[gdal-dev] Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

Even Rouault-2
Hi,

I call this pre-RFC as I don't have yet any code implementing the following
ideas (so there might be holes in the analysis), but I wanted to get some
feedback before investing too much time into this substantial work.

Let's start with the easiest topic:

1) Use of OGRSpatialReference in GDAL

Currently GDAL datasets accept and return a WKT 1 string to describe the SRS.
To be more independant of the actual encoding, and for example allowing a
GeoPackage raster dataset to be able to use WKT 2, I propose we switch from a
string representation to the more abstract representation offered by
OGRSpatialReference. In the current state of GDAL, OGRSpatialReference is just
a direct mapping of WKT 1, but in my in-progress gdalbarn branch at
https://github.com/rouault/gdal/tree/gdalbarn , OGRSpatialReference is now
more or less encapsulating a osgeo::proj::crs object (ie a C++ object
representing a CRS as modelled by OGC Topic 2 / ISO-19111), that can have
several representations: OGC WKT1 (GDAL WKT1), WKT2, ESRI WKT1, PROJ strings.

So the proposal would be to change:

GDALDataset:
    virtual const char* GetProjectionRef();

to
    virtual OGRSpatialReference *GetSpatialRef();

    // compatibility layer (conceptually)
    const char* GetProjectionRef() {
                return GetSpatialRef()->ExportToWkt1(); }

Similarly on the set side,

Change
    virtual CPLErr SetProjection(const char* pszWKT);
to
    virtual CPLErr SetSpatialRef(OGRSpatialReference*);

    // compatibility layer (conceptually)
    CPLErr SetProjection(const char* pszWKT) {
                return SetSpatialRef(OGRSpatialReference::importFromWkt(pszWKT)); }

Regarding ownership, I guess we would do similarly to the vector side, that is
GetSpatialRef() returns an an object that it owns and doesn't increment the
refcount. And SetSpatialRef() would potentially increment the reference count
on the passed object or clone it (like ICreateLayer() does)

To ease the transition that will require going through all the ~150 raster
drivers, I may add in GDALDataset:
  protected:
    OGRSpatialReference* GetSpatialRefFromOldGetProjectionRef()
                // conceptual implementation
                { return OGRSpatialReference::importFromWkt(OldGetProjectionRef()); }

    virtual const char* OldGetProjectionRef() { return ""; }

Drivers would rename their existing GetProjectionRef() to
OldGetProjectionRef() and implement a GetSpatialRef() that would just call
GetSpatialRefFromOldGetProjectionRef(). This is just a technical detail to be
able to complete the transition of API easier in a semi-automatic way. Over
time we may remove the hack and implement natively GetSpatialRef(). And
similarly for the write side.


Now, the harder one:

2) Axis order revamp

Ever recurring headache [1], not sure my below proposal will solve it in a
definitive way

Let's look at the current situation with a common example

EPSG:4326 returns
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"]]

whereas EPSGA:4326 returns (not the A for Axis)
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"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]]

In EPSG:4326 expansion, there is no AXIS mention. This is technically
conformant with OGC WKT1 grammar where AXIS is optional. However, GDAL uses
this omission to implicitely mean it doesn't honour the axis order that the
EPSG authority has set for this geographic CRS, which is latitude first,
longitude second, but instead use longitude, latitude. Note also that EPSGA:
4326 is not really honoured by GDAL: if you use it with gdaltransform, it will
still expect/return coordinates in longitude, latitude order.

In WKT2, AXIS is a required element, so we cannot hide ourselves behind the
omission we did in WKT1, and we can't either modify the definition to
advertize longitude, latitude order, as it would be a great source of
confusion (using an authority code but putting a contradictory element into it
is a no-no)

I propose to add a "data axis to SRS axis mapping" concept, which is a bit
similar to what is done in WCS DescribeCoverage response to explain how the
SRS axis map to the grid axis of a coverage

Extract from
https://docs.geoserver.org/stable/en/user/extensions/wcs20eo/index.html
for a coverage that uses EPSG:4326
         <gml:coverageFunction>
        <gml:GridFunction>
          <gml:sequenceRule axisOrder="+2 +1">Linear</gml:sequenceRule>
          <gml:startPoint>0 0</gml:startPoint>
        </gml:GridFunction>
      </gml:coverageFunction>

