[Geotiff] geos projection sweep parameter

[Even Rouault]

On jeudi 6 juillet 2017 17:12:44 CEST David Hoese wrote:
> Hi,
> 
> I work with data from the GOES-16 satellite ABI instrument. The official
> ABI Level 1B product consists of data mapped to the GEOS projection with
> a sweep angle axis parameter set to 'x' instead of the default 'y'. More
> information on this projection can be found here:
> http://proj4.org/projections/geos.html
> 
> As of GDAL v2.1+ the "sweep" parameter is handled internally and I can't
> remember what version of PROJ.4 added support for this. As far as I can
> tell there is no standard way to store the "sweep" parameter in a
> GeoTIFF file. I'm not sure if this is due to a limitation of the WKT
> specification or if it is something that only requires changes in the
> GeoTIFF library. Does anyone know what needs to happen to get "sweep"
> added to a GeoTIFF? If this requires changes to libgeotiff I am willing
> to add them or provide patches to add the functionality. Thanks for any
> information you can provide.
> 

David,

The issue is more fundamental than the sweep parameter being stripped off.
The main issue is that the GEOS projection does not exist at all in GeoTIFF representation
(GeoTIFF encoding is not nominally based on WKT, but consists of keys with codified values)
The list of  GeoTIFFprojection methods is :
http://web.archive.org/web/20160407200550/http://www.remotesensing.org:80/geotiff/spec/geotiff6.html#6.3.3.3

What makes things seem to work (except the sweep parameter I mean) is that GDAL
embeds a form of WKT in the PCSCitationGeoKey as a fallback when there is no proper
GeoTIFF way of representing the SRS

Given a source netCDF file with :

PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Geostationary_Satellite"],
    PARAMETER["central_meridian",-89.5],
    PARAMETER["satellite_height",35785831],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    EXTENSION["PROJ4","+proj=geos +lon_0=-89.5 +h=35785831 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs  +sweep=x"]]
Origin = (-5434865.725631997920573,-721442.352958934963681)
Pixel Size = (6012.019607998958236,6012.019607998958236)



gdal_translat'ing to TIF gives a GeoTIFF whose listgeo dump is

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTransformationTag (4,4):
         6012.01960799896  0                 0                 -5434865.725632  
         0                 6012.01960799896  0                 -721442.352958935
         0                 0                 0                 0                
         0                 0                 0                 1                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): User-Defined
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,24): "Geostationary_Satellite"
      GeographicTypeGeoKey (Short,1): GCS_WGS_84
      GeogCitationGeoKey (Ascii,13): "GCS_WGS_1984"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogInvFlatteningGeoKey (Double,1): 298.257223563    
      PCSCitationGeoKey (Ascii,383): "ESRI PE String = PROJCS["Geostationary_Satellite",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Geostationary_Satellite"],PARAMETER["central_meridian",-89.5],PARAMETER["satellite_height",35785831],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]"
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.


The WKT form used is the ESRI variant of WKT.  This is the way GDAL tries to fit a SRS in GeoTIFF when
there is no proper GeoTIFF encoding for it (for supported GeoTIFF SRS, this ESRI WKT
is not set). This is what ArcGIS is supposed to use too (not
completely sure if the massaging of OGC WKT to ESRI WKT completely matches ArcGIS though,
especially for geos which is quite new).
The issue is that ESRI WKT doesn't support the EXTENSION PROJ4 node, and hence the sweep
parameter is lost.
If there is an official way in ESRI WKT to represent the sweep parameter, then we could use it in GDAL.

So actions could be to try investigating with ESRI how they would represent that. Or as a fallback
preserve the EXTENSION node in that case, hoping that this doesn't cause too much issues to
other software (but given that the current generated ESRI WKT already lacks an important information,
interoperability is already non existent...)

A proper solution would be of course to map GEOS to GeoTIFF, but unfortunately there has not
been activity in years to make official evolutions to the specification (OGC created a charter for a GeoTIFF
SWG , http://www.opengeospatial.org/projects/groups/geotiffswg , but we haven't seen any
public activity coming out from that)


For the sake of completness, you can workaround the issue by creating a GDAL PAM .aux.xml
sidecar file

gdal_translate in.nc out.tif
gdal_translate in.nc out.png -of PNG (just to get a out.png.aux.xml with the GDAL WKT)
mv out.png.aux.xml out.tif.aux.xml
gdalinfo out.tif


Even


-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/geotiff/attachments/20170707/6b3a699b/attachment-0001.htm