[GRASS-dev] datums not recognized by g.proj?
Michael Barton
Michael.Barton at asu.edu
Mon Oct 1 11:54:15 PDT 2012
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 at osgeo.org>
wrote:
> (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!]
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, 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.
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
More information about the grass-dev
mailing list