error in projection definition

[vincent]

Setting up a new location via location wizard results in an error when it's done by specifying the epsg code 3857.

This (not recognized) geodetic system uses a spherical ellipsoid (+a=+b), but the created location's projection parameters show a flattened ellipsoid (+rf=298,...).

[neteler]

I cannot reproduce the problem with 6.4.svn r60122. Screenshot attached. My system is:

System Info
GRASS version: 6.4.4svn
GRASS SVN Revision: 60122M
GDAL/OGR: 1.10.1
PROJ4: Rel. 4.8.0, 6 March 2012
Python: 2.7.5
wxPython: 2.8.12.0
Platform: Linux-3.14.3-200.fc20.x86_64-x86_64-with-fedora-20-Heisenbug

I also tested on the same machine GRASS 7.1.svn:

grass71 -c epsg:3857 ~/grassdatum/test
Cleaning up temporary files...
Creating new GRASS GIS location/mapset...
...
Welcome to GRASS 7.1.svn (r60176M)
...

GRASS 7.1.svn (test):~ > g.proj -w
PROJCS["Mercator",
    GEOGCS["wgs84",
        DATUM["WGS_1984",
            SPHEROID["WGS_1984",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

All fine, perhaps your PROJ.4 version is outdated?

[neteler]

Successfully defining the virtual globe epsg code 3857

[hamish]

I suspect the fact that epsg 3857 is a new entry for the epsg file (replaces 900973 from extra.esri) may have something to do with it. How does it act for codes which do not exist?

Vincent, what version of proj4 and gdal are you using? what platform? did you try to then specify the definition manually by +proj4 terms or by manually choosing ellips to be sphere, etc?

Hamish

ps- please try to only use pseudo-mercator if data import/export demands it, it is the cause of much pain. :-)

[vincent]

Replying to hamish:

Vincent, what version of proj4 and gdal are you using? what platform?

Quite old ones :

proj4 : Rel. 4.7.1, 23 September 2009

gdal : GDAL 1.9.0, released 2011/12/29

uname --all : Linux pc1 3.2.0-4-amd64 Debian 3.2.57-3 x86_64 GNU/Linux

did you try to then specify the definition manually by +proj4 terms or by manually choosing ellips to be sphere, etc?

Yes, passing the proj4 chain gives the right projection, see this ​thread for history.

ps- please try to only use pseudo-mercator if data import/export demands it, it is the cause of much pain. :-)

I keep it in mind from then...

V.

[vincent]

Replying to neteler:

I cannot reproduce the problem with 6.4.svn r60122. Screenshot attached.

Markus, sorry I missed your reply 6 weeks ago. I don't really understand your attachment : to me it seems you have the same issue i.e. a non-null rf parameter.

V.

[hamish]

note what shows up on the loc'n wizard summary screen is *not* the terms used to create the location! Check your $MAPSET/PERMANENT/PROJ_INFO file (g.proj -p) to see what really made it in.

vincent:

proj4 : Rel. 4.7.1, 23 September 2009

gdal : GDAL 1.9.0, released 2011/12/29

Debian/wheezy then? same here. the epsg file which comes with the proj-data package there does include epsg 3857, as does wkt/csv versions of the epsg database in grass-core, libgdal1, and libgeotiff-epsg packages. So all is new enough.

I'm not sure how grass's /usr/lib/grass64/etc/ogr_csv/pcs.csv handles the @null part though, or if that is even used here, it may just be the version picked up by the OGR library functions which matters.

are you using the standard grass package from debian wheezy? (ver 6.4.2-2)

Hamish

[vincent]

Replying to hamish:

Debian/wheezy then?

yes

I'm not sure how grass's /usr/lib/grass64/etc/ogr_csv/pcs.csv handles the @null part though, or if that is even used here, it may just be the version picked up by the OGR library functions which matters.

are you using the standard grass package from debian wheezy? (ver 6.4.2-2)

no, ver 6.4.4svn (last compiled about 2 weeks ago)

Vincent.

[neteler]

As per r68131 (trunk) the PROJCS (coordinate system name) is now maintained properly:

g.proj epsg=3857 datumtrans=-1 -w 
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["wgs84",
        DATUM["WGS_1984",
            SPHEROID["WGS_1984",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

[vincent]

To me the issue is still the same, except that it may not be a grass-specific one. On a fresh grass7 (r68135) install :

• Create a new pseudo-mercator location with loc'n wiz by selecting epsg code 3857
• Pick "0" from datum transformations list ("do not apply any datum transformations")
• Then execute :

m.proj -o -d coordinates=720000,5740000
6.46787005|45.93986903|0.00000000

• regarding this result execute :

echo "6.46787005 45.93986903" | cs2cs +init=epsg:4326 -f "%.8f" +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
720000.00048307	5770718.42069077 0.00000000

There is a north-south error of approx. 30 km; It looks like m.proj sends "+rf=298.25..." to cs2cs, where it shoud not...

[martinl]

Replying to vincent:

There is a north-south error of approx. 30 km; It looks like m.proj sends "+rf=298.25..." to cs2cs, where it shoud not...

right, g.proj -jf returns

+proj=merc +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +no_defs +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000 +to_meter=1

while proj.4 epsg file says

<3857> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs <>

[neteler]

Re: EPSG 3857 Pseudomercator: there are both WGS84 and spherical in use...: see ​https://en.wikipedia.org/wiki/Web_Mercator#Spherical_or_ellipsoidal.3F

[martinl]

See also related GDAL tickets: ​https://trac.osgeo.org/gdal/ticket/3962, ​https://trac.osgeo.org/gdal/ticket/5924

[martinl]

GRASS is using GDAL library to convert EPSG code to projection definition, so that's the reason why GRASS reports WGS84 PseudoMercator even PROJ.4 uses Spherical PseudoMercator for EPSG 3857. See source:grass/trunk/general/g.proj/input.c#L171​. The issue should be solved on GDAL side I would say.

[vincent]

Or should GRASS be aware of the parameter EXTENSION ?

EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"]

[martinl]

Replying to vincent:

Or should GRASS be aware of the parameter EXTENSION ?

EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"]

I think GDAL should be aware of that not GRASS (which is using GDAL function to translate EPSG code).

[neteler]

FYI: New details about lib/proj/ have been added to

​https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/README.txt

along with the proj related CSV removal in r69211. Please test again.

[vincent]

still the same issue with r69289

[mlennert]

As discussed during the code sprint, the issue is because of the very weird nature of this "projection" system:

As mentioned on ​epsg.io:

"Uses spherical development of ellipsoidal coordinates. Relative to WGS 84 / World Mercator (CRS code 3395) errors of 0.7 percent in scale and differences in northing of up to 43km in the map (equivalent to 21km on the ground) may arise."

It uses the WGS 84 datum (which normally implies the WGS 84 ellipsoid), but actually develops the coordinates on a sphere. AFAIU, in GRASS, we have not foreseen the option to declare one datum, but using a different ellipsoid than the one normally linked to that datum.

The proj.4 definition of this projection system seems to be somewhat of a hack.

In the discussions, one idea that came up was to provide alternative datum transformation parameters for this system, one on the sphere and one on the wgs84 ellipsoid.

The most reasonable conclusion, however, was to not use this system at all because of it's weird nature. Several official administrations forbid its use...

[vincent]

Replying to mlennert:

As discussed during the code sprint, the issue is because of the very weird nature of this "projection" system [...]

I fully agree with these comments

The most reasonable conclusion, however, was to not use this system at all because of it's weird nature. Several official administrations forbid its use...

Ok bu sooner or later, users may have to manage data in this "projection" system, and they might be puzzled (angry) to see it does not work despite a proper location definition via the epsg code.

In the discussions, one idea that came up was to provide alternative datum transformation parameters for this system, one on the sphere and one on the wgs84 ellipsoid.

If there is no solution to automatically "patch" this issue, meanwhile we could edit a wiki page describing a solution (like passing a proj.4 chain) to the location's creation.

[mlennert]

Replying to vincent:

Replying to mlennert:

As discussed during the code sprint, the issue is because of the very weird nature of this "projection" system [...]

I fully agree with these comments

The most reasonable conclusion, however, was to not use this system at all because of it's weird nature. Several official administrations forbid its use...

Ok bu sooner or later, users may have to manage data in this "projection" system, and they might be puzzled (angry) to see it does not work despite a proper location definition via the epsg code.

In the discussions, one idea that came up was to provide alternative datum transformation parameters for this system, one on the sphere and one on the wgs84 ellipsoid.

If there is no solution to automatically "patch" this issue, meanwhile we could edit a wiki page describing a solution (like passing a proj.4 chain) to the location's creation.

At this stage, I'm not sure what the solution would be. I created a location with the proj.4 chain you gave:

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs <>

and then attempted to project a file from EPSG 31370 (Belgian Lambert 1972) to that location, but I got the following error:

WARNING: pj_transform() a échoué: failed to load datum shift file
ERROR: Unable to re-project vector map <rues@urbis> from <Belgique31370>

[neteler]

Can you please try the debugging tools described in ​https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/README.txt

--> Testing: comparing results with PROJ_DEBUG=ON CPL_DEBUG=ON ?

[mlennert]

Replying to neteler:

Can you please try the debugging tools described in ​https://trac.osgeo.org/grass/browser/grass/trunk/lib/proj/README.txt

--> Testing: comparing results with PROJ_DEBUG=ON CPL_DEBUG=ON ?

> v.proj location=Belgique31370 mapset=urbis input=rues
Reprojecting primitives ...
pj_open_lib(/data/home/mlennert/SRC/GRASS/grass_trunk/dist.x86_64-pc-linux-gnu/etc/proj/nad/@null): call fopen(/data/home/mlennert/SRC/GRASS/grass_trunk/dist.x86_64-pc-linux-gnu/etc/proj/nad/@null) - failed

WARNING: pj_transform() failed: failed to load datum shift file
ERROR: Unable to re-project vector map <rues@urbis> from <Belgique31370>
GDAL: In GDALDestroy - unloading GDAL shared library.

The problem seems to be the +nadgrids=@null in the proj.4 string.

[neteler]

Can you restore it locally from the last rev and test?

[mlennert]

Replying to neteler:

Can you restore it locally from the last rev and test?

Restore what ?

At the code sprint, Even said that this nadgrids=@null was a hack that he didn't think was the best approach...

[hellik]

Replying to mlennert:

Replying to neteler:

Can you restore it locally from the last rev and test?

Restore what ?

At the code sprint, Even said that this nadgrids=@null was a hack that he didn't think was the best approach...

Just some recent experiences with EPSG:3857:

• created location with EPSG:3857
• some coordinates processing along rivers in this location
• feeded these coordinates to gdal_translate to download some tiles of a WMTS service which delivers data in EPSG:3857
• then gdalwarp'ed EPSG:3857-data to an Austrian local projection
• reprojected rasters are ok

[martinl]

No activity.

[mmetz]

Fixed in trunk r71523