point2las: propagate SRS?

[jlaura]

Describe the bug point2las appears to be stripping rich projection information and replacing it with a truncated proj4 string. Is this intended?

To Reproduce Steps to reproduce the behavior. Take a lunar PC that one has lying around (I used Kaguya data) and run point2las with it a la:

point2las PC.tif"   \
--compressed   \
-o out  \
--t_srs "GEOGCRS[\"Moon(2015)-Sphere/Ocentric\",DATUM[\"Moon(2015)-Sphere\",ELLIPSOID[\"Moon(2015)-Sphere\",1737400,0,LENGTHUNIT[\"metre\",1]],ID[\"IAU\",30100,2015]],PRIMEM[\"ReferenceMeridian\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"IAU\",30100,2015]],CS[ellipsoidal,3],AXIS[\"geodeticlatitude(Lat)\",north,ORDER[1],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"geodeticlongitude(Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]],AXIS[\"ellipsoidalheight(h)\",up,ORDER[3],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],REMARK[\"Promotedto3DfromIAU_2015:30100.SourceofIAUCoordinatesystems:https://doi.org/10.1007/s10569-017-9805-5\"]]"

Take the resulting LAZ file and run PROJ_IGNORE_CELESTIAL_BODY=YES pdal info <the laz file>.

The resulting SRS is reported as:

"GEOGCS[\"unknown\",DATUM[\"unknown\",SPHEROID[\"unknown\",1737400,0]],PRIMEM[\"Reference meridian\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST]]",

and

    "spatialreference": "GEOGCS[\"unknown\",DATUM[\"unknown\",SPHEROID[\"unknown\",1737400,0]],PRIMEM[\"Reference meridian\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST]]".

Note that I tested passing the projection fully escaped (as in the command above) and with all \" escape/quotes removed. The resulting projection reported from proj info was the same.

One thing that I am concerned about opening this bug report is that the issue could be with pdal either in the write or in the info call. I figured I would start here, to confirm that the rich SRS is being written and then, as needed link this issue to a PDAL issue.

Expected behavior I would expect the passed SRS to be present / reported in the resulting LAZ file.

Your Environment (please complete the following information):

• OS: Linux
• ASP Version: 1-23-24 nightly
  • Any other environment information that might be helpful?

Additional context I can upload (and email a download link) as needed. The DTMs that I have are ~55MB and git caps uploads on issues at 25MB.

[dshean]

Oof. PDAL should use latest PROJ9 and respect CRS definitions, but I have not tested for planetary bodies. My guess is the ASP PDAL call is responsible, maybe here: https://github.com/NeoGeographyToolkit/StereoPipeline/blob/f2bcbd219d6787961f189e9a16c5f5620c353b0d/src/asp/Tools/point2las.cc#L280

This brings up the larger issue to update ASP CRS support to PROJ9.

Just the other day, I noticed that pc_align is still relying on proj4 strings to define CRS of points in csv file.
I did a quick search and many ASP utilities currently rely on this: https://stereopipeline.readthedocs.io/en/latest/search.html?q=proj4

I think anywhere ASP currently expects proj4, we should accept any supported CRS definition (EPSG codes, WKT2, proj4, whatever...see list of options here https://proj.org/en/9.2/apps/projinfo.html) and leverage the latest PROJ9 parser.

[jlaura]

I think anywhere ASP currently expects proj4, we should accept any supported CRS definition (EPSG codes, WKT2, proj4, whatever...see list of options here https://proj.org/en/9.2/apps/projinfo.html) and leverage the latest PROJ9 parser.

Agreed. I had (erroneously) assumed that that was the case. Basically, pass anything PROJ will ingest and then forget about it.

Ultimately, I am trying to create a huge merged point cloud using entwine and I suspect that PROJ issues exist through the entire stack. The issues that may be in ASP are just the first step of getting a large, merged point cloud in a cloud native form.

[oleg-alexandrov]

This will require some work to sort out. Is there no suitable PROJ.4 string that may be good enough?

[thareUSGS]

Here is what I posted the other day to a OpenPlanetary channel. But the hope is that you don't need to deal with it. Let PROJ deal with it (which means it will require a newer version -- if that is the issue).

WKT2 and/or PROJJSON are likely the methods to use going forward.

• both are very verbose (unlike the short-hand proj string).
  • WKT2 supported by OGC: https://docs.ogc.org/is/12-063r5/12-063r5.html
  • And yes I know that a WKT2 string is ugly (almost too verbose).
• WKT2 has good support in PROJ/GDAL/QGIS/...
• WKT2 I think is used in GeoPackage
• WKT2 can be used in LAZ/COPC lidar format (perhaps some limitations).
• While still in the works by OGC, GeoTIFF may be updated to support a WKT2 string. This will take a longtime to take over current GeoTIFF keywords, but shows what others are thinking for pushing WKT2 forward.
• WKT2 is the default for the IAU 2015 codes (if anything, these codes can be used for examples to be manually updated for custom WKT2 strings).
• Our STAC metadata is using WKT2, e.g.,: DTEPD_032625_1120_041144_1120_A01
• Lastly, you can see what a WKT2 string might look like by using PROJ or GDAL: gdalsrsinfo -o WKT2 +proj=stere +lon_0=31.6218 +lat_0=-85.42088 +R=1737400 +units=m
  • Note the many "unknown" description tags. Hence the issue with PROJ strings (and this conversion is also making a best guess for defining meter and degree lengths).

So I guess until JSON takes over the world, I would recommend using WKT2 when possible. The main complaint for WKT2 is that it has a somewhat odd formatting method.

[jlaura]

I'll second what Trent has posted. Ideally, this is using anything PROJ supported and data producers are using either short codes or WKT2 (or the wrapped WKT2 that is PROJJSON).

In the short term, I can think about how viable it is to debug the projection issues down stream when working with truncated proj strings. This is less of an issue when looking, for example, at orthos or intersection error images in a GIS. It becomes a bit trickier working with point cloud data that should have a proper vertical datum assigned.

All of that is to say, not a fire, but definitely something to consider, from my perspective, as important.

[dshean]

+1 for WKT2 as best current option. It's cumbersome, but the only way to unambiguously define 3D CRS. We can recommend people use WKT2 when possible, and point users to PROJ and other external doc on how to define CRS. But we also want to support legacy CRS definitions and the familiar EPSG codes. In practice, I think we can just replace any ASP utility options formerly requiring proj4 to be more generic srs, mirroring the GDAL syntax.

As @oleg-alexandrov mentioned, this is a larger task, requiring additional research and dev time to identify and update instances throughout the ASP/VW code base. But it will absolutely be worth the trouble in the long run, especially now that the bulk of the implementation is done in PROJ/GDAL.

Can refer to existing implementations in GDAL and PROJ for generic parsing of user-defined SRS input. A few examples from utilities that do this (I haven't dug into API for this functionality):

• https://github.com/OSGeo/gdal/blob/a3f8f3dddfefc516c3ba86d195196d28cb673690/apps/gdalsrsinfo.cpp#L315
• https://github.com/OSGeo/PROJ/blob/master/src/apps/projinfo.cpp

Beyond supporting unambiguous CRS definitions, at least on Earth, we also need to start specifying and tracking dataset epochs (e.g., stereo image acquisition date). And assuming we have a proper CRS definition (with epoch) for our ASP dataset, we can then rely on PROJ to correct for plate motion when the user specifies a point2dem/point2las output CRS with a fixed epoch like NAD83(2011) and the forthcoming updated NGS datums (https://geodesy.noaa.gov/datums/newdatums/index.shtml) for the US.

Most vendor camera position metadata will use the latest ITRF realization, and most users will specify a dynamic output CRS (like "WGS84", or hopefully, a specific WGS84 realization, like "WGS84 (G2139)"). This is fine, but we need to record the epoch of the observation so that these PC/DEM products can be combined with other datasets acquired at different epochs in the same dynamic CRS. For example, when combining a WV-1 DEM from 2007 with a WV-3 DEM of the same site from 2024, features on the ground will have moved >30-50 cm horizontally due to plate motion (over a pixel), depending on the location, even up to ~2 m for fastest plate motion of 10 cm/yr. See Figure 11 (https://agupubs.onlinelibrary.wiley.com/cms/asset/0224905c-ac15-4961-b2dc-b52da01bc8e6/jgrb51713-fig-0011-m.jpg) from https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2016JB013098. This may not sound like a big deal, but the geolocation errors introduced by ignoring these issues with current ASP CRS handling/assumptions greatly exceeds accuracy of most reference lidar and GCP datasets.

In short, we need improved CRS support anywhere we integrate external datasets in ASP (e.g., bundle_adjust, pc_align) and anywhere we export products with a CRS definition (e.g., point2dem, point2las, dem_geoid). I think next step is to identify instances throughout the code base where this needs attention, so we can assess scope and dev time required.

[oleg-alexandrov]

I made ASP accept any string passed via --t_srs. It goes directly to GDAL. I tested with Jay's Moon string, GeoJSON, and a couple of Earth strings. This required a big overhaul of our 10-year old projection handling logic.

This will be available in tomorrow's build.

Here is what pdal info will output:

"gtiff": "Geotiff_Information:\n Version: 1\n Key_Revision: 1.0\n Tagged_Information:\n End_Of_Tags.\n Keyed_Information:\n GTModelTypeGeoKey (Short,1): ModelTypeGeographic\n GTRasterTypeGeoKey (Short,1): RasterPixelIsArea\n GeographicTypeGeoKey (Short,1): User-Defined\n GeogCitationGeoKey (Ascii,122): \"GCS Name = Moon(2015)-Sphere/Ocentric|Datum = Moon(2015)-Sphere|Ellipsoid = Moon(2015)-Sphere|Primem = ReferenceMeridian|\"\n GeogGeodeticDatumGeoKey (Short,1): User-Defined\n GeogAngularUnitsGeoKey (Short,1): Angular_Degree\n GeogEllipsoidGeoKey (Short,1): User-Defined\n GeogSemiMajorAxisGeoKey (Double,1): 1737400 \n GeogSemiMinorAxisGeoKey (Double,1): 1737400 \n GeogPrimeMeridianLongGeoKey (Double,1): 0 \n End_Of_Keys.\n End_Of_Geotiff.\n",