libproj: use new PROJ 5+ API if available

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

include/defs/gprojects.h

@@ -2,14 +2,23 @@
 #define GRASS_GPROJECTSDEFS_H
 
 /* do_proj.c */
+/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
 int pj_do_proj(double *, double *, const struct pj_info *, const struct pj_info *);
 int pj_do_transform(int, double *, double *, double *,
 		    const struct pj_info *, const struct pj_info *);
 
+/* do_proj_ll.c */
+int GPJ_do_proj_ll(double *, double *, const struct pj_info *, int);
+int GPJ_do_transform_ll(int, double *, double *, double *,
+		       const struct pj_info *, int);
+
 /* get_proj.c */
+/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
 int pj_get_kv(struct pj_info *, const struct Key_Value *, const struct Key_Value *);
 int pj_get_string(struct pj_info *, char *);
+#ifndef HAVE_PROJ_H
 int GPJ_get_equivalent_latlong(struct pj_info *, struct pj_info *);
+#endif
 const char *set_proj_lib(const char *);
 int pj_print_proj_params(const struct pj_info *, const struct pj_info *);
 
@@ -42,12 +51,13 @@ int GPJ__get_ellipsoid_params(const struct Key_Value *,
 			      double *, double *, double *);
 void GPJ_free_ellps(struct gpj_ellps *);
 
-
+#ifndef HAVE_PROJ_H
 /* PROJ.4's private datastructures copied from projects.h as removed
    from upstream; pending better solution. see:
    http://trac.osgeo.org/proj/ticket/98 */
 
 int pj_factors(LP, void *, double, struct FACTORS *);
 /* end of copy */
+#endif
 
 #endif

include/gprojects.h

@@ -18,7 +18,14 @@
 #define GRASS_GPROJECTS_H
 
 #include <grass/config.h>
+/* TODO: clean up support for PROJ 5+ */
+#ifdef HAVE_PROJ_H
+#include <proj.h>
+#define RAD_TO_DEG    57.295779513082321
+#define DEG_TO_RAD   .017453292519943296
+#else
 #include <proj_api.h>
+#endif
 #ifdef HAVE_OGR
 #    include <ogr_srs_api.h>
 #endif
@@ -30,12 +37,17 @@
 /* GRASS relative location of datum conversion lookup tables */
 #define GRIDDIR "/etc/proj/nad"
 
+/* TODO: rename pj_ to gpj_ to avoid symbol clash with PROJ lib */
 struct pj_info
 {
+#ifdef HAVE_PROJ_H
+    PJ *pj;
+#else
     projPJ pj;
+#endif
     double meters;
     int zone;
-    char proj[100];
+    char proj[256];
 };
 
 struct gpj_datum
@@ -66,6 +78,7 @@ struct gpj_ellps
     double a, es, rf;
 };
 
+#ifndef HAVE_PROJ_H
 /* PROJ.4's private datastructures copied from projects.h as removed
    from upstream; pending better solution. see:
    http://trac.osgeo.org/proj/ticket/98 */
@@ -87,6 +100,7 @@ struct FACTORS {
 	int code;		/* info as to analytics, see following */
 };
 /* end of copy */
+#endif
 
 #include <grass/defs/gprojects.h>
 

lib/proj/convert.c

@@ -33,6 +33,58 @@
 static void DatumNameMassage(char **);
 #endif
 
