uw-cryo / 3d_crs_transformation_resources Goto Github PK

projinfo EPSG:4326
ensemble
2.0 m uncertainty
What happens when you have a clear definition for one dataset (ITRF2014) combined with one in WGS84 ensemble?
Can introduce an apparent offset!

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....

check versions!
PROJ_NETWORK = 'ON'
horizontal coordinate system (lat/lon) in a WGS84 realization for a WGS84 ellipsoid in a ITRF realization
WGS w/ GPS week definitions
Local plane vs. geodetic surveying
Definitions for center of coordinate system (mass, shape)
EPSG db updated with latest PROJ
projsync --bbox
WGS84 ensemble description
Make sure all epsg codes have links: https://epsg.io/7912
geopandas 0.13
projinfo bbox