libproj: +new GRASS API for coordinate transformation

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

include/defs/gprojects.h

@@ -2,16 +2,19 @@
 #define GRASS_GPROJECTSDEFS_H
 
 /* do_proj.c */
-/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
+int GPJ_init_transform(const struct pj_info *, const struct pj_info *,
+		       struct pj_info *);
+int GPJ_transform(const struct pj_info *, const struct pj_info *,
+                  const struct pj_info *, int, double *, double *, double *);
+int GPJ_transform_array(const struct pj_info *, const struct pj_info *,
+                        const struct pj_info *, int,
+		        double *, double *, double *, int);
+
+/* old API, to be removed */
 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 *);

lib/proj/convert.c

@@ -5,12 +5,12 @@
    \brief GProj Library - Functions for manipulating co-ordinate
    system representations
 
-   (C) 2003-2017 by the GRASS Development Team
+   (C) 2003-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.
 
-   \author Paul Kelly, Frank Warmerdam 
+   \author Paul Kelly, Frank Warmerdam, Markus Metz 
 */
 
 #include <grass/config.h>

lib/proj/do_proj.c

@@ -5,9 +5,9 @@
    \brief GProj library - Functions for re-projecting point data
 
    \author Original Author unknown, probably Soil Conservation Service
-   Eric Miller, Paul Kelly
+   Eric Miller, Paul Kelly, Markus Metz
 
