Specifying exact reprojection with proj pipeline in gdalwarp

[bpurinton]

My approach

Let's say I have a tif in UTM Zone 3N (EPSG:32603) and I want to reproject into NAD83(2011) / Alaska Albers (EPSG:6393). It's a trivial example because I could easily enough just:

gdalwarp -t_srs EPSG:6393 source_epsg32603.tif dest_epsg6393.tif

But, let's say I was interested in the gory, technical details and wanted a PROJ pipeline to do the conversion:

$ projinfo -s EPSG:32603 -t EPSG:6393 -o PROJ --hide-ballpark --spatial-test intersects

Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of UTM zone 3N + Inverse of NAD83(2011) to WGS 84 (1) + Alaska Albers (meters), 2 m, Puerto Rico -
 onshore and offshore. United States (USA) onshore and offshore - Alabama; Alaska; Arizona; Arkansas; California; Colorado; 
Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; 
Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New 
Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South 
Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands - onshore 
and offshore.

PROJ string:
+proj=pipeline
  +step +inv +proj=utm +zone=3 +ellps=WGS84
  +step +proj=aea +lat_0=50 +lon_0=-154 +lat_1=55 +lat_2=65 +x_0=0 +y_0=0
        +ellps=GRS80

There's only one pipeline to choose from, and I would like to use that fully qualified pipeline (despite the 2 m uncertainty.. yikes) to do the reprojection.

I then run:

gdalwarp -ct "+proj=pipeline +step +inv +proj=utm +zone=3 +ellps=WGS84 +step +proj=aea +lat_0=50 +lon_0=-154 +lat_1=55 +lat_2=65 +x_0=0 +y_0=0 +ellps=GRS80" source_epsg32603.tif dest_epsg6393.tif

That works and produces the same result as the shorter gdalwarp comand with EPSG codes. But the benefit of this second approach is that I can (a) see the uncertainty associated with the conversion and potentially factor that into an analysis; (b) in the case of multiple pipeline options, I could actually choose one to pass to -ct in gdalwarp

My Question

I'm curious: is this a reasonable thing to do? Or should I just trust EPSG codes and do my conversions the way I always have in the past with the first gdalwarp command?

Sidenote

It appears that the second gdalwarp command with the -ct pipeline option does not write the CRS information into the GeoTiff metadata. A gdalinfo on dest_epsg6393.tif returns the old UTM 3N CRS info. Possible bug, opened an issue: OSGeo/gdal#8242

Edit

This CRS issue is resolved on OSGeo/gdal#8243. Turns out you still need to use the t_srs option to write out the metadata even though the parameters of the proj pipeline will be respected and override the default params, so the updated, second gdalwarp command should have been:

gdalwarp -ct "+proj=pipeline +step +inv +proj=utm +zone=3 +ellps=WGS84 +step +proj=aea +lat_0=50 +lon_0=-154 +lat_1=55 +lat_2=65 +x_0=0 +y_0=0 +ellps=GRS80" -t_srs EPSG:6393 source_epsg32603.tif dest_epsg6393.tif

Now the commands are really redundant, but at least you gain access control over the exact transformation parameters when using the ct flag....

[dshean]

Thanks for detailed question!

The short answer is that explicitly defining your PROJ pipeline is the safest approach, but it can be cumbersome.

If your t_srs and s_srs definitions are unambiguous, then your results should be identical to those obtained with the pipeline. If your EPSG codes unambiguously define your CRS, then you can use those for convenience.

But in this case, your EPSG codes are only 2D CRS, with no definition of specific ITRF/WGS84 realizations. You can explicitly define your 3D CRS using well-known text (WKT2), adding explicit information about which WGS84 you're using with your UTM projection.

The 2m uncertainty is due to the ambiguity of the WGS84 realization with the EPSG:326XX codes. See the output for the ENSEMBLE.

gdalsrsinfo EPSG:32603

PROJ.4 : +proj=utm +zone=3 +datum=WGS84 +units=m +no_defs

OGC WKT2:2019 :
PROJCRS["WGS 84 / UTM zone 3N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 3N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-165,
            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["Navigation and medium accuracy spatial referencing."],
        AREA["Between 168°W and 162°W, northern hemisphere between equator and 84°N, onshore and offshore. United States (USA) - Alaska (AK)."],
        BBOX[0,-168,84,-162]],
    ID["EPSG",32603]]

Note the ENSEMBLEACCURACY. The way around this is to explicitly define the WGS84/ITRF realization of your dataset. I don't think there are existing EPSG codes for UTM projection with specific realiztion (like "World Geodetic System 1984 (G1762)"), but you can do this with WKT2 (can save as text file and pass to t_srs). That should significantly reduce the uncertainty of the transformation.

Hopefully that makes sense. There are some good discussions of the WGS84 Ensemble issues in the GDAL and PROJ mailing list archives.

[dshean]

Need to add some text on this WGS84 Ensemble issue to the notes: https://github.com/ICESAT-2HackWeek/3D_CRS_Transformation_Resources/issues/2

[bpurinton]

Note the ENSEMBLEACCURACY. The way around this is to explicitly define the WGS84/ITRF realization of your dataset. I don't think there are existing EPSG codes for UTM projection with specific realiztion (like "World Geodetic System 1984 (G1762)"), but you can do this with WKT2 (can save as text file and pass to t_srs). That should significantly reduce the uncertainty of the transformation.

Woah! Not the answer I was expecting here. Thanks for the detailed response. WGS84 is not WGS84 indeed...