-
Notifications
You must be signed in to change notification settings - Fork 8
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add support to crs "local" (recently defined in terra) #64
Comments
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 infosessioninfo::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)
#>
#>
#> ------------------------------------------------------------------------------ |
Thank you for including it.
According to rspatial/terra#797 (comment) |
The latest version on CRAN is from 10-09-2022, so this change is not included. |
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 A potential workaroundThis can be implemented in tidyterra. Basically it consists on removing the CRS for plotting and adding a # 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 |
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.
This seems to come from
It also gives an error with an empty CRS, but ggplot2 appears to check for that
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?) |
|
When reading files with |
@edzer. It seems odd to me that you would remove perfectly legal WKT. Why not keep it and instead have
That could also be a separate function such that |
There can be many reasons why |
And that is why it would be good to do this test? |
We discussed in r-spatial/sf#1776 whether to substitute (and keep) |
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? |
The test doesn't tell you whether the coordinates are to be treated as unknown, Cartesian. |
The test would tell you that you cannot transform the data. That is equivalent to having a crs that is NA in the context of making a map. Otherwise, While returning NA could solve the problem for some cases; you still would not be able to use ggplot for mapping an area on Mars? (Unless you set the crs to NA)? |
That is not how things work currently work with
Of course, it is just that we have treated them identical over the last 2 decades.
Once PROJ gets equipped with datums for other bodies, I see no problems there. Maybe it already has, or someone has added them locally to proj.db ... |
of course -
|
My assumption is that ggplot calls |
See r-spatial/sf#2049:
as a missing planar representation. Vector GPKG mis-assigns a GEOGCRS with unknown datum even if the coordinates are far outside the valid ranges for geographical coordinates. |
# version 1.0-13 * `gdal_utils()` adds `"ogrinfo"` utility (requires GDAL >= 3.7.0); #2160 * `st_as_sf()` catches errors when setting invalid crs values, raised by Jon Skøien * add `rename_with.sf()` method; #1472 * use GEOS' overlayNG routines for (GEOS) Intersection, Difference, Union and SymDifference; #2143 * added `duplicated.sf()`; #2138, #2140, thanks to @bart1 * `select.sf()` allows selecting the same column twice under different names; #1886 * `st_as_sf.ppplist()` is deprecated; #1926 * `st_cast()` handles empty geometries; #1961 * don't repeat longlat messages in `summarise.sf()`; #1519 * fix random sampling on the sphere; #2133 # version 1.0-12 * update NAMESPACE to `useDynLib(sf, .registration=TRUE)`; #2127 thanks to @eddelbuettel * fix call in `gdal_addo()`; #2124 * fix issues that came up with older GDAL version, use `GDAL_VERSION_NUM` consistently; #2123 #2121 #2119 # version 1.0-11 * `st_make_grid()` also accepts area units for `cellsize`, for square and hexagonal grids; #1505 * add `st_concave_hull()`, for concave hulls, if GEOS >= 3.11.0; #1964 * add `st_triangulate_constrained()`, for constrained Delaunay triangulation, if GEOS >= 3.10.0; #1964 * clean up the retrieval of length or angle units from WKT representations; https://lists.osgeo.org/pipermail/gdal-dev/2023-March/056994.html * conversion to GEOS uses the `GEOS_PREC_VALID_OUTPUT` flag, which makes sure that the "[o]utput is always valid. Collapsed geometry elements (including both polygons and lines) are removed." # version 1.0-10 * `gdal_utils()` has a `config_options` argument to set further GDAL options, just like `st_write()`; #2003 * fix slow writing of logical vectors in `st_write()`; #1409; #1689 * `st_drivers()` has an argument `regex` to filter on driver (long) name; #2090 * drop C++11 as a system requirement * `c.sfc()` (and, consequently, `dplyr::bind_rows()`) gives an error if components have different CRS; #1884 * data imported from `maps` are associated with the Clark 1866 ellipsoid; #2080 * fix importing legacy `SpatialPolygon` objects without comments; #2063, #2069, rstudio/leaflet#833 * `st_read()` no longer errors on mixes of `XY` and `XYZ` geometries; #2046 #1592 * in `plot.sf()`, when numeric `breaks` are given a legend key is always plotted; #2065 * `st_crs()$axes` returns a `data.frame` with axes properties (name, orientation, conversion factor) when GDAL >= 3.0.0 * clean up unit handling for geometry measures (length, area, distance) and crs; * `st_crs(x)$ud_unit` returns `NULL` if units are unknown; #2049 * `st_write()` substitutes an `NA` crs with `ENGCRS["Undefined Cartesian SRS with unknown unit"]`; #2049, #2054 * `st_can_transform()` checks whether a transformation between two crs exists; see dieghernan/tidyterra#64; #2049 * the matrix returned by `st_coordinates()` has no row names, to reduce output size # version 1.0-9 * adjust for changes how R-devel handles `POSIXlt`; #2028 * add `st_break_antimeridian()`; #1983, #1991 by Roger Bivand * add `Fibonacci` as a spatial sampling type in `st_sample()` * use the global `options("sf_use_s2")` to determine whether to use s2, rather than a value in a local environment; #1977 * fix utils `mdiminfo` and `mdimtranslate` in `gdal_utils()` * extend arguments of `gdal_read_mdim()` needed by `stars::read_mdim()` if `stars` >= 0.5-7; add `gdal_write_mdim()` * add `drop_na()` method for `sf` objects; #1975 # version 1.0-8 * `st_drop_geometry.default()` returns `x` unmodified; * `sf_project()` accepts 3- or 4-column matrices, containing z and t values; * optimizations for `st_sfc()` by @paleolimbot; #1938, #1925 * `[<-.sfc()` recomputes the bounding box; `st_sfc()` gets parameter `compute_bbox`; #1965 * add new algorithm and drop option to `st_make_valid()` when using GEOS and GEOS >= 3.10.1; #1655 * add `st_minimum_rotated_rectangle()`, available when GEOS >= 3.9.0; #1953 * fix `st_sample()` with `type = "hexagonal"` for corner case (n=1), add a `progress` argument for a progress bar; #1945 * add package `pbapply` to Suggests; #1945 * add pdf driver to windows build; #1942 * clarify `pipeline` argument in `st_transform()` when axis order is ambiguous; #1934 * handle argument `xpd` in calls to `plot.sfc_POLYGON()` and `plot.sfc_MULTIPOLYGON()` * add `pivot_wider()` method, by Henning Teickner; #1915 * add `gdal_addo()` to add or remove overviews from raster images; #1921 * `st_layers()` returns `crs` of each layer in a `crs` list of `crs` objects * restore `st_graticule()` behaviour to pre-sf 1.0-0; tidyverse/ggplot2#4571 * `gdal_metadata()` sets metadata item names properly * `st_read()` gains an argument `optional` passed on to `as.data.frame` to avoid changing column names; #1916 * GPX files are autodetected by `st_read()`; #1917 * unnecessary coordinate names are not returned in `st_sample()`, making the output size smaller; #1879 # version 1.0-7 * `st_drop_geometry()` is a generic; #1914 * `st_crs(x)$ud_unit` returns the unit of the coordinate reference system of `x` * geometric predicates return `sgbp` objects omitting self-intersections etc. by passing `remove_self = TRUE` and unique symmetric relationship by passing `retain_unique = TRUE` (to `...` if needed); this simplifies identifying (and removing) duplicated geometries; duplicates are identified by e.g. by `st_equals(x, retain_unique = TRUE) |> unlist() |> unique()`; #1893 * fix compile issue against GDAL < 2.5.0 introduced in 1.0-6; #1899 # version 1.0-6 * adapt to new `spatstat.random` package; #1892 * `st_geometry<-()` also allows to rename a geometry column in an `sf` object; #1890 * for `sf` objects, the `st_as_sfc()` method is an alias for `st_geometry()`; #1882 * `st_make_grid()` speeded up; #1579 thanks to Krzysztof Dyba * remove direct and indirect dependencies on `rgeos` and `rgdal`; #1869 * use `stats::dist` rather than GEOS for symmetric point-point Euclidian distance computation; #1874 # version 1.0-5 * package startup message reports status of `sf_use_s2()`; #1782 * `sf_use_s2()` uses `message()` to report a change; #1782 * `st_sample()` using regular sampling for ellipsoidal coordinates "works" as if coordinates were Cartesian; #1837 # version 1.0-4 * new function `st_delete()` deletes a data source, or layer(s) within a data source; #1828 * fix memory leak in `WKT1_ESRI` retrieval; #1690 # version 1.0-3 * cope with how GEOS >= 3.10.0 handles illegal geometries (e.g., non-closed rings); #1807 * `crs` objects have a `$srid` method to extract the SRID (as authority "name:code"); #1804 * `st_as_grob()` methods for `sfc_*` objects correctly handle empty geometries; #1789 with help from Hiroaki Yutani * when writing objects with `NA` as CRS to GeoPackage, assign "Unknown Cartesian CRS" first - this is in line with using Cartesian geometry operations for objects with `NA` as CRS; #1776 * add coerce method from `sgbp` to `sparseMatrix`; #1750 * fix `st_cast()` for `GEOMETRYCOLLECTIONS` containing empty geometries; #1767 * fix `st_is_valid()` for bogus polygons and projected coordinates; #1666, #1760; #1761
Please add support to crs "local"
rspatial/terra#797
Given:
The following works as expected:
But:
The text was updated successfully, but these errors were encountered: