diff --git a/docs/index.html b/docs/index.html index 9936d7b..e630874 100644 --- a/docs/index.html +++ b/docs/index.html @@ -167,7 +167,9 @@

On this page

@@ -195,384 +197,609 @@

Land Cover Analysis

Loading data

-

We load the possible sites (quiet = TRUE is for not displaying verbose loading information). We also load the two halves of the far north land cover dataset, along with the the attribute table of land cover classes.

+

We load the possible sites (quiet = TRUE is for not displaying verbose loading information). We also load the two halves of the far north land cover dataset, along with the the attribute table of land cover classes. We load the ecodistrict data and select for the relevant lowlands disctrict, coded as 1028.

sites_possible <- sf::st_read(
   "data/Peat/GRTS_PossibleCaARU_sample_draw_base.shp", 
   quiet = TRUE)
-lu_16 <- raster::raster("data/land_use/FarNorth_LandCover_Class_UTM16.tif")
-lu_17 <- raster::raster("data/land_use/FarNorth_LandCover_Class_UTM17.tif")
-lu_dat <- readr::read_csv("data/land_use/attr_table_northen_ont_lc.txt") |>
-  dplyr::mutate(cats = as.factor(code))
+ecodistrict <- sf::st_read( + "data/ecodistrict_shp/Ecodistricts/ecodistricts.shp", + quiet = TRUE) |> + dplyr::filter(ECODISTRIC == 1028) +lu_16 <- raster::raster("data/land_use/FarNorth_LandCover_Class_UTM16.tif") +lu_17 <- raster::raster("data/land_use/FarNorth_LandCover_Class_UTM17.tif") +lu_dat <- readr::read_csv("data/land_use/attr_table_northen_ont_lc.txt") |> + dplyr::mutate(cats = as.factor(code))
+
+

Plotting spatial data

+

It is always a good idea to try and plot spatial data before any processing.

+
+
ggplot() +
+  geom_sf(data = ecodistrict) +
+  geom_sf(data = sf::st_transform(sites_possible, 
+                                  sf::st_crs(ecodistrict))) +
+  theme_bw()
+
+
+
+

+
+
+
+
+

Plotting the land cover data is difficult because it is provided is two different UTMs.

+

Extracting Land Cover data

-

The following functions will take care of land cover extraction.

+

The following functions will take care of land cover extraction for sites.

-
extract_from_points <- function(scale_m, sites, lu) {
-  
-  sites_buffer <- sites |>
-    sf::st_transform(sf::st_crs(lu)) |> 
-    sf::st_buffer(dist = scale_m) |> 
-    dplyr::select(siteID)
-  
-  extr <- exactextractr::exact_extract(lu, sites_buffer,
-                                       progress = FALSE)
-  
-  extr <- mapply(extr, 1:length(extr), 
-                 FUN = \(x, y) dplyr::mutate(x, id = y),
-                 SIMPLIFY = F)
-  
-  extr_df <- do.call(rbind, extr) |>
-    dplyr::filter(!is.na(value)) |>
-    dplyr::relocate(id)
-  
-  return(extr_df)
-   
-}
-
-compute_land_cover <- function(scale_m, sites, 
-                               lu_16, lu_17, lu_dat) {
-  
-  extr_16_df <- extract_from_points(scale_m, sites, lu_16)
-  extr_17_df <- extract_from_points(scale_m, sites, lu_17)
-  
-  stopifnot(all(!(extr_16_df$siteID %in% extr_17_df$siteID)))
-  
-  extr <- rbind(extr_16_df, extr_17_df) |>
-    dplyr::arrange(id, value)
-  
-  extr_table <- extr |>
-    dplyr::group_by(id, value) |>
-    dplyr::summarise(coverage_fraction_sum = sum(coverage_fraction)) |>
-    dplyr::mutate(props = 
-                    coverage_fraction_sum/sum(coverage_fraction_sum)) |>
-    dplyr::ungroup()  |>
-    dplyr::mutate(value = as.factor(value)) |>
-    dplyr::left_join(lu_dat, by = c("value" = "cats")) |>
-    dplyr::select(id, category_code, props, label)
-  
-  extr_table[is.na(extr_table)] <- 0
-  
-  extr_table_sum <- extr_table |>
-    dplyr::group_by(category_code, label) |>
-    dplyr::summarise(prop_sum = sum(props)) |>
-    dplyr::ungroup()
-  
-  return(extr_table_sum)
-
-}
+
extract_from_points <- function(scale_m, sites, lu) {
+  
+  sites_buffer <- sites |>
+    sf::st_transform(sf::st_crs(lu)) |> 
+    sf::st_buffer(dist = scale_m) |> 
+    dplyr::select(siteID)
+  
+  extr <- exactextractr::exact_extract(lu, sites_buffer,
+                                       progress = FALSE)
+  
+  extr <- mapply(extr, 1:length(extr), 
+                 FUN = \(x, y) dplyr::mutate(x, id = y),
+                 SIMPLIFY = F)
+  
+  extr_df <- do.call(rbind, extr) |>
+    dplyr::filter(!is.na(value)) |>
+    dplyr::relocate(id)
+  
+  return(extr_df)
+  
+}
+
+compute_land_cover <- function(scale_m, sites, 
+                               lu_16, lu_17, lu_dat) {
+  
+  extr_16_df <- extract_from_points(scale_m, sites, lu_16)
+  extr_17_df <- extract_from_points(scale_m, sites, lu_17)
+  
+  stopifnot(all(!(extr_16_df$siteID %in% extr_17_df$siteID)))
+  
+  extr <- rbind(extr_16_df, extr_17_df) |>
+    dplyr::arrange(id, value)
+  
+  extr_table <- extr |>
+    dplyr::group_by(value) |>
+    dplyr::summarise(coverage_fraction_sum = sum(coverage_fraction)) |>
+    dplyr::mutate(prop = 
+                    coverage_fraction_sum/sum(coverage_fraction_sum)) |>
+    dplyr::ungroup()  |>
+    dplyr::mutate(value = as.factor(value)) |>
+    dplyr::left_join(lu_dat, by = c("value" = "cats")) |>
+    dplyr::select(category_code, prop, label)
+  
+  extr_table[is.na(extr_table)] <- 0
+  
+  return(extr_table)
+  
+}

