[GRASS-dev] PROJ6 support in GRASS

GRASS GRASS developer
1

Hi all,

there is a new PR for PROJ6 support in GRASS
https://github.com/OSGeo/grass/pull/118

There are two important changes when using PROJ6:

First, reprojection with v.proj and r.proj is no longer always possible without the user making informed decisions. The reason is that there can be several different operations available to reproject coordinates from one CRS to another CRS. These different operations are listed and the user has to provide the appropriate operation with the pipeline option, taking care of any axisswap.

Second, axis order is no longer always easting, northing, e.g. for EPSG:4326 it is northing, easting, and an axisswap might need to be removed from operations provided by PROJ.

There are many more changes (see details in the PR), but these are the two most important ones.

Feedback welcome!

2

Hi,

On Tue, Sep 3, 2019 at 4:22 AM Markus Metz <markus.metz.giswork@gmail.com> wrote:

Hi all,

there is a new PR for PROJ6 support in GRASS
https://github.com/OSGeo/grass/pull/118

There are two important changes when using PROJ6:

First, reprojection with v.proj and r.proj is no longer always possible without the user making informed decisions. The reason is that there can be several different operations available to reproject coordinates from one CRS to another CRS. These different operations are listed and the user has to provide the appropriate operation with the pipeline option, taking care of any axisswap.

first, thanks for all the work. Second, I don’t see how most users are supposed to know what to pick. Is there perhaps a way to pick a good default? I just can’t imagine not having r.import/v.import…

Anna

Second, axis order is no longer always easting, northing, e.g. for EPSG:4326 it is northing, easting, and an axisswap might need to be removed from operations provided by PROJ.

There are many more changes (see details in the PR), but these are the two most important ones.

Feedback welcome!

grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

3

Markus,
Thank you for working on this!

Perhaps a flag to default to the transformation with the least error ( via something like the projinfo -s x -t y --summary query and with details captured about the option selected in the metadata)? https://media.ccc.de/v/bucharest-198-revamp-of-coordinate-reference-system-management-in-the-osgeo-c-c-stack-with-proj-and-gdal#t=666 around 10:33 minutes in.

Looking at the reprojection options in QGIS 3.8.2 the menu for picking the reprojection option lists the accuracy, the number of contributing stations, and a preferred alternative. If this data can be gleaned from a projinfo query, a preferred alternative can be established. I would not make that the default.

A sample workflow would be that someone tries r.proj with no flags. If there are multiple alternatives, it does not proceed and gives a message along the lines of " Multiple projection options detected, please rerun with --info flag to see alternatives.

User reruns with --info flag ( which supersedes all other flags ) and sees 5 alternatives with accuracy and number of stations listed with an asterisk next to the “preferred” alternative

User runs r.proj with a --prefer flag that takes the option with the least error and most stations. Alternatively, --file flag that points to a text file with a proj6 pipeline to enforce exactly what is desired.

The preferred option may not always be the same over time, but if the pipeline is captured in the metadata you will have a record of how it was done.

Doug

On Tue, Sep 3, 2019 at 10:31 PM Anna Petrášová <kratochanna@gmail.com> wrote:

Hi,

On Tue, Sep 3, 2019 at 4:22 AM Markus Metz <markus.metz.giswork@gmail.com> wrote:

Hi all,

there is a new PR for PROJ6 support in GRASS
https://github.com/OSGeo/grass/pull/118

There are two important changes when using PROJ6:

First, reprojection with v.proj and r.proj is no longer always possible without the user making informed decisions. The reason is that there can be several different operations available to reproject coordinates from one CRS to another CRS. These different operations are listed and the user has to provide the appropriate operation with the pipeline option, taking care of any axisswap.

first, thanks for all the work. Second, I don’t see how most users are supposed to know what to pick. Is there perhaps a way to pick a good default? I just can’t imagine not having r.import/v.import…

Anna

Second, axis order is no longer always easting, northing, e.g. for EPSG:4326 it is northing, easting, and an axisswap might need to be removed from operations provided by PROJ.

There are many more changes (see details in the PR), but these are the two most important ones.

Feedback welcome!

grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

–

Doug Newcomb - Cartographer
USFWS
551F Pylon Dr
Raleigh, NC
919-856-4520 ext. 14 doug_newcomb@fws.gov

NOTE: This email correspondence and any attachments to and from this sender is subject to the Freedom of Information Act (FOIA) and may be disclosed to third parties.

