The rSDP
package provides functions for finding,
accessing, querying, wrangling, and sub-setting a set of curated raster
data products. This vignette provides some examples of how to perform
the first two steps (finding and accessing), along with some example
workflows for exploring and visualizing the data.
To reproduce the examples in this article, you will need to install
the latest versions of the rSDP
,terra
and
leaflet
packages. For the web maps, you will need to
install the development version of the leaflet package from GitHub.
remotes::install_github("rmbl-sdp/rSDP")
remotes::install_github("rstudio/leaflet")
install.packages(c("terra","tidyterra","sf"),type="binary") ## `binary` install prevents dependency issues on Mac and Windows
Loading the libraries makes the functions available.
Finding datasets
The RMBL Spatial Data Platform provides more than 100 spatial data products, so the first challenge is finding data that is relevant to your work. There are currently two ways to find data:
- Searching and browsing the web-based Data Catalog on RMBL’s website.
- Searching and filtering datasets using the built-in catalog in the
rSDP
package.
Searching the SDP Catalog from R
You can get a data frame with the catalog information for all datasets with
cat <- sdp_get_catalog()
head(cat[,1:4])
#> CatalogID Release Type Product
#> 1 R1D001 Release1 Hydro Large Stream Flowlines (Multi-flow Direction)
#> 2 R1D002 Release1 Hydro Large Stream Flowlines (Single-flow Direction)
#> 3 R1D003 Release1 Hydro Large Watersheds (Multi-flow Direction)
#> 4 R1D004 Release1 Hydro Multi-direction Flow Accumulation
#> 5 R1D005 Release1 Hydro Small Watersheds (Single Flow Direction)
#> 6 R1D006 Release1 Hydro Single Direction Flow Accumulation
If you are working in Rstudio, the best way to explore the catalog is by opening it up in the Viewer pane.
View(cat)
This allows you to use the built-in filtering and searching facilities of the RStudio Viewer.
If you want to filter the catalog programmatically, you can specify
arguments to sdp_get_catalog()
that return matching rows.
For example, if you wan to return only the products that have to do with
vegetation, you can specify:
veg_cat <- sdp_get_catalog(types="Vegetation")
head(veg_cat[,1:4])
#> CatalogID Release Type Product
#> 25 R1D025 Release1 Vegetation Basic Landcover
#> 26 R1D026 Release1 Vegetation Leaf-on Digital Surface Model
#> 27 R1D027 Release1 Vegetation Vegetation Canopy Height
#> 40 R2D013 Release2 Vegetation High-res Canopy Height
#> 41 R2D014 Release2 Vegetation Lidar Vegetation Metrics
#> 54 R3D013 Release3 Vegetation Understory Cover
For more advanced filtering, you can use regular expressions to match particular products.
snow_cat <- cat[grepl("Snowpack Persist",cat$Product),]
head(snow_cat[,1:4])
#> CatalogID Release Type
#> 75 R4D001 Release4 Snow
#> 135 R4D061 Release4 Snow
#> 136 R4D062 Release4 Snow
#> 137 R4D063 Release4 Snow
#> Product
#> 75 Snowpack Persistence Day of Year Yearly Timeseries
#> 135 Snowpack Persistence Day of Year Mean (1993 - 2022)
#> 136 Snowpack Persistence Day of Year Standard Deviation (1993-2022)
#> 137 Snowpack Persistence Uncertainty Yearly Timeseries
The “CatalogID” field provides a concise, unique identifier for each data product, and is the preferred way to specify an individual product. Here, the CatalogID “R4D001” represents the data product representing annual time-series of snow persistence across the Upper Gunnison domain.
Getting detailed metadata
To get detailed metadata for each item (beyond the basic information
that is provided in the Catalog), we can use the
sdp_get_metadata()
function.
snow_meta <- sdp_get_metadata("R4D001")
print(snow_meta$qgis$extent$spatial)
#> list()
#> attr(,"minx")
#> [1] "305073"
#> attr(,"miny")
#> [1] "4256064"
#> attr(,"maxz")
#> [1] "0"
#> attr(,"crs")
#> [1] "EPSG:32613"
#> attr(,"maxx")
#> [1] "388098"
#> attr(,"dimensions")
#> [1] "2"
#> attr(,"maxy")
#> [1] "4328667"
#> attr(,"minz")
#> [1] "0"
Visualizing the SDP domains
We currently provide SDP data products in three different spatial domains. To visualize the boundaries of these domains, we can construct a web map that shows them. First we need to connect to these datasets stored in the cloud.
GT_bound <- sf::st_read("https://rmbl-sdp.s3.us-east-2.amazonaws.com/data_products/supplemental/GT_region_vect_1m.geojson")
#> Reading layer `GT_region_vect_1m' from data source
#> `https://rmbl-sdp.s3.us-east-2.amazonaws.com/data_products/supplemental/GT_region_vect_1m.geojson'
#> using driver `GeoJSON'
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 327105 ymin: 4313359 xmax: 328312 ymax: 4314793
#> Projected CRS: WGS 84 / UTM zone 13N
UER_bound <- sf::st_read("https://rmbl-sdp.s3.us-east-2.amazonaws.com/data_products/supplemental/UER_region_vect_1m.geojson")
#> Reading layer `UER_region_vect_1m' from data source
#> `https://rmbl-sdp.s3.us-east-2.amazonaws.com/data_products/supplemental/UER_region_vect_1m.geojson'
#> using driver `GeoJSON'
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 315965 ymin: 4298054 xmax: 337062 ymax: 4322641
#> Projected CRS: WGS 84 / UTM zone 13N
UG_bound <- sf::st_read("https://rmbl-sdp.s3.us-east-2.amazonaws.com/data_products/supplemental/UG_region_vect_1m.geojson")
#> Reading layer `UG_region_vect_1m' from data source
#> `https://rmbl-sdp.s3.us-east-2.amazonaws.com/data_products/supplemental/UG_region_vect_1m.geojson'
#> using driver `GeoJSON'
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 307895 ymin: 4258017 xmax: 385679 ymax: 4326087
#> Projected CRS: WGS 84 / UTM zone 13N
domain_bounds <- rbind(UG_bound,UER_bound,GT_bound)
domain_bounds$Domain <- c("Upper Gunnison (UG)","Upper East River (UER)","Gothic Townsite (GT)")
Once we’ve got them in, we can use the terra::plet()
function to plot the domains on a web map.