[GRASS-SVN] r71523 - grass/trunk/lib/proj
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Oct 1 13:44:08 PDT 2017
Author: mmetz
Date: 2017-10-01 13:44:08 -0700 (Sun, 01 Oct 2017)
New Revision: 71523
Modified:
grass/trunk/lib/proj/convert.c
grass/trunk/lib/proj/get_proj.c
Log:
libproj: use PROJ4 EXTENSION in WKT (TODO: PROJ4_GRIDS)
Modified: grass/trunk/lib/proj/convert.c
===================================================================
--- grass/trunk/lib/proj/convert.c 2017-10-01 13:55:35 UTC (rev 71522)
+++ grass/trunk/lib/proj/convert.c 2017-10-01 20:44:08 UTC (rev 71523)
@@ -388,6 +388,7 @@
const char *pszProjCS = NULL;
char *datum = NULL;
struct gpj_datum dstruct;
+ const char *ograttr;
if (hSRS == NULL)
goto default_to_xy;
@@ -399,6 +400,66 @@
* to start with... */
OSRMorphFromESRI(hSRS);
+ *projinfo = G_create_key_value();
+
+ /* use proj4 definition from EXTENSION attribute if existing */
+ ograttr = OSRGetAttrValue(hSRS, "EXTENSION", 0);
+ if (ograttr && *ograttr && strcmp(ograttr, "PROJ4") == 0) {
+ ograttr = OSRGetAttrValue(hSRS, "EXTENSION", 1);
+ G_debug(3, "proj4 extension:");
+ G_debug(3, "%s", ograttr);
+
+ if (ograttr && *ograttr) {
+ char *proj4ext;
+ OGRSpatialReferenceH hSRS2;
+
+ proj4ext = G_store(ograttr);
+ hSRS2 = OSRNewSpatialReference(NULL);
+ /* test */
+ if (OSRImportFromProj4(hSRS2, proj4ext) != OGRERR_NONE) {
+ G_warning(_("Updating spatial reference with embedded proj4 definition failed. "
+ "Proj4 definition: <%s>"), proj4ext);
+ }
+ else {
+ OSRDestroySpatialReference(hSRS2);
+ /* update OGR spatial reference with proj4 string */
+ /* TODO: replace warning with important_message once confirmed working */
+ G_warning(_("Updating spatial reference with embedded proj4 definition"));
+
+ /* -------------------------------------------------------------------- */
+ /* Derive the user name for the coordinate system. */
+ /* -------------------------------------------------------------------- */
+ pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
+ if (!pszProjCS)
+ pszProjCS = OSRGetAttrValue(hSRS, "GEOGCS", 0);
+
+ if (pszProjCS) {
+ G_set_key_value("name", pszProjCS, *projinfo);
+ }
+ else if (pszProj) {
+ char path[4095];
+ char name[80];
+
+ /* use name of the projection as name for the coordinate system */
+
+ sprintf(path, "%s/etc/proj/projections", G_gisbase());
+ if (G_lookup_key_value_from_file(path, pszProj, name, sizeof(name)) >
+ 0)
+ G_set_key_value("name", name, *projinfo);
+ else
+ G_set_key_value("name", pszProj, *projinfo);
+ }
+
+ OSRDestroySpatialReference(hSRS);
+ hSRS = OSRNewSpatialReference(NULL);
+ if (OSRImportFromProj4(hSRS, proj4ext) != OGRERR_NONE) {
+ G_fatal_error(_("Can't parse PROJ.4-style parameter string"));
+ }
+ }
+ G_free(proj4ext);
+ }
+ }
+
/* -------------------------------------------------------------------- */
/* Set cellhd for well known coordinate systems. */
/* -------------------------------------------------------------------- */
@@ -474,6 +535,11 @@
pszProj = pszValue;
}
+ /* discard @null nadgrids */
+ if (G_strcasecmp(pszToken, "nadgrids") == 0 &&
+ G_strcasecmp(pszValue, "@null") == 0)
+ continue;
+
/* Ellipsoid and datum handled separately below */
if (G_strcasecmp(pszToken, "ellps") == 0
|| G_strcasecmp(pszToken, "a") == 0
@@ -493,33 +559,32 @@
if (!pszProj)
G_warning(_("No projection name! Projection parameters likely to be meaningless."));
- *projinfo = G_create_key_value();
-
/* -------------------------------------------------------------------- */
/* Derive the user name for the coordinate system. */
/* -------------------------------------------------------------------- */
- pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
- if (!pszProjCS)
- pszProjCS = OSRGetAttrValue(hSRS, "GEOGCS", 0);
+ if (!G_find_key_value("name", *projinfo)) {
+ pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
+ if (!pszProjCS)
+ pszProjCS = OSRGetAttrValue(hSRS, "GEOGCS", 0);
- if (pszProjCS) {
- G_set_key_value("name", pszProjCS, *projinfo);
- }
- else if (pszProj) {
- char path[4095];
- char name[80];
-
- /* use name of the projection as name for the coordinate system */
+ if (pszProjCS) {
+ G_set_key_value("name", pszProjCS, *projinfo);
+ }
+ else if (pszProj) {
+ char path[4095];
+ char name[80];
+
+ /* use name of the projection as name for the coordinate system */
- sprintf(path, "%s/etc/proj/projections", G_gisbase());
- if (G_lookup_key_value_from_file(path, pszProj, name, sizeof(name)) >
- 0)
- G_set_key_value("name", name, *projinfo);
- else
- G_set_key_value("name", pszProj, *projinfo);
+ sprintf(path, "%s/etc/proj/projections", G_gisbase());
+ if (G_lookup_key_value_from_file(path, pszProj, name, sizeof(name)) >
+ 0)
+ G_set_key_value("name", name, *projinfo);
+ else
+ G_set_key_value("name", pszProj, *projinfo);
+ }
}
-
/* -------------------------------------------------------------------- */
/* Find the GRASS datum name and choose parameters either */
/* interactively or not. */
@@ -581,18 +646,18 @@
}
if (paramsets > 0) {
- struct gpj_datum_transform_list *list, *old;
+ struct gpj_datum_transform_list *tlist, *old;
- list = GPJ_get_datum_transform_by_name(datum);
+ tlist = GPJ_get_datum_transform_by_name(datum);
- if (list != NULL) {
+ if (tlist != NULL) {
do {
- if (list->count == datumtrans)
- chosenparams = G_store(list->params);
- old = list;
- list = list->next;
+ if (tlist->count == datumtrans)
+ chosenparams = G_store(tlist->params);
+ old = tlist;
+ tlist = tlist->next;
GPJ_free_datum_transform(old);
- } while (list != NULL);
+ } while (tlist != NULL);
}
}
@@ -653,8 +718,10 @@
* accept first one that matches. These numbers were found
* by trial and error and could be fine-tuned, or possibly
* a direct comparison of IEEE floating point values used. */
- if ((a == list->a || fabs(a - list->a) < 0.1 || fabs(1 - a / list->a) < 0.0000001) && ((es == 0 && list->es == 0) || /* Special case for sphere */
- (invflat == list->rf || fabs(invflat - list->rf) < 0.0000001))) {
+ if ((a == list->a || fabs(a - list->a) < 0.1 || fabs(1 - a / list->a) < 0.0000001) &&
+ ((es == 0 && list->es == 0) ||
+ /* Special case for sphere */
+ (invflat == list->rf || fabs(invflat - list->rf) < 0.0000001))) {
ellps = G_store(list->name);
break;
}
Modified: grass/trunk/lib/proj/get_proj.c
===================================================================
--- grass/trunk/lib/proj/get_proj.c 2017-10-01 13:55:35 UTC (rev 71522)
+++ grass/trunk/lib/proj/get_proj.c 2017-10-01 20:44:08 UTC (rev 71523)
@@ -64,7 +64,7 @@
double a, es, rf;
int returnval = 1;
char buffa[300], factbuff[50];
- char proj_in[50], *datum, *params;
+ char proj_in[250], *datum, *params;
projPJ *pj;
proj_in[0] = '\0';
More information about the grass-commit
mailing list