PROJ6+ / WKT2 support (#1240)

Co-authored-by: nilason <n_larsson@yahoo.com>

* adjust lib/gis, lib/proj, g.proj, r.proj, v.proj, r.in.gdal, r.external, r.out.gdal, v.in.ogr, v.external, v.out.ogr

general/g.proj/create.c

@@ -11,7 +11,7 @@ void create_location(const char *location)
     int ret;
 
     ret = G_make_location_crs(location, &cellhd, projinfo, projunits,
-                               projepsg, projwkt, projsrid);
+                              projsrid, projwkt);
     if (ret == 0)
 	G_message(_("Location <%s> created"), location);
     else if (ret == -1)
@@ -32,28 +32,28 @@ 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];
-	
+
     /* Write out the PROJ_INFO, PROJ_UNITS, and PROJ_EPSG if available. */
     if (projinfo != NULL) {
 	G_file_name(path, "", "PROJ_INFO", "PERMANENT");
 	G_write_key_value_file(path, projinfo);
     }
-    
+
     if (projunits != NULL) {
 	G_file_name(path, "", "PROJ_UNITS", "PERMANENT");
 	G_write_key_value_file(path, projunits);
     }
-    
+
     if (projepsg != NULL) {
 	G_file_name(path, "", "PROJ_EPSG", "PERMANENT");
 	G_write_key_value_file(path, projepsg);

general/g.proj/input.c

@@ -1,7 +1,7 @@
-/*  
+/*
  ****************************************************************************
  *
- * MODULE:       g.proj 
+ * MODULE:       g.proj
  * AUTHOR(S):    Paul Kelly - paul-grass@stjohnspoint.co.uk
  *               Frank Warmerdam
  *               Radim Blazek
@@ -36,7 +36,7 @@ static void set_default_region(void);
 
 /**
  * \brief Read projection and region information from current location
- * 
+ *
  * Reads projection and region information from current location and
  * stores in global structs cellhd, projinfo and projunits.
  **/
@@ -49,7 +49,7 @@ void input_currloc(void)
 	projwkt = G_get_projwkt();
 	projinfo = G_get_projinfo();
 	projunits = G_get_projunits();
-        projepsg = G_get_projepsg();
+        /* projepsg = G_get_projepsg(); */
     }
 
     return;
@@ -90,14 +90,14 @@ static void set_authnamecode(OGRSpatialReferenceH);
 
 /**
  * \brief Read projection information in WKT format from stdin or a file
- * 
- * Reads projection information from a file or stdin and stores in global 
+ *
+ * Reads projection information from a file or stdin and stores in global
  * structs projinfo and projunits.
  * Populates global cellhd with default region information.
- * 
+ *
  * \param wktfile File to read WKT co-ordinate system description from; -
  *                for stdin
- * 
+ *
  * \return        2 if a projected or lat/long co-ordinate system has been
  *                defined; 1 if an unreferenced XY co-ordinate system has
  *                been defined
@@ -108,7 +108,7 @@ int input_wkt(char *wktfile)
     FILE *infd;
     char buff[8000], *tmpwkt;
     OGRSpatialReferenceH hSRS;
-    const char *papszOptions[3];
+    char *papszOptions[3];
     int ret;
 
     if (strcmp(wktfile, "-") == 0)
@@ -177,11 +177,13 @@ int input_wkt(char *wktfile)
     /* find authname and authcode */
     hSRS = OSRNewSpatialReference(buff);
     /* get clean WKT definition */
+#if GDAL_VERSION_MAJOR >= 3
     papszOptions[0] = G_store("MULTILINE=YES");
     papszOptions[1] = G_store("FORMAT=WKT2");
     papszOptions[2] = NULL;
-#if GDAL_VERSION_MAJOR >= 3
-    OSRExportToWktEx(hSRS, &tmpwkt, papszOptions);
+    OSRExportToWktEx(hSRS, &tmpwkt, (const char **)papszOptions);
+    G_free(papszOptions[0]);
+    G_free(papszOptions[1]);
 #else
     OSRExportToPrettyWkt(hSRS, &tmpwkt, FALSE);
 #endif
@@ -194,16 +196,16 @@ int input_wkt(char *wktfile)
 }
 
 /**
- * \brief Read projection information in PROJ.4 format from a string 
+ * \brief Read projection information in PROJ.4 format from a string
  *        or stdin
- * 
- * Reads projection information from a string or stdin and stores in global 
+ *
+ * Reads projection information from a string or stdin and stores in global
  * structs projinfo and projunits.
  * Populates global cellhd with default region information.
- * 
- * \param proj4params String representation of PROJ.4 co-ordinate system 
+ *
+ * \param proj4params String representation of PROJ.4 co-ordinate system
  *                    description, or - to read it from stdin
- * 
+ *
  * \return        2 if a projected or lat/long co-ordinate system has been
  *                defined; 1 if an unreferenced XY co-ordinate system has
  *                been defined
@@ -278,7 +280,7 @@ int input_srid(char *srid)
 #if PROJ_VERSION_MAJOR >= 6
     OGRSpatialReferenceH hSRS;
     int ret = 0;
-    const char *papszOptions[3];
+    char *papszOptions[3];
     PJ *obj;
     const char *tmpwkt;
 
@@ -287,7 +289,7 @@ int input_srid(char *srid)
     if (!obj)
 	G_fatal_error(_("SRID <%s> not recognized by PROJ"), srid);
 
-    tmpwkt = proj_as_wkt(NULL, obj, PJ_WKT2_2019, NULL);
+    tmpwkt = proj_as_wkt(NULL, obj, PJ_WKT2_LATEST, NULL);
     hSRS = OSRNewSpatialReference(tmpwkt);
     if (!hSRS)
 	G_fatal_error(_("WKT for SRID <%s> not recognized by GDAL"), srid);
@@ -298,9 +300,10 @@ int input_srid(char *srid)
     papszOptions[0] = G_store("MULTILINE=YES");
     papszOptions[1] = G_store("FORMAT=WKT2");
     papszOptions[2] = NULL;
-    if (OSRExportToWktEx(hSRS, &projwkt, papszOptions) != OGRERR_NONE)
+    if (OSRExportToWktEx(hSRS, &projwkt, (const char **)papszOptions) != OGRERR_NONE)
 	G_warning(_("Unable to convert srid to WKT"));
-
+    G_free(papszOptions[0]);
+    G_free(papszOptions[1]);
     /* GRASS proj info + units */
     ret = GPJ_osr_to_grass(&cellhd, &projinfo, &projunits, hSRS, 0);
 
@@ -317,15 +320,15 @@ int input_srid(char *srid)
 }
 
 /**
- * \brief Read projection information corresponding to an EPSG co-ordinate 
+ * \brief Read projection information corresponding to an EPSG co-ordinate
  *        system number
- * 
- * Determines projection information corresponding to an EPSG co-ordinate 
+ *
+ * Determines projection information corresponding to an EPSG co-ordinate
  * system number and stores in global structs projinfo and projunits.
  * Populates global cellhd with default region information.
- * 
+ *
  * \param epsg_num    EPSG number for co-ordinate system
- * 
+ *
  * \return        2 if a projected or lat/long co-ordinate system has been
  *                defined; 1 if an unreferenced XY co-ordinate system has
  *                been defined
@@ -336,7 +339,7 @@ int input_epsg(int epsg_num)
     OGRSpatialReferenceH hSRS;
     char epsgstr[100];
     int ret = 0;
-    const char *papszOptions[3];
+    char *papszOptions[3];
 
     /* Set finder function for locating OGR csv co-ordinate system tables */
     /* SetCSVFilenameHook(GPJ_set_csv_loc); */
@@ -360,8 +363,11 @@ int input_epsg(int epsg_num)
     papszOptions[0] = G_store("MULTILINE=YES");
     papszOptions[1] = G_store("FORMAT=WKT2");
     papszOptions[2] = NULL;
-    if (OSRExportToWktEx(hSRS, &projwkt, papszOptions) != OGRERR_NONE)
+    if (OSRExportToWktEx(hSRS, &projwkt, (const char **)papszOptions) != OGRERR_NONE)
 	G_warning(_("Unable to convert EPSG code to WKT"));
+
+    G_free(papszOptions[0]);
+    G_free(papszOptions[1]);
 #endif
 
     OSRDestroySpatialReference(hSRS);
@@ -372,19 +378,19 @@ int input_epsg(int epsg_num)
 }
 
 /**
- * \brief Read projection and region information associated with a 
+ * \brief Read projection and region information associated with a
  *        georeferenced file
- * 
+ *
  * Reads projection information associated with a georeferenced file, if
  * available. Attempts are made to open the file with OGR and GDAL in turn.
  * (GDAL conventionally supports raster formats, and OGR vector formats.)
- * 
+ *
  * If successful, projection and region information are read from the file
  * using GDAL/OGR functions and stored in global structs cellhd, projinfo
  * and projunits.
- * 
+ *
  * \param geofile Path to georeferenced file
- * 
+ *
  * \return        2 if a projected or lat/long co-ordinate system has been
  *                defined; 1 if an unreferenced XY co-ordinate system has
  *                been defined
@@ -392,6 +398,93 @@ int input_epsg(int epsg_num)
 
 int input_georef(char *geofile)
 {
+/* GDAL >= 3 */
+#if GDAL_VERSION_MAJOR >= 3
+    GDALDatasetH hDS = NULL;
+    OGRSpatialReferenceH hSRS = NULL;
+    int ret = 0;
+
+    GDALAllRegister();
+
+    /* Try opening file as vector first because it doesn't output a
+     * (potentially confusing) error message if it can't open the file */
+
+    /* Try opening as vector */
+    G_debug(1, "Trying to open <%s> as vector...", geofile);
+    if ((hDS = GDALOpenEx(geofile, GDAL_OF_VECTOR, NULL, NULL, NULL))
+        && GDALDatasetGetLayerCount(hDS) > 0) {
+
+	OGRLayerH ogr_layer;
+
+	ogr_layer = GDALDatasetGetLayer(hDS, 0);
+	hSRS = OGR_L_GetSpatialRef(ogr_layer);
+
+	if (hSRS)
+	    set_default_region();
+    }
+    else {
+	/* Try opening as raster */
+	G_debug(1, "Trying to open <%s> as raster...", geofile);
+
+	if ((hDS = GDALOpen(geofile, GA_ReadOnly))) {
+	    char **sds;
+
+	    /* does the dataset include subdatasets? */
+	    sds = GDALGetMetadata(hDS, "SUBDATASETS");
+	    if (sds && *sds) {
+		G_warning(_("Input dataset <%s> contains subdatasets. "
+		            "Please select a subdataset."), geofile);
+	    }
+	    else {
+		hSRS = GDALGetSpatialRef(hDS);
+
+		if (hSRS)
+		    set_gdal_region(hDS);
+	    }
+	}
+	else {
+	    int namelen;
+
+	    namelen = strlen(geofile);
+	    if (namelen > 4 && G_strcasecmp(geofile + (namelen - 4), ".prj") == 0) {
+		G_warning(_("<%s> is not a GDAL dataset, trying to open it as ESRI WKT"),
+			  geofile);
+
+		return input_wkt(geofile);
+	    }
+	    else {
+		G_fatal_error(_("Unable to read georeferenced file <%s> using "
+				"GDAL library"), geofile);
+	    }
+	}
+    }
+    if (hSRS) {
+	char **papszOptions = NULL;
+
+	ret = GPJ_osr_to_grass(&cellhd, &projinfo, &projunits, hSRS, 0);
+
+	if (cellhd.proj == PROJECTION_XY)
+	    G_warning(_("Read of file %s was successful, but it did not contain "
+			"projection information. 'XY (unprojected)' will be used"),
+		      geofile);
+
+	papszOptions = G_calloc(3, sizeof(char *));
+	papszOptions[0] = G_store("MULTILINE=YES");
+	papszOptions[1] = G_store("FORMAT=WKT2");
+	OSRExportToWktEx(hSRS, &projwkt, (const char **)papszOptions);
+	G_free(papszOptions[0]);
+	G_free(papszOptions[1]);
+	G_free(papszOptions);
+
+	set_authnamecode(hSRS);
+    }
+    if (hDS)
+	GDALClose(hDS);
+
+    return ret;
+
+/* GDAL < 3 */
+#else
     OGRDataSourceH ogr_ds;
     OGRSpatialReferenceH hSRS;
     int ret = 0;
@@ -407,20 +500,13 @@ int input_georef(char *geofile)
     if ((ogr_ds = OGROpen(geofile, FALSE, NULL))
 	&& (OGR_DS_GetLayerCount(ogr_ds) > 0)) {
 	OGRLayerH ogr_layer;
-	const char **papszOptions = NULL;
 
 	G_debug(1, "...succeeded.");
 	/* Get the first layer */
 	ogr_layer = OGR_DS_GetLayer(ogr_ds, 0);
 	hSRS = OGR_L_GetSpatialRef(ogr_layer);
 	ret = GPJ_osr_to_grass(&cellhd, &projinfo, &projunits, hSRS, 0);
-#if GDAL_VERSION_MAJOR >= 3
-	papszOptions = G_calloc(2, sizeof(char *));
-	papszOptions[0] = G_store("FORMAT=WKT2");
-	OSRExportToWktEx(hSRS, &projwkt, papszOptions);
-#else
 	OSRExportToWkt(hSRS, &projwkt);
-#endif
 	set_default_region();
     }
     else {
@@ -434,6 +520,7 @@ int input_georef(char *geofile)
 	    char *wktstring;
 
 	    G_debug(1, "...succeeded.");
+	    /* TODO: change for GDAL 3+ */
 	    wktstring = (char *)GDALGetProjectionRef(gdal_ds);
 	    projwkt = G_store(wktstring);
 	    ret =
@@ -441,11 +528,7 @@ int input_georef(char *geofile)
 				 0);
 
 	    set_gdal_region(gdal_ds);
-#if GDAL_VERSION_MAJOR >= 3
-	    hSRS = OSRClone(GDALGetSpatialRef(gdal_ds));
-#else
 	    hSRS = OSRNewSpatialReference(projwkt);
-#endif
 	    GDALClose(gdal_ds);
 	}
 	else {
@@ -478,12 +561,13 @@ int input_georef(char *geofile)
 	OSRDestroySpatialReference(hSRS);
 
     return ret;
+#endif
 }
 
 /**
- * \brief Populates global cellhd with region settings based on 
+ * \brief Populates global cellhd with region settings based on
  *        georeferencing information in a GDAL dataset
- * 
+ *
  * \param hDS GDAL dataset object to retrieve georeferencing information from
  **/
 
@@ -548,7 +632,7 @@ void set_authnamecode(OGRSpatialReferenceH hSRS)
 	    authcode = OSRGetAuthorityCode(hSRS, authkey);
 	    if (authcode && *authcode) {
 		G_asprintf(&projsrid, "%s:%s", authname, authcode);
-		/* for backwards compatibility, remove ? */
+		/* for backwards compatibility; remove ? */
 		if (strcmp(authname, "EPSG") == 0) {
 		    projepsg = G_create_key_value();
 		    G_set_key_value("epsg", authcode, projepsg);

general/g.proj/output.c

@@ -1,7 +1,7 @@
-/*  
+/*
  ****************************************************************************
  *
- * MODULE:       g.proj 
+ * MODULE:       g.proj
  * AUTHOR(S):    Paul Kelly - paul-grass@stjohnspoint.co.uk
  * PURPOSE:      Provides a means of reporting the contents of GRASS
  *               projection information files and creating
@@ -48,21 +48,18 @@ void print_projinfo(int shell)
 	    fprintf(stdout, "%-11s: %s\n", projinfo->key[i], projinfo->value[i]);
     }
 
-    if (projepsg) {
-        const char *epsg_value, *epsg_key;
-
-	epsg_key = projepsg->key[0];
-	epsg_value = projepsg->value[0];
+    /* TODO: use projsrid instead */
+    if (projsrid) {
 
 	if (!shell) {
 	    fprintf(stdout,
-		"-PROJ_EPSG-------------------------------------------------\n");
-	    fprintf(stdout, "%-11s: %s\n", epsg_key, epsg_value);
+		"-PROJ_SRID-------------------------------------------------\n");
+	    fprintf(stdout, "%-11s: %s\n", "SRID", projsrid);
 	}
 	else
-	    fprintf(stdout, "%s=%s\n", epsg_key, epsg_value);
+	    fprintf(stdout, "%s=%s\n", "srid", projsrid);
     }
- 
+
     if (projunits) {
 	if (!shell)
 	    fprintf(stdout,
@@ -76,10 +73,11 @@ void print_projinfo(int shell)
 			projunits->key[i], projunits->value[i]);
 	}
     }
-    
+
     return;
 }
 
+/* DEPRECATED: datum transformation is handled by PROJ */
 void print_datuminfo(void)
 {
     char *datum, *params;
@@ -212,7 +210,7 @@ void print_wkt(int esristyle, int dontprettify)
     {
 	OGRSpatialReferenceH hSRS;
 	const char *tmpwkt;
-	const char **papszOptions = NULL;
+	char **papszOptions = NULL;
 
 	papszOptions = G_calloc(3, sizeof(char *));
 	if (dontprettify)
@@ -233,13 +231,13 @@ void print_wkt(int esristyle, int dontprettify)
 	    obj = proj_create(NULL, projsrid);
 	    if (!obj)
 		G_fatal_error(_("Unable to create PROJ definition from srid <%s>"), projsrid);
-	    tmpwkt = proj_as_wkt(NULL, obj, PJ_WKT2_2019, NULL);
+	    tmpwkt = proj_as_wkt(NULL, obj, PJ_WKT2_LATEST, NULL);
 	    hSRS = OSRNewSpatialReference(tmpwkt);
-	    OSRExportToWktEx(hSRS, &outwkt, papszOptions);
+	    OSRExportToWktEx(hSRS, &outwkt, (const char **)papszOptions);
 	}
 	if (!outwkt && projwkt) {
 	    hSRS = OSRNewSpatialReference(projwkt);
-	    OSRExportToWktEx(hSRS, &outwkt, papszOptions);
+	    OSRExportToWktEx(hSRS, &outwkt, (const char **)papszOptions);
 	}
 	if (!outwkt && projepsg) {
 	    int epsg_num;
@@ -248,15 +246,18 @@ void print_wkt(int esristyle, int dontprettify)
 
 	    hSRS = OSRNewSpatialReference(NULL);
 	    OSRImportFromEPSG(hSRS, epsg_num);
-	    OSRExportToWktEx(hSRS, &outwkt, papszOptions);
+	    OSRExportToWktEx(hSRS, &outwkt, (const char **)papszOptions);
 	}
 	if (!outwkt) {
 	    /* use GRASS proj info + units */
 	    projwkt = GPJ_grass_to_wkt2(projinfo, projunits, projepsg, esristyle,
 				      !(dontprettify));
 	    hSRS = OSRNewSpatialReference(projwkt);
-	    OSRExportToWktEx(hSRS, &outwkt, papszOptions);
+	    OSRExportToWktEx(hSRS, &outwkt, (const char **)papszOptions);
 	}
+	G_free(papszOptions[0]);
+	G_free(papszOptions[1]);
+	G_free(papszOptions);
 	if (hSRS)
 	    OSRDestroySpatialReference(hSRS);
     }

general/g.region/printwindow.c

@@ -264,7 +264,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 	    if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
 		G_fatal_error(_("Unable to initialize coordinate transformation"));
 
-	    /* 
+	    /*
 	     *  1 ------ 2
 	     *  |        |  map corners
 	     *  |        |
@@ -276,7 +276,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 	    /* get lat/long w/ same datum/ellipsoid as input */
 	    if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 			      &longitude, &latitude, NULL) < 0)
-		G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 		               "GPJ_transform()");
 
 	    lo1 = longitude;
@@ -287,7 +287,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 	    /* get lat/long w/ same datum/ellipsoid as input */
 	    if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 			      &longitude, &latitude, NULL) < 0)
-		G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 		               "GPJ_transform()");
 
 	    lo2 = longitude;
@@ -298,7 +298,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 	    /* get lat/long w/ same datum/ellipsoid as input */
 	    if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 			      &longitude, &latitude, NULL) < 0)
-		G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 		               "GPJ_transform()");
 
 	    lo3 = longitude;
@@ -309,7 +309,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 	    /* get lat/long w/ same datum/ellipsoid as input */
 	    if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 			      &longitude, &latitude, NULL) < 0)
-		G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 		               "GPJ_transform()");
 
 	    lo4 = longitude;
@@ -322,7 +322,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 	    /* get lat/long w/ same datum/ellipsoid as input */
 	    if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 			      &longitude, &latitude, NULL) < 0)
-		G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 		               "GPJ_transform()");
 
 	    loc = longitude;
@@ -507,7 +507,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 	    /* get lat/long w/ same datum/ellipsoid as input */
 	    if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 			      &longitude, &latitude, NULL) < 0)
-		G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 		               "GPJ_transform()");
 
 #ifdef HAVE_PROJ_H
@@ -534,7 +534,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
     }
 
     /* flag.bbox
-       Calculate the largest bounding box in lat-lon coordinates 
+       Calculate the largest bounding box in lat-lon coordinates
        and print it to stdout
      */
     if (print_flag & PRINT_MBBOX) {
@@ -556,7 +556,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 		G_fatal_error(_("Can't get projection info of current location"));
 	    /* do not wrap to -180, 180, otherwise east can be < west */
 	    /* TODO: for PROJ 6+, the +over switch must be added to the
-	     * transformation pipeline if authority:name or WKt are used 
+	     * transformation pipeline if authority:name or WKt are used
 	     * as crs definition */
 	    G_set_key_value("over", "defined", in_proj_info);
 
@@ -605,7 +605,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 	    longitude = (window->west + window->east) / 2.;
 	    if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 			      &longitude, &latitude, NULL) < 0)
-		G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 		               "GPJ_transform()");
 
 	    sh_ll_w = sh_ll_e = longitude;
@@ -619,7 +619,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 		longitude = window->west;
 		if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 				  &longitude, &latitude, NULL) < 0)
-		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 				   "GPJ_transform()");
 
 		if (sh_ll_n < latitude)
@@ -638,7 +638,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 		longitude = window->east;
 		if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 				  &longitude, &latitude, NULL) < 0)
-		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 				   "GPJ_transform()");
 
 		if (sh_ll_n < latitude)
@@ -658,7 +658,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 		longitude = window->west + c * window->ew_res;
 		if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 				  &longitude, &latitude, NULL) < 0)
-		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 				   "GPJ_transform()");
 
 		if (sh_ll_n < latitude)
@@ -675,7 +675,7 @@ void print_window(struct Cell_head *window, int print_flag, int flat_flag)
 		longitude = window->west + c * window->ew_res;
 		if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
 				  &longitude, &latitude, NULL) < 0)
-		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 				   "GPJ_transform()");
 
 		if (sh_ll_n < latitude)

gui/wxpython/core/render.py

@@ -843,7 +843,7 @@ class Map(object):
 
         for line in ret.splitlines():
             if ':' in line:
-                key, val = map(lambda x: x.strip(), line.split(':'))
+                key, val = map(lambda x: x.strip(), line.split(':', 1))
                 if key in ['units']:
                     val = val.lower()
                 projinfo[key] = val

include/defs/gis.h

@@ -41,9 +41,9 @@
 #   include <sys/param.h>
 #  endif
 #  if (defined(BSD))
-    /* no malloc.h, no alloca.h ? 
-     * TODO: better 
-     * check if alloca.h exists, 
+    /* no malloc.h, no alloca.h ?
+     * TODO: better
+     * check if alloca.h exists,
      * if not, check if malloc.h exists,
      * if not use stdlib.h */
 #   include <stdlib.h>
@@ -72,11 +72,11 @@
 #define RELDIR "?"
 #endif
 
-/* GDAL < 2.3 does not define HAVE_LONG_LONG when compiled with 
+/* GDAL < 2.3 does not define HAVE_LONG_LONG when compiled with
  * Visual Studio as for OSGeo4W, even though long long is available,
  * and GIntBig falls back to long which is on Windows always 4 bytes.
  * This patch ensures that GIntBig is defined as long long (8 bytes)
- * if GDAL is compiled with Visual Studio and GRASS is compiled with 
+ * if GDAL is compiled with Visual Studio and GRASS is compiled with
  * MinGW. This patch must be applied before other GDAL/OGR headers are
  * included, as done by gprojects.h and vector.h */
 #if defined(__MINGW32__) && HAVE_GDAL
@@ -149,10 +149,10 @@ void G_ascii_check(char *);
 /* Do it better if you know how */
 /* asprintf is not found on MINGW but exists */
 
-/* 
+/*
  *  Because configure script in GDAL test is G_asprintf exists in gis lib
  *  the G_asprintf macro is disabled until a stable version of GDAL
- *  with a different function becomes widely used 
+ *  with a different function becomes widely used
  */
 int G_vasprintf(char **, const char *, va_list);
 int G_asprintf(char **, const char *, ...)
@@ -510,7 +510,7 @@ int G_make_location(const char *, struct Cell_head *, const struct Key_Value *,
 		    const struct Key_Value *);
 int G_make_location_epsg(const char *, struct Cell_head *, const struct Key_Value *,
 			 const struct Key_Value *, const struct Key_Value *);
-int G_make_location_crs(const char *, struct Cell_head *, const struct Key_Value *,
+int G_make_location_crs(const char *, struct Cell_head *,
 			const struct Key_Value *, const struct Key_Value *,
 			const char *, const char *);
 int G_write_projsrid(const char *, const char *);

include/gprojects.h

@@ -23,6 +23,27 @@
 #include <proj.h>
 #define RAD_TO_DEG    57.295779513082321
 #define DEG_TO_RAD   .017453292519943296
+
+/* adapted from gdal_version.h */
+#ifndef PROJ_COMPUTE_VERSION
+#define PROJ_COMPUTE_VERSION(maj,min,rev) ((maj)*1000000+(min)*10000+(rev)*100)
+#endif
+
+/* just in case PROJ introduces PROJ_VERSION_NUM in a future version */
+#ifdef PROJ_VERSION_NUM
+#undef PROJ_VERSION_NUM
+#endif
+#define PROJ_VERSION_NUM      (PROJ_COMPUTE_VERSION(PROJ_VERSION_MAJOR,PROJ_VERSION_MINOR,PROJ_VERSION_PATCH))
+
+#ifndef PJ_WKT2_LATEST
+/* update if new PROJ versions support new WKT2 standards */
+#if PROJ_VERSION_NUM < 6030000
+#define PJ_WKT2_LATEST PJ_WKT2_2018
+#else
+#define PJ_WKT2_LATEST PJ_WKT2_2019
+#endif
+#endif
+
 #else
 #include <proj_api.h>
 #define PJ_FWD 	 1
@@ -57,6 +78,7 @@ struct pj_info
     char proj[100];
     char *def;
     char *srid;
+    char *wkt;
 };
 
 struct gpj_datum

lib/gis/adj_cellhd.c

@@ -34,7 +34,7 @@ static int ll_check_ew(struct Cell_head *cellhd);
  * region). It also makes projection-specific adjustments. The
  * <i>cellhd</i> structure must have its <i>north, south, east,
  * west</i>, and <i>proj</i> fields set.
- * 
+ *
  * If <i>row_flag</i> is true, then the north-south resolution is
  * computed from the number of <i>rows</i> in the <i>cellhd</i>
  * structure. Otherwise the number of <i>rows</i> is computed from the
@@ -127,16 +127,16 @@ void G_adjust_Cell_head(struct Cell_head *cellhd, int row_flag, int col_flag)
  * region).  It also makes projection-specific adjustments. The
  * <i>cellhd</i> structure must have its <i>north, south, east,
  * west</i>, and <i>proj</i> fields set.
- * 
- * If <i>row_flag</i> is true, then the north-south resolution is computed 
- * from the number of <i>rows</i> in the <i>cellhd</i> structure. 
- * Otherwise the number of <i>rows</i> is computed from the north-south 
- * resolution in the structure, similarly for <i>col_flag</i> and the 
- * number of columns and the east-west resolution. 
  *
- * If <i>depth_flag</i> is true, top-bottom resolution is calculated 
+ * If <i>row_flag</i> is true, then the north-south resolution is computed
+ * from the number of <i>rows</i> in the <i>cellhd</i> structure.
+ * Otherwise the number of <i>rows</i> is computed from the north-south
+ * resolution in the structure, similarly for <i>col_flag</i> and the
+ * number of columns and the east-west resolution.
+ *
+ * If <i>depth_flag</i> is true, top-bottom resolution is calculated
  * from depths.
- * If <i>depth_flag</i> are false, number of depths is calculated from 
+ * If <i>depth_flag</i> are false, number of depths is calculated from
  * top-bottom resolution.
  *
  * \warning This function can cause segmentation fault without any warning
@@ -341,7 +341,7 @@ static int ll_check_ns(struct Cell_head *cellhd)
 
     G_debug(3, "ll_check_ns: epsilon: %g", llepsilon);
 
-    /* North, South: allow a half cell spill-over */ 
+    /* North, South: allow a half cell spill-over */
 
     diff = (cellhd->north - cellhd->south) / cellhd->ns_res;
     ncells = (int) (diff + 0.5);
