[GRASS-dev] [GRASS GIS] #2456: read CSV from GDAL data directory

GDAL · April 24, 2016, 5:07pm

#2456: read CSV from GDAL data directory
---------------------------------+---------------------------------
  Reporter: martinl | Owner: grass-dev@…
      Type: task | Status: new
  Priority: blocker | Milestone: 7.1.0
Component: Projections/Datums | Version: svn-releasebranch70
Resolution: | Keywords: gdal
       CPU: Unspecified | Platform: Unspecified
---------------------------------+---------------------------------

Comment (by neteler):

Also EPSG 3358 (NC) is an issue:

{{{
grass71 -c epsg:3358 ~/grassdata/nc_data
GRASS 7.1.svn (nc_data):~ > g.proj -w
PROJCS["NAD83(HARN) / North Carolina",
     GEOGCS["grs80",
         DATUM["unknown",
SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101],
             TOWGS84[0,0,0,0,0,0,0]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Lambert_Conformal_Conic_2SP"],
     PARAMETER["standard_parallel_1",36.16666666666666],
     PARAMETER["standard_parallel_2",34.33333333333334],
     PARAMETER["latitude_of_origin",33.75],
     PARAMETER["central_meridian",-79],
     PARAMETER["false_easting",609601.22],
     PARAMETER["false_northing",0],
     UNIT["Meter",1]]
}}}

Expected: http://epsg.io/3358 :

{{{
PROJCS["NAD83(HARN) / North Carolina",
     GEOGCS["NAD83(HARN)",
         DATUM["NAD83_High_Accuracy_Reference_Network",
             SPHEROID["GRS 1980",6378137,298.257222101,
                 AUTHORITY["EPSG","7019"]],
             TOWGS84[0,0,0,0,0,0,0],
             AUTHORITY["EPSG","6152"]],
         PRIMEM["Greenwich",0,
             AUTHORITY["EPSG","8901"]],
         UNIT["degree",0.0174532925199433,
             AUTHORITY["EPSG","9122"]],
         AUTHORITY["EPSG","4152"]],
     PROJECTION["Lambert_Conformal_Conic_2SP"],
     PARAMETER["standard_parallel_1",36.16666666666666],
     PARAMETER["standard_parallel_2",34.33333333333334],
     PARAMETER["latitude_of_origin",33.75],
     PARAMETER["central_meridian",-79],
     PARAMETER["false_easting",609601.22],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]],
     AXIS["X",EAST],
     AXIS["Y",NORTH],
     AUTHORITY["EPSG","3358"]]
}}}

The recognition of datum name "NAD83_High_Accuracy_Reference_Network"
fails.

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/2456#comment:16>
GRASS GIS <https://grass.osgeo.org>

Paul_Kelly · April 25, 2016, 8:02pm

Hi Markus, Even,

On 04/24/2016 10:29 PM, Markus Neteler wrote:

Hi Paul,

you are probably the person with most insights into lib/proj/convert.c.

