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

GRASS GIS GRASS GIS developer
1

Markus,

It looks like some datums are not being recognized by g.proj. When this happens, datum transform information is not returned. Take a look at this.

wgs84 is recognized by g.proj

GRASS 7.0.svn (newLocation):~ > g.proj -d proj4=‘+proj=utm +datum=wgs84’
GRASS datum code: wgs84
WKT Name: WGS_1984
Datum transformation parameters (PROJ.4 format):
towgs84=0,0,0,0,0,0,0

but European Datum 1995 (eur50) is not recognized. This datum IS listed in …/etc/proj/datums.table

GRASS 7.0.svn (newLocation):~ > g.proj -d proj4=‘+proj=utm +datum=eur50’
WARNING: Datum not recognised by GRASS and no parameters found
Datum name not present
Datum parameters not present

Unless I misunderstand something, this seems to be a serious bug because it affects location creation. I will file a bug report on this if you can verify.

Thanks
Michael

C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

voice: 480-965-6262 (SHESC), 480-727-9746 (CSDC)
fax: 480-965-7671 (SHESC), 480-727-0709 (CSDC)

www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu

2

Hello all,
Markus has prompted me to try and give some insight into what is going on here, which I will gladly try to do:

Michael Barton wrote:

Markus,

It looks like some datums are not being recognized by g.proj. When this
happens, datum transform information is not returned. Take a look at this.

wgs84 is recognized by g.proj

GRASS 7.0.svn (newLocation):~ > g.proj -d proj4='+proj=utm +datum=wgs84'
GRASS datum code: wgs84
WKT Name: WGS_1984
Datum transformation parameters (PROJ.4 format):
towgs84=0,0,0,0,0,0,0

but European Datum 1995 (eur50) is not recognized. This datum IS listed
in ../etc/proj/datums.table

As far as I can think, the problem here is that +datum=eur50 is not a PROJ.4-compatible datum name. The only hard-coded datum names supported in PROJ.4 strings are those output by the "proj -ld" command, i.e. WGS84, GGRS87, NAD83, NAD27, potsdam, carthage, hermannskogel, ire65, nzgd49 and OSB36.

The string used for the proj4= option in g.proj needs to be compatible with any installation of PROJ.4, so it can't include GRASS-specific datum names such as eur50 and the others in the GRASS datum.table file. Internally to g.proj, the string from the proj4= option is actually parsed by GDAL (using its OSRImportFromProj4() function) - and GDAL is not aware of the internal GRASS datum names so this is the result:

GRASS 7.0.svn (newLocation):~ > g.proj -d proj4='+proj=utm +datum=eur50'
WARNING: Datum <unknown> not recognised by GRASS and no parameters found
Datum name not present
Datum parameters not present

I am slightly confused though, as to why +datum=wgs84 works fine - because the GDAL name is WGS84, i.e. in capitals. The GRASS version is wgs84 (i.e. lower case), so this shouldn't work either - but maybe GDAL is doing a case-insensitive comparison somewhere.

I am very out-of-touch with GRASS development and know very little about the new location wizard, but I wonder if it handles PROJ-4 strings differently and parses them itself, allowing GRASS-specific datum names for the +datum option, and if this is leading to some confusion?

For compatibility with other applications, I think it's best not to use GRASS-specific datum names in PROJ.4 strings as they won't be valid if the PROJ.4 string is exported for use with another GIS application.

Hope this sheds some light on things.

Paul

3

Thanks for the information Paul,

This is directly a g.proj problem rather than a proj4 problem per se--although I don't know how g.proj uses proj4. So it may be indirectly something to do with proj4 too.

I do know that this USED to work fine. That is, eur50 (and all other datums recognized by GRASSS) was recognized by g.proj and spawned a datum transforms list. The datums list and datum transforms list files that g.proj uses (or used to use) are the ones maintained by the GRASS project, not the official ones of proj4.