A similar mapping would be used to define how the 'x' and 'y' components in
the geotransform matrix or in a OGRGeometry map to the axis defined by the CRS
definition.

Such mapping would be given by a new method in OGRSpatialReference

    std::vector<int> GetDataAxisToSRSAxisMapping() const

To explain its semantics, imagine that it return 2,-1,3.
That would be interpreted as:
* 2: the first axis of the CRS maps to the second axis of the data
* -1: the second axis of the CRS maps to the first axis of the data, with
values negated
* 3: the third axis of the CRS maps to the third axis of the data

This would be similar to the PROJ axisswap operation:
        https://proj4.org/operations/conversions/axisswap.html

By default, GetDataAxisToSRSAxisMapping() would return the identity 1,2[,3],
that is be conformant with the axis order defined by the authority.

As all GDAL and a vast majority of OGR drivers depend on using the "GIS axis
mapping", a method
        SetDataAxisToSRSAxisMappingStrategy(
                        TRADITIONAL_GIS_ORDER or AUTHORITY_COMPLIANT )
would be added to make their job of specifying the axis mapping easier;

TRADITIONAL_GIS_ORDER would mean:
* for geographic 2D CRS,
    - for Latitude NORTH, Longitude EAST (such as EPSG:4326),
GetDataAxisToSRSAxisMapping() would return {2,1}, meaning that the data order
is longitude, latitude
    - for Longitude EAST, Latitude NORTH (such as OGC:CRS84), would return
{1,2}

* for projected CRS,
    - for EAST, NORTH (ie most projected CRS), would return {1,2}
    - for NORTH, EAST, would return {2,1}
    - for North Pole CRS, with East/SOUTH, North/SOUTH, such as EPSG:5041
("WGS 84 / UPS North (E,N)"), would return {1,2}
    - for North Pole CRS, with northing/SOUTH, easting/SOUTH, such as EPSG:
32661 ("WGS 84 / UPS North (N,E)"), would return {2,1}
    - similarly for South Pole CRS
    - Krovak East North, with EAST, NORTH, such as EPSG:5221 "S-JTSK (Ferro) /
Krovak", would return {1,2}

        Cases where I'm less sure about what the 'GIS order' is and suspect we
might not do the right thing currently:
    - Transverse Mercator South Oriented, with WEST, SOUTH, would return {1,2}  
(I assume that this system is not used for rasters)
    - traditional Krovak, with SOUTH, WEST, such as EPSG:2065 "S-JTSK (Ferro)
/ Krovak", would return {1,2}   (I assume that this system is not used for
rasters)
    - other cases, return {1,2}

OGRCreateCoordinateTransformation(OGRSpatialReference* source,
OGRSpatialReference* target) will honour the data axis to srs axis mapping
That means that gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32631 will expect
latitude, longitude order (so backward incompatible)
(PROJ6 will also be axis order compliant when using EPSG codes)
We may add a -s_axis gis switch (and/or -s_axis long,lat for more explicit) /
-t_axis syntax to have some easier transition.

Raster datasets will all set
SetDataAxisToSRSAxisMappingStrategy(TRADITIONAL_GIS_ORDER) on the
OGRSpatialReference* they return, and will also assume it in SetSpatialRef()
(will probably be just assumed and unchecked for now)

Vector layers will mostly all call
SetDataAxisToSRSAxisMappingStrategy(TRADITIONAL_GIS_ORDER) on the
OGRSpatialReference* returned by GetSpatialRef(). We may decide that some
drivers like GML no longer do coordinate axis swapping for GML3 and then use
the AUTHORITY_COMPLIANT strategy.
ICreateLayer() when receiving a OGRSpatialReference* may decide to change the
axis mapping strategy. That is: if it receives a OGRSpatialReference with
AUTHORITY_COMPLIANT order, it may decide to switch to TRADITIONAL_GIS_ORDER
and GetSpatialRef()::GetDataAxisToSRSAxisMapping() will reflect that. Tools
like ogr2ogr will need to do the geometry axis swapping in that case