@@ -372,7 +372,7 @@ static int ll_check_ns(struct Cell_head *cellhd)
 	if (diff <= 0.5 + llepsilon) {
 	    G_important_message(_("90 degree north is exceeded by %g cells"),
 		      diff);
-	    
+
 	    if (diff < llepsilon && diff > fpepsilon) {
 		G_verbose_message(_("Subtle input data rounding error of north boundary (%g)"),
 			  cellhd->north - 90.0);
@@ -422,7 +422,7 @@ static int ll_check_ns(struct Cell_head *cellhd)
 	if (diff <= 0.5 + llepsilon) {
 	    G_important_message(_("90 degree south is exceeded by %g cells"),
 		      diff);
-	    
+
 	    if (diff < llepsilon && diff > fpepsilon) {
 		G_verbose_message(_("Subtle input data rounding error of south boundary (%g)"),
 			  cellhd->south + 90);
@@ -451,7 +451,7 @@ static int ll_check_ns(struct Cell_head *cellhd)
 	else
 	    G_fatal_error(_("Illegal latitude for South"));
     }
-    
+
     if (lladjust)
 	cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
 
@@ -500,7 +500,7 @@ static int ll_check_ew(struct Cell_head *cellhd)
 /*!
  * \brief Adjust window for lat/lon.
  *
- * This function tries to automatically fix fp precision issues and 
+ * This function tries to automatically fix fp precision issues and
  * adjust rounding errors for lat/lon.
  *
  * <b>Note:</b> 3D values are not adjusted.

lib/gis/get_projinfo.c

@@ -1,10 +1,10 @@
 /*!
   \file lib/gis/get_projinfo.c
-  
+
   \brief GIS Library - Get projection info
-  
+
   (C) 1999-2014 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.
 */
@@ -23,7 +23,7 @@
 
   Note: Allocated Key_Value structure should be freed by
   G_free_key_value().
-  
+
   Prints a warning if no units information available.
 
   \return pointer to Key_Value structure with key/value pairs
@@ -49,10 +49,10 @@ struct Key_Value *G_get_projunits(void)
 
 /*!
   \brief Gets projection information for location
-  
+
   Note: Allocated Key_Value structure should be freed by
   G_free_key_value().
-  
+
   Prints a warning if no projection information available.
 
   \return pointer to Key_Value structure with key/value pairs
@@ -72,12 +72,12 @@ struct Key_Value *G_get_projinfo(void)
 	return NULL;
     }
     in_proj_keys = G_read_key_value_file(path);
-    
+
     /* TODO: do not restrict to EPSG as the only authority */
     if ((in_epsg_keys = G_get_projepsg()) != NULL) {
 	const char *epsgstr = G_find_key_value("epsg", in_epsg_keys);
 	char buf[4096];
-	
+
 	sprintf(buf, "EPSG:%s", epsgstr);
 	G_set_key_value("init", buf, in_proj_keys);
 	G_free_key_value(in_epsg_keys);
@@ -88,15 +88,17 @@ struct Key_Value *G_get_projinfo(void)
 
 /*!
   \brief Gets EPSG information for the current location
-  
+
+  DEPRECATED: Use G_get_projsrid() instead.
+
   Note: Allocated Key_Value structure should be freed by
   G_free_key_value().
-  
+
   \return pointer to Key_Value structure with key/value pairs
   \return NULL when EPSG code is not defined for location
 */
 
-/* superseded by G_get_projsrid(), keep for backwards compatibility */ 
+/* superseded by G_get_projsrid(), keep for backwards compatibility */
 struct Key_Value *G_get_projepsg(void)
 {
     struct Key_Value *in_epsg_keys;
@@ -232,8 +234,24 @@ char *G_get_projsrid(void)
     G_file_name(path, "", SRID_FILE, "PERMANENT");
     if (access(path, 0) != 0) {
 	if (G_projection() != PROJECTION_XY) {
+	    struct Key_Value *projepsg;
+	    const char *epsg_num;
+
 	    G_debug(1, "<%s> file not found for location <%s>",
 		      SRID_FILE, G_location());
+
+	    /* for backwards compatibility, check if PROJ_EPSG exists */
+	    if ((projepsg = G_get_projepsg()) != NULL) {
+		epsg_num = G_find_key_value("epsg", projepsg);
+		if (*epsg_num) {
+		    G_debug(1, "Using <%s> file instead for location <%s>",
+			    EPSG_FILE, G_location());
+		    G_asprintf(&sridstring, "EPSG:%s", epsg_num);
+		    G_free_key_value(projepsg);
+
+		    return sridstring;
+		}
+	    }
 	}
 	return NULL;
     }

lib/gis/make_loc.c

@@ -26,7 +26,7 @@
 
 /*!
  * \brief Create a new location
- * 
+ *
  * This function creates a new location in the current database,
  * initializes the projection, default window and current window.
  *
@@ -48,7 +48,7 @@
  * \return 0 on success
  * \return -1 to indicate a system error (check errno).
  * \return -2 failed to create projection file (currently not used)
- * \return -3 illegal name 
+ * \return -3 illegal name
  */
 int G_make_location(const char *location_name,
                     struct Cell_head *wind,
@@ -96,7 +96,7 @@ int G_make_location(const char *location_name,
 
 /*!
  * \brief Create a new location
- * 
+ *
  * This function creates a new location in the current database,
  * initializes the projection, default window and current window,
  * and sets the EPSG code if present
@@ -115,14 +115,14 @@ int G_make_location(const char *location_name,
  *
  * \param proj_units    projection units suitable to write to the PROJ_UNITS
  *                      file, or NULL.
- * 
+ *
  * \param proj_epsg     EPSG code suitable to write to the PROJ_EPSG
  *                      file, or NULL.
  *
  * \return 0 on success
  * \return -1 to indicate a system error (check errno).
  * \return -2 failed to create projection file (currently not used)
- * \return -3 illegal name 
+ * \return -3 illegal name
  */
 int G_make_location_epsg(const char *location_name,
 			 struct Cell_head *wind,
@@ -149,7 +149,7 @@ int G_make_location_epsg(const char *location_name,
 
 /*!
  * \brief Create a new location
- * 
+ *
  * This function creates a new location in the current database,
  * initializes the projection, default window and current window,
  * and sets WKT, srid, and EPSG code if present
@@ -168,7 +168,7 @@ int G_make_location_epsg(const char *location_name,
  *
  * \param proj_units    projection units suitable to write to the PROJ_UNITS
  *                      file, or NULL.
- * 
+ *
  * \param proj_epsg     EPSG code suitable to write to the PROJ_EPSG
  *                      file, or NULL.
  *
@@ -181,15 +181,14 @@ int G_make_location_epsg(const char *location_name,
  * \return 0 on success
  * \return -1 to indicate a system error (check errno).
  * \return -2 failed to create projection file (currently not used)
- * \return -3 illegal name 
+ * \return -3 illegal name
  */
 int G_make_location_crs(const char *location_name,
 			struct Cell_head *wind,
 			const struct Key_Value *proj_info,
 			const struct Key_Value *proj_units,
-			const struct Key_Value *proj_epsg,
-			const char *proj_wkt,
-			const char *proj_srid)
+			const char *proj_srid,
+			const char *proj_wkt)
 {
     int ret;
     char path[GPATH_MAX];
@@ -199,10 +198,9 @@ int G_make_location_crs(const char *location_name,
     if (ret != 0)
 	return ret;
 
-    /* Write out PROJ_EPSG if epsg code is available. */
-    if (proj_epsg != NULL) {
-	G_file_name(path, "", "PROJ_EPSG", "PERMANENT");
-	G_write_key_value_file(path, proj_epsg);
+    /* Write out PROJ_SRID if srid is available. */
+    if (proj_srid != NULL) {
+	G_write_projsrid(location_name, proj_srid);
     }
 
     /* Write out PROJ_WKT if WKT is available. */
@@ -210,11 +208,6 @@ int G_make_location_crs(const char *location_name,
 	G_write_projwkt(location_name, proj_wkt);
     }
 
-    /* Write out PROJ_SRID if srid is available. */
-    if (proj_srid != NULL) {
-	G_write_projsrid(location_name, proj_srid);
-    }
-
     return 0;
 }
 
@@ -469,7 +462,7 @@ int G_compare_projections(const struct Key_Value *proj_info1,
 	    /* lat_1 differ */
 	    /* check for swapped lat_1, lat_2 */
 	    l_2 = G_find_key_value("lat_2", proj_info2);
-	    
+
 	    if (!l_2)
 		return -11;
 	    if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
@@ -488,7 +481,7 @@ int G_compare_projections(const struct Key_Value *proj_info1,
 	    /* lat_2 differ */
 	    /* check for swapped lat_1, lat_2 */
 	    l_2 = G_find_key_value("lat_1", proj_info2);
-	    
+
 	    if (!l_2)
 		return -11;
 	    if (l_1 && l_2 && (fabs(atof(l_1) - atof(l_2)) > 0.000001)) {
@@ -508,7 +501,7 @@ int G_compare_projections(const struct Key_Value *proj_info1,
 
 /*!
    \brief Write WKT definition to file
-   
+
    Any WKT string and version recognized by PROJ is supported.
 
    \param location_name name of the location to write the WKT definition

lib/proj/get_proj.c

@@ -8,7 +8,7 @@
    Eric Miller, Paul Kelly, Markus Metz
 
    (C) 2003-2008, 2018 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.
@@ -38,18 +38,18 @@ static int nopt;
 /**
  * \brief Create a pj_info struct Co-ordinate System definition from a set of
  *        PROJ_INFO / PROJ_UNITS-style key-value pairs
- * 
+ *
  * This function takes a GRASS-style co-ordinate system definition as stored
  * in the PROJ_INFO and PROJ_UNITS files and processes it to create a pj_info
  * representation for use in re-projecting with pj_do_proj(). In addition to
- * the parameters passed to it it may also make reference to the system 
+ * the parameters passed to it it may also make reference to the system
  * ellipse.table and datum.table files if necessary.
- * 
- * \param info Pointer to a pj_info struct (which must already exist) into 
+ *
+ * \param info Pointer to a pj_info struct (which must already exist) into
  *        which the co-ordinate system definition will be placed
  * \param in_proj_keys PROJ_INFO-style key-value pairs
  * \param in_units_keys PROJ_UNITS-style key-value pairs
- * 
+ *
  * \return -1 on error (unable to initialise PROJ.4)
  *          2 if "default" 3-parameter datum shift values from datum.table
  *            were used
@@ -82,6 +82,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
     info->def = NULL;
     info->pj = NULL;
     info->srid = NULL;
+    info->wkt = NULL;
 
     str = G_find_key_value("meters", in_units_keys);
     if (str != NULL) {
@@ -120,7 +121,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
 	else if (strcmp(in_proj_keys->key[i], "zone") == 0) {
 	    continue;
 
-	    /* Datum and ellipsoid-related parameters will be handled 
+	    /* Datum and ellipsoid-related parameters will be handled
 	     * separately after end of this loop PK */
 
 	}
@@ -192,7 +193,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
     else {
 	sprintf(buffa, "a=%.16g", a);
 	alloc_options(buffa);
-	/* Cannot use es directly because the OSRImportFromProj4() 
+	/* Cannot use es directly because the OSRImportFromProj4()
 	 * function in OGR only accepts b or rf as the 2nd parameter */
 	if (es == 0)
 	    sprintf(buffa, "b=%.16g", a);
@@ -214,7 +215,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
 	alloc_options(buffa);
 	G_free(params);
 
-	/* else if a datum name is present take it and look up the parameters 
+	/* else if a datum name is present take it and look up the parameters
 	 * from the datum.table file */
     }
     else if (datum != NULL) {
@@ -225,7 +226,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
 	    returnval = 2;
 	    G_free(params);
 
-	    /* else just pass the datum name on and hope it is recognised by 
+	    /* else just pass the datum name on and hope it is recognised by
 	     * PROJ.4 even though it isn't recognised by GRASS */
 	}
 	else {
@@ -242,7 +243,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
 
 #ifdef HAVE_PROJ_H
 #if PROJ_VERSION_MAJOR >= 6
-    /* without type=crs, PROJ6 does not recognize what this is, 
+    /* without type=crs, PROJ6 does not recognize what this is,
      * a crs or some kind of coordinate operation, falling through to
      * PJ_TYPE_OTHER_COORDINATE_OPERATION */
     alloc_options("type=crs");
@@ -276,7 +277,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
     if (perr)
 	G_fatal_error("PROJ 5 error %d", perr);
 #endif
-    
+
     info->pj = pj;
 
     deflen = 0;
@@ -309,20 +310,20 @@ static void alloc_options(char *buffa)
 }
 
 /**
- * \brief Create a pj_info struct Co-ordinate System definition from a 
+ * \brief Create a pj_info struct Co-ordinate System definition from a
  *        string with a sequence of key=value pairs
- * 
- * This function takes a GRASS- or PROJ style co-ordinate system definition 
- * and processes it to create a pj_info representation for use in 
- * re-projecting with pj_do_proj(). In addition to the parameters passed 
- * to it it may also make reference to the system ellipse.table and 
+ *
+ * This function takes a GRASS- or PROJ style co-ordinate system definition
+ * and processes it to create a pj_info representation for use in
+ * re-projecting with pj_do_proj(). In addition to the parameters passed
+ * to it it may also make reference to the system ellipse.table and
  * datum.table files if necessary.
- * 
- * \param info Pointer to a pj_info struct (which must already exist) into 
+ *
+ * \param info Pointer to a pj_info struct (which must already exist) into
  *        which the co-ordinate system definition will be placed
  * \param str input string with projection definition
  * \param in_units_keys PROJ_UNITS-style key-value pairs
- * 
+ *
  * \return -1 on error (unable to initialise PROJ.4)
  *          1 on success
  **/
@@ -346,12 +347,12 @@ int pj_get_string(struct pj_info *info, char *str)
     info->def = NULL;
     info->srid = NULL;
     info->pj = NULL;
-    
+
     nopt = 0;
 
     if ((str == NULL) || (str[0] == '\0')) {
-	/* Null Pointer or empty string is supplied for parameters, 
-	 * implying latlong projection; just need to set proj 
+	/* Null Pointer or empty string is supplied for parameters,
+	 * implying latlong projection; just need to set proj
 	 * parameter and call pj_init PK */
 	sprintf(info->proj, "ll");
 	sprintf(buffa, "proj=latlong ellps=WGS84");
@@ -360,7 +361,7 @@ int pj_get_string(struct pj_info *info, char *str)
     else {
 	/* Parameters have been provided; parse through them but don't
 	 * bother with most of the checks in pj_get_kv; assume the
-	 * programmer knows what he / she is doing when using this 
+	 * programmer knows what he / she is doing when using this
 	 * function rather than reading a PROJ_INFO file       PK */
 	s = str;
 	while (s = strtok(s, " \t\n"), s) {
@@ -405,7 +406,7 @@ int pj_get_string(struct pj_info *info, char *str)
 
 #ifdef HAVE_PROJ_H
 #if PROJ_VERSION_MAJOR >= 6
-    /* without type=crs, PROJ6 does not recognize what this is, 
+    /* without type=crs, PROJ6 does not recognize what this is,
      * a crs or some kind of coordinate operation, falling through to
      * PJ_TYPE_OTHER_COORDINATE_OPERATION */
     alloc_options("type=crs");
@@ -454,16 +455,16 @@ int pj_get_string(struct pj_info *info, char *str)
 /**
  * \brief Define a latitude / longitude co-ordinate system with the same
  *        ellipsoid and datum parameters as an existing projected system
- * 
+ *
  * This function is useful when projected co-ordinates need to be simply
  * converted to and from latitude / longitude.
- * 
+ *
  * \param pjnew Pointer to pj_info struct for geographic co-ordinate system
  *        that will be created
  * \param pjold Pointer to pj_info struct for existing projected co-ordinate
  *        system
- * 
- * \return 1 on success; -1 if there was an error (i.e. if the PROJ.4 
+ *
+ * \return 1 on success; -1 if there was an error (i.e. if the PROJ.4
  *         pj_latlong_from_proj() function returned NULL)
  **/
 
@@ -520,11 +521,11 @@ const char *set_proj_share(const char *name)
 /**
  * \brief Print projection parameters as used by PROJ.4 for input and
  *        output co-ordinate systems
- * 
+ *
  * \param iproj 'Input' co-ordinate system
  * \param oproj 'Output' co-ordinate system
- * 
- * \return 1 on success, -1 on error (i.e. if the PROJ-style definition 
+ *
+ * \return 1 on success, -1 on error (i.e. if the PROJ-style definition
  *         is NULL for either co-ordinate system)
  **/
 

raster/r.external/proj.c

@@ -11,22 +11,72 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
 		      int check_only)
 {
     struct Cell_head loc_wind;
-    struct Key_Value *proj_info, *proj_units, *proj_epsg;
-    struct Key_Value *loc_proj_info, *loc_proj_units;
-    const char *wkt;
+    struct Key_Value *proj_info = NULL,
+                     *proj_units = NULL;
+    struct Key_Value *loc_proj_info = NULL,
+                     *loc_proj_units = NULL;
+    char *wkt = NULL, *srid = NULL;
     char error_msg[8096];
     int proj_trouble;
 
     /* -------------------------------------------------------------------- */
-    /*      Fetch the projection in GRASS form.                             */
+    /*      Fetch the projection in GRASS form, SRID, and WKT.              */
     /* -------------------------------------------------------------------- */
-    proj_info = NULL;
-    proj_units = NULL;
-    proj_epsg = NULL;
-    loc_proj_info = NULL;
-    loc_proj_units = NULL;
 
-    wkt = GDALGetProjectionRef(hDS);
+#if GDAL_VERSION_NUM >= 3000000
+    OGRSpatialReferenceH hSRS;
+
+    hSRS = GDALGetSpatialRef(hDS);
+    if (hSRS) {
+	/* get WKT2 definition */
+	char **papszOptions;
+
+	papszOptions = G_calloc(3, sizeof(char *));
+	papszOptions[0] = G_store("MULTILINE=YES");
+	papszOptions[1] = G_store("FORMAT=WKT2");
+	OSRExportToWktEx(hSRS, &wkt, (const char **)papszOptions);
+	G_free(papszOptions[0]);
+	G_free(papszOptions[1]);
+	G_free(papszOptions);
+    }
+    /* proj_trouble:
+     * 0: valid srs
+     * 1: no srs, default to xy
+     * 2: unreadable srs, default to xy
+     */
+
+    /* Projection only required for checking so convert non-interactively */
+    proj_trouble = 0;
+    if (wkt && *wkt) {
+	if (hSRS != NULL)
+	    GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0);
+
+	if (!hSRS || (!OSRIsProjected(hSRS) && !OSRIsGeographic(hSRS))) {
+	    G_important_message(_("Input contains an invalid SRS. "
+	                          "WKT definition:\n%s"), wkt);
+
+	    proj_trouble = 2;
+	}
+	else{
+	    const char *authkey, *authname, *authcode;
+
+	    if (OSRIsProjected(hSRS))
+		authkey = "PROJCS";
+	    else /* is geographic */
+		authkey = "GEOGCS";
+
+	    authname = OSRGetAuthorityName(hSRS, authkey);
+	    if (authname && *authname) {
+		authcode = OSRGetAuthorityCode(hSRS, authkey);
+		if (authcode && *authcode) {
+		    G_asprintf(&srid, "%s:%s", authname, authcode);
+		}
+	    }
+	}
+    }
+
+#else
+    wkt = G_store(GDALGetProjectionRef(hDS));
     /* proj_trouble:
      * 0: valid srs
      * 1: no srs, default to xy
@@ -43,7 +93,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
 	    GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0);
 
 	if (!hSRS || (!OSRIsProjected(hSRS) && !OSRIsGeographic(hSRS))) {
-	    G_important_message(_("Input contains an invalid SRS. " 
+	    G_important_message(_("Input contains an invalid SRS. "
 	                          "WKT definition:\n%s"), wkt);
 
 	    proj_trouble = 2;
@@ -57,18 +107,17 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
 		authkey = "GEOGCS";
 
 	    authname = OSRGetAuthorityName(hSRS, authkey);
-	    if (authname && *authname && strcmp(authname, "EPSG") == 0) {
+	    if (authname && *authname) {
 		authcode = OSRGetAuthorityCode(hSRS, authkey);
 		if (authcode && *authcode) {
-		    G_debug(1, "found EPSG:%s", authcode);
-		    proj_epsg = G_create_key_value();
-		    G_set_key_value("epsg", authcode, proj_epsg);
+		    G_asprintf(&srid, "%s:%s", authname, authcode);
 		}
 	    }
 	}
 	if (hSRS)
 	    OSRDestroySpatialReference(hSRS);
     }
+#endif
     else {
 	G_important_message(_("No projection information available"));
 	cellhd->proj = PROJECTION_XY;
@@ -80,14 +129,14 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
     /*      Do we need to create a new location?                            */
     /* -------------------------------------------------------------------- */
     if (outloc != NULL) {
-	/* do not create a xy location if an existing SRS was unreadable */ 
+	/* do not create a xy location if an existing SRS was unreadable */
 	if (proj_trouble == 2) {
 	    G_fatal_error(_("Unable to convert input map projection to GRASS "
 			    "format; cannot create new location."));
 	}
 	else {
-            if (0 != G_make_location_epsg(outloc, cellhd, proj_info,
-	                                  proj_units, proj_epsg)) {
+            if (0 != G_make_location_crs(outloc, cellhd, proj_info,
+	                                 proj_units, srid, wkt)) {
                 G_fatal_error(_("Unable to create new location <%s>"),
                               outloc);
             }
@@ -282,7 +331,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
                        _("Consider generating a new location from the input dataset using "
                          "the 'location' parameter.\n"));
             }
-            
+
 	    if (check_only)
 		msg_fn = G_message;
 	    else
@@ -297,7 +346,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
 	    if (check_only)
 		msg_fn = G_message;
 	    else
-		msg_fn = G_verbose_message;            
+		msg_fn = G_verbose_message;
 	    msg_fn(_("Projection of input dataset and current location "
 		     "appear to match"));
 	    if (check_only) {

raster/r.in.gdal/main.c

@@ -2,7 +2,7 @@
 /****************************************************************************
  *
  * MODULE:       r.in.gdal
- *               
+ *
  * AUTHOR(S):    Frank Warmerdam (copyright of this file)
  *               Added optional GCP transformation: Markus Neteler 10/2001
  *
@@ -63,7 +63,7 @@ static GDALDatasetH opends(char *dsname, const char **doo, GDALDriverH *hDriver)
     if (hDS == NULL)
         G_fatal_error(_("Unable to open datasource <%s>"), dsname);
     G_add_error_handler(error_handler_ds, hDS);
-    
+
     *hDriver = GDALGetDatasetDriver(hDS);	/* needed for AVHRR data */
     /* L1B - NOAA/AVHRR data must be treated differently */
     /* for hDriver names see gdal/frmts/gdalallregister.cpp */
@@ -134,7 +134,7 @@ int main(int argc, char *argv[])
     /* -------------------------------------------------------------------- */
     parm.input = G_define_standard_option(G_OPT_F_BIN_INPUT);
     parm.input->description = _("Name of raster file to be imported");
-    
+
     parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
 
     parm.band = G_define_option();
@@ -159,7 +159,7 @@ int main(int argc, char *argv[])
 	_("Name of location to create or to read projection from for GCPs transformation");
     parm.target->key_desc = "name";
     parm.target->guisection = _("Projection");
-    
+
     parm.title = G_define_option();
     parm.title->key = "title";
     parm.title->key_desc = "phrase";
@@ -228,7 +228,7 @@ int main(int argc, char *argv[])
     flag_o->description =
 	_("Assume that the dataset has same projection as the current location");
     flag_o->guisection = _("Projection");
- 
+
     flag_j = G_define_flag();
     flag_j->key = 'j';
     flag_j->description =
@@ -236,20 +236,20 @@ int main(int argc, char *argv[])
     flag_j->suppress_required = YES;
     flag_j->guisection = _("Projection");
     G_option_requires(flag_j, parm.input, NULL);
-    
+
     flag_e = G_define_flag();
     flag_e->key = 'e';
     flag_e->label = _("Extend region extents based on new dataset");
     flag_e->description =
 	_("Also updates the default region if in the PERMANENT mapset");
     flag_e->guisection = _("Region");
- 
+
     flag_f = G_define_flag();
     flag_f->key = 'f';
     flag_f->description = _("List supported formats and exit");
     flag_f->guisection = _("Print");
     flag_f->suppress_required = YES;
-    
+
     flag_l = G_define_flag();
     flag_l->key = 'l';
     flag_l->description =
@@ -287,7 +287,7 @@ int main(int argc, char *argv[])
     /* 1. list supported formats */
     /* 2. print input bands */
     /* 3. open input */
-    /* 4. check projection / create location */ 
+    /* 4. check projection / create location */
     /* 5. open output */
 
     /* The parser checks if the map already exists in current mapset, this is
@@ -345,7 +345,7 @@ int main(int argc, char *argv[])
     offset = atoi(parm.offset->answer);
 
     num_digits = atoi(parm.num_digits->answer);
-    
+
     if ((title = parm.title->answer))
 	G_strip(title);
 
@@ -379,7 +379,7 @@ int main(int argc, char *argv[])
     croptoregion = flag_r->answer;
     if (flag_r->answer && parm.outloc->answer) {
 	G_warning(_("Disabling '-r' flag for new location"));
-	croptoregion = 0; 
+	croptoregion = 0;
     }
 
     /* default GDAL memory cache size appears to be only 40 MiB, slowing down r.in.gdal */
@@ -436,7 +436,7 @@ int main(int argc, char *argv[])
 
     /* does the driver support subdatasets? */
     /* test for capability GDAL_DMD_SUBDATASETS */
-    
+
     /* does the dataset include subdatasets? */
     {
 	char **sds = GDALGetMetadata(hDS, "SUBDATASETS");
@@ -535,7 +535,7 @@ int main(int argc, char *argv[])
 
     /* zero cell header */
     G_zero(&cellhd, sizeof(struct Cell_head));
-    
+
     /* -------------------------------------------------------------------- */
     /*      Set up the window representing the data we have.                */
     /* -------------------------------------------------------------------- */
@@ -858,7 +858,7 @@ int main(int argc, char *argv[])
 	    int create_target;
 	    struct Cell_head gcpcellhd;
 	    double emin, emax, nmin, nmax;
-	    
+
 	    G_zero(&gcpcellhd, sizeof(struct Cell_head));
 
 	    sPoints.count = GDALGetGCPCount(hDS);
@@ -884,7 +884,7 @@ int main(int argc, char *argv[])
 	    create_target = 0;
 	    if (parm.target->answer) {
 		char target_mapset[GMAPSET_MAX];
-		
+
 		/* does the target location exist? */
 		G_create_alt_env();
 		G_setenv_nogisrc("LOCATION_NAME", parm.target->answer);
@@ -927,7 +927,7 @@ int main(int argc, char *argv[])
 				      &(sPoints.e2[iGCP]),
 				      &(sPoints.n2[iGCP]), NULL) < 0)
 			G_fatal_error(_("Error in %s (can't "
-					"re-project GCP %i)"), 
+					"re-project GCP %i)"),
 				       "GPJ_transform()", iGCP);
 		}
 
@@ -946,13 +946,30 @@ int main(int argc, char *argv[])
 
 	    I_put_control_points(output, &sPoints);
 	    if (create_target) {
+		char *gdalsrid = NULL, *gdalwkt = NULL;
+		OGRSpatialReferenceH hSRS = NULL;
+
+		/* GDAL >= 3 */
+#if GDAL_VERSION_MAJOR >= 3
+		hSRS = GDALGetGCPSpatialRef(hDS);
+		char **papszOptions = NULL;
+#else
+		gdalwkt = G_store(GDALGetGCPProjection(hDS));
+		hSRS = OSRNewSpatialReference(NULL);
+		if (OSRImportFromWkt(hSRS, &gdalwkt) != OGRERR_NONE) {
+		    OSRDestroySpatialReference(hSRS);
+		    hSRS = NULL;
+		}
+#endif
 		/* create target location */
-		if (GPJ_wkt_to_grass(&gcpcellhd, &proj_info,
-				     &proj_units, GDALGetGCPProjection(hDS), 0) < 0) {
+		if (!hSRS || GPJ_osr_to_grass(&gcpcellhd, &proj_info,
+				              &proj_units, hSRS, 0) == 1) {
 		    G_warning(_("Unable to convert input map projection to GRASS "
 				    "format; cannot create new location."));
 		}
 		else {
+		    const char *authkey, *authname, *authcode;
+
 		    gcpcellhd.west = emin;
 		    gcpcellhd.east = emax;
 		    gcpcellhd.south = nmin;
@@ -967,12 +984,41 @@ int main(int argc, char *argv[])
 		    gcpcellhd.bottom = 0.;
 		    gcpcellhd.tb_res = 1.;
 		    gcpcellhd.depths = 1;
-		    
+
 		    G_adjust_Cell_head(&gcpcellhd, 1, 1);
 
+		    /* get id of spatial reference (srid) */
+		    authkey = NULL;
+		    if (OSRIsProjected(hSRS))
+			authkey = "PROJCS";
+		    else if (OSRIsGeographic(hSRS))
+			authkey = "GEOGCS";
+
+		    if (authkey) {
+			authname = OSRGetAuthorityName(hSRS, authkey);
+			if (authname && *authname) {
+			    authcode = OSRGetAuthorityCode(hSRS, authkey);
+			    if (authcode && *authcode) {
+				G_asprintf(&gdalsrid, "%s:%s", authname, authcode);
+			    }
+			}
+		    }
+
+		    /* get WKT of spatial reference */
+#if GDAL_VERSION_MAJOR >= 3
+		    papszOptions = G_calloc(3, sizeof(char *));
+		    papszOptions[0] = G_store("MULTILINE=YES");
+		    papszOptions[1] = G_store("FORMAT=WKT2");
+		    OSRExportToWktEx(hSRS, &gdalwkt, (const char **)papszOptions);
+		    G_free(papszOptions[0]);
+		    G_free(papszOptions[1]);
+		    G_free(papszOptions);
+#endif
 		    G_create_alt_env();
-		    if (0 != G_make_location(parm.target->answer, &gcpcellhd,
-					     proj_info, proj_units)) {
+		    if (0 != G_make_location_crs(parm.target->answer,
+						 &gcpcellhd,
+					         proj_info, proj_units,
+						 gdalsrid, gdalwkt)) {
 			G_fatal_error(_("Unable to create new location <%s>"),
 				      parm.target->answer);
 		    }
@@ -1120,7 +1166,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
     GDALRasterAttributeTableH gdal_rat;
 
     G_message(_("Importing raster map <%s>..."), output);
-    
+
     /* -------------------------------------------------------------------- */
     /*      Select a cell type for the new cell.                            */
     /* -------------------------------------------------------------------- */
@@ -1240,8 +1286,10 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 
 	    G_percent(row, nrows, 2);
 
-	    GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
-			 bufComplex, ncols_gdal, 1, eGDT, 0, 0);
+	    if (GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
+			     bufComplex, ncols_gdal, 1, eGDT, 0, 0) != CE_None) {
+		G_fatal_error(_("Unable to read row %d"), row);
+	    }
 
 	    for (indx = ncols - 1; indx >= 0; indx--) {	/* CEOS: flip east-west during import - MN */
 		if (eGDT == GDT_Int32) {
@@ -1287,7 +1335,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
     }				/* end of complex */
     else if (l1bdriver) {		/* AVHRR */
 	/* AVHRR - read from south to north to match GCPs */
-	/* AVHRR - as for other formats, read from north to south to match GCPs 
+	/* AVHRR - as for other formats, read from north to south to match GCPs
 	 * MM 2013 with gdal 1.10 */
 	for (row = 0; row < nrows; row++) {
 	    if (rowmap[row] < 0)
@@ -1295,8 +1343,10 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 
 	    G_percent(row, nrows, 2);
 
-	    GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
-			 cell_gdal, ncols_gdal, 1, eGDT, 0, 0);
+	    if (GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
+			     cell_gdal, ncols_gdal, 1, eGDT, 0, 0) != CE_None) {
+		G_fatal_error(_("Unable to read row %d"), row);
+	    }
 
 	    if (nullFlags != NULL) {
 		memset(nullFlags, 0, ncols);
@@ -1305,7 +1355,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		    for (indx = 0; indx < ncols; indx++) {
 			if (colmap[indx] < 0)
 			    nullFlags[indx] = 1;
-			else if (bNoDataEnabled && 
+			else if (bNoDataEnabled &&
 			         ((CELL *) cell_gdal)[colmap[indx]] == (GInt32) dfNoData) {
 			    nullFlags[indx] = 1;
 			}
@@ -1317,7 +1367,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		    for (indx = 0; indx < ncols; indx++) {
 			if (colmap[indx] < 0)
 			    nullFlags[indx] = 1;
-			else if (bNoDataEnabled && 
+			else if (bNoDataEnabled &&
 			         ((FCELL *)cell_gdal)[colmap[indx]] == (float)dfNoData) {
 			    nullFlags[indx] = 1;
 			}
@@ -1329,7 +1379,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		    for (indx = 0; indx < ncols; indx++) {
 			if (colmap[indx] < 0)
 			    nullFlags[indx] = 1;
-			else if (bNoDataEnabled && 
+			else if (bNoDataEnabled &&
 			         ((DCELL *)cell_gdal)[colmap[indx]] == dfNoData) {
 			    nullFlags[indx] = 1;
 			}
@@ -1369,8 +1419,10 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 
 	    G_percent(row, nrows, 2);
 
-	    GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
-			 cell_gdal, ncols_gdal, 1, eGDT, 0, 0);
+	    if (GDALRasterIO(hBand, GF_Read, 0, rowmap[row], ncols_gdal, 1,
+			     cell_gdal, ncols_gdal, 1, eGDT, 0, 0) != CE_None) {
+		G_fatal_error(_("Unable to read row %d"), row);
+	    }
 
 	    if (nullFlags != NULL) {
 		memset(nullFlags, 0, ncols);
@@ -1379,7 +1431,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		    for (indx = 0; indx < ncols; indx++) {
 			if (colmap[indx] < 0)
 			    nullFlags[indx] = 1;
-			else if (bNoDataEnabled && 
+			else if (bNoDataEnabled &&
 			         ((CELL *) cell_gdal)[colmap[indx]] == (GInt32) dfNoData) {
 			    nullFlags[indx] = 1;
 			}
@@ -1391,7 +1443,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		    for (indx = 0; indx < ncols; indx++) {
 			if (colmap[indx] < 0)
 			    nullFlags[indx] = 1;
-			else if (bNoDataEnabled && 
+			else if (bNoDataEnabled &&
 			         ((FCELL *)cell_gdal)[colmap[indx]] == (float)dfNoData) {
 			    nullFlags[indx] = 1;
 			}
@@ -1403,7 +1455,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		    for (indx = 0; indx < ncols; indx++) {
 			if (colmap[indx] < 0)
 			    nullFlags[indx] = 1;
-			else if (bNoDataEnabled && 
+			else if (bNoDataEnabled &&
 			         ((DCELL *)cell_gdal)[colmap[indx]] == dfNoData) {
 			    nullFlags[indx] = 1;
 			}
@@ -1482,22 +1534,22 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 
     /* GRASS color rules in metadata? */
     GDALmetadata = GDALGetMetadata(hBand, "");
-    
+
     if (GDALmetadata) {
 	struct Colors colors;
 	DCELL val1, val2;
 	int r1, g1, b1, r2, g2, b2;
 
 	Rast_init_colors(&colors);
-	
+
 	while (GDALmetadata && GDALmetadata[0]) {
 	    G_debug(2, "%s", GDALmetadata[0]);
 
 	    if (!strncmp("COLOR_TABLE_RULE_RGB_", GDALmetadata[0], 21)) {
 		char *p;
-		
+
 		for (p = GDALmetadata[0]; *p != '=' && *p != '\0'; p++);
-		
+
 		if (*p == '=') {
 		    p++;
 		}
@@ -1525,7 +1577,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
     if (!have_colors && gdal_rat != NULL) {
 	nrows = GDALRATGetRowCount(gdal_rat);
 	ncols = GDALRATGetColumnCount(gdal_rat);
-	
+
 	if (nrows > 0 && ncols > 0) {
 	    int minc, maxc, minmaxc;
 	    int rc, gc, bc, rminc, rmaxc, gminc, gmaxc, bminc, bmaxc;
@@ -1536,13 +1588,13 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 	    int cf;
 
 	    Rast_init_colors(&colors);
-	    
+
 	    minc = maxc = minmaxc = -1;
 	    rc = gc = bc = rminc = rmaxc = gminc = gmaxc = bminc = bmaxc = -1;
 
 	    for (indx = 0; indx < ncols; indx++) {
 		 field_use = GDALRATGetUsageOfCol(gdal_rat, indx);
-		 
+
 		 if (field_use == GFU_Min)
 		    minc = indx;
 		 else if (field_use == GFU_Max)
@@ -1568,14 +1620,14 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		 else if (field_use == GFU_BlueMax)
 		    bmaxc = indx;
 	    }
-	    
+
 	    /* guess color range 0, 1 or 0, 255 */
 
 	    if (minc >= 0 && maxc >= 0 && rminc >= 0 && rmaxc >= 0 &&
 		gminc >= 0 && gmaxc >= 0 && bminc >= 0 && bmaxc >= 0) {
 
 		cf = 1;
-		
+
 		/* analyze color rules */
 		for (indx = 0; indx < nrows; indx++) {
 		    val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minc);
@@ -1655,7 +1707,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		}
 	    }
 	    else if (minmaxc >= 0 && rc >= 0 && gc >= 0 && bc >= 0) {
-		
+
 		cf = 1;
 
 		/* analyze color table */
@@ -1699,14 +1751,14 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 			r1 = GDALRATGetValueAsDouble(gdal_rat, indx, rc);
 			g1 = GDALRATGetValueAsDouble(gdal_rat, indx, gc);
 			b1 = GDALRATGetValueAsDouble(gdal_rat, indx, bc);
-			
+
 			Rast_set_d_color(val1, r1 * cf, g1 * cf, b1 * cf, &colors);
 		    }
 		}
 	    }
-	    
+
 	    have_colors = Rast_colors_count(&colors) > 0;
-	    
+
 	    if (have_colors)
 		Rast_write_colors((char *)output, G_mapset(), &colors);
 
@@ -1772,13 +1824,13 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 	    Rast_free_colors(&colors);
 	}
     }
-    
+
     /* categories in raster attribute table? */
-    
+
     if (gdal_rat != NULL) {
 	nrows = GDALRATGetRowCount(gdal_rat);
 	ncols = GDALRATGetColumnCount(gdal_rat);
-	
+
 	if (nrows > 0 && ncols > 0) {
 	    int minc, maxc, minmaxc, namec;
 	    GDALRATFieldUsage field_use;
@@ -1789,7 +1841,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 	    minc = maxc = minmaxc = namec = -1;
 	    for (indx = 0; indx < ncols; indx++) {
 		 field_use = GDALRATGetUsageOfCol(gdal_rat, indx);
-		 
+
 		 if (field_use == GFU_Min)
 		    minc = indx;
 		 else if (field_use == GFU_Max)
@@ -1808,7 +1860,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		    val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minmaxc);
 		    val2 = val1;
 		    label = GDALRATGetValueAsString(gdal_rat, indx, namec);
-		    
+
 		    if (label)
 			Rast_set_d_cat(&val1, &val2, label, &cats);
 		}
@@ -1824,7 +1876,7 @@ static void ImportBand(GDALRasterBandH hBand, const char *output,
 		    val1 = GDALRATGetValueAsDouble(gdal_rat, indx, minc);
 		    val2 = GDALRATGetValueAsDouble(gdal_rat, indx, maxc);
 		    label = GDALRATGetValueAsString(gdal_rat, indx, namec);
-		    
+
 		    if (label)
 			Rast_set_d_cat(&val1, &val2, label, &cats);
 		}
@@ -1847,7 +1899,7 @@ static int dump_rat(GDALRasterBandH hBand, char *outrat, int nBand)
     GDALRasterAttributeTableH gdal_rat;
     FILE *fp;
     char fname[GNAME_MAX];
-    
+
     if ((gdal_rat = GDALGetDefaultRAT(hBand)) == NULL)
 	return 0;
 
@@ -1858,25 +1910,25 @@ static int dump_rat(GDALRasterBandH hBand, char *outrat, int nBand)
 	return 0;
 
     field_type = G_malloc(ncols * sizeof(GDALRATFieldType));
-    
+
     G_snprintf(fname, GNAME_MAX, "%s_%d.csv", outrat, nBand);
     if (!(fp = fopen(fname, "w"))) {
 	int err = errno;
-	
+
 	G_fatal_error(_("Unable to open file <%s>: %s."),
 	              fname, strerror(err));
     }
 
     /* dump column names and usage */
     for (col = 0; col < ncols; col++) {
-	
+
 	if (col)
 	    fprintf(fp, "|");
-	
+
 	field_name = GDALRATGetNameOfCol(gdal_rat, col);
 	fprintf(fp, "%s", field_name);
 	field_use = GDALRATGetUsageOfCol(gdal_rat, col);
-	
+
 	if (field_use == GFU_Generic)
 	    fprintf(fp, " (General purpose field)");
 	else if (field_use == GFU_PixelCount)
@@ -1922,7 +1974,7 @@ static int dump_rat(GDALRasterBandH hBand, char *outrat, int nBand)
 	field_type[col] = GDALRATGetTypeOfCol(gdal_rat, col);
     }
     fprintf(fp, "\n");
-    
+
     /* dump entries */
     for (row = 0; row < nrows; row++) {
 

raster/r.in.gdal/proj.c

@@ -11,21 +11,71 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
 		      int check_only)
 {
     struct Cell_head loc_wind;
-    struct Key_Value *proj_info, *proj_units, *proj_epsg;
-    struct Key_Value *loc_proj_info, *loc_proj_units;
-    const char *wkt;
+    struct Key_Value *proj_info = NULL,
+                     *proj_units = NULL;
+    struct Key_Value *loc_proj_info = NULL,
+                     *loc_proj_units = NULL;
+    char *wkt = NULL, *srid = NULL;
     char error_msg[8096];
     int proj_trouble;
 
     /* -------------------------------------------------------------------- */
-    /*      Fetch the projection in GRASS form.                             */
+    /*      Fetch the projection in GRASS form, SRID, and WKT.              */
     /* -------------------------------------------------------------------- */
-    proj_info = NULL;
-    proj_units = NULL;
-    proj_epsg = NULL;
-    loc_proj_info = NULL;
-    loc_proj_units = NULL;
 
+#if GDAL_VERSION_NUM >= 3000000
+    OGRSpatialReferenceH hSRS;
+
+    hSRS = GDALGetSpatialRef(hDS);
+    if (hSRS) {
+	char **papszOptions;
+
+	/* get WKT2 definition */
+	papszOptions = G_calloc(3, sizeof(char *));
+	papszOptions[0] = G_store("MULTILINE=YES");
+	papszOptions[1] = G_store("FORMAT=WKT2");
+	OSRExportToWktEx(hSRS, &wkt, (const char **)papszOptions);
+	G_free(papszOptions[0]);
+	G_free(papszOptions[1]);
+	G_free(papszOptions);
+    }
+    /* proj_trouble:
+     * 0: valid srs
+     * 1: no srs, default to xy
+     * 2: unreadable srs, default to xy
+     */
+
+    /* Projection only required for checking so convert non-interactively */
+    proj_trouble = 0;
+    if (wkt && *wkt) {
+	if (hSRS != NULL)
+	    GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0);
+
+	if (!hSRS || (!OSRIsProjected(hSRS) && !OSRIsGeographic(hSRS))) {
+	    G_important_message(_("Input contains an invalid SRS. "
+	                          "WKT definition:\n%s"), wkt);
+
+	    proj_trouble = 2;
+	}
+	else{
+	    const char *authkey, *authname, *authcode;
+
+	    if (OSRIsProjected(hSRS))
+		authkey = "PROJCS";
+	    else /* is geographic */
+		authkey = "GEOGCS";
+
+	    authname = OSRGetAuthorityName(hSRS, authkey);
+	    if (authname && *authname) {
+		authcode = OSRGetAuthorityCode(hSRS, authkey);
+		if (authcode && *authcode) {
+		    G_asprintf(&srid, "%s:%s", authname, authcode);
+		}
+	    }
+	}
+    }
+
+#else
     wkt = GDALGetProjectionRef(hDS);
     /* proj_trouble:
      * 0: valid srs
@@ -43,7 +93,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
 	    GPJ_osr_to_grass(cellhd, &proj_info, &proj_units, hSRS, 0);
 
 	if (!hSRS || (!OSRIsProjected(hSRS) && !OSRIsGeographic(hSRS))) {
-	    G_important_message(_("Input contains an invalid SRS. " 
+	    G_important_message(_("Input contains an invalid SRS. "
 	                          "WKT definition:\n%s"), wkt);
 
 	    proj_trouble = 2;
@@ -57,18 +107,17 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
 		authkey = "GEOGCS";
 
 	    authname = OSRGetAuthorityName(hSRS, authkey);
-	    if (authname && *authname && strcmp(authname, "EPSG") == 0) {
+	    if (authname && *authname) {
 		authcode = OSRGetAuthorityCode(hSRS, authkey);
 		if (authcode && *authcode) {
-		    G_debug(1, "found EPSG:%s", authcode);
-		    proj_epsg = G_create_key_value();
-		    G_set_key_value("epsg", authcode, proj_epsg);
+		    G_asprintf(&srid, "%s:%s", authname, authcode);
 		}
 	    }
 	}
 	if (hSRS)
 	    OSRDestroySpatialReference(hSRS);
     }
+#endif
     else {
 	G_important_message(_("No projection information available"));
 	cellhd->proj = PROJECTION_XY;
@@ -80,14 +129,14 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
     /*      Do we need to create a new location?                            */
     /* -------------------------------------------------------------------- */
     if (outloc != NULL) {
-	/* do not create a xy location if an existing SRS was unreadable */ 
+	/* do not create a xy location if an existing SRS was unreadable */
 	if (proj_trouble == 2) {
 	    G_fatal_error(_("Unable to convert input map projection to GRASS "
 			    "format; cannot create new location."));
 	}
 	else {
-            if (0 != G_make_location_epsg(outloc, cellhd, proj_info,
-	                                  proj_units, proj_epsg)) {
+            if (0 != G_make_location_crs(outloc, cellhd, proj_info,
+	                                 proj_units, srid, wkt)) {
                 G_fatal_error(_("Unable to create new location <%s>"),
                               outloc);
             }
@@ -282,7 +331,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
                        _("Consider generating a new location from the input dataset using "
                          "the 'location' parameter.\n"));
             }
-            
+
 	    if (check_only)
 		msg_fn = G_message;
 	    else
@@ -297,7 +346,7 @@ void check_projection(struct Cell_head *cellhd, GDALDatasetH hDS,
 	    if (check_only)
 		msg_fn = G_message;
 	    else
-		msg_fn = G_verbose_message;            
+		msg_fn = G_verbose_message;
 	    msg_fn(_("Projection of input dataset and current location "
 		     "appear to match"));
 	    if (check_only) {

raster/r.out.gdal/Makefile

@@ -2,7 +2,7 @@ MODULE_TOPDIR = ../..
 
 PGM  = r.out.gdal
 
-LIBES = $(GPROJLIB) $(IMAGERYLIB) $(DBMILIB) $(RASTERLIB) $(GISLIB) $(GDALLIBS) $(MATHLIB)
+LIBES = $(GPROJLIB) $(IMAGERYLIB) $(DBMILIB) $(RASTERLIB) $(GISLIB) $(GDALLIBS) $(PROJLIB) $(MATHLIB)
 DEPENDENCIES = $(GPROJDEP) $(IMAGERYDEP) $(DBMIDEP) $(RASTERDEP) $(GISDEP)
 EXTRA_INC = $(PROJINC) $(GDALCFLAGS)
 EXTRA_CFLAGS = -DGRASS_VERSION_NUMBER=\"'$(GRASS_VERSION_NUMBER)'\" \

raster/r.out.gdal/main.c

@@ -195,7 +195,7 @@ int main(int argc, char *argv[])
 #endif
     format->answer = "GTiff";
     format->required = YES;
-    
+
     type = G_define_option();
     type->key = "type";
     type->type = TYPE_STRING;
@@ -204,7 +204,7 @@ int main(int argc, char *argv[])
 	"Byte,Int16,UInt16,Int32,UInt32,Float32,Float64,CInt16,CInt32,CFloat32,CFloat64";
     type->required = NO;
     type->guisection = _("Creation");
- 
+
     createopt = G_define_option();
     createopt->key = "createopt";
     createopt->type = TYPE_STRING;
@@ -239,7 +239,7 @@ int main(int argc, char *argv[])
     nodataopt->multiple = NO;
     nodataopt->required = NO;
     nodataopt->guisection = _("Creation");
-    
+
     overviewopt = G_define_option();
     overviewopt->key = "overviews";
     overviewopt->type = TYPE_INTEGER;
@@ -264,7 +264,7 @@ int main(int argc, char *argv[])
 	supported_formats(&gdal_formats);
 	exit(EXIT_SUCCESS);
     }
-    
+
     /* Find input GRASS raster.. */
     mapset = G_find_raster2(input->answer, "");
 
@@ -286,7 +286,71 @@ int main(int argc, char *argv[])
     struct Key_Value *projinfo = G_get_projinfo();
     struct Key_Value *projunits = G_get_projunits();
     struct Key_Value *projepsg = G_get_projepsg();
-    char *srswkt = GPJ_grass_to_wkt2(projinfo, projunits, projepsg, 0, 0);
+    char *srswkt = NULL;
+
+#if GDAL_VERSION_MAJOR >= 3 && PROJ_VERSION_MAJOR >= 6
+    char *indef;
+
+    if ((indef = G_get_projsrid())) {
+	PJ *obj = NULL;
+
+	if ((obj = proj_create(NULL, indef))) {
+	    srswkt = G_store(proj_as_wkt(NULL, obj, PJ_WKT2_LATEST, NULL));
+
+	    if (srswkt && !*srswkt) {
+		G_free(srswkt);
+		srswkt = NULL;
+	    }
+	    proj_destroy(obj);
+	}
+    }
+    if (!srswkt && (indef = G_get_projwkt())) {
+	OGRSpatialReferenceH hSRS;
+	char **papszOptions = NULL;
+
+	hSRS = OSRNewSpatialReference(indef);
+	papszOptions = G_calloc(2, sizeof(char *));
+	papszOptions[0] = G_store("FORMAT=WKT2");
+	OSRExportToWktEx(hSRS, &srswkt, (const char **)papszOptions);
+	G_free(papszOptions[0]);
+	G_free(papszOptions);
+	srswkt = G_store(srswkt);
+	if (srswkt && !*srswkt) {
+	    G_free(srswkt);
+	    srswkt = NULL;
+	}
+	OSRDestroySpatialReference(hSRS);
+    }
+#endif
+    if (!srswkt) {
+	srswkt = GPJ_grass_to_wkt2(projinfo, projunits, projepsg, 0, 0);
+
+#if GDAL_VERSION_MAJOR >= 3 && PROJ_VERSION_MAJOR >= 6
+	/* convert bound CRS */
+	{
+	    PJ *obj = NULL;
+
+	    indef = srswkt;
+	    if ((obj = proj_create(NULL, indef))) {
+		if (proj_get_type(obj) == PJ_TYPE_BOUND_CRS) {
+		    PJ *source_crs;
+
+		    G_debug(1, "found bound crs");
+		    source_crs = proj_get_source_crs(NULL, obj);
+		    if (source_crs) {
+			srswkt = G_store(proj_as_wkt(NULL, source_crs, PJ_WKT2_LATEST, NULL));
+			if (srswkt && !*srswkt) {
+			    G_free(srswkt);
+			    srswkt = NULL;
+			}
+			proj_destroy(source_crs);
+		    }
+		}
+		proj_destroy(obj);
+	    }
+	}
+#endif
+    }
 
     G_get_window(&cellhead);
 
@@ -298,7 +362,7 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Unable to get <%s> driver"), format->answer);
     /* Does driver support GDALCreate ? */
     if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, NULL) == NULL) {
-	/* If not - create MEM driver for intermediate dataset. 
+	/* If not - create MEM driver for intermediate dataset.
 	 * Check if raster can be created at all (with GDALCreateCopy) */
 	if (GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATECOPY, NULL)) {
 	    G_message(_("Driver <%s> does not support direct writing. "
@@ -454,7 +518,7 @@ int main(int argc, char *argv[])
     /* if GDAL datatype set by user, do checks */
     if (type->answer) {
 
-	/* Check if raster data range is outside of the range of 
+	/* Check if raster data range is outside of the range of
 	 * given GDAL datatype, not even overlapping */
 	if (range_check(export_min, export_max, datatype))
 	    G_fatal_error(_("Raster export would result in complete data loss, aborting."));
@@ -668,7 +732,7 @@ int main(int argc, char *argv[])
 	int i, oi, *ol;
 
 	G_message(_("Building overviews ..."));
-	
+
 	ol = G_malloc(n_overviews * sizeof(int));
 	oi = 2;
 	for (i = 0; i < n_overviews; i++) {
@@ -681,7 +745,7 @@ int main(int argc, char *argv[])
 	}
     }
 
-    /* Finally create user requested raster format from memory raster 
+    /* Finally create user requested raster format from memory raster
      * if in-memory driver was used */
     if (hMEMDS) {
 	hDstDS =

raster/r.proj/bordwalk.c

@@ -88,8 +88,11 @@ static void proj_update(const struct pj_info *from_pj, const struct pj_info *to_
 			struct Cell_head *to_hd, double hx, double hy)
 {
 	if (GPJ_transform(from_pj, to_pj, trans_pj, dir,
-			  &hx, &hy, NULL) < 0)
+			  &hx, &hy, NULL) < 0) {
+	    G_fatal_error(_("unable to transform coordinates %g, %g"),
+			  hx, hy);
 	    return;
+	}
 	update(to_hd, hx, hy);
 }
 
@@ -132,8 +135,12 @@ static int proj_inside(const struct pj_info *from_pj, const struct pj_info *to_p
 		       const struct pj_info *trans_pj, int dir, 
 		       const struct Cell_head *ref_hd, double hx, double hy)
 {
-    if (GPJ_transform(from_pj, to_pj, trans_pj, -dir, &hx, &hy, NULL) < 0)
+    if (GPJ_transform(from_pj, to_pj, trans_pj, -dir, &hx, &hy, NULL) < 0) {
+	G_fatal_error(_("unable to transform coordinates %g, %g"),
+		      hx, hy);
+
 	return 0;
+    }
     return inside(ref_hd, hx, hy);
 }
 
@@ -247,8 +254,9 @@ void bordwalk_edge(const struct Cell_head *from_hd, struct Cell_head *to_hd,
     hy = (from_hd->north + from_hd->south) / 2.0;
 
     if (GPJ_transform(from_pj, to_pj, trans_pj, dir,
-		      &hx, &hy, NULL) < 0)
-	G_fatal_error(_("Unable to reproject map center"));
+		      &hx, &hy, NULL) < 0) {
+	G_fatal_error(_("Unable to reproject map center %g, %g"), hx, hy);
+    }
 
     to_hd->east  = hx;
     to_hd->west  = hx;

raster/r.proj/main.c

@@ -8,8 +8,8 @@
 *		 Dept. of Geography
 *		 emes@geo0.geog.uni-heidelberg.de
 *
-* 		 (With the help of a lot of existing GRASS sources, in 
-*		  particular v.proj) 
+* 		 (With the help of a lot of existing GRASS sources, in
+*		  particular v.proj)
 *
 * PURPOSE:      r.proj converts a map to a new geographic projection. It reads a
 *	        map from a different location, projects it and write it out
@@ -26,27 +26,27 @@
 * Changes
 *		 Morten Hulden <morten@untamo.net>, Aug 2000:
 *		 - aborts if input map is outside current location.
-*		 - can handle projections (conic, azimuthal etc) where 
-*		 part of the map may fall into areas where south is 
+*		 - can handle projections (conic, azimuthal etc) where
+*		 part of the map may fall into areas where south is
 *		 upward and east is leftward.
 *		 - avoids passing location edge coordinates to PROJ
 *		 (they may be invalid in some projections).
 *		 - output map will be clipped to borders of the current region.
-*		 - output map cell edges and centers will coinside with those 
+*		 - output map cell edges and centers will coinside with those
 *		 of the current region.
 *		 - output map resolution (unless changed explicitly) will
 *		 match (exactly) the resolution of the current region.
-*		 - if the input map is smaller than the current region, the 
+*		 - if the input map is smaller than the current region, the
 *		 output map will only cover the overlapping area.
 *                - if the input map is larger than the current region, only the
 *		 needed amount of memory will be allocated for the projection
-*	
+*
 *		 Bugfixes 20050328: added floor() before (int) typecasts to in avoid
-*		 asymmetrical rounding errors. Added missing offset outcellhd.ew_res/2 
-*		 to initial xcoord for each row in main projection loop (we want to  project 
+*		 asymmetrical rounding errors. Added missing offset outcellhd.ew_res/2
+*		 to initial xcoord for each row in main projection loop (we want to  project
 *		 center of cell, not border).
 *
-*                Glynn Clements 2006: Use G_interp_* functions, modified      
+*                Glynn Clements 2006: Use G_interp_* functions, modified
 *                  version of r.proj which uses a tile cache instead  of loading the
 *                  entire map into memory.
 *                Markus Metz 2010: lanczos and lanczos fallback interpolation methods
@@ -165,7 +165,7 @@ int main(int argc, char **argv)
     inmap->description = _("Name of input raster map to re-project");
     inmap->required = NO;
     inmap->guisection = _("Source");
-    
+
     indbase = G_define_standard_option(G_OPT_M_DBASE);
     indbase->label = _("Path to GRASS database of input location");
 
@@ -175,7 +175,7 @@ int main(int argc, char **argv)
     outmap->guisection = _("Target");
 
     ipolname = make_ipol_list();
-    
+
     interpol = G_define_option();
     interpol->key = "method";
     interpol->type = TYPE_STRING;
@@ -207,7 +207,7 @@ int main(int argc, char **argv)
     list->key = 'l';
     list->description = _("List raster maps in input mapset and exit");
     list->guisection = _("Print");
-    
+
     nocrop = G_define_flag();
     nocrop->key = 'n';
     nocrop->description = _("Do not perform region cropping optimization");
@@ -217,7 +217,7 @@ int main(int argc, char **argv)
     print_bounds->description =
 	_("Print input map's bounds in the current projection and exit");
     print_bounds->guisection = _("Print");
-    
+
     gprint_bounds = G_define_flag();
     gprint_bounds->key = 'g';
     gprint_bounds->description =
@@ -274,6 +274,9 @@ int main(int argc, char **argv)
     if (pj_get_kv(&oproj, out_proj_info, out_unit_info) < 0)
 	G_fatal_error(_("Unable to get projection key values of output raster map"));
 
+    oproj.srid = G_get_projsrid();
+    oproj.wkt = G_get_projwkt();
+
     /* Change the location           */
     G_create_alt_env();
     G_setenv_nogisrc("GISDBASE", indbase->answer ? indbase->answer : G_gisdbase());
@@ -330,6 +333,10 @@ int main(int argc, char **argv)
     if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
 	G_fatal_error(_("Unable to get projection key values of input map"));
 
+    iproj.srid = G_get_projsrid();
+    iproj.wkt = G_get_projwkt();
+
+    tproj.pj = NULL;
     tproj.def = NULL;
 #ifdef HAVE_PROJ_H
     if (pipeline->answer) {
@@ -337,11 +344,6 @@ int main(int argc, char **argv)
     }
 #endif
 
-    /* switch back to current location */
-    G_switch_env();
-    if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
-	G_fatal_error(_("Unable to initialize coordinate transformation"));
-
     G_free_key_value(in_proj_info);
     G_free_key_value(in_unit_info);
     G_free_key_value(out_proj_info);
@@ -349,14 +351,9 @@ int main(int argc, char **argv)
     if (G_verbose() > G_verbose_std())
 	pj_print_proj_params(&iproj, &oproj);
 
-    /* switch to input location */
-    G_switch_env();
-
-    /* this call causes r.proj to read the entire map into memeory */
+    /* this call causes r.proj to read the entire map into memory */
     Rast_get_cellhd(inmap->answer, setname, &incellhd);
 
-    Rast_set_input_window(&incellhd);
-
     if (G_projection() == PROJECTION_XY)
 	G_fatal_error(_("Unable to work with unprojected data (xy location)"));
 
@@ -380,6 +377,17 @@ int main(int argc, char **argv)
 	G_message(_("Input map <%s@%s> in location <%s>:"),
 	    inmap->answer, setname, inlocation->answer);
 
+	/* reproject input raster extents from input to output */
+	G_set_window(&incellhd);
+
+	G_debug(1, "input window north: %.8f", incellhd.north);
+	G_debug(1, "input window south: %.8f", incellhd.south);
+	G_debug(1, "input window east: %.8f", incellhd.east);
+	G_debug(1, "input window west: %.8f", incellhd.west);
+
+	if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
+	    G_fatal_error(_("Unable to initialize coordinate transformation"));
+
 	outcellhd.north = -1e9;
 	outcellhd.south =  1e9;
 	outcellhd.east  = -1e9;
@@ -411,10 +419,30 @@ int main(int argc, char **argv)
 	exit(EXIT_SUCCESS);
     }
 
-
     /* Cut non-overlapping parts of input map */
-    if (!nocrop->answer)
-	bordwalk(&outcellhd, &incellhd, &iproj, &oproj, &tproj, PJ_INV);
+    if (!nocrop->answer) {
+	/* reproject current region from output to input */
+	/* switch back to current location,
+	 * initialize transformation pipeline */
+	G_switch_env();
+	G_unset_window();
+	G_set_window(&outcellhd);
+	tproj.def = NULL;
+	tproj.pj = NULL;
+#ifdef HAVE_PROJ_H
+	if (pipeline->answer) {
+	    tproj.def = G_store(pipeline->answer);
+	}
+#endif
+	if (GPJ_init_transform(&oproj, &iproj, &tproj) < 0)
+	    G_fatal_error(_("Unable to initialize coordinate transformation"));
+
+	/* switch to input location */
+	G_switch_env();
+
+	/* update cellhead of input map */
+	bordwalk(&outcellhd, &incellhd, &oproj, &iproj, &tproj, PJ_FWD);
+    }
 
     /* Add 2 cells on each side for bilinear/cubic & future interpolation methods */
     /* (should probably be a factor based on input and output resolution) */
@@ -431,16 +459,35 @@ int main(int argc, char **argv)
     if (incellhd.west < iwest)
 	incellhd.west = iwest;
 
-    Rast_set_input_window(&incellhd);
-
-    /* And switch back to original location */
+    /* switch to current location */
 
     G_switch_env();
 
     /* Adjust borders of output map */
 
-    if (!nocrop->answer)
+    if (!nocrop->answer) {
+	/* reproject from input to output */
+	/* switch input location,
+	 * initialize transformation pipeline */
+	G_switch_env();
+	G_unset_window();
+	G_set_window(&incellhd);
+	tproj.def = NULL;
+	tproj.pj = NULL;
+#ifdef HAVE_PROJ_H
+	if (pipeline->answer) {
+	    tproj.def = G_store(pipeline->answer);
+	}
+#endif
+	if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
+	    G_fatal_error(_("Unable to initialize coordinate transformation"));
+
+	/* switch to output location */
+	G_switch_env();
+
+	/* reduce output region */
 	bordwalk(&incellhd, &outcellhd, &iproj, &oproj, &tproj, PJ_FWD);
+    }
 
 #if 0
     outcellhd.west = outcellhd.south = HUGE_VAL;
@@ -501,9 +548,23 @@ int main(int argc, char **argv)
     ibuffer = readcell(fdi, memory->answer);
     Rast_close(fdi);
 
+    /* And switch back to original location */
     G_switch_env();
     Rast_set_output_window(&outcellhd);
 
+    /* reproject from output to input */
+    G_unset_window();
+    G_set_window(&outcellhd);
+    tproj.def = NULL;
+    tproj.pj = NULL;
+#ifdef HAVE_PROJ_H
+    if (pipeline->answer) {
+	tproj.def = G_store(pipeline->answer);
+    }
+#endif
+    if (GPJ_init_transform(&oproj, &iproj, &tproj) < 0)
+	G_fatal_error(_("Unable to initialize coordinate transformation"));
+
     if (strcmp(interpol->answer, "nearest") == 0) {
 	fdo = Rast_open_new(mapname, cell_type);
 	obuffer = (CELL *) Rast_allocate_output_buf(cell_type);
@@ -527,7 +588,7 @@ int main(int argc, char **argv)
 
 #if 0
 	/* parallelization does not always work,
-	 * segfaults in the interpolation functions 
+	 * segfaults in the interpolation functions
 	 * can happen */
         #pragma omp parallel for schedule (static)
 #endif
@@ -540,9 +601,9 @@ int main(int argc, char **argv)
 
 	    /* project coordinates in output matrix to       */
 	    /* coordinates in input matrix                   */
-	    if (GPJ_transform(&iproj, &oproj, &tproj, PJ_INV,
+	    if (GPJ_transform(&oproj, &iproj, &tproj, PJ_FWD,
 			      &xcoord1, &ycoord1, NULL) < 0) {
-		G_warning(_("Error in %s"), "GPJ_transform()");
+		G_fatal_error(_("Error in %s"), "GPJ_transform()");
 		Rast_set_null_value(obufptr, 1, cell_type);
 	    }
 	    else {
@@ -612,10 +673,10 @@ char *make_ipol_desc(void)
 
     for (i = 0; menu[i].name; i++)
 	size += strlen(menu[i].name) + strlen(menu[i].text) + 2;
-    
+
     buf = G_malloc(size);
     *buf = '\0';
-    
+
     for (i = 0; menu[i].name; i++) {
 	if (i)
 	    strcat(buf, ";");

vector/v.external/proj.c

@@ -12,16 +12,16 @@
 int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
 		   struct Key_Value **proj_info,
 		   struct Key_Value **proj_units,
-		   struct Key_Value **proj_epsg,
+		   char **proj_srid, char **proj_wkt,
 		   char *geom_col, int verbose)
 {
     OGRSpatialReferenceH hSRS;
-    char *pszProj4 = NULL;
 
     hSRS = NULL;
     *proj_info = NULL;
     *proj_units = NULL;
-    *proj_epsg = NULL;
+    *proj_srid = NULL;
+    *proj_wkt = NULL;
 
     /* Fetch input layer projection in GRASS form. */
 #if GDAL_VERSION_NUM >= 1110000
@@ -29,7 +29,7 @@ int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
 	int igeom;
         OGRGeomFieldDefnH Ogr_geomdefn;
 	OGRFeatureDefnH Ogr_featuredefn;
-        
+
         Ogr_featuredefn = OGR_L_GetLayerDefn(Ogr_layer);
         igeom = OGR_FD_GetGeomFieldIndex(Ogr_featuredefn, geom_col);
         if (igeom < 0)
@@ -66,7 +66,7 @@ int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
 
 	return 2;
     }
-    /* custom checks because if in doubt GPJ_osr_to_grass() returns a 
+    /* custom checks because if in doubt GPJ_osr_to_grass() returns a
      * xy CRS */
     if (hSRS == NULL) {
 	if (verbose) {
@@ -96,6 +96,20 @@ int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
     }
     else{
 	const char *authkey, *authname, *authcode;
+#if GDAL_VERSION_NUM >= 3000000
+	char **papszOptions;
+
+	/* get WKT2 definition */
+	papszOptions = G_calloc(3, sizeof(char *));
+	papszOptions[0] = G_store("MULTILINE=YES");
+	papszOptions[1] = G_store("FORMAT=WKT2");
+	OSRExportToWktEx(hSRS, proj_wkt, (const char **)papszOptions);
+	G_free(papszOptions[0]);
+	G_free(papszOptions[1]);
+	G_free(papszOptions);
+#else
+	OSRExportToWkt(hSRS, proj_wkt);
+#endif
 
 	if (OSRIsProjected(hSRS))
 	    authkey = "PROJCS";
@@ -103,31 +117,12 @@ int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
 	    authkey = "GEOGCS";
 
 	authname = OSRGetAuthorityName(hSRS, authkey);
-	if (authname && *authname && strcmp(authname, "EPSG") == 0) {
+	if (authname && *authname) {
 	    authcode = OSRGetAuthorityCode(hSRS, authkey);
 	    if (authcode && *authcode) {
-		*proj_epsg = G_create_key_value();
-		G_set_key_value("epsg", authcode, *proj_epsg);
-	    }
-	}
-    }
-
-    if (OSRExportToProj4(hSRS, &pszProj4) != OGRERR_NONE) {
-	G_important_message(_("Projection for layer <%s> can not be converted to proj4"),
-			    OGR_L_GetName(Ogr_layer));
-
-	if (verbose) {
-	    char *wkt = NULL;
-
-	    if (OSRExportToPrettyWkt(hSRS, &wkt, FALSE) != OGRERR_NONE) {
-		G_important_message(_("Can't get WKT-style parameter string"));
-	    }
-	    else if (wkt) {
-		G_important_message(_("WKT-style definition:\n%s"), wkt);
+		G_asprintf(proj_srid, "%s:%s", authname, authcode);
 	    }
 	}
-
-	return 2;
     }
 
     return 0;
@@ -139,8 +134,11 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
 		      int check_only)
 {
     struct Cell_head loc_wind;
-    struct Key_Value *proj_info, *proj_units, *proj_epsg;
-    struct Key_Value *loc_proj_info, *loc_proj_units;
+    struct Key_Value *proj_info = NULL,
+                     *proj_units = NULL;
+    struct Key_Value *loc_proj_info = NULL,
+                     *loc_proj_units = NULL;
+    char *wkt = NULL, *srid = NULL;
     char error_msg[8096];
     int proj_trouble;
     OGRLayerH Ogr_layer;
@@ -149,13 +147,8 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
     Ogr_layer = ds_getlayerbyindex(hDS, layer);
 
     /* -------------------------------------------------------------------- */
-    /*      Fetch the projection in GRASS form.                             */
+    /*      Fetch the projection in GRASS form, SRID, and WKT.              */
     /* -------------------------------------------------------------------- */
-    proj_info = NULL;
-    proj_units = NULL;
-    proj_epsg = NULL;
-    loc_proj_info = NULL;
-    loc_proj_units = NULL;
 
     /* proj_trouble:
      * 0: valid srs
@@ -165,21 +158,21 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
 
     /* Projection only required for checking so convert non-interactively */
     proj_trouble = get_layer_proj(Ogr_layer, cellhd, &proj_info, &proj_units,
-		                  &proj_epsg, geom_col, 1);
+		                  &srid, &wkt, geom_col, 1);
 
     /* -------------------------------------------------------------------- */
     /*      Do we need to create a new location?                            */
     /* -------------------------------------------------------------------- */
     if (outloc != NULL) {
 	/* do not create a xy location because this can mean that the
-	 * real SRS has not been recognized or is missing */ 
+	 * real SRS has not been recognized or is missing */
 	if (proj_trouble) {
 	    G_fatal_error(_("Unable to convert input map projection to GRASS "
 			    "format; cannot create new location."));
 	}
 	else {
-            if (0 != G_make_location_epsg(outloc, cellhd, proj_info,
-	                                  proj_units, proj_epsg)) {
+            if (0 != G_make_location_crs(outloc, cellhd, proj_info,
+	                                 proj_units, srid, wkt)) {
                 G_fatal_error(_("Unable to create new location <%s>"),
                               outloc);
             }
@@ -375,7 +368,7 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
                        _("Consider generating a new location from the input dataset using "
                          "the 'location' parameter.\n"));
             }
-            
+
 	    if (check_only)
 		msg_fn = G_message;
 	    else
@@ -390,7 +383,7 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
 	    if (check_only)
 		msg_fn = G_message;
 	    else
-		msg_fn = G_verbose_message;            
+		msg_fn = G_verbose_message;
 	    msg_fn(_("Projection of input dataset and current location "
 		     "appear to match"));
 	    if (check_only) {

vector/v.in.ogr/proj.c

@@ -3,6 +3,7 @@
 #include <grass/glocale.h>
 
 #include <ogr_srs_api.h>
+#include <cpl_conv.h>
 #include "global.h"
 
 /* get projection info of OGR layer in GRASS format
@@ -12,16 +13,16 @@
 int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
 		   struct Key_Value **proj_info,
 		   struct Key_Value **proj_units,
-		   struct Key_Value **proj_epsg,
+		   char **proj_srid, char **proj_wkt,
 		   char *geom_col, int verbose)
 {
     OGRSpatialReferenceH hSRS;
-    char *pszProj4 = NULL;
 
     hSRS = NULL;
     *proj_info = NULL;
     *proj_units = NULL;
-    *proj_epsg = NULL;
+    *proj_srid = NULL;
+    *proj_wkt = NULL;
 
     /* Fetch input layer projection in GRASS form. */
 #if GDAL_VERSION_NUM >= 1110000
@@ -29,7 +30,7 @@ int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
 	int igeom;
         OGRGeomFieldDefnH Ogr_geomdefn;
 	OGRFeatureDefnH Ogr_featuredefn;
-        
+
         Ogr_featuredefn = OGR_L_GetLayerDefn(Ogr_layer);
         igeom = OGR_FD_GetGeomFieldIndex(Ogr_featuredefn, geom_col);
         if (igeom < 0)
@@ -66,7 +67,7 @@ int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
 
 	return 2;
     }
-    /* custom checks because if in doubt GPJ_osr_to_grass() returns a 
+    /* custom checks because if in doubt GPJ_osr_to_grass() returns a
      * xy CRS */
     if (hSRS == NULL) {
 	if (verbose) {
@@ -96,6 +97,20 @@ int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
     }
     else{
 	const char *authkey, *authname, *authcode;
+#if GDAL_VERSION_NUM >= 3000000
+	char **papszOptions;
+
+	/* get WKT2 definition */
+	papszOptions = G_calloc(3, sizeof(char *));
+	papszOptions[0] = G_store("MULTILINE=YES");
+	papszOptions[1] = G_store("FORMAT=WKT2");
+	OSRExportToWktEx(hSRS, proj_wkt, (const char **)papszOptions);
+	G_free(papszOptions[0]);
+	G_free(papszOptions[1]);
+	G_free(papszOptions);
+#else
+	OSRExportToWkt(hSRS, proj_wkt);
+#endif
 
 	if (OSRIsProjected(hSRS))
 	    authkey = "PROJCS";
@@ -103,33 +118,14 @@ int get_layer_proj(OGRLayerH Ogr_layer, struct Cell_head *cellhd,
 	    authkey = "GEOGCS";
 
 	authname = OSRGetAuthorityName(hSRS, authkey);
-	if (authname && *authname && strcmp(authname, "EPSG") == 0) {
+	if (authname && *authname) {
 	    authcode = OSRGetAuthorityCode(hSRS, authkey);
 	    if (authcode && *authcode) {
-		*proj_epsg = G_create_key_value();
-		G_set_key_value("epsg", authcode, *proj_epsg);
+		G_asprintf(proj_srid, "%s:%s", authname, authcode);
 	    }
 	}
     }
 
-    if (OSRExportToProj4(hSRS, &pszProj4) != OGRERR_NONE) {
-	G_important_message(_("Projection for layer <%s> can not be converted to proj4"),
-			    OGR_L_GetName(Ogr_layer));
-
-	if (verbose) {
-	    char *wkt = NULL;
-
-	    if (OSRExportToPrettyWkt(hSRS, &wkt, FALSE) != OGRERR_NONE) {
-		G_important_message(_("Can't get WKT-style parameter string"));
-	    }
-	    else if (wkt) {
-		G_important_message(_("WKT-style definition:\n%s"), wkt);
-	    }
-	}
-
-	return 2;
-    }
-
     return 0;
 }
 
@@ -141,16 +137,20 @@ int cmp_layer_srs(ds_t Ogr_ds, int nlayers, int *layers,
 		  char **layer_names, char *geom_col)
 {
     int layer;
-    struct Key_Value *proj_info1, *proj_units1, *proj_epsg1;
-    struct Key_Value *proj_info2, *proj_units2, *proj_epsg2;
+    struct Key_Value *proj_info1, *proj_units1;
+    char *proj_srid1, *proj_wkt1;
+    struct Key_Value *proj_info2, *proj_units2;
+    char *proj_srid2, *proj_wkt2;
     struct Cell_head cellhd1, cellhd2;
     OGRLayerH Ogr_layer;
 
     if (nlayers == 1)
 	return 0;
 
-    proj_info1 = proj_units1 = proj_epsg1 = NULL;
-    proj_info2 = proj_units2 = proj_epsg2 = NULL;
+    proj_info1 = proj_units1 = NULL;
+    proj_srid1 = proj_wkt1 = NULL;
+    proj_info2 = proj_units2 = NULL;
+    proj_srid2 = proj_wkt2 = NULL;
 
     G_get_window(&cellhd1);
     layer = 0;
@@ -159,7 +159,7 @@ int cmp_layer_srs(ds_t Ogr_ds, int nlayers, int *layers,
 	Ogr_layer = ds_getlayerbyindex(Ogr_ds, layers[layer]);
 
 	if (get_layer_proj(Ogr_layer, &cellhd1, &proj_info1, &proj_units1,
-			   &proj_epsg1, geom_col, 0) == 0) {
+			   &proj_srid1, &proj_wkt1, geom_col, 0) == 0) {
 	    break;
 	}
 	layer++;
@@ -173,8 +173,10 @@ int cmp_layer_srs(ds_t Ogr_ds, int nlayers, int *layers,
 	    G_free_key_value(proj_info1);
 	if (proj_units1)
 	    G_free_key_value(proj_units1);
-	if (proj_epsg1)
-	    G_free_key_value(proj_epsg1);
+	if (proj_srid1)
+	    G_free(proj_srid1);
+	if (proj_wkt1)
+	    CPLFree(proj_wkt1);
 
 	return 0;
     }
@@ -187,8 +189,10 @@ int cmp_layer_srs(ds_t Ogr_ds, int nlayers, int *layers,
 	    G_free_key_value(proj_info1);
 	if (proj_units1)
 	    G_free_key_value(proj_units1);
-	if (proj_epsg1)
-	    G_free_key_value(proj_epsg1);
+	if (proj_srid1)
+	    G_free(proj_srid1);
+	if (proj_wkt1)
+	    CPLFree(proj_wkt1);
 
 	return 1;
     }
@@ -198,13 +202,15 @@ int cmp_layer_srs(ds_t Ogr_ds, int nlayers, int *layers,
 	Ogr_layer = ds_getlayerbyindex(Ogr_ds, layers[layer]);
 	G_get_window(&cellhd2);
 	if (get_layer_proj(Ogr_layer, &cellhd2, &proj_info2, &proj_units2,
-			   &proj_epsg2, geom_col, 0) != 0) {
+			   &proj_srid2, &proj_wkt2, geom_col, 0) != 0) {
 	    if (proj_info1)
 		G_free_key_value(proj_info1);
 	    if (proj_units1)
 		G_free_key_value(proj_units1);
-	    if (proj_epsg1)
-		G_free_key_value(proj_epsg1);
+	    if (proj_srid1)
+		G_free(proj_srid1);
+	    if (proj_wkt1)
+		CPLFree(proj_wkt1);
 
 	    return 1;
 	}
@@ -216,15 +222,19 @@ int cmp_layer_srs(ds_t Ogr_ds, int nlayers, int *layers,
 		G_free_key_value(proj_info1);
 	    if (proj_units1)
 		G_free_key_value(proj_units1);
-	    if (proj_epsg1)
-		G_free_key_value(proj_epsg1);
+	    if (proj_srid1)
+		G_free(proj_srid1);
+	    if (proj_wkt1)
+		CPLFree(proj_wkt1);
 	    if (proj_info2)
 		G_free_key_value(proj_info2);
 	    if (proj_units2)
 		G_free_key_value(proj_units2);
-	    if (proj_epsg2)
-		G_free_key_value(proj_epsg2);
-	    
+	    if (proj_srid2)
+		G_free(proj_srid2);
+	    if (proj_wkt2)
+		CPLFree(proj_wkt2);
+
 	    G_warning(_("Projection of layer <%s> is different from "
 			"projection of layer <%s>"),
 			layer_names[layer], layer_names[layer - 1]);
@@ -235,15 +245,19 @@ int cmp_layer_srs(ds_t Ogr_ds, int nlayers, int *layers,
 	    G_free_key_value(proj_info2);
 	if (proj_units2)
 	    G_free_key_value(proj_units2);
-	if (proj_epsg2)
-	    G_free_key_value(proj_epsg2);
+	if (proj_srid2)
+	    G_free(proj_srid2);
+	if (proj_wkt2)
+	    CPLFree(proj_wkt2);
     }
     if (proj_info1)
 	G_free_key_value(proj_info1);
     if (proj_units1)
 	G_free_key_value(proj_units1);
-    if (proj_epsg1)
-	G_free_key_value(proj_epsg1);
+    if (proj_srid1)
+	G_free(proj_srid1);
+    if (proj_wkt1)
+	CPLFree(proj_wkt1);
 
     return 0;
 }
@@ -254,8 +268,11 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
 		      int check_only)
 {
     struct Cell_head loc_wind;
-    struct Key_Value *proj_info, *proj_units, *proj_epsg;
-    struct Key_Value *loc_proj_info, *loc_proj_units;
+    struct Key_Value *proj_info = NULL,
+                     *proj_units = NULL;
+    struct Key_Value *loc_proj_info = NULL,
+                     *loc_proj_units = NULL;
+    char *wkt = NULL, *srid = NULL;
     char error_msg[8096];
     int proj_trouble;
     OGRLayerH Ogr_layer;
@@ -264,13 +281,8 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
     Ogr_layer = ds_getlayerbyindex(hDS, layer);
 
     /* -------------------------------------------------------------------- */
-    /*      Fetch the projection in GRASS form.                             */
+    /*      Fetch the projection in GRASS form, SRID, and WKT.              */
     /* -------------------------------------------------------------------- */
-    proj_info = NULL;
-    proj_units = NULL;
-    proj_epsg = NULL;
-    loc_proj_info = NULL;
-    loc_proj_units = NULL;
 
     /* proj_trouble:
      * 0: valid srs
@@ -280,21 +292,21 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
 
     /* Projection only required for checking so convert non-interactively */
     proj_trouble = get_layer_proj(Ogr_layer, cellhd, &proj_info, &proj_units,
-		                  &proj_epsg, geom_col, 1);
+		                  &srid, &wkt, geom_col, 1);
 
     /* -------------------------------------------------------------------- */
     /*      Do we need to create a new location?                            */
     /* -------------------------------------------------------------------- */
     if (outloc != NULL) {
 	/* do not create a xy location because this can mean that the
-	 * real SRS has not been recognized or is missing */ 
+	 * real SRS has not been recognized or is missing */
 	if (proj_trouble) {
 	    G_fatal_error(_("Unable to convert input map projection to GRASS "
 			    "format; cannot create new location."));
 	}
 	else {
-            if (0 != G_make_location_epsg(outloc, cellhd, proj_info,
-	                                  proj_units, proj_epsg)) {
+            if (0 != G_make_location_crs(outloc, cellhd, proj_info,
+	                                 proj_units, srid, wkt)) {
                 G_fatal_error(_("Unable to create new location <%s>"),
                               outloc);
             }
@@ -490,7 +502,7 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
                        _("Consider generating a new location from the input dataset using "
                          "the 'location' parameter.\n"));
             }
-            
+
 	    if (check_only)
 		msg_fn = G_message;
 	    else
@@ -505,7 +517,7 @@ void check_projection(struct Cell_head *cellhd, ds_t hDS, int layer, char *geom_
 	    if (check_only)
 		msg_fn = G_message;
 	    else
-		msg_fn = G_verbose_message;            
+		msg_fn = G_verbose_message;
 	    msg_fn(_("Projection of input dataset and current location "
 		     "appear to match"));
 	    if (check_only) {

vector/v.out.ogr/Makefile

@@ -6,7 +6,7 @@ EXTRA_CFLAGS = $(GDALCFLAGS)
 
 OBJS = main.o
 
-LIBES = $(GPROJLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(GDALLIBS)
+LIBES = $(GPROJLIB) $(VECTORLIB) $(DBMILIB) $(GISLIB) $(GDALLIBS) $(PROJLIB)
 DEPENDENCIES = $(GPROJDEP) $(VECTORDEP) $(DBMIDEP) $(GISDEP)
 EXTRA_INC = $(VECT_INC) $(PROJINC)
 EXTRA_CFLAGS = $(VECT_CFLAGS)

vector/v.out.ogr/main.c

@@ -5,7 +5,7 @@
  *
  * AUTHOR(S):    Radim Blazek
  *               Some extensions: Markus Neteler, Benjamin Ducke
- *               Multi-features support by Martin Landa <landa.martin gmail.com> 
+ *               Multi-features support by Martin Landa <landa.martin gmail.com>
  *
  * PURPOSE:      Converts GRASS vector to one of supported OGR vector formats.
  *
@@ -35,7 +35,7 @@ int main(int argc, char *argv[])
     int num_to_export;
     int field;
     int overwrite, found;
-    
+
     struct GModule *module;
     struct Options options;
     struct Flags flags;
@@ -75,7 +75,7 @@ int main(int argc, char *argv[])
     int num_types;
     char *dsn;
     int outer_ring_ccw;
-    
+
     G_gisinit(argv[0]);
 
     /* Module options */
@@ -90,7 +90,7 @@ int main(int argc, char *argv[])
 	_("Exports a vector map layer to any of the supported OGR vector formats.");
     module->description = _("By default a vector map layer is exported to OGC GeoPackage format.");
     module->overwrite = TRUE;
-    
+
     /* parse & read options */
     parse_args(argc, argv, &options, &flags);
 
@@ -104,7 +104,7 @@ int main(int argc, char *argv[])
         strcmp(options.format->answer, "PostgreSQL") != 0)
         G_warning(_("Data source starts with \"PG:\" prefix, expecting \"PostgreSQL\" "
                     "format (\"%s\" given)"), options.format->answer);
-    
+
     /* parse dataset creation options */
     i = 0;
     while (options.dsco->answers[i]) {
@@ -268,7 +268,71 @@ int main(int argc, char *argv[])
 	projinfo = G_get_projinfo();
 	projunits = G_get_projunits();
 	projepsg = G_get_projepsg();
-	Ogr_projection = GPJ_grass_to_osr2(projinfo, projunits, projepsg);
+
+#if GDAL_VERSION_MAJOR >= 3 && PROJ_VERSION_MAJOR >= 6
+	char *indef = NULL, *inwkt = NULL;
+
+	if ((indef = G_get_projsrid())) {
+	    PJ *obj = NULL;
+
+	    if ((obj = proj_create(NULL, indef))) {
+		inwkt = G_store(proj_as_wkt(NULL, obj, PJ_WKT2_LATEST, NULL));
+
+		if (inwkt && !*inwkt) {
+		    G_free(inwkt);
+		    inwkt = NULL;
+		}
+	    }
+	}
+	if (!inwkt) {
+	    inwkt = G_get_projwkt();
+	}
+	if (inwkt) {
+	    Ogr_projection = OSRNewSpatialReference(inwkt);
+	}
+#endif
+	if (!Ogr_projection)
+	    Ogr_projection = GPJ_grass_to_osr2(projinfo, projunits, projepsg);
+
+#if GDAL_VERSION_MAJOR >= 3 && PROJ_VERSION_MAJOR >= 6
+	if (Ogr_projection) {
+	    /* convert bound CRS */
+	    PJ *obj = NULL;
+	    char *papszOptions[2];
+
+	    inwkt = NULL;
+
+	    papszOptions[0] = G_store("FORMAT=WKT2");
+	    papszOptions[1] = NULL;
+	    OSRExportToWktEx(Ogr_projection, &inwkt, (const char **)papszOptions);
+	    G_free(papszOptions[0]);
+
+	    if ((obj = proj_create(NULL, inwkt))) {
+		if (proj_get_type(obj) == PJ_TYPE_BOUND_CRS) {
+		    PJ *source_crs;
+
+		    G_debug(1, "found bound crs");
+		    source_crs = proj_get_source_crs(NULL, obj);
+		    if (source_crs) {
+			inwkt = G_store(proj_as_wkt(NULL, source_crs, PJ_WKT2_LATEST, NULL));
+			if (inwkt && !*inwkt) {
+			    G_free(inwkt);
+			    inwkt = NULL;
+			}
+			proj_destroy(source_crs);
+			Ogr_projection = NULL;
+			if (inwkt) {
+			    char *inwkttmp = inwkt;
+			    OSRImportFromWkt(Ogr_projection, &inwkttmp);
+			}
+		    }
+		}
+		proj_destroy(obj);
+	    }
+	    if (inwkt)
+		CPLFree(inwkt);
+	}
+#endif
 
 	if (Ogr_projection ==  NULL)
 	    G_fatal_error(_("Unable to create OGR spatial reference"));
@@ -496,7 +560,7 @@ int main(int argc, char *argv[])
 	    G_debug(1, "Create OGR data source");
 	    hDS = OGR_Dr_CreateDataSource(hDriver, dsn, papszDSCO);
 	}
-#endif	
+#endif
     }
     else {
 #if GDAL_VERSION_NUM >= 2020000
@@ -519,12 +583,12 @@ int main(int argc, char *argv[])
 	}
 #endif
     }
-	
+
     CSLDestroy(papszDSCO);
     if (hDS == NULL)
 	G_fatal_error(_("Unable to open OGR data source '%s'"),
 		      options.dsn->answer);
-    
+
     /* check if OGR layer exists */
     overwrite = G_check_overwrite(argc, argv);
     found = FALSE;
@@ -539,7 +603,7 @@ int main(int argc, char *argv[])
 	Ogr_field = OGR_L_GetLayerDefn(Ogr_layer);
 	if (G_strcasecmp(OGR_FD_GetName(Ogr_field), options.layer->answer))
 	    continue;
-	
+
 	found = TRUE;
 	if (!overwrite && !flags.append->answer) {
 	    G_fatal_error(_("Layer <%s> already exists in OGR data source '%s'"),
@@ -563,11 +627,11 @@ int main(int argc, char *argv[])
 		  options.layer->answer);
 	flags.append->answer = FALSE;
     }
-    
+
     /* Automatically append driver options for 3D output to layer
        creation options if '2' is not given.*/
     if (!flags.force2d->answer &&
-        Vect_is_3d(&In) && 
+        Vect_is_3d(&In) &&
         strcmp(options.format->answer, "ESRI_Shapefile") == 0) {
 	/* find right option */
 	char shape_geom[20];
@@ -592,7 +656,7 @@ int main(int argc, char *argv[])
 	    papszLCO = CSLSetNameValue(papszLCO, "SHPT", shape_geom);
 	}
     }
-    
+
     /* check if the map is 3d */
     if (Vect_is_3d(&In)) {
 	/* specific check for ESRI ShapeFile */
@@ -631,19 +695,19 @@ int main(int argc, char *argv[])
 #if GDAL_VERSION_NUM >= 2020000
     if (flags.append->answer)
 	Ogr_layer = GDALDatasetGetLayerByName(hDS, options.layer->answer);
-    else 
+    else
 	Ogr_layer = GDALDatasetCreateLayer(hDS, options.layer->answer,
 	                                   Ogr_projection, wkbtype,
 					   papszLCO);
 #else
     if (flags.append->answer)
 	Ogr_layer = OGR_DS_GetLayerByName(hDS, options.layer->answer);
-    else 
+    else
 	Ogr_layer = OGR_DS_CreateLayer(hDS, options.layer->answer,
 	                               Ogr_projection, wkbtype,
 				       papszLCO);
 #endif
-    
+
     CSLDestroy(papszLCO);
     if (Ogr_layer == NULL) {
 	if (flags.append->answer)
@@ -651,7 +715,7 @@ int main(int argc, char *argv[])
 	else
 	    G_fatal_error(_("Unable to create OGR layer"));
     }
-    
+
     db_init_string(&dbstring);
 
     /* Vector attributes -> OGR fields */
@@ -668,16 +732,16 @@ int main(int argc, char *argv[])
 		G_warning(_("Exporting 'cat' anyway, as it is the only attribute table field"));
 		flags.nocat->answer = FALSE;
 	    }
-	    
+
 	    if (flags.append->answer) {
 		Ogr_field = OGR_L_GetLayerDefn(Ogr_layer);
 		if (OGR_FD_GetFieldIndex(Ogr_field, GV_KEY_COLUMN) > -1)
 		    create_field = FALSE;
-		else 
+		else
 		    G_warning(_("New attribute column <%s> added to the table"),
 			      GV_KEY_COLUMN);
 	    }
-	    
+
 	    if (create_field) {
 		Ogr_field = OGR_Fld_Create(GV_KEY_COLUMN, OFTInteger);
 		if (OGR_L_CreateField(Ogr_layer, Ogr_field, 0) != OGRERR_NONE)
@@ -685,7 +749,7 @@ int main(int argc, char *argv[])
 		                  GV_KEY_COLUMN);
 		OGR_Fld_Destroy(Ogr_field);
 	    }
-	    
+
 	    doatt = 0;
 	}
 	else {
@@ -711,7 +775,7 @@ int main(int argc, char *argv[])
 		colwidth = db_get_column_length(Column);
 		G_debug(3, "col %d: %s sqltype=%d ctype=%d width=%d",
 			i, colname[i], colsqltype, colctype[i], colwidth);
-		 
+
 		switch (colctype[i]) {
 		case DB_C_TYPE_INT:
 		    ogr_ftype = OFTInteger;
@@ -755,27 +819,30 @@ int main(int argc, char *argv[])
 			G_warning(_("New attribute column <%s> added to the table"),
 				   colname[i]);
 		}
-		 
+
 		Ogr_field = OGR_Fld_Create(colname[i], ogr_ftype);
 		if (ogr_ftype == OFTString && colwidth > 0)
 		    OGR_Fld_SetWidth(Ogr_field, colwidth);
 		if (OGR_L_CreateField(Ogr_layer, Ogr_field, 0) != OGRERR_NONE)
 		    G_fatal_error(_("Unable to create column <%s>"),
 		                  colname[i]);
-		 
+
 		OGR_Fld_Destroy(Ogr_field);
 	    }
 	    if (keycol == -1)
 		G_fatal_error(_("Key column <%s> not found"), Fi->key);
 	}
     }
-    
+
     Ogr_featuredefn = OGR_L_GetLayerDefn(Ogr_layer);
 
     n_feat = n_nocat = n_noatt = 0;
 
-    if (OGR_L_TestCapability(Ogr_layer, OLCTransactions))
-	OGR_L_StartTransaction(Ogr_layer);
+    if (OGR_L_TestCapability(Ogr_layer, OLCTransactions)) {
+	if (OGR_L_StartTransaction(Ogr_layer) != OGRERR_NONE) {
+	    G_warning(_("Unable to start OGR transaction"));
+	}
+    }
 
     /* export polygons oriented according to OGC simple features standard 1.2.1
      * outer rings are oriented counter-clockwise (CCW)
@@ -802,7 +869,7 @@ int main(int argc, char *argv[])
         n_feat += export_lines(&In, field, otype, flags.multi->answer ? TRUE : FALSE,
                                donocat, ftype == GV_BOUNDARY ? TRUE : FALSE,
                                Ogr_featuredefn, Ogr_layer,
-                               Fi, Driver, ncol, colctype, 
+                               Fi, Driver, ncol, colctype,
                                colname, doatt, flags.nocat->answer ? TRUE : FALSE,
                                &n_noatt, &n_nocat);
     }
@@ -814,9 +881,9 @@ int main(int argc, char *argv[])
                      Vect_get_num_areas(&In)),
 		  Vect_get_num_areas(&In));
 
-        n_feat += export_areas(&In, field, flags.multi->answer ? TRUE : FALSE, donocat, 
+        n_feat += export_areas(&In, field, flags.multi->answer ? TRUE : FALSE, donocat,
                                Ogr_featuredefn, Ogr_layer,
-                               Fi, Driver, ncol, colctype, 
+                               Fi, Driver, ncol, colctype,
                                colname, doatt, flags.nocat->answer ? TRUE : FALSE,
                                &n_noatt, &n_nocat, outer_ring_ccw);
     }
@@ -835,8 +902,11 @@ int main(int argc, char *argv[])
 	G_warning(_("Export of volumes not implemented yet. Skipping."));
     }
 
-    if (OGR_L_TestCapability(Ogr_layer, OLCTransactions))
-	OGR_L_CommitTransaction(Ogr_layer);
+    if (OGR_L_TestCapability(Ogr_layer, OLCTransactions)) {
+	if (OGR_L_CommitTransaction(Ogr_layer) != OGRERR_NONE) {
+	    G_warning(_("Unable to commit OGR transaction"));
+	}
+    }
 
     ds_close(hDS);
 

vector/v.proj/main.c

@@ -5,12 +5,12 @@
  * AUTHOR(S):    Irina Kosinovsky, US ARMY CERL,
  *               M.L. Holko, USDA, SCS, NHQ-CGIS,
  *               R.L. Glenn, USDA, SCS, NHQ-CGIS (original contributors
- *               Update to GRASS 6: Radim Blazek <radim.blazek gmail.com> 
+ *               Update to GRASS 6: Radim Blazek <radim.blazek gmail.com>
  *               Huidae Cho <grass4u gmail.com>, Hamish Bowman <hamish_b yahoo.com>,
  *               Jachym Cepicky <jachym les-ejk.cz>, Markus Neteler <neteler itc.it>,
  *               Paul Kelly <paul-grass stjohnspoint.co.uk>
  *               Markus Metz
- * PURPOSE:      
+ * PURPOSE:
  * COPYRIGHT:    (C) 1999-2008, 2018 by the GRASS Development Team
  *
  *               This program is free software under the GNU General Public
@@ -81,7 +81,7 @@ int main(int argc, char *argv[])
     ilocopt->required = YES;
     ilocopt->label = _("Location containing input vector map");
     ilocopt->guisection = _("Source");
-    
+
     isetopt = G_define_standard_option(G_OPT_M_MAPSET);
     isetopt->label = _("Mapset containing input vector map");
     isetopt->description = _("Default: name of current mapset");
@@ -92,10 +92,10 @@ int main(int argc, char *argv[])
     mapopt->label = _("Name of input vector map to re-project");
     mapopt->description = NULL;
     mapopt->guisection = _("Source");
-    
+
     ibaseopt = G_define_standard_option(G_OPT_M_DBASE);
     ibaseopt->label = _("Path to GRASS database of input location");
-    
+
     smax = G_define_option();
     smax->key = "smax";
     smax->type = TYPE_DOUBLE;
@@ -183,9 +183,40 @@ int main(int argc, char *argv[])
     Out_proj = G_projection();
     if (Out_proj == PROJECTION_LL && flag.wrap->answer)
 	nowrap = 1;
-    
+
     G_begin_distance_calculations();
 
+    /****** get the output projection parameters ******/
+    out_proj_keys = G_get_projinfo();
+    if (out_proj_keys == NULL)
+	exit(EXIT_FAILURE);
+
+    out_unit_keys = G_get_projunits();
+    if (out_unit_keys == NULL)
+	exit(EXIT_FAILURE);
+
+    if (pj_get_kv(&info_out, out_proj_keys, out_unit_keys) < 0)
+	exit(EXIT_FAILURE);
+
+    info_out.srid = G_get_projsrid();
+    info_out.wkt = G_get_projwkt();
+
+    if (G_verbose() == G_verbose_max()) {
+	pj_print_proj_params(&info_in, &info_out);
+    }
+
+    info_trans.def = NULL;
+#ifdef HAVE_PROJ_H
+    if (pipeline->answer) {
+	info_trans.def = G_store(pipeline->answer);
+    }
+#endif
+
+    /* Initialize the Point / Cat structure */
+    Points = Vect_new_line_struct();
+    Points2 = Vect_new_line_struct();
+    Cats = Vect_new_cats_struct();
+
     /* Change the location here and then come back */
 
     select_target_env();
@@ -193,7 +224,7 @@ int main(int argc, char *argv[])
     G_setenv_nogisrc("LOCATION_NAME", iloc_name);
     G_setenv_nogisrc("MAPSET", iset_name);
     stat = G_mapset_permissions(iset_name);
-    
+
     if (stat >= 0) {		/* yes, we can access the mapset */
 	/* if requested, list the vector maps in source location - MN 5/2001 */
 	if (flag.list->answer) {
@@ -233,7 +264,7 @@ int main(int argc, char *argv[])
 	/* apparently the +over switch must be set in the input projection,
 	 * not the output latlon projection
 	 * TODO: for PROJ 6+, the +over switch must be added to the
-	 * transformation pipeline if authority:name or WKt are used as
+	 * transformation pipeline if authority:name or WKT are used as
 	 * crs definition */
 	if (Out_proj == PROJECTION_LL && nowrap == 1)
 	    G_set_key_value("over", "defined", in_proj_keys);
@@ -245,14 +276,88 @@ int main(int argc, char *argv[])
 	if (pj_get_kv(&info_in, in_proj_keys, in_unit_keys) < 0)
 	    exit(EXIT_FAILURE);
 
+	info_in.srid = G_get_projsrid();
+	info_in.wkt = G_get_projwkt();
+
 	Vect_set_open_level(1);
 	G_debug(1, "Open old: location: %s mapset : %s", G_location_path(),
 		G_mapset());
 	if (Vect_open_old(&Map, map_name, mapset) < 0)
 	    G_fatal_error(_("Unable to open vector map <%s>"), map_name);
+
+#if PROJ_VERSION_MAJOR >= 6
+	/* need to set the region to the input vector
+	 * for PROJ to select the appropriate pipeline */
+	{
+	    int first = 1, counter = 0;
+	    struct Cell_head inwindow;
+
+	    G_unset_window();
+	    G_get_window(&inwindow);
+
+	    /* Cycle through all lines */
+	    Vect_rewind(&Map);
+	    while (1) {
+		type = Vect_read_next_line(&Map, Points, Cats);	/* read line */
+		if (type == 0)
+		    continue;		/* Dead */
+
+		if (type == -1)
+		    G_fatal_error(_("Reading input vector map"));
+		if (type == -2)
+		    break;
+
+		if (first && Points->n_points > 0) {
+		    first = 0;
+		    inwindow.east = inwindow.west = Points->x[0];
+		    inwindow.north = inwindow.south = Points->y[0];
+		    inwindow.top = inwindow.bottom = Points->z[0];
+		}
+		for (i = 0; i < Points->n_points; i++) {
+		    if (inwindow.east < Points->x[i])
+			inwindow.east = Points->x[i];
+		    if (inwindow.west > Points->x[i])
+			inwindow.west = Points->x[i];
+		    if (inwindow.north < Points->y[i])
+			inwindow.north = Points->y[i];
+		    if (inwindow.south > Points->y[i])
+			inwindow.south = Points->y[i];
+		}
+		counter++;
+	    }
+	    if (counter == 0) {
+		G_warning(_("Input vector map <%s> is empty"), omap_name);
+		exit(EXIT_SUCCESS);
+	    }
+	    inwindow.ns_res = inwindow.ew_res = 1;
+	    /* align to resolution */
+	    inwindow.east = ceil(inwindow.east);
+	    inwindow.west = floor(inwindow.west);
+	    inwindow.north = ceil(inwindow.north);
+	    inwindow.south = floor(inwindow.south);
+	    if (inwindow.east == inwindow.west) {
+		inwindow.east += 0.5;
+		inwindow.west -= 0.5;
+	    }
+	    if (inwindow.north == inwindow.south) {
+		inwindow.north += 0.5;
+		inwindow.south -= 0.5;
+	    }
+	    G_debug(1, "input map north: %.8f", inwindow.north);
+	    G_debug(1, "input map south: %.8f", inwindow.south);
+	    G_debug(1, "input map east: %.8f", inwindow.east);
+	    G_debug(1, "input map west: %.8f", inwindow.west);
+
+	    G_set_window(&inwindow);
+	}
+	/* GPJ_init_transform() must be called only after the region has been set */
+#endif
+	if (GPJ_init_transform(&info_in, &info_out, &info_trans) < 0)
+	    G_fatal_error(_("Unable to initialize coordinate transformation"));
+
     }
-    else if (stat < 0)
-    {				/* allow 0 (i.e. denied permission) */
+    else if (stat < 0) {
+    	/* allow 0 (i.e. denied permission) */
 	/* need to be able to read from others */
 	if (stat == 0)
 	    G_fatal_error(_("Mapset <%s> in input location <%s> - permission denied"),
@@ -264,46 +369,16 @@ int main(int argc, char *argv[])
 
     select_current_env();
 
-    /****** get the output projection parameters ******/
-    out_proj_keys = G_get_projinfo();
-    if (out_proj_keys == NULL)
-	exit(EXIT_FAILURE);
-
-    out_unit_keys = G_get_projunits();
-    if (out_unit_keys == NULL)
-	exit(EXIT_FAILURE);
-
-    if (pj_get_kv(&info_out, out_proj_keys, out_unit_keys) < 0)
-	exit(EXIT_FAILURE);
-
     G_free_key_value(in_proj_keys);
     G_free_key_value(in_unit_keys);
     G_free_key_value(out_proj_keys);
     G_free_key_value(out_unit_keys);
 
-    if (G_verbose() == G_verbose_max()) {
-	pj_print_proj_params(&info_in, &info_out);
-    }
-
-    info_trans.def = NULL;
-#ifdef HAVE_PROJ_H
-    if (pipeline->answer) {
-	info_trans.def = G_store(pipeline->answer);
-    }
-#endif
-    if (GPJ_init_transform(&info_in, &info_out, &info_trans) < 0)
-	G_fatal_error(_("Unable to initialize coordinate transformation"));
-
-    /* Initialize the Point / Cat structure */
-    Points = Vect_new_line_struct();
-    Points2 = Vect_new_line_struct();
-    Cats = Vect_new_cats_struct();
-
     /* test if latlon wrapping to -180,180 should be disabled */
     if (Out_proj == PROJECTION_LL && nowrap == 0) {
 	int first = 1, counter = 0;
 	double x, y;
-	
+
 	/* Cycle through all lines */
 	Vect_rewind(&Map);
 	while (1) {
@@ -315,7 +390,7 @@ int main(int argc, char *argv[])
 		G_fatal_error(_("Reading input vector map"));
 	    if (type == -2)
 		break;
-		
+
 	    if (first && Points->n_points > 0) {
 		first = 0;
 		src_box.E = src_box.W = Points->x[0];
@@ -343,7 +418,7 @@ int main(int argc, char *argv[])
 	y = src_box.N;
 	if (GPJ_transform(&info_in, &info_out, &info_trans, PJ_FWD,
 			  &x, &y, NULL) < 0)
-	    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+	    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 			   "GPJ_transform()");
 
 	tgt_box.E = x;
@@ -355,7 +430,7 @@ int main(int argc, char *argv[])
 	y = src_box.S;
 	if (GPJ_transform(&info_in, &info_out, &info_trans, PJ_FWD,
 			  &x, &y, NULL) < 0)
-	    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+	    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 			   "GPJ_transform()");
 	if (tgt_box.W > x)
 	    tgt_box.W = x;
@@ -370,7 +445,7 @@ int main(int argc, char *argv[])
 	y = src_box.N;
 	if (GPJ_transform(&info_in, &info_out, &info_trans, PJ_FWD,
 			  &x, &y, NULL) < 0)
-	    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+	    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 			   "GPJ_transform()");
 	if (tgt_box.W > x) {
 	    tgt_box.E = x + 360;
@@ -385,7 +460,7 @@ int main(int argc, char *argv[])
 	y = src_box.S;
 	if (GPJ_transform(&info_in, &info_out, &info_trans, PJ_FWD,
 			  &x, &y, NULL) < 0)
-	    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+	    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 			   "GPJ_transform()");
 	if (tgt_box.W > x) {
 	    if (tgt_box.E < x + 360)
@@ -405,7 +480,7 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Unable to create vector map <%s>"), omap_name);
 
     Vect_set_error_handler_io(NULL, &Out_Map); /* register standard i/o error handler */
-    
+
     Vect_copy_head_data(&Map, &Out_Map);
     Vect_hist_copy(&Map, &Out_Map);
     Vect_hist_command(&Out_Map);
@@ -465,12 +540,12 @@ int main(int argc, char *argv[])
 
 		if (GPJ_transform(&info_in, &info_out, &info_trans, PJ_FWD,
 				  &x1, &y1, flag.transformz->answer ? &z1 : NULL) < 0)
-		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 				   "GPJ_transform()");
 
 		if (GPJ_transform(&info_in, &info_out, &info_trans, PJ_FWD,
 				  &x2, &y2, flag.transformz->answer ? &z2 : NULL) < 0)
-		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"), 
+		    G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
 				   "GPJ_transform()");
 
 		Vect_append_point(Points2, x1, y1, z1);