I try to debug why the North Carolina Location cannot be generated
properly (EPSG 3358) but comes out like this:

  grass71 -c epsg:3358 ~/grassdata/nc_data
  GRASS 7.1.svn (nc_data):~ > g.proj -w
  PROJCS["NAD83(HARN) / North Carolina",
      GEOGCS["grs80",
          DATUM["unknown",

The datum name "unknown" is coming through from GDAL. Here is a simple test program to illustrate it (using latest SVN GDAL from this evening):

#include <stdio.h>
#include <ogr_srs_api.h>

main()
{
OGRSpatialReferenceH hSRS;
const char *datum;

hSRS = OSRNewSpatialReference(NULL);
OSRImportFromProj4(hSRS, "+init=epsg:3358");

datum = OSRGetAttrValue(hSRS, "DATUM", 0);
printf("Datum name is %s\n", datum);

return 0;
}

Compile as:
gcc -Wall -o test test.c -lgdal
and run:
./test
Datum name is unknown

Perhaps there is a better / more correct way to get the co-ordinate system definition from GDAL than by requesting it as a "+init=epsg:XXXX" PROJ string?

Paul

rouault · April 25, 2016, 8:12pm

Le lundi 25 avril 2016 22:02:20, Paul Kelly a écrit :

Hi Markus, Even,

On 04/24/2016 10:29 PM, Markus Neteler wrote:
> Hi Paul,
>
> you are probably the person with most insights into lib/proj/convert.c.
>
> I try to debug why the North Carolina Location cannot be generated
>
> properly (EPSG 3358) but comes out like this:
> grass71 -c epsg:3358 ~/grassdata/nc_data
> GRASS 7.1.svn (nc_data):~ > g.proj -w
> PROJCS["NAD83(HARN) / North Carolina",
>
> GEOGCS["grs80",
>
> DATUM["unknown",

The datum name "unknown" is coming through from GDAL. Here is a simple
test program to illustrate it (using latest SVN GDAL from this evening):

#include <stdio.h>
#include <ogr_srs_api.h>

main()
{
     OGRSpatialReferenceH hSRS;
     const char *datum;

     hSRS = OSRNewSpatialReference(NULL);
     OSRImportFromProj4(hSRS, "+init=epsg:3358");

     datum = OSRGetAttrValue(hSRS, "DATUM", 0);
     printf("Datum name is %s\n", datum);

     return 0;
}

Compile as:
gcc -Wall -o test test.c -lgdal
and run:
./test
Datum name is unknown

Perhaps there is a better / more correct way to get the co-ordinate
system definition from GDAL than by requesting it as a "+init=epsg:XXXX"
PROJ string?

Yes, if you use OSRImportFromProj4(), the building of the SRS WKT will be done
by "reverse engineering" the proj.4 string that gets expended by proj.4 from
+init=epsg:3358 to "+proj=lcc +lat_1=36.16666666666666
+lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs". Which has lost any
datum name.

The solution is simple : use OSRImportFromEPSG(hSRS, code)

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

Paul_Kelly · April 25, 2016, 8:37pm

On 04/25/2016 09:12 PM, Even Rouault wrote:

Le lundi 25 avril 2016 22:02:20, Paul Kelly a écrit :

[...]

Perhaps there is a better / more correct way to get the co-ordinate
system definition from GDAL than by requesting it as a "+init=epsg:XXXX"
PROJ string?

Yes, if you use OSRImportFromProj4(), the building of the SRS WKT will be done
by "reverse engineering" the proj.4 string that gets expended by proj.4 from
+init=epsg:3358 to "+proj=lcc +lat_1=36.16666666666666
+lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs". Which has lost any
datum name.

The solution is simple : use OSRImportFromEPSG(hSRS, code)

Ah, thanks for the hint - but I had got confused: g.proj has actually been using OSRImportFromEPSG() for a long time (since r31652, in 2008 actually). Faulty memory on my part.

The problem is actually even simpler (datum simply not recognised in GRASS) and there are two possible solutions:

1) If we want NAD83_High_Accuracy_Reference_Network to be treated as a completely separate datum from North_American_Datum_1983, then we add it as a new line in lib/gis/datum.table as per the first attached patch.

2) If we want NAD83_High_Accuracy_Reference_Network to be treated as equivalent to North_American_Datum_1983 (this means once the location has been created, 'g.proj -w' will report the datum name as North_American_Datum_1983), then we add two new lines to the equivalent pairs array in lib/proj/convert.c as per the second attached patch.

I don't really feel qualified to decide which is the most desirable behaviour, so I'll leave that up to someone else to decide, if that's OK.

Paul

(attachments)

datum.table.patch (512 Bytes)
convert.c.patch (295 Bytes)

neteler · April 25, 2016, 8:47pm

(cc Helena)

First of all, thanks to Even and Paul to shed some light on this!

On Mon, Apr 25, 2016 at 10:37 PM, Paul Kelly
<paul-grass@stjohnspoint.co.uk> wrote:
...

The problem is actually even simpler (datum simply not recognised in GRASS)
and there are two possible solutions:

1) If we want NAD83_High_Accuracy_Reference_Network to be treated as a
completely separate datum from North_American_Datum_1983, then we add it as
a new line in lib/gis/datum.table as per the first attached patch.

2) If we want NAD83_High_Accuracy_Reference_Network to be treated as
equivalent to North_American_Datum_1983 (this means once the location has
been created, 'g.proj -w' will report the datum name as
North_American_Datum_1983), then we add two new lines to the equivalent
pairs array in lib/proj/convert.c as per the second attached patch.

I don't really feel qualified to decide which is the most desirable
behaviour, so I'll leave that up to someone else to decide, if that's OK.

I also don't feel too qualified here but from a user's point of
view a change of name is confusing.

So option 1) looks better to me (i.e., " datum.table.patch") which is
tested ok here.
Maybe Helena as our NC expert has an opinion here?

Paul, still struggling with "SIRGAS2000"
https://trac.osgeo.org/grass/ticket/2456#comment:15
(g.proj versus testepsg output). Needs a different trick?

Markus

--
Markus Neteler
http://www.mundialis.de - free data with free software
http://courses.neteler.org/blog/

Paul_Kelly · April 25, 2016, 9:15pm

On 04/25/2016 09:47 PM, Markus Neteler wrote:

Paul, still struggling with "SIRGAS2000"
https://trac.osgeo.org/grass/ticket/2456#comment:15
(g.proj versus testepsg output). Needs a different trick?

Hi Markus,
It looks to me like the canonical and equivalent names are the wrong way round. I will attempt a brief summary.

datum.table is the oldest file, going back to the earliest days of datums in GRASS. It should contain:

1. Short code for datum, used only in GRASS
2. Full EPSG name of datum, in quotation marks
3. ellipsoid used with this datum (GRASS ellipsoid code, from ellipse.table)
4. A single 3-parameter datum transformation that should cover the whole area the datum is used in

datumtransform.table is newer and allows multiple different sets of datum transformation parameters to be supplied for any datum. These can be in any format accepted by PROJ.4, so towgs84= 3 parameters or 7 parameters, or nadgrids=
A datum does not have to have an entry in datumtransform.table, and there can be multiple entries for each datum if necessary.

Some datums have different names (i.e. different from the EPSG name) that are used in WKT descriptions. These should be added as a pair of names in the papszDatumEquiv array lib/proj/convert.c. The equivalent name comes first in the pair and the EPSG name (as used in the GRASS datum.table file) comes second.

So for the SIRGAS2000 datum, according to the testepsg program "Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000" is the EPSG name, so this is what should go in datum.table. I don't think you even need the equivalent "Sistema_de_Referencia_Geocentrico_para_America_del_Sur_2000" in convert.c, unless there have been problems with this?

I would suggest the attached patch to solve both problems (taking into account Helena's comments that NAD83 HARN should be separate).

Paul

(attachments)

datum-changes.patch (1.33 KB)

neteler · April 25, 2016, 10:27pm

Hi Paul,

thanks for your patience and the explanations!

n Mon, Apr 25, 2016 at 11:15 PM, Paul Kelly
<paul-grass@stjohnspoint.co.uk> wrote:

I would suggest the attached patch to solve both problems (taking into
account Helena's comments that NAD83 HARN should be separate).

Paul

I have now submitted the changes along with the addition of some code
comments as per your emails in:
https://trac.osgeo.org/grass/changeset/68308

Hope I got it right! If yes, I'll backport that.

Now I have recompiled trunk and made a new test with SIRGAS2000 (seems
to be an interesting test case besides Krovak):

grass71 -c epsg:4674 ~/grassdata/SIRGAS2000

GRASS 7.1.svn (SIRGAS2000):~ > g.proj -w
PROJCS["SIRGAS 2000",
    GEOGCS["grs80",
        DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
            SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]]]

... which still differs from ...

GRASS 7.1.svn (SIRGAS2000):~ > testepsg epsg:4674
Validate Succeeds.
WKT[epsg:4674] =
GEOGCS["SIRGAS 2000",
    DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6674"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4674"]]
...

The main difference is that GRASS actually generates it as PROJCS
while GDAL generates it as GEOGCS. Hence a reverse test of the g.proj
output checked in testepsg fails at time.

thanks
Markus

Paul_Kelly · April 25, 2016, 11:34pm

Hi Markus,

On 04/25/2016 11:27 PM, Markus Neteler wrote:

I have now submitted the changes along with the addition of some code
comments as per your emails in:
https://trac.osgeo.org/grass/changeset/68308

Hope I got it right! If yes, I'll backport that.

Yes, that all looks correct to me.

Now I have recompiled trunk and made a new test with SIRGAS2000 (seems
to be an interesting test case besides Krovak):

grass71 -c epsg:4674 ~/grassdata/SIRGAS2000

GRASS 7.1.svn (SIRGAS2000):~ > g.proj -w
PROJCS["SIRGAS 2000",
     GEOGCS["grs80",
         DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
             SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]]]