4

On Wed, Sep 4, 2019 at 4:31 AM Anna Petrášová <kratochanna@gmail.com> wrote:

On Tue, Sep 3, 2019 at 4:22 AM Markus Metz <markus.metz.giswork@gmail.com> wrote:

Hi all,

there is a new PR for PROJ6 support in GRASS
https://github.com/OSGeo/grass/pull/118

There are two important changes when using PROJ6:

First, reprojection with v.proj and r.proj is no longer always possible without the user making informed decisions. The reason is that there can be several different operations available to reproject coordinates from one CRS to another CRS. These different operations are listed and the user has to provide the appropriate operation with the pipeline option, taking care of any axisswap.

first, thanks for all the work. Second, I don't see how most users are supposed to know what to pick. Is there perhaps a way to pick a good default? I just can't imagine not having r.import/v.import...

I see the same problem: users won't know what to select if defaults
are gone with PROJ 6.

Anna

Second, axis order is no longer always easting, northing, e.g. for EPSG:4326 it is northing, easting, and an axisswap might need to be removed from operations provided by PROJ.

There are many more changes (see details in the PR), but these are the two most important ones.

Feedback welcome!

Thanks for your hard work.

I talked to Robert Bivand here at the GeoSTAT 2019 in Münster who gave
a talk about the topic:

"Not just R-spatial: sustaining open source geospatial software stacks"
https://github.com/rsbivand/geostat19_talk

You can quick-view the talk in rendered HTML like this (maybe there
are better ways):
http://htmlpreview.github.io/?https://raw.githubusercontent.com/rsbivand/geostat19_talk/master/docs/geostat_bivand_19.html

It contains a section about PROJ (6).

markusN

5

On Wed, Sep 4, 2019 at 6:01 PM Markus Neteler <neteler@osgeo.org> wrote:

On Wed, Sep 4, 2019 at 4:31 AM Anna Petrášová <kratochanna@gmail.com> wrote:

On Tue, Sep 3, 2019 at 4:22 AM Markus Metz <markus.metz.giswork@gmail.com> wrote:

Hi all,

there is a new PR for PROJ6 support in GRASS
https://github.com/OSGeo/grass/pull/118

There are two important changes when using PROJ6:

First, reprojection with v.proj and r.proj is no longer always possible without the user making informed decisions. The reason is that there can be several different operations available to reproject coordinates from one CRS to another CRS. These different operations are listed and the user has to provide the appropriate operation with the pipeline option, taking care of any axisswap.

first, thanks for all the work. Second, I don’t see how most users are supposed to know what to pick. Is there perhaps a way to pick a good default? I just can’t imagine not having r.import/v.import…

I see the same problem: users won’t know what to select if defaults
are gone with PROJ 6.

We can provide information that enables users to make an informed decision. The options listed by PROJ6 are ordered by guessed best choice, i.e. the first one is, given the provided information, the best choice. But users need to review the list of options.

Reprojection from one CRS to another CRS was never easy. For both raster and vector data, some preparation was always needed to decide on appropriate settings. PROJ6 provides methods to improve the accuracy of reprojected coordinates, depending on the actual use case. The user decides (must decide).

Anna

Second, axis order is no longer always easting, northing, e.g. for EPSG:4326 it is northing, easting, and an axisswap might need to be removed from operations provided by PROJ.

There are many more changes (see details in the PR), but these are the two most important ones.

Feedback welcome!

Thanks for your hard work.

I talked to Robert Bivand here at the GeoSTAT 2019 in Münster who gave
a talk about the topic:

“Not just R-spatial: sustaining open source geospatial software stacks”
https://github.com/rsbivand/geostat19_talk

You can quick-view the talk in rendered HTML like this (maybe there
are better ways):
http://htmlpreview.github.io/?https://raw.githubusercontent.com/rsbivand/geostat19_talk/master/docs/geostat_bivand_19.html

It contains a section about PROJ (6).

Not so random citation:
“…we need to manipulate the CRS read in with the file to insert our choice of how to make the transformation…”

I will try to restrict the number of options based on the current region in a modified PR.

Markus M

6

Several different methods might be available to reproject from one CRS to another CRS. PROJ6 can select the most appropriate method if a bounding box is provided. This means that the selected method depends on the bounding box and that results of a reprojection can differ depending on the provided bounding box.

In GRASS, this bounding box can be obtained from the current region. That means if the current region changes, the method and thus the results of coordinate reprojection might change. This effect can not be underestimated.

