datums not recognized by g.proj?

classic Classic list List threaded Threaded
15 messages Options
Reply | Threaded
Open this post in threaded view
|

datums not recognized by g.proj?

Michael Barton
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 <unknown> 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)












_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Paul Kelly
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
_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Michael Barton
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 <[hidden email]>
 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

_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Paul Kelly
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
_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Michael Barton
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 <[hidden email]>
 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

_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Paul Kelly
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
_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

hamish-2
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)
_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Markus Neteler
In reply to this post by Paul Kelly
(back to list)

On Mon, Oct 1, 2012 at 5:55 PM, Michael Barton <[hidden email]> 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,@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 :)

Markus
_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Michael Barton
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 <[hidden email]>
 wrote:

> (back to list)
>
> On Mon, Oct 1, 2012 at 5:55 PM, Michael Barton <[hidden email]> 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!]

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 :)
>
> Markus

_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Markus Neteler
In reply to this post by Markus Neteler
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
_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Markus Neteler
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Paul Kelly
In reply to this post by Markus Neteler
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
_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Michael Barton
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 <[hidden email]>
 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

_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Michael Barton
In reply to this post by Paul Kelly
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 <[hidden email]> 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

_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev
Reply | Threaded
Open this post in threaded view
|

Re: datums not recognized by g.proj?

Markus Neteler
In reply to this post by Paul Kelly
Hallo Paul,

On Mon, Oct 1, 2012 at 11:44 PM, Paul Kelly
<[hidden email]> 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 :)

Thanks much, Paul!

Markus
_______________________________________________
grass-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/grass-dev