Introduction to the rasterbc package

Loading/merging blocks

To display the elevation data for the entire district, we need to combine the three blocks downloaded earlier. This can be done using opendata_bc, which loads all required blocks, merges them into a single layer, crops and masks as needed, and then loads into memory the returned SpatRaster object:

example.tif = opendata_bc(example.sf, collection='dem', varname='dem')
#> creating mosaic of 3 block(s)
#> clipping layer...masking layer...done
print(example.tif)
#> class       : SpatRaster 
#> dimensions  : 622, 893, 1  (nrow, ncol, nlyr)
#> resolution  : 100, 100  (x, y)
#> extent      : 1426588, 1515888, 539188, 601388  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs 
#> source      : memory 
#> name        :     dem 
#> min value   :  340.43 
#> max value   : 2153.29
plot(example.tif, main=paste(example.name, ': elevation'))
plot(st_geometry(example.sf), add=TRUE)
plot(st_geometry(blocks), add=TRUE, border='red')

If you’re opening the downloaded data right away, you can skip the getdata_bc call. The argument dl (by default TRUE) in opendata_bc will cause it to call getdata_bc automatically to download any missing blocks.

Any simple features object of class sf or sfc can be used for the argument geo, provided its geometry intersects with the provincial boundary of BC. The intended usage is for the user to delineate their region of interest as a (MULTI)POLYGON object (here, example.sf is a MULTIPOLYGON). Geometries of other classes (such SpatialPolygons, as defined by sp; or data frames containing coordinates of vertices) can often be coerced to sf using a command like sf::st_as_sf(other_geometry_class_object). Alternatively, users can directly download individual blocks by specifying their NTS/SNRC codes, eg. the next chunk plots the DEM again, using the codes:

example.codes = findblocks_bc(example.sf)
example.tif = opendata_bc(example.codes, collection='dem', varname='dem')
#> creating mosaic of 3 block(s)
#> done
plot(example.tif, main=paste('NTS/SNRC mapsheets ', paste(example.codes, collapse=', '), ': elevation'))
plot(st_geometry(blocks), add=TRUE, border='red')
plot(st_geometry(example.sf), add=TRUE)
text(st_coordinates(st_centroid(st_geometry(blocks))), labels=blocks$NTS_SNRC, cex=0.5)

File management

If you forget which files have been downloaded, you can either check the directory data.dir using your file browser (subfolder ‘/dem/blocks’, in this case), or use listfiles_bc to get a logical vector indicating which files are currently found in your local storage directory:

is.downloaded = listdata_bc(collection='dem', varname='dem', simple=TRUE)
paste('downloaded: ', sum(is.downloaded), '/', length(is.downloaded)) |> print()
#> [1] "downloaded:  3 / 89"

This shows that of the 89 blocks for the variable name ‘dem’ (in the collection ‘dem’), we have downloaded three so far. Notice the return value of listdata_bc is a named vector, with names indicating the destination filenames and paths. This shows where they will be written by getdata_bc. All filenames are either of the form ‘varname_mapsheet.tif’ (as in this example) or else varname_year_mapsheet.tif (for time-series data).

By default, the listdata_bc function prints a list of all available layers. eg. in the ‘dem’ collection we also have ‘aspect’ and ‘slope’:

listdata_bc(collection='dem', verbose=2)
#>        year                        description                                  unit tiles
#> dem      NA              digital elevation map              (metres above sea level)  3/89
#> slope    NA derived from digital elevation map            (degrees above horizontal)  3/89
#> aspect   NA derived from digital elevation map (degrees counterclockwise from north)  0/89

Notice the ‘dem’ blocks that were downloaded earlier. We merged these blocks in the opendata_bc function call that created example.tif. Currently, this layer resides in memory and can be accessed via the R object example.tif. To save a copy, one can use the terra::writeRaster function:

dem.path = file.path(getOption('rasterbc.data.dir'), 'dem', 'example_dem.tif')
terra::writeRaster(example.tif, dem.path, overwrite=TRUE)

getdata_bc writes all of its data inside a ‘blocks’ subdirectory (in this case ‘/dem/blocks’), and the subfolder of the data directory corresponding to the collection (in this case ‘/dem’) is, by default, left empty. So it is a good place to store and organize such derivative files, where they can be loaded more quickly (in future), eg.

terra::rast(dem.path)
#> class       : SpatRaster 
#> dimensions  : 2459, 2983, 1  (nrow, ncol, nlyr)
#> resolution  : 100, 100  (x, y)
#> extent      : 1286588, 1584888, 450888, 696788  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs 
#> source      : example_dem.tif 
#> name        :      dem 
#> min value   : 7.653875 
#> max value   : 2918.951

If you’re finished with rasterbc and want to remove all of the stored data, or if you simply want to free up space, the entire data directory or any of its contents can be deleted using your file browser. This will not break the rasterbc installation. However, all downloaded data will be erased and you will need to run datadir_bc again before using the other package functions.

Integer codes

Note that the bgcz collection data are factors, which are then encoded in the geotiff files as integer codes. opendata_bc returns a plot-ready categorical raster in this case, with coded levels replaced by factor names. The complete lookup tables are also stored in the lazy loaded list object metadata_bc.

lookup.list = rasterbc::metadata_bc$bgcz$metadata$coding
print(lookup.list$zone)
#>  [1] "BAFA"
#>  [2] "BG"  
#>  [3] "BWBS"
#>  [4] "CDF" 
#>  [5] "CMA" 
#>  [6] "CWH" 
#>  [7] "ESSF"
#>  [8] "ICH" 
#>  [9] "IDF" 
#> [10] "IMA" 
#> [11] "MH"  
#> [12] "MS"  
#> [13] "PP"  
#> [14] "SBPS"
#> [15] "SBS" 
#> [16] "SWB" 
#> [17] NA

For example we have zone, 1 = Boreal Altai Fescue Alpine (BAFA), 2 = Bunchgrass, etc. See the documentation for the bgcz source script for links to a complete description of all codes.

The code below plots this data for the example region:

# open the biogeoclimatic zone raster
bgcz.raster = opendata_bc(geo=example.sf, collection='bgcz', varname='zone', quiet=TRUE)

# set up a colour palette and plot with legend defined manually
plot(bgcz.raster, col=rainbow(5), main='Biogeoclimatic zones')
plot(st_geometry(example.sf), add=TRUE)