Opened 8 years ago

Closed 3 years ago

#2229 closed defect (fixed)

error in projection definition

Reported by: vincent Owned by: grass-dev@…
Priority: normal Milestone: 7.0.7
Component: Projections/Datums Version: svn-releasebranch70
Keywords: pseudomercator Cc:
CPU: Unspecified Platform: Unspecified

Description

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,...).

Attachments (1)

grass644_loc_wiz_epsg3857.png (212.5 KB) - added by neteler 8 years ago.
Successfully defining the virtual globe epsg code 3857

Download all attachments as: .zip

Change History (32)

comment:1 Changed 8 years ago by 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?

Changed 8 years ago by neteler

Successfully defining the virtual globe epsg code 3857

comment:2 Changed 8 years ago by 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. :-)

comment:3 in reply to:  2 Changed 8 years ago by 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.

comment:4 in reply to:  1 Changed 8 years ago by 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.

comment:5 Changed 8 years ago by 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

comment:6 in reply to:  5 Changed 8 years ago by 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.

comment:7 Changed 6 years ago by neteler

Milestone: 6.4.47.0.4
Version: unspecifiedsvn-releasebranch70

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]]

comment:8 Changed 6 years ago by 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...

comment:9 Changed 6 years ago by martinl

Milestone: 7.0.47.0.5

comment:10 in reply to:  8 Changed 5 years ago by 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 <>

comment:11 Changed 5 years ago by 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

comment:13 Changed 5 years ago by 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.

comment:14 Changed 5 years ago by martinl

Keywords: pseudomercator added

comment:15 in reply to:  13 ; Changed 5 years ago by 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"]

comment:16 in reply to:  15 Changed 5 years ago by 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).

comment:17 Changed 5 years ago by 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.

comment:18 Changed 5 years ago by vincent

still the same issue with r69289

comment:19 Changed 5 years ago by 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...

comment:20 in reply to:  19 ; Changed 5 years ago by 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.

comment:21 in reply to:  20 Changed 5 years ago by 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>

comment:22 Changed 5 years ago by 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 ?

comment:23 in reply to:  22 Changed 5 years ago by 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.

comment:24 Changed 5 years ago by neteler

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

comment:25 in reply to:  24 ; Changed 5 years ago by 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...

comment:26 in reply to:  25 Changed 5 years ago by 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

comment:27 Changed 5 years ago by neteler

Milestone: 7.0.57.0.6

comment:28 Changed 4 years ago by neteler

Milestone: 7.0.67.0.7

comment:29 Changed 3 years ago by martinl

Resolution: wontfix
Status: newclosed

No activity.

comment:30 Changed 3 years ago by mmetz

Resolution: wontfix
Status: closedreopened

Fixed in trunk r71523

comment:31 Changed 3 years ago by mmetz

Resolution: fixed
Status: reopenedclosed
Note: See TracTickets for help on using tickets.