[GRASS-dev] PROJ6 support in GRASS

[Markus Metz]

On Thu, Nov 7, 2019 at 6:47 PM Martin Landa <landa.martin at gmail.com> wrote:
>
> Hi,
>
> čt 7. 11. 2019 v 15:17 odesílatel Markus Metz
> <markus.metz.giswork at 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=<value>

    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
-->
WarningUse 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
> :-)
>
> > 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20191107/a8a72b64/attachment-0001.html>