So there has been some kind of change in the way g.proj, proj4, and the various GRASS projection files interact.

As I said, this is causing incorrect location projections and is thus quite serious.

Michael
____________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

voice: 480-965-6262 (SHESC), 480-727-9746 (CSDC)
fax: 480-965-7671 (SHESC), 480-727-0709 (CSDC)
www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu

On Sep 30, 2012, at 6:19 AM, Paul Kelly <paul-grass@stjohnspoint.co.uk>
wrote:

Hello all,
Markus has prompted me to try and give some insight into what is going on here, which I will gladly try to do:

Michael Barton wrote:

Markus,

It looks like some datums are not being recognized by g.proj. When this
happens, datum transform information is not returned. Take a look at this.

wgs84 is recognized by g.proj

GRASS 7.0.svn (newLocation):~ > g.proj -d proj4='+proj=utm +datum=wgs84'
GRASS datum code: wgs84
WKT Name: WGS_1984
Datum transformation parameters (PROJ.4 format):
towgs84=0,0,0,0,0,0,0

but European Datum 1995 (eur50) is not recognized. This datum IS listed
in ../etc/proj/datums.table

As far as I can think, the problem here is that +datum=eur50 is not a PROJ.4-compatible datum name. The only hard-coded datum names supported in PROJ.4 strings are those output by the "proj -ld" command, i.e. WGS84, GGRS87, NAD83, NAD27, potsdam, carthage, hermannskogel, ire65, nzgd49 and OSB36.

The string used for the proj4= option in g.proj needs to be compatible with any installation of PROJ.4, so it can't include GRASS-specific datum names such as eur50 and the others in the GRASS datum.table file. Internally to g.proj, the string from the proj4= option is actually parsed by GDAL (using its OSRImportFromProj4() function) - and GDAL is not aware of the internal GRASS datum names so this is the result:

GRASS 7.0.svn (newLocation):~ > g.proj -d proj4='+proj=utm +datum=eur50'
WARNING: Datum <unknown> not recognised by GRASS and no parameters found
Datum name not present
Datum parameters not present

I am slightly confused though, as to why +datum=wgs84 works fine - because the GDAL name is WGS84, i.e. in capitals. The GRASS version is wgs84 (i.e. lower case), so this shouldn't work either - but maybe GDAL is doing a case-insensitive comparison somewhere.

I am very out-of-touch with GRASS development and know very little about the new location wizard, but I wonder if it handles PROJ-4 strings differently and parses them itself, allowing GRASS-specific datum names for the +datum option, and if this is leading to some confusion?

For compatibility with other applications, I think it's best not to use GRASS-specific datum names in PROJ.4 strings as they won't be valid if the PROJ.4 string is exported for use with another GIS application.

Hope this sheds some light on things.

Paul

4

Michael Barton wrote:

Thanks for the information Paul,

This is directly a g.proj problem rather than a proj4 problem per se--although I don't know how g.proj uses proj4. So it may be indirectly something to do with proj4 too.

g.proj does not use PROJ.4 for importing projection definitions at all. All the PROJ.4 and WKT string parsing is done by calls to GDAL functions.

I do know that this USED to work fine. That is, eur50 (and all other datums recognized by GRASSS) was recognized by g.proj and spawned a datum transforms list. The datums list and datum transforms list files that g.proj uses (or used to use) are the ones maintained by the GRASS project, not the official ones of proj4.

Hmmm, I really can't see any way it can have ever worked. It certainly isn't how it was supposed to work when it was originally written, although I'm not sure if anybody has changed anything in the meantime.

The way g.proj knows which datum transformations to include in the list that it spawns is by checking the official EPSG name of the datum. This is included as the second field in the GRASS datum.table, and is matched up against the EPSG name passed back from GDAL if a co-ordinate system is specified using a WKT description, georeferenced file or an EPSG code.

But when specifying a co-ordinate system using a GRASS-style datum name in a PROJ.4 string the official EPSG name for the datum is not available, since GDAL (which g.proj uses to parse the PROJ.4 string) knows nothing of the GRASS datum names.

> As I said, this is causing incorrect location projections and is thus quite serious.

Where are the incorrect PROJ.4 strings with GRASS datum names in them being generated? I think that is worth sorting out, since as I said it's not a good idea to be circulating PROJ.4 strings with non-standard datum names that won't work with other GIS.

So there has been some kind of change in the way g.proj, proj4, and the various GRASS projection files interact.

If you can come up with a combination of GRASS and GDAL versions that result in this behaviour (i.e. g.proj accepts GRASS datum names for the +datum option in a PROJ.4 string) then I promise to look into it and debug.

Paul

5

Paul,

See below.

Michael
____________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

voice: 480-965-6262 (SHESC), 480-727-9746 (CSDC)
fax: 480-965-7671 (SHESC), 480-727-0709 (CSDC)
www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu

On Sep 30, 2012, at 11:29 AM, Paul Kelly <paul-grass@stjohnspoint.co.uk>
wrote:

Michael Barton wrote:

Thanks for the information Paul,

This is directly a g.proj problem rather than a proj4 problem per se--although I don't know how g.proj uses proj4. So it may be indirectly something to do with proj4 too.

g.proj does not use PROJ.4 for importing projection definitions at all. All the PROJ.4 and WKT string parsing is done by calls to GDAL functions.

I do know that this USED to work fine. That is, eur50 (and all other datums recognized by GRASSS) was recognized by g.proj and spawned a datum transforms list. The datums list and datum transforms list files that g.proj uses (or used to use) are the ones maintained by the GRASS project, not the official ones of proj4.

Hmmm, I really can't see any way it can have ever worked. It certainly isn't how it was supposed to work when it was originally written, although I'm not sure if anybody has changed anything in the meantime.

The way g.proj knows which datum transformations to include in the list that it spawns is by checking the official EPSG name of the datum. This is included as the second field in the GRASS datum.table, and is matched up against the EPSG name passed back from GDAL if a co-ordinate system is specified using a WKT description, georeferenced file or an EPSG code.

