Use terra/sf to replace raster/sp/rgdal

[brownag]

Preliminary investigations into replacing rgdal, sp and raster with terra and sf internally

Currently SDA_spatialQuery() will return Spatial* object if you use Spatial as input, and sf otherwise. Other functions that previously returned RasterLayer now return SpatRaster. I like the idea of returning the type the user provided, but most spatial functions do not (always) take spatial objects as input.

• Considering an option that can be used to control preferred output types

  • Functions like seriesExtent, taxaExtent, ISSR800.wcs/mukey.wcs, fetchSDA_spatial would need either an argument or a package option to set the preferred raster and vector types (raster options: "terra", "raster"; vector options: "sf", "terra", "sp")
    • Need raster back in Suggests to call raster::raster(<SpatRaster>)
    • Use sf to convert sf to sp sf::as_Spatial(<sf>)
    • Currently vector data default to sf, but could use terra SpatVector in all cases and provide it as an option

• can the CRS of the taxaExtent/seriesExtent rasters be set on the server side?

[dylanbeaudette]

Interesting, thanks for attempting this big switch. I'll do some testing after Christmas. Overall, I like the idea of swapping sp for sf in every place that it is used. I'm on the fence about terra / raster, mostly because of the major differences in raster attribute table management but also due to the fluid nature of terra at this time. Testing with a couple of different workflows will be helpful. I do like the option of returning as terra objects.

As for dependency management, I think that any package that uses raster or terra will automatically "pull in" the other as a hard dependency.

[brownag]

Cool thanks. I'll be interested to see what you find. From my testing things are working reasonably well, and existing workflows that "must" use raster should be able to convert the results to RasterLayer equivalents. I would be interested to see a specific example of issues with the RAT. I think terra has been in the process of stabilizing-- and it seems to have leveled out now considering that raster uses it internally (has it in imports). terra still suggests raster. I think this would be a significant change and warrants a minor version bump to 2.7.x if/when implemented

[brownag]

Re-running CI on this PR to see if switch to terra resolves issues with current CRAN status of raster and terra.

At time of writing this raster has not been updated on CRAN for a couple weeks, whereas terra was updated yesterday. More methods and functions have been duplicated / moved from raster to terra, resulting in some namespace issues. I believe these errors should go away when raster is updated on CRAN.

[brownag]

It has been reported that the current namespace conflict due to CRAN raster v3.5-11 and terra 1.5-13 interferes with loading the soilDB package. The CI has been breaking for a week--not sure what the timeline is on a raster release. Hopefully raster 3.5-12 will be sent to CRAN soon.

[brownag]

raster v3.5-12 was released on Saturday (2022-01-22), which should resolve warnings and issues with loading namespaces:

Result: WARN
    Found the following significant warnings:
     Warning: multiple methods tables found for ‘direction’
     Warning: multiple methods tables found for ‘gridDistance’ 

https://github.com/ncss-tech/soilDB/actions/runs/1731472001

[dylanbeaudette]

Doing some testing of taxaExtent and seriesExtent, using code + output from https://ncss-tech.github.io/AQP/soilDB/taxa-extent.html

R session info:

 [1] sf_1.0-5           spData_2.0.1       mapview_2.10.0     rgdal_1.5-28       sp_1.4-6           viridis_0.6.2     
 [7] viridisLite_0.4.0  rasterVis_0.51.1   lattice_0.20-45    terra_1.5-19       SoilTaxonomy_0.1.4 soilDB_2.6.13     

Made a couple of minor changes, and noticed the following:

• terra::aggregate() needs to have na.rm = TRUE explicitly set, this was not the case with raster::aggregate()
• terra::values(r) <- terra::values(r) is essential, not sure why

This isn't a big deal, but it means that examples / docs will need to be edited accordingly.

Will test WCS stuff next.

[dylanbeaudette]

Testing ISSR800 WCS, notes:

• setting / retrieving factor levels is markedly different in terra, many examples and existing code will break as-is
• rasterVis::levelplot() does not appear to work with categorical SpatRaster objects, conversion via raster::raster() works but will have to be added to a lot of examples / production code
• setting all factor levels vs. just those present (as in previous version) can cause problems when the number of categories is very large (e.g. series names) → I'd prefer to keep the original approach