With the current pull request #118 the current region is used to help PROJ6 select the best method, and [r|v].import should work as before. Results can be different compared to PROJ5 or earlier.

IMHO, there is no way around that users become more familiar with details of coordinate reprojection and need to read the output of [r|v].proj carefully regarding different methods known to PROJ6.

PROJ6 will not use the best method if any required datum transformation grid is not available. Users will need to obtain the corresponding grid themselves, which raises another question: where to put this grid? On Linux, this would be e.g. /usr/share/proj, but if the user can not write to /usr/share/proj, the grids need to be saved somewhere in $HOME, and PROJ6 must look at that place. This problem, i.e. the place where PROJ6 should look for in $HOME is not yet solved AFAICT (proj-6.2.0).

Markus M

On Wed, Sep 4, 2019 at 9:46 PM Markus Metz <markus.metz.giswork@gmail.com> wrote:

On Wed, Sep 4, 2019 at 6:01 PM Markus Neteler <neteler@osgeo.org> wrote:

On Wed, Sep 4, 2019 at 4:31 AM Anna Petrášová <kratochanna@gmail.com> wrote:

On Tue, Sep 3, 2019 at 4:22 AM Markus Metz <markus.metz.giswork@gmail.com> wrote:

Hi all,

there is a new PR for PROJ6 support in GRASS
https://github.com/OSGeo/grass/pull/118

There are two important changes when using PROJ6:

First, reprojection with v.proj and r.proj is no longer always possible without the user making informed decisions. The reason is that there can be several different operations available to reproject coordinates from one CRS to another CRS. These different operations are listed and the user has to provide the appropriate operation with the pipeline option, taking care of any axisswap.

first, thanks for all the work. Second, I don’t see how most users are supposed to know what to pick. Is there perhaps a way to pick a good default? I just can’t imagine not having r.import/v.import…

I see the same problem: users won’t know what to select if defaults
are gone with PROJ 6.

We can provide information that enables users to make an informed decision. The options listed by PROJ6 are ordered by guessed best choice, i.e. the first one is, given the provided information, the best choice. But users need to review the list of options.

Reprojection from one CRS to another CRS was never easy. For both raster and vector data, some preparation was always needed to decide on appropriate settings. PROJ6 provides methods to improve the accuracy of reprojected coordinates, depending on the actual use case. The user decides (must decide).

Anna

Second, axis order is no longer always easting, northing, e.g. for EPSG:4326 it is northing, easting, and an axisswap might need to be removed from operations provided by PROJ.

There are many more changes (see details in the PR), but these are the two most important ones.

Feedback welcome!

Thanks for your hard work.

I talked to Robert Bivand here at the GeoSTAT 2019 in Münster who gave
a talk about the topic:

“Not just R-spatial: sustaining open source geospatial software stacks”
https://github.com/rsbivand/geostat19_talk

You can quick-view the talk in rendered HTML like this (maybe there
are better ways):
http://htmlpreview.github.io/?https://raw.githubusercontent.com/rsbivand/geostat19_talk/master/docs/geostat_bivand_19.html

It contains a section about PROJ (6).

Not so random citation:
“…we need to manipulate the CRS read in with the file to insert our choice of how to make the transformation…”

I will try to restrict the number of options based on the current region in a modified PR.

Markus M

7

On mercredi 25 septembre 2019 22:47:21 CEST Markus Metz wrote:

Several different methods might be available to reproject from one CRS to
another CRS. PROJ6 can select the most appropriate method if a bounding box
is provided. This means that the selected method depends on the bounding
box and that results of a reprojection can differ depending on the provided
bounding box.

In GRASS, this bounding box can be obtained from the current region. That
means if the current region changes, the method and thus the results of
coordinate reprojection might change. This effect can not be underestimated.

With the current pull request #118 the current region is used to help PROJ6
select the best method, and [r|v].import should work as before. Results can
be different compared to PROJ5 or earlier.

IMHO, there is no way around that users become more familiar with details
of coordinate reprojection and need to read the output of [r|v].proj
carefully regarding different methods known to PROJ6.

PROJ6 will not use the best method if any required datum transformation
grid is not available. Users will need to obtain the corresponding grid
themselves, which raises another question: where to put this grid? On
Linux, this would be e.g. /usr/share/proj, but if the user can not write to
/usr/share/proj, the grids need to be saved somewhere in $HOME, and PROJ6
must look at that place. This problem, i.e. the place where PROJ6 should
look for in $HOME is not yet solved AFAICT (proj-6.2.0).

