[gdal-dev] Issue defining compound projection with a local projection for GeoTIFF

[Chris Smemoe]

I'm not sure whether this is a bug or some issue with how I'm doing things
(like maybe my WKT is bogus), but when I try to define a compound projection
for a GeoTIFF with a local projection, the vertical projection information
does not get saved to the GeoTIFF and I get an error reported from PROJ.  I
am trying to upgrade the version of GDAL I'm using and I think this used to
work in GDAL version 2.x but no longer works in GDAL 3.x/PROJ 7.x.  What I'm
trying to do is to update the projection of a GeoTIFF from a local
projection with horizontal units only to a compound projection with both
horizontal and vertical units.  When I close the GDAL dataset, It seems like
PROJ is trying to look for an ellipsoid node for the spatial reference
object but I don't have one since I'm using a local projection.  Then when I
reopen the dataset, nothing seems to have been saved.  I'm running GDAL
3.4.2.  Here's an example of what I'm trying to do:

 

-------------------begin-------------------

python

>>> from osgeo import gdal

>>> gd = gdal.Open("C:\\temp\\grid_meter_exported.tif", gdal.GA_Update)

>>> gd

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at
0x000001B6BC578630> >

>>> gd.GetProjectionRef()

'LOCAL_CS["None",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAS
T],AXIS["Northing",NORTH]]'

>>>
gd.SetProjection('COMPD_CS["None",LOCAL_CS["None",LOCAL_DATUM["None",32767],
UNIT["Foot_US",0.304800609601219],AXIS["None",OTHER]],VERT_CS["Local
Meter",VERT_DATUM["Local",2005],UNIT["Meter",1,AUTHORITY["EPSG","9001"]],AXI
S["Gravity-related height",UP]]]')

0

>>> gd.GetProjection()

'COMPD_CS["None",LOCAL_CS["None",LOCAL_DATUM["None",32767],UNIT["Foot_US",0.
304800609601219],AXIS["None",OTHER]],VERT_CS["Local
Meter",VERT_DATUM["Local",2005],UNIT["Meter",1,AUTHORITY["EPSG","9001"]],AXI
S["Gravity-related height",UP]]]'

>>> gd = None

ERROR 1: PROJ: proj_get_ellipsoid: CRS has no geodetic CRS

ERROR 1: PROJ: proj_get_ellipsoid: Object is not a CRS or
GeodeticReferenceFrame

ERROR 1: PROJ: proj_get_ellipsoid: CRS has no geodetic CRS

ERROR 1: PROJ: proj_get_ellipsoid: Object is not a CRS or
GeodeticReferenceFrame

ERROR 1: PROJ: proj_crs_get_geodetic_crs: CRS has no geodetic CRS

>>> gd = gdal.Open("C:\\temp\\grid_meter_exported.tif", gdal.GA_Update)

>>> gd

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at
0x000001B6BC578720> >

>>> gd.GetProjectionRef()

'LOCAL_CS["None",UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAS
T],AXIS["Northing",NORTH]]'

 

-------------------end-------------------

 

The dataset I'm working with can be downloaded here:

 

https://www.dropbox.com/s/2jn7vkcl690e90r/grid_meter_exported.tif?dl=1

 

I'd appreciate any guidance.

 

Thanks,

 

Chris

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20220608/d02615b0/attachment.htm>