[GRASS-user] projection issue

Vincent_Bain · March 24, 2014, 10:33am

Hello,

the trouble I am currently experiencing is dealing partly with cs2cs,
but partly with grass, so I post my question on this list...

Being within a Mercator location, defined with the location wizard
(calling epsg:3857), I type :

        echo "6.46311951 45.78037144" | m.proj -idg
        input parameters=[+proj=longlat +datum=WGS84]
        output parameters=[+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]

Which results in :

719471.17278927|5714587.81400429|0.00000000

The thing is this set of geographic coordinates should return something
else...

        echo "6.46311951 45.78037144" | cs2cs +proj=longlat -f "%.8f"
        +to +init=epsg:3857

gives :

719471.17278927 5745223.17292760 0.00000000

I guess this is not a m.proj issue, because passing projection argument
directly to cs2cs gives the same error :

curproj=`eval g.proj -jf`

        echo "6.46311951 45.78037144" | cs2cs +proj=longlat -f "%.8f"
        +to $curproj

returns :

719471.17278927 5714587.81400429 0.00000000

I can't figure what is going wrong. Would anyone have an idea ?
Thanks,

Vincent

Vincent_Bain · March 24, 2014, 12:10pm

Comment on my previous message :

if we have a closer look at parameters, we notice a difference between
the current (grass location) projection parameters and epsg:3857
parameters :

Current :
        a (semi-major axis) = 6378137
        rf (flattening factor inverse) =298.257223563
am I wrong or as a consequence
        b (semi-minor axis) = 6356752 ?

epsg:3857 :
        a (semi-major axis) = 6378137
        b (semi-minor axis) = 6378137 (a sphere ?)

Now if I run cs2cs with these distinct params :

         echo "6.46311951 45.78037144" | 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
returns
        719471.17278927 5745223.17292760 0.00000000

and
        echo "6.46311951 45.78037144" | cs2cs +init=epsg:4326 -f "%.8f"
        +to +proj=merc +a=6378137 +b=6356752 +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
returns
        719471.17278927 5714587.36406099 0.00000000
that is what m.proj initially said.

It looks as if grass was taking into account the flatten-shape of the
ellipsoid, where epsg registry considers it is a sphere.

Sorry for the noise, but I would actually enjoy understanding where I am
wrong ;-(

Bye,
Vincent

mlennert · March 24, 2014, 12:49pm

On 24/03/14 13:10, Vincent Bain wrote:

Comment on my previous message :

if we have a closer look at parameters, we notice a difference between
the current (grass location) projection parameters and epsg:3857
parameters :

Current :
         a (semi-major axis) = 6378137
         rf (flattening factor inverse) =298.257223563
am I wrong or as a consequence
         b (semi-minor axis) = 6356752 ?

epsg:3857 :
         a (semi-major axis) = 6378137
         b (semi-minor axis) = 6378137 (a sphere ?)

Now if I run cs2cs with these distinct params :

          echo "6.46311951 45.78037144" | 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
returns
         719471.17278927 5745223.17292760 0.00000000

and
         echo "6.46311951 45.78037144" | cs2cs +init=epsg:4326 -f "%.8f"
         +to +proj=merc +a=6378137 +b=6356752 +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
returns
         719471.17278927 5714587.36406099 0.00000000
that is what m.proj initially said.

It looks as if grass was taking into account the flatten-shape of the
ellipsoid, where epsg registry considers it is a sphere.

Sorry for the noise, but I would actually enjoy understanding where I am
wrong ;-(

You're not wrong. The difference comes from the GRASS file /etc/ogr_csv/pcs.csv in which there is a specific rule set for 3857.

Don't know why, though.

Moritz

Vincent_Bain · March 24, 2014, 1:03pm

Le lundi 24 mars 2014 à 13:49 +0100, Moritz Lennert a écrit :

You're not wrong. The difference comes from the GRASS file
/etc/ogr_csv/pcs.csv in which there is a specific rule set for 3857.

Don't know why, though.

Moritz

Thank you Moritz for your reply. Yes, I finally tried to make a new
location with a proj.4 chain rather than the epsg code picker, and it
works fine.

The source of this "error" is very surprising, as much as the few
reactions it's been triggering !

Bye,
Vincent

mlennert · March 24, 2014, 2:00pm

On 24/03/14 14:03, Vincent Bain wrote:

Le lundi 24 mars 2014 à 13:49 +0100, Moritz Lennert a écrit :

You're not wrong. The difference comes from the GRASS file
/etc/ogr_csv/pcs.csv in which there is a specific rule set for 3857.

Don't know why, though.

Moritz

Thank you Moritz for your reply. Yes, I finally tried to make a new
location with a proj.4 chain rather than the epsg code picker, and it
works fine.

The source of this "error" is very surprising, as much as the few
reactions it's been triggering !

Well, not sure this is really an "error", or rather an issue about how exactly to define this "pseudo-mercator" projection system.

In any case, it seems to me to be more a gdal issue than a grass issue.

I guess https://trac.osgeo.org/gdal/ticket/3962 is related to this.

Moritz

Vincent_Bain · March 24, 2014, 2:18pm

Le lundi 24 mars 2014 à 15:00 +0100, Moritz Lennert a écrit :

I guess https://trac.osgeo.org/gdal/ticket/3962 is related to this.

interesting discussions

confusion/debate on this particular point is weird : there should be
only one way to choose between a sphere and an ellipsoid ?

As far as I'm concerned, I keep record of the hand-made location
definition for further work with epsg:3857

Yours,
V.

Hamish · March 24, 2014, 8:24pm

Vincent:

the trouble I am currently experiencing is dealing partly with
cs2cs, but partly with grass, so I post my question on this
list...

Being within a Mercator location, defined with the location
wizard (calling epsg:3857), I type :

...

I can't figure what is going wrong. Would anyone have an idea ?

Hi,

First of all, avoid using Google's "spherical mercator" projection for anything other than necessary data import/export if you can possibly help it. But sometimes we don't have a choice, so..

Google's spherical Mercator projection (formerly the unofficial
epsg-ish code 900913 from the esri.extra file) isn't really WGS84, it just borrows the major axis of the Earth used in WGS84 for both the major and minor axes of the ellipsoid, making it a sphere larger than Earth really is. So it's just a simple Mercator on a sphere, which happens to have a particular radius.

Different ellipsoids (and a sphere being a type of ellipsoid) squash the Earth north-south, so the northing value changes the most. Since Mercator doesn't deform east-west (lines of longitude are all parallel and perfectly vertical), in that projection changing the ellipsoid does not change the easting value at all.

Next thing to know is the +nadgrids=@null hack to get around the datum transform. See
http://trac.osgeo.org/proj/wiki/GenParms#ThenullGrid

and this basically explains the rest:

http://trac.osgeo.org/proj/wiki/FAQ#ChangingEllipsoidWhycantIconvertfromWGS84toGoogleEarthVirtualGlobeMercator

fyi, +wktext has to do with "well known text" .prj file text, and requests to software in-the-know that the non-standard "+nadgrids=@null" projection term gets recorded when you export in that way. The default (without +wktext) is to only export official spec. terms, for compatibility with brittle software which might have to import it later and can't deal with anything beyond the official WKT spec (or ESRI's diversions from it, but there's another g.proj flag for that..).

good luck,
Hamish

Vincent_Bain · March 24, 2014, 9:47pm

Le lundi 24 mars 2014 à 13:24 -0700, Hamish a écrit :

Hi,

First of all, avoid using Google's "spherical mercator" projection for anything other than necessary data import/export if you can possibly help it. But sometimes we don't have a choice, so..

I wish I could never have heard about it... till I had to cope with
french geographical web services : their wmts service is based on this
projection system, so I need to define a spherical mercator location to
correctly retrieve tilesets.

Google's spherical Mercator projection (formerly the unofficial
epsg-ish code 900913 from the esri.extra file) isn't really WGS84, it just borrows the major axis of the Earth used in WGS84 for both the major and minor axes of the ellipsoid, making it a sphere larger than Earth really is. So it's just a simple Mercator on a sphere, which happens to have a particular radius.

Different ellipsoids (and a sphere being a type of ellipsoid) squash the Earth north-south, so the northing value changes the most. Since Mercator doesn't deform east-west (lines of longitude are all parallel and perfectly vertical), in that projection changing the ellipsoid does not change the easting value at all.

Yes, in my case I noticed a n-s 30 km gap

Next thing to know is the +nadgrids=@null hack to get around the datum transform. See
http://trac.osgeo.org/proj/wiki/GenParms#ThenullGrid

and this basically explains the rest:

http://trac.osgeo.org/proj/wiki/FAQ#ChangingEllipsoidWhycantIconvertfromWGS84toGoogleEarthVirtualGlobeMercator

I believed the nadgrid was a shifting matrix applying to ellipsoid
center (in cartesian geocentric coordinates). Does this mean
+nadgrids=@null implicitly tells cs2cs : keep the from-ellipsoid (ie
ignore the to-sphere), and shift it (of 0,0,0) ?

good luck,
Hamish

Your explanations and links on this topic were very welcome !
Thank you,

Vincent

Glynn_Clements1 · March 26, 2014, 3:48pm

Vincent Bain wrote:

I believed the nadgrid was a shifting matrix applying to ellipsoid
center (in cartesian geocentric coordinates). Does this mean
+nadgrids=@null implicitly tells cs2cs : keep the from-ellipsoid (ie
ignore the to-sphere), and shift it (of 0,0,0) ?

+nadgrids=@null tells it not to perform any datum transformation.

The default behaviour is to perform a matrix-based datum
transformation based upon the source and destination datums.

The nadgrids parameter overrides this. Normally, it is used to specify
a data file containing the datum-transformed positions of a set of
points; this allows for a more accurate transformation than a
matrix-based transformation.

nadgrids=@null is a "hack" which allows the datum transformation to be
suppressed entirely. This is required when the lat/lon coordinates are
based upon an ellipsoid but the cartographic transformation is based
upon a sphere, as is the case for the Google Earth projection.

Without that parameter, whatever ellipsoid was used for the
cartographic transformation would also be used for the datum
transformation.

--
Glynn Clements <glynn@gclements.plus.com>

Vincent_Bain · March 30, 2014, 8:37am

Le mercredi 26 mars 2014 à 15:48 +0000, Glynn Clements a écrit :

Thank you Glynn for the clarification

The default behaviour is to perform a matrix-based datum
transformation based upon the source and destination datums.

The nadgrids parameter overrides this. Normally, it is used to specify
a data file containing the datum-transformed positions of a set of
points; this allows for a more accurate transformation than a
matrix-based transformation.

yes, here in France we often need to switch between two geodetic
systems, the elder being based on the Clarke 1880 ellipsoid. One can
choose between a standard datum shift and a grid-based transformation.
Many gis users still use the old geodetic system and don't pay attention
to the transformation method when it comes to retrieve data from a
gps...
In the place I live in the Alps, the error observed between standard-
and grid-based transformation is about 40 meters in EW direction.

VB