Applications may decide for an appropriate user directory and set it with
proj_context_set_search_paths() (in that case this overrides PROJ_LIB or
hardcoded directories, so you have to add in the search paths if this is the
desired behaviour)

QGIS has code to do exactly that:
https://github.com/qgis/QGIS/blob/65359bc7eafbfe967c669d1428eeedb87c5cd2a1/
src/core/qgsapplication.cpp#L312
and
https://github.com/qgis/QGIS/blob/d45c3dd4f1e8fe3642d42d84aa978b28dba913aa/
src/core/qgsprojutils.cpp#L258

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

8

On Wed, Sep 25, 2019 at 10:55 PM Even Rouault <even.rouault@spatialys.com> wrote:

On mercredi 25 septembre 2019 22:47:21 CEST Markus Metz wrote:

PROJ6 will not use the best method if any required datum transformation
grid is not available. Users will need to obtain the corresponding grid
themselves, which raises another question: where to put this grid? On
Linux, this would be e.g. /usr/share/proj, but if the user can not write to
/usr/share/proj, the grids need to be saved somewhere in $HOME, and PROJ6
must look at that place. This problem, i.e. the place where PROJ6 should
look for in $HOME is not yet solved AFAICT (proj-6.2.0).

Applications may decide for an appropriate user directory and set it with
proj_context_set_search_paths() (in that case this overrides PROJ_LIB or
hardcoded directories, so you have to add in the search paths if this is the
desired behaviour)

QGIS has code to do exactly that:
https://github.com/qgis/QGIS/blob/65359bc7eafbfe967c669d1428eeedb87c5cd2a1/
src/core/qgsapplication.cpp#L312
and
https://github.com/qgis/QGIS/blob/d45c3dd4f1e8fe3642d42d84aa978b28dba913aa/
src/core/qgsprojutils.cpp#L258

The (my) idea is that proj suggests reasonable paths to store grids, search paths should be available in PJ_INFO. Otherwise applications may decide for a directory specific to that application, not generally valid for proj used by different applications. Applications can then check the existing proj search paths and install grids to a path where the current user has write access. Thus no need to use proj_context_set_search_paths() and no need to install the same grid(s) repeatedly.

However, there are no default search paths in PJ_INFO. Can you give a hint where the default search paths accessible in the C API?

Markus M

9

The (my) idea is that proj suggests reasonable paths to store grids, search
paths should be available in PJ_INFO. Otherwise applications may decide for
a directory specific to that application, not generally valid for proj used
by different applications. Applications can then check the existing proj
search paths and install grids to a path where the current user has write
access. Thus no need to use proj_context_set_search_paths() and no need to
install the same grid(s) repeatedly.

That's not adressed indeed. Currently each application has to handle that for
itself.
I may work at a later point on something related, which is that PROJ would
automatically download the grids (actually the chunk of the grid it needs) it
needs if they are not already available.
That's the discussion in
https://lists.osgeo.org/pipermail/proj/2019-September/008858.html
We're currently half of that effort funded; we're looking for the other half

However, there are no default search paths in PJ_INFO. Can you give a hint
where the default search paths accessible in the C API?

Hum, proj_info().searchpath ? But I may not understand your question

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

10

PROJ 6 support is not yet finished in GRASS. Currently it sort of works in master and relbr78, but I expect troubles because GRASS is heavily relying on deprecated proj strings as CRS definitions. This is becoming problematic e.g. when creating locations from GDAL/OGR datasets because GRASS goes from OGR spatialreference through proj string to GRASS definition.

For GDAL 3 + PROJ 6, GRASS should switch to WKT as main CRS definition format, ideally also storing a WKT definition next to PROJ_INFO. When creating locations from GDAL/OGR datasets, GRASS should go from OGR spatialreference through WKT to GRASS definition, and probably use the new PROJ6+ methods as much as possible.

When getting the projection info of a GRASS location, GRASS should then check in that order 1) WKT, 2) EPSG or other SRID, 3) GRASS native definition (which deviated from proj definition some decades ago).

That means, g.proj, r.in.gdal, r.external, v.in.ogr, v.external need updates, and the mechanism of conversion between different formats of CRS definitions in lib/proj needs to be largely rewritten.

Markus M

11

Hi,

