|
| 1 | +# Examples |
| 2 | + |
| 3 | +## Habitat Corridor Mapping |
| 4 | + |
| 5 | +First, install the necessary packages and import them: |
| 6 | + |
| 7 | +```julia |
| 8 | +using Pkg |
| 9 | +Pkg.add(["Rasters", "Plots", "Graphs", "SimpleWeightedGraphs", "SpatialGraphs"]) |
| 10 | +using SpatialGraphs, Rasters, Plots, Graphs, SimpleWeightedGraphs |
| 11 | +``` |
| 12 | + |
| 13 | +```@setup corridors |
| 14 | +using Pkg |
| 15 | +Pkg.add(["Rasters", "Graphs", "SimpleWeightedGraphs"]) |
| 16 | +using SpatialGraphs, Rasters, Graphs, SimpleWeightedGraphs |
| 17 | +
|
| 18 | +url_base = "https://raw.githubusercontent.com/Circuitscape/datasets/main/" |
| 19 | +download(string(url_base, "data/nlcd_2016_frederick_md.tif"), |
| 20 | + "nlcd_2016_frederick_md.tif") |
| 21 | +
|
| 22 | +land_cover = Raster("nlcd_2016_frederick_md.tif") |
| 23 | +nothing |
| 24 | +``` |
| 25 | + |
| 26 | +Now, let's download the land cover map that we'll use as the basis for defining |
| 27 | +the weights of our graph, and plot it. Graph weights correspond to the cost of |
| 28 | +an animal moving from one graph vertex to another. This is also commonly refered |
| 29 | +to as resistance in the context of connectivity modeling. |
| 30 | + |
| 31 | +```julia |
| 32 | +url_base = "https://raw.githubusercontent.com/Circuitscape/datasets/main/" |
| 33 | +download(string(url_base, "nlcd_2016_frederick_md.tif"), |
| 34 | + "nlcd_2016_frederick_md.tif") |
| 35 | +values = [11, 21, 22, 23, 24, 31, 41, 42, 43, 52, 71, 81, 82, 90, 95] |
| 36 | +palette = ["#476BA0", "#DDC9C9", "#D89382", "#ED0000", "#AA0000", |
| 37 | + "#b2b2b2", "#68AA63", "#1C6330", "#B5C98E", "#CCBA7C", |
| 38 | + "#E2E2C1", "#DBD83D", "#AA7028", "#BAD8EA", "#70A3BA"] |
| 39 | + |
| 40 | +plot(Raster("nlcd_2016_frederick_md.tif"), title = "Land Cover Type", |
| 41 | + xlabel = "Easting", ylabel = "Northing", |
| 42 | + seriescolor = cgrad(palette, (values .- 12) ./ 84, |
| 43 | + categorical = true), size = (700, 640)) |
| 44 | +``` |
| 45 | +```@raw html |
| 46 | +<img src='../figs/mdlc.png' width=500><br> |
| 47 | +``` |
| 48 | + |
| 49 | +Load the land cover data as a Raster, convert it to Float64, and read it into |
| 50 | +memory. |
| 51 | + |
| 52 | +```@example corridors |
| 53 | +landcover = convert.( |
| 54 | + Float64, |
| 55 | + read( |
| 56 | + Raster( |
| 57 | + "nlcd_2016_frederick_md.tif", |
| 58 | + missingval = -9999.0 |
| 59 | + ) |
| 60 | + ) |
| 61 | +) |
| 62 | +
|
| 63 | +``` |
| 64 | + |
| 65 | +Now, we need to convert land cover into resistance by reclassifying the values |
| 66 | +in the land cover raster to their corresponding costs/weights. |
| 67 | + |
| 68 | +```@example corridors |
| 69 | +# Copy the landcover object to use as a place holder for resistance |
| 70 | +resistance = deepcopy(landcover) |
| 71 | +
|
| 72 | +# Define how values should be reclassified. Original values are on the left, |
| 73 | +# and new values are on the right. |
| 74 | +reclass_table = [ |
| 75 | + 11. 100; # Water |
| 76 | + 21 500; # Developed, open space |
| 77 | + 22 1000; # Developed, low intensity |
| 78 | + 23 -9999; # Developed, medium intensity |
| 79 | + 24 -9999; # Developed, high intensity |
| 80 | + 31 100; # Barren land |
| 81 | + 41 1; # Deciduous forest |
| 82 | + 42 1; # Evergreen forest |
| 83 | + 43 1; # Mixed forest |
| 84 | + 52 20; # Shrub/scrub |
| 85 | + 71 30; # Grassland/herbaceous |
| 86 | + 81 200; # Pasture/hay |
| 87 | + 82 300; # Cultivated crops |
| 88 | + 90 20; # Woody wetlands |
| 89 | + 95 30; # Emergent herbaceous wetlands |
| 90 | +] |
| 91 | +
|
| 92 | +# Redefine the values in resistance based on the corresponding land cover |
| 93 | +# values, using reclass_table |
| 94 | +for i in 1:(size(reclass_table)[1]) |
| 95 | + resistance[landcover .== reclass_table[i, 1]] .= reclass_table[i, 2] |
| 96 | +end |
| 97 | +``` |
| 98 | + |
| 99 | +Now that we have a resistance raster, we can start using SpatialGraphs.jl! |
| 100 | + |
| 101 | +Convert the resistance raster into a `WeightedRasterGraph`: |
| 102 | + |
| 103 | +```@example corridors |
| 104 | +wrg = weightedrastergraph(resistance) |
| 105 | +``` |
| 106 | + |
| 107 | +Now that we have a `WeightedRasterGraph`, we can make use of the functions in |
| 108 | +Graphs.jl. |
| 109 | + |
| 110 | +Let's calculate cost distances from two different vertices. Vertex IDs, |
| 111 | +the second argument to `dijkstra_shortest_paths`, correspond to the values in |
| 112 | +wrg.vertex_raster. |
| 113 | + |
| 114 | +First, we'll use `dijkstra_shortest_paths` to calculate to total cost distance |
| 115 | +from each vertex of interest to all other vertices in the graph |
| 116 | + |
| 117 | +```@example corridors |
| 118 | +cd_1 = dijkstra_shortest_paths(wrg, 1) |
| 119 | +cd_2 = dijkstra_shortest_paths(wrg, 348021) |
| 120 | +``` |
| 121 | + |
| 122 | +Now, create placeholder rasters to store the cost distance values for each |
| 123 | +vertex, and populate them using the cost distances calculated |
| 124 | +above. Finally, we'll sum the |
| 125 | + |
| 126 | +```@example corridors |
| 127 | +# Create rasters |
| 128 | +cd_1_raster = Raster( # Start with an empty raster that will be populated |
| 129 | + fill(float(resistance.missingval), size(wrg.vertex_raster)), |
| 130 | + dims(resistance), |
| 131 | + missingval=resistance.missingval |
| 132 | +) |
| 133 | +cd_2_raster = Raster( |
| 134 | + fill(float(wrg.vertex_raster.missingval), size(wrg.vertex_raster)), |
| 135 | + dims(resistance), |
| 136 | + missingval=resistance.missingval |
| 137 | +) |
| 138 | +
|
| 139 | +# Populate with cost distances |
| 140 | +cd_1_raster[wrg.vertex_raster .!= wrg.vertex_raster.missingval] .= cd_1.dists[ |
| 141 | + wrg.vertex_raster[wrg.vertex_raster.!= wrg.vertex_raster.missingval] |
| 142 | +] |
| 143 | +cd_2_raster[wrg.vertex_raster .!= wrg.vertex_raster.missingval] .= cd_2.dists[ |
| 144 | + wrg.vertex_raster[wrg.vertex_raster.!= wrg.vertex_raster.missingval] |
| 145 | +] |
| 146 | +
|
| 147 | +# Sum the rasters to enalbe corridor mapping |
| 148 | +cwd_sum = cd_1_raster + cd_2_raster |
| 149 | +``` |
| 150 | + |
| 151 | +Setting some land cover classes to NoData resulted in some vertices being |
| 152 | +completely isolated from the vertices from which we calculated cost distances |
| 153 | +above. In these cases, `dijkstra_shortest_paths` returns a value of `Inf` for |
| 154 | +the cost distance. Let's get rid of those by setting them to the NoData value |
| 155 | +of our raster. |
| 156 | + |
| 157 | +```@example corridors |
| 158 | +cwd_sum[cwd_sum .== Inf] .= -9999. |
| 159 | +``` |
| 160 | + |
| 161 | +Finally, let's plot the result! We'll use a threshold the set pixels above |
| 162 | +a certain value to NoData, so we will only retain the pixels within our |
| 163 | +"corridor". |
| 164 | + |
| 165 | +```julia |
| 166 | +corridor = copy(cwd_sum) |
| 167 | +corridor[cwd_sum .>= (minimum(cwd_sum[cwd_sum .>= 0]) + 250)] .= -9999 |
| 168 | + |
| 169 | +plot( |
| 170 | + cwd_sum, |
| 171 | + title = "Cost Distance Corridor", |
| 172 | + xlabel = "Easting", ylabel = "Northing", |
| 173 | + colorbar_title = " \nCost Distance Sum", |
| 174 | + size = (650, 475) |
| 175 | +) |
| 176 | +``` |
| 177 | + |
| 178 | +```@raw html |
| 179 | +<img src='../figs/corridor_plot.png' width=500><br> |
| 180 | +``` |
0 commit comments