Examples

Forest connectivity in central Maryland

Land cover datasets are commonly used to parameterize resistance for connectivity modeling. This example uses the National Land Cover Dataset for the United States to model forest connectivity in central Maryland. Each value in the categorical land cover dataset is assigned a resistance score. We can have Omniscape.jl assign these values internally by providing a reclassification table (see Resistance Reclassification).

First, install the necessary packages and import them:

using Pkg; Pkg.add(["Omniscape", "GeoData", "Plots"])
using Omniscape, GeoData, Plots

Next, download the landcover data we'll use in this example, and plot it:

url_base = "https://raw.githubusercontent.com/Circuitscape/datasets/main/"
# Download the NLCD tile used to create the resistance surface and load it
download(string(url_base, "data/nlcd_2016_frederick_md.tif"),
         "nlcd_2016_frederick_md.tif")

# Plot the landcover data
values = [11, 21, 22, 23, 24, 31, 41, 42, 43, 52, 71, 81, 82, 90, 95]
palette = ["#476BA0", "#DDC9C9", "#D89382", "#ED0000", "#AA0000",
           "#b2b2b2", "#68AA63", "#1C6330", "#B5C98E", "#CCBA7C",
           "#E2E2C1", "#DBD83D", "#AA7028", "#BAD8EA", "#70A3BA"]

plot(GeoArray(GDALarray("nlcd_2016_frederick_md.tif")),
     title = "Land Cover Type", xlabel = "Easting", ylabel = "Northing",
     seriescolor = cgrad(palette, (values .- 12) ./ 84, categorical = true))

Now, load the array using Omniscape's internal read_raster() function or a function from a GIS Julia package of your choice. read_raster() returns a tuple with the data array, a wkt string containing geographic projection info, and an array containing geotransform values. We'll use the wkt and geotransform later.

land_cover, wkt, transform = Omniscape.read_raster("nlcd_2016_frederick_md.tif", Float64)
(Union{Missing, Float64}[82.0 82.0 … 43.0 43.0; 82.0 82.0 … 43.0 43.0; … ; 82.0 81.0 … 82.0 82.0; 81.0 81.0 … 22.0 82.0], "PROJCS[\"Albers_Conical_Equal_Area\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Albers_Conic_Equal_Area\"],PARAMETER[\"latitude_of_center\",0],PARAMETER[\"longitude_of_center\",0],PARAMETER[\"standard_parallel_1\",29.5],PARAMETER[\"standard_parallel_2\",45.5],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]", [1.56369e6, 30.0, 0.0, 1.98666e6, 0.0, -30.0])

The next step is to create a resistance reclassification table that defines a resistance value for each land cover value. Land cover values go in the left column, and resistance values go in the right column. In this case, we are modeling forest connectivity, so forest classes receive the lowest resistance score of one. Other "natural" land cover types are assigned moderate values, and human-developed land cover types are assigned higher values. Medium- to high-intensity development are given a value of missing, which denotes infinite resistance (absolute barriers to movement).

# Create the reclassification table used to translate land cover into resistance
reclass_table = [
    11.	100; # Water
    21	500; # Developed, open space
    22	1000; # Developed, low intensity
    23	missing; # Developed, medium intensity
    24	missing; # Developed, high intensity
    31	100; # Barren land
    41	1; # Deciduous forest
    42	1; # Evergreen forest
    43	1; # Mixed forest
    52	20; # Shrub/scrub
    71	30; # Grassland/herbaceous
    81	200; # Pasture/hay
    82	300; # Cultivated crops
    90	20; # Woody wetlands
    95	30; # Emergent herbaceous wetlands
]
15×2 Array{Union{Missing, Float64},2}:
 11.0   100.0
 21.0   500.0
 22.0  1000.0
 23.0      missing
 24.0      missing
 31.0   100.0
 41.0     1.0
 42.0     1.0
 43.0     1.0
 52.0    20.0
 71.0    30.0
 81.0   200.0
 82.0   300.0
 90.0    20.0
 95.0    30.0

Next, we define the configuration options for this model run. See the Arguments section in the User Guide for more information about each of the configuration options.

# Specify the configuration options
config = Dict{String, String}(
    "radius" => "100",
    "block_size" => "21",
    "project_name" => "md_nlcd_omniscape_output",
    "source_from_resistance" => "true",
    "r_cutoff" => "1", # Only forest pixels should be sources
    "reclassify_resistance" => "true"
)
Dict{String,String} with 6 entries:
  "radius"                 => "100"
  "source_from_resistance" => "true"
  "block_size"             => "21"
  "r_cutoff"               => "1"
  "project_name"           => "md_nlcd_omniscape_output"
  "reclassify_resistance"  => "true"

Finally, compute connectivity using run_omniscape(), feeding in the configuration dictionary, the resistance array, the reclass table, as well as the wkt and geotransform information loaded earlier. Passing in the wkt and geotransform, along with true for the write_outputs argument, will allow Omniscape to write the outputs as properly projected rasters. run_omniscape will print some information to the console and show progress, along with an ETA, in the form of a progress bar.

output = run_omniscape(config,
                       land_cover,
                       reclass_table = reclass_table,
                       wkt = wkt,
                       geotransform = transform,
                       write_outputs = true)
694×725 Array{Union{Missing, Float64},2}:
 0.0027259   0.00949946   0.0171301  …   0.967888  0.770175  0.558886
 0.0087916   0.0168562    0.0276003      1.47789   1.19632   0.810082
 0.0126561   0.023666     0.0364483      1.81501   1.56552   1.0461
 0.0128648   0.0248707    0.042003       2.47138   2.0677    1.17163
 0.162925    0.902127     2.75012        3.09767   3.05185   0.330843
 1.85257     3.87142      4.94455    …   4.08826   1.72065   0.846308
 3.64995     5.62095      6.90531        4.31015   3.03374   1.50058
 5.35248     7.9939       8.57694        5.68669   4.61169   0.48815
 7.08709    10.514       11.0659         9.93759   2.62505   0.0954267
 9.00027    13.0939      13.5202        10.6486    2.12736   0.073892
 ⋮                                   ⋱                       
 0.826616    1.19082      1.23115    …   1.12573   0.864408  4.12477
 0.740572    1.03662      1.06037        1.53066   2.46731   1.36385
 0.733328    0.859761     0.902033       1.66945   1.71243   1.29806
 0.640941    0.759785     0.776764       1.11693   1.66565   1.16405
 0.554711    0.642528     0.662308       1.33137   1.42498   1.0173
 0.467123    0.545933     0.539674   …   1.32905   1.30457   0.867044
 0.351557    0.495893     0.428266       1.41289   1.19599   0.706723
 0.223074    0.621662     0.339845       1.77895   0.983626  0.443013
 0.0682433   0.124514     1.31383        0.67803   0.375597  0.137477

You'll see that outputs are written to a new folder called "md_nlcd_omniscape_output". This is specified by the "project_name" value in config above. The cumulative current map will always be called "cum_currmap.tif", and it will be located in the output folder.

Now, load the current map back into Julia as spatial data and plot it:

current = GDALarray("md_nlcd_omniscape_output/cum_currmap.tif")
plot(current,
     title = "Cumulative Current Flow", xlabel = "Easting", ylabel = "Northing",
     seriescolor = cgrad(:inferno, [0, 0.005, 0.03, 0.06, 0.1, 0.15]))

Cumulative current flow representing foret connectivity. Note that areas in white correspond to built up areas (NLCD values of 23 and 24) that act as absolute barriers to movement.