libgis: add read/write functions for WKT and srid (#976)

GRASS native projection information is reaching its limits in representing CRS's, while PROJ and GDAL have been substantially enhanced.Use WKT or authority code if possible to convert to different formats

+ libgis: read/write functions for WKT and spatial reference id as recognized by PROJ

+ g.proj: use WKT or authority code if possible

general/g.proj/create.c

@@ -10,8 +10,8 @@ void create_location(const char *location)
 {
     int ret;
 
-    ret = G_make_location_epsg(location, &cellhd, projinfo, projunits,
-                               projepsg);
+    ret = G_make_location_crs(location, &cellhd, projinfo, projunits,
+                               projepsg, projwkt, projsrid);
     if (ret == 0)
 	G_message(_("Location <%s> created"), location);
     else if (ret == -1)
@@ -59,6 +59,14 @@ void modify_projinfo()
 	G_write_key_value_file(path, projepsg);
     }
 
+    if (projwkt != NULL) {
+	G_write_projwkt(NULL, projwkt);
+    }
+
+    if (projsrid != NULL) {
+	G_write_projsrid(NULL, projsrid);
+    }
+
     if ((old_cellhd.zone != cellhd.zone) ||
 	(old_cellhd.proj != cellhd.proj)) {
 	/* Recreate the default, and current window files if projection

general/g.proj/input.c

@@ -34,10 +34,6 @@
 
 static void set_default_region(void);
 
-#ifdef HAVE_OGR
-static void set_gdal_region(GDALDatasetH);
-#endif
-
 /**
  * \brief Read projection and region information from current location
  * 
@@ -49,6 +45,8 @@ void input_currloc(void)
 {
     G_get_default_window(&cellhd);
     if (cellhd.proj != PROJECTION_XY) {
+	projsrid = G_get_projsrid();
+	projwkt = G_get_projwkt();
 	projinfo = G_get_projinfo();
 	projunits = G_get_projunits();
         projepsg = G_get_projepsg();
@@ -57,7 +55,38 @@ void input_currloc(void)
     return;
 }
 
+/**
+ * \brief Populates global cellhd with "default" region settings
+ **/
+
+static void set_default_region(void)
+{
+    /* If importing projection there will be no region information, so
+     * set some default values */
+    cellhd.rows = 1;
+    cellhd.rows3 = 1;
+    cellhd.cols = 1;
+    cellhd.cols3 = 1;
+    cellhd.depths = 1.;
+    cellhd.north = 1.;
+    cellhd.ns_res = 1.;
+    cellhd.ns_res3 = 1.;
+    cellhd.south = 0.;
+    cellhd.west = 0.;
+    cellhd.ew_res = 1.;
+    cellhd.ew_res3 = 1.;
+    cellhd.east = 1.;
+    cellhd.top = 1.;
+    cellhd.tb_res = 1.;
+    cellhd.bottom = 0.;
+
+    return;
+}
+
+
 #ifdef HAVE_OGR
+static void set_gdal_region(GDALDatasetH);
+static void set_authnamecode(OGRSpatialReferenceH);
 
 /**
  * \brief Read projection information in WKT format from stdin or a file
@@ -77,8 +106,9 @@ void input_currloc(void)
 int input_wkt(char *wktfile)
 {
     FILE *infd;
-    char buff[8000];
+    char buff[8000], *tmpwkt;
     OGRSpatialReferenceH hSRS;
+    const char *papszOptions[3];
     int ret;
 
     if (strcmp(wktfile, "-") == 0)
@@ -89,7 +119,7 @@ int input_wkt(char *wktfile)
     if (infd) {
 	fread(buff, sizeof(buff), 1, infd);
 	if (ferror(infd))
-	    G_fatal_error(_("Error reading WKT projection description"));
+	    G_fatal_error(_("Error reading WKT definition"));
 	else
 	    fclose(infd);
 	/* Get rid of newlines */
@@ -98,30 +128,67 @@ int input_wkt(char *wktfile)
     else
 	G_fatal_error(_("Unable to open file '%s' for reading"), wktfile);
 
+    projwkt = G_store(buff);
+
+#if PROJ_VERSION_MAJOR >= 6
+    /* validate input WKT */
+    {
+	PROJ_STRING_LIST wkt_warnings, wkt_grammar_errors;
+	PJ *obj;
+
+	wkt_warnings = NULL;
+	wkt_grammar_errors = NULL;
+
+	/* no strict validation */
+	obj = proj_create_from_wkt(NULL, buff, NULL, &wkt_warnings,
+	                           &wkt_grammar_errors);
+
+	if (wkt_warnings) {
+	    int i;
+
+	    G_warning(_("WKT validation warnings:"));
+	    for (i = 0; wkt_warnings[i]; i++)
+		G_warning("%s", wkt_warnings[i]);
+
+	    proj_string_list_destroy(wkt_warnings);
+	}
+
+	if (wkt_grammar_errors) {
+	    int i;
+
+	    G_warning(_("WKT validation grammar errors:"));
+	    for (i = 0; wkt_grammar_errors[i]; i++)
+		G_warning("%s", wkt_grammar_errors[i]);
+
+	    proj_string_list_destroy(wkt_grammar_errors);
+	}
+	proj_destroy(obj);
+    }
+#endif
+
+    /* get GRASS proj info + units */
+    /* NOTE: GPJ_wkt_to_grass() converts any WKT version to WKT1 */
     ret = GPJ_wkt_to_grass(&cellhd, &projinfo, &projunits, buff, 0);
+    if (ret < 2)
+	G_fatal_error(_("WKT not recognized: %s"), buff);
+
     set_default_region();
 
+    /* find authname and authcode */
     hSRS = OSRNewSpatialReference(buff);
-    if (hSRS) {
-	const char *authkey, *authname, *authcode;
-
-	authkey = NULL;
-	if (OSRIsProjected(hSRS))
-	    authkey = "PROJCS";
-	else if (OSRIsGeographic(hSRS))
-	    authkey = "GEOGCS";
-
-	if (authkey) {
-	    authname = OSRGetAuthorityName(hSRS, authkey);
-	    if (authname && *authname && strcmp(authname, "EPSG") == 0) {
-		authcode = OSRGetAuthorityCode(hSRS, authkey);
-		if (authcode && *authcode) {
-		    projepsg = G_create_key_value();
-		    G_set_key_value("epsg", authcode, projepsg);
-		}
-	    }
-	}
-    }
+    /* get clean WKT definition */
+    papszOptions[0] = G_store("MULTILINE=YES");
+    papszOptions[1] = G_store("FORMAT=WKT2");
+    papszOptions[2] = NULL;
+#if GDAL_VERSION_MAJOR >= 3
+    OSRExportToWktEx(hSRS, &tmpwkt, papszOptions);
+#else
+    OSRExportToPrettyWkt(hSRS, &tmpwkt, FALSE);
+#endif
+    projwkt = G_store(tmpwkt);
+    CPLFree(tmpwkt);
+    set_authnamecode(hSRS);
+    OSRDestroySpatialReference(hSRS);
 
     return ret;
 }
@@ -150,13 +217,23 @@ int input_proj4(char *proj4params)
     OGRSpatialReferenceH hSRS;
     int ret = 0;
 
+    /* TEST: use PROJ proj_create(), convert to WKT,
+     *       OSRImportFromWkt  */
     if (strcmp(proj4params, "-") == 0) {
 	infd = stdin;
 	fgets(buff, sizeof(buff), infd);
-	G_asprintf(&proj4string, "%s +no_defs", buff);
     }
     else
-	G_asprintf(&proj4string, "%s +no_defs", proj4params);
+	strcpy(buff, proj4params);
+
+#if PROJ_VERSION_MAJOR >= 6
+	if (!strstr(buff, "+type=crs"))
+	    G_asprintf(&proj4string, "%s +no_defs +type=crs", buff);
+	else
+	    G_asprintf(&proj4string, "%s +no_defs", buff);
+#else
+	G_asprintf(&proj4string, "%s +no_defs", buff);
+#endif
 
     /* Set finder function for locating OGR csv co-ordinate system tables */
     /* SetCSVFilenameHook(GPJ_set_csv_loc); */
@@ -169,6 +246,8 @@ int input_proj4(char *proj4params)
 
     ret = GPJ_osr_to_grass(&cellhd, &projinfo, &projunits, hSRS, 0);
 
+    /* authority name and code not available in PROJ definitions */
+
     OSRDestroySpatialReference(hSRS);
 
     set_default_region();
@@ -177,6 +256,67 @@ int input_proj4(char *proj4params)
 }
 
 /**
+ * \brief Read projection information corresponding to a spatial
+ *        reference id (srid)
+ *
+ * Determines projection information corresponding to a srid
+ * composed of authority name and code and stores in global structs
+ * projinfo and projunits. Populates global cellhd with default region
+ * information.
+ *
+ * Examples: "EPSG:4326", "urn:ogc:def:crs:EPSG::4326"
+ *
+ * \param srid    spatial reference id
+ *
+ * \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
+ **/
+
+int input_srid(char *srid)
+{
+#if PROJ_VERSION_MAJOR >= 6
+    OGRSpatialReferenceH hSRS;
+    int ret = 0;
+    const char *papszOptions[3];
+    PJ *obj;
+    const char *tmpwkt;
+
+    /* GDAL alternative: OSRSetFromUserInput() */
+    obj = proj_create(NULL, srid);
+    if (!obj)
+	G_fatal_error(_("SRID <%s> not recognized by PROJ"), srid);
+
+    tmpwkt = proj_as_wkt(NULL, obj, PJ_WKT2_2019, NULL);
+    hSRS = OSRNewSpatialReference(tmpwkt);
+    if (!hSRS)
+	G_fatal_error(_("WKT for SRID <%s> not recognized by GDAL"), srid);
+
+    projsrid = G_store(srid);
+
+    /* WKT */
+    papszOptions[0] = G_store("MULTILINE=YES");
+    papszOptions[1] = G_store("FORMAT=WKT2");
+    papszOptions[2] = NULL;
+    if (OSRExportToWktEx(hSRS, &projwkt, papszOptions) != OGRERR_NONE)
+	G_warning(_("Unable to convert srid to WKT"));
+
+    /* GRASS proj info + units */
+    ret = GPJ_osr_to_grass(&cellhd, &projinfo, &projunits, hSRS, 0);
+
+    proj_destroy(obj);
+    OSRDestroySpatialReference(hSRS);
+
+    set_default_region();
+    return ret;
+#else
+    G_fatal_error(_("Input srid requires GDAL 3+ and PROJ 6+"));
+
+    return 0;
+#endif
+}
+
+/**
  * \brief Read projection information corresponding to an EPSG co-ordinate 
  *        system number
  * 
@@ -196,6 +336,7 @@ int input_epsg(int epsg_num)
     OGRSpatialReferenceH hSRS;
     char epsgstr[100];
     int ret = 0;
+    const char *papszOptions[3];
 
     /* Set finder function for locating OGR csv co-ordinate system tables */
     /* SetCSVFilenameHook(GPJ_set_csv_loc); */
