Different results from m.proj in different locations
|Reported by:||wenzeslaus||Owned by:|
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
-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
-78.60835872 35.80956421 -78.77456513 35.68755254
Which is something close but different. Compare the points at OSM (you can change the
mlon in URL):
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
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.
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
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.