[PROJ] question about EPSG:3003

[Even Rouault]

Hi Sandro,

First, if you instanciate with proj_create() / proj_create_from_database() and 
export with proj_as_proj_string(), you'll never get a towgs84 clause.
Since the CRS by itself has no transformation clause.
However if you look at the output of the projinfo utility or 
OSRImportFromEPSG(), you may, sometimes, get a towgs84 clause when exporting 
to PROJ string (e.g with EPSG:4322 / WGS72), because they use the 
proj_crs_create_bound_crs_to_WGS84() function which tries to find the most 
precise transformation from the CRS to EPSG:4326, and, that's the important 
point, which is valid for the whole area of use of the CRS.

For EPSG:3003, there are however 3 possibilities (you can ignore the last one 
in the below output), none of which cover the whole area of use of the 
projected CRS.

$ src/projinfo -s EPSG:3003 -t EPSG:4326 --spatial-test intersects --summary
Candidate operations found: 4
unknown id, Inverse of Italy zone 1 + Monte Mario to WGS 84 (4), 4 m, Italy - 
mainland
unknown id, Inverse of Italy zone 1 + Monte Mario to WGS 84 (2), 4 m, Italy - 
Sardinia onshore
unknown id, Inverse of Italy zone 1 + Monte Mario to WGS 84 (11), 10 m, Italy 
- Sicily Strait west of 13°E
unknown id, Inverse of Italy zone 1 + Ballpark geographic offset from Monte 
Mario to WGS 84, unknown accuracy, World, has ballpark transformation

Previous PROJ version hard-coded the transformation that has the largest area 
of use, but that's not necessarily ideal in all use cases (Sardinia or Sicily)

So, this is intended behaviour. I'm campaigning here and on all relevants 
channels: "With PROJ 6, PROJ strings to express CRS are very suboptimal. They 
should be restricted to express a transformation pipeline"

If you want Spatialite to have the same behaviour as before regarding 
generation of PROJ strings, you could possibly use the 
proj_create_operations() API to list the possible transformations, by setting 
proj_operation_factory_context_set_spatial_criterion(PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION) 
on the factory_ctx, pickup the one that has the largest area of use, create 
the final bound CRS with that transformation with proj_crs_create_bound_crs(), 
and export it to a PROJ string.
Admitedly, we could have a releaxed proj_crs_create_bound_crs_to_WGS84() 
implementation in PROJ itself to do that, but I'm not sure if we really need 
to replicate/encourage past suboptimal behaviour.

Even

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