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

[Even Rouault]

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
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20220608/cf939ca8/attachment.htm>