WIP PROJ6 support #118

metzm · 2019-09-03T08:19:08Z

Convert lowercase epsg to uppercase EPSG
Axis order of a CRS is no longer always easting, northing, it can also be northing, easting, e.g. EPSG:4326. If source and target CRS are available, axis order is adjusted using proj_normalize_for_visualization() to easting, northing if needed.
Previously, GRASS did conversions of degress to/from radians and meters from/to map units. Now PROJ6+ might do these conversions itself, or not. GRASS is now doing the conversions only if needed.
In PROJ6+, there can be several different operations to transform coordinates from one CRS to another CRS. If more than one operation is available, information is provided about these different operations and the user has to provide one of these as pipeline option to r.proj or v.proj. That means that automated reprojection on the fly with r.import or v.import no longer works if more than one operation is available. Further on, the user must take care about axis order and remove any axisswap step from the pipeline if needed.

Example trying to convert from EPSG:2264 to EPSG:4326:

WARNING: Found 2 possible transformations
************************
Operation 1:
Description: Inverse of SPCS83 North Carolina zone (US Survey feet) + NAD83
to WGS 84 (55)

Area of use: USA - North Carolina

Accuracy within area of use: 1 m

PROJ string:
+proj=pipeline +step +proj=unitconvert +xy_in=us-ft +z_in=us-ft +xy_out=m
+z_out=m +step +inv +proj=lcc +lat_0=33.75 +lon_0=-79
+lat_1=36.1666666666667 +lat_2=34.3333333333333 +x_0=609601.219202438
+y_0=0 +ellps=GRS80 +step +proj=hgridshift +grids=nchpgn.gsb +step
+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
************************
Operation 2:
Description: Inverse of SPCS83 North Carolina zone (US Survey feet) + NAD83
to WGS 84 (1)

Area of use: North America - Canada and USA (CONUS, Alaska mainland)

Accuracy within area of use: 4 m

PROJ string:
+proj=pipeline +step +proj=unitconvert +xy_in=us-ft +z_in=us-ft +xy_out=m
+z_out=m +step +inv +proj=lcc +lat_0=33.75 +lon_0=-79
+lat_1=36.1666666666667 +lat_2=34.3333333333333 +x_0=609601.219202438
+y_0=0 +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step
+proj=axisswap +order=2,1
************************
See also output of:
projinfo -o PROJ -s "EPSG:2264" -t "EPSG:4326"
Please provide the appropriate PROJ string with the pipeline option
ERROR: Unable to initialize coordinate transformation

The problems here are 1) there is more than one operation available, 2) there is an axisswap that needs to be removed. Both problems need to be resolved by the user.

adjust for PROJ6

- Convert lowercase epsg to uppercase EPSG - Axis order of a CRS is no longer always easting, northing, it can also be northing, easting, e.g. EPSG:4326. If source and target CRS are available, axis order is adjusted using proj_normalize_for_visualization() to easting, northing if needed. - Previously, GRASS did conversions of degress to/from radians and meters from/to map units. Now PROJ6+ might do these conversions itself, or not. GRASS is now doing the conversions only if needed. - In PROJ6+, there can be several different operations to transform coordinates from one CRS to another CRS. If more than one operation is available, information is provided about these different operations and the user has to provide one of these as pipeline option to r.proj or v.proj. That means that automated reprojection on the fly with r.import or v.import no longer works if more than one operation is available. Further on, the user must take care about axis order and remove any axisswap step from the pipeline if needed.

kbevers · 2019-09-04T07:39:02Z

I am happy to review this, but please bear in mind that I am not particularly familiar with GRASS so some of my comments may not be valid. I'll start with a few comments to the description of the PR.

Convert lowercase epsg to uppercase EPSG

Is that really needed? cs2cs has no problems initializing both forms:

echo 55 12 | cs2cs EPSG:4326 EPSG:25832
691875.63       6098907.83 0.00

echo 55 12 | cs2cs epsg:4326 epsg:25832
691875.63       6098907.83 0.00

If more than one operation is available, information is provided about these different operations and the user has to provide one of these as pipeline option to r.proj or v.proj. That means that automated reprojection on the fly with r.import or v.import no longer works if more than one operation is available.

Seems unnecessary. Generally PROJ will provide the best transformation available between two CRS. I would use that by default and the give the user the option to select another one if more than one candidate is found. When using proj_create_crs_to_crs without specifying an area of use all possible transformations are initialized and the best one is used depending on the input coordinate.

Further on, the user must take care about axis order and remove any axisswap step from the pipeline if needed.

