gdal2tiles.py: add capability to specify coordinate transformation (-ct switch)

[ghemann]

gdal2tiles.py generates tiles that can be displayed on various map engines, and many coordinate frames work, but there is a shift for all CRS's that are in the NAD83(HARN) coordinate frame. The image attached is the outline of a building in EPSG 2874 in the correct location, but the tiles of the image are shifted north east. This has been verified on 8 other EPSG all in NAD83(HARN).

All of my examples are in the north eastern quadrant, so they always shift in this direction, but vary based on true distance from 0,0.

gdal 3.2.2 and proj 8.0.0, and rendered on mapbox.

[rouault]

I doubt this is a GDAL bug. This might rather be a lack of PROJ grids (like the ones in https://github.com/OSGeo/PROJ-data/tree/master/us_noaa but they are only for NAD83 to NAD83(HARN)), or perhaps a lack in PROJ database of information to do datum transformation from NAD83(HARN) to your target datum (which you didn't specify), or the data on which you overlay isn't correctly referenced. In any case, this ticket is not actionable with just the information you provided.

[ghemann]

Sorry for my lack of knowledge here, trying to wrap my head around where the core problem is.

I assume that gdal2tiles.py always converts to 4326, is that not the case? At least I dont see any arguments to specify the target. When I transform that green square overlay from the original EPSG 2874 to 4326 using postgis, it converts to the correct location. And both that overlay and the original image in other software line up fine.

Are there ways I can validate its not a gdal issue?

[rouault]

I assume that gdal2tiles.py always converts to 4326, is that not the case?

not necessarily, but in the default mode that uses WebMercator, yes, it is WGS 84.

Are there ways I can validate its not a gdal issue?

You could use "gdaltransform -s_srs EPSG:2874 -t_srs EPSG:4326" with the coordinates of your rectangle.

But actually using the PROJ projinfo utility, I see that "projinfo -s EPSG:2874 -t EPSG:4326 --spatial-test intersects" will use a "NAD83(HARN) to WGS 84 (3)" transformation that introduces a shift between NAD83(2011) and EPSG:4326; There's also another possible transformation that use "NAD83(HARN) to WGS 84 (1)" transformation, which is a no-op transformation. Both transformations have a 1 metre accuracy, so PROJ has to select between them. It selects "NAD83(HARN) to WGS 84 (3)" because its area of use is CONUS, whereas the one for "NAD83(HARN) to WGS 84 (1)" also include off-shore areas. I suspect in your use case the basemap has been created with the assumption that NAD83(2011) and EPSG:4326 are identical (that is using "NAD83(HARN) to WGS 84 (1)")

As a workaround you could probably add a "-s EPSG:3498" option to gdal2tiles to force the use of "NAD83(NSRS2007) / California zone 5 (ftUS)" instead, since the transformation between NAD83(NSRS2007) and WGS 84 is a null one (at least currently with the geodetic database used by PROJ).

This is clearly a hack though. gdal2tiles.py should perhaps offer an option to specify which precise coordinate transformation is needed.

[ghemann]

@rouault Wow, you're suggestion to force it to the other source worked! So does this mean there is a shift in the proj database for all HARN and whenever that exists I should use the NSRS2007 instead? Sorry, didn't entirely follow all the projections you mentioned. I guess I'm just curious if this is an issue because (a) HARN just has a shift to wgs84 in proj (b) im interpretting the wrong crs when looking at the input image (c) or something else entirely. Thanks for the help!

[rouault]

So does this mean there is a shift in the proj database for all HARN and whenever that exists I should use the NSRS2007 instead? Sorry, didn't entirely follow all the projections you mentioned. I guess I'm just curious if this is an issue because (a) HARN just has a shift to wgs84 in proj (b) im interpretting the wrong crs when looking at the input image (c) or something else entirely. Thanks for the help!

This is admitedly a tricky subject if you're not familiar with geodesy concepts. But basically when transforming between 2 datums (HARN and WGS 84 here), there is not always a single universal way of doing that. The PROJ database has 2 potential entries for that: one that consider both datums equivalent with an accuracy of 1m, and another one that applies a shift, still with an accuracy of 1m. Which one to select highly depends on the context. And using WGS 84 itself might add an inaccuracy of 2 meters for various reasons... The trick I used is just based on the fact that I observed that NSRS2007 had just one entry registered in the database for a transformation from it to WGS84 which was a null transformation, and as a null transformation was apparently what was used for your basemap, the hack does the trick.

[reyemtm]

This depends on whether or not you want a null transformation from NAD83<>WGS84 or not. I have been using the following method to implement a virtual transformation for our imagery (also in NAD83 HARN) into WGS84 using the same calculations that Arc uses. The latest versions of the Arc software (Pro) use the ITRF00 transformations by default and it is very difficult to get it to use a null transformation. Also, all of our GPS data is transformed using this transformation. So, assuming your data is in either HARN or just NAD83 (so there very little shift between these), you can create a virtual raster into EPSG:3857 (not 4326 as this will lead to blurry tiles) then create your tiles from this virtual raster.

The key is to provide the desired shift to WGS84 in the -s_srs option since the default shift in proj is NULL. I am still getting a very slight shift when creating my tiles as compared to how it looks in QGIS and how they are exported if I transform the actual raster in Arc* into Web Mercator vs using the virtual raster method, but for these building outlines it would not be noticeable.

gdalwarp -t_srs EPSG:3857 -s_srs "+proj=lcc +lat_0=38 +lon_0=-82.5 +lat_1=40.0333333333333 +lat_2=38.7333333333333 +x_0=600000 +y_0=0 +ellps=GRS80 +towgs84=-0.9956,1.9013,0.5215,0.025915,0.009246,0.011599,-0.00062 +units=us-ft +no_defs" -srcnodata 255 -of VRT INPUT_TIFF.tif OUTPUT_TIFF.vrt

If you really want to do a deep dive into all the various projection transformations supported by Arc*, take a look here. It provides all the coordinate frame values for their transformations so you can use those in the gdal command above.