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

[Even Rouault]

Hi,

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)" 
/ 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.

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