[GRASS-dev] datums not recognized by g.proj?

[Markus Neteler]

(back to list)

On Mon, Oct 1, 2012 at 5:55 PM, Michael Barton <Michael.Barton at asu.edu> wrote:
...
> The worse problem is that locations are created ignoring the selected datum.

Thanks to the indication below I was able to understand what you mean.

So
1. Run GRASS in an existing location, set
    g.gisenv set=DEBUG=3
   leave GRASS.
2. Restart GRASS, enter Location Wizard, select "Select coordinate system"
3. Enter "utm"
4. Select "[x] Datum with associated ellipsoid"
    Select Zone 29
    Southern Hemisphere No
5. Datum code: eur50   [it is not ed50!]

Now find in the terminal:
D1/3: grass.script.core.start_command(): g.proj proj4=+proj=utm
+datum=eur50 datumtrans=-1
D1/3: grass.script.core.start_command(): g.proj -jf proj4=+proj=utm
+zone=29 +ellps=international  +a=6378388.000  +rf=297.0 +datum=eur50
+dx=-87.0 +dy=-98 +dz=-121 +no_defs location=newLocation2

So, the datum= string is well there. I am using proj-4.7.0-5.fc17.x86_64.

> I assume that g.proj will simply use a default datum (WGS84??) during location creation.

No, see above.
But, yes, the resulting location does not contain the datum string:


GRASS 6.4.3svn (newLocation2):~ > g.proj -w
...
PROJCS["UTM Zone 29, Northern Hemisphere",
    GEOGCS["international",
        DATUM["unknown",
            SPHEROID["International_1924",6378388,297]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

GRASS 6.4.3svn (newLocation2):~ > g.proj -p
...
-PROJ_INFO-------------------------------------------------
name       : Universal Transverse Mercator
proj       : utm
ellps      : international
zone       : 29
no_defs    : defined
-PROJ_UNITS------------------------------------------------
unit       : Meter
units      : Meters
meters     : 1

This is a problem. But I don't know if this ever worked (I never created a
location this way). Michael, please indicate the last GRASS version which
you found it working like this.

> However, running g.proj -d afterwards returns a "Datum parameters not present" error. So I
> can't tell what it is doing. This suggests to me that datum transform problem aside, many
> (most???) locations created by choosing projection and datum are wrong.

No, ONLY, if created in this style while being OK when using the EPSG code.
It is important to distinguish this.

...
> For example, if you have some data from a map that says it is in "Universal
> Transverse Mercator, Zone 29, European Datum 1950", searching on almost
> any of those terms in the projection and datum lists will give you a hit on the
> correct information. But doing the same search in the EPSG window gives you
> nothing, except for searching on "mercator" or "transverse"--and these give
> the wrong information. The best way to find the correct projection via the
> EPSG listing is if you already know that European Datum 1950 has a
> code name of "ED50".

(the code name is "eur50")

As Paul pointed out, PROJ4 doesn't seem to know much (anything?) about EUR50:

cs2cs -ld
__datum_id__ __ellipse___ __definition/comments______________________________
       WGS84 WGS84        towgs84=0,0,0
      GGRS87 GRS80        towgs84=-199.87,74.79,246.62
                          Greek_Geodetic_Reference_System_1987
       NAD83 GRS80        towgs84=0,0,0
                          North_American_Datum_1983
       NAD27 clrk66       nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat
                          North_American_Datum_1927
     potsdam bessel       towgs84=606.0,23.0,413.0
                          Potsdam Rauenberg 1950 DHDN
    carthage clark80      towgs84=-263.0,6.0,431.0
                          Carthage 1934 Tunisia
hermannskogel bessel       towgs84=653.0,-212.0,449.0
                          Hermannskogel
       ire65 mod_airy
towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15
                          Ireland 1965
      nzgd49 intl         towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993
                          New Zealand Geodetic Datum 1949
      OSGB36 airy
towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894
                          Airy 1830

When creating locations using the EPSG way (my preferred method), it
uses GDAL/OGR instead of PROJ4 which knows about EUR50. Hence the
difference.

I do agree that getting the datum lost is bad. Perhaps we could enhance
g.proj to write it to PROJ_INFO when using the "Select coordinate system" way
of creating locations. Proof of concept:

GRASS 6.4.3svn (newLocation2):~ > eval `g.gisenv`

GRASS 6.4.3svn (newLocation2):~ > echo "datum: eur50" >>
$GISDBASE/$LOCATION_NAME/$MAPSET/PROJ_INFO

GRASS 6.4.3svn (newLocation2):~ > g.proj -w
PROJCS["UTM Zone 29, Northern Hemisphere",
    GEOGCS["international",
        DATUM["European_Datum_1950",
            SPHEROID["International_1924",6378388,297]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

GRASS 6.4.3svn (newLocation2):~ > g.proj -p
-PROJ_INFO-------------------------------------------------
name       : Universal Transverse Mercator
proj       : utm
ellps      : international
zone       : 29
no_defs    : defined
datum      : eur50
-PROJ_UNITS------------------------------------------------
unit       : Meter
units      : Meters
meters     : 1

... solved :)

Markus