čt 7. 11. 2019 v 15:17 odesílatel Markus Metz
<markus.metz.giswork@gmail.com> napsal:

PROJ 6 support is not yet finished in GRASS. Currently it sort of works in master and relbr78, but I expect troubles because GRASS is heavily relying on deprecated proj strings as CRS definitions. This is becoming problematic e.g. when creating locations from GDAL/OGR datasets because GRASS goes from OGR spatialreference through proj string to GRASS definition.

OSGeo4W packages are planned to be switched to GDAL3/PROJ6 soon, are
you aware of any known bugs?

For GDAL 3 + PROJ 6, GRASS should switch to WKT as main CRS definition format, ideally also storing a WKT definition next to PROJ_INFO. When creating locations from GDAL/OGR datasets, GRASS should go from OGR spatialreference through WKT to GRASS definition, and probably use the new PROJ6+ methods as much as possible.

Yes, good topic for
https://grasswiki.osgeo.org/wiki/GRASS_GIS_Community_Sprint_Prague_2019
:slight_smile:

When getting the projection info of a GRASS location, GRASS should then check in that order 1) WKT, 2) EPSG or other SRID, 3) GRASS native definition (which deviated from proj definition some decades ago).

That means, g.proj, r.in.gdal, r.external, v.in.ogr, v.external need updates, and the mechanism of conversion between different formats of CRS definitions in lib/proj needs to be largely rewritten.

This task for next stable GRASS release (7.10?) I guess.

Ma

--
Martin Landa
http://geo.fsv.cvut.cz/gwiki/Landa
http://gismentors.cz/mentors/landa

12

On Thu, Nov 7, 2019 at 6:47 PM Martin Landa <landa.martin@gmail.com> wrote:

Hi,

čt 7. 11. 2019 v 15:17 odesílatel Markus Metz
<markus.metz.giswork@gmail.com> napsal:

PROJ 6 support is not yet finished in GRASS. Currently it sort of works in master and relbr78, but I expect troubles because GRASS is heavily relying on deprecated proj strings as CRS definitions. This is becoming problematic e.g. when creating locations from GDAL/OGR datasets because GRASS goes from OGR spatialreference through proj string to GRASS definition.

OSGeo4W packages are planned to be switched to GDAL3/PROJ6 soon, are
you aware of any known bugs?

yes, EPSG:3857 is causing trouble. This is the WGS 84 / Pseudo-Mercator projection on a sphere with radius 6378137 that needs special treatment which no longer works with GRASS + PROJ 6 + GDAL 3.

g.proj epsg=3857 -p gives me
-PROJ_INFO-------------------------------------------------
name : WGS 84 / Pseudo-Mercator
datum : wgs84
ellps : wgs84
proj : merc
lat_ts : 0
lon_0 : 0
x_0 : 0
y_0 : 0
k : 1
wktext : defined
no_defs : defined
-PROJ_EPSG-------------------------------------------------
epsg : 3857
-PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1

which is wrong because it should be a Mercator projection on a sphere, not an ellipsoid

from https://proj.org/operations/projections/webmerc.html:
“This is a variant of the regular Mercator projection, except that the computation is done on a sphere, using the semi-major axis of the ellipsoid.”
and
“All parameters for the projection are optional, except the ellipsoid definition, which is WGS84 for the typical use case of EPSG:3857.”
and
"+ellps=

See proj -le for a list of available ellipsoids.

Defaults to “GRS80”."

projinfo -o PROJ EPSG:3857
gives me
PROJ.4 string:
+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs

that is correct, but unfortunately not returned by GDAL’s function OSRExportToProj4() which is now deprecated. That means GRASS needs to avoid OSRExportToProj4() with GDAL 3+:
https://gdal.org/doxygen/ogr__srs__api_8h.html#ae410c0b0bb86b7d7903193828db8a1f5

-->
Warning
Use of this function is discouraged. Its behaviour in GDAL >= 3 / PROJ >= 6 is significantly different from earlier versions. In particular +datum will only encode WGS84, NAD27 and NAD83, and +towgs84/+nadgrids terms will be missing most of the time. PROJ strings to encode CRS should be considered as a legacy solution. Using a AUTHORITY:CODE or WKT representation is the recommended way.
<--
Thus GRASS should use an SRID or WKT instead of proj4 strings.

Markus M

For GDAL 3 + PROJ 6, GRASS should switch to WKT as main CRS definition format, ideally also storing a WKT definition next to PROJ_INFO. When creating locations from GDAL/OGR datasets, GRASS should go from OGR spatialreference through WKT to GRASS definition, and probably use the new PROJ6+ methods as much as possible.