+#ifdef HAVE_PROJ_H
+char *gpj_get_def(PJ *P)
+{
+    char *pjdef;
+    PJ_PROJ_INFO pjinfo;
+
+    if (P == NULL)
+	G_fatal_error("Invalid PJ pointer");
+
+    pjinfo = proj_pj_info(P);
+
+    pjdef = G_store(pjinfo.definition);
+
+    return pjdef;
+}
+#endif
+
+
+/* from proj-5.0.0/src/pj_units.c */
+struct gpj_units {
+    char    *id;        /* units keyword */
+    char    *to_meter;  /* multiply by value to get meters */
+    char    *name;      /* comments */
+    double   factor;       /* to_meter factor in actual numbers */
+};
+
+struct gpj_units
+gpj_units[] = {
+    {"km",      "1000.",                "Kilometer",                    1000.0},
+    {"m",       "1.",                   "Meter",                        1.0},
+    {"dm",      "1/10",                 "Decimeter",                    0.1},
+    {"cm",      "1/100",                "Centimeter",                   0.01},
+    {"mm",      "1/1000",               "Millimeter",                   0.001},
+    {"kmi",     "1852.0",               "International Nautical Mile",  1852.0},
+    {"in",      "0.0254",               "International Inch",           0.0254},
+    {"ft",      "0.3048",               "International Foot",           0.3048},
+    {"yd",      "0.9144",               "International Yard",           0.9144},
+    {"mi",      "1609.344",             "International Statute Mile",   1609.344},
+    {"fath",    "1.8288",               "International Fathom",         1.8288},
+    {"ch",      "20.1168",              "International Chain",          20.1168},
+    {"link",    "0.201168",             "International Link",           0.201168},
+    {"us-in",   "1./39.37",             "U.S. Surveyor's Inch",         0.0254},
+    {"us-ft",   "0.304800609601219",    "U.S. Surveyor's Foot",         0.304800609601219},
+    {"us-yd",   "0.914401828803658",    "U.S. Surveyor's Yard",         0.914401828803658},
+    {"us-ch",   "20.11684023368047",    "U.S. Surveyor's Chain",        20.11684023368047},
+    {"us-mi",   "1609.347218694437",    "U.S. Surveyor's Statute Mile", 1609.347218694437},
+    {"ind-yd",  "0.91439523",           "Indian Yard",                  0.91439523},
+    {"ind-ft",  "0.30479841",           "Indian Foot",                  0.30479841},
+    {"ind-ch",  "20.11669506",          "Indian Chain",                 20.11669506},
+    {NULL,      NULL,                   NULL,                           0.0}
+};
+
 static char *grass_to_wkt(const struct Key_Value *proj_info,
                           const struct Key_Value *proj_units,
                           const struct Key_Value *proj_epsg,
@@ -160,11 +212,19 @@ OGRSpatialReferenceH GPJ_grass_to_osr(const struct Key_Value * proj_info,
 	return NULL;
     }
 
+#ifdef HAVE_PROJ_H
+    if ((proj4 = gpj_get_def(pjinfo.pj)) == NULL) {
+	G_warning(_("Unable get PROJ.4-style parameter string"));
+	return NULL;
+    }
+    proj_destroy(pjinfo.pj);
+#else
     if ((proj4 = pj_get_def(pjinfo.pj, 0)) == NULL) {
 	G_warning(_("Unable get PROJ.4-style parameter string"));
 	return NULL;
     }
     pj_free(pjinfo.pj);
+#endif
 
     unit = G_find_key_value("unit", proj_units);
     unfact = G_find_key_value("meters", proj_units);
@@ -172,8 +232,11 @@ OGRSpatialReferenceH GPJ_grass_to_osr(const struct Key_Value * proj_info,
 	G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
     else
 	proj4mod = G_store(proj4);
+#ifdef HAVE_PROJ_H
+    G_free(proj4);
+#else
     pj_dalloc(proj4);
-
+#endif
     if ((errcode = OSRImportFromProj4(hSRS, proj4mod)) != OGRERR_NONE) {
 	G_warning(_("OGR can't parse PROJ.4-style parameter string: "
 		    "%s (OGR Error code was %d)"), proj4mod, errcode);
@@ -337,7 +400,7 @@ OGRSpatialReferenceH GPJ_grass_to_osr2(const struct Key_Value * proj_info,
 
         OSRImportFromEPSG(hSRS, epsgcode);
 
-        /* take +towgs84 from projinfo file if defined) */
+        /* take +towgs84 from projinfo file if defined */
         towgs84 = G_find_key_value("towgs84", proj_info);
         if (towgs84) {
             char **tokens;
@@ -387,6 +450,7 @@ int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
     char *pszProj = NULL;
     const char *pszProjCS = NULL;
     char *datum = NULL;
+    char *proj4_unit = NULL;
     struct gpj_datum dstruct;
     const char *ograttr;
     OGRSpatialReferenceH hSRS;
@@ -554,9 +618,13 @@ int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
 	    continue;
 
 	/* We will handle units separately */
-	if (G_strcasecmp(pszToken, "to_meter") == 0
-	    || G_strcasecmp(pszToken, "units") == 0)
+	if (G_strcasecmp(pszToken, "to_meter") == 0)
+	    continue;
+
+	if (G_strcasecmp(pszToken, "units") == 0) {
+	    proj4_unit = G_store(pszValue);
 	    continue;
+	}
 
 	G_set_key_value(pszToken, pszValue, temp_projinfo);
     }
@@ -789,6 +857,13 @@ int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
 
 	dfToMeters = OSRGetLinearUnits(hSRS, &pszUnitsName);
 
+	/* the unit name can be arbitrary: the following can be the same
+	 * us-ft			(proj.4 keyword)
+	 * U.S. Surveyor's Foot		(proj.4 name)
+	 * US survey foot		(WKT)
+	 * Foot_US			(WKT)
+	 */
+
 	/* Workaround for the most obvious case when unit name is unknown */
 	if ((G_strcasecmp(pszUnitsName, "unknown") == 0) &&
 	    (dfToMeters == 1.))
@@ -798,6 +873,21 @@ int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
 	    G_asprintf(&pszUnitsName, "meter");
 	if ((G_strcasecmp(pszUnitsName, "kilometre") == 0))
 	    G_asprintf(&pszUnitsName, "kilometer");
+	
+	if (dfToMeters != 1. && proj4_unit) {
+	    int i;
+	    
+	    i = 0;
+	    while (gpj_units[i].id != NULL) {
+		if (strcmp(proj4_unit, gpj_units[i].id) == 0) {
+		    if (pszUnitsName)
+			G_free(pszUnitsName);
+		    G_asprintf(&pszUnitsName, gpj_units[i].name);
+		    break;
+		}
+		i++;
+	    }
+	}
 
 	G_set_key_value("unit", pszUnitsName, *projunits);
 

lib/proj/do_proj.c

@@ -43,6 +43,20 @@ do {\
 
 static double METERS_in = 1.0, METERS_out = 1.0;
 
+#ifdef HAVE_PROJ_H
+static char *gpj_get_def(PJ *P)
+{
+    char *pjdef;
+    PJ_PROJ_INFO pj_proj_info = proj_pj_info(P);
+
+    pjdef = G_store(pj_proj_info.definition);
+
+    return pjdef;
+}
+#endif
+
+/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
+
 /** 
  * \brief Re-project a point between two co-ordinate systems
  * 
@@ -57,13 +71,78 @@ static double METERS_in = 1.0, METERS_out = 1.0;
  * \param info_in pointer to pj_info struct for input co-ordinate system
  * \param info_out pointer to pj_info struct for output co-ordinate system
  * 
- * \return Return value from PROJ pj_transform() function
+ * \return Return value from PROJ proj_trans() function
  **/
 
 int pj_do_proj(double *x, double *y,
 	       const struct pj_info *info_in, const struct pj_info *info_out)
 {
     int ok;
+#ifdef HAVE_PROJ_H
+    PJ *P;
+    char *projdef, *projdefin, *projdefout;
+    PJ_COORD c;
+
+    projdefin = gpj_get_def(info_in->pj);
+    projdefout = gpj_get_def(info_out->pj);
+    projdef = NULL;
+    G_asprintf(&projdef, "+proj=pipeline +step %s +inv +step %s", projdefin, projdefout);
+    P = proj_create(PJ_DEFAULT_CTX, projdef);
+    if (P == NULL)
+	G_fatal_error(_("proj_create() failed"));
+    G_free(projdefin);
+    G_free(projdefout);
+    G_free(projdef);
+
+    METERS_in = info_in->meters;
+    METERS_out = info_out->meters;
+
+    if (strncmp(info_in->proj, "ll", 2) == 0) {
+	/* convert to radians */
+	c.lp.lam = (*x) / RAD_TO_DEG;
+	c.lp.phi = (*y) / RAD_TO_DEG;
+	c.lp.z = 0;
+	c.lp.t = 0;
+	c = proj_trans(P, PJ_FWD, c);
+	ok = proj_errno(P);
+
+	if (strncmp(info_out->proj, "ll", 2) == 0) {
+	    /* convert to degrees */
+	    *x = c.lp.lam * RAD_TO_DEG;
+	    *y = c.lp.phi * RAD_TO_DEG;
+	}
+	else {
+	    /* convert to map units */
+	    *x = c.xy.x / METERS_out;
+	    *y = c.xy.y / METERS_out;
+	}
+    }
+    else {
+	/* convert to meters */
+	c.xy.x = *x * METERS_in;
+	c.xy.y = *y * METERS_in;
+	c.xy.z = 0;
+	c.xy.t = 0;
+	c = proj_trans(P, PJ_FWD, c);
+	ok = proj_errno(P);
+
+	if (strncmp(info_out->proj, "ll", 2) == 0) {
+	    /* convert to degrees */
+	    *x = c.lp.lam * RAD_TO_DEG;
+	    *y = c.lp.phi * RAD_TO_DEG;
+	}
+	else {
+	    /* convert to map units */
+	    *x = c.xy.x / METERS_out;
+	    *y = c.xy.y / METERS_out;
+	}
+    }
+    proj_destroy(P);
+
+    if (ok < 0) {
+	G_warning(_("proj_trans() failed: %d"), ok);
+    }
+#else
     double u, v;
     double h = 0.0;
 
@@ -105,6 +184,7 @@ int pj_do_proj(double *x, double *y,
     if (ok < 0) {
 	G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
     }
+#endif
     return ok;
 }
 
@@ -128,7 +208,7 @@ int pj_do_proj(double *x, double *y,
  * \param info_in pointer to pj_info struct for input co-ordinate system
  * \param info_out pointer to pj_info struct for output co-ordinate system
  * 
- * \return Return value from PROJ pj_transform() function
+ * \return Return value from PROJ proj_trans() function
  **/
 
 int pj_do_transform(int count, double *x, double *y, double *h,
@@ -136,7 +216,103 @@ int pj_do_transform(int count, double *x, double *y, double *h,
 {
     int ok;
     int has_h = 1;
+#ifdef HAVE_PROJ_H
+    int i;
+    PJ *P;
+    char *projdef, *projdefin, *projdefout;
+    PJ_COORD c;
 
+    projdefin = gpj_get_def(info_in->pj);
+    projdefout = gpj_get_def(info_out->pj);
+    projdef = NULL;
+    G_asprintf(&projdef, "+proj=pipeline +step %s +inv +step %s", projdefin, projdefout);
+    P = proj_create(PJ_DEFAULT_CTX, projdef);
+    if (P == NULL)
+	G_fatal_error(_("proj_create() failed"));
+    G_free(projdefin);
+    G_free(projdefout);
+    G_free(projdef);
+
+    METERS_in = info_in->meters;
+    METERS_out = info_out->meters;
+
+    if (h == NULL) {
+	h = G_malloc(sizeof *h * count);
+	/* they say memset is only guaranteed for chars ;-( */
+	for (i = 0; i < count; ++i)
+	    h[i] = 0.0;
+	has_h = 0;
+    }
+    ok = 0;
+    if (strncmp(info_in->proj, "ll", 2) == 0) {
+	c.lp.t = 0;
+	if (strncmp(info_out->proj, "ll", 2) == 0) {
+	    for (i = 0; i < count; i++) {
+		/* convert to radians */
+		c.lp.lam = x[i] / RAD_TO_DEG;
+		c.lp.phi = y[i] / RAD_TO_DEG;
+		c.lp.z = h[i];
+		c = proj_trans(P, PJ_FWD, c);
+		if ((ok = proj_errno(P)) < 0)
+		    break;
+		/* convert to degrees */
+		x[i] = c.lp.lam * RAD_TO_DEG;
+		y[i] = c.lp.phi * RAD_TO_DEG;
+	    }
+	}
+	else {
+	    for (i = 0; i < count; i++) {
+		/* convert to radians */
+		c.lp.lam = x[i] / RAD_TO_DEG;
+		c.lp.phi = y[i] / RAD_TO_DEG;
+		c.lp.z = h[i];
+		c = proj_trans(P, PJ_FWD, c);
+		if ((ok = proj_errno(P)) < 0)
+		    break;
+		/* convert to map units */
+		x[i] = c.xy.x / METERS_out;
+		y[i] = c.xy.y / METERS_out;
+	    }
+	}
+    }
+    else {
+	c.xy.t = 0;
+	if (strncmp(info_out->proj, "ll", 2) == 0) {
+	    for (i = 0; i < count; i++) {
+		/* convert to meters */
+		c.xy.x = x[i] * METERS_in;
+		c.xy.y = y[i] * METERS_in;
+		c.xy.z = h[i];
+		c = proj_trans(P, PJ_FWD, c);
+		if ((ok = proj_errno(P)) < 0)
+		    break;
+		/* convert to degrees */
+		x[i] = c.lp.lam * RAD_TO_DEG;
+		y[i] = c.lp.phi * RAD_TO_DEG;
+	    }
+	}
+	else {
+	    for (i = 0; i < count; i++) {
+		/* convert to meters */
+		c.xy.x = x[i] * METERS_in;
+		c.xy.y = y[i] * METERS_in;
+		c.xy.z = h[i];
+		c = proj_trans(P, PJ_FWD, c);
+		if ((ok = proj_errno(P)) < 0)
+		    break;
+		/* convert to map units */
+		x[i] = c.xy.x / METERS_out;
+		y[i] = c.xy.y / METERS_out;
+	    }
+	}
+    }
+    if (!has_h)
+	G_free(h);
+    proj_destroy(P);
+
+    if (ok < 0) {
+	G_warning(_("proj_trans() failed: %d"), ok);
+#else
     METERS_in = info_in->meters;
     METERS_out = info_out->meters;
 
@@ -178,6 +354,7 @@ int pj_do_transform(int count, double *x, double *y, double *h,
 
     if (ok < 0) {
 	G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
+#endif
     }
     return ok;
 }

lib/proj/do_proj_ll.c

@@ -0,0 +1,221 @@
+
+/**
+   \file do_proj.c
+
+   \brief GProj library - Functions for re-projecting point data
+
+   \author Original Author unknown, probably Soil Conservation Service
+   Eric Miller, Paul Kelly
+
+   (C) 2003-2008 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.
+**/
+
+#include <stdio.h>
+#include <string.h>
+#include <ctype.h>
+
+#include <grass/gis.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+
+static double METERS_in = 1.0;
+
+#ifdef HAVE_PROJ_H
+static char *gpj_get_def(PJ *P)
+{
+    char *pjdef;
+    PJ_PROJ_INFO pj_proj_info = proj_pj_info(P);
+
+    pjdef = G_store(pj_proj_info.definition);
+
+    return pjdef;
+}
+
+/** 
+ * \brief Re-project a point between a co-ordinate system and its 
+ *        latitude / longitude equivalent, using the same datum, 
+ *        ellipsoid and datum transformation
+ * 
+ * This function takes a pointer to a pj_info structure as argument, 
+ * and projects a point between the co-ordinate system and its 
+ * latitude / longitude equivalent. With forward projection, the point 
+ * is projected from the co-ordinate system to its latitude / longitude 
+ * equivalent, while backward projection does the reverse.
+ * The easting and northing of the point are contained in two pointers 
+ * passed to the function; these will be overwritten by the co-ordinates
+ * of the re-projected point.
+ * 
+ * \param x Pointer to a double containing easting or longitude
+ * \param y Pointer to a double containing northing or latitude
+ * \param info_in pointer to pj_info struct for input co-ordinate system
+ * \param direction direction of re-projection (PJ_FWD or PJ_INV)
+ * 
+ * \return Return value from PROJ proj_trans() function
+ **/
+
+int GPJ_do_proj_ll(double *x, double *y,
+	          const struct pj_info *info_in, int direction)
+{
+    int ok;
+    PJ *P;
+    char *projdef;
+    PJ_COORD c;
+
+    ok = 0;
+    if (strncmp(info_in->proj, "ll", 2) == 0) {
+	/* nothing to do */
+	return ok;
+    }
+
+    projdef = gpj_get_def(info_in->pj);
+    P = proj_create(PJ_DEFAULT_CTX, projdef);
+    if (P == NULL)
+	G_fatal_error(_("proj_create() failed"));
+    G_free(projdef);
+
+    METERS_in = info_in->meters;
+
+    if (direction == PJ_FWD) {
+	/* from ll to projected */
+
+	/* convert to radians */
+	c.lp.lam = (*x) / RAD_TO_DEG;
+	c.lp.phi = (*y) / RAD_TO_DEG;
+
+	c = proj_trans(P, PJ_FWD, c);
+	ok = proj_errno(P);
+
+	/* convert to map units */
+	*x = c.xy.x / METERS_in;
+	*y = c.xy.y / METERS_in;
+    }
+    else {
+	/* from projected to ll */
+
+	/* convert to meters */
+	c.xy.x = *x * METERS_in;
+	c.xy.y = *y * METERS_in;
+
+	c = proj_trans(P, PJ_INV, c);
+	ok = proj_errno(P);
+
+	/* convert to degrees */
+	*x = c.lp.lam * RAD_TO_DEG;
+	*y = c.lp.phi * RAD_TO_DEG;
+    }
+    proj_destroy(P);
+
+    if (ok < 0) {
+	G_warning(_("proj_trans() failed: %d"), ok);
+    }
+    return ok;
+}
+
+/** 
+ * \brief Re-project an array of points between a co-ordinate system 
+ *        and its latitude / longitude equivalent, using the same 
+ *        datum, ellipsoid and datum transformation
+ * 
+ * This function takes a pointer to a pj_info structure as argument, 
+ * and projects an array of points between the co-ordinate system and its 
+ * latitude / longitude equivalent. With forward projection, the points 
+ * are projected from the co-ordinate system to its latitude / longitude 
+ * equivalent, while backward projection does the reverse.
+ * The easting and northing of the points are contained in two pointers passed 
+ * to the function; these will be overwritten by the co-ordinates of the 
+ * re-projected point.
+ * 
+ * \param count Number of points in the arrays to be transformed
+ * \param x Pointer to an array of type double containing easting or longitude
+ * \param y Pointer to an array of type double containing northing or latitude
+ * \param h Pointer to an array of type double containing ellipsoidal height. 
+ *          May be null in which case a two-dimensional re-projection will be 
+ *          done
+ * \param info_in pointer to pj_info struct for input co-ordinate system
+ * \param direction direction of re-projection (PJ_FWD or PJ_INV)
+ * 
+ * \return Return value from PROJ proj_trans() function
+ **/
+
+int GPJ_do_transform_ll(int count, double *x, double *y, double *h,
+		       const struct pj_info *info_in, int direction)
+{
+    int ok;
+    int has_h = 1;
+    int i;
+    PJ *P;
+    char *projdef;
+    PJ_COORD c;
+
+    ok = 0;
+    if (strncmp(info_in->proj, "ll", 2) == 0) {
+	/* nothing to do */
+	return ok;
+    }
+
+    projdef = gpj_get_def(info_in->pj);
+    P = proj_create(PJ_DEFAULT_CTX, projdef);
+    if (P == NULL)
+	G_fatal_error(_("proj_create() failed"));
+    G_free(projdef);
+
+    METERS_in = info_in->meters;
+
+    if (h == NULL) {
+	h = G_malloc(sizeof *h * count);
+	/* they say memset is only guaranteed for chars ;-( */
+	for (i = 0; i < count; ++i)
+	    h[i] = 0.0;
+	has_h = 0;
+    }
+
+    if (direction == PJ_FWD) {
+	/* from ll to projected */
+
+	for (i = 0; i < count; i++) {
+	    /* convert to radians */
+	    c.lp.lam = x[i] / RAD_TO_DEG;
+	    c.lp.phi = y[i] / RAD_TO_DEG;
+
+	    c = proj_trans(P, PJ_FWD, c);
+	    if ((ok = proj_errno(P)) < 0)
+		break;
+
+	    /* convert to map units */
+	    x[i] = c.xy.x / METERS_in;
+	    y[i] = c.xy.y / METERS_in;
+	}
+    }
+    else {
+	/* from projected to ll */
+
+	for (i = 0; i < count; i++) {
+	    /* convert to meters */
+	    c.xy.x = x[i] * METERS_in;
+	    c.xy.y = y[i] * METERS_in;
+
+	    c = proj_trans(P, PJ_INV, c);
+	    if ((ok = proj_errno(P)) < 0)
+		break;
+
+	    /* convert to degrees */
+	    x[i] = c.lp.lam * RAD_TO_DEG;
+	    y[i] = c.lp.phi * RAD_TO_DEG;
+	}
+    }
+    if (!has_h)
+	G_free(h);
+    proj_destroy(P);
+
+    if (ok < 0) {
+	G_warning(_("proj_trans() failed: %d"), ok);
+    }
+    return ok;
+}
+
+#endif /* HAVE_PROJ_H */
+

lib/proj/get_proj.c

@@ -31,7 +31,26 @@
 static void alloc_options(char *);
 
 static char *opt_in[MAX_PARGS];
-static int nopt1;
+static int nopt;
+
+#ifdef HAVE_PROJ_H
+static char *gpj_get_def(PJ *P)
+{
+    char *pjdef;
+    PJ_PROJ_INFO pjinfo;
+
+    if (P == NULL)
+	G_fatal_error("Invalid PJ pointer");
+
+    pjinfo = proj_pj_info(P);
+
+    pjdef = G_store(pjinfo.definition);
+
+    return pjdef;
+}
+#endif
+
+/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
 
 /**
  * \brief Create a pj_info struct Co-ordinate System definition from a set of
@@ -65,7 +84,12 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
     int returnval = 1;
     char buffa[300], factbuff[50];
     char proj_in[250], *datum, *params;
+#ifdef HAVE_PROJ_H
+    PJ *pj;
+    PJ_CONTEXT *pjc;
+#else
     projPJ *pj;
+#endif
 
     proj_in[0] = '\0';
     info->zone = 0;
@@ -89,7 +113,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
     if (strlen(info->proj) <= 0)
 	sprintf(info->proj, "ll");
 
-    nopt1 = 0;
+    nopt = 0;
     for (i = 0; i < in_proj_keys->nitems; i++) {
 	/* the name parameter is just for grasses use */
 	if (strcmp(in_proj_keys->key[i], "name") == 0) {
@@ -220,24 +244,42 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
     }
     G_free(datum);
 
+#ifdef HAVE_PROJ_H
+    pjc = proj_context_create();
+    if (!(pj = proj_create_argv(pjc, nopt, opt_in))) {
+#else
     /* Set finder function for locating datum conversion tables PK */
     pj_set_finder(FINDERFUNC);
 
-    if (!(pj = pj_init(nopt1, opt_in))) {
+    if (!(pj = pj_init(nopt, opt_in))) {
+#endif
 	strcpy(buffa,
-	       _("Unable to initialise PROJ.4 with the following parameter list:"));
-	for (i = 0; i < nopt1; i++) {
+	       _("Unable to initialise PROJ with the following parameter list:"));
+	for (i = 0; i < nopt; i++) {
 	    char err[50];
 
 	    sprintf(err, " +%s", opt_in[i]);
 	    strcat(buffa, err);
 	}
 	G_warning("%s", buffa);
-	G_warning(_("The error message: %s"), pj_strerrno(pj_errno));
+#ifndef HAVE_PROJ_H
+	G_warning(_("The PROJ error message: %s"), pj_strerrno(pj_errno));
+#endif
 	return -1;
     }
+
+#ifdef HAVE_PROJ_H
+    int perr = proj_errno(pj);
+
+    if (perr)
+	G_fatal_error("PROJ 5 error %d", perr);
+#endif
+    
     info->pj = pj;
 
+    for (i = 0; i < nopt; i++)
+	G_free(opt_in[i]);
+
     return returnval;
 }
 
@@ -246,23 +288,28 @@ static void alloc_options(char *buffa)
     int nsize;
 
     nsize = strlen(buffa);
-    opt_in[nopt1++] = (char *)G_malloc(nsize + 1);
-    sprintf(opt_in[nopt1 - 1], "%s", buffa);
+    opt_in[nopt++] = (char *)G_malloc(nsize + 1);
+    sprintf(opt_in[nopt - 1], "%s", buffa);
     return;
 }
 
 int pj_get_string(struct pj_info *info, char *str)
 {
-    char *opt_in[MAX_PARGS];
     char *s;
-    int nopt = 0;
-    int nsize;
+    int i, nsize;
     char zonebuff[50], buffa[300];
+#ifdef HAVE_PROJ_H
+    PJ *pj;
+    PJ_CONTEXT *pjc;
+#else
     projPJ *pj;
+#endif
 
     info->zone = 0;
     info->proj[0] = '\0';
     info->meters = 1.0;
+    
+    nopt = 0;
 
     if ((str == NULL) || (str[0] == '\0')) {
 	/* Null Pointer or empty string is supplied for parameters, 
@@ -270,9 +317,7 @@ int pj_get_string(struct pj_info *info, char *str)
 	 * parameter and call pj_init PK */
 	sprintf(info->proj, "ll");
 	sprintf(buffa, "proj=latlong ellps=WGS84");
-	nsize = strlen(buffa);
-	opt_in[nopt] = (char *)G_malloc(nsize + 1);
-	sprintf(opt_in[nopt++], "%s", buffa);
+	alloc_options(buffa);
     }
     else {
 	/* Parameters have been provided; parse through them but don't
@@ -309,15 +354,20 @@ int pj_get_string(struct pj_info *info, char *str)
 		    else {
 			sprintf(buffa, "%s", s);
 		    }
-		    nsize = strlen(buffa);
-		    opt_in[nopt] = (char *)G_malloc(nsize + 1);
-		    sprintf(opt_in[nopt++], "%s", buffa);
+		    alloc_options(buffa);
 		}
 	    }
 	    s = 0;
 	}
     }
 
+#ifdef HAVE_PROJ_H
+    pjc = proj_context_create();
+    if (!(pj = proj_create_argv(pjc, nopt, opt_in))) {
+	G_warning(_("Unable to initialize pj"));
+	return -1;
+    }
+#else
     /* Set finder function for locating datum conversion tables PK */
     pj_set_finder(FINDERFUNC);
 
@@ -326,8 +376,12 @@ int pj_get_string(struct pj_info *info, char *str)
 		  pj_strerrno(pj_errno));
 	return -1;
     }
+#endif
     info->pj = pj;
 
+    for (i = 0; i < nopt; i++)
+	G_free(opt_in[i]);
+
     return 1;
 }
 
@@ -349,6 +403,10 @@ int pj_get_string(struct pj_info *info, char *str)
 
 int GPJ_get_equivalent_latlong(struct pj_info *pjnew, struct pj_info *pjold)
 {
+#ifdef HAVE_PROJ_H
+    G_fatal_error(_("GPJ_get_equivalent_latlong(): with the new PROJ 5+ API "
+                    "use the old pj directly with PJ_FWD/PJ_INV transformation"));
+#else
     pjnew->meters = 1.;
     pjnew->zone = 0;
     sprintf(pjnew->proj, "ll");
@@ -356,6 +414,7 @@ int GPJ_get_equivalent_latlong(struct pj_info *pjnew, struct pj_info *pjold)
 	return -1;
     else
 	return 1;
+#endif
 }
 
 /* set_proj_lib()
@@ -396,11 +455,19 @@ int pj_print_proj_params(const struct pj_info *iproj, const struct pj_info *opro
     char *str;
 
     if (iproj) {
+#ifdef HAVE_PROJ_H
+	str = gpj_get_def(iproj->pj);
+#else
 	str = pj_get_def(iproj->pj, 1);
+#endif
 	if (str != NULL) {
 	    fprintf(stderr, "%s: %s\n", _("Input Projection Parameters"),
 		    str);
+#ifdef HAVE_PROJ_H
+	    G_free(str);
+#else
 	    pj_dalloc(str);
+#endif
 	    fprintf(stderr, "%s: %.16g\n", _("Input Unit Factor"),
 		    iproj->meters);
 	}
@@ -409,11 +476,19 @@ int pj_print_proj_params(const struct pj_info *iproj, const struct pj_info *opro
     }
 
     if (oproj) {
+#ifdef HAVE_PROJ_H
+	str = gpj_get_def(oproj->pj);
+#else
 	str = pj_get_def(oproj->pj, 1);
+#endif
 	if (str != NULL) {
 	    fprintf(stderr, "%s: %s\n", _("Output Projection Parameters"),
 		    str);
+#ifdef HAVE_PROJ_H
+	    G_free(str);
+#else
 	    pj_dalloc(str);
+#endif
 	    fprintf(stderr, "%s: %.16g\n", _("Output Unit Factor"),
 		    oproj->meters);
 	}