Opened 13 years ago

Last modified 13 years ago

#3672 new bug

region error when creating a GRASS location?

Reported by: lutra Owned by:
Priority: major: does not work as expected Milestone: Version 1.7.0
Component: GRASS Version: Trunk
Keywords: Cc: pcav, neteler, rblazek
Must Fix for Release: Yes Platform: All
Platform Version: Awaiting user input: no

Description

Hi, I have seen this on both Linux and Windows.

If I open a "old" mapset (don't remember exactly when it was created), that was created with EPSG 40000, g.region -p will return as expected

projection: 99 (Transverse Mercator)
zone:       0
datum:      rome40
ellipsoid:  international
north:      4834913.78536118
south:      4803364.45415656
west:       1650370.47056664
east:       1689398.62729587
nsres:      24.99947005
ewres:      25.00202225
rows:       1262
cols:       1561
cells:      1969982

if I create new location with EPSG 40000 then g.region -p will return

projection: 99 (Transverse Mercator)
zone:       0
datum:      towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
ellipsoid:  international
north:      5336010
south:      4004460
west:       1240010
east:       2430820
nsres:      1191.01073345
ewres:      1192.002002
rows:       1118
cols:       999
cells:      1116882

that I don't know if it is correct.

Now if I create a location in EPSG 3003 (that should be 40000==3003+towgs) the result of g.region -p is

projection: 99 (Transverse Mercator)
zone:       0
datum:      ** unknown (default: WGS84) **
ellipsoid:  international
north:      5336090
south:      4004520
west:       1239980
east:       2430820
nsres:      1191.02862254
ewres:      1190.84
rows:       1118
cols:       1000
cells:      1118000

that seems wrong.

The same result happens also with other projection, for example epsg 3763 (or all other CRS system I use normally)

projection: 99 (Transverse Mercator)
zone:       0
datum:      ** unknown (default: WGS84) **
ellipsoid:  grs80
north:      -497.84044931
south:      -292994.94849781
west:       -98340.22283752
east:       105157.76515519
nsres:      4.99995056
ewres:      4.99995056
rows:       58500
cols:       40700
cells:      2380950000

Change History (5)

comment:1 by lutra, 13 years ago

FYI

creating the location (ex. epsg 3763) with wingrass (6.4.0) returns the expected g.region result

g.region -p                                                                     
projection: 99 (Transverse Mercator)
zone:       0
datum:      etrs89
ellipsoid:  grs80

opening with wingrass a different location created with QGIS (with the same epsg) r.region returns

g.region -p                                                                     
projection: 99 (Transverse Mercator)
zone:       0
datum:      ** unknown (default: WGS84) **
ellipsoid:  international

comment:2 by neteler, 13 years ago

The reason will be that several datums are associated to EPSG 3003:

g.proj -c epsg=3003 datumtrans=-1 location=loc_epsg_3003
---
1
Used in whole rome40 region
towgs84=-225.000,-65.000,9.000
Default 3-Parameter Transformation (May not be optimum for older datums; use this only if no more appropriate options are available.)
---
2
Used in Italy (Peninsular Part)
towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
Accuracy 3-4m
---
3
Used in Italy (Sardinia)
towgs84=-168.6,-34.0,38.6,-0.374,-0.679,-1.379,-9.48
Accuracy 3-4m
---
4
Used in Italy (Sicily)
towgs84=-50.2,-50.4,84.8,-0.690,-2.012,0.459,-28.08
Accuracy 3-4m

When creating such a new location (the same applies to all projections with multiple datum choices), the user must be presented with a related datum dialog to choose from.

comment:3 by lutra, 13 years ago

The same applies, for example, for epsg 3763?

comment:4 by neteler, 13 years ago

Here the test with 3763:

g.proj -c epsg=3763 datumtrans=-1 location=loc_epsg_3763
Location loc_epsg_3763 created!
exit

grass64 ~/grassdata/loc_epsg_3763/PERMANENT/
GRASS 6.4.1svn (loc_epsg_3763):~ > g.proj -w
PROJCS["Transverse Mercator",
    GEOGCS["grs80",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",39.66825833333333],
    PARAMETER["central_meridian",-8.133108333333334],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

Compared to http://spatialreference.org/ref/epsg/3763/prettywkt/ it looks as expected in GRASS.

Concerning the original report:

projection: 99 (Transverse Mercator)
zone:       0
datum:      towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
ellipsoid:  international
...

that is (see http://trac.osgeo.org/grass/browser/grass/branches/releasebranch_6_4/lib/gis/datumtransform.table) the first choice here:

63	rome40  "towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68" "Italy (Peninsular Part)" "Accuracy 3-4m"
64	rome40  "towgs84=-168.6,-34.0,38.6,-0.374,-0.679,-1.379,-9.48" "Italy (Sardinia)" "Accuracy 3-4m"
65	rome40  "towgs84=-50.2,-50.4,84.8,-0.690,-2.012,0.459,-28.08" "Italy (Sicily)" "Accuracy 3-4m"

Something in the QGIS method to create a new location with multiple datum choices seems to be broken.

comment:5 by pcav, 13 years ago

In QGIS we used a different approach: we defined additional ad hoc projections (e.g. code 40000) with the additional datum. It used to work, so I think it has been broken by some changes recently.

Note: See TracTickets for help on using tickets.