Yes, good topic for
https://grasswiki.osgeo.org/wiki/GRASS_GIS_Community_Sprint_Prague_2019
:slight_smile:

When getting the projection info of a GRASS location, GRASS should then check in that order 1) WKT, 2) EPSG or other SRID, 3) GRASS native definition (which deviated from proj definition some decades ago).

That means, g.proj, r.in.gdal, r.external, v.in.ogr, v.external need updates, and the mechanism of conversion between different formats of CRS definitions in lib/proj needs to be largely rewritten.

This task for next stable GRASS release (7.10?) I guess.

Ma

–
Martin Landa
http://geo.fsv.cvut.cz/gwiki/Landa
http://gismentors.cz/mentors/landa

13

On jeudi 7 novembre 2019 19:50:33 CET Markus Metz wrote:

On Thu, Nov 7, 2019 at 6:47 PM Martin Landa <landa.martin@gmail.com> wrote:
> Hi,
>
> čt 7. 11. 2019 v 15:17 odesílatel Markus Metz
>
> <markus.metz.giswork@gmail.com> napsal:
> > PROJ 6 support is not yet finished in GRASS. Currently it sort of works

in master and relbr78, but I expect troubles because GRASS is heavily
relying on deprecated proj strings as CRS definitions. This is becoming
problematic e.g. when creating locations from GDAL/OGR datasets because
GRASS goes from OGR spatialreference through proj string to GRASS
definition.

> OSGeo4W packages are planned to be switched to GDAL3/PROJ6 soon, are
> you aware of any known bugs?

yes, EPSG:3857 is causing trouble. This is the WGS 84 / Pseudo-Mercator
projection on a sphere with radius 6378137 that needs special treatment
which no longer works with GRASS + PROJ 6 + GDAL 3.

g.proj epsg=3857 -p gives me
-PROJ_INFO-------------------------------------------------
name : WGS 84 / Pseudo-Mercator
datum : wgs84
ellps : wgs84
proj : merc
lat_ts : 0
lon_0 : 0
x_0 : 0
y_0 : 0
k : 1
wktext : defined
no_defs : defined
-PROJ_EPSG-------------------------------------------------
epsg : 3857
-PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1

which is wrong because it should be a Mercator projection on a sphere, not
an ellipsoid

Well, this is a messy topic. EPSG:3857 uses the WebMercator projection, which
is not the standard Mercator projection since it uses spherical formulas even
if the datum and ellipsoid is WGS84

So neither proj=merc + ellps=wgs84 or proj=merc + ellps=sphere are
technically correct

Nothing new with PROJ 6 for that case

projinfo -o PROJ EPSG:3857
gives me
PROJ.4 string:
+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1
+units=m +nadgrids=@null +wktext +no_defs +type=crs

Yes, this is a hack of PROJ.4 still supported by PROJ 6 to more or less get
the correct result of WebMercator using the standard Mercator, by lying on the
ellipsoid definition (using the Sphere) and applying a nadgrids=null, so that
no datum transformation between WGS84 and that Sphere is attempted !

that is correct, but unfortunately not returned by GDAL's function
OSRExportToProj4() which is now deprecated.

Hum, I do ot agree with this. This *is* still what is returned by GDAL
OSRExportToProj4() with GDAL + PROJ master :

$ gdalsrsinfo EPSG:3857

PROJ.4 : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0
+k=1 +units=m +nadgrids=@null +wktext +no_defs

And you also actually use (PROJ 6 required here)

$ gdalsrsinfo "+proj=webmerc +datum=WGS84"
(or $ projinfo "+proj=webmerc +datum=WGS84 +type=crs" -o PROJ )

PROJ.4 : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0
+k=1 +units=m +nadgrids=@null +wktext +no_defs

That means GRASS needs to avoid
OSRExportToProj4() with GDAL 3+:

In the general case, I agree with your conclusion, but not for the above one
:slight_smile:

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

14

On Thu, Nov 7, 2019 at 8:04 PM Even Rouault <even.rouault@spatialys.com> wrote:

On jeudi 7 novembre 2019 19:50:33 CET Markus Metz wrote:

On Thu, Nov 7, 2019 at 6:47 PM Martin Landa <landa.martin@gmail.com> wrote:

Hi,

čt 7. 11. 2019 v 15:17 odesílatel Markus Metz

