geos projection sweep parameter

Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

geos projection sweep parameter

David Hoese
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.

Dave
_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: geos projection sweep parameter

Even Rouault-2

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


_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: geos projection sweep parameter

David Hoese
Even,

So to summarize, a typical geotiff reading client should look for the
standard geotiff projection information, if not found look for "ESRI PE
String". Right? You made it sound like ArcGIS uses its own form of WKT?
Does it also use the "ESRI PE String" key?

As you stated the two "possible" solutions are to get the "geos"
projection (and others) added to the geotiff specification or to get
"sweep" added to the ESRI WKT specification. Both of which would require
client libraries to be updated to properly parse and use this
information which could take years to reach a majority of the GIS
community; maybe less given the right documentation and google search
results.

Perhaps I will start with contacting some of the working group and see
what path forward they might recommend. Thank you for your help.

Dave

On 7/6/17 7:18 PM, Even Rouault wrote:

> 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
>
_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: geos projection sweep parameter

Even Rouault-2

On lundi 10 juillet 2017 10:10:48 CEST David Hoese wrote:

> Even,

>

> So to summarize, a typical geotiff reading client should look for the

> standard geotiff projection information, if not found look for "ESRI PE

> String".

 

Yes. But ESRI PE string is vendor specific, so a strict GeoTIFF reader will generally not understand it.

 

> You made it sound like ArcGIS uses its own form of WKT?

 

WKT v1 was not strongly standardized, so various vendors have come with their own variations,

to decide how to name projection methods, parameters, etc...

Some preview :

https://web.archive.org/web/20130728081442/http://home.gdal.org:80/projects/opengis/wktproblems.html

 

> Does it also use the "ESRI PE String" key?

 

I guess so ! I don't recall the details, but I believe that support for GeoTIFF ESRI PE string in GDAL

probably originates from merging this part of the ESRI fork of GDAL to GDAL mainline.

 

>

> As you stated the two "possible" solutions are to get the "geos"

> projection (and others) added to the geotiff specification or to get

> "sweep" added to the ESRI WKT specification. Both of which would require

> client libraries to be updated to properly parse and use this

> information which could take years to reach a majority of the GIS

> community; maybe less given the right documentation and google search

> results.

 

For GDAL-based software, potentially no change on the reading side would be required.

If the ESRI PE string contains the EXTENSION PROJ4 node (which is rather hacky), GDAL will use it.

I've just tested it by manually crafting such a GeoTIFF. Works at least with GDAL 1.11 or later.

 

>

> Perhaps I will start with contacting some of the working group and see

> what path forward they might recommend.

 

I actually forwarded your email to one of the members of the OGC GeoTIFF SWG

 

Even

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff