Chris,
Please file a ticket about that at https://github.com/OSGeo/gdal/issues/new
Even
Le 08/06/2022 à 09:11, Chris Smemoe via gdal-dev a écrit :
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",EAST],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"]],AXIS["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"]],AXIS["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",EAST],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
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
