towgs84 datum transform options not maintained with shapefile export

[voncasec]

Perhaps this issue is related to the GDAL ticket 4345 (​http://trac.osgeo.org/gdal/ticket/4345).

Including it here in case GRASS handles the data differently.

In general, if I am using a projection defined with towgs84 datum transform parameters and I want to export the layer to a shapefile, the +towgs84 datum transform parameters are not included.

Often times, I will just manually modify the .prj file after the fact to include the TOWGS84 flag to ensure that there is no shift when viewing the locations in a GIS that allow for on the fly reprojections.

Example of g.proj -e

GRASS 6.4.3svn (Bayaguana_NAD27_Caribbean):/usr/local/grass-6.4.3svn/etc > g.proj -e PROJCS["UTM_Zone_19_Northern_Hemisphere",

GEOGCS["GCS_clark66",

DATUM["D_North_American_1927",

SPHEROID["Clarke_1866",6378206.4,294.9786982]],

PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]],

PROJECTIONTransverse_Mercator, PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-69], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["Meter",1]]

Desired output of g.proj -e + .prj file from v.out.ogr

GRASS 6.4.3svn (Bayaguana_NAD27_Caribbean):/usr/local/grass-6.4.3svn/etc > g.proj -e PROJCS["UTM_Zone_19_Northern_Hemisphere",

GEOGCS["GCS_clark66",

DATUM["D_North_American_1927",

SPHEROID["Clarke_1866",6378206.4,294.9786982], TOWGS84[-3.0,142.0,183.0,0,0,0,0]],

PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]],

PROJECTIONTransverse_Mercator, PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-69], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["Meter",1]]

[hamish]

Yes, it's a known issue, and it needs to be solved in GDAL/OGR first. If I recall correctly, WKT doesn't have a defined way to include datum transforms, ISTR some time ago (older than that ticket) Frank suggested it could be packed into something like a comment term, but then that would only be known to other osgeo-stack software which knew to look for it.

Also, the proj4 epsg file includes explicit +towgs84 terms for some SRSs, but others not (several options for the user to choose from under the umbrella of the generalized code); the one referred to in the gdal ticket (epsg:4632) is one of the ones with explicit +towgs84 terms in the epsg def'n.

?, Hamish

[hamish]

Casey wrote: (grass-dev ML)

Thanks Hamish,

I assumed it needed to be looked at with GDAL/OGR first when exporting to a shapefile, I just wasn't sure if this was also the case when using g.proj -e.

Hi,

I did a little digging. You get the TOWGS84 in the WKT if you use 'g.proj -w' and your PERMANENT/PROJ_INFO file has a towgs84 line, but not 'g.proj -we' (-e modifies -w, and triggers it because it implies it). Grid datum transform files are not handled by any variant of Well Known Text afaik (a slight correction to my last comment), except perhaps by using a custom osgeo-family "EXTENSION" keyword, which I'm not sure has been commonly defined yet.

So the good news is that for generic Shapefile export you could leave off the v.out.ogr/g.prok '-e' flag to use generic WKT, which supports TOWGS84. For export to ESRI products it would be lost though, as ESRI apparently (according to the OGR code) doesn't support TOWGS84 statements in the WKT.

specifically 'g.proj -we' calls GPJ_grass_to_wkt() which if the ESRI-compatibility flag is given calls OGR's OSRMorphToESRI() (aka morphToESRI()).

morphToESRI() has this comment as it strips away any TOWGS84:

/* -------------------------------------------------------------------- */
/*      Strip all CT parameters (AXIS, AUTHORITY, TOWGS84, etc).        */
/* -------------------------------------------------------------------- */

OGR's StripCTParms() is called in to do that:

 * This method will remove all components of the coordinate system
 * that are specific to the OGC CT Specification.  That is it will attempt
 * to strip it down to being compatible with the Simple Features 1.0 
 * specification.
...
    poCurrent->StripNodes( "AUTHORITY" );
    poCurrent->StripNodes( "TOWGS84" );
    poCurrent->StripNodes( "AXIS" );
    poCurrent->StripNodes( "EXTENSION" );

    return OGRERR_NONE;
}

see grass's lib/proj/convert.c and GDAL/OGR's ogr/ogr_srs_esri.cpp + ogrspatialreference.cpp.

so in the case of 'g.proj -w -e' it behaves as expected, TOWGS84 is removed on purpose.

hope it helps, Hamish

[neteler]

Please try r68136 in the 7.0.x release branch.

[martinl]

No activity, closing, feel free to reopen if needed.