I'm midly satisfied with adding the axis mapping in the OGRSpatialReference
class itself, but I feel it would be overkill to have an extra object to
capture both (and can't come with a name for it). The coupling of both will
require careful explanation of what some methods do. For example, I think
OGRSpatialReference::IsSame() would ignore the data axis mapping setting when
comparing 2 OGRSpatialReference object.

I would also lean on making WKT 1 export now always return the AXIS element,
and EPSG:xxxx would then behave like EPSGA:xxxx

So a summary view of this approach is that in the formal SRS definition, we no
longer do derogations regarding axis order, but we add an additional interface
to describe how we actually make our match match with the SRS definition.

In case you reached that point, thoughts... ?

Even

[1] Fresh ticket about such confusion:
https://github.com/OSGeo/gdal/issues/1155

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

Andrew C Aitchison-2

Do we need to consider (real and false) precision when converting between
OGRSpatialReference and string ?

On Tue, 18 Dec 2018, Even Rouault wrote:

> Hi,
>
> I call this pre-RFC as I don't have yet any code implementing the following
> ideas (so there might be holes in the analysis), but I wanted to get some
> feedback before investing too much time into this substantial work.
>
> Let's start with the easiest topic:
>
> 1) Use of OGRSpatialReference in GDAL
>
> Currently GDAL datasets accept and return a WKT 1 string to describe the SRS.
> To be more independant of the actual encoding, and for example allowing a
> GeoPackage raster dataset to be able to use WKT 2, I propose we switch from a
> string representation to the more abstract representation offered by
> OGRSpatialReference. In the current state of GDAL, OGRSpatialReference is just
> a direct mapping of WKT 1, but in my in-progress gdalbarn branch at
> https://github.com/rouault/gdal/tree/gdalbarn , OGRSpatialReference is now
> more or less encapsulating a osgeo::proj::crs object (ie a C++ object
> representing a CRS as modelled by OGC Topic 2 / ISO-19111), that can have
> several representations: OGC WKT1 (GDAL WKT1), WKT2, ESRI WKT1, PROJ strings.
>
> So the proposal would be to change:
>
> GDALDataset:
>    virtual const char* GetProjectionRef();
>
> to
>    virtual OGRSpatialReference *GetSpatialRef();
>
>    // compatibility layer (conceptually)
>    const char* GetProjectionRef() {
> return GetSpatialRef()->ExportToWkt1(); }
>
> Similarly on the set side,
>
> Change
>    virtual CPLErr SetProjection(const char* pszWKT);
> to
>    virtual CPLErr SetSpatialRef(OGRSpatialReference*);
>
>    // compatibility layer (conceptually)
>    CPLErr SetProjection(const char* pszWKT) {
> return SetSpatialRef(OGRSpatialReference::importFromWkt(pszWKT)); }
>
> Regarding ownership, I guess we would do similarly to the vector side, that is
> GetSpatialRef() returns an an object that it owns and doesn't increment the
> refcount. And SetSpatialRef() would potentially increment the reference count
> on the passed object or clone it (like ICreateLayer() does)
>
> To ease the transition that will require going through all the ~150 raster
> drivers, I may add in GDALDataset:
>  protected:
>    OGRSpatialReference* GetSpatialRefFromOldGetProjectionRef()
> // conceptual implementation
> { return OGRSpatialReference::importFromWkt(OldGetProjectionRef()); }
>
>    virtual const char* OldGetProjectionRef() { return ""; }
>
> Drivers would rename their existing GetProjectionRef() to
> OldGetProjectionRef() and implement a GetSpatialRef() that would just call
> GetSpatialRefFromOldGetProjectionRef(). This is just a technical detail to be
> able to complete the transition of API easier in a semi-automatic way. Over
> time we may remove the hack and implement natively GetSpatialRef(). And
> similarly for the write side.

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

Even Rouault-2
On mercredi 19 décembre 2018 06:30:30 CET Andrew C Aitchison wrote:
> Do we need to consider (real and false) precision when converting between
> OGRSpatialReference and string ?

Can you explain more what you mean ? I don't understand.

--
Spatialys - Geospatial professional services
http://www.spatialys.com

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