@@ -204,11 +345,24 @@ int input_epsg(int epsg_num)
     if (OSRImportFromEPSG(hSRS, epsg_num) != OGRERR_NONE)
 	G_fatal_error(_("Unable to translate EPSG code"));
 
+    /* GRASS proj info + units */
     ret = GPJ_osr_to_grass(&cellhd, &projinfo, &projunits, hSRS, 0);
 
+    /* EPSG code */
     sprintf(epsgstr, "%d", epsg_num);
     projepsg = G_create_key_value();
     G_set_key_value("epsg", epsgstr, projepsg);
+    /* srid as AUTHORITY:CODE */
+    G_asprintf(&projsrid, "EPSG:%s", epsgstr);
+
+#if GDAL_VERSION_MAJOR >= 3
+    /* WKT */
+    papszOptions[0] = G_store("MULTILINE=YES");
+    papszOptions[1] = G_store("FORMAT=WKT2");
+    papszOptions[2] = NULL;
+    if (OSRExportToWktEx(hSRS, &projwkt, papszOptions) != OGRERR_NONE)
+	G_warning(_("Unable to convert EPSG code to WKT"));
+#endif
 
     OSRDestroySpatialReference(hSRS);
 
@@ -249,15 +403,24 @@ int input_georef(char *geofile)
 
     ogr_ds = NULL;
     hSRS = NULL;
