[gdal-dev] Discrepancy when tiling on MAC and Linux

classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Discrepancy when tiling on MAC and Linux

Grégory Bataille
Hey all,

So I have a strange issue that I somewhat isolated, but I don't know how to move further.
I have dataset in EPSG 31468-15949

My goal is to tile it, render it with openlayers, and add a marker on the GCP location.
Trouble is that the GCP location was falling a bit away from the marker image.

I had a hard time investigating that, and I still do but I realized that tiling in other conditions did work perfectly.

So in Production
- I have a docker image which is basically osgeo/gdal:alpine-normal-3.0.2
  - this means GDAL 3.0.2
  - PROJ 6.2.1

If I use this docker to do the tiling, then I get the "badly aligned" tile.

On my local MAC
- I have gdal 3.0.1 installed from the brew repository osgeo/osgeo4mac/osgeo-gdal
- I have proj 6.2.0
- I have the GDAL 3.0.1 pip package installed

If I use my MAC to tile, I get the "aligned" tile

And that's it, I'm stuck.
I tried in a python shell to do coordinates transform of the GCP from the source WKT to EPSG 3857 and I get the same value whether in the docker or on my MAC

Any idea how I can move further?

I have put a zip with some data here

Thanks

---
Gregory Bataille

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Even Rouault-2
On vendredi 8 novembre 2019 08:18:11 CET Grégory Bataille wrote:
> Hey all,
>
> So I have a strange issue that I somewhat isolated, but I don't know how to
> move further.
> I have dataset in EPSG 31468-15949

I would check if that wouldn't be related to availability of the BETA2007.gsb
grid. osgeo/gdal:alpine-normal-3.0.2 has it. Does your MAC install have it ?

That said, your TIFF file has a hardcoded TOWGS84=598.1,73.7,418.2,0.202          
,0.045,-2.455,6.7 . When using

gdalwarp Investigation_tile_alignment/ClippedProject.tif out.tif -overwrite -
t_srs EPSG:4326 --debug on

I see locally that BETA2007.gsb is used. This is not so bad, but a bit
unexpected since I'd thought that forced TOWGS84 would have take the
precedence.

Checking further in ogrct.cpp, I see that as the GeoTIFF has also the EPSG:
31468 code, and when instanciating the PROJ coordinate transformation, this is
the preferred initializer, and thus the hardcoded TOWGS84 is ignored. So this
explains what I observe. I'm not completely clear of what would be the desired
behaviour here. In some use cases, one might prefer to honour the possibly
suboptimal TOWGS84 from the geotiff info (the simple fact of doing
gdal_translate -a_srs EPSG:31468 hardcodes it. Even with GDAL 3, which is
probably no longer a good idea...). In other cases, just trust PROJ to use the
best transformation method.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Grégory Bataille
That’s a problem for us. With gdal 3 we started to force values for the towgs84 so that proj would not decide.
This is wanted because our front end will place markers and it does coordinates transform. Without access to proj 6. Meaning it needs to have the transform Value.

I’ll check what you mentioned (on Monday though I think)

Thanks for having a look

Greg

On Fri, 8 Nov 2019 at 12:30, Even Rouault <[hidden email]> wrote:
On vendredi 8 novembre 2019 08:18:11 CET Grégory Bataille wrote:
> Hey all,
>
> So I have a strange issue that I somewhat isolated, but I don't know how to
> move further.
> I have dataset in EPSG 31468-15949

I would check if that wouldn't be related to availability of the BETA2007.gsb
grid. osgeo/gdal:alpine-normal-3.0.2 has it. Does your MAC install have it ?

That said, your TIFF file has a hardcoded TOWGS84=598.1,73.7,418.2,0.202           
,0.045,-2.455,6.7 . When using

gdalwarp Investigation_tile_alignment/ClippedProject.tif out.tif -overwrite -
t_srs EPSG:4326 --debug on

I see locally that BETA2007.gsb is used. This is not so bad, but a bit
unexpected since I'd thought that forced TOWGS84 would have take the
precedence.

Checking further in ogrct.cpp, I see that as the GeoTIFF has also the EPSG:
31468 code, and when instanciating the PROJ coordinate transformation, this is
the preferred initializer, and thus the hardcoded TOWGS84 is ignored. So this
explains what I observe. I'm not completely clear of what would be the desired
behaviour here. In some use cases, one might prefer to honour the possibly
suboptimal TOWGS84 from the geotiff info (the simple fact of doing
gdal_translate -a_srs EPSG:31468 hardcodes it. Even with GDAL 3, which is
probably no longer a good idea...). In other cases, just trust PROJ to use the
best transformation method.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
--
Greg

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Grégory Bataille
Hey Even,

So I can confirm that this is the issue
On Mac
OGRCT: Selecting transformation +proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 (Inverse of 3-degree Gauss-Kruger zone 4 + DHDN to WGS 84 (2) + Inverse of Transformation to WGS84)