Andrew C Aitchison-2
On Wed, 19 Dec 2018, Even Rouault wrote:

> On mercredi 19 décembre 2018 06:30:30 CET Andrew C Aitchison wrote:
>> Do we need to consider (real and false) precision when converting between
>> OGRSpatialReference and string ?

My concern is that since since a WKT1 string encodes floating point
values as strings and an OGRSpatialReference stores them as doubles,
a round trip (eg saving to file and reading back, or vice versa)
will sometimes not preserve the last digit(s).

My reference to "false precision" was something like 1.1 becoming
eg. 1.1000000000001

Looking in more detail at your proposal, I guess this wont bring up
anything new, since we are using the existing routines
OGRSpatialReference::importFromWkt() and GetSpatialRef()->ExportToWkt1()
where this issue will have already been considered.

--
Andrew C. Aitchison Cambridge, UK
  [hidden email]
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

Even Rouault-2
On mercredi 19 décembre 2018 11:29:09 CET Andrew C Aitchison wrote:

> On Wed, 19 Dec 2018, Even Rouault wrote:
> > On mercredi 19 décembre 2018 06:30:30 CET Andrew C Aitchison wrote:
> >> Do we need to consider (real and false) precision when converting between
> >> OGRSpatialReference and string ?
>
> My concern is that since since a WKT1 string encodes floating point
> values as strings and an OGRSpatialReference stores them as doubles,
> a round trip (eg saving to file and reading back, or vice versa)
> will sometimes not preserve the last digit(s).
>
> My reference to "false precision" was something like 1.1 becoming
> eg. 1.1000000000001
>
> Looking in more detail at your proposal, I guess this wont bring up
> anything new, since we are using the existing routines
> OGRSpatialReference::importFromWkt() and GetSpatialRef()->ExportToWkt1()
> where this issue will have already been considered.

Indeed, nothing new on this. Except that the import/export from/to WKT will be
done by PROJ 6 now. It uses %.15g formatting

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

jratike80
In reply to this post by Even Rouault-2
Even Rouault-2 wrote
> ICreateLayer() when receiving a OGRSpatialReference* may decide to change
> the
> axis mapping strategy. That is: if it receives a OGRSpatialReference with
> AUTHORITY_COMPLIANT order, it may decide to switch to
> TRADITIONAL_GIS_ORDER
> and GetSpatialRef()::GetDataAxisToSRSAxisMapping() will reflect that.
> Tools
> like ogr2ogr will need to do the geometry axis swapping in that case

Does this guarantee that it is impossible or at least pretty difficult to
create for example crazy shapefiles with flipped axis and rotated geometries
even if user tries to force that with AUTHORITY_COMPLIANT mode switch?

-Jukka Rahkonen-




--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

Even Rouault-2
On mercredi 19 décembre 2018 05:46:26 CET jratike80 wrote:

> Even Rouault-2 wrote
>
> > ICreateLayer() when receiving a OGRSpatialReference* may decide to change
> > the
> > axis mapping strategy. That is: if it receives a OGRSpatialReference with
> > AUTHORITY_COMPLIANT order, it may decide to switch to
> > TRADITIONAL_GIS_ORDER
> > and GetSpatialRef()::GetDataAxisToSRSAxisMapping() will reflect that.
> > Tools
> > like ogr2ogr will need to do the geometry axis swapping in that case
>
> Does this guarantee that it is impossible or at least pretty difficult to
> create for example crazy shapefiles with flipped axis and rotated geometries
> even if user tries to force that with AUTHORITY_COMPLIANT mode switch?

Hi Jukka,

<jokingly>
"""THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
OR IMPLIED"""
</jokingly>

More seriously, the goal is indeed to try to make things work as currently for
ogr2ogr, gdalwarp workflows by making drivers to report properly how they map
the data axis to the CRS axis.
For now, the only backward incompatible change I can think of is the use of
the low-level OGRCoordinateTransformation() interface (and gdaltransform),
that will be axis compliant w.r.t the authoritative CRS definition, but that's
a wanted feature I think.
With the caveat that this is only words for now, and the work remains to be
done and will certainly bring a few surprises.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev