Different results from m.proj in different locations

[wenzeslaus]

I wanted to get the latitude and longitude extent for a map, so I applied m.proj to the map extent. However, I've noticed that I get different results for the same map when running m.proj in NC sample dataset (Lambert Conformal Conic) and in my custom location (WGS 84 Pseudo-Mercator).

This is related with the #2100 ticket but I'm not sure if this is really the same. It it something different at least from user point of view.

I cannot tell about the sample elevation map but for some other data it seemed to me that the LL from Merkator is correct and the the other are wrong.

How to reproduce

To reproduce, run in NC sample dataset:

eval $(r.info map=elevation -g)
m.proj -o -d input=- separator=space <<EOF
$east $north
$west $south
EOF

I got:

-78.60830318 35.80918894
-78.77462049 35.68792712

Then create a new location, I've used EPSG:3857 (WGS 84 Pseudo-Mercator).

In new location run:

eval $(r.proj -g input=elevation location=nc_spm_08_grass7 mapset=PERMANENT)
m.proj -o -d input=- separator=space <<EOF
$e $n
$w $s
EOF

I got:

-78.60835872 35.80956421
-78.77456513 35.68755254

Which is something close but different. Compare the points at OSM (you can change the mlat and mlon in URL):

• ​http://www.openstreetmap.org/?mlat=35.80918894&mlon=-78.60830318#map=17/35.80850/-78.60932
• ​http://www.openstreetmap.org/?mlat=35.80956421&mlon=-78.60835872#map=17/35.80850/-78.60932

Relation to the r.proj bounding box ticket

When I apply the (first) patch from the possibly related ticket #2100, I get different results but still does not match with the others:

-78.60830339 35.80960928
-78.77462031 35.68750748

OSM link:

• ​http://www.openstreetmap.org/?mlat=35.80960928&mlon=-78.60830339#map=17/35.80850/-78.60932

What g.region -l says

During writing of this ticket I also find the g.region -l which I can actually use instead of m.proj, however I still get the different results.

LCC:

north-east corner: long: 78:36:29.891436W lat: 35:48:33.080191N
south-west corner: long: 78:46:28.633781W lat: 35:41:16.537646N

Mercator:

north-east corner: long: 78:36:30.091395W lat: 35:48:34.431172N
south-west corner: long: 78:46:28.434454W lat: 35:41:15.189156N

Creating the new location

I suggest to create the new location using Location wizard instead of g.proj because I've tried to write a testing code for this but the g.proj didn't worked for me. The commands:

g.proj epsg=3857 location=test_location_3857
g.mapset mapset=PERMANENT location=test_location_3857

create and switch the location but than the r.proj fails with:

WARNING: pj_transform() failed: failed to load datum shift file
...

If I'm using it right this might be another ticket. I've used it recently through Python grass.script.core.create_location() function and it worked well.

Warning about the LL format

When trying to do something related to this ticket and with LL either by hand or programmatically be warned about different order and formats (speaking about GRASS only). I spent some time with this. Maybe this inconsistency is topic for different ticket.

[neteler]

Probably EPSG:3857 (WGS 84 Pseudo-Mercator) comes with too many problems anyway:

• ​http://earth-info.nga.mil/GandG/wgs84/web_mercator/index.html
• ​http://blog.geogarage.com/2014/09/advisory-notice-on-web-mercator.html
• ​https://blogs.esri.com/esri/arcgis/2014/09/25/what-does-the-nga-web-mercator-advisory-mean-for-esri-defense-and-intelligence-users/

[neteler]

I just made a new test here:

# part 1:
GRASS 7.4.0svn (nc_spm_08_grass7):~/grassdata > eval $(r.info map=elevation -g)
GRASS 7.4.0svn (nc_spm_08_grass7):~/grassdata > m.proj -o -d input=- separator=space <<EOF
> $east $north
> $west $south
> EOF
-78.60830318 35.80918894 0.00000000
-78.77462049 35.68792712 0.00000000

g.version -reg
version=7.4.0svn
date=2018
revision=r72020
build_date=2018-01-04
build_platform=x86_64-pc-linux-gnu
build_off_t_size=8
libgis_revision=70829 
libgis_date="2017-04-04 09:43:02 +0200 (Tue, 04 Apr 2017) "
proj4=4.9.3
gdal=2.1.3
geos=3.6.1
sqlite=3.20.1

Then:

# part 2:
GRASS 7.4.0svn (TEST):~/grassdata > g.proj -w
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["unnamed",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,"inf"]],
        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]]

eval $(r.proj -g input=elevation location=nc_spm_08_grass7 mapset=PERMANENT)

m.proj -o -d input=- separator=space <<EOF
> $e $n
> $w $s
> EOF
-78.60830318 35.80960938 0.00000000
-78.77462049 35.68750730 0.00000000

# expected (?)
-78.60830318 35.80918894 0.00000000
-78.77462049 35.68792712 0.00000000

[mmetz]

Replying to neteler:

I just made a new test here:

# part 1:
GRASS 7.4.0svn (nc_spm_08_grass7):~/grassdata > eval $(r.info map=elevation -g)
GRASS 7.4.0svn (nc_spm_08_grass7):~/grassdata > m.proj -o -d input=- separator=space <<EOF
> $east $north
> $west $south
> EOF
-78.60830318 35.80918894 0.00000000
-78.77462049 35.68792712 0.00000000

g.version -reg
version=7.4.0svn
date=2018
revision=r72020
build_date=2018-01-04
build_platform=x86_64-pc-linux-gnu
build_off_t_size=8
libgis_revision=70829 
libgis_date="2017-04-04 09:43:02 +0200 (Tue, 04 Apr 2017) "
proj4=4.9.3
gdal=2.1.3
geos=3.6.1
sqlite=3.20.1

Then:

# part 2:
GRASS 7.4.0svn (TEST):~/grassdata > g.proj -w
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["unnamed",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,"inf"]],
        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]]

eval $(r.proj -g input=elevation location=nc_spm_08_grass7 mapset=PERMANENT)

m.proj -o -d input=- separator=space <<EOF
> $e $n
> $w $s
> EOF
-78.60830318 35.80960938 0.00000000
-78.77462049 35.68750730 0.00000000

# expected (?)
-78.60830318 35.80918894 0.00000000
-78.77462049 35.68792712 0.00000000

No, not expected: first you reproject the upper-right corner and the lower-left corner with m.proj. Then you run r.proj -g which reprojects all outer corners of all border cells to figure out the new extents. This test assumes that a rectangle stays a rectangle after reprojection and that north of the upper-left and upper-right corner are the same after reprojection (not to mention northing of points in between the upper-left and upper-right corner). This assumption is wrong.

EPSG:3857 (WGS 84 Pseudo-Mercator) is now handled reasonably well in 7.4 and 7.5, as long as the location is newly created to get correct projection information.

[martinl]

Can we close this issue?

[neteler]

Ticket retargeted after milestone closed

[neteler]

Replying to martinl:

Can we close this issue?

Ping wenzeslaus :)

[neteler]

Fixed in 7.4.x and later, closing.