We extract at different scales (buffer radius around points): 1 m, 50 m, 100 m and 1 km.

-
res <- mapply(FUN = compute_land_cover, c(`1 m` = 1, `50 m` = 50,
-                                          `100 m` = 100, `1 km` = 1000),
-              MoreArgs = list(
-                sites = sites_possible,
-                lu_16 = lu_16, lu_17 = lu_17, lu_dat = lu_dat),
-              SIMPLIFY = F) |>
-  dplyr::bind_rows(.id = 'scale') |>
-  dplyr::mutate(scale = forcats::fct_relevel(scale, "1 m", "50 m",
-                                             "100 m", "1 km"),
-                label =  forcats::fct_reorder(label, prop_sum)) |> 
-  dplyr::arrange(scale, dplyr::desc(prop_sum))
-
-knitr::kable(res)
+
res_points <- mapply(FUN = compute_land_cover, 
+                     c(`1 m` = 1, `50 m` = 50,
+                       `100 m` = 100, `1 km` = 1000),
+                     MoreArgs = list(
+                       sites = sites_possible,
+                       lu_16 = lu_16, lu_17 = lu_17, lu_dat = lu_dat),
+                     SIMPLIFY = F) |>
+  dplyr::bind_rows(.id = 'scale') |>
+  dplyr::mutate(scale = forcats::fct_relevel(scale, "1 m", "50 m",
+                                             "100 m", "1 km"),
+                label =  forcats::fct_reorder(label, prop)) |> 
+  dplyr::arrange(scale, dplyr::desc(prop))
+
+knitr::kable(res_points)
+ - + - + - + - + - + - + - - - - + + + - - - + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + +
scale category_codeprop labelprop_sum
1 m TrBOG0.26 Treed Bog13.00
1 m TrFEN0.23 Treed Fen11.34
1 m OBOG0.20 Open Bog10.00
1 m ConSWA0.11 Coniferous Swamp5.66
1 m ConTRE0.10 Coniferous Treed5.00
1 m OFEN0.06 Open Fen3.00
1 mMixTREMixed Treed1.00SpTRE0.02Sparse Treed
1 mSpTRESparse Treed1.00MixTRE0.02Mixed Treed
50 m TrBOG0.23 Treed Bog11.69
50 m TrFEN0.22 Treed Fen10.93
50 m OBOG0.21 Open Bog10.56
50 m ConSWA0.12 Coniferous Swamp5.81
50 m ConTRE0.09 Coniferous Treed4.50
50 m OFEN0.07 Open Fen3.71
50 m SpTRE0.03 Sparse Treed1.38
50 m MixTRE0.02 Mixed Treed0.76
50 m WAT0.01 Clear Open Water0.61
50 m XWAT0.00 Turbid Water0.05
100 m TrBOG0.23 Treed Bog11.32
100 m TrFEN0.22 Treed Fen10.83
100 m OBOG0.20 Open Bog9.82
100 m ConSWA0.12 Coniferous Swamp6.03
100 m ConTRE0.10 Coniferous Treed4.84
100 m OFEN0.08 Open Fen4.07
100 m WAT0.03 Clear Open Water1.36
100 m SpTRE0.02 Sparse Treed0.80
100 m MixTRE0.01 Mixed Treed0.67
100 m XWAT0.01 Turbid Water0.25
1 km TrFEN0.21 Treed Fen10.69
1 km TrBOG0.21 Treed Bog10.59
1 km OBOG0.19 Open Bog9.50
1 km ConSWA0.17 Coniferous Swamp8.33
1 km OFEN0.07 Open Fen3.74
1 km ConTRE0.06 Coniferous Treed3.21
1 km WAT0.04 Clear Open Water2.21
1 km MixTRE0.01 Mixed Treed0.66
1 km XWAT0.01 Turbid Water0.41
1 km ThSWA0.01 Thicket Swamp0.35
1 km SpTRE0.00 Sparse Treed0.14
1 km NSWood0.00 Disturbance - Non and Sparse Woody0.12
1 km DecTRE0.00 Deciduous Treed0.02
1 km TrOrSHr0.00 Disturbance - Treed and/or Shrub0.01
1 km BED0.00 Bedrock
+
+
+

We also want to do the same operation for the ecodistrict to allow for comparison. We don’t need to use exact extraction, insteadt the crop and mask each raster. This operation is costly so we write out the rasters and load them again (see unrendered code).

