[Proj] Integration of RFC2 work with existing PROJ code: axis order and unit issues
even.rouault at spatialys.com
Wed Nov 14 15:15:48 EST 2018
Now that RFC2 implementation is ready to be 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)"
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.
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.
Spatialys - Geospatial professional services
More information about the Proj