It is important to keep all columns from the RAT, as there may be additional data like colors in there.

cols <- cats(texture_2550cm)[[1]]$hex
levelplot(raster::raster(texture_2550cm), att = 'names', margin = FALSE, col.regions = cols)

I'll add some more tests to the main branch, and then test here. Mainly, are raster dimensions as expected.

[dylanbeaudette]

I'm not sure how to "correctly" return the mukey WCS object. The categorical ISSR-800 grids seem fine, but those do not use terra::categories().

Stuck here:

# make a bounding box and assign a CRS (4326: GCS, WGS84)
a <- st_bbox(
  c(xmin = -114.16, xmax = -114.08, ymin = 47.65, ymax = 47.68), 
  crs = st_crs(4326)
)

# convert bbox to sf geometry
a <- st_as_sfc(a)

# fetch gSSURGO map unit keys at native resolution (~30m)
x <- mukey.wcs(aoi = a, db = 'gssurgo')

# error
plot(x)

# ok this works
plot(x, type = 'classes')

# this does not
levelplot(
  x, 
  att = 'gSSURGO map unit keys',
  main = 'gSSURGO map unit keys',
  sub = 'Albers Equal Area Projection',
  margin = FALSE, 
  colorkey = FALSE, 
  col.regions = viridis
)

[brownag]

Before (soilDB 2.6.14 -- main branch)

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(raster)
#> Loading required package: sp
library(soilDB)

# make a bounding box and assign a CRS (4326: GCS, WGS84)
a <- st_bbox(
  c(xmin = -114.16, xmax = -114.08, ymin = 47.65, ymax = 47.68), 
  crs = st_crs(4326)
)

# fetch gSSURGO map unit keys at native resolution (~30m)
x <- mukey.wcs(aoi = a, db = 'gssurgo')

x
#> class      : RasterLayer 
#> dimensions : 147, 219, 32193  (nrow, ncol, ncell)
#> resolution : 29.99083, 29.92466  (x, y)
#> extent     : -1365483, -1358915, 2869259, 2873658  (xmin, xmax, ymin, ymax)
#> crs        : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
#> source     : memory
#> names      : gSSURGO.map.unit.keys 
#> values     : 144983, 1716001  (min, max)
#> attributes :
#>             ID
#>  from:  144983
#>   to : 1716001

# raster plot
plot(x)

Using development rasterVis from github (>0.51.2)

library(rasterVis)
#> Loading required package: lattice
levelplot(x,
          att = 'ID',
          main = 'gSSURGO map unit keys',
          sub = 'Albers Equal Area Projection')

Latest on this PR (69d694c)

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.5.21
library(soilDB)
library(rasterVis)
#> Loading required package: lattice

# make a bounding box and assign a CRS (4326: GCS, WGS84)
a <- st_bbox(
  c(xmin = -114.16, xmax = -114.08, ymin = 47.65, ymax = 47.68), 
  crs = st_crs(4326)
)

# fetch gSSURGO map unit keys at native resolution (~30m)
x <- mukey.wcs(aoi = a, db = 'gssurgo')

x
#> class       : SpatRaster 
#> dimensions  : 147, 219, 1  (nrow, ncol, nlyr)
#> resolution  : 29.99541, 29.92466  (x, y)
#> extent      : -1365484, -1358915, 2869259, 2873658  (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83 / Conus Albers (EPSG:5070) 
#> source      : memory 
#> categories  : gSSURGO.map.unit.keys, ID 
#> name        : gSSURGO.map.unit.keys 
#> min value   :                144983 
#> max value   :               1716001

# terra plot
plot(x)

levelplot(x,
          att = 'ID',
          main = 'gSSURGO map unit keys',
          sub = 'Albers Equal Area Projection')

Full options:

levelplot(
  x, 
  att = 'gSSURGO.map.unit.keys',
  main = 'gSSURGO map unit keys',
  sub = 'Albers Equal Area Projection',
  margin = FALSE, 
  colorkey = FALSE, 
  col.regions = viridisLite::viridis
)

[brownag]

Last commit fixes plot(x, type="classes") and makes terra::values() output on resulting SpatRaster consistent with what it was with raster.

[dylanbeaudette]

Thanks for all of the recent additions. I'm 100% on board with the switch and ready to convert all of my related code over to using / expecting sf + terra.