+
+
ecodistrict_16 <- sf::st_transform(ecodistrict, sf::st_crs(lu_16))
+ecodistrict_17 <- sf::st_transform(ecodistrict, sf::st_crs(lu_17))
+
+lu_16_crop <- raster::crop(lu_16, ecodistrict_16)
+lu_16_crop_mask <- raster::mask(lu_16_crop, ecodistrict_16)
+
+lu_17_crop <- raster::crop(lu_17, ecodistrict_17)
+lu_17_crop_mask <- raster::mask(lu_17_crop, ecodistrict_17)
+
+

We can then get the frequencies of values. This operation is also costly so we write out the objects and load them again (see unrendered code).

+
+
lu_16_freq <- raster::freq(lu_16_crop_mask)
+lu_17_freq <- raster::freq(lu_17_crop_mask)
+
+

We combine the results of both UTMs.

+
+
res_ecodistrict <- rbind(lu_16_freq, lu_17_freq) |>
+  as.data.frame() |> 
+  dplyr::group_by(value) |> 
+  dplyr::summarise(count = sum(count)) |> 
+  dplyr::ungroup() |> 
+  dplyr::filter(!is.na(value)) |> 
+  dplyr::mutate(prop = count/sum(count)) |> 
+  dplyr::mutate(value = as.factor(value)) |>
+  dplyr::left_join(lu_dat, by = c("value" = "cats")) |>
+  dplyr::filter(!is.na(label)) |> 
+  dplyr::select(category_code, prop, label) |> 
+  dplyr::mutate(scale = "Ecodistrict") |> 
+  dplyr::relocate(scale) |>
+  dplyr::arrange(scale, dplyr::desc(prop))
+
+knitr::kable(res_ecodistrict)
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
scalecategory_codeproplabel
EcodistrictTrFEN0.25Treed Fen
EcodistrictOBOG0.20Open Bog
EcodistrictTrBOG0.19Treed Bog
EcodistrictConSWA0.12Coniferous Swamp
EcodistrictOFEN0.09Open Fen
EcodistrictWAT0.07Clear Open Water
EcodistrictConTRE0.04Coniferous Treed
EcodistrictNSWood0.01Disturbance - Non and Sparse Woody
EcodistrictTrOrSHr0.01Disturbance - Treed and/or Shrub
EcodistrictMixTRE0.01Mixed Treed
EcodistrictThSWA0.01Thicket Swamp
EcodistrictSpTRE0.00Sparse Treed
EcodistrictXWAT0.00Turbid Water
EcodistrictDecTRE0.00Deciduous Treed
EcodistrictMIN0.00Sand/Gravel/Mine Tailings
EcodistrictFrMAR0.00Freshwater Marsh
EcodistrictBED0.00Bedrock
EcodistrictURB0.00Community/Infrastructure
EcodistrictDecSWA 0.00Deciduous Swamp
EcodistrictInMAR0.00Intertidal Marsh
+

And then combine the results between scales and utm.

+
+
res <- rbind(res_points, res_ecodistrict) |> 
+  tidyr::complete(scale, label) |> 
+  tidyr::replace_na(list(prop = 0)) |> 
+  dplyr::mutate(label =  forcats::fct_reorder(label, prop))
+
+
+
+

Results

We can plot the results with “dodged” ggplot2 barplots.

-
library(ggplot2)
-
-my_pal <- c('#c7e9b4','#7fcdbb','#41b6c4','#2c7fb8','#253494')
-
-ggplot(res) +
-  geom_bar(aes(x = label, y = prop_sum, fill = scale, colour = scale), 
-           alpha = 0.8,
-           stat = "identity",
-           position = "dodge") +
-  theme_bw() +
-  theme(axis.text.x = element_text(angle = 20, vjust = 1, hjust = 1)) +
-  labs(x = "Land Use Class", y = "Sum of Proportions",
-       fill = "Scale", colour = "Scale") +
-  scale_fill_manual(values = my_pal) +
-  scale_color_manual(values = my_pal)
+
my_pal <- c('#c7e9b4','#7fcdbb','#41b6c4','#1d91c0','#225ea8','#0c2c84')
+
+ggplot(res) +
+  geom_bar(aes(x = label, y = prop, fill = scale, colour = scale), 
+           alpha = 0.8,
+           stat = "identity",
+           position = "dodge") +
+  theme_bw() +
+  theme(axis.text.x = element_text(angle = 20, vjust = 1, hjust = 1)) +
+  labs(x = "Land Use Class", y = "Sum of Proportions",
+       fill = "Scale", colour = "Scale") +
+  scale_fill_manual(values = my_pal) +
+  scale_color_manual(values = my_pal)
+
+
+
+

+
+
+
+
+

Removing the land use classes than are not present around sites, we get a slightly easier graph to read.

+
+
only_at_sites <- res |> 
+  dplyr::filter(prop > 0) |> 
+  dplyr::filter(scale != "Ecodistrict") |> 
+  dplyr::pull(label)
+
+res_filt <- res |> 
+  dplyr::filter(label %in% only_at_sites)
+
+ggplot(res_filt) +
+  geom_bar(aes(x = label, y = prop, fill = scale, colour = scale), 
+           alpha = 0.8,
+           stat = "identity",
+           position = "dodge") +
+  theme_bw() +
+  theme(axis.text.x = element_text(angle = 20, vjust = 1, hjust = 1)) +
+  labs(x = "Land Use Class", y = "Sum of Proportions",
+       fill = "Scale", colour = "Scale") +
+  scale_fill_manual(values = my_pal) +
+  scale_color_manual(values = my_pal)
-