+    /* Try opening with OGR */
     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 {
@@ -272,12 +435,18 @@ int input_georef(char *geofile)
 
 	    G_debug(1, "...succeeded.");
 	    wktstring = (char *)GDALGetProjectionRef(gdal_ds);
+	    projwkt = G_store(wktstring);
 	    ret =
-		GPJ_wkt_to_grass(&cellhd, &projinfo, &projunits, wktstring,
+		GPJ_wkt_to_grass(&cellhd, &projinfo, &projunits, projwkt,
 				 0);
 
 	    set_gdal_region(gdal_ds);
-	    hSRS = OSRNewSpatialReference(wktstring);
+#if GDAL_VERSION_MAJOR >= 3
+	    hSRS = OSRClone(GDALGetSpatialRef(gdal_ds));
+#else
+	    hSRS = OSRNewSpatialReference(projwkt);
+#endif
+	    GDALClose(gdal_ds);
 	}
 	else {
 	    int namelen;
@@ -301,63 +470,14 @@ int input_georef(char *geofile)
 		    "projection information. 'XY (unprojected)' will be used"),
 		  geofile);
 
-    if (hSRS) {
-	const char *authkey, *authname, *authcode;
-
-	authkey = NULL;
-	if (OSRIsProjected(hSRS))
-	    authkey = "PROJCS";
-	else if (OSRIsGeographic(hSRS))
-	    authkey = "GEOGCS";
+    set_authnamecode(hSRS);
+    OSRDestroySpatialReference(hSRS);
 
-	if (authkey) {
-	    authname = OSRGetAuthorityName(hSRS, authkey);
-	    if (authname && *authname && strcmp(authname, "EPSG") == 0) {
-		authcode = OSRGetAuthorityCode(hSRS, authkey);
-		if (authcode && *authcode) {
-		    projepsg = G_create_key_value();
-		    G_set_key_value("epsg", authcode, projepsg);
-		}
-	    }
-	}
-    }
     if (ogr_ds)
 	OGR_DS_Destroy(ogr_ds);
 
     return ret;
-
 }