On my linux docker
OGRCT: Selecting transformation +proj=pipeline +step +proj=axisswap +order=2,1 +step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +step +proj=hgridshift +grids=BETA2007.gsb +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=598.1 +y=73.7 +z=418.2 +rx=0.202 +ry=0.045 +rz=-2.455 +s=6.7 +convention=position_vector +step +inv +proj=cart +ellps=bessel +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1 (Inverse of 3-degree Gauss-Kruger zone 4 + DHDN to WGS 84 (4) + Inverse of Transformation to WGS84)

Now this is a problem to me and I don't quite know what to do with it.

So I have a fully specified dataset, that gives a TOWGS84, that my user has chosen (properly or not)
Once the dataset is using this representation, I can't have auto projection because those will/might depend on the application using the dataset.
Currently (most of) my backend processes are in Gdal3 / Proj6 and therefore (per my understanding), the TOWGS84 might be overriden
However my front end processes only have some proj4 equivalence and will trust exactly the dataset.

I'm a bit at a loss to know what I should do.
Should I do a gdalwarp step for gdal to "chose" the best projection and get a new dataset (with another TOWGS84) that I use internally for alignment, while giving the client the one he has asked for when he wants to get the raw dataset?

Do you have any ideas/recommandations? Is there a way to disable the PROJ6 auto-selection feature?

Thanks

---
Gregory Bataille


On Fri, Nov 8, 2019 at 12:36 PM Grégory Bataille <[hidden email]> wrote:
That’s a problem for us. With gdal 3 we started to force values for the towgs84 so that proj would not decide.
This is wanted because our front end will place markers and it does coordinates transform. Without access to proj 6. Meaning it needs to have the transform Value.

I’ll check what you mentioned (on Monday though I think)

Thanks for having a look

Greg

On Fri, 8 Nov 2019 at 12:30, Even Rouault <[hidden email]> wrote:
On vendredi 8 novembre 2019 08:18:11 CET Grégory Bataille wrote:
> Hey all,
>
> So I have a strange issue that I somewhat isolated, but I don't know how to
> move further.
> I have dataset in EPSG 31468-15949

I would check if that wouldn't be related to availability of the BETA2007.gsb
grid. osgeo/gdal:alpine-normal-3.0.2 has it. Does your MAC install have it ?

That said, your TIFF file has a hardcoded TOWGS84=598.1,73.7,418.2,0.202           
,0.045,-2.455,6.7 . When using

gdalwarp Investigation_tile_alignment/ClippedProject.tif out.tif -overwrite -
t_srs EPSG:4326 --debug on

I see locally that BETA2007.gsb is used. This is not so bad, but a bit
unexpected since I'd thought that forced TOWGS84 would have take the
precedence.

Checking further in ogrct.cpp, I see that as the GeoTIFF has also the EPSG:
31468 code, and when instanciating the PROJ coordinate transformation, this is
the preferred initializer, and thus the hardcoded TOWGS84 is ignored. So this
explains what I observe. I'm not completely clear of what would be the desired
behaviour here. In some use cases, one might prefer to honour the possibly
suboptimal TOWGS84 from the geotiff info (the simple fact of doing
gdal_translate -a_srs EPSG:31468 hardcodes it. Even with GDAL 3, which is
probably no longer a good idea...). In other cases, just trust PROJ to use the
best transformation method.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
--
Greg

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Even Rouault-2
On mercredi 13 novembre 2019 13:18:52 CET Grégory Bataille wrote:

> Hey Even,
>
> So I can confirm that this is the issue
> On Mac
> OGRCT: Selecting transformation +proj=pipeline +step +proj=axisswap
> +order=2,1 +step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000
> +y_0=0 +ellps=bessel +step +proj=unitconvert +xy_in=rad +xy_out=deg +step
> +proj=axisswap +order=2,1 (Inverse of 3-degree Gauss-Kruger zone 4 + DHDN
> to WGS 84 (2) + Inverse of Transformation to WGS84)
>
> On my linux docker
> OGRCT: Selecting transformation +proj=pipeline +step +proj=axisswap
> +order=2,1 +step +inv +proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000
> +y_0=0 +ellps=bessel +step +proj=hgridshift +grids=BETA2007.gsb +step
> +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert
> +x=598.1 +y=73.7 +z=418.2 +rx=0.202 +ry=0.045 +rz=-2.455 +s=6.7
> +convention=position_vector +step +inv +proj=cart +ellps=bessel +step
> +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step
> +proj=axisswap +order=2,1 (Inverse of 3-degree Gauss-Kruger zone 4 + DHDN
> to WGS 84 (4) + Inverse of Transformation to WGS84)
>
> Now this is a problem to me and I don't quite know what to do with it.
>
> So I have a fully specified dataset, that gives a TOWGS84, that my user has
> chosen (properly or not)
> Once the dataset is using this representation, I can't have auto projection
> because those will/might depend on the application using the dataset.
> Currently (most of) my backend processes are in Gdal3 / Proj6 and therefore
> (per my understanding), the TOWGS84 might be overriden
> However my front end processes only have some proj4 equivalence and will
> trust exactly the dataset.
>
> I'm a bit at a loss to know what I should do.
> Should I do a gdalwarp step for gdal to "chose" the best projection and get
> a new dataset (with another TOWGS84) that I use internally for alignment,
> while giving the client the one he has asked for when he wants to get the
> raw dataset?
>
> Do you have any ideas/recommandations? Is there a way to disable the PROJ6
> auto-selection feature?

