Add new datum= option allowing the datum in the input co-ordinate

system to be overridden (or added if missing). Should be used very
sparingly, as it doesn't check that the correct ellipsoid is used
which could potentially result in ellipsoid/datum mis-matches.


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@53297 15284696-431f-4ddb-bdfa-cd5b030d7da7

general/g.proj/datumtrans.c

@@ -18,6 +18,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include <grass/gis.h>
+#include <grass/glocale.h>
 #include <grass/gprojects.h>
 
 #include "local_proto.h"
@@ -27,6 +28,9 @@
  * \brief Add or replace datum transformation parameters to the current
  *        co-ordinate system definition
  * 
+ * \param datum       Use this datum (overrides any datum found in the
+ *                    current co-ordinate system definition).
+ * 
  * \param datumtrans  Index number of parameter set to use, 0 to leave
  *                    unspecifed (or remove specific parameters, leaving just
  *                    the datum name), -1 to list the available parameter
@@ -38,17 +42,24 @@
  * \return            1 if a change was made, 0 if not.
  **/
 
-int set_datumtrans(int datumtrans, int force)
+int set_datumtrans(char *datum, int datumtrans, int force)
 {
-    char *params, *datum = NULL;
     int paramsets, status;
 
     if (cellhd.proj == PROJECTION_XY)
 	return 0;
 
-    status = GPJ__get_datum_params(projinfo, &datum, &params);
-    G_debug(3, "set_datumtrans(): GPJ__get_datum_params() status=%d", status);
-    G_free(params);
+    if (datum)
+	/* datum name provided; set status to 1 to indicate no parameters available */
+	status = 1;
+    else {
+	/* extract datum name and parameter status from input co-ordinate system */
+	char *params;
+
+	status = GPJ__get_datum_params(projinfo, &datum, &params);
+	G_debug(3, "set_datumtrans(): GPJ__get_datum_params() status=%d", status);
+	G_free(params);
+    }
 
     if (datum) {
 	/* A datum name is specified; need to determine if
@@ -136,7 +147,8 @@ int set_datumtrans(int datumtrans, int force)
 	/* Copy old PROJ_INFO, skipping out any keys related
 	 * to datum parameters */
 	for (i = 0; i < projinfo->nitems; i++) {
-	    if (strcmp(projinfo->key[i], "dx") == 0
+	    if (strcmp(projinfo->key[i], "datum") == 0
+		|| strcmp(projinfo->key[i], "dx") == 0
 		|| strcmp(projinfo->key[i], "dy") == 0
 		|| strcmp(projinfo->key[i], "dz") == 0
 		|| strcmp(projinfo->key[i], "datumparams") == 0
@@ -148,7 +160,9 @@ int set_datumtrans(int datumtrans, int force)
 			    temp_projinfo);
 	}
 
-	/* Finally add new parameters (if we have them) */
+	/* Finally (re-)add datum name and add new parameters (if we have them) */
+	G_set_key_value("datum", datum, temp_projinfo);
+	G_message(_("Datum set to <%s>"), datum);
 	if (chosenparams != NULL) {
 	    /* Now split 'chosenparams' into key/value format */
 	    paramkey = strtok(chosenparams, "=");

general/g.proj/g.proj.html

@@ -7,8 +7,8 @@ functionality is limited to:
 <ul>
 <li>Reporting the projection information for the current location, 
 either in conventional GRASS (-p flag) or PROJ.4 (-j flag) format</li>
-<li>Reporting and modifying the datum transformation parameters for
-the current location</li>
+<li>Changing the datum, or reporting and modifying the datum transformation 
+parameters, for the current location</li>
 </ul>
 
 <p>When compiled with OGR, functionality is increased and allows output of 
@@ -33,7 +33,7 @@ file.</dd>
 
 <dt>proj4=<em>description</em> or <em>-</em></dt>
 <dd><em>description</em> should be a projection description in 
-<a href="http://remotesensing.org/proj/">PROJ.4</a> format, enclosed in
+<a href="http://proj.osgeo.org/">PROJ.4</a> format, enclosed in
 quotation marks if there are any spaces. If <em>-</em> is given for
 <em>description</em>, the PROJ.4 description will be read from stdin rather 
 than as a directly-supplied command-line parameter.</dd>
@@ -47,6 +47,18 @@ co-ordinate system and datum information files, stored in the directory
 to support future revisions of the EPSG database.</dd>
 </dl>
 
+<p>If datum information is incorrect or missing in the input
+co-ordinate system definition (e.g. PROJ.4 descriptions have very limited
+support for specifying datum names), a GRASS datum abbreviation can instead be
+supplied using the <em>datum</em> parameter. This will override any
+datum contained in the input co-ordinate system, and discard
+any datum transformation parameters. Enter datum=<em>list</em> to return a
+list of all the datums supported by GRASS. Note that the -t flag is
+implicit when a value for <em>datum</em> is specified. Since any
+existing datum transformation parameters will have been discarded, the
+<em>datumtrans</em> parameter should in general always be used in
+conjunction with <em>datum</em>.
+
 <p>The -p, -j, -w, etc. flags are all functional when importing projection
 information from an external source, meaning that <em>g.proj</em> can be
 used to convert between representations of the information. It is
@@ -167,6 +179,14 @@ stored in a text file:<br>
 g.proj -c wkt=irish_grid.prj location=irish_grid
 </pre></div>
 
+<p>Create a new location from a PROJ.4 description, explicitly
+specifying a datum and using the default datum transformation
+parameters:<br>
+
+<div class="code"><pre>
+g.proj -c location=spain proj4="+proj=utm +zone=30 +ellps=intl" datum=eur50 datumtrans=0
+</pre></div>
+
 <p>Reproject external raster map to current GRASS projection (does not always make sense!)
 using the GDAL 'gdalwarp' tool. We recommend to use the ERDAS/Img format and not
 to use the ESRI style of WKT:<br>

general/g.proj/local_proto.h

@@ -21,7 +21,7 @@ void print_wkt(int, int);
 #endif
 
 /* datumtrans.c */
-int set_datumtrans(int, int);
+int set_datumtrans(char *, int, int);
 
 /* create.c */
 void create_location(char *);

general/g.proj/main.c

@@ -48,6 +48,7 @@ int main(int argc, char *argv[])
 	*ingeo,			/* Input geo-referenced file readable by 
 				 * GDAL or OGR                              */
 #endif
+	*datum,			/* datum to add (or replace existing datum) */
 	*dtrans;		/* index to datum transform option          */
     struct GModule *module;
     
@@ -156,6 +157,17 @@ int main(int argc, char *argv[])
     inepsg->description = _("EPSG projection code");
 #endif
 
+    datum = G_define_option();
+    datum->key = "datum";
+    datum->type = TYPE_STRING;
+    datum->key_desc = "name";
+    datum->required = NO;
+    datum->guisection = _("Datum");
+    datum->label =
+	_("Datum (overrides any datum specified in input co-ordinate system)");
+    datum->description = 
+	_("Accepts standard GRASS datum codes, or \"list\" to list and exit");
+
     dtrans = G_define_option();
     dtrans->key = "datumtrans";
     dtrans->type = TYPE_INTEGER;
@@ -204,6 +216,20 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Only one of '%s', '%s', '%s' or '%s' options may be specified"),
 		      ingeo->key, inwkt->key, inproj4->key, inepsg->key);
 
+    /* List supported datums if requested; code originally 
+     * from G_ask_datum_name() (formerly in libgis) */
+    if (datum->answer && strcmp(datum->answer, "list") == 0) {
+	const char *dat;
+	int i;
+
+	for (i = 0; (dat = G_datum_name(i)); i++) {	
+	    fprintf(stdout, "---\n%d\n%s\n%s\n%s ellipsoid\n",
+		    i, dat, G_datum_description(i), G_datum_ellipsoid(i));
+	}
+
+	exit(EXIT_SUCCESS);
+    }
+
     /* Input */
     /* We can only have one input source, hence if..else construct */
 
@@ -233,7 +259,7 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Projection files missing"));
 
     /* Set Datum Parameters if necessary or requested */
-    set_datumtrans(atoi(dtrans->answer), forcedatumtrans->answer);
+    set_datumtrans(datum->answer, atoi(dtrans->answer), forcedatumtrans->answer);
 
 
     /* Output */