[GRASS-SVN] r53453 - in grass/branches/develbranch_6: general/g.proj include lib/proj

svn_grass at osgeo.org svn_grass at osgeo.org
Wed Oct 17 13:20:25 PDT 2012


Author: neteler
Date: 2012-10-17 13:20:24 -0700 (Wed, 17 Oct 2012)
New Revision: 53453

Added:
   grass/branches/develbranch_6/general/g.proj/create.c
Modified:
   grass/branches/develbranch_6/general/g.proj/datumtrans.c
   grass/branches/develbranch_6/general/g.proj/description.html
   grass/branches/develbranch_6/general/g.proj/input.c
   grass/branches/develbranch_6/general/g.proj/local_proto.h
   grass/branches/develbranch_6/general/g.proj/main.c
   grass/branches/develbranch_6/general/g.proj/output.c
   grass/branches/develbranch_6/include/gprojects.h
   grass/branches/develbranch_6/lib/proj/datum.c
Log:
backport with GRASS 6 adaptations of datum handling patches by Paul Kelly: r53297, r53305, r53310 from trunk

Added: grass/branches/develbranch_6/general/g.proj/create.c
===================================================================
--- grass/branches/develbranch_6/general/g.proj/create.c	                        (rev 0)
+++ grass/branches/develbranch_6/general/g.proj/create.c	2012-10-17 20:20:24 UTC (rev 53453)
@@ -0,0 +1,68 @@
+#include <errno.h>
+#include <string.h>
+
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "local_proto.h"
+
+void create_location(char *location)
+{
+    int ret;
+
+    ret = G__make_location(location, &cellhd, projinfo, projunits, NULL);
+    if (ret == 0)
+	G_message(_("Location <%s> created"), location);
+    else if (ret == -1)
+	G_fatal_error(_("Unable to create location <%s>: %s"),
+		    location, strerror(errno));
+    else if (ret == -2)
+	G_fatal_error(_("Unable to create projection files: %s"),
+		    strerror(errno));
+    else
+	/* Shouldn't happen */
+	G_fatal_error(_("Unspecified error while creating new location"));
+
+    G_message(_("You can switch to the new location by\n`%s=%s`"),
+	      "g.mapset mapset=PERMANENT location", location);
+}
+
+void modify_projinfo()
+{
+    const char *mapset = G_mapset();
+    struct Cell_head old_cellhd;
+    
+    if (strcmp(mapset, "PERMANENT") != 0)
+	G_fatal_error(_("You must select the PERMANENT mapset before updating the "
+			"current location's projection (current mapset is <%s>)."),
+		      mapset);
+    
+    /* Read projection information from current location first */
+    G_get_default_window(&old_cellhd);
+    
+    char path[GPATH_MAX];
+    int stat;
+	
+    /* Write out the PROJ_INFO, and PROJ_UNITS if available. */
+    if (projinfo != NULL) {
+	G__file_name(path, "", "PROJ_INFO", "PERMANENT");
+	G_write_key_value_file(path, projinfo, &stat);
+    }
+    
+    if (projunits != NULL) {
+	G__file_name(path, "", "PROJ_UNITS", "PERMANENT");
+	G_write_key_value_file(path, projunits, &stat);
+    }
+    
+    if ((old_cellhd.zone != cellhd.zone) ||
+	(old_cellhd.proj != cellhd.proj)) {
+	/* Recreate the default, and current window files if projection
+	 * number or zone have changed */
+	G__put_window(&cellhd, "", "DEFAULT_WIND");
+	G__put_window(&cellhd, "", "WIND");
+	G_message(_("Default region was updated to the new projection, but if you have "
+		    "multiple mapsets `g.region -d` should be run in each to update the "
+		    "region from the default"));
+    }
+    G_important_message(_("Projection information updated"));
+}

Modified: grass/branches/develbranch_6/general/g.proj/datumtrans.c
===================================================================
--- grass/branches/develbranch_6/general/g.proj/datumtrans.c	2012-10-17 20:15:48 UTC (rev 53452)
+++ grass/branches/develbranch_6/general/g.proj/datumtrans.c	2012-10-17 20:20:24 UTC (rev 53453)
@@ -16,13 +16,77 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 #include <grass/gis.h>
+#include <grass/glocale.h>
 #include <grass/gprojects.h>
 
 #include "local_proto.h"
 
 /**
  * 
+ * \brief Add or replace datum to the current co-ordinate system definition
+ * 
+ * \param datum       Use this datum (overrides any datum found in the
+ *                    current co-ordinate system definition).
+ * 
+ * \return            1 if a change was made, 0 if not.
+ **/
+
+int set_datum(char *datum)
+{
+    struct gpj_datum dstruct;
+    struct Key_Value *temp_projinfo;
+    int i;
+
+    if (cellhd.proj == PROJECTION_XY)
+	return 0;
+
+    if (!datum || GPJ_get_datum_by_name(datum, &dstruct) < 0)
+    {
+	G_fatal_error(_("Invalid datum code <%s>"), datum);
+	return 0;
+    }
+
+    temp_projinfo = G_create_key_value();
+
+    /* Copy old PROJ_INFO, skipping out any keys related
+     * to datum or ellipsoid parameters */
+    for (i = 0; i < projinfo->nitems; i++) {
+        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
+	    || strcmp(projinfo->key[i], "nadgrids") == 0
+	    || strcmp(projinfo->key[i], "towgs84") == 0
+	    || strcmp(projinfo->key[i], "ellps") == 0
+	    || strcmp(projinfo->key[i], "a") == 0
+	    || strcmp(projinfo->key[i], "b") == 0
+	    || strcmp(projinfo->key[i], "es") == 0
+	    || strcmp(projinfo->key[i], "f") == 0
+	    || strcmp(projinfo->key[i], "rf") == 0)
+	    continue;
+
+	G_set_key_value(projinfo->key[i], projinfo->value[i],
+			temp_projinfo);
+    }
+
+    /* Finally add datum and ellipsoid names */
+    G_set_key_value("datum", dstruct.name, temp_projinfo);
+    G_message(_("Datum set to <%s>"), dstruct.name);
+    G_set_key_value("ellps", dstruct.ellps, temp_projinfo);
+    G_message(_("Ellipsoid set to <%s>"), dstruct.ellps);
+
+    /* Destroy original key/value structure and replace with new one */
+    G_free_key_value(projinfo);
+    projinfo = temp_projinfo;
+
+    return 1;
+}
+
+/**
+ * 
  * \brief Add or replace datum transformation parameters to the current
  *        co-ordinate system definition
  * 
@@ -34,27 +98,25 @@
  * \param force       Force editing of parameters even if current co-ordinate
  *                    system already contains fully specified parameters.
  * 
- * \param interactive Interactively prompt user for parameters, ignoring
- *                    value of datumtrans variable.
- * 
  * \return            1 if a change was made, 0 if not.
  **/
 
-int set_datumtrans(int datumtrans, int force, int interactive)
+int set_datumtrans(int datumtrans, int force)
 {
     char *params, *datum = NULL;
     int paramsets, status;
-    struct gpj_datum dstruct;
 
     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) {
 	/* A datum name is specified; need to determine if
 	 * there are parameters to choose from for this datum */
+	struct gpj_datum dstruct;
 
 	if (GPJ_get_datum_by_name(datum, &dstruct) > 0) {
 	    char *defparams;
@@ -62,6 +124,8 @@
 	    paramsets =
 		GPJ_get_default_datum_params_by_name(dstruct.name,
 						     &defparams);
+	    G_free(defparams);
+	    GPJ_free_datum(&dstruct);
 
 	    G_debug(3, "set_datumtrans(): datum transform terms found "
 		    "with %d options", paramsets);
@@ -89,15 +153,9 @@
 	struct Key_Value *temp_projinfo;
 	int i;
 
-	/* First of all obtain the new parameters either interactively
-	 * or through the supplied transform number index */
-	if (interactive) {
-	    /* Print to stderr here as the function call will also */
-	    fprintf(stderr, "Datum name: %s (%s)", dstruct.name,
-		    dstruct.longname);
-	    GPJ_ask_datum_params(datum, &chosenparams);
-	}
-	else {
+	/* First of all obtain the new parameters 
+	 * through the supplied transform number index */
+	{
 	    struct gpj_datum_transform_list *list;
 
 	    if (datumtrans > paramsets)
@@ -113,21 +171,23 @@
 	    if (list != NULL) {
 		if (datumtrans == -1) {
 		    do {
+			struct gpj_datum_transform_list *old = list;
 			fprintf(stdout, "---\n%d\nUsed in %s\n%s\n%s\n",
 				list->count, list->where_used,
 				list->params, list->comment);
 			list = list->next;
+		        GPJ_free_datum_transform(old);
 		    } while (list != NULL);
 
 		    exit(EXIT_SUCCESS);
 		}
 		else {
 		    do {
-			if (list->count == datumtrans) {
+			struct gpj_datum_transform_list *old = list;
+			if (list->count == datumtrans)
 			    chosenparams = G_store(list->params);
-			    break;
-			}
 			list = list->next;
+		        GPJ_free_datum_transform(old);
 		    } while (list != NULL);
 		}
 	    }
@@ -165,6 +225,7 @@
 	projinfo = temp_projinfo;
 
     }
+    G_free(datum);
 
     return force;
 }

Modified: grass/branches/develbranch_6/general/g.proj/description.html
===================================================================
--- grass/branches/develbranch_6/general/g.proj/description.html	2012-10-17 20:15:48 UTC (rev 53452)
+++ grass/branches/develbranch_6/general/g.proj/description.html	2012-10-17 20:20:24 UTC (rev 53453)
@@ -1,18 +1,24 @@
 <h2>DESCRIPTION</h2>
 
-<p><em>g.proj</em> provides a means of converting a co-ordinate system
+<em>g.proj</em> provides a means of converting a co-ordinate system
 description (i.e. projection information) between various formats.
-It requires <a href="http://www.gdal.org/ogr/">OGR</a> to compile. The basic
-functionality of the module is to report the projection information for the
-current location, either in conventional GRASS (-p flag) or PROJ.4 (-j flag) 
-format.
+If compiled without <a href="http://www.gdal.org/ogr/">OGR</a> present, the 
+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>Changing the datum, or reporting and modifying the datum transformation 
+parameters, for the current location</li>
+</ul>
 
-<p>Projection information may also be output in the Well-Known Text (WKT) 
-format popularised
+<p>
+When compiled with OGR, functionality is increased and allows output of 
+the projection information in the Well-Known Text (WKT) format popularised 
 by proprietary GIS. In addition, if one of the parameters <em>georef</em>, 
 <em>wkt</em>, <em>proj4</em> or <em>epsg</em> is specified, rather than the 
 projection information being read from the current location it is imported 
 from an external source as follows:
+
 <dl>
 <dt>georef=<em>filename</em></dt>
 <dd><em>g.proj</em> attempts to invoke GDAL and OGR in turn to read a
@@ -28,7 +34,7 @@
 
 <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>
@@ -38,36 +44,51 @@
 co-ordinate system in the <a href="http://www.epsg.org/CurrentDB.html">EPSG 
 database</a>. EPSG code support is based upon a local copy of the GDAL CSV 
 co-ordinate system and datum information files, stored in the directory 
-${GISBASE}/etc/ogr_csv. These can be updated if necessary to support future
-revisions of the EPSG database.</dd>
+<tt>$GISBASE/etc/proj/ogr_csv</tt>. These can be updated if necessary
+to support future revisions of the EPSG database.</dd>
 </dl>
 
-<p>The -p, -j, -w, etc. flags are all functional when importing projection
+<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. 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
 <strong>not</strong> required that either the input or output be in GRASS
 format.
 
-<p>In addition however, if the -c flag is specified, <em>g.proj</em> will 
+<p>
+In addition however, if the -c flag is specified, <em>g.proj</em> will 
 create new GRASS projection files (PROJ_INFO, PROJ_UNITS, WIND and 
 DEFAULT_WIND) based on the imported information. If the <em>location</em> 
 parameter is specified in addition to -c, then a new location will be created. 
 Otherwise the projection information files in the current location will be
-overwritten. The program will warn before doing this only if command-line
-interactive mode (<em>-i</em> flag) is selected.
+overwritten. The program will <strong>not</strong> warn before doing this.
 
-<p>The final mode of operation of g.proj is to report on the datum
+<p>
+The final mode of operation of g.proj is to report on the datum
 information and datum transformation parameters associated with the
 co-ordinate system. The -d flag will report a human-readable summary of
 this.
 
 <h2>NOTES</h2>
 
-<p>If the input co-ordinate system contains a datum name but no
+<p>
+If the input co-ordinate system contains a datum name but no
 transformation parameters, and there is more than one suitable parameter set
 available (according to the files datum.table and datumtransform.table in
-${GISBASE}/etc), g.proj will check the value of the <em>datumtrans</em>
-option and act according to the following:<br>
+<tt>$GISBASE/etc/proj</tt>), <em>g.proj</em> will check the value of
+the <em>datumtrans</em> option and act according to the following:<br>
 <strong>-1:</strong> List available parameter sets in a GUI-parsable (but also
 human-readable) format and exit.<br>
 <strong>0 (default):</strong> Continue without specifying parameters - if 
@@ -76,24 +97,24 @@
 <strong>Any other number less than or equal to the number of parameter sets
 available for this datum:</strong> Choose this parameter set and add it to the
 co-ordinate system description.<br>
-If the module is being used from the command-line through an interactive
-terminal, the <em>-i</em> flag can be specified to enable interactive
-selection of the parameter set, and the value of <em>datumtrans</em> (if 
-specified) is ignored.<br>
 If the <em>-t</em> flag is specified, the module will attempt to change the
 datum transformation parameters using one of the above two methods 
 <strong>even if</strong> a valid parameter set is already specified in the 
-input co-ordinate system.
+input co-ordinate system. This can be useful to change the datum information
+for an existing location.
 
-<p>Output is simply based on the input projection information. g.proj does 
+<p>
+Output is simply based on the input projection information. g.proj does 
 <strong>not</strong> attempt to verify that the co-ordinate system thus 
 described matches an existing system in use in the world. In particular,
 this means there are no EPSG Authority codes in the WKT output.
 
-<p>WKT format shows the false eastings and northings in the projected unit
+<p>
+WKT format shows the false eastings and northings in the projected unit
 (e.g. meters, feet) but in PROJ format it should always be given in meters.
 
-<p>The maximum size of input WKT or PROJ.4 projection descriptions is
+<p>
+The maximum size of input WKT or PROJ.4 projection descriptions is
 limited to 8000 bytes.
 
 <h2>EXAMPLES</h2>
@@ -105,6 +126,7 @@
 </pre></div>
 
 <p>
+
 Create a '.prj' file in ESRI format corresponding to the current location:<br>
 
 <div class="code"><pre>
@@ -154,18 +176,17 @@
 
 <p>
 Create a new location with the same co-ordinate system as the current
-location:<br>
+location, but forcing a change to datum transformation parameter set no. 1:<br>
 
 <div class="code"><pre>
-g.proj -c location=newloc
+g.proj -c location=newloc -t datumtrans=1
 </pre></div>
 
 <p>
-Interactively change/update the datum transformation parameters for the
-current location:<br>
+List the possible datum transformation parameters for the current location:<br>
 
 <div class="code"><pre>
-g.proj -itc
+g.proj -t datumtrans=-1
 </pre></div>
 
 <p>
@@ -177,6 +198,15 @@
 </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>
@@ -217,4 +247,5 @@
 
 Paul Kelly
 
-<p><i>Last changed: $Date$</i>
+<p>
+<i>Last changed: $Date$</i>

Modified: grass/branches/develbranch_6/general/g.proj/input.c
===================================================================
--- grass/branches/develbranch_6/general/g.proj/input.c	2012-10-17 20:15:48 UTC (rev 53452)
+++ grass/branches/develbranch_6/general/g.proj/input.c	2012-10-17 20:20:24 UTC (rev 53453)
@@ -22,15 +22,22 @@
 #include <grass/gis.h>
 #include <grass/gprojects.h>
 #include <grass/glocale.h>
-#include <gdal.h>
-#include <ogr_api.h>
-#include <cpl_csv.h>
+#include <grass/config.h>
 
+#ifdef HAVE_OGR
+#  include <gdal.h>
+#  include <ogr_api.h>
+#  include <cpl_csv.h>
+#endif
+
 #include "local_proto.h"
 
 static void set_default_region(void);
+
+#ifdef HAVE_OGR
 static void set_gdal_region(GDALDatasetH);
 static void set_ogr_region(OGRLayerH);
+#endif
 
 /**
  * \brief Read projection and region information from current location
@@ -50,6 +57,8 @@
     return;
 }
 
+#ifdef HAVE_OGR
+
 /**
  * \brief Read projection information in WKT format from stdin or a file
  * 
@@ -86,7 +95,7 @@
 	G_squeeze(buff);
     }
     else
-	G_fatal_error(_("Unable to open file [%s] for reading"), wktfile);
+	G_fatal_error(_("Unable to open file '%s' for reading"), wktfile);
 
     ret = GPJ_wkt_to_grass(&cellhd, &projinfo, &projunits, buff, 0);
     set_default_region();
@@ -254,6 +263,7 @@
     return ret;
 
 }
+#endif /* HAVE_OGR */
 
 /**
  * \brief Populates global cellhd with "default" region settings
@@ -283,6 +293,8 @@
     return;
 }
 
+#ifdef HAVE_OGR
+
 /**
  * \brief Populates global cellhd with region settings based on 
  *        georeferencing information in a GDAL dataset
@@ -365,3 +377,4 @@
 
     return;
 }
+#endif /* HAVE_OGR */

Modified: grass/branches/develbranch_6/general/g.proj/local_proto.h
===================================================================
--- grass/branches/develbranch_6/general/g.proj/local_proto.h	2012-10-17 20:15:48 UTC (rev 53452)
+++ grass/branches/develbranch_6/general/g.proj/local_proto.h	2012-10-17 20:20:24 UTC (rev 53453)
@@ -1,25 +1,29 @@
-#ifdef G_PROJ_MAIN
-#  define G_PROJ_GLOBAL
-#else
-#  define G_PROJ_GLOBAL extern
-#endif
+#include <grass/config.h>
 
-G_PROJ_GLOBAL struct Key_Value *projinfo, *projunits;
-G_PROJ_GLOBAL struct Cell_head cellhd;
+extern struct Key_Value *projinfo, *projunits;
+extern struct Cell_head cellhd;
 
 /* input.c */
 void input_currloc(void);
+#ifdef HAVE_OGR
 int input_wkt(char *);
 int input_proj4(char *);
 int input_epsg(int);
 int input_georef(char *);
+#endif
 
 /* output.c */
-void print_projinfo(void);
+void print_projinfo(int);
 void print_datuminfo(void);
 void print_proj4(int);
+#ifdef HAVE_OGR
 void print_wkt(int, int);
-void create_location(char *, int);
+#endif
 
 /* datumtrans.c */
-int set_datumtrans(int, int, int);
+int set_datum(char *);
+int set_datumtrans(int, int);
+
+/* create.c */
+void create_location(char *);
+void modify_projinfo();

Modified: grass/branches/develbranch_6/general/g.proj/main.c
===================================================================
--- grass/branches/develbranch_6/general/g.proj/main.c	2012-10-17 20:15:48 UTC (rev 53452)
+++ grass/branches/develbranch_6/general/g.proj/main.c	2012-10-17 20:20:24 UTC (rev 53453)
@@ -3,14 +3,15 @@
  *
  * MODULE:       g.proj 
  * AUTHOR(S):    Paul Kelly - paul-grass at stjohnspoint.co.uk
+ *               Shell script style by Martin Landa <landa.martin gmail.com>
  * PURPOSE:      Provides a means of reporting the contents of GRASS
  *               projection information files and creating
  *               new projection information files.
- * COPYRIGHT:    (C) 2003-2007 by the GRASS Development Team
+ * COPYRIGHT:    (C) 2003-2007, 2011 by the GRASS Development Team
  *
- *               This program is free software under the GNU General Public
- *               License (>=v2). Read the file COPYING that comes with GRASS
- *               for details.
+ *               This program is free software under the GNU General
+ *               Public License (>=v2). Read the file COPYING that
+ *               comes with GRASS for details.
  *
  *****************************************************************************/
 
@@ -18,33 +19,41 @@
 #include <string.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
+#include <grass/config.h>
 
-#define G_PROJ_MAIN
 #include "local_proto.h"
 
+struct Key_Value *projinfo, *projunits;
+struct Cell_head cellhd;
+
 int main(int argc, char *argv[])
 {
     struct Flag *printinfo,	/* Print contents of PROJ_INFO & PROJ_UNITS */
-     *printproj4,		/* Print projection in PROJ.4 format        */
-     *datuminfo,		/* Check if datum information is present    */
-     *create,			/* Create new projection files              */
-     *printwkt,			/* Print projection in WKT format           */
-     *esristyle,		/* Use ESRI-style WKT format                */
-     *dontprettify,		/* Print 'flat' output (no linebreaks)      */
-     *forcedatumtrans,		/* Force override of datumtrans parameters  */
-     *interprompt;		/* Enable interactive prompting             */
-
+	*shellinfo,             /* Print in shell script style              */
+	*printproj4,		/* Print projection in PROJ.4 format        */
+	*datuminfo,		/* Check if datum information is present    */
+	*create,		/* Create new projection files              */
+#ifdef HAVE_OGR
+	*printwkt,		/* Print projection in WKT format           */
+	*esristyle,		/* Use ESRI-style WKT format                */
+#endif
+	*dontprettify,		/* Print 'flat' output (no linebreaks)      */
+	*forcedatumtrans;	/* Force override of datumtrans parameters  */
+    
     struct Option *location,	/* Name of new location to create           */
-     *inepsg,			/* EPSG projection code                     */
-     *dtrans,			/* index to datum transform option          */
-     *inwkt,			/* Input file with projection in WKT format */
-     *inproj4,			/* Projection in PROJ.4 format              */
-     *ingeo;			/* Input geo-referenced file readable by 
+#ifdef HAVE_OGR
+	*inepsg,		/* EPSG projection code                     */
+	*inwkt,			/* Input file with projection in WKT format */
+	*inproj4,		/* Projection in PROJ.4 format              */
+	*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;
+    
+    int formats;
 
-    int importformats;
-
     G_set_program_name(argv[0]);
     G_no_gisinit();		/* We don't call G_gisinit() here because it validates the
 				 * mapset, whereas this module may legitmately be used 
@@ -52,18 +61,29 @@
 
     module = G_define_module();
     module->keywords = _("general, projection");
+#ifdef HAVE_OGR
     module->label =
-	_("Converts co-ordinate system descriptions (i.e. projection "
-	  "information) between various formats (including GRASS format).");
+	_("Prints and manipulates GRASS projection information files "
+	  "(in various co-ordinate system descriptions).");
     module->description =
-	_("Can also be used to create GRASS locations.");
+	_("Can also be used to create new GRASS locations.");
+#else
+    module->description =
+	_("Prints and manipulates GRASS projection information files.");
+#endif
 
     printinfo = G_define_flag();
     printinfo->key = 'p';
     printinfo->guisection = _("Print");
     printinfo->description =
-	_("Print projection information (in conventional GRASS format)");
+	_("Print projection information in conventional GRASS format");
 
+    shellinfo = G_define_flag();
+    shellinfo->key = 'g';
+    shellinfo->guisection = _("Print");
+    shellinfo->description =
+	_("Print projection information in shell script style");
+
     datuminfo = G_define_flag();
     datuminfo->key = 'd';
     datuminfo->guisection = _("Print");
@@ -76,6 +96,17 @@
     printproj4->description =
 	_("Print projection information in PROJ.4 format");
 
+    dontprettify = G_define_flag();
+    dontprettify->key = 'f';
+    dontprettify->guisection = _("Print");
+    dontprettify->description =
+	_("Print 'flat' output with no linebreaks (applies to "
+#ifdef HAVE_OGR
+	  "WKT and "
+#endif
+	  "PROJ.4 output)");
+
+#ifdef HAVE_OGR
     printwkt = G_define_flag();
     printwkt->key = 'w';
     printwkt->guisection = _("Print");
@@ -86,21 +117,14 @@
     esristyle->guisection = _("Print");
     esristyle->description =
 	_("Use ESRI-style format (applies to WKT output only)");
-
-    dontprettify = G_define_flag();
-    dontprettify->key = 'f';
-    dontprettify->guisection = _("Print");
-    dontprettify->description =
-	_("Print 'flat' output with no linebreaks (applies to WKT and "
-	  "PROJ.4 output)");
-
+    
     ingeo = G_define_option();
     ingeo->key = "georef";
     ingeo->type = TYPE_STRING;
     ingeo->key_desc = "file";
     ingeo->required = NO;
-    ingeo->guisection = _("Input");
-    ingeo->description = _("Georeferenced data file to read projection "
+    ingeo->guisection = _("Specification");
+    ingeo->description = _("Name of georeferenced data file to read projection "
 			   "information from");
 
     inwkt = G_define_option();
@@ -108,29 +132,45 @@
     inwkt->type = TYPE_STRING;
     inwkt->key_desc = "file";
     inwkt->required = NO;
-    inwkt->guisection = _("Input");
-    inwkt->description = _("ASCII file containing a WKT projection "
-			   "description (- for stdin)");
+    inwkt->guisection = _("Specification");
+    inwkt->label = _("Name of ASCII file containing a WKT projection "
+		     "description");
+    inwkt->description = _("'-' for standard input");
 
     inproj4 = G_define_option();
     inproj4->key = "proj4";
     inproj4->type = TYPE_STRING;
     inproj4->key_desc = "params";
     inproj4->required = NO;
-    inproj4->guisection = _("Input");
-    inproj4->description = _("PROJ.4 projection description (- for stdin)");
+    inproj4->guisection = _("Specification");
+    inproj4->label = _("PROJ.4 projection description");
+    inproj4->description = _("'-' for standard input");
 
     inepsg = G_define_option();
     inepsg->key = "epsg";
     inepsg->type = TYPE_INTEGER;
+    inepsg->key_desc = "code";
     inepsg->required = NO;
     inepsg->options = "1-1000000";
-    inepsg->guisection = _("Input");
+    inepsg->guisection = _("Specification");
     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;
+    dtrans->key_desc = "index";
     dtrans->required = NO;
     dtrans->options = "-1-100";
     dtrans->answer = "0";
@@ -147,48 +187,57 @@
 
     create = G_define_flag();
     create->key = 'c';
-    create->guisection = _("Create/Edit");
+    create->guisection = _("Modify");
     create->description = _("Create new projection files (modifies current "
-			    "location unless 'location' option specified)");
+			    "location)");
 
     location = G_define_option();
     location->key = "location";
     location->type = TYPE_STRING;
     location->key_desc = "name";
     location->required = NO;
-    location->guisection = _("Create/Edit");
+    location->guisection = _("Create");
     location->description = _("Name of new location to create");
 
-    interprompt = G_define_flag();
-    interprompt->key = 'i';
-    interprompt->description =
-	_("Enable interactive prompting (for command-line use only)");
-
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
 
     /* Initialisation & Validation */
 
+#ifdef HAVE_OGR
     /* -e implies -w */
     if (esristyle->answer && !printwkt->answer)
 	printwkt->answer = 1;
 
-    projinfo = NULL;
-    projunits = NULL;
-
-    importformats = ((ingeo->answer ? 1 : 0) + (inwkt->answer ? 1 : 0) +
-		     (inproj4->answer ? 1 : 0) + (inepsg->answer ? 1 : 0));
-    if (importformats > 1)
+    formats = ((ingeo->answer ? 1 : 0) + (inwkt->answer ? 1 : 0) +
+	       (inproj4->answer ? 1 : 0) + (inepsg->answer ? 1 : 0));
+    if (formats > 1)
 	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 */
 
-    if (importformats == 0)
+    if (formats == 0)
+#endif
 	/* Input is projection of current location */
 	input_currloc();
+#ifdef HAVE_OGR
     else if (inwkt->answer)
 	/* Input in WKT format */
 	input_wkt(inwkt->answer);
@@ -201,41 +250,93 @@
     else
 	/* Input from georeferenced file */
 	input_georef(ingeo->answer);
+#endif
 
-
     /* Consistency Check */
 
     if ((cellhd.proj != PROJECTION_XY)
 	&& (projinfo == NULL || projunits == NULL))
 	G_fatal_error(_("Projection files missing"));
 
+    /* Override input datum if requested */
+    if(datum->answer)
+	set_datum(datum->answer);
+
     /* Set Datum Parameters if necessary or requested */
-    set_datumtrans(atoi(dtrans->answer), forcedatumtrans->answer,
-		   interprompt->answer);
+    set_datumtrans(atoi(dtrans->answer), forcedatumtrans->answer);
 
 
     /* Output */
-    /* We can output the same information in multiple formats if
-     * necessary, hence multiple if statements */
+    /* Only allow one output format at a time, to reduce confusion */
+    formats = ((printinfo->answer ? 1 : 0) + (shellinfo->answer ? 1 : 0) +
+	       (datuminfo->answer ? 1 : 0) +
+	       (printproj4->answer ? 1 : 0) +
+#ifdef HAVE_OGR
+	       (printwkt->answer ? 1 : 0) +
+#endif
+	       (create->answer ? 1 : 0));
+    if (formats > 1)
+	G_fatal_error(_("Only one of -%c, -%c, -%c, -%c"
+#ifdef HAVE_OGR
+			", -%c"
+#endif
+			" or -%c flags may be specified"),
+		      printinfo->key, shellinfo->key, datuminfo->key, printproj4->key,
+#ifdef HAVE_OGR
+		      printwkt->key,
+#endif
+		      create->key);
 
-    if (printinfo->answer)
-	print_projinfo();
-
-    if (datuminfo->answer)
+    if (printinfo->answer || shellinfo->answer)
+	print_projinfo(shellinfo->answer);
+    else if (datuminfo->answer)
 	print_datuminfo();
-
-    if (printproj4->answer)
+    else if (printproj4->answer)
 	print_proj4(dontprettify->answer);
-
-    if (printwkt->answer)
+#ifdef HAVE_OGR
+    else if (printwkt->answer)
 	print_wkt(esristyle->answer, dontprettify->answer);
+#endif
+    else if (location->answer)
+	create_location(location->answer);
+    else if (create->answer)
+	modify_projinfo();
+    else
+#ifdef HAVE_OGR
+	G_fatal_error(_("No output format specified, define one "
+			"of flags -%c, -%c, -%c, or -%c"),
+		      printinfo->key, shellinfo->key, printproj4->key, printwkt->key);
+#else
+	G_fatal_error(_("No output format specified, define one "
+			"of flags -%c, -%c, or -%c"),
+		      printinfo->key, shellinfo->key, printproj4->key);
+#endif
 
-    if (create->answer)
-	create_location(location->answer, interprompt->answer);
+#ifdef HAVE_OGR
+    if (create->answer && inepsg->answer) {
+#else
+    if (create->answer){ 
+#endif
+	/* preserve epsg code for user records only (not used by grass's pj routines) */
+	FILE *fp;
+	char path[GPATH_MAX];
+	/* if inputs were not clean it should of failed by now */
+	if (location->answer) {
+            snprintf(path, sizeof(path), "%s/%s/%s/%s", G_gisdbase(),
+		     location->answer, "PERMANENT", "PROJ_EPSG");
+	    path[sizeof(path)-1] = '\0';
+	}
+	else
+	    G__file_name(path, "", "PROJ_EPSG", "PERMANENT");
+	fp = fopen(path, "w");
+#ifdef HAVE_OGR
+	fprintf(fp, "epsg: %s\n", inepsg->answer);
+#endif
+	fclose(fp);
+    }
 
 
     /* Tidy Up */
-
     if (projinfo != NULL)
 	G_free_key_value(projinfo);
     if (projunits != NULL)

Modified: grass/branches/develbranch_6/general/g.proj/output.c
===================================================================
--- grass/branches/develbranch_6/general/g.proj/output.c	2012-10-17 20:15:48 UTC (rev 53452)
+++ grass/branches/develbranch_6/general/g.proj/output.c	2012-10-17 20:20:24 UTC (rev 53453)
@@ -15,34 +15,70 @@
  *****************************************************************************/
 
 #include <stdio.h>
+#include <unistd.h>
 #include <string.h>
-#include <errno.h>
+#include <proj_api.h>
+
 #include <grass/gis.h>
 #include <grass/gprojects.h>
 #include <grass/glocale.h>
+#include <grass/config.h>
 
 #include "local_proto.h"
 
-static int check_xy(void);
+static int check_xy(int shell);
 
-void print_projinfo(void)
+void print_projinfo(int shell)
 {
     int i;
+    char path[GPATH_MAX];
 
-    if (check_xy())
+    if (check_xy(shell))
 	return;
 
-    fprintf(stdout,
-	    "-PROJ_INFO-------------------------------------------------\n");
-    for (i = 0; i < projinfo->nitems; i++)
-	fprintf(stdout, "%-11s: %s\n", projinfo->key[i], projinfo->value[i]);
+    if (!shell)
+	fprintf(stdout,
+		"-PROJ_INFO-------------------------------------------------\n");
+    for (i = 0; i < projinfo->nitems; i++) {
+	if (shell)
+	    fprintf(stdout, "%s=%s\n", projinfo->key[i], projinfo->value[i]);
+	else
+	    fprintf(stdout, "%-11s: %s\n", projinfo->key[i], projinfo->value[i]);
+    }
 
-    fprintf(stdout,
-	    "-PROJ_UNITS------------------------------------------------\n");
-    for (i = 0; i < projunits->nitems; i++)
-	fprintf(stdout, "%-11s: %s\n",
-		projunits->key[i], projunits->value[i]);
+    /* EPSG code is preserved for historical metadata interest only:
+	the contents of this file are not used by pj_*() routines at all */
+    G__file_name(path, "", "PROJ_EPSG", "PERMANENT");
+    if (access(path, F_OK) == 0) {
+	struct Key_Value *in_epsg_key;
+	int stat;
+	in_epsg_key = G_read_key_value_file(path, &stat);
+	if (!shell) {
+	    fprintf(stdout,
+		"-PROJ_EPSG-------------------------------------------------\n");
+	    fprintf(stdout, "%-11s: %s\n", in_epsg_key->key[0],
+		    in_epsg_key->value[0]);
+	}
+	else
+	    fprintf(stdout, "%s=%s\n", in_epsg_key->key[0],
+		    in_epsg_key->value[0]);
 
+	if (in_epsg_key != NULL)
+	    G_free_key_value(in_epsg_key);
+    }
+ 
+    if (!shell)
+	fprintf(stdout,
+		"-PROJ_UNITS------------------------------------------------\n");
+    for (i = 0; i < projunits->nitems; i++) {
+	if (shell)
+	    fprintf(stdout, "%s=%s\n",
+		    projunits->key[i], projunits->value[i]);
+	else
+	    fprintf(stdout, "%-11s: %s\n",
+		    projunits->key[i], projunits->value[i]);
+    }
+    
     return;
 }
 
@@ -52,7 +88,7 @@
     struct gpj_datum dstruct;
     int validdatum = 0;
 
-    if (check_xy())
+    if (check_xy(FALSE))
 	return;
 
     GPJ__get_datum_params(projinfo, &datum, &params);
@@ -93,13 +129,15 @@
 void print_proj4(int dontprettify)
 {
     struct pj_info pjinfo;
-    char *proj4, *proj4mod, *i, *unfact;
+    char *proj4, *proj4mod, *i;
+    const char *unfact;
 
-    if (check_xy())
+    if (check_xy(FALSE))
 	return;
 
     pj_get_kv(&pjinfo, projinfo, projunits);
     proj4 = pj_get_def(pjinfo.pj, 0);
+    pj_free(pjinfo.pj);
 
     /* GRASS-style PROJ.4 strings don't include a unit factor as this is
      * handled separately in GRASS - must include it here though */
@@ -107,7 +145,8 @@
     if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
 	G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
     else
-	proj4mod = proj4;
+	proj4mod = G_store(proj4);
+    pj_dalloc(proj4);
 
     for (i = proj4mod; *i; i++) {
 	/* Don't print the first space */
@@ -120,15 +159,17 @@
 	    fputc(*i, stdout);
     }
     fputc('\n', stdout);
+    G_free(proj4mod);
 
     return;
 }
 
+#ifdef HAVE_OGR
 void print_wkt(int esristyle, int dontprettify)
 {
     char *outwkt;
 
-    if (check_xy())
+    if (check_xy(FALSE))
 	return;
 
     outwkt = GPJ_grass_to_wkt(projinfo, projunits, esristyle,
@@ -142,121 +183,13 @@
 
     return;
 }
+#endif
 
-void create_location(char *location, int interactive)
+static int check_xy(int shell)
 {
-    int ret;
-
-    if (location) {
-	ret = G__make_location(location, &cellhd, projinfo, projunits, NULL);
-	if (ret == 0)
-	    G_message(_("Location %s created!"), location);
-	else if (ret == -1)
-	    G_fatal_error(_("Unable to create location: %s"),
-			  strerror(errno));
-	else if (ret == -2)
-	    G_fatal_error(_("Unable to create projection files: %s"),
-			  strerror(errno));
-	else
-	    /* Shouldn't happen */
-	    G_fatal_error(_("Unspecified error while creating new location"));
-    }
-    else {
-	/* Create flag given but no location specified; overwrite
-	 * projection files for current location */
-
-	int go_ahead = 0;
-	char *mapset = G_mapset();
-	struct Key_Value *old_projinfo = NULL, *old_projunits = NULL;
-	struct Cell_head old_cellhd;
-
-	if (strcmp(mapset, "PERMANENT") != 0)
-	    G_fatal_error(_("You must select the PERMANENT mapset before updating the "
-			   "current location's projection. (Current mapset is %s)"),
-			  mapset);
-
-	/* Read projection information from current location first */
-	G_get_default_window(&old_cellhd);
-
-	if (old_cellhd.proj != PROJECTION_XY) {
-	    old_projinfo = G_get_projinfo();
-	    old_projunits = G_get_projunits();
-	}
-
-	if (old_projinfo && old_projunits) {
-	    if (interactive) {
-		/* Warn as in g.setproj before overwriting current location */
-		fprintf(stderr,
-			_("\n\nWARNING!  A projection file already exists for this location\n"));
-		fprintf(stderr,
-			"\nThis file contains all the parameters for the\nlocation's projection: %s\n",
-			G_find_key_value("proj", old_projinfo));
-		fprintf(stderr,
-			"\n    Overriding this information implies that the old projection parameters\n");
-		fprintf(stderr,
-			"    were incorrect.  If you change the parameters, all existing data will be\n");
-		fprintf(stderr,
-			"    interpreted differently by the projection software.\n%c%c%c",
-			7, 7, 7);
-		fprintf(stderr,
-			"    GRASS will not re-project your data automatically\n\n");
-
-		if (G_yes
-		    (_("Would you still like to overwrite the current projection information "),
-		     0))
-		    go_ahead = 1;
-	    }
-	    else
-		go_ahead = 1;
-	}
-	else {
-	    /* Projection files missing for some reason;
-	     * go ahead and update */
-	    go_ahead = 1;
-	}
-
-	if (go_ahead) {
-	    int out_stat;
-	    char path[4096];
-
-	    /* Write out the PROJ_INFO, and PROJ_UNITS if available. */
-	    if (projinfo != NULL) {
-		G__file_name(path, "", "PROJ_INFO", "PERMANENT");
-		G_write_key_value_file(path, projinfo, &out_stat);
-		if (out_stat != 0)
-		    G_fatal_error(_("Error writing PROJ_INFO"));
-	    }
-
-	    if (projunits != NULL) {
-		G__file_name(path, "", "PROJ_UNITS", "PERMANENT");
-		G_write_key_value_file(path, projunits, &out_stat);
-		if (out_stat != 0)
-		    G_fatal_error(_("Error writing PROJ_UNITS"));
-	    }
-
-	    if ((old_cellhd.zone != cellhd.zone) ||
-		(old_cellhd.proj != cellhd.proj)) {
-		/* Recreate the default, and current window files if projection
-		 * number or zone have changed */
-		G__put_window(&cellhd, "", "DEFAULT_WIND");
-		G__put_window(&cellhd, "", "WIND");
-		G_message(_("N.B. The default region was updated to the new projection, but if you have "
-			   "multiple mapsets g.region -d should be run in each to update the region from "
-			   "the default."));
-	    }
-	    G_message(_("Projection information updated!"));
-	}
-	else
-	    G_message(_("The projection information will not be updated."));
-
-    }
-
-    return;
-}
-
-static int check_xy(void)
-{
     if (cellhd.proj == PROJECTION_XY) {
+	if (shell)
+	    fprintf(stdout, "name=");
 	fprintf(stdout, "XY location (unprojected)\n");
 	return 1;
     }

Modified: grass/branches/develbranch_6/include/gprojects.h
===================================================================
--- grass/branches/develbranch_6/include/gprojects.h	2012-10-17 20:15:48 UTC (rev 53452)
+++ grass/branches/develbranch_6/include/gprojects.h	2012-10-17 20:20:24 UTC (rev 53453)
@@ -95,6 +95,7 @@
 int GPJ_get_datum_params(char **, char **);
 int GPJ__get_datum_params(struct Key_Value *, char **, char **);
 void GPJ_free_datum(struct gpj_datum *);
+void GPJ_free_datum_transform(struct gpj_datum_transform_list *);
 int GPJ_ask_datum_params(const char *, char **);
 struct gpj_datum_transform_list *GPJ_get_datum_transform_by_name(const char
 								 *inputname);

Modified: grass/branches/develbranch_6/lib/proj/datum.c
===================================================================
--- grass/branches/develbranch_6/lib/proj/datum.c	2012-10-17 20:15:48 UTC (rev 53452)
+++ grass/branches/develbranch_6/lib/proj/datum.c	2012-10-17 20:20:24 UTC (rev 53453)
@@ -85,7 +85,7 @@
 int GPJ_get_default_datum_params_by_name(const char *name, char **params)
 {
     struct gpj_datum_transform_list *list, *old;
-    int count = 1;
+    int count = 0;
 
     list = GPJ_get_datum_transform_by_name(name);
 
@@ -98,14 +98,13 @@
      * (will normally be a 3-parameter transformation)        */
     *params = G_store(list->params);
 
-    while (list->next != NULL) {
+    while (list != NULL) {
 	count++;
 	old = list;
 	list = list->next;
-	G_free(old);
+	GPJ_free_datum_transform(old);
     }
 
-    G_free(list);
     return count;
 
 }
@@ -459,6 +458,21 @@
 }
 
 /**
+ * \brief Free the memory used by a gpj_datum_transform_list struct
+ *
+ * \param item gpj_datum_transform_list struct to be freed
+ **/
+
+void GPJ_free_datum_transform(struct gpj_datum_transform_list *item)
+{
+    G_free(item->params);
+    G_free(item->where_used);
+    G_free(item->comment);
+    G_free(item);
+    return;
+}
+
+/**
  * \brief Read the current GRASS datum.table from disk and store in
  *        memory
  * 



More information about the grass-commit mailing list