This is not a PROJ6 issue, but a GDAL one related to the logic starting at
https://github.com/OSGeo/gdal/blob/master/gdal/ogr/ogrct.cpp#L843

You could try to add a test that if OGRSpatialReference::GetTOWGS84() returns
OGRERR_NONE then the resolution to the WKT of the EPSG code is skipped

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Even Rouault-2
> This is not a PROJ6 issue, but a GDAL one related to the logic starting at
> https://github.com/OSGeo/gdal/blob/master/gdal/ogr/ogrct.cpp#L843
>
> You could try to add a test that if OGRSpatialReference::GetTOWGS84()
> returns OGRERR_NONE then the resolution to the WKT of the EPSG code is
> skipped

I've implemented this in https://github.com/OSGeo/gdal/pull/2016

The last commit is the one that is of interest for you. I've kept the current
behaviour as default, as I think it is what most people would want. For your
use case, you'd define OSR_CT_USE_DEFAULT_EPSG_TOWGS84=YES

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Grégory Bataille
wow, thanks.

I'm not quite sure I'm following it all though. I'm reaching a bit my knowledge limit here with those transforms. I'll need to spend a bit more time on it to get more understanding and see what I can do, how I might use your update.
But thanks for the quick turnaround

Cheers

---
Gregory Bataille


On Thu, Nov 14, 2019 at 1:24 AM Even Rouault <[hidden email]> wrote:
> This is not a PROJ6 issue, but a GDAL one related to the logic starting at
> https://github.com/OSGeo/gdal/blob/master/gdal/ogr/ogrct.cpp#L843
>
> You could try to add a test that if OGRSpatialReference::GetTOWGS84()
> returns OGRERR_NONE then the resolution to the WKT of the EPSG code is
> skipped

I've implemented this in https://github.com/OSGeo/gdal/pull/2016

The last commit is the one that is of interest for you. I've kept the current
behaviour as default, as I think it is what most people would want. For your
use case, you'd define OSR_CT_USE_DEFAULT_EPSG_TOWGS84=YES

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Even Rouault-2
On jeudi 14 novembre 2019 05:54:09 CET Grégory Bataille wrote:
> wow, thanks.
>
> I'm not quite sure I'm following it all though. I'm reaching a bit my
> knowledge limit here with those transforms. I'll need to spend a bit more
> time on it to get more understanding and see what I can do, how I might use
> your update.
> But thanks for the quick turnaround

Another option if you want full control on the transformation would be to
provide the PROJ pipeline directly. gdalwarp has the -ct switch for that,
which ultimately uses OGRCoordinateTransformation::SetCoordinateOperation()

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Grégory Bataille
Oh, I could gdalwarp the dataset myself to WGS84 and then tile it, so as to consistently "ignore" grids.
That's an idea indeed.

Thanks for the pointer, I'll take it into consideration

---
Gregory Bataille


On Thu, Nov 14, 2019 at 12:46 PM Even Rouault <[hidden email]> wrote:
On jeudi 14 novembre 2019 05:54:09 CET Grégory Bataille wrote:
> wow, thanks.
>
> I'm not quite sure I'm following it all though. I'm reaching a bit my
> knowledge limit here with those transforms. I'll need to spend a bit more
> time on it to get more understanding and see what I can do, how I might use
> your update.
> But thanks for the quick turnaround

Another option if you want full control on the transformation would be to
provide the PROJ pipeline directly. gdalwarp has the -ct switch for that,
which ultimately uses OGRCoordinateTransformation::SetCoordinateOperation()

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Discrepancy when tiling on MAC and Linux

Grégory Bataille
Hey Even,

Actually a "simpler" solution for now for us is to remove the grid from the docker image we use to run dataset tiling. This way, the projection follows TOWGS84.
It's suboptimal in term of absolute accuracy of the tiles, but it gives us consistency in our data which is more important to us.

Thanks for the support. I don't think I would have found this grid thing on my own.

---
Gregory Bataille


On Thu, Nov 14, 2019 at 1:03 PM Grégory Bataille <[hidden email]> wrote:
Oh, I could gdalwarp the dataset myself to WGS84 and then tile it, so as to consistently "ignore" grids.
That's an idea indeed.

Thanks for the pointer, I'll take it into consideration

---
Gregory Bataille


On Thu, Nov 14, 2019 at 12:46 PM Even Rouault <[hidden email]> wrote:
On jeudi 14 novembre 2019 05:54:09 CET Grégory Bataille wrote:
> wow, thanks.
>
> I'm not quite sure I'm following it all though. I'm reaching a bit my
> knowledge limit here with those transforms. I'll need to spend a bit more
> time on it to get more understanding and see what I can do, how I might use
> your update.
> But thanks for the quick turnaround

Another option if you want full control on the transformation would be to
provide the PROJ pipeline directly. gdalwarp has the -ct switch for that,
which ultimately uses OGRCoordinateTransformation::SetCoordinateOperation()

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev