0

I am unable to define a local projection with no units for a GeoTIFF file using GDAL. Here is a sample of what I'm trying to do in Python:

python
>>> from osgeo import gdal
>>> gd = gdal.Open("C:\\temp\\57898_local_units.tif", gdal.GA_Update)
>>> gd
<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x00000223C4818690> >
>>> gd.SetProjection('LOCAL_CS["None",LOCAL_DATUM["None",32767],UNIT["None",0],AXIS["None",OTHER]]')
0
>>> gd.GetProjectionRef()
'LOCAL_CS["None",LOCAL_DATUM["None",32767],UNIT["None",0],AXIS["None",OTHER]]'
>>> gd = None
>>> gd = gdal.Open("C:\\temp\\57898_local_units.tif", gdal.GA_Update)
>>> gd.GetProjectionRef()
'LOCAL_CS["None",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'

The problem is that when I close and re-open the dataset the WKT for the dataset changes from having no units with a conversion factor of 0 to having metric units with a conversion factor of 1. I am using GDAL 3.4.1 with python, but I have similar issues with newer versions of GDAL in C++. In older versions of GDAL (version 2.3.1 maybe), I used to be able to get a local projection with no units from the projection string described above, though admittedly the projection string was not exactly the same as the input. In those versions, I got the following for the projection string after updating the raster data with the same input projection as above:

LOCAL_CS["unnamed",UNIT["unknown",0]]

Is there any way to fix this issue of not keeping a local projection with no units after exporting it to the GeoTIFF or is this something that is no longer supported in GDAL?

1 Answer 1

2

Your issue may be connected to this GDAL commit https://github.com/OSGeo/PROJ/commit/fbc654ab0c3dc6d707bdee226121275736c2bfb5

EngineeringCRS: when exporting to WKT1_GDAL, output unit and axis (fixes OSGeo/gdal#2347)

The bug it fixed is https://github.com/OSGeo/gdal/issues/2347, Unit of Local CS is removed when doing Dataset.GetProjection()

There seems to be a test case for a custom unit createVerticalCRS()->alterCSLinearUnit(UnitOfMeasure("my unit", 2)); but maybe the changed behavior comes from the GeoTIFF driver.

Maybe you should write mail to gdal-dev mailing list and ask directly from the developers. Tell also why do you want to use an undefined unit that correspond to 0 meters in that way.

2
  • I'll send a message to the GDAL dev mailing list or just report a bug on github. I think I have a workaround where I can use some value besides 0 for the metric unit conversion value and GDAL appears to save the WKT correctly. We're just trying to save out to the GeoTIFF that no units are defined so we can load this information back into our software. Thanks.
    – Burrow
    Commented Jun 2, 2022 at 15:54
  • I do not know very well the OGC CRS standard and GeoTIFF specifications. It will be interesting to see it there is a bug in GDAL now or if your software is made to depend on a bug that is now corrected.
    – user30184
    Commented Jun 2, 2022 at 16:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.