... which still differs from ...

GRASS 7.1.svn (SIRGAS2000):~ > testepsg epsg:4674
Validate Succeeds.
WKT[epsg:4674] =
GEOGCS["SIRGAS 2000",
     DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
         SPHEROID["GRS 1980",6378137,298.257222101,
             AUTHORITY["EPSG","7019"]],
         TOWGS84[0,0,0,0,0,0,0],
         AUTHORITY["EPSG","6674"]],
     PRIMEM["Greenwich",0,
         AUTHORITY["EPSG","8901"]],
     UNIT["degree",0.0174532925199433,
         AUTHORITY["EPSG","9122"]],
     AUTHORITY["EPSG","4674"]]
...

The main difference is that GRASS actually generates it as PROJCS
while GDAL generates it as GEOGCS. Hence a reverse test of the g.proj
output checked in testepsg fails at time.

The code handling this was changed pretty recently: https://trac.osgeo.org/grass/changeset/68131/

Simply deleting three lines from that revision fixes it for me (see below) - the PROJCS name seems to be being set unconditionally even when it's not a projected co-ordinate system. Deleting this doesn't seem to cause any regression as the PROJCS name is already set later in the GPJ_grass_to_osr() function.

I haven't committed it as it should probably be tested first with a few different types of co-ordinate system, and I don't have time for that at the minute unfortunately.

Paul

--- lib/proj/convert.c (revision 68310)
+++ lib/proj/convert.c (working copy)
@@ -141,9 +141,6 @@
    return NULL;
      }
      G_free(proj4mod);
- sysname = G_find_key_value("name", proj_info);
- if (sysname)
- OSRSetProjCS(hSRS, sysname);

if ((errcode = OSRExportToWkt(hSRS, &wkt)) != OGRERR_NONE) {
G_warning(_("OGR can't get WKT-style parameter string "