If you need the conventionel GIS axis ordering then just run your PJ through proj_normalize_for_visualization(). No need for users to faff about with pipeline strings.

rouault · 2019-09-04T07:55:01Z

lib/proj/do_proj.c

+	G_debug(1, "trying %s", info_trans->def);
+
+	/* check number of operations */
+	source_crs = proj_create(NULL, indef);


as mentionned by Kristian, proj_create_crs_to_crs() + proj_trans() should behave reasonably well if several candidate operations are available.
If keeping that check, source_crs, target_crs, operation_ctx, op_list and op should be freed

rouault · 2019-09-04T07:56:56Z

lib/proj/do_proj.c

-		insrid = G_store(info_in->srid);
+	    return -1;
+	}
+	projstr = proj_as_proj_string(NULL, info_trans->pj, PJ_PROJ_5, NULL);


coding consideration: the value of projstr only seems to be used in that part of the code, so better move the top level declaration to the beginning of this block for clarity. applies to other similar instances

rouault · 2019-09-04T07:58:01Z

lib/proj/do_proj.c

+    /* PROJ6+: EPSG must uppercase EPSG */
+    if (info_in->srid) {
+	if (strncmp(info_in->srid, "epsg", 4) == 0)
+	    insrid = G_store_upper(info_in->srid);


I'm not familiar with GRASS memory management, but assuming G_store() and G_store_upper() allocates dynamic memory, insrid doesn't seem to be freed. Same for outsrid

rouault · 2019-09-04T08:03:37Z

lib/proj/do_proj.c

+	/* what about axis order?
+	 * is it always enu? 
+	 * probably yes, as long as there is no +proj=axisswap step */
+	G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s",


in the case where info_in->def would contain a +towgs84/+nadgrids clause, I believe this pipeline would apply it, whereas you probably only want the reverse projection, and no datum shift.
The easiest would probably to mess up with the PROJ string. Otherwise with the PROJ API, you could instanciate a PJ object from the string, check if it is a BoundCRS with proj_get_source_crs(), and in that case, take the source CRS with proj_get_source_crs(), and do the inverse transform on it

rouault · 2019-09-04T08:05:52Z

Convert lowercase epsg to uppercase EPSG

Is that really needed?

At least I found recently there can be performance issues if not using EPSG uppercase, since when not finding the 'EPSG' match, PROJ will query the database for all available authority names, and this turns out to be quite slow, before doing case-insensitive checks.

support PROJ6, use area to select appropriate operation

metzm · 2019-09-13T16:37:59Z

Thanks a lot for taking time to review this PR! I followed your suggestions as much as possible. Some issues remain, however

If more than one operation is available, information is provided about these different operations and the user has to provide one of these as pipeline option to r.proj or v.proj. That means that automated reprojection on the fly with r.import or v.import no longer works if more than one operation is available.

Seems unnecessary. Generally PROJ will provide the best transformation available between two CRS. I would use that by default and the give the user the option to select another one if more than one candidate is found. When using proj_create_crs_to_crs without specifying an area of use all possible transformations are initialized and the best one is used depending on the input coordinate.

Are you really sure? See
https://github.com/OSGeo/PROJ/blob/master/src/4D_api.cpp#L1192
if there is more than one coordinate operation, "The returned P is rather dummy"

Providing a PJ_AREA helps in my local tests, a single operation is selected from the available operations.

If you need the conventionel GIS axis ordering then just run your PJ through proj_normalize_for_visualization(). No need for users to faff about with pipeline strings.

proj_normalize_for_visualization() apparently works only for a PJ created with proj_create_crs_to_crs(). If I select an operation obtained with proj_create_operations() and use the selected pipeline with proj_create(), proj_normalize_for_visualization() does not work because source and/or target crs are unknown to proj_normalize_for_visualization(). But in GRASS GIS we need to enforce an axis order of easting, northing. I think it is dangerous to expect users to manually modify the proj pipeline.

rouault · 2019-09-13T16:48:58Z

Are you really sure? See
https://github.com/OSGeo/PROJ/blob/master/src/4D_api.cpp#L1192
if there is more than one coordinate operation, "The returned P is rather dummy"

Yes, that's what the cs2cs utility uses for example (https://github.com/OSGeo/PROJ/blob/master/src/apps/cs2cs.cpp#L598). "dummy" in the sense that has it is a container for multiple internal transformations a number of query functions on PJ* object will return NULL/error, but the proj_trans() family of functions will work perfectly fine on it.

proj_normalize_for_visualization() apparently works only for a PJ created with proj_create_crs_to_crs(). If I select an operation obtained with proj_create_operations() and use the selected pipeline with proj_create(), proj_normalize_for_visualization() does not work

Indeed, but you should be able to apply proj_create_operations() on the PJ* returned by proj_list_get() when iterating on the results of proj_create_operations(). (Once you've extracted a PROJ pipeline string from a PJ* transformation and re-build a transformation with proj_create(), indeed there's no way to figure out which CRS was source and target)

use proj_normalize_for_visualization() for operations

metzm · 2019-09-14T11:22:07Z

[...] you should be able to apply proj_create_operations() on the PJ* returned by proj_list_get() when iterating on the results of proj_create_operations().

Thanks for the tip, applying proj_normalize_for_visualization() to each operation returned by proj_list_get() works for me.

support for PROJ6 continued

fix typos fix debug level

PROJ6 support - Convert lowercase epsg to uppercase EPSG - Axis order of a CRS is no longer always easting, northing, it can also be northing, easting, e.g. EPSG:4326. If source and target CRS are available, axis order is adjusted using proj_normalize_for_visualization() to easting, northing if needed. - Previously, GRASS did conversions of degress to/from radians and meters from/to map units. Now PROJ6+ might do these conversions itself, or not. GRASS is now doing the conversions only if needed. - In PROJ6+, there can be several different operations to transform coordinates from one CRS to another CRS. If more than one operation is available, information is provided about these different operations. The current region is used by PROJ to select an appropriate operation.

neteler · 2019-10-29T10:23:41Z

Backported to relbranch78 in 66b4ea2

PROJ6 support - Convert lowercase epsg to uppercase EPSG - Axis order of a CRS is no longer always easting, northing, it can also be northing, easting, e.g. EPSG:4326. If source and target CRS are available, axis order is adjusted using proj_normalize_for_visualization() to easting, northing if needed. - Previously, GRASS did conversions of degress to/from radians and meters from/to map units. Now PROJ6+ might do these conversions itself, or not. GRASS is now doing the conversions only if needed. - In PROJ6+, there can be several different operations to transform coordinates from one CRS to another CRS. If more than one operation is available, information is provided about these different operations. The current region is used by PROJ to select an appropriate operation.

metzm added 3 commits May 28, 2019 21:32

proj6 support 1

fd9fa3c

libproj

eee4662

adjust for PROJ6

metzm added enhancement New feature or request backport_needed labels Sep 3, 2019

metzm requested review from neteler, wenzeslaus, landam and mlennert September 3, 2019 08:19

metzm self-assigned this Sep 3, 2019

neteler requested review from kbevers and rouault September 4, 2019 04:07

rouault reviewed Sep 4, 2019

View reviewed changes

libproj

4c96cd1

support PROJ6, use area to select appropriate operation

rsbivand mentioned this pull request Sep 11, 2019

PROJ 6 CRS handling changes impact workflows using +datum= (and other fun) r-spatial/discuss#28

Open

lib/proj

5e853cb

use proj_normalize_for_visualization() for operations

metzm added 2 commits September 18, 2019 09:54

lib/proj

3300428

support for PROJ6 continued

lib/proj

c05298c

fix typos fix debug level

metzm mentioned this pull request Sep 25, 2019

osgeo4w: build fixes and use proj 7 and gdal 3 #136

Closed

metzm merged commit b32d8be into OSGeo:master Oct 29, 2019

metzm deleted the proj6 branch October 29, 2019 07:51

neteler removed the backport_needed label Aug 22, 2020

neteler added this to the 7.8.1 milestone Dec 9, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

WIP PROJ6 support #118

WIP PROJ6 support #118

metzm commented Sep 3, 2019

kbevers commented Sep 4, 2019

rouault Sep 4, 2019

rouault Sep 4, 2019

rouault Sep 4, 2019

rouault Sep 4, 2019

rouault commented Sep 4, 2019

metzm commented Sep 13, 2019

rouault commented Sep 13, 2019

metzm commented Sep 14, 2019

neteler commented Oct 29, 2019

WIP PROJ6 support #118

WIP PROJ6 support #118

Conversation

metzm commented Sep 3, 2019

kbevers commented Sep 4, 2019

rouault Sep 4, 2019

Choose a reason for hiding this comment

rouault Sep 4, 2019

Choose a reason for hiding this comment

rouault Sep 4, 2019

Choose a reason for hiding this comment

rouault Sep 4, 2019

Choose a reason for hiding this comment

rouault commented Sep 4, 2019

metzm commented Sep 13, 2019

rouault commented Sep 13, 2019

metzm commented Sep 14, 2019

neteler commented Oct 29, 2019