Issue in spTransform() with rgdal >= 1.5 / PROJ 3

Summary

This draft shows an issue occurring using spTransform() with rgdal >= 1.5, which sometimes causes reprojection failing. In particular, we encountered an issue when one (or both) among input and output CRS are associated to a WKT2 in which the axes order is Y – X instead than X – Y (these two situations were hereinafter called “YX CRS” and “XY CRS”). Issues do not present using the “old” method involving PROJ.4 strings. We debugged the code in order to understand where the error is, finding that it occur in "transform_ng" method called by spTransform(). We were not able to fix it, being in a C method and not in the R code chunks.

Case study

We tested the output reprojections by building an input point in different CRSs, using sf methods (which are not affected by the analysed issues). The analysed CRSs are:

WGS 84 geographic coordinates (EPSG:4326, hereinafter called “LONLAT”);
ETRS89-extended / LAEA Europe (EPSG:3035, hereinafter called “LAEA89”);
WGS 84 / UTM zone 32N (EPSG:32632, hereinafter called “WGS32N”);
WGS 84 / Pseudo-Mercator (EPSG:3857, hereinafter called “PSMERC”).

In this report only these four CRSs were considered, so to analyse 12 cases (\(4^2-4\) cases) of the combinations of the four CRSs. Authors tested a higher number of cases, including two more CRS, in order to check \(6^2-6=30\) cases:

ETRS89 / Austria Lambert (EPSG:3416, hereinafter called “ETRSAT”);
ETRS89 / UTM zone 32N (EPSG:25832, hereinafter called “ETRS32”).

Results of these tests are available here.

# Load libraries
library(rgdal)

## Loading required package: sp

## rgdal: version: 1.5-2, (SVN revision (unknown))
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 3.0.2, released 2019/10/28
##  Path to GDAL shared files: /usr/share/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 6.2.1, November 1st, 2019, [PJ_VERSION: 621]
##  Path to PROJ.4 shared files: /usr/share/proj
##  Linking to sp version: 1.3-3

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1

sf_extSoftVersion()

##           GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
##        "3.8.0"        "3.0.2"        "6.2.1"         "true"         "true"

packageVersion("rgdal")

## [1] '1.5.2'

packageVersion("sf")

## [1] '0.8.1'

