[gdal-dev] Transforming from EPSG:29902 to EPSG:4326 at the limit of the area of use.

[Sam Hall]

I tried using both:

gdalwarp -ct "$(projinfo -s EPSG:29902 -t EPSG:4326 -o proj -q)" input.tif output.tif
gdalwarp -ct "$(projinfo -s EPSG:29902 -t EPSG:4326 -o proj -q --bbox -10.562,52.08,-10.561,52.1)" input.tif output.tif

and they both produce the desired result of a smooth transition across the area of use.

Thanks so much for your help.

Sam

From: Sam Hall <Sam.Hall at jbarisk.com>
Sent: Wednesday, November 8, 2023 14:57
To: gdal-dev at lists.osgeo.org <gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Transforming from EPSG:29902 to EPSG:4326 at the limit of the area of use.
 
I didn't mean to just respond to you Javier, sorry about that.

That looks promising. I'll work out which is the one is using which transformation and see if I can use the right one.

Sam

From: Javier Jimenez Shaw <j1 at jimenezshaw.com>
Sent: Wednesday, November 8, 2023 14:39
To: Sam Hall <Sam.Hall at jbarisk.com>
Subject: Re: [gdal-dev] Transforming from EPSG:29902 to EPSG:4326 at the limit of the area of use.

projinfo -s EPSG:29902 -t EPSG:4326 -o proj
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of Irish Grid + TM65 to WGS 84 (2), 1 m, Ireland - onshore. United Kingdom (UK) - Northern Ireland (Ulster) - onshore.

PROJ string:
+proj=pipeline
  +step +inv +proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000
        +y_0=250000 +a=6377340.189 +rf=299.3249646
  +step +inv +proj=longlat +a=6377340.189 +rf=299.3249646
  +step +proj=push +v_3
  +step +proj=cart +a=6377340.189 +rf=299.3249646
  +step +proj=helmert +x=482.5 +y=-130.6 +z=564.6 +rx=-1.042 +ry=-0.214
        +rz=-0.631 +s=8.15 +convention=position_vector
  +step +inv +proj=cart +ellps=WGS84
  +step +proj=pop +v_3
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1



projinfo -s EPSG:29902 -t EPSG:4326 -o proj --bbox "-10.562,52.08,-10.561,52.1"
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of Irish Grid + Ballpark geographic offset from TM65 to WGS 84, unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=pipeline
  +step +inv +proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000
        +y_0=250000 +a=6377340.189 +rf=299.3249646
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

So I guess that somehow one is taking the datum transformation and the other doesn't. Which one is doing what you should check.

BTW, you are writing only to me, not to the mailing list. Are you aware of that?

Best,
Javier.

On Wed, 8 Nov 2023 at 15:23, Sam Hall <Sam.Hall at jbarisk.com> wrote:
Thanks Javier,

It doesn't happen on the Irish mainland.

Is it possible to define the transformation to use outside the area of usage of the datum transformation using a proj string with something like this "+proj=hgridshift +grids=first,second" in it? Or are you suggesting that gdalwarp doesn't take the area of usage into account and therefore wouldn't look for the second datum shift and use a ballpark transform.

Sam


From: Javier Jimenez Shaw <j1 at jimenezshaw.com>
Sent: Wednesday, November 8, 2023 13:47
To: Sam Hall <Sam.Hall at jbarisk.com>
Cc: gdal-dev at lists.osgeo.org <gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Transforming from EPSG:29902 to EPSG:4326 at the limit of the area of use.

Assuming that it does not happen in the main area of Ireland...I guess it is because it is outside of the area of usage of the datum transformation https://epsg.org/transformation_1641/TM65-to-WGS-84-2.html
Probably (this is just a guess) gdalwarp is not taking it into consideration, making a ballpark transformation between TM65 and WGS84, while ogr2ogr doesn't.

On Wed, 8 Nov 2023 at 13:35, Sam Hall via gdal-dev <gdal-dev at lists.osgeo.org> wrote:
I have some data that covers the whole of the Republic of Ireland, including the Blasket Islands off the West Coast, situated around 24439, 95540 in EPSG:29902. These islands straddle the western boundary of the area of use of the coordinate system -10.56° W.

When transforming a raster dataset from EPSG:29902 to EPSG:4326 using gdalwarp at gdal >= 3.5.2, there seems to be a roughly 50 m shift in the data at the area of use boundary, while if I use ogr2ogr on a vector copy of this data, I don't see this shift.

I have a reproducible example of this at https://gitlab.com/Sam.Hall1/29902_gdalwarp_vs_ogr2ogr.git. In this example I created a rectangular polygon with interpolated points and saved it to a shapefile, then used gdal_rasterize to create a GeoTiff. gdalwarp gives the desired result using a gdal 3.5.1 docker image. The full commands are in the docker-compose.yml in the repository but the parameter's I'm using are:

gdalwarp -t_srs EPSG:4326 -tr 0.000045 0.000045 -tap -overwrite input.tif output.tif
ogr2ogr -t_srs EPSG:4326 output.shp input.shp

I've also seen similar results when transforming from EPSG:27700 to EPSG:4326 in Jersey and Guernsey in the English Channel but I'd have to use a much earlier version of gdal to reproduce the desired result.

Is there a gridshift file or any additional parameters that I could use to get the output to look more continuous across the area of use line?

Sam
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev