Add crs argument to the three cded functions

[dfilatow]

add an argument to be able to specify the projection and have the output reprojected. By default it should be bcalbers and cded is not. Otherwise the output of the cded functions do not play nicely with the other bcmaps spatial objects.

ws <- bcdc_query_geodata("51f20b1a-ab75-42de-809d-bf415a0f9c62", crs = bcalb) %>%
  filter( WATERSHED_GROUP_CODE == "PARS") %>%
  collect() %>% st_as_sf()
dem.ws <- cded_raster(aoi = ws)
crs(dem.ws)
st_crs(ws)

[stephhazlitt]

Should this be the default behaviour (maybe w/o an argument, users can wrote code to re-project from bcalbers if needed)?

[ateucher]

We considered this, but decided that it was very computationally expensive and so left it up to the user

library(bcdata)
#> 
#> Attaching package: 'bcdata'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(bcmaps)
library(raster)

ws <- bcdc_query_geodata("51f20b1a-ab75-42de-809d-bf415a0f9c62") %>%
  filter( WATERSHED_GROUP_CODE == "PARS") %>%
  collect() %>% st_as_sf()

dem.ws <- cded_raster(aoi = ws)
#> checking your existing tiles for mapsheet 93i are up to date
#> checking your existing tiles for mapsheet 93o are up to date
#> checking your existing tiles for mapsheet 93j are up to date
#> checking your existing tiles for mapsheet 93p are up to date

system.time(
  newdem.ws <- projectRaster(dem.ws, crs = 3005)
)
#>    user  system elapsed 
#>  77.458  84.469 262.218

^{Created on 2021-02-18 by the reprex package (v1.0.0)}

[bevingtona]

It would be super fast if you don't touch the tif files, and just warp the vrt in the cded functions.

sf::gdal_utils(
  util = "warp", 
  source = "out.vrt", 
  destination = "out3005.vrt", 
  options = c("-t_srs", "EPSG:3005",
              "-r", "near",
              "-of", "VRT"))

[ateucher]

So would the actual warping then happen if you pulled it into R as a regular raster/stars object?

[bevingtona]

library(bcdata)
#> Warning: package 'bcdata' was built under R version 4.0.3
#> 
#> Attaching package: 'bcdata'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(bcmaps)
#> Warning: package 'bcmaps' was built under R version 4.0.3
#> Loading required package: sf
#> Warning: package 'sf' was built under R version 4.0.3
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)
#> Warning: package 'raster' was built under R version 4.0.3
#> Loading required package: sp
#> Warning: package 'sp' was built under R version 4.0.3
#> 
#> Attaching package: 'raster'
#> The following object is masked from 'package:bcdata':
#> 
#>     select

ws <- bcdc_query_geodata("51f20b1a-ab75-42de-809d-bf415a0f9c62") %>%
  filter( WATERSHED_GROUP_CODE == "PARS") %>%
  collect() %>% st_as_sf()
#> Note: method with signature 'DBIConnection#character' chosen for function 'dbQuoteIdentifier',
#>  target signature 'wfsConnection#ident'.
#>  "wfsConnection#ANY" would also be valid
#> Note: method with signature 'DBIConnection#character' chosen for function 'dbQuoteIdentifier',
#>  target signature 'wfsConnection#character'.
#>  "wfsConnection#ANY" would also be valid
#> Note: method with signature 'DBIConnection#character' chosen for function 'dbQuoteString',
#>  target signature 'wfsConnection#character'.
#>  "wfsConnection#ANY" would also be valid

dem.ws <- cded_raster(aoi = ws, dest_vrt = "temp.vrt")
#> checking your existing tiles for mapsheet 93i are up to date
#> checking your existing tiles for mapsheet 93o are up to date
#> checking your existing tiles for mapsheet 93j are up to date
#> checking your existing tiles for mapsheet 93p are up to date

sf::gdal_utils(
  util = "warp", 
  source = "temp.vrt", 
  destination = "out3005.vrt", 
  options = c("-t_srs", "EPSG:3005",
              "-r", "near",
              "-of", "VRT"))

raster::raster("out3005.vrt") %>% st_crs() == st_crs(ws)
#> [1] TRUE

^{Created on 2021-02-18 by the reprex package (v0.3.0)}

[bevingtona]

This is a work arround, but pretty sure this would be easy to add into the cded functions... so it is only the VRT that is being warped

[boshek]

One thing to consider is that it would be a breaking change. One other option would be to leverage transform_bc_albers and provide stars and raster methods for that function. Then this would just work:

cded_stars(aoi) %>%
   transform_bc_albers()

[ateucher]

I guess the other question is the choice of interpolation method... is there one that's best for dem data? Is the best choice different for different scenarios?

[ateucher]

I.e., @bevingtona it looks like you've used nearest-neighbour, but I would have expected cubic or bilinear would make sense?

[bevingtona]

busted @ateucher ! yes, bilinear would be best.