(errors analysed in this report appeared using both an ArchLinux OS and on the Docker machine jakubnowosad/geocompr_proj6, with the following package versions installed:

devtools::install_github("rsbivand/sp")
install.packages("rgdal", repos = "https://r-forge.r-project.org")
devtools::install_github("r-spatial/sf", ref = "wkt2")

Points were created:

pt0_lonlat_sf <- st_sfc(st_point(c(10,46)), crs = 4326)
pt0_laea89_sf <- st_transform(pt0_lonlat_sf, 3035)
pt0_wgs32n_sf <- st_transform(pt0_lonlat_sf, 32632)
pt0_psmerc_sf <- st_transform(pt0_lonlat_sf, 3857)

pt0_lonlat <- as(pt0_lonlat_sf, "Spatial")
pt0_laea89 <- as(pt0_laea89_sf, "Spatial")

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown_based_on_GRS80_ellipsoid in CRS definition

pt0_wgs32n <- as(pt0_wgs32n_sf, "Spatial")
pt0_psmerc <- as(pt0_psmerc_sf, "Spatial")

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

pt0_lonlat@coords

##      coords.x1 coords.x2
## [1,]        10        46

pt0_laea89@coords

##      coords.x1 coords.x2
## [1,]   4321000   2543009

pt0_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

pt0_psmerc@coords

##      coords.x1 coords.x2
## [1,]   1113195   5780349

library(mapview)
mapview(pt0_lonlat) + mapview(pt0_laea89) + mapview(pt0_wgs32n) + mapview(pt0_psmerc)

(the disallineation of pt0_laea89 is due to the dischargement of GRS80 ellipsoid)

st_transform(st_as_sfc(pt0_laea89), 4326)[[1]]

## POINT (10 46)

st_transform(st_as_sfc(pt0_wgs32n), 4326)[[1]]

## POINT (10 46)

st_transform(st_as_sfc(pt0_psmerc), 4326)[[1]]

## POINT (10 46)

Then, the points were reprojected in the 12 combinations of CRSs:

pt1_lonlat_laea89 <- spTransform(pt0_lonlat, CRS(SRS_string = "EPSG:3035"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum European_Terrestrial_Reference_System_1989 in CRS definition

pt1_lonlat_wgs32n <- spTransform(pt0_lonlat, CRS(SRS_string = "EPSG:32632"))
pt1_lonlat_psmerc <- spTransform(pt0_lonlat, CRS(SRS_string = "EPSG:3857"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

pt1_laea89_lonlat <- spTransform(pt0_laea89, CRS(SRS_string = "EPSG:4326"))
pt1_laea89_wgs32n <- spTransform(pt0_laea89, CRS(SRS_string = "EPSG:32632"))
pt1_laea89_psmerc <- spTransform(pt0_laea89, CRS(SRS_string = "EPSG:3857"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

pt1_wgs32n_lonlat <- spTransform(pt0_wgs32n, CRS(SRS_string = "EPSG:4326"))
pt1_wgs32n_laea89 <- spTransform(pt0_wgs32n, CRS(SRS_string = "EPSG:3035"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum European_Terrestrial_Reference_System_1989 in CRS definition

pt1_wgs32n_psmerc <- spTransform(pt0_wgs32n, CRS(SRS_string = "EPSG:3857"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

pt1_psmerc_lonlat <- spTransform(pt0_psmerc, CRS(SRS_string = "EPSG:4326"))
pt1_psmerc_laea89 <- spTransform(pt0_psmerc, CRS(SRS_string = "EPSG:3035"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum European_Terrestrial_Reference_System_1989 in CRS definition

pt1_psmerc_wgs32n <- spTransform(pt0_psmerc, CRS(SRS_string = "EPSG:32632"))

The reprojection using spTransform() failed in some cases. Here all the combinations:

Reprojections to LONLAT:

pt0_lonlat@coords # reference

##      coords.x1 coords.x2
## [1,]        10        46

pt1_laea89_lonlat@coords

##      coords.x1 coords.x2
## [1,]   58.3622 -21.52435

pt1_wgs32n_lonlat@coords

##      coords.x1 coords.x2
## [1,]        46        10

pt1_psmerc_lonlat@coords

##      coords.x1 coords.x2
## [1,]        46        10

All reprojections failed.

Reprojections to LAEA:

pt0_laea89@coords # reference

##      coords.x1 coords.x2
## [1,]   4321000   2543009

pt1_lonlat_laea89@coords

##      coords.x1 coords.x2
## [1,] -461198.4   8410923

pt1_wgs32n_laea89@coords

##      coords.x1 coords.x2
## [1,]   2543009   4321000

pt1_psmerc_laea89@coords

##      coords.x1 coords.x2
## [1,]   2543009   4321000

All reprojections failed.

Reprojections to WGS32N:

pt0_wgs32n@coords # reference

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

pt1_lonlat_wgs32n@coords

##      coords.x1 coords.x2
## [1,]   4849477   1378793

pt1_laea89_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  -1244836   6879742

pt1_psmerc_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

Only reprojection PSMERC -> WGS32N was corrected.

Reprojections to PSMERC:

pt0_psmerc@coords # reference

##      coords.x1 coords.x2
## [1,]   1113195   5780349

pt1_lonlat_psmerc@coords

##      coords.x1 coords.x2
## [1,]   5120697   1118890

pt1_laea89_psmerc@coords

##      coords.x1 coords.x2
## [1,]  -2396080   8043792

pt1_wgs32n_psmerc@coords

##      coords.x1 coords.x2
## [1,]   1113195   5780349

Only reprojection WGS32N -> PSMERC was corrected.

In the analysed 12 cases, only in two situations the result was correct: from/to WGS32N to/from PSMERC. These two CRS have are both XY CRSs (see summary for the definitions of XY / YX CRS), while the other two cases are YX CRSs.

cat(showSRID(comment(CRS(SRS_string = "EPSG:4326")), multiline = "YES"))

## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

cat(showSRID(comment(CRS(SRS_string = "EPSG:3035")), multiline = "YES"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum European_Terrestrial_Reference_System_1989 in CRS definition

## PROJCRS["ETRS89-extended / LAEA Europe",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["Europe Equal Area 2001",
##         METHOD["Lambert Azimuthal Equal Area",
##             ID["EPSG",9820]],
##         PARAMETER["Latitude of natural origin",52,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",10,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",4321000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",3210000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (Y)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (X)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Europe - LCC & LAEA"],
##         BBOX[24.6,-35.58,84.17,44.83]],
##     ID["EPSG",3035]]

cat(showSRID(comment(CRS(SRS_string = "EPSG:32632")), multiline = "YES"))

## PROJCRS["WGS 84 / UTM zone 32N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 32N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",9,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World - N hemisphere - 6°E to 12°E - by country"],
##         BBOX[0,6,84,12]],
##     ID["EPSG",32632]]

cat(showSRID(comment(CRS(SRS_string = "EPSG:3857")), multiline = "YES"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## PROJCRS["WGS 84 / Pseudo-Mercator",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["unnamed",
##         METHOD["Popular Visualisation Pseudo Mercator",
##             ID["EPSG",1024]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World - 85°S to 85°N"],
##         BBOX[-85.06,-180,85.06,180]],
##     ID["EPSG",3857]]

Also in the additional 18 analysed cases, the reprojections resulted correct only in four more cases, all involving both input and output XY CRSs (ETRS32 <-> WGS32N and ETRS32 <-> PSMERC, in addition to the exposed PSMERC <-> WGS32N).

Further analysis

Estabilished that the issue appears involving input or output YX CRSs, we checked that the reason was actually the wrong interpretation of these two components. In other words, we checked that the coordinates were inverted before the reprojection in case of input YX CRS and output XY CRS, after the reprojection in case of input XY CRS and output YX CRS, and both before and after the reprojection in case of both input and output YX CRS.

Reprojections from / to XY CRSs:

# Output WGS32N
pt1_psmerc_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

pt0_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

# Output PSMERC
pt1_wgs32n_psmerc@coords

##      coords.x1 coords.x2
## [1,]   1113195   5780349

pt0_psmerc@coords

##      coords.x1 coords.x2
## [1,]   1113195   5780349

As already seen, these reprojections are correct.

Reprojections from XY to YX CRSs:

# Output LONLAT
pt1_wgs32n_lonlat@coords

##      coords.x1 coords.x2
## [1,]        46        10

pt1_psmerc_lonlat@coords

##      coords.x1 coords.x2
## [1,]        46        10

pt0_lonlat@coords

##      coords.x1 coords.x2
## [1,]        10        46

# Output LAEA89
pt1_wgs32n_laea89@coords

##      coords.x1 coords.x2
## [1,]   2543009   4321000

pt1_psmerc_laea89@coords

##      coords.x1 coords.x2
## [1,]   2543009   4321000

pt0_laea89@coords

##      coords.x1 coords.x2
## [1,]   4321000   2543009

In all the cases, the output coordinates corresponds to the correct coordinates with the two X-Y components inverted, confirming that the inversion occurred after the reprojection.

Reprojections from YX to XY CRSs:

# Output WGS32N
pt1_lonlat_wgs32n@coords

##      coords.x1 coords.x2
## [1,]   4849477   1378793

pt1_laea89_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  -1244836   6879742

pt0_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

# Output PSMERC
pt1_lonlat_psmerc@coords

##      coords.x1 coords.x2
## [1,]   5120697   1118890

pt1_laea89_psmerc@coords

##      coords.x1 coords.x2
## [1,]  -2396080   8043792

pt0_psmerc@coords

##      coords.x1 coords.x2
## [1,]   1113195   5780349

The input / output coordinates do not corresponds (the inversion is not done after the reprojection). Instead, inverting the coordinates before the reprojection led to the correct output coordinates:

# Inverted input points
pt1_lonlat_inv_sf <- st_sfc(st_point(pt0_lonlat_sf[[1]][2:1]), crs = st_crs(pt0_lonlat_sf))
pt1_laea89_inv_sf <- st_sfc(st_point(pt0_laea89_sf[[1]][2:1]), crs = st_crs(pt0_laea89_sf))
pt1_lonlat_inv <- as(pt1_lonlat_inv_sf, "Spatial")
pt1_laea89_inv <- as(pt1_laea89_inv_sf, "Spatial")

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown_based_on_GRS80_ellipsoid in CRS definition

# Inverted output points
pt2_lonlat_inv_wgs32n <- spTransform(pt1_lonlat_inv, CRS(SRS_string = "EPSG:32632"))
pt2_laea89_inv_wgs32n <- spTransform(pt1_laea89_inv, CRS(SRS_string = "EPSG:32632"))
pt2_lonlat_inv_psmerc <- spTransform(pt1_lonlat_inv, CRS(SRS_string = "EPSG:3857"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

pt2_laea89_inv_psmerc <- spTransform(pt1_laea89_inv, CRS(SRS_string = "EPSG:3857"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

# Output WGS32N
pt2_lonlat_inv_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

pt2_laea89_inv_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

pt0_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

# Output PSMERC
pt2_lonlat_inv_psmerc@coords

##      coords.x1 coords.x2
## [1,]   1113195   5780349

pt2_laea89_inv_psmerc@coords

##      coords.x1 coords.x2
## [1,]   1113195   5780349

pt0_psmerc@coords

##      coords.x1 coords.x2
## [1,]   1113195   5780349

Reprojections from / to YX CRSs:

# Inverted output points
pt2_laea89_inv_lonlat <- spTransform(pt1_laea89_inv, CRS(SRS_string = "EPSG:4326"))
pt2_lonlat_inv_laea89 <- spTransform(pt1_lonlat_inv, CRS(SRS_string = "EPSG:3035"))

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum European_Terrestrial_Reference_System_1989 in CRS definition

# Output LONLAT
pt2_laea89_inv_lonlat@coords

##      coords.x1 coords.x2
## [1,]        46        10

pt0_lonlat@coords

##      coords.x1 coords.x2
## [1,]        10        46

# Output LAEA89
pt2_lonlat_inv_laea89@coords

##      coords.x1 coords.x2
## [1,]   2543009   4321000

pt0_laea89@coords

##      coords.x1 coords.x2
## [1,]   4321000   2543009

In this cases, inverting the coordinates both before and after the reprojection led to the reference values.

Code debugging

rgdal code was debugged, it resulted that the issue occurred in the transform_ng C method (called in file "R/project.R"):

pt2_lonlat_wgs32n_raw <- .Call(
  "transform_ng", 
  comment(slot(pt0_lonlat, "proj4string")), 
  comment(CRS(SRS_string = "EPSG:32632")),
  NULL, as.integer(1),
  as.double(10), as.double(46), NULL,
  PACKAGE="rgdal"
)
pt2_lonlat_wgs32n <- SpatialPoints(
  matrix(c(pt2_lonlat_wgs32n_raw[[1]], pt2_lonlat_wgs32n_raw[[2]]), ncol=2), 
  proj4string = CRS(SRS_string = "EPSG:32632")
)
pt2_lonlat_wgs32n@coords

##      coords.x1 coords.x2
## [1,]   4849477   1378793

pt0_wgs32n@coords # reference

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

This error occurs also manually specifying the input projection with EPSG code, while it does not occur passing the PROJ.4 string:

pt3_lonlat_wgs32n_raw <- .Call(
  "transform_ng", 
  comment(CRS(SRS_string = "EPSG:4326")),
  comment(CRS(SRS_string = "EPSG:32632")),
  NULL, as.integer(1),
  as.double(10), as.double(46), NULL,
  PACKAGE="rgdal"
)
pt3_lonlat_wgs32n <- SpatialPoints(
  matrix(c(pt3_lonlat_wgs32n_raw[[1]], pt3_lonlat_wgs32n_raw[[2]]), ncol=2), 
  proj4string = CRS(SRS_string = "EPSG:32632")
)
pt3_lonlat_wgs32n@coords

##      coords.x1 coords.x2
## [1,]   4849477   1378793

pt4_lonlat_wgs32n_raw <- .Call(
  "transform_ng", 
  comment(CRS("+init=epsg:4326")),
  comment(CRS(SRS_string = "EPSG:32632")),
  NULL, as.integer(1),
  as.double(10), as.double(46), NULL,
  PACKAGE="rgdal"
)
pt4_lonlat_wgs32n <- SpatialPoints(
  matrix(c(pt4_lonlat_wgs32n_raw[[1]], pt4_lonlat_wgs32n_raw[[2]]), ncol=2), 
  proj4string = CRS(SRS_string = "EPSG:32632")
)
pt4_lonlat_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

pt5_lonlat_wgs32n_raw <- .Call(
  "transform_ng", 
  comment(CRS("+proj=longlat +datum=WGS84 +no_defs")),
  comment(CRS(SRS_string = "EPSG:32632")),
  NULL, as.integer(1),
  as.double(10), as.double(46), NULL,
  PACKAGE="rgdal"
)
pt5_lonlat_wgs32n <- SpatialPoints(
  matrix(c(pt5_lonlat_wgs32n_raw[[1]], pt5_lonlat_wgs32n_raw[[2]]), ncol=2), 
  proj4string = CRS(SRS_string = "EPSG:32632")
)
pt5_lonlat_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

In particular, the first error appeared also using the CRAN version of the package sf, while the previous ones (using spTransform() function and applying "transform_ng" method in the first case) do not (probably, in the CRAN version, spTransform() makes use of the PROJ.4 string).

Manually correcting the axis order, the error disappeared:

wkt2_4326 <- comment(CRS(SRS_string = "EPSG:4326"))
cat(showSRID(wkt2_4326, multiline = "YES"))

## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

wkt2_4326_mod <- paste0(
  "GEOGCRS[\"WGS 84\",",
  "    DATUM[\"World Geodetic System 1984\",",
  "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,",
  "            LENGTHUNIT[\"metre\",1]]],",
  "    PRIMEM[\"Greenwich\",0,",
  "        ANGLEUNIT[\"degree\",0.0174532925199433]],",
  "    CS[ellipsoidal,2],",
  "        AXIS[\"geodetic longitude (Lon)\",east,",
  "            ORDER[1],",
  "            ANGLEUNIT[\"degree\",0.0174532925199433]],",
  "        AXIS[\"geodetic latitude (Lat)\",north,",
  "            ORDER[2],",
  "            ANGLEUNIT[\"degree\",0.0174532925199433]],",
  "    USAGE[",
  "        SCOPE[\"unknown\"],",
  "        AREA[\"World\"],",
  "        BBOX[-90,-180,90,180]],",
  "    ID[\"EPSG\",4326]]"
)
pt6_lonlat_wgs32n_raw <- .Call(
  "transform_ng", 
  comment(CRS(SRS_string = wkt2_4326_mod)),
  comment(CRS(SRS_string = "EPSG:32632")),
  NULL, as.integer(1),
  as.double(10), as.double(46), NULL,
  PACKAGE="rgdal"
)
pt6_lonlat_wgs32n <- SpatialPoints(
  matrix(c(pt6_lonlat_wgs32n_raw[[1]], pt6_lonlat_wgs32n_raw[[2]]), ncol=2), 
  proj4string = CRS(SRS_string = "EPSG:32632")
)
pt6_lonlat_wgs32n@coords

##      coords.x1 coords.x2
## [1,]  577432.2   5094534

The issue does not occur using standalone GDAL:

st_write(pt0_lonlat_sf, pt0_lonlat_sf_path <- tempfile(fileext=".gpkg"))

## Updating layer `file18d3a532800a7' to data source `/tmp/RtmpdDM13v/file18d3a532800a7.gpkg' using driver `GPKG'
## Writing 1 features with 0 fields and geometry type Point.

system(paste0("gdalinfo --version"), intern = TRUE)

## [1] "GDAL 3.0.2, released 2019/10/28"

system(paste0(
  "ogr2ogr -s_srs \"",gsub("\\\"", "\\\\\"", wkt2_4326_mod),"\" ",
  "-t_srs EPSG:32632 ",
  pt7_lonlat_wgs32n_path <- tempfile(fileext=".gpkg")," ",
  pt0_lonlat_sf_path
))
pt7_lonlat_wgs32n <- st_read(pt7_lonlat_wgs32n_path)

## Reading layer `file18d3a532800a7' from data source `/tmp/RtmpdDM13v/file18d3a5504e038.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 0 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 577432.2 ymin: 5094534 xmax: 577432.2 ymax: 5094534
## epsg (SRID):    32632
## proj4string:    +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

pt7_lonlat_wgs32n[[1]]

## Geometry set for 1 feature 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 577432.2 ymin: 5094534 xmax: 577432.2 ymax: 5094534
## epsg (SRID):    32632
## proj4string:    +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

## POINT (577432.2 5094534)

system(paste0(
  "ogr2ogr -s_srs EPSG:4326 ",
  "-t_srs EPSG:32632 ",
  pt8_lonlat_wgs32n_path <- tempfile(fileext=".gpkg")," ",
  pt0_lonlat_sf_path
))
pt8_lonlat_wgs32n <- st_read(pt8_lonlat_wgs32n_path)

## Reading layer `file18d3a532800a7' from data source `/tmp/RtmpdDM13v/file18d3ae885550.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 0 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 577432.2 ymin: 5094534 xmax: 577432.2 ymax: 5094534
## epsg (SRID):    32632
## proj4string:    +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

pt8_lonlat_wgs32n[[1]]

## Geometry set for 1 feature 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 577432.2 ymin: 5094534 xmax: 577432.2 ymax: 5094534
## epsg (SRID):    32632
## proj4string:    +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

## POINT (577432.2 5094534)

system(paste0(
  "ogr2ogr -s_srs \"",gsub("\\\"", "\\\\\"", wkt2_4326),"\" ",
  "-t_srs EPSG:32632 ",
  pt9_lonlat_wgs32n_path <- tempfile(fileext=".gpkg")," ",
  pt0_lonlat_sf_path
))
pt9_lonlat_wgs32n <- st_read(pt9_lonlat_wgs32n_path)

## Reading layer `file18d3a532800a7' from data source `/tmp/RtmpdDM13v/file18d3a2d736b6a.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 0 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 577432.2 ymin: 5094534 xmax: 577432.2 ymax: 5094534
## epsg (SRID):    32632
## proj4string:    +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

pt9_lonlat_wgs32n[[1]]

## Geometry set for 1 feature 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 577432.2 ymin: 5094534 xmax: 577432.2 ymax: 5094534
## epsg (SRID):    32632
## proj4string:    +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

## POINT (577432.2 5094534)

Issue in `spTransform()` with rgdal >= 1.5 / PROJ 3

Luigi Ranghetti, Lorenzo Busetto

2019-12-05

Summary

Case study

Further analysis

Code debugging