sf
sf copied to clipboard
r-spatial
`st_transform(, ...pipeline=)` issue?
I'm trying to check that rgdal, terra @rhijmans and sf not only find pipelines properly, but that they apply them properly:
f <- system.file("ex/lux.shp", package="terra")
v <- terra::vect(f)
crs <- "EPSG:2169"
p <- terra::project(v, crs)
terra::ext(p)
s1 <- rgdal::readOGR(f)
sp::bbox(sp::spTransform(s1, sp::CRS("EPSG:2169")))
rgdal::get_last_coordOp()
o1 <- rgdal::list_coordOps("EPSG:4326", "EPSG:2169")
o1$definition[1]
sp::bbox(sp::spTransform(s1, sp::CRS("EPSG:2169"), pipeline=o1$definition[1]))
s2 <- sf::st_read(f)
sf::st_bbox(sf::st_transform(s2, "EPSG:2169"))
o2 <- sf::sf_proj_pipelines("EPSG:4326", "EPSG:2169")
o2$definition[1]
sf::st_bbox(sf::st_transform(s2, "EPSG:2169", pipeline=o2$definition[1])) # bbox and geoms wrong
So sf::st_transform() (using GDAL not PROJ) finds the best available transformation (as do terra::project() and sp::spTransform()), but passing the pipeline defining this transformation does something wierd (sf::st_transform() C++ function at https://github.com/r-spatial/sf/blob/898b51e3f7f8a56150468d6a208b290859a313e0/src/gdal.cpp#L532-L594).
This matters as I'd like to answer https://github.com/r-spatial/sf/issues/1933 properly by showing that grids are accessed from the CDN in sf::st_transform(), terra::project() and sp::spTransform() in the same way cross-platform.
I got stuck on this simple non-grid case. Since neither sf nor terra report the pipeline anyway, I have to compare actions across packages. It may be worth finding out how to report applied pipelines from GDAL; for example, can an OGRCoordinateTransformation object be reported as a string PROJ pipeline like the definition component @rouault?
As things stand, sp::spTransform() only uses PROJ for transformation, it reports the last pipeline used, and when the pipeline is given as as argument, makes the same transformation if the given pipeline is the best available. sf::st_transform() and terra::project() use GDAL, find an unreported pipeline, and here make the same transformation. sf::st_transform() can take a pipeline argument (from for example sf::sf_proj_pipelines()), but passing the best available pipeline gives, here, a weird transformation.
It is caused by GDAL requiring "The pipeline must take into account the axis order of the source and target SRS." in SetCoordinateOperation():
od1 <- o2$definition[1]
aswp <- "+step +proj=axisswap +order=2,1"
od1a <- paste0(substring(od1, 1, 15), aswp, substring(od1, 15), " ", aswp)
sf::st_bbox(sf::st_transform(s2, "EPSG:2169", pipeline=od1)) # wrong
sf::st_bbox(sf::st_transform(s2, "EPSG:2169", pipeline=od1a)) # correct
but with:
Warning message:
In st_transform.sfc(st_geometry(x), crs, ...) :
pipeline not found in PROJ-suggested candidate transformations
because the PROJ-generated pipeline in sf seems to apply traditional GIS order (rgdal applies it symmetrically), but:
> system("projinfo -s EPSG:4326 -t EPSG:2169")
Candidate operations found: 2
-------------------------------------
Operation No. 1:
unknown id, Inverse of LUREF to WGS 84 (4) + Luxembourg TM, 1 m, Luxembourg.
PROJ string:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +proj=push +v_3
+step +proj=cart +ellps=WGS84
+step +inv +proj=molobadekas +x=-265.8867 +y=76.9851 +z=20.2667 +rx=0.33746
+ry=3.09264 +rz=-2.53861 +s=0.4598 +px=4103620.3943 +py=440486.4235
+pz=4846923.4558 +convention=coordinate_frame
+step +inv +proj=cart +ellps=intl
+step +proj=pop +v_3
+step +proj=tmerc +lat_0=49.8333333333333 +lon_0=6.16666666666667 +k=1
+x_0=80000 +y_0=100000 +ellps=intl
+step +proj=axisswap +order=2,1
...
sf has sf::st_axis_order() default FALSE. @edzer : curiously terra and sp agree on the x/y names of the bbox, but sf reverses them?
can an
OGRCoordinateTransformationobject be reported as a string PROJ pipeline like the definition component @rouault?
the issue is that there isn't necessarily a single pipeline for a given OGRCoordinateTransformation. GDAL will defer to PROJ to select the "best" pipeline for each coordinate tuple to transform. At least when the user doesn't set an explicit pipeline on OGRCoordinateTransformation but just the source and target CRS
Can one retrieve the pipeline actually applied (the best available used)?
not at the GDAL level.
at the PROJ level, using the following functions could give you what you want (well, to be exact, very close to what proj_trans() will do. There are exceptional situations where the transformation used by proj_trans() might not match the one returned by proj_get_suggested_operation(). Like the one mentioned at https://github.com/OSGeo/PROJ/blob/058756c697f5a53b958455916018043cecbcade5/src/4D_api.cpp#L286)
PJ_OBJ_LIST PROJ_DLL *proj_create_operations(
PJ_CONTEXT *ctx,
const PJ *source_crs,
const PJ *target_crs,
const PJ_OPERATION_FACTORY_CONTEXT *operationContext);
int PROJ_DLL proj_get_suggested_operation(PJ_CONTEXT *ctx,
PJ_OBJ_LIST *operations,
PJ_DIRECTION direction,
PJ_COORD coord);
PJ PROJ_DLL *proj_list_get(PJ_CONTEXT *ctx,
const PJ_OBJ_LIST *result,
int index);
void PROJ_DLL proj_list_destroy(PJ_OBJ_LIST *result);
In rgdal I use proj_pj_info(pj_transform).definition) as a string from either pj_transform = proj_create(PJ_DEFAULT_CTX, CHAR(STRING_ELT(coordOp, 0))) if coordOp is given, or pj_transform = proj_create_crs_to_crs_from_pj(PJ_DEFAULT_CTX, source_crs, target_crs, area_of_interest, NULL) if only source, target and opionally AOI. However, in rgdal by default pj_transform is passed through proj_normalize_for_visualization(). If I set the non-default value:
> o1 <- rgdal::list_coordOps("EPSG:4326", "EPSG:2169", visualization_order=FALSE)
> o1$definition[1]
[1] "+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=molobadekas +x=-265.8867 +y=76.9851 +z=20.2667 +rx=0.33746 +ry=3.09264 +rz=-2.53861 +s=0.4598 +px=4103620.3943 +py=440486.4235 +pz=4846923.4558 +convention=coordinate_frame +step +inv +proj=cart +ellps=intl +step +proj=pop +v_3 +step +proj=tmerc +lat_0=49.8333333333333 +lon_0=6.16666666666667 +k=1 +x_0=80000 +y_0=100000 +ellps=intl +step +proj=axisswap +order=2,1"
I get the same output as projinfo. The odd behaviour of sf::st_transform() is an interaction with axis policy that we need to think through. rgdal gives the insight needed (hopefully):
> sp::bbox(sp::spTransform(s1, sp::CRS("EPSG:2169"), enforce_xy=TRUE))
min max
x 49540.31 105922.0
y 57009.53 138631.1
> rgdal::get_last_coordOp()
[1] "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=molobadekas +x=-265.8867 +y=76.9851 +z=20.2667 +rx=0.33746 +ry=3.09264 +rz=-2.53861 +s=0.4598 +px=4103620.3943 +py=440486.4235 +pz=4846923.4558 +convention=coordinate_frame +step +inv +proj=cart +ellps=intl +step +proj=pop +v_3 +step +proj=tmerc +lat_0=49.8333333333333 +lon_0=6.16666666666667 +k=1 +x_0=80000 +y_0=100000 +ellps=intl"
> sp::bbox(sp::spTransform(s1, sp::CRS("EPSG:2169"), enforce_xy=FALSE))
min max
x 5391667 5505686
y -4544033 -4425804
> rgdal::get_last_coordOp()
[1] "+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=molobadekas +x=-265.8867 +y=76.9851 +z=20.2667 +rx=0.33746 +ry=3.09264 +rz=-2.53861 +s=0.4598 +px=4103620.3943 +py=440486.4235 +pz=4846923.4558 +convention=coordinate_frame +step +inv +proj=cart +ellps=intl +step +proj=pop +v_3 +step +proj=tmerc +lat_0=49.8333333333333 +lon_0=6.16666666666667 +k=1 +x_0=80000 +y_0=100000 +ellps=intl"
with the odd output being caused by unmatched axisswaps.
So, the problem was the use of "EPSG:4326" when data was eastings, northings:
f <- system.file("ex/lux.shp", package="terra")
v <- terra::vect(f)
cat(terra::crs(v), "\n")
terra::crs(v) <- "OGC:CRS84"
cat(terra::crs(v), "\n")
p <- terra::project(v, crs)
terra::ext(p)
s1 <- rgdal::readOGR(f)
slot(s1, "proj4string")
slot(s1, "proj4string") <- sp::CRS("OGC:CRS84")
slot(s1, "proj4string")
sp::bbox(sp::spTransform(s1, sp::CRS("EPSG:2169"), enforce_xy=TRUE))
rgdal::get_last_coordOp()
sp::bbox(sp::spTransform(s1, sp::CRS("EPSG:2169"), enforce_xy=FALSE))
rgdal::get_last_coordOp()
s2 <- sf::st_read(f)
sf::st_crs(s2)
sf::st_crs(s2) <- sf::st_crs("OGC:CRS84")
sf::st_crs(s2)
o3 <- sf::sf_proj_pipelines("OGC:CRS84", "EPSG:2169")
o3$definition[1]
sf::st_bbox(sf::st_transform(s2, "EPSG:2169", pipeline=o3$definition[1]))
system("projinfo -s OGC:CRS84 -t EPSG:2169")
The northing, easting axis order of EPSG:4326 is certainly confusing, and further, lux.shp prevents the storage of the correct axis order, because WKT1 doesn't write it. Use of GPKG preserves axis order. @rhijmans do you want a PR for terra?
curiously terra and sp agree on the x/y names of the bbox, but sf reverses them?
f <- system.file("ex/lux.shp", package="terra")
v <- terra::vect(f)
crs <- "EPSG:2169"
p <- terra::project(v, crs)
terra::ext(p)
# SpatExtent : 49540.3062979379, 105922.009918318, 57009.5323630146, 138631.128143484 (xmin, xmax, ymin, ymax)
s1 <- rgdal::readOGR(f)
# OGR data source with driver: ESRI Shapefile
# Source: "/home/edzer/R/x86_64-pc-linux-gnu-library/4.0/terra/ex/lux.shp", layer: "lux"
# with 12 features
# It has 6 fields
# Integer64 fields read as strings: POP
bb = sp::bbox(sp::spTransform(s1, sp::CRS("EPSG:2169")))
# Warning message:
# In showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
# Discarded datum Luxembourg 1930 in Proj4 definition
as.vector(bb)
# [1] 49540.31 57009.53 105922.01 138631.13
s2 <- sf::read_sf(f)
sf::st_crs(s2) = "OGC:CRS84"
sf::st_bbox(sf::st_transform(s2, "EPSG:2169"))
# xmin ymin xmax ymax
# 49540.31 57009.53 105922.01 138631.13
sp bounding boxes are matrices with unnamed elements; their order matches that of sf bbox vector.
I guess the st_transform needs to have a bit of text in the Details section added on the axis order mess when pipeline is specified (as of: do it yourself).
Add https://github.com/r-spatial/sf/commit/2f08ef1d9d66a8bdd4894aac6d61366bc8702a5d manually. For reasons unknown, github insists on tracking r-spatial-sf, not my fork, so PRs are not possible, because they include all the mess I make trying to keep up with your repo. Combining roxygen and github is too much for me. I'll pen a blog in my time on this, it is serious and as we saw with the Toronto tower will cause plenty of user errors.
How do you generate the Rd files from the doxygen markup? Can one generate incrementally (only from a modified file)?
I tried that, it just dumped all the *.o etc in source but didn't update the man file (I did try to limit to only Rd files with their roclet, but maybe that was the problem). Without setting the roclet, it polluted the src/ directory, but updated the only modified target Rd output. Can it be told not to install in-tree?