Reproducible workflow using R
The Cumberland Plain Woodland in the Sydney Basin Bioregion is listed as a critically endangered ecological community according to the New South Wales Threatened Species Conservation Act (NSW Scientific Committee 2009). This community is included within the critically endangered ecological community listed under the Environment Protection and Biodiversity Conservation Act as “Cumberland Plain Shale Woodlands and Shale-Gravel Transition Forest” (Department of the Environment 2024).
For this workflow we will load datasets from the Central Resource for Sharing and Enabling Environmental Data in NSW (SEED NSW) and other portals.
We identified the following relevant datasets for this workflow:
Data can be downloaded manually or programmatically into a data folder. An example for use in Linux or MacOSX is given in the next code snippets.
mkdir -p sandbox/
cd sandbox/
wget --continue https://datasets.seed.nsw.gov.au/dataset/b2220d40-fd3c-4404-ba7b-b8a867c7b8de/resource/99a64722-0803-48b1-8011-310a04cc9fff/download/vegetation-condition-benchmarks-cover-and-richness-raw-data-v1-2.xlsx
wget --continue https://www.environment.nsw.gov.au/-/media/OEH/Corporate-Site/Documents/BioNet/bionet-plant-community-type-data.xlsxmkdir -p sandbox/
cd sandbox/
wget --continue https://datasets.seed.nsw.gov.au/dataset/a55ea79e-967e-42cc-b6a7-2f2627a606be/resource/5426a635-fa31-48d3-bb0a-ed8ca4a0d60e/download/vegetationsouthcoastsciviv14e2230.zip
unzip -U vegetationsouthcoastsciviv14e2230.zip
wget --continue https://datasets.seed.nsw.gov.au/dataset/d872b07e-253d-47fa-8e27-fea0d2522b9d/resource/c2cb2276-7706-40ad-8922-63cca8520d37/download/ageonetworkzipscumberlandplainupdate3785.zip
unzip -U ageonetworkzipscumberlandplainupdate3785.zip This reproducible workflow is based on R and uses packages sf, leaflet and mapview for visualisation of the spatial data, plus additional packages for data import (readxl, httr) and data wrangling (dplyr, tidyr).
library(sf)
library(dplyr)
library(mapview)
library(leaflet)
library(readxl)
library(httr)
library(knitr)
library(tidyr)
library(units)We read the BioNet PCT data directly from the downloaded file and filter the PCTs in the Cumberland region and belonging to the Grassy Woodlands formation.
here::i_am("CPW-workflow-mapview.qmd")
data_dir <- here::here("sandbox")
BioNet_PCT <- read_excel(here::here(data_dir,
"bionet-plant-community-type-data.xlsx"), "PCT Data PQ")
BioNet_PCT |>
filter(vegetationFormation %in% "Grassy Woodlands", grepl("Cumberland",PCTName)) |>
select(PCTID,PCTName) |>
kable()| PCTID | PCTName |
|---|---|
| 830 | Cumberland moist shale woodland |
| 849 | Cumberland shale plains woodland |
| 850 | Cumberland shale hills woodland |
| 1395 | Cumberland shale - sandstone Ironbark forest |
| 3318 | Cumberland Moist Shale Woodland |
| 3319 | Cumberland Shale Hills Woodland |
| 3320 | Cumberland Shale Plains Woodland |
| 3321 | Cumberland Shale-Sandstone Ironbark Forest |
Cumberland Plain Woodland belongs to the Coastal Valley Grassy Woodlands vegetation class (Keith 2004) and includes: ‘Shale Hills Woodland’ (map unit 9) and ‘Shale Plains Woodland’ (map unit 10) of M. Tozer (2003); equivalent to the ‘Cumberland Shale Hills Woodland’ (map unit GW p28) and ‘Cumberland Shale Plains Woodland’ (map unit GW p29) of M. G. Tozer et al. (2010).
Now we can query spatial data from the NSW BioNet Flora Survey Data Collection for these PCT IDs using a REST API call.
# PCT codes to query:
PCT_ids <- BioNet_PCT |>
filter(vegetationFormation %in% "Grassy Woodlands", grepl("Cumberland",PCTName)) |>
pull(PCTID)
# Constructing url request for REST API:
mapserver_url <- paste("https://mapprod3.environment.nsw.gov.au/arcgis/rest/services",
"EDP/SystematicFloraSurvey/MapServer/0/query", sep = "/")
url <- parse_url(mapserver_url)
url$query <- list(where = sprintf("currentClassification IN (%s)",
paste(PCT_ids, collapse = ",")),
outFields = "*",
returnGeometry = "true",
f = "geoJSON")
request <- build_url(url)
# read spatial data using simple features package:
FloraSurveyQuery <- st_read(request)Reading layer `OGRGeoJSON' from data source
`https://mapprod3.environment.nsw.gov.au/arcgis/rest/services/EDP/SystematicFloraSurvey/MapServer/0/query?where=currentClassification%20IN%20%28830%2C849%2C850%2C1395%2C3318%2C3319%2C3320%2C3321%29&outFields=%2A&returnGeometry=true&f=geoJSON'
using driver `GeoJSON'
Simple feature collection with 768 features and 51 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 150.51 ymin: -34.38476 xmax: 151.0866 ymax: -33.24606
Geodetic CRS: WGS 84This query returns 768 features (point locations) and many columns, but we want to match this information with benchmark data.
The vegetation condition benchmark raw data is also available as a downloaded file, we read this file and join it with the downloaded Flora Survey data using the site id column, then we select the columns we want to display in the map. Only a subset of survey sites has cover and richness benchmark data.
benchmarks <- read_excel(here::here(data_dir,
"vegetation-condition-benchmarks-cover-and-richness-raw-data-v1-2.xlsx")) |>
filter(!is.na(Longitude)) |>
st_as_sf(coords = c("Longitude","Latitude"), crs = 4326) |>
filter(`Class IBRA` %in% "Coastal Valley Grassy Woodlands/Sydney Basin")
benchmarks_surveys <-
FloraSurveyQuery |>
inner_join(st_drop_geometry(benchmarks), by = c("siteID" = "Site no")) |>
select(`PCT id` = currentClassification,
`PCT name` = currentClassificationDesc,
elevationInMeters, DateTxt, FrCov:Rain12)
benchmarks_surveys |>
st_drop_geometry() |>
group_by(`PCT id`, `PCT name`) |>
summarise(`Nr. of sites` = n(), .groups = 'drop') |>
kable()| PCT id | PCT name | Nr. of sites |
|---|---|---|
| 3318 | Cumberland Moist Shale Woodland | 5 |
| 3319 | Cumberland Shale Hills Woodland | 52 |
| 3320 | Cumberland Shale Plains Woodland | 168 |
| 3321 | Cumberland Shale-Sandstone Ironbark Forest | 77 |
We will load the vegetation map from a web mapping service (WMS). The labels of the coloured polygons are only visible when zooming in to the map. The labels refers to class codes from M. G. Tozer et al. (2010). We can look up the codes in the raster attribute table that we have downloaded.
SciVi <- read_sf(here::here(data_dir,
"Vegetation_SouthCoast_SCIVI_v14_E_2230/Vegetation_SouthCoast_SCIVI_v14_E_2230/SouthCoast_SCIVI_v14_E_2230/Data/SouthCoast_SCIVI_v14_E_2230.tif.vat.dbf"))
SciVi |> filter(CLASS %in% c("p28", "p29")) |> select(CLASS, MAPUNIT) |> distinct()# A tibble: 2 × 2
CLASS MAPUNIT
<chr> <chr>
1 p28 Cumberland Shale Hills Woodland
2 p29 Cumberland Shale Plains WoodlandWe will also visualise a layer with the recent changes in the distribution of CPW. First we import the downloaded vector data into R and then filter the corresponding communities based on their description. To optimise visualisation in an on-line document, we group the original polygons by community description and status.
CPW_dist <- read_sf(here::here(data_dir, "cumberland_V2_2008_VISmap_3785.shp")) |>
st_transform(crs = '+proj=longlat +datum=WGS84 +no_defs') |>
filter(COMM_DESC %in%
c('14 - Moist Shale Woodland',
'103 - Shale/Gravel Transition Forest',
'9 - Shale Hills Woodland',
'10 - Shale Plains Woodland')) |>
st_make_valid()
CPW_changes <- CPW_dist |>
group_by(COMM_DESC, STATUS2007) |>
summarise(geometry = st_union(geometry),
.groups = "drop")We can summarise the information from the vector data into a table. Here we calculate the area (in hectars) under each category for each community.
options(knitr.kable.NA = '--')
CPW_changes |>
mutate(area=st_area(geometry) |> set_units('ha') |> round(2)) |>
st_drop_geometry() |>
pivot_wider(names_from=STATUS2007, values_from=area) |>
knitr::kable(digits = 2)| COMM_DESC | Cleared 1998 - 2007 | Cleared pre 1998 | Extant 2007 | Possibly cleared 1998 - 2007 |
|---|---|---|---|---|
| 10 - Shale Plains Woodland | 306.74 [ha] | 0.44 [ha] | 5695.60 [ha] | 45.78 [ha] |
| 103 - Shale/Gravel Transition Forest | 38.17 [ha] | 1.87 [ha] | 1639.08 [ha] | 16.02 [ha] |
| 14 - Moist Shale Woodland | 0.20 [ha] | – | 46.41 [ha] | 0.15 [ha] |
| 9 - Shale Hills Woodland | 89.11 [ha] | – | 2505.82 [ha] | 15.76 [ha] |
This datasets cover ca. 80% of the total distribution of Cumberland Plain Woodland (NSW Scientific Committee and Simpson 2008).
Now we can show these layers in a interactive map.
wms_src <- "https://mapprod3.environment.nsw.gov.au/arcgis/services/VIS/Vegetation_SouthCoast_SCIVI_V14_E_2230/MapServer/WMSServer?"
mapviewOptions(fgb=FALSE)
m <-
mapview(benchmarks_surveys,
zcol = 'PCT name',
layer.name = "Benchmark raw data v1.2", col.regions = sf.colors) +
mapview(CPW_changes, zcol = 'STATUS2007',
layer.name = "CPW distribution changes (Simpson 2008)")
m@map <- m@map |>
addWMSTiles(group = 'Vegetation map (Tozer 2010)',
wms_src,
layers = "0,1",
options = WMSTileOptions(format = "image/png", transparent = TRUE),
attribution = "") |>
mapview:::mapViewLayersControl(names = c('Vegetation map (Tozer 2010)'))
mR version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.4.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Australia/Sydney
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] units_0.8-5 tidyr_1.3.0 knitr_1.45 httr_1.4.7 readxl_1.4.3
[6] leaflet_2.2.1 mapview_2.11.2 dplyr_1.1.4 sf_1.0-14
loaded via a namespace (and not attached):
[1] xfun_0.41 raster_3.6-26 htmlwidgets_1.6.4
[4] lattice_0.22-5 leaflet.providers_2.0.0 vctrs_0.6.5
[7] tools_4.3.1 crosstalk_1.2.1 generics_0.1.3
[10] stats4_4.3.1 curl_5.2.0 tibble_3.2.1
[13] proxy_0.4-27 fansi_1.0.6 pkgconfig_2.0.3
[16] KernSmooth_2.23-22 satellite_1.0.4 uuid_1.1-1
[19] lifecycle_1.0.4 farver_2.1.1 compiler_4.3.1
[22] munsell_0.5.0 terra_1.7-55 codetools_0.2-19
[25] htmltools_0.5.7 class_7.3-22 yaml_2.3.8
[28] jquerylib_0.1.4 pillar_1.9.0 ellipsis_0.3.2
[31] classInt_0.4-10 wk_0.9.1 brew_1.0-8
[34] tidyselect_1.2.0 digest_0.6.33 purrr_1.0.2
[37] rprojroot_2.0.4 fastmap_1.1.1 grid_4.3.1
[40] here_1.0.1 colorspace_2.1-0 cli_3.6.2
[43] magrittr_2.0.3 base64enc_0.1-3 utf8_1.2.4
[46] leafem_0.2.3 e1071_1.7-14 withr_2.5.2
[49] scales_1.3.0 sp_2.1-2 rmarkdown_2.25
[52] cellranger_1.1.0 png_0.1-8 evaluate_0.23
[55] s2_1.1.5 rlang_1.1.2 Rcpp_1.0.11
[58] leafpop_0.1.0 glue_1.6.2 DBI_1.1.3
[61] svglite_2.1.3 rstudioapi_0.15.0 jsonlite_1.8.8
[64] R6_2.5.1 systemfonts_1.0.5 @online{ferrer-paris2024,
author = {Ferrer-Paris, José R.},
title = {Reproducible Workflow for Visualisation of Geospatial Data of
the {Cumberland} {Plain} {Woodland}},
date = {2024-05-16},
url = {https://jrfep.quarto.pub/cpw-workflow-mapview/},
langid = {en}
}