+

diff --git a/docs/index_files/figure-html/unnamed-chunk-14-1.png b/docs/index_files/figure-html/unnamed-chunk-14-1.png new file mode 100644 index 0000000..9e19fc5 Binary files /dev/null and b/docs/index_files/figure-html/unnamed-chunk-14-1.png differ diff --git a/docs/index_files/figure-html/unnamed-chunk-15-1.png b/docs/index_files/figure-html/unnamed-chunk-15-1.png new file mode 100644 index 0000000..9e97453 Binary files /dev/null and b/docs/index_files/figure-html/unnamed-chunk-15-1.png differ diff --git a/docs/index_files/figure-html/unnamed-chunk-3-1.png b/docs/index_files/figure-html/unnamed-chunk-3-1.png new file mode 100644 index 0000000..0a07a6c Binary files /dev/null and b/docs/index_files/figure-html/unnamed-chunk-3-1.png differ diff --git a/docs/search.json b/docs/search.json index e6131d0..df6043b 100644 --- a/docs/search.json +++ b/docs/search.json @@ -4,7 +4,7 @@ "href": "index.html", "title": "Land Cover Analysis", "section": "", - "text": "We load the possible sites (quiet = TRUE is for not displaying verbose loading information). We also load the two halves of the far north land cover dataset, along with the the attribute table of land cover classes.\n\nsites_possible <- sf::st_read(\n \"data/Peat/GRTS_PossibleCaARU_sample_draw_base.shp\", \n quiet = TRUE)\nlu_16 <- raster::raster(\"data/land_use/FarNorth_LandCover_Class_UTM16.tif\")\nlu_17 <- raster::raster(\"data/land_use/FarNorth_LandCover_Class_UTM17.tif\")\nlu_dat <- readr::read_csv(\"data/land_use/attr_table_northen_ont_lc.txt\") |>\n dplyr::mutate(cats = as.factor(code))", + "text": "We load the possible sites (quiet = TRUE is for not displaying verbose loading information). We also load the two halves of the far north land cover dataset, along with the the attribute table of land cover classes. We load the ecodistrict data and select for the relevant lowlands disctrict, coded as 1028.\n\nsites_possible <- sf::st_read(\n \"data/Peat/GRTS_PossibleCaARU_sample_draw_base.shp\", \n quiet = TRUE)\necodistrict <- sf::st_read(\n \"data/ecodistrict_shp/Ecodistricts/ecodistricts.shp\", \n quiet = TRUE) |> \n dplyr::filter(ECODISTRIC == 1028)\nlu_16 <- raster::raster(\"data/land_use/FarNorth_LandCover_Class_UTM16.tif\")\nlu_17 <- raster::raster(\"data/land_use/FarNorth_LandCover_Class_UTM17.tif\")\nlu_dat <- readr::read_csv(\"data/land_use/attr_table_northen_ont_lc.txt\") |>\n dplyr::mutate(cats = as.factor(code))", "crumbs": [ "Site Analysis", "Land Cover Analysis" @@ -15,7 +15,7 @@ "href": "index.html#loading-data", "title": "Land Cover Analysis", "section": "", - "text": "We load the possible sites (quiet = TRUE is for not displaying verbose loading information). We also load the two halves of the far north land cover dataset, along with the the attribute table of land cover classes.\n\nsites_possible <- sf::st_read(\n \"data/Peat/GRTS_PossibleCaARU_sample_draw_base.shp\", \n quiet = TRUE)\nlu_16 <- raster::raster(\"data/land_use/FarNorth_LandCover_Class_UTM16.tif\")\nlu_17 <- raster::raster(\"data/land_use/FarNorth_LandCover_Class_UTM17.tif\")\nlu_dat <- readr::read_csv(\"data/land_use/attr_table_northen_ont_lc.txt\") |>\n dplyr::mutate(cats = as.factor(code))", + "text": "We load the possible sites (quiet = TRUE is for not displaying verbose loading information). We also load the two halves of the far north land cover dataset, along with the the attribute table of land cover classes. We load the ecodistrict data and select for the relevant lowlands disctrict, coded as 1028.\n\nsites_possible <- sf::st_read(\n \"data/Peat/GRTS_PossibleCaARU_sample_draw_base.shp\", \n quiet = TRUE)\necodistrict <- sf::st_read(\n \"data/ecodistrict_shp/Ecodistricts/ecodistricts.shp\", \n quiet = TRUE) |> \n dplyr::filter(ECODISTRIC == 1028)\nlu_16 <- raster::raster(\"data/land_use/FarNorth_LandCover_Class_UTM16.tif\")\nlu_17 <- raster::raster(\"data/land_use/FarNorth_LandCover_Class_UTM17.tif\")\nlu_dat <- readr::read_csv(\"data/land_use/attr_table_northen_ont_lc.txt\") |>\n dplyr::mutate(cats = as.factor(code))", "crumbs": [ "Site Analysis", "Land Cover Analysis" @@ -26,7 +26,7 @@ "href": "index.html#extracting-land-cover-data", "title": "Land Cover Analysis", "section": "Extracting Land Cover data", - "text": "Extracting Land Cover data\nThe following functions will take care of land cover extraction.\n\nextract_from_points <- function(scale_m, sites, lu) {\n \n sites_buffer <- sites |>\n sf::st_transform(sf::st_crs(lu)) |> \n sf::st_buffer(dist = scale_m) |> \n dplyr::select(siteID)\n \n extr <- exactextractr::exact_extract(lu, sites_buffer,\n progress = FALSE)\n \n extr <- mapply(extr, 1:length(extr), \n FUN = \\(x, y) dplyr::mutate(x, id = y),\n SIMPLIFY = F)\n \n extr_df <- do.call(rbind, extr) |>\n dplyr::filter(!is.na(value)) |>\n dplyr::relocate(id)\n \n return(extr_df)\n \n}\n\ncompute_land_cover <- function(scale_m, sites, \n lu_16, lu_17, lu_dat) {\n \n extr_16_df <- extract_from_points(scale_m, sites, lu_16)\n extr_17_df <- extract_from_points(scale_m, sites, lu_17)\n \n stopifnot(all(!(extr_16_df$siteID %in% extr_17_df$siteID)))\n \n extr <- rbind(extr_16_df, extr_17_df) |>\n dplyr::arrange(id, value)\n \n extr_table <- extr |>\n dplyr::group_by(id, value) |>\n dplyr::summarise(coverage_fraction_sum = sum(coverage_fraction)) |>\n dplyr::mutate(props = \n coverage_fraction_sum/sum(coverage_fraction_sum)) |>\n dplyr::ungroup() |>\n dplyr::mutate(value = as.factor(value)) |>\n dplyr::left_join(lu_dat, by = c(\"value\" = \"cats\")) |>\n dplyr::select(id, category_code, props, label)\n \n extr_table[is.na(extr_table)] <- 0\n \n extr_table_sum <- extr_table |>\n dplyr::group_by(category_code, label) |>\n dplyr::summarise(prop_sum = sum(props)) |>\n dplyr::ungroup()\n \n return(extr_table_sum)\n\n}\n\nWe extract at different scales (buffer radius around points): 1 m, 50 m, 100 m and 1 km.\n\nres <- mapply(FUN = compute_land_cover, c(`1 m` = 1, `50 m` = 50,\n `100 m` = 100, `1 km` = 1000),\n MoreArgs = list(\n sites = sites_possible,\n lu_16 = lu_16, lu_17 = lu_17, lu_dat = lu_dat),\n SIMPLIFY = F) |>\n dplyr::bind_rows(.id = 'scale') |>\n dplyr::mutate(scale = forcats::fct_relevel(scale, \"1 m\", \"50 m\",\n \"100 m\", \"1 km\"),\n label = forcats::fct_reorder(label, prop_sum)) |> \n dplyr::arrange(scale, dplyr::desc(prop_sum))\n\nknitr::kable(res)\n\n\n\n\nscale\ncategory_code\nlabel\nprop_sum\n\n\n\n\n1 m\nTrBOG\nTreed Bog\n13.00\n\n\n1 m\nTrFEN\nTreed Fen\n11.34\n\n\n1 m\nOBOG\nOpen Bog\n10.00\n\n\n1 m\nConSWA\nConiferous Swamp\n5.66\n\n\n1 m\nConTRE\nConiferous Treed\n5.00\n\n\n1 m\nOFEN\nOpen Fen\n3.00\n\n\n1 m\nMixTRE\nMixed Treed\n1.00\n\n\n1 m\nSpTRE\nSparse Treed\n1.00\n\n\n50 m\nTrBOG\nTreed Bog\n11.69\n\n\n50 m\nTrFEN\nTreed Fen\n10.93\n\n\n50 m\nOBOG\nOpen Bog\n10.56\n\n\n50 m\nConSWA\nConiferous Swamp\n5.81\n\n\n50 m\nConTRE\nConiferous Treed\n4.50\n\n\n50 m\nOFEN\nOpen Fen\n3.71\n\n\n50 m\nSpTRE\nSparse Treed\n1.38\n\n\n50 m\nMixTRE\nMixed Treed\n0.76\n\n\n50 m\nWAT\nClear Open Water\n0.61\n\n\n50 m\nXWAT\nTurbid Water\n0.05\n\n\n100 m\nTrBOG\nTreed Bog\n11.32\n\n\n100 m\nTrFEN\nTreed Fen\n10.83\n\n\n100 m\nOBOG\nOpen Bog\n9.82\n\n\n100 m\nConSWA\nConiferous Swamp\n6.03\n\n\n100 m\nConTRE\nConiferous Treed\n4.84\n\n\n100 m\nOFEN\nOpen Fen\n4.07\n\n\n100 m\nWAT\nClear Open Water\n1.36\n\n\n100 m\nSpTRE\nSparse Treed\n0.80\n\n\n100 m\nMixTRE\nMixed Treed\n0.67\n\n\n100 m\nXWAT\nTurbid Water\n0.25\n\n\n1 km\nTrFEN\nTreed Fen\n10.69\n\n\n1 km\nTrBOG\nTreed Bog\n10.59\n\n\n1 km\nOBOG\nOpen Bog\n9.50\n\n\n1 km\nConSWA\nConiferous Swamp\n8.33\n\n\n1 km\nOFEN\nOpen Fen\n3.74\n\n\n1 km\nConTRE\nConiferous Treed\n3.21\n\n\n1 km\nWAT\nClear Open Water\n2.21\n\n\n1 km\nMixTRE\nMixed Treed\n0.66\n\n\n1 km\nXWAT\nTurbid Water\n0.41\n\n\n1 km\nThSWA\nThicket Swamp\n0.35\n\n\n1 km\nSpTRE\nSparse Treed\n0.14\n\n\n1 km\nNSWood\nDisturbance - Non and Sparse Woody\n0.12\n\n\n1 km\nDecTRE\nDeciduous Treed\n0.02\n\n\n1 km\nTrOrSHr\nDisturbance - Treed and/or Shrub\n0.01\n\n\n1 km\nBED\nBedrock\n0.00\n\n\n\n\n\nWe can plot the results with “dodged” ggplot2 barplots.\n\nlibrary(ggplot2)\n\nmy_pal <- c('#c7e9b4','#7fcdbb','#41b6c4','#2c7fb8','#253494')\n\nggplot(res) +\n geom_bar(aes(x = label, y = prop_sum, fill = scale, colour = scale), \n alpha = 0.8,\n stat = \"identity\",\n position = \"dodge\") +\n theme_bw() +\n theme(axis.text.x = element_text(angle = 20, vjust = 1, hjust = 1)) +\n labs(x = \"Land Use Class\", y = \"Sum of Proportions\",\n fill = \"Scale\", colour = \"Scale\") +\n scale_fill_manual(values = my_pal) +\n scale_color_manual(values = my_pal)", + "text": "Extracting Land Cover data\nThe following functions will take care of land cover extraction for sites.\n\nextract_from_points <- function(scale_m, sites, lu) {\n \n sites_buffer <- sites |>\n sf::st_transform(sf::st_crs(lu)) |> \n sf::st_buffer(dist = scale_m) |> \n dplyr::select(siteID)\n \n extr <- exactextractr::exact_extract(lu, sites_buffer,\n progress = FALSE)\n \n extr <- mapply(extr, 1:length(extr), \n FUN = \\(x, y) dplyr::mutate(x, id = y),\n SIMPLIFY = F)\n \n extr_df <- do.call(rbind, extr) |>\n dplyr::filter(!is.na(value)) |>\n dplyr::relocate(id)\n \n return(extr_df)\n \n}\n\ncompute_land_cover <- function(scale_m, sites, \n lu_16, lu_17, lu_dat) {\n \n extr_16_df <- extract_from_points(scale_m, sites, lu_16)\n extr_17_df <- extract_from_points(scale_m, sites, lu_17)\n \n stopifnot(all(!(extr_16_df$siteID %in% extr_17_df$siteID)))\n \n extr <- rbind(extr_16_df, extr_17_df) |>\n dplyr::arrange(id, value)\n \n extr_table <- extr |>\n dplyr::group_by(value) |>\n dplyr::summarise(coverage_fraction_sum = sum(coverage_fraction)) |>\n dplyr::mutate(prop = \n coverage_fraction_sum/sum(coverage_fraction_sum)) |>\n dplyr::ungroup() |>\n dplyr::mutate(value = as.factor(value)) |>\n dplyr::left_join(lu_dat, by = c(\"value\" = \"cats\")) |>\n dplyr::select(category_code, prop, label)\n \n extr_table[is.na(extr_table)] <- 0\n \n return(extr_table)\n \n}\n\nWe extract at different scales (buffer radius around points): 1 m, 50 m, 100 m and 1 km.\n\nres_points <- mapply(FUN = compute_land_cover, \n c(`1 m` = 1, `50 m` = 50,\n `100 m` = 100, `1 km` = 1000),\n MoreArgs = list(\n sites = sites_possible,\n lu_16 = lu_16, lu_17 = lu_17, lu_dat = lu_dat),\n SIMPLIFY = F) |>\n dplyr::bind_rows(.id = 'scale') |>\n dplyr::mutate(scale = forcats::fct_relevel(scale, \"1 m\", \"50 m\",\n \"100 m\", \"1 km\"),\n label = forcats::fct_reorder(label, prop)) |> \n dplyr::arrange(scale, dplyr::desc(prop))\n\nknitr::kable(res_points)\n\n\n\n\nscale\ncategory_code\nprop\nlabel\n\n\n\n\n1 m\nTrBOG\n0.26\nTreed Bog\n\n\n1 m\nTrFEN\n0.23\nTreed Fen\n\n\n1 m\nOBOG\n0.20\nOpen Bog\n\n\n1 m\nConSWA\n0.11\nConiferous Swamp\n\n\n1 m\nConTRE\n0.10\nConiferous Treed\n\n\n1 m\nOFEN\n0.06\nOpen Fen\n\n\n1 m\nSpTRE\n0.02\nSparse Treed\n\n\n1 m\nMixTRE\n0.02\nMixed Treed\n\n\n50 m\nTrBOG\n0.23\nTreed Bog\n\n\n50 m\nTrFEN\n0.22\nTreed Fen\n\n\n50 m\nOBOG\n0.21\nOpen Bog\n\n\n50 m\nConSWA\n0.12\nConiferous Swamp\n\n\n50 m\nConTRE\n0.09\nConiferous Treed\n\n\n50 m\nOFEN\n0.07\nOpen Fen\n\n\n50 m\nSpTRE\n0.03\nSparse Treed\n\n\n50 m\nMixTRE\n0.02\nMixed Treed\n\n\n50 m\nWAT\n0.01\nClear Open Water\n\n\n50 m\nXWAT\n0.00\nTurbid Water\n\n\n100 m\nTrBOG\n0.23\nTreed Bog\n\n\n100 m\nTrFEN\n0.22\nTreed Fen\n\n\n100 m\nOBOG\n0.20\nOpen Bog\n\n\n100 m\nConSWA\n0.12\nConiferous Swamp\n\n\n100 m\nConTRE\n0.10\nConiferous Treed\n\n\n100 m\nOFEN\n0.08\nOpen Fen\n\n\n100 m\nWAT\n0.03\nClear Open Water\n\n\n100 m\nSpTRE\n0.02\nSparse Treed\n\n\n100 m\nMixTRE\n0.01\nMixed Treed\n\n\n100 m\nXWAT\n0.01\nTurbid Water\n\n\n1 km\nTrFEN\n0.21\nTreed Fen\n\n\n1 km\nTrBOG\n0.21\nTreed Bog\n\n\n1 km\nOBOG\n0.19\nOpen Bog\n\n\n1 km\nConSWA\n0.17\nConiferous Swamp\n\n\n1 km\nOFEN\n0.07\nOpen Fen\n\n\n1 km\nConTRE\n0.06\nConiferous Treed\n\n\n1 km\nWAT\n0.04\nClear Open Water\n\n\n1 km\nMixTRE\n0.01\nMixed Treed\n\n\n1 km\nXWAT\n0.01\nTurbid Water\n\n\n1 km\nThSWA\n0.01\nThicket Swamp\n\n\n1 km\nSpTRE\n0.00\nSparse Treed\n\n\n1 km\nNSWood\n0.00\nDisturbance - Non and Sparse Woody\n\n\n1 km\nDecTRE\n0.00\nDeciduous Treed\n\n\n1 km\nTrOrSHr\n0.00\nDisturbance - Treed and/or Shrub\n\n\n1 km\nBED\n0.00\nBedrock\n\n\n\n\n\nWe also want to do the same operation for the ecodistrict to allow for comparison. We don’t need to use exact extraction, insteadt the crop and mask each raster. This operation is costly so we write out the rasters and load them again (see unrendered code).\n\necodistrict_16 <- sf::st_transform(ecodistrict, sf::st_crs(lu_16))\necodistrict_17 <- sf::st_transform(ecodistrict, sf::st_crs(lu_17))\n\nlu_16_crop <- raster::crop(lu_16, ecodistrict_16)\nlu_16_crop_mask <- raster::mask(lu_16_crop, ecodistrict_16)\n\nlu_17_crop <- raster::crop(lu_17, ecodistrict_17)\nlu_17_crop_mask <- raster::mask(lu_17_crop, ecodistrict_17)\n\nWe can then get the frequencies of values. This operation is also costly so we write out the objects and load them again (see unrendered code).\n\nlu_16_freq <- raster::freq(lu_16_crop_mask)\nlu_17_freq <- raster::freq(lu_17_crop_mask)\n\nWe combine the results of both UTMs.\n\nres_ecodistrict <- rbind(lu_16_freq, lu_17_freq) |>\n as.data.frame() |> \n dplyr::group_by(value) |> \n dplyr::summarise(count = sum(count)) |> \n dplyr::ungroup() |> \n dplyr::filter(!is.na(value)) |> \n dplyr::mutate(prop = count/sum(count)) |> \n dplyr::mutate(value = as.factor(value)) |>\n dplyr::left_join(lu_dat, by = c(\"value\" = \"cats\")) |>\n dplyr::filter(!is.na(label)) |> \n dplyr::select(category_code, prop, label) |> \n dplyr::mutate(scale = \"Ecodistrict\") |> \n dplyr::relocate(scale) |>\n dplyr::arrange(scale, dplyr::desc(prop))\n\nknitr::kable(res_ecodistrict)\n\n\n\n\nscale\ncategory_code\nprop\nlabel\n\n\n\n\nEcodistrict\nTrFEN\n0.25\nTreed Fen\n\n\nEcodistrict\nOBOG\n0.20\nOpen Bog\n\n\nEcodistrict\nTrBOG\n0.19\nTreed Bog\n\n\nEcodistrict\nConSWA\n0.12\nConiferous Swamp\n\n\nEcodistrict\nOFEN\n0.09\nOpen Fen\n\n\nEcodistrict\nWAT\n0.07\nClear Open Water\n\n\nEcodistrict\nConTRE\n0.04\nConiferous Treed\n\n\nEcodistrict\nNSWood\n0.01\nDisturbance - Non and Sparse Woody\n\n\nEcodistrict\nTrOrSHr\n0.01\nDisturbance - Treed and/or Shrub\n\n\nEcodistrict\nMixTRE\n0.01\nMixed Treed\n\n\nEcodistrict\nThSWA\n0.01\nThicket Swamp\n\n\nEcodistrict\nSpTRE\n0.00\nSparse Treed\n\n\nEcodistrict\nXWAT\n0.00\nTurbid Water\n\n\nEcodistrict\nDecTRE\n0.00\nDeciduous Treed\n\n\nEcodistrict\nMIN\n0.00\nSand/Gravel/Mine Tailings\n\n\nEcodistrict\nFrMAR\n0.00\nFreshwater Marsh\n\n\nEcodistrict\nBED\n0.00\nBedrock\n\n\nEcodistrict\nURB\n0.00\nCommunity/Infrastructure\n\n\nEcodistrict\nDecSWA\n0.00\nDeciduous Swamp\n\n\nEcodistrict\nInMAR\n0.00\nIntertidal Marsh\n\n\n\n\n\nAnd then combine the results between scales and utm.\n\nres <- rbind(res_points, res_ecodistrict) |> \n tidyr::complete(scale, label) |> \n tidyr::replace_na(list(prop = 0)) |> \n dplyr::mutate(label = forcats::fct_reorder(label, prop))", "crumbs": [ "Site Analysis", "Land Cover Analysis" @@ -38,5 +38,27 @@ "title": "About", "section": "", "text": "About this site\n\n1 + 1\n\n[1] 2" + }, + { + "objectID": "index.html#plotting-spatial-data", + "href": "index.html#plotting-spatial-data", + "title": "Land Cover Analysis", + "section": "Plotting spatial data", + "text": "Plotting spatial data\nIt is always a good idea to try and plot spatial data before any processing.\n\nggplot() +\n geom_sf(data = ecodistrict) +\n geom_sf(data = sf::st_transform(sites_possible, \n sf::st_crs(ecodistrict))) +\n theme_bw()\n\n\n\n\n\n\n\n\nPlotting the land cover data is difficult because it is provided is two different UTMs.", + "crumbs": [ + "Site Analysis", + "Land Cover Analysis" + ] + }, + { + "objectID": "index.html#results", + "href": "index.html#results", + "title": "Land Cover Analysis", + "section": "Results", + "text": "Results\nWe can plot the results with “dodged” ggplot2 barplots.\n\nmy_pal <- c('#c7e9b4','#7fcdbb','#41b6c4','#1d91c0','#225ea8','#0c2c84')\n\nggplot(res) +\n geom_bar(aes(x = label, y = prop, fill = scale, colour = scale), \n alpha = 0.8,\n stat = \"identity\",\n position = \"dodge\") +\n theme_bw() +\n theme(axis.text.x = element_text(angle = 20, vjust = 1, hjust = 1)) +\n labs(x = \"Land Use Class\", y = \"Sum of Proportions\",\n fill = \"Scale\", colour = \"Scale\") +\n scale_fill_manual(values = my_pal) +\n scale_color_manual(values = my_pal)\n\n\n\n\n\n\n\n\nRemoving the land use classes than are not present around sites, we get a slightly easier graph to read.\n\nonly_at_sites <- res |> \n dplyr::filter(prop > 0) |> \n dplyr::filter(scale != \"Ecodistrict\") |> \n dplyr::pull(label)\n\nres_filt <- res |> \n dplyr::filter(label %in% only_at_sites)\n\nggplot(res_filt) +\n geom_bar(aes(x = label, y = prop, fill = scale, colour = scale), \n alpha = 0.8,\n stat = \"identity\",\n position = \"dodge\") +\n theme_bw() +\n theme(axis.text.x = element_text(angle = 20, vjust = 1, hjust = 1)) +\n labs(x = \"Land Use Class\", y = \"Sum of Proportions\",\n fill = \"Scale\", colour = \"Scale\") +\n scale_fill_manual(values = my_pal) +\n scale_color_manual(values = my_pal)", + "crumbs": [ + "Site Analysis", + "Land Cover Analysis" + ] } ] \ No newline at end of file diff --git a/index.qmd b/index.qmd index d2f42f8..586331e 100644 --- a/index.qmd +++ b/index.qmd @@ -183,6 +183,8 @@ And then combine the results between scales and utm. ```{r} res <- rbind(res_points, res_ecodistrict) |> + tidyr::complete(scale, label) |> + tidyr::replace_na(list(prop = 0)) |> dplyr::mutate(label = forcats::fct_reorder(label, prop)) ``` @@ -206,3 +208,27 @@ ggplot(res) + scale_fill_manual(values = my_pal) + scale_color_manual(values = my_pal) ``` + +Removing the land use classes than are not present around sites, we get a slightly easier graph to read. + +```{r fig.width=13, fig.height=8} +only_at_sites <- res |> + dplyr::filter(prop > 0) |> + dplyr::filter(scale != "Ecodistrict") |> + dplyr::pull(label) + +res_filt <- res |> + dplyr::filter(label %in% only_at_sites) + +ggplot(res_filt) + + geom_bar(aes(x = label, y = prop, fill = scale, colour = scale), + alpha = 0.8, + stat = "identity", + position = "dodge") + + theme_bw() + + theme(axis.text.x = element_text(angle = 20, vjust = 1, hjust = 1)) + + labs(x = "Land Use Class", y = "Sum of Proportions", + fill = "Scale", colour = "Scale") + + scale_fill_manual(values = my_pal) + + scale_color_manual(values = my_pal) +```