-   (C) 2003-2008 by the GRASS Development Team
+   (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
@@ -43,48 +43,464 @@ do {\
 
 static double METERS_in = 1.0, METERS_out = 1.0;
 
-/* TODO: add to gprojects.h */
-/* Create a PROJ transformation object 
- * to transform coordinates from an input SRS to an output SRS */
-/* PROJ 5+:
- *   info_in, info_trans must not be null
- *   if info_out is null, assume info_out to be ll equivalent of info_in
- *   create info_trans from info_in to info_out
- *   NOTE: this is against the logic of PROJ 5 which by default 
- *         converts from ll to a given SRS
+/**
+ * \brief Create a PROJ transformation object to transform coordinates 
+ *        from an input SRS to an output SRS 
+ * 
+ * After the transformation has been initialized with this function, 
+ * coordinates can be transformed from input SRS to output SRS with 
+ * GPJ_transform() and direction = PJ_FWD, and back from output SRS to
+ * input SRS with direction = OJ_INV.
+ * If coordinates should be transformed between the input SRS and its 
+ * latlong equivalent, an uninitialized info_out with 
+ * info_out->pj = NULL can be passed to the function. In this case,
+ * coordinates will be transformed between the input SRS and its 
+ * latlong equivalent, and for PROJ 5+, the transformation object is 
+ * created accordingly, while for PROJ 4, the output SRS is created as 
+ * latlong equivalent of the input SRS
+ * 
+* PROJ 5+:
+ *   info_in->pj must not be null
+ *   if info_out->pj is null, assume info_out to be the ll equivalent 
+ *   of info_in
+ *   create info_trans as conversion from info_in to its ll equivalent
+ *   NOTE: this is the inverse of the logic of PROJ 5 which by default 
+ *         converts from ll to a given SRS, not from a given SRS to ll
  *         thus PROJ 5+ itself uses an inverse transformation in the
  *         first step of the pipeline for proj_create_crs_to_crs()
  * PROJ 4:
- *   info_in, info_out must not be null
- *   do nothing
- */
-int GPJ_prepare_pjinfo(const struct pj_info *info_in,
+ *   info_in->pj must not be null
+ *   if info_out->pj is null, create info_out as ll equivalent
+ *   else do nothing, info_trans is not used
+ * 
+ * \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
+ * \param info_trans pointer to pj_info struct for a transformation object (PROJ 5+)
+ * 
+ * \return 1 on success, -1 on failure
+ **/
+int GPJ_init_transform(const struct pj_info *info_in,
                        const struct pj_info *info_out,
 		       struct pj_info *info_trans)
 {
-    if (info_in == NULL)
+    if (info_in->pj == NULL)
 	G_fatal_error(_("Input coordinate system is NULL"));
 
 #ifdef HAVE_PROJ_H
-    if (info_trans == NULL)
-	G_fatal_error(_("Transformation object for PROJ is NULL"));
+    if (info_out->pj != NULL && info_out->def != NULL)
+	G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s +step %s",
+		   info_in->def, info_out->def);
+    else
+	/* assume info_out to be ll equivalent of info_in */
+	G_asprintf(&(info_trans->def), "+proj=pipeline +step +inv %s",
+		   info_in->def);
 
-    G_asprintf(&(info_trans->def), "+proj=pipeline +step %s +inv +step %s",
-                         info_in->def, info_out->def);
     info_trans->pj = proj_create(PJ_DEFAULT_CTX, info_trans->def);
-    if (info_trans->pj == NULL)
-	G_fatal_error(_("proj_create() failed for '%s'"), info_trans->def);
+    if (info_trans->pj == NULL) {
+	G_warning(_("proj_create() failed for '%s'"), info_trans->def);
+	return -1;
+    }
+
+    info_trans->meters = 1.;
+    info_trans->zone = 0;
+    sprintf(info_trans->proj, "pipeline");
+#else
+    if (info_out->pj == NULL) {
+	if (GPJ_get_equivalent_latlong(info_out, info_in) < 0) {
+	    G_warning(_("Unable to create latlong equivalent for '%s'"),
+	              info_in->def);
+	    return -1;
+	}
+    }
+#endif
 
     return 1;
+}
+
+/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
+
+/** 
+ * \brief Re-project a point between two co-ordinate systems using a 
+ *        transformation object prepared with GPJ_prepare_pj()
+ * 
+ * This function takes pointers to three pj_info structures as arguments, 
+ * and projects a point between the input and output co-ordinate system.
+ * The pj_info structure info_trans must have been initialized with 
+ * GPJ_init_transform().
+ * The direction determines if a point is projected from input CRS to 
+ * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
+ * The easting, northing, and height of the point are contained in the 
+ * pointers passed to the function; these will be overwritten by the 
+ * coordinates of the transformed point.
+ * 
+ * \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
+ * \param info_trans pointer to pj_info struct for a transformation object (PROJ 5+)
+ * \param dir direction of the transformation (PJ_FWD or PJ_INV)
+ * \param x Pointer to a double containing easting or longitude
+ * \param y Pointer to a double containing northing or latitude
+ * \param z Pointer to a double containing height, or NULL
+ * 
+ * \return Return value from PROJ proj_trans() function
+ **/
+
+int GPJ_transform(const struct pj_info *info_in,
+                  const struct pj_info *info_out,
+                  const struct pj_info *info_trans, int dir,
+		  double *x, double *y, double *z)
+{
+    int ok = 0;
+
+#ifdef HAVE_PROJ_H
+    /* PROJ 5+ variant */
+    int in_is_ll, out_is_ll;
+    PJ_COORD c;
+
+    if (info_in->pj == NULL)
+	G_fatal_error(_("No input projection"));
+
+    if (info_trans->pj == NULL)
+	G_fatal_error(_("No transformation object"));
+
+    if (dir == PJ_FWD) {
+	/* info_in -> info_out */
+	METERS_in = info_in->meters;
+	in_is_ll = !strncmp(info_in->proj, "ll", 2);
+	if (info_out->pj) {
+	    METERS_out = info_out->meters;
+	    out_is_ll = !strncmp(info_out->proj, "ll", 2);
+	}
+	else {
+	    METERS_out = 1.0;
+	    out_is_ll = 1;
+	}
+    }
+    else {
+	/* info_out -> info_in */
+	METERS_out = info_in->meters;
+	out_is_ll = !strncmp(info_in->proj, "ll", 2);
+	if (info_out->pj) {
+	    METERS_in = info_out->meters;
+	    in_is_ll = !strncmp(info_out->proj, "ll", 2);
+	}
+	else {
+	    METERS_in = 1.0;
+	    in_is_ll = 1;
+	}
+    }
+
+    /* prepare */
+    if (in_is_ll) {
+	/* convert to radians */
+	c.lpzt.lam = (*x) / RAD_TO_DEG;
+	c.lpzt.phi = (*y) / RAD_TO_DEG;
+	c.lpzt.z = 0;
+	if (z)
+	    c.lpzt.z = *z;
+	c.lpzt.t = 0;
+    }
+    else {
+	/* convert to meters */
+	c.xyzt.x = *x * METERS_in;
+	c.xyzt.y = *y * METERS_in;
+	c.xyzt.z = 0;
+	if (z)
+	    c.xyzt.z = *z;
+	c.xyzt.t = 0;
+    }
+
+    /* transform */
+    c = proj_trans(info_trans->pj, dir, c);
+    ok = proj_errno(info_trans->pj);
+
+    if (ok < 0) {
+	G_warning(_("proj_trans() failed: %d"), ok);
+	return ok;
+    }
+
+    /* output */
+    if (out_is_ll) {
+	/* convert to degrees */
+	*x = c.lpzt.lam * RAD_TO_DEG;
+	*y = c.lpzt.phi * RAD_TO_DEG;
+	if (z)
+	    *z = c.lpzt.z;
+    }
+    else {
+	/* convert to map units */
+	*x = c.xyzt.x / METERS_out;
+	*y = c.xyzt.y / METERS_out;
+	if (z)
+	    *z = c.xyzt.z;
+    }
 #else
+    /* PROJ 4 variant */
+    double u, v;
+    double h = 0.0;
+    const struct pj_info *p_in, *p_out;
+
     if (info_out == NULL)
-	G_fatal_error(_("Output coordinate system is NULL"));
+	G_fatal_error(_("No output projection"));
 
-    return 1;
+    if (dir == FWD) {
+	p_in = info_in;
+	p_out = info_out;
+    }
+    else {
+	p_in = info_out;
+	p_out = info_in;
+    }
+
+    METERS_in = p_in->meters;
+    METERS_out = p_out->meters;
+    
+    if (z)
+	h = *z;
+
+    if (strncmp(p_in->proj, "ll", 2) == 0) {
+	u = (*x) / RAD_TO_DEG;
+	v = (*y) / RAD_TO_DEG;
+    }
+    else {
+	u = *x * METERS_in;
+	v = *y * METERS_in;
+    }
+
+    ok = pj_transform(p_in->pj, p_out->pj, 1, 0, &u, &v, &h);
+
+    if (ok < 0) {
+	G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
+	return ok;
+    }
+
+    if (strncmp(p_out->proj, "ll", 2) == 0) {
+	*x = u * RAD_TO_DEG;
+	*y = v * RAD_TO_DEG;
+    }
+    else {
+	*x = u / METERS_out;
+	*y = v / METERS_out;
+    }
+    if (z)
+	*z = h;
 #endif
+
+    return ok;
 }
 
-/* TODO: rename pj_ to GPJ_ to avoid symbol clash with PROJ lib */
+/** 
+ * \brief Re-project an array of points between two co-ordinate systems 
+ *        using a transformation object prepared with GPJ_prepare_pj()
+ * 
+ * This function takes pointers to three pj_info structures as arguments, 
+ * and projects an array of pointd between the input and output 
+ * co-ordinate system. The pj_info structure info_trans must have been 
+ * initialized with GPJ_init_transform().
+ * The direction determines if a point is projected from input CRS to 
+ * output CRS (PJ_FWD) or from output CRS to input CRS (PJ_INV).
+ * The easting, northing, and height of the point are contained in the 
+ * pointers passed to the function; these will be overwritten by the 
+ * coordinates of the transformed point.
+ * 
+ * \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
+ * \param info_trans pointer to pj_info struct for a transformation object (PROJ 5+)
+ * \param dir direction of the transformation (PJ_FWD or PJ_INV)
+ * \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 z pointer to an array of type double containing height, or NULL
+ * \param n number of points in the arrays to be transformed
+ * 
+ * \return Return value from PROJ proj_trans() function
+ **/
+
+int GPJ_transform_array(const struct pj_info *info_in,
+                        const struct pj_info *info_out,
+                        const struct pj_info *info_trans, int dir,
+		        double *x, double *y, double *z, int n)
+{
+    int ok;
+    int i;
+    int has_z = 1;
+
+#ifdef HAVE_PROJ_H
+    /* PROJ 5+ variant */
+    int in_is_ll, out_is_ll;
+    PJ_COORD c;
+
+    if (info_trans->pj == NULL)
+	G_fatal_error(_("No transformation object"));
+
+    if (dir == PJ_FWD) {
+	/* info_in -> info_out */
+	METERS_in = info_in->meters;
+	in_is_ll = !strncmp(info_in->proj, "ll", 2);
+	if (info_out->pj) {
+	    METERS_out = info_out->meters;
+	    out_is_ll = !strncmp(info_out->proj, "ll", 2);
+	}
+	else {
+	    METERS_out = 1.0;
+	    out_is_ll = 1;
+	}
+    }
+    else {
+	/* info_out -> info_in */
+	METERS_out = info_in->meters;
+	out_is_ll = !strncmp(info_in->proj, "ll", 2);
+	if (info_out->pj) {
+	    METERS_in = info_out->meters;
+	    in_is_ll = !strncmp(info_out->proj, "ll", 2);
+	}
+	else {
+	    METERS_in = 1.0;
+	    in_is_ll = 1;
+	}
+    }
+
+    if (z == NULL) {
+	z = G_malloc(sizeof(double) * n);
+	/* they say memset is only guaranteed for chars ;-( */
+	for (i = 0; i < n; i++)
+	    z[i] = 0.0;
+	has_z = 0;
+    }
+    ok = 0;
+    if (in_is_ll) {
+	c.lpzt.t = 0;
+	if (out_is_ll) {
+	    /* what is more costly ?
+	     * calling proj_trans for each point
+	     * or having three loops over all points ?
+	     * proj_trans_array() itself calls proj_trans() in a loop
+	     * -> one loop over all points is better than 
+	     *    three loops over all points 
+	     */
+	    for (i = 0; i < n; i++) {
+		/* convert to radians */
+		c.lpzt.lam = x[i] / RAD_TO_DEG;
+		c.lpzt.phi = y[i] / RAD_TO_DEG;
+		c.lpzt.z = z[i];
+		c = proj_trans(info_trans->pj, dir, c);
+		if ((ok = proj_errno(info_trans->pj)) < 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 < n; i++) {
+		/* convert to radians */
+		c.lpzt.lam = x[i] / RAD_TO_DEG;
+		c.lpzt.phi = y[i] / RAD_TO_DEG;
+		c.lpzt.z = z[i];
+		c = proj_trans(info_trans->pj, dir, c);
+		if ((ok = proj_errno(info_trans->pj)) < 0)
+		    break;
+		/* convert to map units */
+		x[i] = c.xy.x / METERS_out;
+		y[i] = c.xy.y / METERS_out;
+	    }
+	}
+    }
+    else {
+	c.xyzt.t = 0;
+	if (out_is_ll) {
+	    for (i = 0; i < n; i++) {
+		/* convert to meters */
+		c.xyzt.x = x[i] * METERS_in;
+		c.xyzt.y = y[i] * METERS_in;
+		c.xyzt.z = z[i];
+		c = proj_trans(info_trans->pj, dir, c);
+		if ((ok = proj_errno(info_trans->pj)) < 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 < n; i++) {
+		/* convert to meters */
+		c.xyzt.x = x[i] * METERS_in;
+		c.xyzt.y = y[i] * METERS_in;
+		c.xyzt.z = z[i];
+		c = proj_trans(info_trans->pj, dir, c);
+		if ((ok = proj_errno(info_trans->pj)) < 0)
+		    break;
+		/* convert to map units */
+		x[i] = c.xy.x / METERS_out;
+		y[i] = c.xy.y / METERS_out;
+	    }
+	}
+    }
+    if (!has_z)
+	G_free(z);
+
+    if (ok < 0) {
+	G_warning(_("proj_trans() failed: %d"), ok);
+    }
+#else
+    /* PROJ 4 variant */
+    const struct pj_info *p_in, *p_out;
+
+    if (dir == FWD) {
+	p_in = info_in;
+	p_out = info_out;
+    }
+    else {
+	p_in = info_out;
+	p_out = info_in;
+    }
+
+    METERS_in = p_in->meters;
+    METERS_out = p_out->meters;
+
+    if (z == NULL) {
+	z = G_malloc(sizeof(double) * n);
+	/* they say memset is only guaranteed for chars ;-( */
+	for (i = 0; i < n; ++i)
+	    z[i] = 0.0;
+	has_z = 0;
+    }
+    if (strncmp(p_in->proj, "ll", 2) == 0) {
+	if (strncmp(p_out->proj, "ll", 2) == 0) {
+	    DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
+	    ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
+	    MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
+	}
+	else {
+	    DIVIDE_LOOP(x, y, n, RAD_TO_DEG);
+	    ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
+	    DIVIDE_LOOP(x, y, n, METERS_out);
+	}
+    }
+    else {
+	if (strncmp(p_out->proj, "ll", 2) == 0) {
+	    MULTIPLY_LOOP(x, y, n, METERS_in);
+	    ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
+	    MULTIPLY_LOOP(x, y, n, RAD_TO_DEG);
+	}
+	else {
+	    MULTIPLY_LOOP(x, y, n, METERS_in);
+	    ok = pj_transform(info_in->pj, info_out->pj, n, 1, x, y, z);
+	    DIVIDE_LOOP(x, y, n, METERS_out);
+	}
+    }
+    if (!has_z)
+	G_free(z);
+
+    if (ok < 0)
+	G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
+#endif
+
+    return ok;
+}
+
+/*
+ * old API, to be deleted
+ */
 
 /** 
  * \brief Re-project a point between two co-ordinate systems
@@ -108,10 +524,12 @@ int pj_do_proj(double *x, double *y,
 {
     int ok;
 #ifdef HAVE_PROJ_H
-    struct pj_info info_new;
+    struct pj_info info_trans;
     PJ_COORD c;
 
-    GPJ_prepare_pjinfo(info_in, info_out, &info_new);
+    if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
+	return -1;
+    }
 
     METERS_in = info_in->meters;
     METERS_out = info_out->meters;
@@ -122,8 +540,8 @@ int pj_do_proj(double *x, double *y,
 	c.lpzt.phi = (*y) / RAD_TO_DEG;
 	c.lpzt.z = 0;
 	c.lpzt.t = 0;
-	c = proj_trans(info_new.pj, PJ_FWD, c);
-	ok = proj_errno(info_new.pj);
+	c = proj_trans(info_trans.pj, PJ_FWD, c);
+	ok = proj_errno(info_trans.pj);
 
 	if (strncmp(info_out->proj, "ll", 2) == 0) {
 	    /* convert to degrees */
@@ -142,8 +560,8 @@ int pj_do_proj(double *x, double *y,
 	c.xyzt.y = *y * METERS_in;
 	c.xyzt.z = 0;
 	c.xyzt.t = 0;
-	c = proj_trans(info_new.pj, PJ_FWD, c);
-	ok = proj_errno(info_new.pj);
+	c = proj_trans(info_trans.pj, PJ_FWD, c);
+	ok = proj_errno(info_trans.pj);
 
 	if (strncmp(info_out->proj, "ll", 2) == 0) {
 	    /* convert to degrees */
@@ -156,7 +574,7 @@ int pj_do_proj(double *x, double *y,
 	    *y = c.xy.y / METERS_out;
 	}
     }
-    proj_destroy(info_new.pj);
+    proj_destroy(info_trans.pj);
 
     if (ok < 0) {
 	G_warning(_("proj_trans() failed: %d"), ok);
@@ -237,10 +655,12 @@ int pj_do_transform(int count, double *x, double *y, double *h,
     int has_h = 1;
 #ifdef HAVE_PROJ_H
     int i;
-    struct pj_info info_new;
+    struct pj_info info_trans;
     PJ_COORD c;
 
-    GPJ_prepare_pjinfo(info_in, info_out, &info_new);
+    if (GPJ_init_transform(info_in, info_out, &info_trans) < 0) {
+	return -1;
+    }
 
     METERS_in = info_in->meters;
     METERS_out = info_out->meters;
@@ -261,8 +681,8 @@ int pj_do_transform(int count, double *x, double *y, double *h,
 		c.lpzt.lam = x[i] / RAD_TO_DEG;
 		c.lpzt.phi = y[i] / RAD_TO_DEG;
 		c.lpzt.z = h[i];
-		c = proj_trans(info_new.pj, PJ_FWD, c);
-		if ((ok = proj_errno(info_new.pj)) < 0)
+		c = proj_trans(info_trans.pj, PJ_FWD, c);
+		if ((ok = proj_errno(info_trans.pj)) < 0)
 		    break;
 		/* convert to degrees */
 		x[i] = c.lp.lam * RAD_TO_DEG;
@@ -275,8 +695,8 @@ int pj_do_transform(int count, double *x, double *y, double *h,
 		c.lpzt.lam = x[i] / RAD_TO_DEG;
 		c.lpzt.phi = y[i] / RAD_TO_DEG;
 		c.lpzt.z = h[i];
-		c = proj_trans(info_new.pj, PJ_FWD, c);
-		if ((ok = proj_errno(info_new.pj)) < 0)
+		c = proj_trans(info_trans.pj, PJ_FWD, c);
+		if ((ok = proj_errno(info_trans.pj)) < 0)
 		    break;
 		/* convert to map units */
 		x[i] = c.xy.x / METERS_out;
@@ -292,8 +712,8 @@ int pj_do_transform(int count, double *x, double *y, double *h,
 		c.xyzt.x = x[i] * METERS_in;
 		c.xyzt.y = y[i] * METERS_in;
 		c.xyzt.z = h[i];
-		c = proj_trans(info_new.pj, PJ_FWD, c);
-		if ((ok = proj_errno(info_new.pj)) < 0)
+		c = proj_trans(info_trans.pj, PJ_FWD, c);
+		if ((ok = proj_errno(info_trans.pj)) < 0)
 		    break;
 		/* convert to degrees */
 		x[i] = c.lp.lam * RAD_TO_DEG;
@@ -306,8 +726,8 @@ int pj_do_transform(int count, double *x, double *y, double *h,
 		c.xyzt.x = x[i] * METERS_in;
 		c.xyzt.y = y[i] * METERS_in;
 		c.xyzt.z = h[i];
-		c = proj_trans(info_new.pj, PJ_FWD, c);
-		if ((ok = proj_errno(info_new.pj)) < 0)
+		c = proj_trans(info_trans.pj, PJ_FWD, c);
+		if ((ok = proj_errno(info_trans.pj)) < 0)
 		    break;
 		/* convert to map units */
 		x[i] = c.xy.x / METERS_out;
@@ -317,10 +737,11 @@ int pj_do_transform(int count, double *x, double *y, double *h,
     }
     if (!has_h)
 	G_free(h);
-    proj_destroy(info_new.pj);
+    proj_destroy(info_trans.pj);
 
     if (ok < 0) {
 	G_warning(_("proj_trans() failed: %d"), ok);
+    }
 #else
     METERS_in = info_in->meters;
     METERS_out = info_out->meters;
@@ -363,7 +784,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
     }
+#endif
     return ok;
 }

lib/proj/do_proj_ll.c

@@ -1,236 +0,0 @@
-
-/**
-   \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;
-
-/** 
- * \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 equivalent latitude / longitude system to the 
- * projected co-ordinate system, 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;
-
-#ifdef HAVE_PROJ_H
-    PJ_COORD c;
-
-    ok = 0;
-    if (strncmp(info_in->proj, "ll", 2) == 0) {
-	/* nothing to do */
-	return ok;
-    }
-
-    METERS_in = info_in->meters;
-
-    if (direction == PJ_FWD) {
-	/* from ll to projected */
-
-	/* convert to radians */
-	c.lpzt.lam = (*x) / RAD_TO_DEG;
-	c.lpzt.phi = (*y) / RAD_TO_DEG;
-	c.lpzt.z = 0;
-	c.lpzt.t = 0;
-
-	c = proj_trans(info_in->pj, PJ_FWD, c);
-	ok = proj_errno(info_in->pj);
-
-	/* 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.xyzt.x = *x * METERS_in;
-	c.xyzt.y = *y * METERS_in;
-	c.xyzt.z = 0;
-	c.xyzt.t = 0;
-
-	c = proj_trans(info_in->pj, PJ_INV, c);
-	ok = proj_errno(info_in->pj);
-
-	/* convert to degrees */
-	*x = c.lp.lam * RAD_TO_DEG;
-	*y = c.lp.phi * RAD_TO_DEG;
-    }
-
-    if (ok < 0) {
-	G_warning(_("proj_trans() failed: %d"), ok);
-    }
-#else
-    struct pj_info info_out;
-
-    if ((ok = GPJ_get_equivalent_latlong(&info_out, info_in)) != 1)
-	return ok;
-
-    if (direction == PJ_FWD) {
-	/* from ll to projected */
-	ok = pj_do_proj(x, y, &info_out, info_in);
-    }
-    else {
-	/* from projected to ll */
-	ok = pj_do_proj(x, y, info_in, &info_out);
-    }
-    pj_free(info_out.pj);
-#endif
-
-    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 equivalent latitude / longitude system to the 
- * projected co-ordinate system, 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;
-
-#ifdef HAVE_PROJ_H
-    int has_h = 1;
-    int i;
-    PJ_COORD c;
-
-    ok = 0;
-    if (strncmp(info_in->proj, "ll", 2) == 0) {
-	/* nothing to do */
-	return ok;
-    }
-
-    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.lpzt.lam = x[i] / RAD_TO_DEG;
-	    c.lpzt.phi = y[i] / RAD_TO_DEG;
-	    c.lpzt.z = h[i];
-	    c.lpzt.t = 0;
-
-	    c = proj_trans(info_in->pj, PJ_FWD, c);
-	    if ((ok = proj_errno(info_in->pj)) < 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.xyzt.x = x[i] * METERS_in;
-	    c.xyzt.y = y[i] * METERS_in;
-	    c.xyzt.z = h[i];
-	    c.xyzt.t = 0;
-
-	    c = proj_trans(info_in->pj, PJ_INV, c);
-	    if ((ok = proj_errno(info_in->pj)) < 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);
-
-    if (ok < 0) {
-	G_warning(_("proj_trans() failed: %d"), ok);
-    }
-#else
-    struct pj_info info_out;
-
-    if ((ok = GPJ_get_equivalent_latlong(&info_out, info_in)) != 1)
-	return ok;
-
-    if (direction == PJ_FWD) {
-	/* from ll to projected */
-	ok = pj_do_transform(count, x, y, h, &info_out, info_in);
-    }
-    else {
-	/* from projected to ll */
-	ok = pj_do_transform(count, x, y, h, info_in, &info_out);
-    }
-    pj_free(info_out.pj);
-#endif
-    return ok;
-}
-
-

lib/proj/get_proj.c

@@ -5,9 +5,9 @@
    \brief GProj library - Functions for re-projecting point data
 
    \author Original Author unknown, probably Soil Conservation Service,
-   Eric Miller, Paul Kelly
+   Eric Miller, Paul Kelly, Markus Metz
 
-   (C) 2003-2008 by the GRASS Development Team
+   (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
@@ -80,6 +80,7 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
     info->meters = 1.0;
     info->proj[0] = '\0';
     info->def = NULL;
+    info->pj = NULL;
 
     str = G_find_key_value("meters", in_units_keys);
     if (str != NULL) {
@@ -272,8 +273,8 @@ int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys,
     strcpy(info->def, buffa);
     G_free(opt_in[0]);
 
-    for (i = 0; i < nopt; i++) {
-	sprintf(buffa,  "+%s ", opt_in[0]);
+    for (i = 1; i < nopt; i++) {
+	sprintf(buffa,  "+%s ", opt_in[i]);
 	strcat(info->def, buffa);
 	G_free(opt_in[i]);
     }
@@ -327,6 +328,7 @@ int pj_get_string(struct pj_info *info, char *str)
     info->proj[0] = '\0';
     info->meters = 1.0;
     info->def = NULL;
+    info->pj = NULL;
     
     nopt = 0;
 
@@ -408,8 +410,8 @@ int pj_get_string(struct pj_info *info, char *str)
     strcpy(info->def, buffa);
     G_free(opt_in[0]);
 
-    for (i = 0; i < nopt; i++) {
-	sprintf(buffa,  "+%s ", opt_in[0]);
+    for (i = 1; i < nopt; i++) {
+	sprintf(buffa,  "+%s ", opt_in[i]);
 	strcat(info->def, buffa);
 	G_free(opt_in[i]);
     }
@@ -485,8 +487,8 @@ const char *set_proj_lib(const char *name)
  * \param iproj 'Input' co-ordinate system
  * \param oproj 'Output' co-ordinate system
  * 
- * \return 1 on success, -1 on error (i.e. if PROJ.4 pj_get_def() function
- *         returned NULL for either co-ordinate system)
+ * \return 1 on success, -1 on error (i.e. if the PROJ-style definition 
+ *         is NULL for either co-ordinate system)
  **/
 
 int pj_print_proj_params(const struct pj_info *iproj, const struct pj_info *oproj)