GRASS GIS 7 Programmer's Manual  7.2.1(2017)-exported
convert.c
Go to the documentation of this file.
1 
2 /*!
3  \file lib/proj/convert.c
4 
5  \brief GProj Library - Functions for manipulating co-ordinate
6  system representations
7 
8  (C) 2003-2008, 2012 by the GRASS Development Team
9 
10  This program is free software under the GNU General Public License
11  (>=v2). Read the file COPYING that comes with GRASS for details.
12 
13  \author Paul Kelly, Frank Warmerdam
14 */
15 
16 #include <grass/config.h>
17 
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <string.h>
21 #include <math.h>
22 #include <grass/gis.h>
23 #include <grass/gprojects.h>
24 #include <grass/glocale.h>
25 
26 #ifdef HAVE_OGR
27 #include <cpl_csv.h>
28 #include "local_proto.h"
29 
30 /* GRASS relative location of OGR co-ordinate system lookup tables */
31 #define CSVDIR "/etc/proj/ogr_csv"
32 
33 static void DatumNameMassage(char **);
34 #endif
35 
36 /*!
37  * \brief Converts a GRASS co-ordinate system representation to WKT style.
38  *
39  * Takes a GRASS co-ordinate system as specified by two sets of
40  * key/value pairs derived from the PROJ_INFO and PROJ_UNITS files,
41  * and converts it to the 'Well Known Text' format popularised by
42  * proprietary GIS
43  *
44  * \param proj_info Set of GRASS PROJ_INFO key/value pairs
45  * \param proj_units Set of GRASS PROJ_UNIT key/value pairs
46  * \param esri_style boolean Output ESRI-style WKT (Use OSRMorphToESRI()
47  * function provided by OGR library)
48  * \param prettify boolean Use linebreaks and indents to 'prettify' output
49  * WKT string (Use OSRExportToPrettyWkt() function in OGR)
50  *
51  * \return Pointer to a string containing the co-ordinate system in
52  * WKT format
53  * \return NULL on error
54  */
55 char *GPJ_grass_to_wkt(const struct Key_Value *proj_info,
56  const struct Key_Value *proj_units,
57  int esri_style, int prettify)
58 {
59 #ifdef HAVE_OGR
60  OGRSpatialReferenceH hSRS;
61  char *wkt, *local_wkt;
62 
63  hSRS = GPJ_grass_to_osr(proj_info, proj_units);
64 
65  if (hSRS == NULL)
66  return NULL;
67 
68  if (esri_style)
69  OSRMorphToESRI(hSRS);
70 
71  if (prettify)
72  OSRExportToPrettyWkt(hSRS, &wkt, 0);
73  else
74  OSRExportToWkt(hSRS, &wkt);
75 
76  local_wkt = G_store(wkt);
77  CPLFree(wkt);
78  OSRDestroySpatialReference(hSRS);
79 
80  return local_wkt;
81 #else
82  G_warning(_("GRASS is not compiled with OGR support"));
83  return NULL;
84 #endif
85 }
86 
87 #ifdef HAVE_OGR
88 /*!
89  * \brief Converts a GRASS co-ordinate system to an OGRSpatialReferenceH object.
90  *
91  * \param proj_info Set of GRASS PROJ_INFO key/value pairs
92  * \param proj_units Set of GRASS PROJ_UNIT key/value pairs
93  *
94  * \return OGRSpatialReferenceH object representing the co-ordinate system
95  * defined by proj_info and proj_units or NULL if it fails
96  */
97 OGRSpatialReferenceH GPJ_grass_to_osr(const struct Key_Value * proj_info,
98  const struct Key_Value * proj_units)
99 {
100  struct pj_info pjinfo;
101  char *proj4, *proj4mod, *wkt, *modwkt, *startmod, *lastpart;
102  OGRSpatialReferenceH hSRS, hSRS2;
103  OGRErr errcode;
104  struct gpj_datum dstruct;
105  struct gpj_ellps estruct;
106  size_t len;
107  const char *ellpskv, *unit, *unfact;
108  char *ellps, *ellpslong, *datum, *params, *towgs84, *datumlongname,
109  *start, *end;
110  const char *sysname, *osrunit, *osrunfact;
111  double a, es, rf;
112  int haveparams = 0;
113 
114  if ((proj_info == NULL) || (proj_units == NULL))
115  return NULL;
116 
117  hSRS = OSRNewSpatialReference(NULL);
118 
119  if (pj_get_kv(&pjinfo, proj_info, proj_units) < 0) {
120  G_warning(_("Unable parse GRASS PROJ_INFO file"));
121  return NULL;
122  }
123 
124  if ((proj4 = pj_get_def(pjinfo.pj, 0)) == NULL) {
125  G_warning(_("Unable get PROJ.4-style parameter string"));
126  return NULL;
127  }
128  pj_free(pjinfo.pj);
129 
130  unit = G_find_key_value("unit", proj_units);
131  unfact = G_find_key_value("meters", proj_units);
132  if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
133  G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
134  else
135  proj4mod = G_store(proj4);
136  pj_dalloc(proj4);
137 
138  if ((errcode = OSRImportFromProj4(hSRS, proj4mod)) != OGRERR_NONE) {
139  G_warning(_("OGR can't parse PROJ.4-style parameter string: "
140  "%s (OGR Error code was %d)"), proj4mod, errcode);
141  return NULL;
142  }
143  G_free(proj4mod);
144 
145  /* this messes up PROJCS versus GEOGCS!
146  sysname = G_find_key_value("name", proj_info);
147  if (sysname)
148  OSRSetProjCS(hSRS, sysname);
149  */
150 
151  if ((errcode = OSRExportToWkt(hSRS, &wkt)) != OGRERR_NONE) {
152  G_warning(_("OGR can't get WKT-style parameter string "
153  "(OGR Error code was %d)"), errcode);
154  return NULL;
155  }
156 
157  ellpskv = G_find_key_value("ellps", proj_info);
158  GPJ__get_ellipsoid_params(proj_info, &a, &es, &rf);
159  haveparams = GPJ__get_datum_params(proj_info, &datum, &params);
160 
161  if (ellpskv != NULL)
162  ellps = G_store(ellpskv);
163  else
164  ellps = NULL;
165 
166  if ((datum == NULL) || (GPJ_get_datum_by_name(datum, &dstruct) < 0)) {
167  datumlongname = G_store("unknown");
168  if (ellps == NULL)
169  ellps = G_store("unnamed");
170  }
171  else {
172  datumlongname = G_store(dstruct.longname);
173  if (ellps == NULL)
174  ellps = G_store(dstruct.ellps);
175  GPJ_free_datum(&dstruct);
176  }
177  G_debug(3, "GPJ_grass_to_osr: datum: <%s>", datum);
178  G_free(datum);
179  if (GPJ_get_ellipsoid_by_name(ellps, &estruct) > 0) {
180  ellpslong = G_store(estruct.longname);
181  DatumNameMassage(&ellpslong);
182  GPJ_free_ellps(&estruct);
183  }
184  else
185  ellpslong = G_store(ellps);
186 
187  startmod = strstr(wkt, "GEOGCS");
188  lastpart = strstr(wkt, "PRIMEM");
189  len = strlen(wkt) - strlen(startmod);
190  wkt[len] = '\0';
191  if (haveparams == 2) {
192  /* Only put datum params into the WKT if they were specifically
193  * specified in PROJ_INFO */
194  char *paramkey, *paramvalue;
195 
196  paramkey = strtok(params, "=");
197  paramvalue = params + strlen(paramkey) + 1;
198  if (G_strcasecmp(paramkey, "towgs84") == 0)
199  G_asprintf(&towgs84, ",TOWGS84[%s]", paramvalue);
200  else
201  towgs84 = G_store("");
202  G_free(params);
203  }
204  else
205  towgs84 = G_store("");
206 
207  sysname = OSRGetAttrValue(hSRS, "PROJCS", 0);
208  if (sysname == NULL) {
209  /* Not a projected co-ordinate system */
210  start = G_store("");
211  end = G_store("");
212  }
213  else {
214  if ((strcmp(sysname, "unnamed") == 0) &&
215  (G_find_key_value("name", proj_info) != NULL))
216  G_asprintf(&start, "PROJCS[\"%s\",",
217  G_find_key_value("name", proj_info));
218  else
219  start = G_store(wkt);
220 
221  osrunit = OSRGetAttrValue(hSRS, "UNIT", 0);
222  osrunfact = OSRGetAttrValue(hSRS, "UNIT", 1);
223 
224  if ((unfact == NULL) || (G_strcasecmp(osrunit, "unknown") != 0))
225  end = G_store("");
226  else {
227  char *buff;
228  double unfactf = atof(unfact);
229 
230  G_asprintf(&buff, ",UNIT[\"%s\",", osrunit);
231 
232  startmod = strstr(lastpart, buff);
233  len = strlen(lastpart) - strlen(startmod);
234  lastpart[len] = '\0';
235  G_free(buff);
236 
237  if (unit == NULL)
238  unit = "unknown";
239  G_asprintf(&end, ",UNIT[\"%s\",%.16g]]", unit, unfactf);
240  }
241 
242  }
243  OSRDestroySpatialReference(hSRS);
244  G_asprintf(&modwkt,
245  "%sGEOGCS[\"%s\",DATUM[\"%s\",SPHEROID[\"%s\",%.16g,%.16g]%s],%s%s",
246  start, ellps, datumlongname, ellpslong, a, rf, towgs84,
247  lastpart, end);
248  hSRS2 = OSRNewSpatialReference(modwkt);
249  G_free(modwkt);
250 
251  CPLFree(wkt);
252  G_free(start);
253  G_free(ellps);
254  G_free(datumlongname);
255  G_free(ellpslong);
256  G_free(towgs84);
257  G_free(end);
258 
259  return hSRS2;
260 }
261 
262 /*!
263  * \brief Converts an OGRSpatialReferenceH object to a GRASS co-ordinate system.
264  *
265  * \param cellhd Pointer to a GRASS Cell_head structure that will have its
266  * projection-related members populated with appropriate values
267  * \param projinfo Pointer to a pointer which will have a GRASS Key_Value
268  * structure allocated containing a set of GRASS PROJ_INFO values
269  * \param projunits Pointer to a pointer which will have a GRASS Key_Value
270  * structure allocated containing a set of GRASS PROJ_UNITS values
271  * \param hSRS OGRSpatialReferenceH object containing the co-ordinate
272  * system to be converted
273  * \param datumtrans Index number of datum parameter set to use, 0 to leave
274  * unspecified
275  *
276  * \return 2 if a projected or lat/long co-ordinate system has been
277  * defined
278  * \return 1 if an unreferenced XY co-ordinate system has
279  * been defined
280  */
281 int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
282  struct Key_Value **projunits, OGRSpatialReferenceH hSRS,
283  int datumtrans)
284 {
285  struct Key_Value *temp_projinfo;
286  char *pszProj4 = NULL, *pszRemaining;
287  char *pszProj = NULL;
288  const char *pszProjCS = NULL;
289  char *datum = NULL;
290  struct gpj_datum dstruct;
291 
292  if (hSRS == NULL)
293  goto default_to_xy;
294 
295  /* Set finder function for locating OGR csv co-ordinate system tables */
296  /* SetCSVFilenameHook(GPJ_set_csv_loc); */
297 
298  /* Hopefully this doesn't do any harm if it wasn't in ESRI format
299  * to start with... */
300  OSRMorphFromESRI(hSRS);
301 
302  /* -------------------------------------------------------------------- */
303  /* Set cellhd for well known coordinate systems. */
304  /* -------------------------------------------------------------------- */
305  if (!OSRIsGeographic(hSRS) && !OSRIsProjected(hSRS))
306  goto default_to_xy;
307 
308  if (cellhd) {
309  int bNorth;
310 
311  if (OSRIsGeographic(hSRS)) {
312  cellhd->proj = PROJECTION_LL;
313  cellhd->zone = 0;
314  }
315  else if (OSRGetUTMZone(hSRS, &bNorth) != 0) {
316  cellhd->proj = PROJECTION_UTM;
317  cellhd->zone = OSRGetUTMZone(hSRS, &bNorth);
318  if (!bNorth)
319  cellhd->zone *= -1;
320  }
321  else {
322  cellhd->proj = PROJECTION_OTHER;
323  cellhd->zone = 0;
324  }
325  }
326 
327  /* -------------------------------------------------------------------- */
328  /* Get the coordinate system definition in PROJ.4 format. */
329  /* -------------------------------------------------------------------- */
330  if (OSRExportToProj4(hSRS, &pszProj4) != OGRERR_NONE)
331  goto default_to_xy;
332 
333  /* -------------------------------------------------------------------- */
334  /* Parse the PROJ.4 string into key/value pairs. Do a bit of */
335  /* extra work to "GRASSify" the result. */
336  /* -------------------------------------------------------------------- */
337  temp_projinfo = G_create_key_value();
338 
339  /* Create "local" copy of proj4 string so we can modify and free it
340  * using GRASS functions */
341  pszRemaining = G_store(pszProj4);
342  CPLFree(pszProj4);
343  pszProj4 = pszRemaining;
344  while ((pszRemaining = strstr(pszRemaining, "+")) != NULL) {
345  char *pszToken, *pszValue;
346 
347  pszRemaining++;
348 
349  /* Advance pszRemaining to end of this token[=value] pair */
350  pszToken = pszRemaining;
351  while (*pszRemaining != ' ' && *pszRemaining != '\0')
352  pszRemaining++;
353 
354  if (*pszRemaining == ' ') {
355  *pszRemaining = '\0';
356  pszRemaining++;
357  }
358 
359  /* parse token, and value */
360  if (strstr(pszToken, "=") != NULL) {
361  pszValue = strstr(pszToken, "=");
362  *pszValue = '\0';
363  pszValue++;
364  }
365  else
366  pszValue = "defined";
367 
368  /* projection name */
369  if (G_strcasecmp(pszToken, "proj") == 0) {
370  /* The ll projection is known as longlat in PROJ.4 */
371  if (G_strcasecmp(pszValue, "longlat") == 0)
372  pszValue = "ll";
373 
374  pszProj = pszValue;
375  }
376 
377  /* Ellipsoid and datum handled separately below */
378  if (G_strcasecmp(pszToken, "ellps") == 0
379  || G_strcasecmp(pszToken, "a") == 0
380  || G_strcasecmp(pszToken, "b") == 0
381  || G_strcasecmp(pszToken, "es") == 0
382  || G_strcasecmp(pszToken, "rf") == 0
383  || G_strcasecmp(pszToken, "datum") == 0)
384  continue;
385 
386  /* We will handle units separately */
387  if (G_strcasecmp(pszToken, "to_meter") == 0
388  || G_strcasecmp(pszToken, "units") == 0)
389  continue;
390 
391  G_set_key_value(pszToken, pszValue, temp_projinfo);
392  }
393  if (!pszProj)
394  G_warning(_("No projection name! Projection parameters likely to be meaningless."));
395 
396  *projinfo = G_create_key_value();
397 
398  /* -------------------------------------------------------------------- */
399  /* Derive the user name for the coordinate system. */
400  /* -------------------------------------------------------------------- */
401  pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
402  if (!pszProjCS)
403  pszProjCS = OSRGetAttrValue(hSRS, "GEOGCS", 0);
404 
405  if (pszProjCS) {
406  G_set_key_value("name", pszProjCS, *projinfo);
407  }
408  else if (pszProj) {
409  char path[4095];
410  char name[80];
411 
412  /* use name of the projection as name for the coordinate system */
413 
414  sprintf(path, "%s/etc/proj/projections", G_gisbase());
415  if (G_lookup_key_value_from_file(path, pszProj, name, sizeof(name)) >
416  0)
417  G_set_key_value("name", name, *projinfo);
418  else
419  G_set_key_value("name", pszProj, *projinfo);
420  }
421 
422 
423  /* -------------------------------------------------------------------- */
424  /* Find the GRASS datum name and choose parameters either */
425  /* interactively or not. */
426  /* -------------------------------------------------------------------- */
427 
428  {
429  const char *pszDatumNameConst = OSRGetAttrValue(hSRS, "DATUM", 0);
430  struct datum_list *list, *listhead;
431  char *dum1, *dum2, *pszDatumName;
432  int paramspresent =
433  GPJ__get_datum_params(temp_projinfo, &dum1, &dum2);
434 
435  if (pszDatumNameConst) {
436  /* Need to make a new copy of the string so we don't mess
437  * around with the memory inside the OGRSpatialReferenceH? */
438 
439  pszDatumName = G_store(pszDatumNameConst);
440  DatumNameMassage(&pszDatumName);
441  G_debug(3, "GPJ_osr_to_grass: pszDatumNameConst: <%s>", pszDatumName);
442 
443  list = listhead = read_datum_table();
444 
445  while (list != NULL) {
446  if (G_strcasecmp(pszDatumName, list->longname) == 0) {
447  datum = G_store(list->name);
448  break;
449  }
450  list = list->next;
451  }
452  free_datum_list(listhead);
453 
454  if (datum == NULL) {
455  if (paramspresent < 2)
456  /* Only give warning if no parameters present */
457  G_warning(_("Datum <%s> not recognised by GRASS and no parameters found"),
458  pszDatumName);
459  }
460  else {
461  G_set_key_value("datum", datum, *projinfo);
462 
463  if (paramspresent < 2) {
464  /* If no datum parameters were imported from the OSR
465  * object then we should use the set specified by datumtrans */
466  char *params, *chosenparams = NULL;
467  int paramsets;
468 
469  paramsets =
470  GPJ_get_default_datum_params_by_name(datum, &params);
471 
472  if (paramsets < 0)
473  G_warning(_("Datum <%s> apparently recognised by GRASS but no parameters found. "
474  "You may want to look into this."), datum);
475  else if (datumtrans > paramsets) {
476 
477  G_warning(_("Invalid transformation number %d; valid range is 1 to %d. "
478  "Leaving datum transform parameters unspecified."),
479  datumtrans, paramsets);
480  datumtrans = 0;
481  }
482 
483  if (paramsets > 0) {
484  struct gpj_datum_transform_list *list, *old;
485 
486  list = GPJ_get_datum_transform_by_name(datum);
487 
488  if (list != NULL) {
489  do {
490  if (list->count == datumtrans)
491  chosenparams = G_store(list->params);
492  old = list;
493  list = list->next;
495  } while (list != NULL);
496  }
497  }
498 
499  if (chosenparams != NULL) {
500  char *paramkey, *paramvalue;
501 
502  paramkey = strtok(chosenparams, "=");
503  paramvalue = chosenparams + strlen(paramkey) + 1;
504  G_set_key_value(paramkey, paramvalue, *projinfo);
505  G_free(chosenparams);
506  }
507 
508  if (paramsets > 0)
509  G_free(params);
510  }
511 
512  }
513  G_free(pszDatumName);
514  }
515  }
516 
517  /* -------------------------------------------------------------------- */
518  /* Determine an appropriate GRASS ellipsoid name if possible, or */
519  /* else just put a and es values into PROJ_INFO */
520  /* -------------------------------------------------------------------- */
521 
522  if ((datum != NULL) && (GPJ_get_datum_by_name(datum, &dstruct) > 0)) {
523  /* Use ellps name associated with datum */
524  G_set_key_value("ellps", dstruct.ellps, *projinfo);
525  GPJ_free_datum(&dstruct);
526  G_free(datum);
527  }
528  else {
529  /* If we can't determine the ellipsoid from the datum, derive it
530  * directly from "SPHEROID" parameters in WKT */
531  const char *pszSemiMajor = OSRGetAttrValue(hSRS, "SPHEROID", 1);
532  const char *pszInvFlat = OSRGetAttrValue(hSRS, "SPHEROID", 2);
533 
534  if (pszSemiMajor != NULL && pszInvFlat != NULL) {
535  char *ellps = NULL;
536  struct ellps_list *list, *listhead;
537  double a = atof(pszSemiMajor), invflat = atof(pszInvFlat), flat;
538  double es;
539 
540  /* Allow for incorrect WKT describing a sphere where InvFlat
541  * is given as 0 rather than inf */
542  if (invflat > 0)
543  flat = 1 / invflat;
544  else
545  flat = 0;
546 
547  es = flat * (2.0 - flat);
548 
549  list = listhead = read_ellipsoid_table(0);
550 
551  while (list != NULL) {
552  /* Try and match a and es against GRASS defined ellipsoids;
553  * accept first one that matches. These numbers were found
554  * by trial and error and could be fine-tuned, or possibly
555  * a direct comparison of IEEE floating point values used. */
556  if ((a == list->a || fabs(a - list->a) < 0.1 || fabs(1 - a / list->a) < 0.0000001) && ((es == 0 && list->es == 0) || /* Special case for sphere */
557  (invflat == list->rf || fabs(invflat - list->rf) < 0.0000001))) {
558  ellps = G_store(list->name);
559  break;
560  }
561  list = list->next;
562  }
563  if (listhead != NULL)
564  free_ellps_list(listhead);
565 
566  if (ellps == NULL) {
567  /* If we weren't able to find a matching ellps name, set
568  * a and es values directly from WKT-derived data */
569  char es_str[100];
570 
571  G_set_key_value("a", (char *)pszSemiMajor, *projinfo);
572 
573  sprintf(es_str, "%.16g", es);
574  G_set_key_value("es", es_str, *projinfo);
575  }
576  else {
577  /* else specify the GRASS ellps name for readability */
578  G_set_key_value("ellps", ellps, *projinfo);
579  G_free(ellps);
580  }
581 
582  }
583 
584  }
585 
586  /* -------------------------------------------------------------------- */
587  /* Finally append the detailed projection parameters to the end */
588  /* -------------------------------------------------------------------- */
589 
590  {
591  int i;
592 
593  for (i = 0; i < temp_projinfo->nitems; i++)
594  G_set_key_value(temp_projinfo->key[i],
595  temp_projinfo->value[i], *projinfo);
596 
597  G_free_key_value(temp_projinfo);
598  }
599 
600  G_free(pszProj4);
601 
602  /* -------------------------------------------------------------------- */
603  /* Set the linear units. */
604  /* -------------------------------------------------------------------- */
605  *projunits = G_create_key_value();
606 
607  if (OSRIsGeographic(hSRS)) {
608  /* We assume degrees ... someday we will be wrong! */
609  G_set_key_value("unit", "degree", *projunits);
610  G_set_key_value("units", "degrees", *projunits);
611  G_set_key_value("meters", "1.0", *projunits);
612  }
613  else {
614  char szFormatBuf[256];
615  char *pszUnitsName = NULL;
616  double dfToMeters;
617  char *pszUnitsPlural, *pszStringEnd;
618 
619  dfToMeters = OSRGetLinearUnits(hSRS, &pszUnitsName);
620 
621  /* Workaround for the most obvious case when unit name is unknown */
622  if ((G_strcasecmp(pszUnitsName, "unknown") == 0) &&
623  (dfToMeters == 1.))
624  G_asprintf(&pszUnitsName, "meter");
625 
626  if ((G_strcasecmp(pszUnitsName, "metre") == 0))
627  G_asprintf(&pszUnitsName, "meter");
628  if ((G_strcasecmp(pszUnitsName, "kilometre") == 0))
629  G_asprintf(&pszUnitsName, "kilometer");
630 
631  G_set_key_value("unit", pszUnitsName, *projunits);
632 
633  /* Attempt at plural formation (WKT format doesn't store plural
634  * form of unit name) */
635  pszUnitsPlural = G_malloc(strlen(pszUnitsName) + 3);
636  strcpy(pszUnitsPlural, pszUnitsName);
637  pszStringEnd = pszUnitsPlural + strlen(pszUnitsPlural) - 4;
638  if (G_strcasecmp(pszStringEnd, "foot") == 0) {
639  /* Special case for foot - change two o's to e's */
640  pszStringEnd[1] = 'e';
641  pszStringEnd[2] = 'e';
642  }
643  else if (G_strcasecmp(pszStringEnd, "inch") == 0) {
644  /* Special case for inch - add es */
645  pszStringEnd[4] = 'e';
646  pszStringEnd[5] = 's';
647  pszStringEnd[6] = '\0';
648  }
649  else {
650  /* For anything else add an s at the end */
651  pszStringEnd[4] = 's';
652  pszStringEnd[5] = '\0';
653  }
654 
655  G_set_key_value("units", pszUnitsPlural, *projunits);
656  G_free(pszUnitsPlural);
657 
658  sprintf(szFormatBuf, "%.16g", dfToMeters);
659  G_set_key_value("meters", szFormatBuf, *projunits);
660 
661  }
662 
663  return 2;
664 
665  /* -------------------------------------------------------------------- */
666  /* Fallback to returning an ungeoreferenced definition. */
667  /* -------------------------------------------------------------------- */
668  default_to_xy:
669  if (cellhd != NULL) {
670  cellhd->proj = PROJECTION_XY;
671  cellhd->zone = 0;
672  }
673 
674  *projinfo = NULL;
675  *projunits = NULL;
676 
677  return 1;
678 }
679 #endif
680 
681 /*!
682  * \brief Converts a WKT projection description to a GRASS co-ordinate system.
683  *
684  * \param cellhd Pointer to a GRASS Cell_head structure that will have its
685  * projection-related members populated with appropriate values
686  * \param projinfo Pointer to a pointer which will have a GRASS Key_Value
687  * structure allocated containing a set of GRASS PROJ_INFO values
688  * \param projunits Pointer to a pointer which will have a GRASS Key_Value
689  * structure allocated containing a set of GRASS PROJ_UNITS values
690  * \param wkt Well-known Text (WKT) description of the co-ordinate
691  * system to be converted
692  * \param datumtrans Index number of datum parameter set to use, 0 to leave
693  * unspecified
694  *
695  * \return 2 if a projected or lat/long co-ordinate system has been
696  * defined
697  * \return 1 if an unreferenced XY co-ordinate system has
698  * been defined
699  * \return -1 on error
700  */
701 int GPJ_wkt_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
702  struct Key_Value **projunits, const char *wkt,
703  int datumtrans)
704 {
705 #ifdef HAVE_OGR
706  int retval;
707 
708  if (wkt == NULL)
709  retval =
710  GPJ_osr_to_grass(cellhd, projinfo, projunits, NULL, datumtrans);
711  else {
712  OGRSpatialReferenceH hSRS;
713 
714  /* Set finder function for locating OGR csv co-ordinate system tables */
715  /* SetCSVFilenameHook(GPJ_set_csv_loc); */
716 
717  hSRS = OSRNewSpatialReference(wkt);
718  retval =
719  GPJ_osr_to_grass(cellhd, projinfo, projunits, hSRS, datumtrans);
720  OSRDestroySpatialReference(hSRS);
721  }
722 
723  return retval;
724 #else
725  return -1;
726 #endif
727 }
728 
729 #ifdef HAVE_OGR
730 /* GPJ_set_csv_loc()
731  * 'finder function' for use with OGR SetCSVFilenameHook() function */
732 
733 const char *GPJ_set_csv_loc(const char *name)
734 {
735  const char *gisbase = G_gisbase();
736  static char *buf = NULL;
737 
738  if (buf != NULL)
739  G_free(buf);
740 
741  G_asprintf(&buf, "%s%s/%s", gisbase, CSVDIR, name);
742 
743  return buf;
744 }
745 
746 
747 /* The list below is only for files that use a non-standard name for a
748  * datum that is already supported in GRASS. The number of entries must be even;
749  * they are all in pairs. The first one in the pair is the non-standard name;
750  * the second is the GRASS/GDAL name. If a name appears more than once (as for
751  * European_Terrestrial_Reference_System_1989) then it means there was more
752  * than one non-standard name for it that needs to be accounted for.
753  *
754  * N.B. The order of these pairs is different from that in
755  * ogr/ogrfromepsg.cpp in the GDAL source tree! GRASS uses the EPSG
756  * names in its WKT representation except WGS_1984 and WGS_1972 as
757  * these shortened versions seem to be standard.
758  * Below order:
759  * the equivalent name comes first in the pair, and
760  * the EPSG name (as used in the GRASS datum.table file) comes second.
761  *
762  * The datum parameters are stored in
763  * ../gis/datum.table # 3 parameters
764  * ../gis/datumtransform.table # 7 parameters (requires entry in datum.table)
765  *
766  * Hint: use GDAL's "testepsg" to identify the canonical name, e.g.
767  * testepsg epsg:4674
768  */
769 
770 static const char *papszDatumEquiv[] = {
771  "Militar_Geographische_Institute",
772  "Militar_Geographische_Institut",
773  "World_Geodetic_System_1984",
774  "WGS_1984",
775  "World_Geodetic_System_1972",
776  "WGS_1972",
777  "European_Terrestrial_Reference_System_89",
778  "European_Terrestrial_Reference_System_1989",
779  "European_Reference_System_1989",
780  "European_Terrestrial_Reference_System_1989",
781  "ETRS_1989",
782  "European_Terrestrial_Reference_System_1989",
783  "ETRS89",
784  "European_Terrestrial_Reference_System_1989",
785  "ETRF_1989",
786  "European_Terrestrial_Reference_System_1989",
787  "NZGD_2000",
788  "New_Zealand_Geodetic_Datum_2000",
789  "Monte_Mario_Rome",
790  "Monte_Mario",
791  "MONTROME",
792  "Monte_Mario",
793  "Campo_Inchauspe_1969",
794  "Campo_Inchauspe",
795  "S_JTSK",
796  "System_Jednotne_Trigonometricke_Site_Katastralni",
797  "S_JTSK_Ferro",
798  "Militar_Geographische_Institut",
799  "Potsdam_Datum_83",
800  "Deutsches_Hauptdreiecksnetz",
801  "South_American_1969",
802  "South_American_Datum_1969",
803  "ITRF_1992",
804  "ITRF92",
805  NULL
806 };
807 
808 /************************************************************************/
809 /* OGREPSGDatumNameMassage() */
810 /* */
811 /* Massage an EPSG datum name into WMT format. Also transform */
812 /* specific exception cases into WKT versions. */
813 
814 /************************************************************************/
815 
816 static void DatumNameMassage(char **ppszDatum)
817 {
818  int i, j;
819  char *pszDatum = *ppszDatum;
820 
821  G_debug(3, "DatumNameMassage: Raw string found <%s>", (char *)pszDatum);
822  /* -------------------------------------------------------------------- */
823  /* Translate non-alphanumeric values to underscores. */
824  /* -------------------------------------------------------------------- */
825  for (i = 0; pszDatum[i] != '\0'; i++) {
826  if (!(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z')
827  && !(pszDatum[i] >= 'a' && pszDatum[i] <= 'z')
828  && !(pszDatum[i] >= '0' && pszDatum[i] <= '9')) {
829  pszDatum[i] = '_';
830  }
831  }
832 
833  /* -------------------------------------------------------------------- */
834  /* Remove repeated and trailing underscores. */
835  /* -------------------------------------------------------------------- */
836  for (i = 1, j = 0; pszDatum[i] != '\0'; i++) {
837  if (pszDatum[j] == '_' && pszDatum[i] == '_')
838  continue;
839 
840  pszDatum[++j] = pszDatum[i];
841  }
842  if (pszDatum[j] == '_')
843  pszDatum[j] = '\0';
844  else
845  pszDatum[j + 1] = '\0';
846 
847  /* -------------------------------------------------------------------- */
848  /* Search for datum equivalences. Specific massaged names get */
849  /* mapped to OpenGIS specified names. */
850  /* -------------------------------------------------------------------- */
851  G_debug(3, "DatumNameMassage: Search for datum equivalences of <%s>", (char *)pszDatum);
852  for (i = 0; papszDatumEquiv[i] != NULL; i += 2) {
853  if (EQUAL(*ppszDatum, papszDatumEquiv[i])) {
854  G_free(*ppszDatum);
855  *ppszDatum = G_store(papszDatumEquiv[i + 1]);
856  break;
857  }
858  }
859 }
860 
861 #endif /* HAVE_OGR */
int GPJ_wkt_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo, struct Key_Value **projunits, const char *wkt, int datumtrans)
Converts a WKT projection description to a GRASS co-ordinate system.
Definition: convert.c:701
int G_strcasecmp(const char *x, const char *y)
String compare ignoring case (upper or lower)
Definition: strings.c:46
const char * G_find_key_value(const char *key, const struct Key_Value *kv)
Find given key (case sensitive)
Definition: key_value1.c:84
void GPJ_free_datum(struct gpj_datum *dstruct)
Free the memory used for the strings in a gpj_datum struct.
Definition: proj/datum.c:403
int GPJ__get_datum_params(const struct Key_Value *projinfo, char **datumname, char **params)
Extract the datum transformation-related parameters from a set of general PROJ_INFO parameters...
Definition: proj/datum.c:173
void GPJ_free_datum_transform(struct gpj_datum_transform_list *item)
Free the memory used by a gpj_datum_transform_list struct.
Definition: proj/datum.c:326
struct ellps_list * read_ellipsoid_table(int fatal)
Definition: ellipse.c:224
int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo, struct Key_Value **projunits, OGRSpatialReferenceH hSRS, int datumtrans)
Converts an OGRSpatialReferenceH object to a GRASS co-ordinate system.
Definition: convert.c:281
char * GPJ_grass_to_wkt(const struct Key_Value *proj_info, const struct Key_Value *proj_units, int esri_style, int prettify)
Converts a GRASS co-ordinate system representation to WKT style.
Definition: convert.c:55
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:86
int G_asprintf(char **out, const char *fmt,...)
Definition: asprintf.c:70
#define NULL
Definition: ccmath.h:32
struct gpj_datum_transform_list * GPJ_get_datum_transform_by_name(const char *inputname)
Internal function to find all possible sets of transformation parameters for a particular datum...
Definition: proj/datum.c:237
int GPJ_get_datum_by_name(const char *name, struct gpj_datum *dstruct)
Look up a string in datum.table file to see if it is a valid datum name and if so place its informati...
Definition: proj/datum.c:37
void G_free_key_value(struct Key_Value *kv)
Free allocated Key_Value structure.
Definition: key_value1.c:103
int GPJ__get_ellipsoid_params(const struct Key_Value *proj_keys, double *a, double *e2, double *rf)
Get the ellipsoid parameters from proj keys structure.
Definition: ellipse.c:73
int G_lookup_key_value_from_file(const char *file, const char *key, char value[], int n)
Look up for key in file.
Definition: key_value4.c:48
#define EQUAL(a, b)
Definition: gsdrape.c:64
void free_datum_list(struct datum_list *dstruct)
Free the memory used by a datum_list linked list structure.
Definition: proj/datum.c:417
struct list * list
Definition: read_list.c:24
#define CSVDIR
Definition: convert.c:31
const char * GPJ_set_csv_loc(const char *name)
Definition: convert.c:733
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
void GPJ_free_ellps(struct gpj_ellps *estruct)
Free ellipsoid data structure.
Definition: ellipse.c:309
OGRSpatialReferenceH GPJ_grass_to_osr(const struct Key_Value *proj_info, const struct Key_Value *proj_units)
Converts a GRASS co-ordinate system to an OGRSpatialReferenceH object.
Definition: convert.c:97
struct datum_list * read_datum_table(void)
Read the current GRASS datum.table from disk and store in memory.
Definition: proj/datum.c:345
int pj_get_kv(struct pj_info *info, const struct Key_Value *in_proj_keys, const struct Key_Value *in_units_keys)
Create a pj_info struct Co-ordinate System definition from a set of PROJ_INFO / PROJ_UNITS-style key-...
Definition: get_proj.c:59
int GPJ_get_default_datum_params_by_name(const char *name, char **params)
"Last resort" function to retrieve a "default" set of datum parameters for a datum (N...
Definition: proj/datum.c:85
Definition: path.h:16
void G_set_key_value(const char *key, const char *value, struct Key_Value *kv)
Set value for given key.
Definition: key_value1.c:38
int GPJ_get_ellipsoid_by_name(const char *name, struct gpj_ellps *estruct)
Looks up ellipsoid in ellipsoid table and returns the a, e2 parameters for the ellipsoid.
Definition: ellipse.c:160
const char * name
Definition: named_colr.c:7
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition: gisbase.c:41
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
struct Key_Value * G_create_key_value(void)
Allocate and initialize Key_Value structure.
Definition: key_value1.c:23
void free_ellps_list(struct ellps_list *elist)
Definition: ellipse.c:316