<markus.metz.giswork@gmail.com> napsal:

PROJ 6 support is not yet finished in GRASS. Currently it sort of works

in master and relbr78, but I expect troubles because GRASS is heavily
relying on deprecated proj strings as CRS definitions. This is becoming
problematic e.g. when creating locations from GDAL/OGR datasets because
GRASS goes from OGR spatialreference through proj string to GRASS
definition.

OSGeo4W packages are planned to be switched to GDAL3/PROJ6 soon, are
you aware of any known bugs?

yes, EPSG:3857 is causing trouble. This is the WGS 84 / Pseudo-Mercator
projection on a sphere with radius 6378137 that needs special treatment
which no longer works with GRASS + PROJ 6 + GDAL 3.

g.proj epsg=3857 -p gives me
-PROJ_INFO-------------------------------------------------
name : WGS 84 / Pseudo-Mercator
datum : wgs84
ellps : wgs84
proj : merc
lat_ts : 0
lon_0 : 0
x_0 : 0
y_0 : 0
k : 1
wktext : defined
no_defs : defined
-PROJ_EPSG-------------------------------------------------
epsg : 3857
-PROJ_UNITS------------------------------------------------
unit : meter
units : meters
meters : 1

which is wrong because it should be a Mercator projection on a sphere, not
an ellipsoid

Well, this is a messy topic. EPSG:3857 uses the WebMercator projection, which
is not the standard Mercator projection since it uses spherical formulas even
if the datum and ellipsoid is WGS84

Indeed messy. Look at the WKT2 representation of projinfo -o WKT2_2018 EPSG:3857: projection on the ellipsoid “WGS 84”, not on a sphere.

So neither proj=merc + ellps=wgs84 or proj=merc + ellps=sphere are
technically correct

Nothing new with PROJ 6 for that case

projinfo -o PROJ EPSG:3857
gives me
PROJ.4 string:
+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1
+units=m +nadgrids=@null +wktext +no_defs +type=crs

Yes, this is a hack of PROJ.4 still supported by PROJ 6 to more or less get
the correct result of WebMercator using the standard Mercator, by lying on the
ellipsoid definition (using the Sphere) and applying a nadgrids=null, so that
no datum transformation between WGS84 and that Sphere is attempted !

I don’t understand “lying on the ellipsoid definition (using the Sphere)”, it should be a sphere with radius 6378137, no?

that is correct, but unfortunately not returned by GDAL’s function
OSRExportToProj4() which is now deprecated.

Hum, I do ot agree with this. This is still what is returned by GDAL
OSRExportToProj4() with GDAL + PROJ master :

$ gdalsrsinfo EPSG:3857

PROJ.4 : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0
+k=1 +units=m +nadgrids=@null +wktext +no_defs

OK, I will figure out what goes wrong in GRASS regarding the representation of EPSG:3857.

Thanks a lot for your very fast response, Even!

Markus M

And you also actually use (PROJ 6 required here)

$ gdalsrsinfo “+proj=webmerc +datum=WGS84”
(or $ projinfo “+proj=webmerc +datum=WGS84 +type=crs” -o PROJ )

PROJ.4 : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0
+k=1 +units=m +nadgrids=@null +wktext +no_defs

That means GRASS needs to avoid
OSRExportToProj4() with GDAL 3+:

In the general case, I agree with your conclusion, but not for the above one
:slight_smile:

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

15

Hi Markus,

čt 7. 11. 2019 v 19:50 odesílatel Markus Metz
<markus.metz.giswork@gmail.com> napsal:

yes, EPSG:3857 is causing trouble. This is the WGS 84 / Pseudo-Mercator projection on a sphere with radius 6378137 that needs special treatment which no longer works with GRASS + PROJ 6 + GDAL 3.

I will try to reword my question. Are you aware of any related
blockers which should avoid us packaging upcoming GRASS 7.8.1 with
GDAL3/Proj6? Thanks for clarification! Martin

--
Martin Landa
http://geo.fsv.cvut.cz/gwiki/Landa
http://gismentors.cz/mentors/landa

16

Indeed messy. Look at the WKT2 representation of projinfo -o WKT2_2018
EPSG:3857: projection on the ellipsoid "WGS 84", not on a sphere.

as does the authority says:
https://www.epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::3857

I don't understand "lying on the ellipsoid definition (using the Sphere)",
it should be a sphere with radius 6378137, no?

