Add support to crs "local" (recently defined in terra)

[aloboa]

Please add support to crs "local"
rspatial/terra#797

Given:

r <- rast(nrows=4, ncols=4)
values(r) <- 1:ncell(r)
r <- c(r,r,r)

The following works as expected:

crs(r) <- "epsg:25831"
ggplot() + geom_spatraster_rgb(data = r)
ggplot() + geom_spatraster_rgb(data = stretch(r))
crs(r) <- ""
ggplot() + geom_spatraster_rgb(data = stretch(r))
ggplot() + geom_spatraster_rgb(data = stretch(r)) + coord_fixed(ratio = 1)
ggplot() + geom_spatraster_rgb(data = stretch(r))
ggplot() + geom_spatraster_rgb(data = stretch(r)) + coord_fixed(ratio = 1)

But:

 crs(r) <- "local"
 ggplot() + geom_spatraster_rgb(data = stretch(r))
Error in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  : 
  OGRCreateCoordinateTransformation(): transformation not available
In addition: Warning message:
In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 6: Cannot find coordinate operations from `ENGCRS["Cartesian (Meter)",EDATUM["Local Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]' to `ENGCRS["Cartesian (Meter)",EDATUM["Local Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]'

ggplot() + geom_spatraster_rgb(data = stretch(r)) + coord_fixed(ratio = 1)
Error in `reproject_raster_on_stat()`:
! geom_spatraster_*() on SpatRasters with crs must be used with coord_sf().
Run `rlang::last_error()` to see where the error occurred.

[dieghernan]

Thanks for pointing this out. I pressume this functionality has been added to the dev version, right?

library(terra)
#> Warning: package 'terra' was built under R version 4.1.3
#> terra 1.6.17

# Raster
f <- system.file("extdata/cyl_temp.tif", package = "tidyterra")
r <- rast(f)

crs(r) <- "local"
#> Warning: [crs<-] Cannot set raster SRS: empty srs

# Not working
r
#> class       : SpatRaster 
#> dimensions  : 89, 116, 3  (nrow, ncol, nlyr)
#> resolution  : 3856.617, 3856.617  (x, y)
#> extent      : 2893583, 3340950, 2019451, 2362690  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
#> source      : cyl_temp.tif 
#> names       :   tavg_04,   tavg_05,   tavg_06 
#> min values  :  0.565614,  4.294102,  8.817221 
#> max values  : 13.283829, 16.740898, 21.113781

# This works
crs(r) <- ""
r
#> class       : SpatRaster 
#> dimensions  : 89, 116, 3  (nrow, ncol, nlyr)
#> resolution  : 3856.617, 3856.617  (x, y)
#> extent      : 2893583, 3340950, 2019451, 2362690  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source      : cyl_temp.tif 
#> names       :   tavg_04,   tavg_05,   tavg_06 
#> min values  :  0.565614,  4.294102,  8.817221 
#> max values  : 13.283829, 16.740898, 21.113781

^{Created on 2022-11-08 with reprex v2.0.2}

Session info

sessioninfo::session_info()
#> - Session info ---------------------------------------------------------------
#>  setting  value
#>  version  R version 4.1.0 (2021-05-18)
#>  os       Windows 10 x64 (build 19044)
#>  system   x86_64, mingw32
#>  ui       RTerm
#>  language (EN)
#>  collate  Spanish_Spain.1252
#>  ctype    Spanish_Spain.1252
#>  tz       Europe/Paris
#>  date     2022-11-08
#>  pandoc   ***
#> 
#> - Packages -------------------------------------------------------------------
#>  package     * version date (UTC) lib source
#>  cli           3.4.1   2022-09-23 [1] CRAN (R 4.1.3)
#>  codetools     0.2-18  2020-11-04 [2] CRAN (R 4.1.0)
#>  digest        0.6.29  2021-12-01 [1] CRAN (R 4.1.3)
#>  evaluate      0.17    2022-10-07 [1] CRAN (R 4.1.3)
#>  fansi         1.0.3   2022-03-24 [1] CRAN (R 4.1.3)
#>  fastmap       1.1.0   2021-01-25 [2] CRAN (R 4.1.0)
#>  fs            1.5.2   2021-12-08 [1] CRAN (R 4.1.3)
#>  glue          1.6.2   2022-02-24 [1] CRAN (R 4.1.3)
#>  highr         0.9     2021-04-16 [2] CRAN (R 4.1.0)
#>  htmltools     0.5.3   2022-07-18 [1] CRAN (R 4.1.3)
#>  knitr         1.40    2022-08-24 [1] CRAN (R 4.1.3)
#>  lifecycle     1.0.3   2022-10-07 [1] CRAN (R 4.1.0)
#>  magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.1.3)
#>  pillar        1.8.1   2022-08-19 [1] CRAN (R 4.1.3)
#>  pkgconfig     2.0.3   2019-09-22 [2] CRAN (R 4.1.0)
#>  purrr         0.3.5   2022-10-06 [1] CRAN (R 4.1.3)
#>  R.cache       0.16.0  2022-07-21 [1] CRAN (R 4.1.3)
#>  R.methodsS3   1.8.2   2022-06-13 [1] CRAN (R 4.1.3)
#>  R.oo          1.25.0  2022-06-12 [1] CRAN (R 4.1.3)
#>  R.utils       2.12.0  2022-06-28 [1] CRAN (R 4.1.3)
#>  Rcpp          1.0.9   2022-07-08 [1] CRAN (R 4.1.3)
#>  reprex        2.0.2   2022-08-17 [1] CRAN (R 4.1.3)
#>  rlang         1.0.6   2022-09-24 [1] CRAN (R 4.1.3)
#>  rmarkdown     2.17    2022-10-07 [1] CRAN (R 4.1.3)
#>  rstudioapi    0.14    2022-08-22 [1] CRAN (R 4.1.3)
#>  sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.1.2)
#>  stringi       1.7.6   2021-11-29 [1] CRAN (R 4.1.2)
#>  stringr       1.4.1   2022-08-20 [1] CRAN (R 4.1.3)
#>  styler        1.7.0   2022-03-13 [1] CRAN (R 4.1.3)
#>  terra       * 1.6-17  2022-09-10 [1] CRAN (R 4.1.3)
#>  tibble        3.1.8   2022-07-22 [1] CRAN (R 4.1.3)
#>  utf8          1.2.2   2021-07-24 [1] CRAN (R 4.1.3)
#>  vctrs         0.5.0   2022-10-22 [1] CRAN (R 4.1.3)
#>  withr         2.5.0   2022-03-03 [1] CRAN (R 4.1.3)
#>  xfun          0.33    2022-09-12 [1] CRAN (R 4.1.3)
#>  yaml          2.3.5   2022-02-21 [1] CRAN (R 4.1.2)
#> 
#> 
#> ------------------------------------------------------------------------------

[aloboa]

Thank you for including it.

I pressume this functionality has been added to the dev version, right?

According to rspatial/terra#797 (comment)
it should have made it to cran by the end of September,

[kadyb]

The latest version on CRAN is from 10-09-2022, so this change is not included.

[dieghernan]

After digging a bit, I think this is a ggplot2 issue, and IMHO not easy to fix. See:

library(terra)
#> Warning: package 'terra' was built under R version 4.1.3
#> terra 1.6.40
library(tidyterra)
#> Warning: package 'tidyterra' was built under R version 4.1.3
#> 
#> Attaching package: 'tidyterra'
#> The following object is masked from 'package:stats':
#> 
#>     filter

# Mock wkt, see https://github.com/rspatial/terra/blob/11b0b5663a2fe5409b89f3afee1196ffa4b3c765/R/crs.R#L148

local_wkt <- 'LOCAL_CS["Cartesian (Meter)", LOCAL_DATUM["Local Datum",0], UNIT["Meter",1.0], AXIS["X",EAST], AXIS["Y",NORTH]]'

f <- system.file("extdata/volcano2.tif", package = "tidyterra")
r <- rast(f) %>% aggregate(5)


# Add local
crs(r) <- pull_crs(local_wkt)

r
#> class       : SpatRaster 
#> dimensions  : 35, 25, 1  (nrow, ncol, nlyr)
#> resolution  : 25, 25  (x, y)
#> extent      : 1756969, 1757594, 5916998, 5917873  (xmin, xmax, ymin, ymax)
#> coord. ref. : Cartesian (Meter) 
#> source(s)   : memory
#> name        : elevation 
#> min value   :  77.44264 
#> max value   : 192.47299

plot(r)

So far so good. Now I check here how sf behaves with this CRS:

library(sf)
#> Warning: package 'sf' was built under R version 4.1.3
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE
sf_pols <- r %>%
  as.polygons() %>%
  st_as_sf()

sf_pols
#> Simple feature collection with 112 features and 1 field
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: 1756969 ymin: 5917023 xmax: 1757569 ymax: 5917873
#> Projected CRS: Cartesian (Meter)
#> First 10 features:
#>    elevation                       geometry
#> 1         77 POLYGON ((1757344 5917873, ...
#> 2         79 POLYGON ((1757319 5917873, ...
#> 3         80 MULTIPOLYGON (((1757019 591...
#> 4         81 POLYGON ((1757094 5917873, ...
#> 5         82 MULTIPOLYGON (((1757119 591...
#> 6         83 MULTIPOLYGON (((1756994 591...
#> 7         84 MULTIPOLYGON (((1756969 591...
#> 8         85 MULTIPOLYGON (((1757144 591...
#> 9         86 MULTIPOLYGON (((1756969 591...
#> 10        87 MULTIPOLYGON (((1757219 591...



plot(sf_pols["elevation"])

Seems ok too. Problem is when I try to ggplot any of these objects:

# With ggplot

library(ggplot2)

# sf, no tidyterra involved
ggplot(sf_pols) +
  geom_sf()
#> Warning in CPL_transform(x, crs, aoi, pipeline, reverse,
#> desired_accuracy, : GDAL Error 6: Cannot find coordinate
#> operations from `ENGCRS["Cartesian (Meter)",EDATUM["Local
#> Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]'
#> to `ENGCRS["Cartesian (Meter)",EDATUM["Local
#> Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]'
#> Error in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, : OGRCreateCoordinateTransformation(): transformation not available

# Same error than tidyterra
ggplot() +
  geom_spatraster(data = r)
#> Warning in CPL_transform(x, crs, aoi, pipeline, reverse,
#> desired_accuracy, : GDAL Error 6: Cannot find coordinate
#> operations from `ENGCRS["Cartesian (Meter)",EDATUM["Local
#> Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]'
#> to `ENGCRS["Cartesian (Meter)",EDATUM["Local
#> Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]'
#> Error in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, : OGRCreateCoordinateTransformation(): transformation not available

See that the error you mentioned on tidyterra is the same that with the sf object. My guess is that this issue should be located in coord_sf() but don't know how to overcome this neatly.

A potential workaround

This can be implemented in tidyterra. Basically it consists on removing the CRS for plotting and adding a coord_fixed() coordinates. The output is similar to terra::plot() but note that if any other spatial object is added to this map (via geom_sf()/geom_spatraster(), etc.) this layer won't be reprojected to the CRS of these objects, that I think is a cool feature of ggplot2 and may be the root cause of this issue:

# Potential workaround that may be implented?
r_nocrs <- r

crs(r_nocrs) <- NA

ggplot() +
  geom_spatraster(data = r_nocrs) +
  coord_fixed(expand = FALSE) +
  scale_fill_terrain_c()

# Similar to

plot(r)

^{Created on 2022-11-08 with reprex v2.0.2}

[rhijmans]

It seems to me that this should be addressed in ggplot2. That is, there could be a check if a transformation is possible. Such a check exists if the CRS is empty.

library(terra)
library(tidyterra)
r <- rast(system.file("extdata/volcano2.tif", package = "tidyterra")) |> aggregate(5)
crs(r)  <- 'LOCAL_CS["Cartesian (Meter)", LOCAL_DATUM["Local Datum",0], UNIT["Meter",1.0], AXIS["X",EAST], AXIS["Y",NORTH]]'

library(sf)
library(ggplot2)
sf_pols <- st_as_sf(as.polygons(r))

ggplot(sf_pols) + geom_sf()
#Error in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  : 
#  OGRCreateCoordinateTransformation(): transformation not available
#In addition: Warning message:
#In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 6: Cannot find coordinate operations from `ENGCRS["Cartesian (Meter)",EDATUM["Local Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]' to `ENGCRS["Cartesian (Meter)",EDATUM["Local Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]'

This seems to come from sf::st_transform

st_transform(sf_pols, crs(sf_pols))
#Error in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  : 
#  OGRCreateCoordinateTransformation(): transformation not available
#In addition: Warning message:
#In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 6: Cannot find coordinate operations from `ENGCRS["Cartesian (Meter)",EDATUM["Local Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]' to `ENGCRS["Cartesian (Meter)",EDATUM["Local Datum"],CS[Cartesian,2],AXIS["x",east,ORDER[1],LENGTHUNIT["Meter",1]],AXIS["y",north,ORDER[2],LENGTHUNIT["Meter",1]]]'

It also gives an error with an empty CRS, but ggplot2 appears to check for that

st_crs(sf_pols) <- ""
st_transform(sf_pols, "+proj=longlat")
#Error in st_transform.sfc(st_geometry(x), crs, ...) : 
#  cannot transform sfc object with missing crs
ggplot(sf_pols) + geom_sf()

I am guessing that ggplot checks for an empty CRS, but that it should have a more general check whether a transformation is possible (and sf would need to provide that?)

[edzer]

st_crs(NA) (which is what you get when assigning "") is the CRS that is understood by sf as "unknown, Cartesian". Anything else will be assumed as "known, understood, convertible / transformable", and hence leads to the errors above. sf doesn't try to interpret WKT, so has now way to figure out whether a WKT effectively means "unknown, Cartesian". You can see here the WKT that is assigned to objects before writing out as GPKG because GPKG at some stage started to assign 4326 to files with no CRS defined, which is NOT what sf assumes when no CRS is given/known.

geom_sf() calls back to sf for transformations; these are needed to draw the graticules. Setting CRS to NA on the object is the way to go to instruct geom_sf() not to do so.

[edzer]

When reading files with LOCAL_CS["Undefined Cartesian SRS"], sf now assigns CRS NA, with a message.

[rhijmans]

@edzer. It seems odd to me that you would remove perfectly legal WKT. Why not keep it and instead have st_transform check whether the transformation is possible? "terra" uses

poCT = OGRCreateCoordinateTransformation(&source, &target);
if( poCT == NULL ) return FALSE; 

That could also be a separate function such that geom_sf() could call that before deciding whether to call st_transform

[edzer]

There can be many reasons why OGRCreateCoordinateTransformation fails, having a LOCAL_CS["Undefined Cartesian SRS"] is only one of them. You could have for instance ellipsoidal coordinates associated with a different celestial body.

[rhijmans]

And that is why it would be good to do this test?

[edzer]

We discussed in r-spatial/sf#1776 whether to substitute (and keep) LOCAL_CS["Undefined Cartesian SRS"] for NA, as it is more explicit. It just didn't happen; we could still (try to) do it, and have is.na() return TRUE on it, I guess that is what ggplot2 checks.

[edzer]

However,

st_crs('LOCAL_CS["Cartesian (Meter)",LOCAL_DATUM["Local Datum",0],UNIT["Meter",1.0],AXIS["X",EAST],AXIS["Y",NORTH]]') == st_crs('LOCAL_CS["Undefined Cartesian SRS"]')
# FALSE

(this check is done by GDAL SRS->IsSame(), here), so where will this stop?