2

I'm trying to reproject LiDAR point clouds from California zone 2 (ftUS) using spTransform from sp R package to a metric CRS like EPSG 3857 (initial post here) and got non finite transformation detected: error. Here bellow a reproducible example.

library(sp)
X = c(6275235.88, 6275229.48, 6275163.35)
Y = c(2094771.54, 2094815.15, 2094777.31)
WKT = "COMPD_CS[\"NAD83(2011) / California zone 2 (ftUS) + NAVD88 height - Geoid12B (ftUS)\",PROJCS[\"NAD83(2011) / California zone 2 (ftUS)\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",39.83333333333334],PARAMETER[\"standard_parallel_2\",38.33333333333334],PARAMETER[\"latitude_of_origin\",37.66666666666666],PARAMETER[\"central_meridian\",-122],PARAMETER[\"false_easting\",6561666.667],PARAMETER[\"false_northing\",1640416.667],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"6418\"]],VERT_CS[\"NAVD88 height - Geoid12B (ftUS)\",VERT_DATUM[\"North American Vertical Datum 1988\",2005,AUTHORITY[\"EPSG\",\"5103\"]],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AUTHORITY[\"EPSG\",\"6360\"]]]"
XY = cbind(X,Y)
XY = sp::SpatialPoints(XY, sp::CRS(SRS_string = WKT))
sp::spTransform(XY, sp::CRS("+init=epsg:3857"))
#> Warning in sp::spTransform(XY, sp::CRS("+init=epsg:3857")): 3 projected point(s)
#> not finite
#> non finite transformation detected:
#>            X       Y        
#> [1,] 6275236 2094772 Inf Inf
#> [2,] 6275229 2094815 Inf Inf
#> [3,] 6275163 2094777 Inf Inf
#> Error in sp::spTransform(XY, sp::CRS("+init=epsg:3857")): failure in points 1:2:3

Basically, in this project I'm dealing with various CRS across the world and I don't need absolute Z accuracy but I really want to get rid of US feet survey CRS !! The EPSG 3857 is convenient because it is a worldwide metric datum but it is not the cause of the error. I received the same error using local datum like EPSG 6417 and 6418. As suggested by @mkennedy in this post, the Z component of the CRS transformation could be the cause of the error and I am wandering how to use spTransform just for the horizontal CRS.

In other words, is it possible to ask to spTransform to focus on the horizontal component of the targeted CRS? I know this is related to the proj library but I'm looking for a solution within R environment. I'm open to solutions based on other R package like sf

Edit: the problem is reproducible with a Z coordinate or with sf::st_transform(). But sf doesn't fail. It returns NAs instead

X = c(6275235.88, 6275229.48, 6275163.35)
Y = c(2094771.54, 2094815.15, 2094777.31)
Z = c(10,15,20)
WKT = "COMPD_CS[\"NAD83(2011) / California zone 2 (ftUS) + NAVD88 height - Geoid12B (ftUS)\",PROJCS[\"NAD83(2011) / California zone 2 (ftUS)\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",39.83333333333334],PARAMETER[\"standard_parallel_2\",38.33333333333334],PARAMETER[\"latitude_of_origin\",37.66666666666666],PARAMETER[\"central_meridian\",-122],PARAMETER[\"false_easting\",6561666.667],PARAMETER[\"false_northing\",1640416.667],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"6418\"]],VERT_CS[\"NAVD88 height - Geoid12B (ftUS)\",VERT_DATUM[\"North American Vertical Datum 1988\",2005,AUTHORITY[\"EPSG\",\"5103\"]],UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],AUTHORITY[\"EPSG\",\"6360\"]]]"
XYZ = cbind(X,Y,Z)
XYZ = sp::SpatialPoints(XYZ, sp::CRS(SRS_string = WKT))
sp::spTransform(XYZ, sp::CRS("+init=epsg:3857"))

XYZ = sf::st_as_sf(XYZ)
sf::st_transform(XYZ, sf::st_crs("EPSG:3857"))
3
  • 1
    I can't get that WKT string to seem valid to anything in sp or sf. ANd my sp::CRS doesn't have an SRS_string parameter. What version of everything (packages, gdal, proj, as seen by library(sf) output) are you on? Edit these into your question?
    – Spacedman
    Oct 2, 2020 at 17:07
  • 1
    The infinities here might be because you haven't specified a Z (height) for your points in XY. I think if you coordinate system is something like EPSG:8715 then you need coordinates and height values: epsg.io/8715
    – Spacedman
    Oct 2, 2020 at 17:14
  • 1
    @Spacedman. SRS_string is new with recent version of sp and needs proj6, gdal3 and is part of the migration to WTK representation of CRS in R (see also many doc written by Roger Bivand). It should be the preferred way to specify CRS. I tried with Z coordinates as well but it fails too.
    – JRR
    Oct 2, 2020 at 20:45

1 Answer 1

4

You can turn this into a 2D transform by removing the VERT_CS component of the WKT.

WKT2 = "COMPD_CS[\"NAD83(2011) / California zone 2 (ftUS) + NAVD88 height - Geoid12B (ftUS)\",
PROJCS[\"NAD83(2011) / California zone 2 (ftUS)\",
GEOGCS[\"NAD83(2011)\",
DATUM[\"NAD83_National_Spatial_Reference_System_2011\",
SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],
AUTHORITY[\"EPSG\",\"1116\"]],
PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],
UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],
AUTHORITY[\"EPSG\",\"6318\"]],
PROJECTION[\"Lambert_Conformal_Conic_2SP\"],
PARAMETER[\"standard_parallel_1\",39.83333333333334],
PARAMETER[\"standard_parallel_2\",38.33333333333334],
PARAMETER[\"latitude_of_origin\",37.66666666666666],
PARAMETER[\"central_meridian\",-122],
PARAMETER[\"false_easting\",6561666.667],
PARAMETER[\"false_northing\",1640416.667],
UNIT[\"US survey foot\",0.3048006096012192,AUTHORITY[\"EPSG\",\"9003\"]],
AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"6418\"]]]"

and then I only get the warning messages and what looks like a sensible transformation of the X and Y coordinates:

> XYZ = sp::SpatialPoints(XYZ, sp::CRS(SRS_string = WKT2))
Warning message:
In showSRID(SRS_string, format = "PROJ", multiline = "NO") :
  Discarded datum NAD83_National_Spatial_Reference_System_2011 in CRS definition
> sp::spTransform(XYZ, sp::CRS("+init=epsg:3857"))
SpatialPoints:
             X       Y       Z
[1,] -13693038 4708772 2.0e+07
[2,] -13693041 4708789 1.5e+01
[3,] -13693066 4708774 2.0e+01
Coordinate Reference System (CRS) arguments: +proj=merc +a=6378137
+b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m
+nadgrids=@null +wktext +no_defs 
Warning messages:
1: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
2: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum WGS_1984 in CRS definition

the heights being unchanged (so multiply by meter-per-USft to fix).

By sensible I mean if I transform them to lat-long I get something that seems to be around California somewhere...

> sp::spTransform(XYZ, sp::CRS("+init=epsg:4326"))
SpatialPoints:
             X        Y       Z
[1,] -123.0067 38.90989 2.0e+07
[2,] -123.0067 38.91001 1.5e+01
[3,] -123.0069 38.90990 2.0e+01

I don't know why a VERT_CS messes this up, possibly worth a bug report, but I'm not sure who to...

Versions:

> library(sf)
Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

Also I don't know if the vertical coordinate and horizontal components are independent. In other words if transforming (X,Y,Z) to (X', Y', Z') can be written X'=f(X,Y), Y'=g(X,Y), Z'=h(Z). If not then dropping the vertical coordinate system might result in errors because of a difference in the definition of "vertical" in different systems.

To investigate if this is a bug, I used R to convert your 3d CRS to a Proj string, and feed into the proj command line tools to see what it did:

# cs2cs +proj=lcc +lat_0=37.6666666666667 +lon_0=-122 +lat_1=39.8333333333333 +lat_2=38.3333333333333 +x_0=2000000.0001016 +y_0=500000.0001016 +ellps=GRS80 +units=us-ft +vunits=us-ft +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs
6275235.88 2094771.54 100
-13693037.82    4708772.00 30.48

And this looks correct, even down to the change of height from 100 us-ft to 30.48 metres. Feeding this proj string into R still returns empty points on transform:

> WKTc = st_crs(WKT)
> projstring = as(WKTc,"CRS")
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
> projstring
CRS arguments:
 +proj=lcc +lat_0=37.6666666666667 +lon_0=-122 +lat_1=39.8333333333333
+lat_2=38.3333333333333 +x_0=2000000.0001016 +y_0=500000.0001016
+ellps=GRS80 +units=us-ft +vunits=us-ft +no_defs 
> pts = st_as_sf(data.frame(x=6275235.88,y= 2094771.54,z=100), coords=1:3, crs=projstring)
> st_transform(pts, "+init=epsg:3857")
Simple feature collection with 1 feature and 0 fields (with 1 geometry empty)
geometry type:  POINT
dimension:      XY
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
projected CRS:  WGS 84 / Pseudo-Mercator
     geometry
1 POINT EMPTY

which seems to be because R has converted the proj string back to a WKT. Given that (I think) proj internally only deals with proj strings, it might be a bug in the way R is passing the coordinates to Proj, because it seems capable of generating the correct proj string (or my cs2cs command line wouldn't work).

2
  • Thanks !!! To whom should this potential bug be relayed? I think my knowledge is too limited to properly document this. Maybe @Spacedman could post a new post on gis.stackexchange based on last reply?
    – Adrien
    Oct 4, 2020 at 18:34
  • For those who want to remove VERT_CS programmatically from WKT, check this nice explanation: stackoverflow.com/questions/64195441/…
    – Adrien
    Oct 4, 2020 at 18:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.