[Proj] Integration of RFC2 work with existing PROJ code: axis order and unit issues

[Even Rouault]

Hi,

(resending, as apparently my yesterday message wasn't delivered)

Now that the RFC2 implementation has been merged
( https://github.com/OSGeo/proj.4/pull/1175 ), come the next steps. One of
them will be interation in GDAL, and then libgeotiff, so there will be
fixes
and complementary work that will be done as a result of those integrations.

For PROJ itself, one of the main thing I'd like to address is the fact that
we
have now two CRS databases:
- the old one in the data/epsg and data/IGNF files
- the new one in the SQLite3 proj.db database

The new one contains the content of the old one (with some differences,
like
the EPSG database in proj.db is fresher than in the old one). So the plan
is
to kill the old 'epsg' and 'IGNF' file that are used by the +init=file:code
syntax and go instead into the database to fetch their definitions. If no
match is found, so for other authorities than EPSG, IGNF or ESRI currently,
then the existing text file based mechanism would still be used.

All good ? Well, not quite. One of the main issue I foresee before doing
any
implementation is axis order. In the old epsg file, geographic CRS use the
longitude,latitude order, whereas proj.db and the new code I've developped
respect the axis order declared by the authority, so in the case of EPSG,
latitude,longitude for geographic CRS. This also applies for a few
projected
CRS: 'epsg' file uses generally easting,northing order (except for a few
Transverse Mercator South Orientated projected CRS where it uses the
expected
westing,southing order), whereas EPSG might have northing,easting order for
a
few ones, or southing,westing for a few Krovak-based CRS.

What to do ? I can see different options:

- just use the 'correct' axis order, and mark this as a breaking change.
That
is if you use +init=epsg:4326, you'd have to provide / will get latitude
first

- emulate the previous axis order in an attempt of being backward
compatible.
Not necessarily easy to do since that means we have to remember somewhere
the
old axis order and add the inversion when needed in some glue code between
old
and new code. And that would undermine one of the interest of the new work
to
be more axis aware than previously.

- have some switch in proj utilies ('proj', 'cs2cs', 'cct') to specify the
input and output order (like -input_axis=long,lat -output_axis=east,north)
so
that they can adapt user input/output to the expected axis of the CRS used.
That would still be a breaking change since existing scripts would have to
add
the new switches to preserve existing behaviour.

The same question also arises at the API level. For example,
proj_create_crs_to_crs(ctxt, "epsg:4326", "epsg:32632") currently expects
the
input to be in longitude, latitude order.

There is also the question of units. Traditionnaly, in code paths that
would
expect angular values, PROJ expects radians (frontend utilities do the
conversion into/from degrees). The new code assume sthe unit to be the one
declared by the coordinate system in the geographic CRS definition, ie
degrees
in 99% of the case (or grad for a couple odd & old systems like "NTF
(Paris)"
/ EPSG:4807)

Let's consider a practical example. "projinfo -s EPSG:4326 -t EPSG:4807"
outputs a transformation pipeline which is
"
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert
+xy_in=deg +xy_out=rad +step +inv +proj=hgridshift +grids=ntf_r93.gsb +step
+proj=longlat +ellps=clrk80ign +pm=paris +step +proj=unitconvert +xy_in=rad
+xy_out=grad +step +proj=axisswap +order=2,1"

The traduction in plain English of this is:
- take the input coordinate and reverse the first 2 compoments
- convert from degree to radian
- do a horizontal grid shift transformation
- apply a prime meridian shift
- convert from radian to grad
- take the ouput coordinate and reverse the first 2 compoments

That is it expects the input to be lat,long in degrees and output lat,long
in
grads, as in the offical definition of those CRS.

Should the PJ* object returned by proj_create_crs_to_crs(ctxt, "epsg:4326",
"epsg:4807") just use that pipeline, and consequently expect the
coordinates
passed into and received from proj_trans() to also use those orders and
units.

Thoughts ?

Even

PS: the database might also contain geographic CRS definitions that are
longitude,latitude by design. For example I've added an entry for the
OGC:CRS84, and if you do "projinfo -s OGC:CRS84 -t EPSG:4807", then the
pipeline doesn't contain the initial axisswap operations.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20181115/1d45a440/attachment.htm