-#endif /* HAVE_OGR */
-
-/**
- * \brief Populates global cellhd with "default" region settings
- **/
-
-static void set_default_region(void)
-{
-    /* If importing projection there will be no region information, so
-     * set some default values */
-    cellhd.rows = 1;
-    cellhd.rows3 = 1;
-    cellhd.cols = 1;
-    cellhd.cols3 = 1;
-    cellhd.depths = 1.;
-    cellhd.north = 1.;
-    cellhd.ns_res = 1.;
-    cellhd.ns_res3 = 1.;
-    cellhd.south = 0.;
-    cellhd.west = 0.;
-    cellhd.ew_res = 1.;
-    cellhd.ew_res3 = 1.;
-    cellhd.east = 1.;
-    cellhd.top = 1.;
-    cellhd.tb_res = 1.;
-    cellhd.bottom = 0.;
-
-    return;
-}
-
-#ifdef HAVE_OGR
 
 /**
  * \brief Populates global cellhd with region settings based on 
@@ -408,4 +528,33 @@ static void set_gdal_region(GDALDatasetH hDS)
     return;
 }
 
+void set_authnamecode(OGRSpatialReferenceH hSRS)
+{
+    const char *authkey, *authname, *authcode;
+
+    if (!hSRS)
+	return;
+
+    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(&projsrid, "%s:%s", authname, authcode);
+		/* for backwards compatibility, remove ? */
+		if (strcmp(authname, "EPSG") == 0) {
+		    projepsg = G_create_key_value();
+		    G_set_key_value("epsg", authcode, projepsg);
+		}
+	    }
+	}
+    }
+}
+
 #endif /* HAVE_OGR */

general/g.proj/local_proto.h

@@ -1,12 +1,14 @@
 #include <grass/config.h>
 
 extern struct Key_Value *projinfo, *projunits, *projepsg;
+extern char *projsrid, *projwkt;
 extern struct Cell_head cellhd;
 
 /* input.c */
 void input_currloc(void);
 #ifdef HAVE_OGR
 int input_wkt(char *);
+int input_srid(char *);
 int input_proj4(char *);
 int input_epsg(int);
 int input_georef(char *);

general/g.proj/main.c

@@ -23,11 +23,16 @@
 
 #include "local_proto.h"
 
-struct Key_Value *projinfo, *projunits, *projepsg;
+struct Key_Value *projinfo = NULL,
+                 *projunits = NULL,
+		 *projepsg = NULL;
+char *projsrid = NULL, *projwkt = NULL;
 struct Cell_head cellhd;
 
 int main(int argc, char *argv[])
 {
+    /* TODO: replace most of these flags with an option to select the
+     * output format */
     struct Flag *printinfo,	/* Print contents of PROJ_INFO & PROJ_UNITS */
 	*shellinfo,             /* Print in shell script style              */
 	*printproj4,		/* Print projection in PROJ.4 format        */
@@ -42,6 +47,7 @@ int main(int argc, char *argv[])
     
     struct Option *location,	/* Name of new location to create           */
 #ifdef HAVE_OGR
+	*insrid,		/* spatial reference id (auth name + code   */
 	*inepsg,		/* EPSG projection code                     */
 	*inwkt,			/* Input file with projection in WKT format */
 	*inproj4,		/* Projection in PROJ.4 format              */
@@ -141,6 +147,15 @@ int main(int argc, char *argv[])
 		     "description");
     inwkt->description = _("'-' for standard input");
 
+    insrid = G_define_option();
+    insrid->key = "srid";
+    insrid->type = TYPE_STRING;
+    insrid->key_desc = "params";
+    insrid->required = NO;
+    insrid->guisection = _("Specification");
+    insrid->label = _("Spatial reference ID with authority name and code");
+    insrid->description = _("E.g. EPSG:4326 or urn:ogc:def:crs:EPSG::4326");
+
     inproj4 = G_define_option();
     inproj4->key = "proj4";
     inproj4->type = TYPE_STRING;
@@ -228,10 +243,12 @@ int main(int argc, char *argv[])
 	printwkt->answer = 1;
 
     formats = ((ingeo->answer ? 1 : 0) + (inwkt->answer ? 1 : 0) +
-	       (inproj4->answer ? 1 : 0) + (inepsg->answer ? 1 : 0));
+	       (inproj4->answer ? 1 : 0) + (inepsg->answer ? 1 : 0) +
+	       (insrid->answer ? 1 : 0));
     if (formats > 1)
-	G_fatal_error(_("Only one of '%s', '%s', '%s' or '%s' options may be specified"),
-		      ingeo->key, inwkt->key, inproj4->key, inepsg->key);
+	G_fatal_error(_("Only one of '%s', '%s', '%s', '%s' or '%s' options may be specified"),
+		      ingeo->key, inwkt->key, inproj4->key, inepsg->key,
+		      insrid->key);
 
     /* List supported datums if requested; code originally 
      * from G_ask_datum_name() (formerly in libgis) */
@@ -249,6 +266,7 @@ int main(int argc, char *argv[])
 
     epsg = inepsg->answer;
     projinfo = projunits = projepsg = NULL;
+    projsrid = projwkt = NULL;
 
     /* Input */
     /* We can only have one input source, hence if..else construct */
@@ -261,6 +279,9 @@ int main(int argc, char *argv[])
     else if (inwkt->answer)
 	/* Input in WKT format */
 	input_wkt(inwkt->answer);
+    else if (insrid->answer)
+	/* Input as spatial reference ID */
+	input_srid(insrid->answer);
     else if (inproj4->answer)
 	/* Input in PROJ.4 format */
 	input_proj4(inproj4->answer);
@@ -279,13 +300,12 @@ int main(int argc, char *argv[])
 	G_fatal_error(_("Projection files missing"));
 
     /* Override input datum if requested */
-    if(datum->answer)
+    if (datum->answer)
 	set_datum(datum->answer);
 
     /* Set Datum Parameters if necessary or requested */
     set_datumtrans(atoi(dtrans->answer), forcedatumtrans->answer);
 
-
     /* Output */
     /* Only allow one output format at a time, to reduce confusion */
     formats = ((printinfo->answer ? 1 : 0) + (shellinfo->answer ? 1 : 0) +

general/g.proj/output.c

@@ -27,7 +27,8 @@
 
 static int check_xy(int shell);
 
-/* print projection information gathered from one of the possible inputs */
+/* print projection information gathered from one of the possible inputs
+ * in GRASS format */
 void print_projinfo(int shell)
 {
     int i;
@@ -123,43 +124,75 @@ void print_datuminfo(void)
     return;
 }
 
+/* print input projection information in PROJ format */
 void print_proj4(int dontprettify)
 {
     struct pj_info pjinfo;
-    char *proj4, *proj4mod, *i;
+    char *i, *projstrmod;
+    const char *projstr;
     const char *unfact;
 
     if (check_xy(FALSE))
 	return;
 
-    if (pj_get_kv(&pjinfo, projinfo, projunits) == -1)
-        G_fatal_error(_("Unable to convert projection information to PROJ format"));
-    proj4 = pjinfo.def;
-#ifdef HAVE_PROJ_H
-    proj_destroy(pjinfo.pj);
+    projstr = NULL;
+    projstrmod = NULL;
+
+#if PROJ_VERSION_MAJOR >= 6
+    /* PROJ6+: create a PJ object from wkt or srid,
+     * then get PROJ string using PROJ API */
+    {
+	PJ *obj = NULL;
+
+	if (projwkt) {
+	    obj = proj_create_from_wkt(NULL, projwkt, NULL, NULL, NULL);
+	}
+	if (!obj && projsrid) {
+	    obj = proj_create(NULL, projsrid);
+	}
+	if (obj) {
+	    projstr = proj_as_proj_string(NULL, obj, PJ_PROJ_5, NULL);
+
+	    if (projstr)
+		projstr = G_store(projstr);
+	    proj_destroy(obj);
+	}
+    }
+#endif
+
+    if (!projstr) {
+	if (pj_get_kv(&pjinfo, projinfo, projunits) == -1)
+	    G_fatal_error(_("Unable to convert projection information to PROJ format"));
+	projstr = pjinfo.def;
+#if PROJ_VERSION_MAJOR >= 5
+	proj_destroy(pjinfo.pj);
 #else
-    pj_free(pjinfo.pj);
+	pj_free(pjinfo.pj);
 #endif
-    /* GRASS-style PROJ.4 strings don't include a unit factor as this is
-     * handled separately in GRASS - must include it here though */
-    unfact = G_find_key_value("meters", projunits);
-    if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
-	G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
-    else
-	proj4mod = G_store(proj4);
 
-    for (i = proj4mod; *i; i++) {
+	/* GRASS-style PROJ.4 strings don't include a unit factor as this is
+	 * handled separately in GRASS - must include it here though */
+	unfact = G_find_key_value("meters", projunits);
+	if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
+	    G_asprintf(&projstrmod, "%s +to_meter=%s", projstr, unfact);
+	else
+	    projstrmod = G_store(projstr);
+    }
+    if (!projstrmod)
+	projstrmod = G_store(projstr);
+
+    for (i = projstrmod; *i; i++) {
 	/* Don't print the first space */
-	if (i == proj4mod && *i == ' ')
+	if (i == projstrmod && *i == ' ')
 	    continue;
 
-	if (*i == ' ' && *(i+1) == '+' && !(dontprettify))
+	if (*i == ' ' && *(i + 1) == '+' && !(dontprettify))
 	    fputc('\n', stdout);
 	else
 	    fputc(*i, stdout);
     }
     fputc('\n', stdout);
-    G_free(proj4mod);
+    G_free(projstrmod);
 
     return;
 }
@@ -172,8 +205,65 @@ void print_wkt(int esristyle, int dontprettify)
     if (check_xy(FALSE))
 	return;
 
+    outwkt = NULL;
+
+#if PROJ_VERSION_MAJOR >= 6
+    /* print WKT2 using GDAL OSR interface */
+    {
+	OGRSpatialReferenceH hSRS;
+	const char *tmpwkt;
+	const char **papszOptions = NULL;
+
+	papszOptions = G_calloc(3, sizeof(char *));
+	if (dontprettify)
+	    papszOptions[0] = G_store("MULTILINE=NO");
+	else
+	    papszOptions[0] = G_store("MULTILINE=YES");
+	if (esristyle)
+	    papszOptions[1] = G_store("FORMAT=WKT1_ESRI");
+	else
+	    papszOptions[1] = G_store("FORMAT=WKT2");
+	papszOptions[2] = NULL;
+
+	hSRS = NULL;
+
+	if (projsrid) {
+	    PJ *obj;
+
+	    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);
+	    hSRS = OSRNewSpatialReference(tmpwkt);
+	    OSRExportToWktEx(hSRS, &outwkt, papszOptions);
+	}
+	if (!outwkt && projwkt) {
+	    hSRS = OSRNewSpatialReference(projwkt);
+	    OSRExportToWktEx(hSRS, &outwkt, papszOptions);
+	}
+	if (!outwkt && projepsg) {
+	    int epsg_num;
+
+	    epsg_num = atoi(G_find_key_value("epsg", projepsg));
+
+	    OSRImportFromEPSG(hSRS, epsg_num);
+	    OSRExportToWktEx(hSRS, &outwkt, 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);
+	}
+	if (hSRS)
+	    OSRDestroySpatialReference(hSRS);
+    }
+#else
     outwkt = GPJ_grass_to_wkt2(projinfo, projunits, projepsg, esristyle,
 			      !(dontprettify));
+#endif
+
     if (outwkt != NULL) {
 	fprintf(stdout, "%s\n", outwkt);
 	G_free(outwkt);

include/defs/gis.h

@@ -398,6 +398,8 @@ int G_read_ellipsoid_table(int);
 struct Key_Value *G_get_projunits(void);
 struct Key_Value *G_get_projinfo(void);
 struct Key_Value *G_get_projepsg(void);
+char *G_get_projwkt(void);
+char *G_get_projsrid(void);
 
 /* get_window.c */
 void G_get_window(struct Cell_head *);
@@ -522,6 +524,11 @@ 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 *,
+			const struct Key_Value *, const struct Key_Value *,
+			const char *, const char *);
+int G_write_projsrid(const char *, const char *);
+int G_write_projwkt(const char *, const char *);
 int G_compare_projections(const struct Key_Value *, const struct Key_Value *,
 			  const struct Key_Value *, const struct Key_Value *);
 

include/gis.h

@@ -117,6 +117,8 @@ static const char *GRASS_copyright __attribute__ ((unused))
 #define PROJECTION_FILE "PROJ_INFO"
 #define UNIT_FILE       "PROJ_UNITS"
 #define EPSG_FILE       "PROJ_EPSG"
+#define WKT_FILE        "PROJ_WKT"
+#define SRID_FILE       "PROJ_SRID"
 
 #ifdef __MINGW32__
 #define CONFIG_DIR "GRASS7"

include/gprojects.h

@@ -27,9 +27,14 @@
 #include <proj_api.h>
 #define PJ_FWD 	 1
 #define PJ_INV 	-1
+/* PROJ_VERSION_MAJOR is not set in the old PROJ API */
+#define PROJ_VERSION_MAJOR 4
 #endif
 #ifdef HAVE_OGR
 #    include <ogr_srs_api.h>
+#    if PROJ_VERSION_MAJOR >= 6 && GDAL_VERSION_MAJOR < 3
+#        error "PROJ 6+ requires GDAL 3+"
+#    endif
 #endif
 
 /* Data Files */
@@ -60,6 +65,7 @@ struct gpj_datum
     double dx, dy, dz;
 };
 