But the only way to *find* the second field in datum.table is by searching on the first field. So either there is something wrong with the search (e.g., it won't search on "eur50" for some reason) or the 2nd field is wrong (e.g., "European_Datum_1950").

Moreover, the corresponding entries in datum_transform.table match the first field in datum.table. That is, there are entries in datum_transform.table with a first field of "eur50", the same as the entry in datum.table. So again, I don't see how the correct datum transforms can even be found if "eur50" is not recognized by g.proj.

But when specifying a co-ordinate system using a GRASS-style datum name in a PROJ.4 string the official EPSG name for the datum is not available, since GDAL (which g.proj uses to parse the PROJ.4 string) knows nothing of the GRASS datum names.

So what is the purpose of datum.table and datum_transform.table?

> As I said, this is causing incorrect location projections and is thus quite serious.

Where are the incorrect PROJ.4 strings with GRASS datum names in them being generated? I think that is worth sorting out, since as I said it's not a good idea to be circulating PROJ.4 strings with non-standard datum names that won't work with other GIS.

The only way to generate a proj4 string interactively is to pick the appropriate values from a table. When g.proj was still an interactive terminal module, running it would bring up the lists from these tables to generate a proper projection. If all the data tables in ../etc/proj are wrong and are not used for anything anymore, we need to get rid of them and replace them with something that works. A nice, clean set of 4 csv files that includes all relevant data (projections, datums, datum transforms, ellipsoids) would be much easier to parse. The stuff currently in ../etc/proj is currently a real mishmash and difficult to parse.

So there has been some kind of change in the way g.proj, proj4, and the various GRASS projection files interact.

If you can come up with a combination of GRASS and GDAL versions that result in this behaviour (i.e. g.proj accepts GRASS datum names for the +datum option in a PROJ.4 string) then I promise to look into it and debug.

The problem is that it no longer works this way, apparently making it impossible to generate a correct projection based on the information that comes with GRASS.

Michael

Paul

6

Hi Michael,

On Sun, 30 Sep 2012, Michael Barton wrote:
[...]

The only way to generate a proj4 string interactively is to pick the appropriate values from a table. When g.proj was still an interactive terminal module, running it would bring up the lists from these tables to generate a proper projection.

Are you perhaps thinking of g.setproj? It does something fairly similar to what you describe, whereas g.proj has never had a very extensive interactive mode. That might explain a lot of the confusion. g.setproj is a completely different beast, one that I was never really able to tame.

If all the data tables in ../etc/proj are wrong and are not used for anything anymore, we need to get rid of them and replace them with something that works. A nice, clean set of 4 csv files that includes all relevant data (projections, datums, datum transforms, ellipsoids) would be much easier to parse. The stuff currently in ../etc/proj is currently a real mishmash and difficult to parse.

That sounds like a good idea, although it would be a massive amount of work.

If you can come up with a combination of GRASS and GDAL versions that result in this behaviour (i.e. g.proj accepts GRASS datum names for the +datum option in a PROJ.4 string) then I promise to look into it and debug.

The problem is that it no longer works this way

What I was saying was that if you can let me know a certain combination of GRASS and GDAL versions in which g.proj *does* accept GRASS-specific datum names for the +datum option in a PROJ.4 string, let me know which versions and I will investigate. As far as I'm aware up to now it never did work like that.

Best regards

Paul

7

Paul wrote:

Are you perhaps thinking of g.setproj? It does something
fairly similar to what you describe, whereas g.proj has
never had a very extensive interactive mode. That might
explain a lot of the confusion. g.setproj is a completely
different beast, one that I was never really able to tame.

I dunno, AFAIK g.setproj in G6 works extremely well. Extraordinarily well
even. Despite its complications if in doubt it's the guaranteed to work
fallback method.

Michael wrote:

> If all the data tables in ../etc/proj are wrong and are
> not used for anything anymore, we need to get rid of them
> and replace them with something that works.

they are not wrong; they are still used; and they should stay (for now).
Certainly the old method is more robust than the new one is. When the
new one gets to a better state than the old, then let's talk about
replacement. But we're not there yet.. not even close IMO.

Note in GRASS 7 there is a bug in the location wizard. I need the help of
a wxPy expert to fix it..
  http://trac.osgeo.org/grass/ticket/1513#comment:3

this is known to break the ellipsoid (maybe datum too?) setting, so a good
thing to fix before digging too deep.

A nice, clean set of 4 csv files that includes all relevant data

it ain't broke. please don't try to fix it.

If you are not being prompted for datum transform opts, my first question
is --> are you running PROJ4 version 4.7.0 <-- ? (not a svn dev checkout)
If not, rebuild against that version and try again.
(see original report in ongoing http://trac.osgeo.org/grass/ticket/1452)

(projections, datums, datum transforms, ellipsoids) would be
much easier to parse. The stuff currently in ../etc/proj is
currently a real mishmash and difficult to parse.

it's all cleanly formated, just that the logic of how it all works
is perhaps a bit convoluted. I'm pretty sure that between Paul and me
we know how it all fits together though so could add any missing
documentation.

see `proj -ld` for a list of datums known to PROJ.4. If trying to convert
non-epsg coded manually specified strings to +proj4 format, then back
to grass is failing, it seems to me the obvious solution is to only use
the +proj4 parsing for epsg codes, and use the decades worth of field
tested grass tables directly. (~recreate what g.setproj does in python,
or port a non-interactive g.setproj to a new lib fn for 'g.proj -c' to
use, if needed [I'm not sure it is; just don't try to pass GRASS datum
names to +proj=, as PROJ4 only knows its internal offerings])

thanks,
Hamish

(sorry if this is a bit off-base, I haven't seen the full thread yet; I'll
be playing catchup when I can, but I'll still be mostly offline for a wee
while yet/ for anything really important try my personal email address
which I'll check in on more often)

8

(back to list)

On Mon, Oct 1, 2012 at 5:55 PM, Michael Barton <Michael.Barton@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 " 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,@alaska,@ntv2_0.gsb,@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 :slight_smile:

Markus

9

This is encouraging. See below.

Michael
______________________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University
Tempe, AZ 85287-2402
USA

voice: 480-965-6262 (SHESC), 480-727-9746 (CSDC)
fax: 480-965-7671(SHESC), 480-727-0709 (CSDC)
www: http://csdc.asu.edu, http://shesc.asu.edu
    http://www.public.asu.edu/~cmbarton

On Oct 1, 2012, at 11:13 AM, Markus Neteler <neteler@osgeo.org>
wrote:

(back to list)

On Mon, Oct 1, 2012 at 5:55 PM, Michael Barton <Michael.Barton@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 " Datum with associated ellipsoid"
   Select Zone 29
   Southern Hemisphere No
5. Datum code: eur50 [it is not ed50!]

Right. It uses the datum code that comes from GRASS datum.table

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.

Yes. But the datum=eur50 argument seems to be ignored. No error raised. But maybe this is not a problem for accuracy if the other values are set as well?

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:

If the problem is just no datum string but the projection is created appropriately for the given datum, that's a big relief. There are still 2 issues to fix 1) the GRASS datum string has to be put into PROJ.INFO and 2) we still are missing the correct datum transformations. The latter can be fixed in the wxPython code, but I'm not sure how to fix #1.

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.

I cannot say when this worked. What I remember from 5 years back or so is that datum transforms were spawned but now they are not. But perhaps this was because I was luck in testing with a GRASS datum code that worked (i.e., sufficiently similar to a proj4 code) and incorrect datum codes do not raise an error to indicate when it is not working.

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.

This is correct. It is NOT ALL locations that are potentially a problem but locations created by picking projections and datums AND for GRASS datum codes that are not recognized (probably many, but clearly not all of them). With some effort of running g.proj -d with each of the GRASS datum codes, we could ID which are a problem. I can try to do that over the next few days. If someone else can do it quicker or in an automated way, that's even better.

...

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,@alaska,@ntv2_0.gsb,@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.

Actually, EPSG does not know about GRASS datum codes either. EUR50 (or eur50) is not listed anywhere in the EPSG projections list. You must know that it is ED50. I don't know if GDAL/OGR does or does not recognize eur50.

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:

This would certainly solve a big chunk of the problem. Then we need to fix the datum transform issue.

Michael

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 :slight_smile:

Markus

10

For the record, see also this post:

http://fwarmerdam.blogspot.it/2010/03/in-last-few-weeks-i-believe-i-have-made.html
"Slaying the Datum Shift Dragon"

Markus

11

Even more, this has been discussed in 2010:

(this are the most relevant posts)
http://lists.osgeo.org/pipermail/grass-dev/2010-March/049633.html
http://lists.osgeo.org/pipermail/grass-dev/2010-March/049638.html
http://lists.osgeo.org/pipermail/grass-dev/2010-March/049674.html
http://lists.osgeo.org/pipermail/grass-dev/2010-March/049705.html
http://lists.osgeo.org/pipermail/grass-dev/2010-March/049709.html

Markus

12

Markus Neteler wrote:

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

I have just committed r53297 to trunk, which adds a new datum= option to g.proj, which does something very similar to this. If a GRASS datum code is given for the datum= option, it will override any datum in the input co-ordinate system (or add one if it is missing).

So you can now do something like:
g.proj -c loc=spain proj4="+proj=utm +zone=30 +ellps=intl" \
datum=eur50 datumtrans=-1
which will correctly prompt for all the datum transformation options for eur50. You can also do
g.proj datum=list
to get a list of all supported datums like g.setproj does, but in a more easily parseable format similar to the output from datumtrans=-1.

Hope that helps a bit.

Paul

13

GREAT!!!!

I'm in class now (on a break) but will compile and test ASAP. If it works, I'll make required changes in the location wizard.

Michael
____________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

voice: 480-965-6262 (SHESC), 480-727-9746 (CSDC)
fax: 480-965-7671 (SHESC), 480-727-0709 (CSDC)
www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu

On Oct 1, 2012, at 2:44 PM, Paul Kelly <paul-grass@stjohnspoint.co.uk>
wrote:

Markus Neteler wrote:

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

I have just committed r53297 to trunk, which adds a new datum= option to g.proj, which does something very similar to this. If a GRASS datum code is given for the datum= option, it will override any datum in the input co-ordinate system (or add one if it is missing).

So you can now do something like:
g.proj -c loc=spain proj4="+proj=utm +zone=30 +ellps=intl" \
datum=eur50 datumtrans=-1
which will correctly prompt for all the datum transformation options for eur50. You can also do
g.proj datum=list
to get a list of all supported datums like g.setproj does, but in a more easily parseable format similar to the output from datumtrans=-1.

Hope that helps a bit.

Paul

14

Just tested on GRASS 7 and it looks like it works very well.

PROJ.INFO file created using projection and datum selection in location wizard BEFORE updating

name: Universal Transverse Mercator
proj: utm
ellps: international
zone: 29
no_defs: defined

PROJ.INFO file recreated by g.proj AFTER updating

name: Universal Transverse Mercator
proj: utm
ellps: international
zone: 29
no_defs: defined
datum: eur50
towgs84: -131,-100.3,-163.4,-1.244,-0.020,-1.144,9.39

Using the new g.proj with datumtrans=-1 brings up the list of datum transforms for "eur50" as it should.

I will update the location wizard code in trunk shortly to take advantage of this. If others can then test, we should backport to 6.4.3.

Thanks so much. Serious problem identified over the weekend and fixed by Monday night. I can tell my class tomorrow. Wow.

Michael

____________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

voice: 480-965-6262 (SHESC), 480-727-9746 (CSDC)
fax: 480-965-7671 (SHESC), 480-727-0709 (CSDC)
www: http://www.public.asu.edu/~cmbarton, http://csdc.asu.edu

On Oct 1, 2012, at 2:44 PM, Paul Kelly <paul-grass@stjohnspoint.co.uk> wrote:

Markus Neteler wrote:

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

I have just committed r53297 to trunk, which adds a new datum= option to g.proj, which does something very similar to this. If a GRASS datum code is given for the datum= option, it will override any datum in the input co-ordinate system (or add one if it is missing).

So you can now do something like:
g.proj -c loc=spain proj4="+proj=utm +zone=30 +ellps=intl" \
datum=eur50 datumtrans=-1
which will correctly prompt for all the datum transformation options for eur50. You can also do
g.proj datum=list
to get a list of all supported datums like g.setproj does, but in a more easily parseable format similar to the output from datumtrans=-1.

Hope that helps a bit.

Paul

15

Hallo Paul,

On Mon, Oct 1, 2012 at 11:44 PM, Paul Kelly
<paul-grass@stjohnspoint.co.uk> wrote:
...

I have just committed r53297 to trunk, which adds a new datum= option to
g.proj, which does something very similar to this. If a GRASS datum code is
given for the datum= option, it will override any datum in the input
co-ordinate system (or add one if it is missing).

...

You can also do
g.proj datum=list
to get a list of all supported datums like g.setproj does, but in a more
easily parseable format similar to the output from datumtrans=-1.

Hope that helps a bit.

I am very grateful to you for this improvement, especially since you had
to come back to us :slight_smile:

Thanks much, Paul!

Markus