terra objects, RATs, WCS tutorial

[dylanbeaudette]

Found some recent incompatibilities with terra (terra_1.6-7), rasterVis (rasterVis_0.51.2). Same errors with latest available R and RStudio on gov machines.

Using the latest terra (R-Universe) + rasterVis (GH), the following issues warnings and produces no output. This worked as of 2022-05-10.

library(soilDB)
library(rasterVis)
library(sf)
library(terra)
library(viridis)
# 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')

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

Something is wrong / different in terra, as the following no longer works as expected.

a <- st_bbox(
  c(xmin = -96.7696, xmax = -96.6477, ymin = 36.5477, ymax = 36.6139), 
  crs = st_crs(4326)
)

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

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

# extract RAT for thematic mapping
rat <- cats(mu)[[1]]

# variables of interest
vars <- c("dbthirdbar_r", "awc_r", "ph1to1h2o_r")

# get / aggregate specific horizon-level properties from SDA
# be sure to see the manual page for this function
p <-  get_SDA_property(property = vars,
                       method = "Dominant Component (Numeric)", 
                       mukeys = as.integer(rat$mukey),
                       top_depth = 0,
                       bottom_depth = 25)


# convert areasymbol into a factor easy plotting later
p$areasymbol <- factor(p$areasymbol)

# merge aggregate soil data into RAT
rat <- merge(rat, p, by.x = 'mukey', by.y = 'mukey', sort = FALSE, all.x = TRUE)
levels(mu) <- rat

## ERROR HERE

# select 'areasymbol'
activeCat(mu) <- 'areasymbol'
plot(mu, axes = FALSE, maxcell = 1e5, col = viridis(2))
box()

Oddly enough, it is still possible to convert ID + RAT → stand alone grids. The following works as expected.

# convert mukey grid + RAT -> stack of numerical grids
mu.stack <- catalyze(mu)

# keep only properties / remove IDs
mu.stack <- mu.stack[[vars]]

I suspect that this has something to do with changes to terra::levels(), which no longer includes the full RAT...?

[brownag]

I suspect that this has something to do with changes to terra::levels(), which no longer includes the full RAT...?

Negative. It is an issue of recent (1.5-34+) changes to terra::subst(), which is used internally by levelplot(). We rightly use cats() for the full RAT, and this was updated accordingly in rasterVis 0.51.2.

I made an issue in terra to resolve--and it looks like Robert, fastest maintainer in the West, is already hard at work on it!

[dylanbeaudette]

Thanks for checking on this.

Something changed in the behavior of terra::levels() and / or RAT management. For example, I noticed that a column name in an automatically generated RAT change from labels → label with the new version of terra.

Is it wrong to expect levels(x) to return the full RAT? I've noticed that sometimes it does (single attribute RAT), and other times it does not (10+ attributes).

[brownag]

Something changed in the behavior of terra::levels() and / or RAT management. For example, I noticed that a column name in an automatically generated RAT change from labels → label with the new version of terra.

Yes, there have been some changes to levels(), related I think to the way you used to be able to set them with an atomic vector, but I am pretty sure that is not the cause of this issue.

Is it wrong to expect levels(x) to return the full RAT? I've noticed that sometimes it does (single attribute RAT), and other times it does not (10+ attributes).

Both cats and levels return a list of data.frame, one for each layer. If I am not mistaken, levels() does not return the "ID" column in the case of our mukey.wcs() result--but all of the categories should be included.

[brownag]

Both of these problems (unrelated to each other) should now be fixed w/ rasterVis 0.51.4 (after https://github.com/oscarperpinan/rastervis/commit/6eca4c891635e6ad33586acb390edca9c2d6da32) + terra 1.6-27 (after https://github.com/rspatial/terra/commit/27d4e60dde80404778878d5d296bc5f597b49075)

[dylanbeaudette]

Nice detective work, and thanks for submitting issues to both packages. I'm testing the WCS tutorial now and there are still a couple of lines with errors (rasterVis).

[brownag]

I'm testing the WCS tutorial now and there are still a couple of lines with errors (rasterVis).

Ahh yeah, he is still using subst() the old way. I had alluded to this issue but didnt make a specific issue since I wanted it to stabilize in terra first. There are a few options to fix that, described here: https://github.com/oscarperpinan/rastervis/issues/92

[dylanbeaudette]

Great, thanks for the quick issue! At least we can now fall back on the terra plot methods. I've been wanting to extend / update this tutorial for a couple of months now and the plotting bugs kept me away from it.

[brownag]

An update on this issue: all of the {rasterVis} and {terra} bugs have been fixed. Requires terra >=1.6-28 and rasterVis >= 0.51.4, neither of which have been released to CRAN yet.

Since none of this is an issue with soilDB I think this issue could be safely closed, but I plan to leave it open until terra updates on CRAN.

Folks can download the development version of terra from r-universe and should be able to visualize the web coverage service grids effectively. On CCE machines this usually requires the following set of commands to avoid SSL errors:

# update URL as needed for latest R 4.2.x windows binary from: https://rspatial.r-universe.dev/ui#package:terra
curl::curl_download("https://rspatial.r-universe.dev/bin/windows/contrib/4.2/terra_1.6-33.zip", "terra.zip")
install.packages("terra.zip", repos = NULL)
unlink('terra.zip')

[dylanbeaudette]

Thanks. Did some testing last week and over the weekend. Much better after the updates to both packages. There are still some minor errors with rasterVis when layer names contain spaces. I'll bring that up with the author. Otherwise, I've been converting many documents over to terra::plot() when possible. The main reason for using rasterVis is the cleanly merged legend and paneling.

[brownag]

As a heads up last week I did encounter another edge case in the terra plot method for categorical rasters where values in the grid do not have a corresponding record in the levels. Posted it to terra today, so hopefully it will get fixed and included in next CRAN release. EDIT: fixed in terra >=1.6-40

I haven't seen it happening in the wild or any of our documentation materials but I encountered it when I was doing something a little sloppy with a categorical raster.

[dylanbeaudette]

Interesting--those kind of bugs really bother me. I'll keep an eye out. I finally posted the sometimes 0-indexed RAT issue: https://github.com/rspatial/terra/issues/880

[brownag]

I think the terra-based code in soilDB has been stable for some time now, though we have had issues with the visualization methods and examples outside the package.

I would like to create a SSURGO web coverage service based vignette, which will allow us to be regularly assured that everything is working as expected, and then we can close this issue. Building the vignette on the R-devel CI brings in terra from r-universe which should help us catch issues before they make it to CRAN.

[dylanbeaudette]

I think that all of the related topics have been addressed. All of my latest use of terra/categorical grids has been great.

Good idea on the vignette. Feel free to use / adjust / etc. any of the WCS tutorial content. That is due for a re-org. and maybe splint into smaller documents soon.

[brownag]

Good idea on the vignette. Feel free to use / adjust / etc. any of the WCS tutorial content. That is due for a re-org. and maybe splint into smaller documents soon.

I want to use as much of it as I can without dramatically altering dependencies. I was looking over the tutorial the other day and I think there are probably 3 vignettes worth of info in there, so parsing it up will take some thought, and probably there would be overlap/good opportunity for cross-linking topics

[dylanbeaudette]

I think your estimate is about right. Happy to help.