[GRASS-SVN] r72423 - in grass/trunk: include include/defs lib/proj
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Mar 20 11:43:19 PDT 2018
Author: mmetz
Date: 2018-03-20 11:43:19 -0700 (Tue, 20 Mar 2018)
New Revision: 72423
Added:
grass/trunk/lib/proj/do_proj_ll.c
Modified:
grass/trunk/include/defs/gprojects.h
grass/trunk/include/gprojects.h
grass/trunk/lib/proj/convert.c
grass/trunk/lib/proj/do_proj.c
grass/trunk/lib/proj/get_proj.c
Log:
libproj: use new PROJ 5+ API if available
Modified: grass/trunk/include/defs/gprojects.h
===================================================================
--- grass/trunk/include/defs/gprojects.h 2018-03-20 18:40:26 UTC (rev 72422)
+++ grass/trunk/include/defs/gprojects.h 2018-03-20 18:43:19 UTC (rev 72423)
@@ -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,7 +51,7 @@
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 */
@@ -49,5 +58,6 @@
int pj_factors(LP, void *, double, struct FACTORS *);
/* end of copy */
+#endif
#endif
Modified: grass/trunk/include/gprojects.h
===================================================================
--- grass/trunk/include/gprojects.h 2018-03-20 18:40:26 UTC (rev 72422)
+++ grass/trunk/include/gprojects.h 2018-03-20 18:43:19 UTC (rev 72423)
@@ -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 @@
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 @@
int code; /* info as to analytics, see following */
};
/* end of copy */
+#endif
#include <grass/defs/gprojects.h>
Modified: grass/trunk/lib/proj/convert.c
===================================================================
--- grass/trunk/lib/proj/convert.c 2018-03-20 18:40:26 UTC (rev 72422)
+++ grass/trunk/lib/proj/convert.c 2018-03-20 18:43:19 UTC (rev 72423)
@@ -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 @@
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 @@
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 @@
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 @@
char *pszProj = NULL;
const char *pszProjCS = NULL;
char *datum = NULL;
+ char *proj4_unit = NULL;
struct gpj_datum dstruct;
const char *ograttr;
OGRSpatialReferenceH hSRS;
@@ -554,10 +618,14 @@
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);
}
if (!pszProj)
@@ -789,6 +857,13 @@
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 @@
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);
Modified: grass/trunk/lib/proj/do_proj.c
===================================================================
--- grass/trunk/lib/proj/do_proj.c 2018-03-20 18:40:26 UTC (rev 72422)
+++ grass/trunk/lib/proj/do_proj.c 2018-03-20 18:43:19 UTC (rev 72423)
@@ -43,6 +43,20 @@
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,7 +71,7 @@
* \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,
@@ -64,6 +78,71 @@
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 @@
if (ok < 0) {
G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
}
+#endif
return ok;
}
@@ -128,7 +208,7 @@
* \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,11 +216,107 @@
{
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;
+
+ if (h == NULL) {
int i;
h = G_malloc(sizeof *h * count);
@@ -178,6 +354,7 @@
if (ok < 0) {
G_warning(_("pj_transform() failed: %s"), pj_strerrno(ok));
+#endif
}
return ok;
}
Copied: grass/trunk/lib/proj/do_proj_ll.c (from rev 72327, grass/trunk/lib/proj/do_proj.c)
===================================================================
--- grass/trunk/lib/proj/do_proj_ll.c (rev 0)
+++ grass/trunk/lib/proj/do_proj_ll.c 2018-03-20 18:43:19 UTC (rev 72423)
@@ -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 */
+
Modified: grass/trunk/lib/proj/get_proj.c
===================================================================
--- grass/trunk/lib/proj/get_proj.c 2018-03-20 18:40:26 UTC (rev 72422)
+++ grass/trunk/lib/proj/get_proj.c 2018-03-20 18:43:19 UTC (rev 72423)
@@ -31,8 +31,27 @@
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
* PROJ_INFO / PROJ_UNITS-style key-value pairs
@@ -65,7 +84,12 @@
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 @@
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,13 +244,18 @@
}
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]);
@@ -233,11 +262,24 @@
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 @@
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 @@
* 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,9 +354,7 @@
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;
@@ -318,6 +361,13 @@
}
}
+#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 @@
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 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 @@
return -1;
else
return 1;
+#endif
}
/* set_proj_lib()
@@ -396,11 +455,19 @@
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 @@
}
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);
}
More information about the grass-commit
mailing list