+/* no longer needed with PROJ6+ */
 struct gpj_datum_transform_list
 {
 

lib/gis/get_projinfo.c

@@ -9,6 +9,8 @@
   (>=v2). Read the file COPYING that comes with GRASS for details.
 */
 
+#include <string.h>
+#include <errno.h>
 #include <unistd.h>
 #include <stdio.h>
 #include <grass/gis.h>
@@ -93,6 +95,8 @@ struct Key_Value *G_get_projinfo(void)
   \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 */ 
 struct Key_Value *G_get_projepsg(void)
 {
     struct Key_Value *in_epsg_keys;
@@ -110,3 +114,183 @@ struct Key_Value *G_get_projepsg(void)
 
     return in_epsg_keys;
 }
+
+/*!
+  \brief Get WKT information for the current location
+
+  \return pointer to WKT string
+  \return NULL when WKT is not available for the current location
+*/
+
+char *G_get_projwkt(void)
+{
+    char *wktstring = NULL;
+    char path[GPATH_MAX];
+    FILE *fp;
+    int n, nalloc;
+    int c;
+
+    G_file_name(path, "", WKT_FILE, "PERMANENT");
+    if (access(path, 0) != 0) {
+	if (G_projection() != PROJECTION_XY) {
+	    G_debug(1, "<%s> file not found for location <%s>",
+		      WKT_FILE, G_location());
+	}
+	return NULL;
+    }
+
+    fp = fopen(path, "r");
+    if (!fp)
+	G_fatal_error(_("Unable to open input file <%s>: %s"), path, strerror(errno));
+
+    wktstring = G_malloc(1024 * sizeof(char));
+    nalloc = 1024;
+
+    n = 0;
+    while (1) {
+	c = fgetc(fp);
+
+	if (c == EOF) {
+	    break;
+	}
+
+	if (c == '\r') {	/* DOS or MacOS9 */
+	    c = fgetc(fp);
+	    if (c != EOF) {
+		if (c != '\n') {	/* MacOS9 - we have to return the char to stream */
+		    ungetc(c, fp);
+		    c = '\n';
+		}
+	    }
+	    else {	/* MacOS9 - we have to return the char to stream */
+		ungetc(c, fp);
+		c = '\n';
+	    }
+	}
+
+	if (n == nalloc) {
+	    wktstring = G_realloc(wktstring, nalloc + 1024);
+	    nalloc += 1024;
+	}
+
+	wktstring[n] = c;
+
+	n++;
+    }
+
+    if (n > 0) {
+	if (n == nalloc) {
+	    wktstring = G_realloc(wktstring, nalloc + 1);
+	    nalloc += 1;
+	}
+	wktstring[n] = '\0';
+    }
+    else {
+	G_free(wktstring);
+	wktstring = NULL;
+    }
+
+    if (fclose(fp) != 0)
+	G_fatal_error(_("Error closing output file <%s>: %s"), path, strerror(errno));
+
+    return wktstring;
+}
+
+/*!
+  \brief Get srid (spatial reference id) for the current location
+
+  Typically an srid will be of the form authority NAME:CODE,
+  e.g. EPSG:4326
+
+  This srid is passed to proj_create() using PROJ or
+  OSRSetFromUserInput() using GDAL. Therefore various other forms of
+  srid are possible, e.g. in OSRSetFromUserInput():
+
+   1. Well Known Text definition - passed on to importFromWkt().
+   2. "EPSG:n" - number passed on to importFromEPSG().
+   3. "EPSGA:n" - number passed on to importFromEPSGA().
+   4. "AUTO:proj_id,unit_id,lon0,lat0" - WMS auto projections.
+   5. "urn:ogc:def:crs:EPSG::n" - ogc urns
+   6. PROJ.4 definitions - passed on to importFromProj4().
+   7. filename - file read for WKT, XML or PROJ.4 definition.
+   8. well known name accepted by SetWellKnownGeogCS(), such as NAD27, NAD83, WGS84 or WGS72.
+   9. "IGNF:xxxx", "ESRI:xxxx", etc. from definitions from the PROJ database;
+  10. PROJJSON (PROJ >= 6.2)
+
+  \return pointer to srid string
+  \return NULL when srid is not available for the current location
+*/
+
+char *G_get_projsrid(void)
+{
+    char *sridstring = NULL;
+    char path[GPATH_MAX];
+    FILE *fp;
+    int n, nalloc;
+    int c;
+
+    G_file_name(path, "", SRID_FILE, "PERMANENT");
+    if (access(path, 0) != 0) {
+	if (G_projection() != PROJECTION_XY) {
+	    G_debug(1, "<%s> file not found for location <%s>",
+		      SRID_FILE, G_location());
+	}
+	return NULL;
+    }
+
+    fp = fopen(path, "r");
+    if (!fp)
+	G_fatal_error(_("Unable to open input file <%s>: %s"), path, strerror(errno));
+
+    sridstring = G_malloc(1024 * sizeof(char));
+    nalloc = 1024;
+
+    n = 0;
+    while (1) {
+	c = fgetc(fp);
+
+	if (c == EOF) {
+	    break;
+	}
+
+	if (c == '\r') {	/* DOS or MacOS9 */
+	    c = fgetc(fp);
+	    if (c != EOF) {
+		if (c != '\n') {	/* MacOS9 - we have to return the char to stream */
+		    ungetc(c, fp);
+		    c = '\n';
+		}
+	    }
+	    else {	/* MacOS9 - we have to return the char to stream */
+		ungetc(c, fp);
+		c = '\n';
+	    }
+	}
+
+	if (n == nalloc) {
+	    sridstring = G_realloc(sridstring, nalloc + 1024);
+	    nalloc += 1024;
+	}
+
+	sridstring[n] = c;
+
+	n++;
+    }
+
+    if (n > 0) {
+	if (n == nalloc) {
+	    sridstring = G_realloc(sridstring, nalloc + 1);
+	    nalloc += 1;
+	}
+	sridstring[n] = '\0';
+    }
+    else {
+	G_free(sridstring);
+	sridstring = NULL;
+    }
+
+    if (fclose(fp) != 0)
+	G_fatal_error(_("Error closing output file <%s>: %s"), path, strerror(errno));
+
+    return sridstring;
+}

lib/gis/make_loc.c

@@ -18,9 +18,11 @@
 
 #include <stdlib.h>
 #include <string.h>
+#include <errno.h>
 #include <unistd.h>
 #include <sys/stat.h>
 #include <math.h>
+#include <grass/glocale.h>
 
 /*!
  * \brief Create a new location
@@ -146,6 +148,77 @@ 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
+ *
+ * \param location_name Name of the new location. Should not include
+ *                      the full path, the location will be created within
+ *                      the current database.
+ * \param wind          default window setting for the new location.
+ *                      All fields should be set in this
+ *                      structure, and care should be taken to ensure that
+ *                      the proj/zone fields match the definition in the
+ *                      proj_info parameter(see G_set_cellhd_from_projinfo()).
+ *
+ * \param proj_info     projection definition suitable to write to the
+ *                      PROJ_INFO file, or NULL for PROJECTION_XY.
+ *
+ * \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.
+ *
+ * \param proj_wkt      WKT defintion suitable to write to the PROJ_WKT
+ *                      file, or NULL.
+ *
+ * \param proj_srid     Spatial reference ID suitable to write to the PROJ_SRID
+ *                      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 
+ */
+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)
+{
+    int ret;
+    char path[GPATH_MAX];
+
+    ret = G_make_location(location_name, wind, proj_info, proj_units);
+
+    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_WKT if WKT is available. */
+    if (proj_wkt != NULL) {
+	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;
+}
+
+/*!
  * \brief Compare projections including units
  *
  *  \param proj_info1   projection info to compare
@@ -432,3 +505,101 @@ int G_compare_projections(const struct Key_Value *proj_info1,
 
     return 1;
 }
+
+/*!
+   \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
+   \param wktstring pointer to WKT string
+
+   \return 0 success
+   \return -1 error writing
+ */
+
+int G_write_projwkt(const char *location_name, const char *wktstring)
+{
+    FILE *fp;
+    char path[GPATH_MAX];
+    int err, n;
+
+    if (!wktstring)
+	return 0;
+
+    if (location_name && *location_name)
+	sprintf(path, "%s/%s/%s/%s", G_gisdbase(), location_name, "PERMANENT", WKT_FILE);
+    else
+	G_file_name(path, "", WKT_FILE, "PERMANENT");
+
+    fp = fopen(path, "w");
+
+    if (!fp)
+	G_fatal_error(_("Unable to open output file <%s>: %s"), path, strerror(errno));
+
+    err = 0;
+    n = strlen(wktstring);
+    if (wktstring[n - 1] != '\n') {
+	if (n != fprintf(fp, "%s\n", wktstring))
+	    err = -1;
+    }
+    else {
+	if (n != fprintf(fp, "%s", wktstring))
+	    err = -1;
+    }
+
+    if (fclose(fp) != 0)
+	G_fatal_error(_("Error closing output file <%s>: %s"), path, strerror(errno));
+
+    return err;
+}
+
+
+/*!
+   \brief Write srid (spatial reference id) to file
+
+    A srid consists of an authority name and code and must be known to
+    PROJ.
+
+   \param location_name name of the location to write the srid
+   \param sridstring pointer to srid string
+
+   \return 0 success
+   \return -1 error writing
+ */
+
+int G_write_projsrid(const char *location_name, const char *sridstring)
+{
+    FILE *fp;
+    char path[GPATH_MAX];
+    int err, n;
+
+    if (!sridstring)
+	return 0;
+
+    if (location_name && *location_name)
+	sprintf(path, "%s/%s/%s/%s", G_gisdbase(), location_name, "PERMANENT", SRID_FILE);
+    else
+	G_file_name(path, "", SRID_FILE, "PERMANENT");
+
+    fp = fopen(path, "w");
+
+    if (!fp)
+	G_fatal_error(_("Unable to open output file <%s>: %s"), path, strerror(errno));
+
+    err = 0;
+    n = strlen(sridstring);
+    if (sridstring[n - 1] != '\n') {
+	if (n != fprintf(fp, "%s\n", sridstring))
+	    err = -1;
+    }
+    else {
+	if (n != fprintf(fp, "%s", sridstring))
+	    err = -1;
+    }
+
+    if (fclose(fp) != 0)
+	G_fatal_error(_("Error closing output file <%s>: %s"), path, strerror(errno));
+
+    return err;
+}

lib/gis/proj3.c

@@ -11,6 +11,9 @@
   \author Original author CERL
  */
 
+/* TODO: the G_database_*() functions should be renamed to G_location_*()
+ * because they apply to a GRASS location, not to a GRASS database */
+
 #include <string.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>