What I mean is that in the following
"+proj=merc +a=6378137 +b=6378137 +nadgrids=@null"
we use a clever hack to use the standard Mercator projection to simulate the
""Popular Visualisation Pseudo Mercator" by speciyfing the sphere as an
ellipsoid with the a & b parameters, and using +nadgrids=@null which is
interpreted by PROJ as "this sphere is equivalent to WGS84 for the purpose of
datum transformations". Otherwise if you transformed from another datum than
WGS84, you might get a incorrect datum transformation...

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

17

On Thu, Nov 7, 2019 at 9:39 PM Martin Landa <landa.martin@gmail.com> wrote:

Hi Markus,

čt 7. 11. 2019 v 19:50 odesílatel Markus Metz
<markus.metz.giswork@gmail.com> napsal:

yes, EPSG:3857 is causing trouble. This is the WGS 84 / Pseudo-Mercator projection on a sphere with radius 6378137 that needs special treatment which no longer works with GRASS + PROJ 6 + GDAL 3.

I will try to reword my question. Are you aware of any related
blockers which should avoid us packaging upcoming GRASS 7.8.1 with
GDAL3/Proj6? Thanks for clarification! Martin

No, that’s the only bug I found so far.

Markus M

18

Hi,

pá 8. 11. 2019 v 8:22 odesílatel Markus Metz
<markus.metz.giswork@gmail.com> napsal:

> I will try to reword my question. Are you aware of any related
> blockers which should avoid us packaging upcoming GRASS 7.8.1 with
> GDAL3/Proj6? Thanks for clarification! Martin

No, that's the only bug I found so far.

thanks for clarification! In that case let's release 7.8.1 ASAP in
order to allow OSGeo4W switch to GDAL3/PROJ6.

Ma

--
Martin Landa
http://geo.fsv.cvut.cz/gwiki/Landa
http://gismentors.cz/mentors/landa

19

On Fri, Nov 8, 2019 at 10:15 AM Martin Landa <landa.martin@gmail.com> wrote:

Hi,

pá 8. 11. 2019 v 8:22 odesílatel Markus Metz
<markus.metz.giswork@gmail.com> napsal:

I will try to reword my question. Are you aware of any related
blockers which should avoid us packaging upcoming GRASS 7.8.1 with
GDAL3/Proj6? Thanks for clarification! Martin

No, that’s the only bug I found so far.

thanks for clarification! In that case let’s release 7.8.1 ASAP in
order to allow OSGeo4W switch to GDAL3/PROJ6.

the EPSG:3857 (Web Mercator) bug has been fixed in master 4167542. Is there still time to backport to relbr78?

It is pure chance that the previous method worked, this is a fix for the GRASS hack to deal with “the hacky way of encodign webmerc in GDAL WKT1” to cite PROJ 6.

All this trouble because +proj=webmerc exists only since PROJ 5.1.0, and because even GDAL 3 and PROJ 6 (at least PROJ 6 knows +proj=webmerc) still report “+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs” as a workaround.

Markus M

20

On Fri, Nov 8, 2019 at 7:39 PM Markus Metz
<markus.metz.giswork@gmail.com> wrote:

On Fri, Nov 8, 2019 at 10:15 AM Martin Landa <landa.martin@gmail.com> wrote:
>
> Hi,
>
> pá 8. 11. 2019 v 8:22 odesílatel Markus Metz
> <markus.metz.giswork@gmail.com> napsal:
> > > I will try to reword my question. Are you aware of any related
> > > blockers which should avoid us packaging upcoming GRASS 7.8.1 with
> > > GDAL3/Proj6? Thanks for clarification! Martin
> >
> > No, that's the only bug I found so far.
>
> thanks for clarification! In that case let's release 7.8.1 ASAP in
> order to allow OSGeo4W switch to GDAL3/PROJ6.

the EPSG:3857 (Web Mercator) bug has been fixed in master 4167542.

Excellent, thanks a lot.

Is there still time to backport to relbr78?

Yes. I have done so right now.

It is pure chance that the previous method worked, this is a fix for the GRASS hack to deal with "the hacky way of encodign webmerc in GDAL WKT1" to cite PROJ 6.

All this trouble because +proj=webmerc exists only since PROJ 5.1.0, and because even GDAL 3 and PROJ 6 (at least PROJ 6 knows +proj=webmerc) still report "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs" as a workaround.

Markus M

It is still a bit foggy :slight_smile:

markusN