GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71917
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2017 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 static char *grass_to_wkt(const struct Key_Value *proj_info,
37  const struct Key_Value *proj_units,
38  const struct Key_Value *proj_epsg,
39  int esri_style, int prettify)
40 {
41 #ifdef HAVE_OGR
42  OGRSpatialReferenceH hSRS;
43  char *wkt, *local_wkt;
44 
45  hSRS = GPJ_grass_to_osr2(proj_info, proj_units, proj_epsg);
46 
47  if (hSRS == NULL)
48  return NULL;
49 
50  if (esri_style)
51  OSRMorphToESRI(hSRS);
52 
53  if (prettify)
54  OSRExportToPrettyWkt(hSRS, &wkt, 0);
55  else
56  OSRExportToWkt(hSRS, &wkt);
57 
58  local_wkt = G_store(wkt);
59  CPLFree(wkt);
60  OSRDestroySpatialReference(hSRS);
61 
62  return local_wkt;
63 #else
64  G_warning(_("GRASS is not compiled with OGR support"));
65  return NULL;
66 #endif
67 }
68 
69 /*!
70  * \brief Converts a GRASS co-ordinate system representation to WKT style.
71  *
72  * Takes a GRASS co-ordinate system as specified by two sets of
73  * key/value pairs derived from the PROJ_INFO and PROJ_UNITS files,
74  * and converts it to the 'Well Known Text' format.
75  *
76  * \param proj_info Set of GRASS PROJ_INFO key/value pairs
77  * \param proj_units Set of GRASS PROJ_UNIT key/value pairs
78  * \param esri_style boolean Output ESRI-style WKT (Use OSRMorphToESRI()
79  * function provided by OGR library)
80  * \param prettify boolean Use linebreaks and indents to 'prettify' output
81  * WKT string (Use OSRExportToPrettyWkt() function in OGR)
82  *
83  * \return Pointer to a string containing the co-ordinate system in
84  * WKT format
85  * \return NULL on error
86  */
87 char *GPJ_grass_to_wkt(const struct Key_Value *proj_info,
88  const struct Key_Value *proj_units,
89  int esri_style, int prettify)
90 {
91  return grass_to_wkt(proj_info, proj_units, NULL, esri_style, prettify);
92 }
93 
94 /*!
95  * \brief Converts a GRASS co-ordinate system representation to WKT
96  * style. EPSG code is preferred if available.
97  *
98  * Takes a GRASS co-ordinate system as specified key/value pairs
99  * derived from the PROJ_EPSG file. TOWGS84 parameter is scanned
100  * from PROJ_INFO file and appended to co-ordinate system definition
101  * imported from EPSG code by GDAL library. PROJ_UNITS file is
102  * ignored. The function converts it to the 'Well Known Text' format.
103  *
104  * \todo Merge with GPJ_grass_to_wkt() in GRASS 8.
105  *
106  * \param proj_info Set of GRASS PROJ_INFO key/value pairs
107  * \param proj_units Set of GRASS PROJ_UNIT key/value pairs
108  * \param proj_epsg Set of GRASS PROJ_EPSG key/value pairs
109  * \param esri_style boolean Output ESRI-style WKT (Use OSRMorphToESRI()
110  * function provided by OGR library)
111  * \param prettify boolean Use linebreaks and indents to 'prettify' output
112  * WKT string (Use OSRExportToPrettyWkt() function in OGR)
113  *
114  * \return Pointer to a string containing the co-ordinate system in
115  * WKT format
116  * \return NULL on error
117  */
118 char *GPJ_grass_to_wkt2(const struct Key_Value *proj_info,
119  const struct Key_Value *proj_units,
120  const struct Key_Value *proj_epsg,
121  int esri_style, int prettify)
122 {
123  return grass_to_wkt(proj_info, proj_units, proj_epsg, esri_style, prettify);
124 }
125 
126 #ifdef HAVE_OGR
127 /*!
128  * \brief Converts a GRASS co-ordinate system to an OGRSpatialReferenceH object.
129  *
130  * \param proj_info Set of GRASS PROJ_INFO key/value pairs
131  * \param proj_units Set of GRASS PROJ_UNIT key/value pairs
132  *
133  * \return OGRSpatialReferenceH object representing the co-ordinate system
134  * defined by proj_info and proj_units or NULL if it fails
135  */
136 OGRSpatialReferenceH GPJ_grass_to_osr(const struct Key_Value * proj_info,
137  const struct Key_Value * proj_units)
138 {
139  struct pj_info pjinfo;
140  char *proj4, *proj4mod, *wkt, *modwkt, *startmod, *lastpart;
141  OGRSpatialReferenceH hSRS, hSRS2;
142  OGRErr errcode;
143  struct gpj_datum dstruct;
144  struct gpj_ellps estruct;
145  size_t len;
146  const char *ellpskv, *unit, *unfact;
147  char *ellps, *ellpslong, *datum, *params, *towgs84, *datumlongname,
148  *start, *end;
149  const char *sysname, *osrunit, *osrunfact;
150  double a, es, rf;
151  int haveparams = 0;
152 
153  if ((proj_info == NULL) || (proj_units == NULL))
154  return NULL;
155 
156  hSRS = OSRNewSpatialReference(NULL);
157 
158  if (pj_get_kv(&pjinfo, proj_info, proj_units) < 0) {
159  G_warning(_("Unable parse GRASS PROJ_INFO file"));
160  return NULL;
161  }
162 
163  if ((proj4 = pj_get_def(pjinfo.pj, 0)) == NULL) {
164  G_warning(_("Unable get PROJ.4-style parameter string"));
165  return NULL;
166  }
167  pj_free(pjinfo.pj);
168 
169  unit = G_find_key_value("unit", proj_units);
170  unfact = G_find_key_value("meters", proj_units);
171  if (unfact != NULL && (strcmp(pjinfo.proj, "ll") != 0))
172  G_asprintf(&proj4mod, "%s +to_meter=%s", proj4, unfact);
173  else
174  proj4mod = G_store(proj4);
175  pj_dalloc(proj4);
176 
177  if ((errcode = OSRImportFromProj4(hSRS, proj4mod)) != OGRERR_NONE) {
178  G_warning(_("OGR can't parse PROJ.4-style parameter string: "
179  "%s (OGR Error code was %d)"), proj4mod, errcode);
180  return NULL;
181  }
182  G_free(proj4mod);
183 
184  /* this messes up PROJCS versus GEOGCS!
185  sysname = G_find_key_value("name", proj_info);
186  if (sysname)
187  OSRSetProjCS(hSRS, sysname);
188  */
189 
190  if ((errcode = OSRExportToWkt(hSRS, &wkt)) != OGRERR_NONE) {
191  G_warning(_("OGR can't get WKT-style parameter string "
192  "(OGR Error code was %d)"), errcode);
193  return NULL;
194  }
195 
196  ellpskv = G_find_key_value("ellps", proj_info);
197  GPJ__get_ellipsoid_params(proj_info, &a, &es, &rf);
198  haveparams = GPJ__get_datum_params(proj_info, &datum, &params);
199 
200  if (ellpskv != NULL)
201  ellps = G_store(ellpskv);
202  else
203  ellps = NULL;
204 
205  if ((datum == NULL) || (GPJ_get_datum_by_name(datum, &dstruct) < 0)) {
206  datumlongname = G_store("unknown");
207  if (ellps == NULL)
208  ellps = G_store("unnamed");
209  }
210  else {
211  datumlongname = G_store(dstruct.longname);
212  if (ellps == NULL)
213  ellps = G_store(dstruct.ellps);
214  GPJ_free_datum(&dstruct);
215  }
216  G_debug(3, "GPJ_grass_to_osr: datum: <%s>", datum);
217  G_free(datum);
218  if (GPJ_get_ellipsoid_by_name(ellps, &estruct) > 0) {
219  ellpslong = G_store(estruct.longname);
220  DatumNameMassage(&ellpslong);
221  GPJ_free_ellps(&estruct);
222  }
223  else
224  ellpslong = G_store(ellps);
225 
226  startmod = strstr(wkt, "GEOGCS");
227  lastpart = strstr(wkt, "PRIMEM");
228  len = strlen(wkt) - strlen(startmod);
229  wkt[len] = '\0';
230  if (haveparams == 2) {
231  /* Only put datum params into the WKT if they were specifically
232  * specified in PROJ_INFO */
233  char *paramkey, *paramvalue;
234 
235  paramkey = strtok(params, "=");
236  paramvalue = params + strlen(paramkey) + 1;
237  if (G_strcasecmp(paramkey, "towgs84") == 0)
238  G_asprintf(&towgs84, ",TOWGS84[%s]", paramvalue);
239  else
240  towgs84 = G_store("");
241  G_free(params);
242  }
243  else
244  towgs84 = G_store("");
245 
246  sysname = OSRGetAttrValue(hSRS, "PROJCS", 0);
247  if (sysname == NULL) {
248  /* Not a projected co-ordinate system */
249  start = G_store("");
250  end = G_store("");
251  }
252  else {
253  if ((strcmp(sysname, "unnamed") == 0) &&
254  (G_find_key_value("name", proj_info) != NULL))
255  G_asprintf(&start, "PROJCS[\"%s\",",
256  G_find_key_value("name", proj_info));
257  else
258  start = G_store(wkt);
259 
260  osrunit = OSRGetAttrValue(hSRS, "UNIT", 0);
261  osrunfact = OSRGetAttrValue(hSRS, "UNIT", 1);
262 
263  if ((unfact == NULL) || (G_strcasecmp(osrunit, "unknown") != 0))
264  end = G_store("");
265  else {
266  char *buff;
267  double unfactf = atof(unfact);
268 
269  G_asprintf(&buff, ",UNIT[\"%s\",", osrunit);
270 
271  startmod = strstr(lastpart, buff);
272  len = strlen(lastpart) - strlen(startmod);
273  lastpart[len] = '\0';
274  G_free(buff);
275 
276  if (unit == NULL)
277  unit = "unknown";
278  G_asprintf(&end, ",UNIT[\"%s\",%.16g]]", unit, unfactf);
279  }
280 
281  }
282  OSRDestroySpatialReference(hSRS);
283  G_asprintf(&modwkt,
284  "%sGEOGCS[\"%s\",DATUM[\"%s\",SPHEROID[\"%s\",%.16g,%.16g]%s],%s%s",
285  start, ellps, datumlongname, ellpslong, a, rf, towgs84,
286  lastpart, end);
287  hSRS2 = OSRNewSpatialReference(modwkt);
288  G_free(modwkt);
289 
290  CPLFree(wkt);
291  G_free(start);
292  G_free(ellps);
293  G_free(datumlongname);
294  G_free(ellpslong);
295  G_free(towgs84);
296  G_free(end);
297 
298  return hSRS2;
299 }
300 
301 /*!
302  * \brief Converts a GRASS co-ordinate system to an
303  * OGRSpatialReferenceH object. EPSG code is preferred if available.
304  *
305  * The co-ordinate system definition is imported from EPSG (by GDAL)
306  * definition if available. TOWGS84 parameter is scanned from
307  * PROJ_INFO file and appended to co-ordinate system definition. If
308  * EPSG code is not available, PROJ_INFO file is used as
309  * GPJ_grass_to_osr() does.
310 
311  * \todo Merge with GPJ_grass_to_osr() in GRASS 8.
312  *
313  * \param proj_info Set of GRASS PROJ_INFO key/value pairs
314  * \param proj_units Set of GRASS PROJ_UNIT key/value pairs
315  * \param proj_epsg Set of GRASS PROJ_EPSG key/value pairs
316  *
317  * \return OGRSpatialReferenceH object representing the co-ordinate system
318  * defined by proj_info and proj_units or NULL if it fails
319  */
320 OGRSpatialReferenceH GPJ_grass_to_osr2(const struct Key_Value * proj_info,
321  const struct Key_Value * proj_units,
322  const struct Key_Value * proj_epsg)
323 {
324  int epsgcode = 0;
325 
326  if (proj_epsg) {
327  const char *epsgstr = G_find_key_value("epsg", proj_epsg);
328  if (epsgstr)
329  epsgcode = atoi(epsgstr);
330  }
331 
332  if (epsgcode) {
333  const char *towgs84;
334  OGRSpatialReferenceH hSRS;
335 
336  hSRS = OSRNewSpatialReference(NULL);
337 
338  OSRImportFromEPSG(hSRS, epsgcode);
339 
340  /* take +towgs84 from projinfo file if defined) */
341  towgs84 = G_find_key_value("towgs84", proj_info);
342  if (towgs84) {
343  char **tokens;
344  int i;
345  double df[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
346 
347  tokens = G_tokenize(towgs84, ",");
348 
349  for (i = 0; i < G_number_of_tokens(tokens); i++)
350  df[i] = atof(tokens[i]);
351  G_free_tokens(tokens);
352 
353  OSRSetTOWGS84(hSRS, df[0], df[1], df[2], df[3], df[4], df[5], df[6]);
354  }
355 
356  return hSRS;
357  }
358 
359  return GPJ_grass_to_osr(proj_info, proj_units);
360 }
361 
362 /*!
363  * \brief Converts an OGRSpatialReferenceH object to a GRASS co-ordinate system.
364  *
365  * \param cellhd Pointer to a GRASS Cell_head structure that will have its
366  * projection-related members populated with appropriate values
367  * \param projinfo Pointer to a pointer which will have a GRASS Key_Value
368  * structure allocated containing a set of GRASS PROJ_INFO values
369  * \param projunits Pointer to a pointer which will have a GRASS Key_Value
370  * structure allocated containing a set of GRASS PROJ_UNITS values
371  * \param hSRS OGRSpatialReferenceH object containing the co-ordinate
372  * system to be converted
373  * \param datumtrans Index number of datum parameter set to use, 0 to leave
374  * unspecified
375  *
376  * \return 2 if a projected or lat/long co-ordinate system has been
377  * defined
378  * \return 1 if an unreferenced XY co-ordinate system has
379  * been defined
380  */
381 int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
382  struct Key_Value **projunits, OGRSpatialReferenceH hSRS1,
383  int datumtrans)
384 {
385  struct Key_Value *temp_projinfo;
386  char *pszProj4 = NULL, *pszRemaining;
387  char *pszProj = NULL;
388  const char *pszProjCS = NULL;
389  char *datum = NULL;
390  struct gpj_datum dstruct;
391  const char *ograttr;
392  OGRSpatialReferenceH hSRS;
393 
394  *projinfo = NULL;
395  *projunits = NULL;
396 
397  hSRS = hSRS1;
398 
399  if (hSRS == NULL)
400  goto default_to_xy;
401 
402  /* Set finder function for locating OGR csv co-ordinate system tables */
403  /* SetCSVFilenameHook(GPJ_set_csv_loc); */
404 
405  /* Hopefully this doesn't do any harm if it wasn't in ESRI format
406  * to start with... */
407  OSRMorphFromESRI(hSRS);
408 
409  *projinfo = G_create_key_value();
410 
411  /* use proj4 definition from EXTENSION attribute if existing */
412  ograttr = OSRGetAttrValue(hSRS, "EXTENSION", 0);
413  if (ograttr && *ograttr && strcmp(ograttr, "PROJ4") == 0) {
414  ograttr = OSRGetAttrValue(hSRS, "EXTENSION", 1);
415  G_debug(3, "proj4 extension:");
416  G_debug(3, "%s", ograttr);
417 
418  if (ograttr && *ograttr) {
419  char *proj4ext;
420  OGRSpatialReferenceH hSRS2;
421 
422  hSRS2 = OSRNewSpatialReference(NULL);
423  proj4ext = G_store(ograttr);
424 
425  /* test */
426  if (OSRImportFromProj4(hSRS2, proj4ext) != OGRERR_NONE) {
427  G_warning(_("Updating spatial reference with embedded proj4 definition failed. "
428  "Proj4 definition: <%s>"), proj4ext);
429  OSRDestroySpatialReference(hSRS2);
430  }
431  else {
432  /* use new OGR spatial reference defined with embedded proj4 string */
433  /* TODO: replace warning with important_message once confirmed working */
434  G_warning(_("Updating spatial reference with embedded proj4 definition"));
435 
436  /* -------------------------------------------------------------------- */
437  /* Derive the user name for the coordinate system. */
438  /* -------------------------------------------------------------------- */
439  pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
440  if (!pszProjCS)
441  pszProjCS = OSRGetAttrValue(hSRS, "GEOGCS", 0);
442 
443  if (pszProjCS) {
444  G_set_key_value("name", pszProjCS, *projinfo);
445  }
446  else if (pszProj) {
447  char path[4095];
448  char name[80];
449 
450  /* use name of the projection as name for the coordinate system */
451 
452  sprintf(path, "%s/etc/proj/projections", G_gisbase());
453  if (G_lookup_key_value_from_file(path, pszProj, name, sizeof(name)) >
454  0)
455  G_set_key_value("name", name, *projinfo);
456  else
457  G_set_key_value("name", pszProj, *projinfo);
458  }
459 
460  /* the original hSRS1 is left as is, ok? */
461  hSRS = hSRS2;
462  }
463  G_free(proj4ext);
464  }
465  }
466 
467  /* -------------------------------------------------------------------- */
468  /* Set cellhd for well known coordinate systems. */
469  /* -------------------------------------------------------------------- */
470  if (!OSRIsGeographic(hSRS) && !OSRIsProjected(hSRS))
471  goto default_to_xy;
472 
473  if (cellhd) {
474  int bNorth;
475 
476  if (OSRIsGeographic(hSRS)) {
477  cellhd->proj = PROJECTION_LL;
478  cellhd->zone = 0;
479  }
480  else if (OSRGetUTMZone(hSRS, &bNorth) != 0) {
481  cellhd->proj = PROJECTION_UTM;
482  cellhd->zone = OSRGetUTMZone(hSRS, &bNorth);
483  if (!bNorth)
484  cellhd->zone *= -1;
485  }
486  else {
487  cellhd->proj = PROJECTION_OTHER;
488  cellhd->zone = 0;
489  }
490  }
491 
492  /* -------------------------------------------------------------------- */
493  /* Get the coordinate system definition in PROJ.4 format. */
494  /* -------------------------------------------------------------------- */
495  if (OSRExportToProj4(hSRS, &pszProj4) != OGRERR_NONE)
496  goto default_to_xy;
497 
498  /* -------------------------------------------------------------------- */
499  /* Parse the PROJ.4 string into key/value pairs. Do a bit of */
500  /* extra work to "GRASSify" the result. */
501  /* -------------------------------------------------------------------- */
502  temp_projinfo = G_create_key_value();
503 
504  /* Create "local" copy of proj4 string so we can modify and free it
505  * using GRASS functions */
506  pszRemaining = G_store(pszProj4);
507  CPLFree(pszProj4);
508  pszProj4 = pszRemaining;
509  while ((pszRemaining = strstr(pszRemaining, "+")) != NULL) {
510  char *pszToken, *pszValue;
511 
512  pszRemaining++;
513 
514  /* Advance pszRemaining to end of this token[=value] pair */
515  pszToken = pszRemaining;
516  while (*pszRemaining != ' ' && *pszRemaining != '\0')
517  pszRemaining++;
518 
519  if (*pszRemaining == ' ') {
520  *pszRemaining = '\0';
521  pszRemaining++;
522  }
523 
524  /* parse token, and value */
525  if (strstr(pszToken, "=") != NULL) {
526  pszValue = strstr(pszToken, "=");
527  *pszValue = '\0';
528  pszValue++;
529  }
530  else
531  pszValue = "defined";
532 
533  /* projection name */
534  if (G_strcasecmp(pszToken, "proj") == 0) {
535  /* The ll projection is known as longlat in PROJ.4 */
536  if (G_strcasecmp(pszValue, "longlat") == 0)
537  pszValue = "ll";
538 
539  pszProj = pszValue;
540  }
541 
542  /* discard @null nadgrids */
543  if (G_strcasecmp(pszToken, "nadgrids") == 0 &&
544  G_strcasecmp(pszValue, "@null") == 0)
545  continue;
546 
547  /* Ellipsoid and datum handled separately below */
548  if (G_strcasecmp(pszToken, "ellps") == 0
549  || G_strcasecmp(pszToken, "a") == 0
550  || G_strcasecmp(pszToken, "b") == 0
551  || G_strcasecmp(pszToken, "es") == 0
552  || G_strcasecmp(pszToken, "rf") == 0
553  || G_strcasecmp(pszToken, "datum") == 0)
554  continue;
555 
556  /* We will handle units separately */
557  if (G_strcasecmp(pszToken, "to_meter") == 0
558  || G_strcasecmp(pszToken, "units") == 0)
559  continue;
560 
561  G_set_key_value(pszToken, pszValue, temp_projinfo);
562  }
563  if (!pszProj)
564  G_warning(_("No projection name! Projection parameters likely to be meaningless."));
565 
566  /* -------------------------------------------------------------------- */
567  /* Derive the user name for the coordinate system. */
568  /* -------------------------------------------------------------------- */
569  if (!G_find_key_value("name", *projinfo)) {
570  pszProjCS = OSRGetAttrValue(hSRS, "PROJCS", 0);
571  if (!pszProjCS)
572  pszProjCS = OSRGetAttrValue(hSRS, "GEOGCS", 0);
573 
574  if (pszProjCS) {
575  G_set_key_value("name", pszProjCS, *projinfo);
576  }
577  else if (pszProj) {
578  char path[4095];
579  char name[80];
580 
581  /* use name of the projection as name for the coordinate system */
582 
583  sprintf(path, "%s/etc/proj/projections", G_gisbase());
584  if (G_lookup_key_value_from_file(path, pszProj, name, sizeof(name)) >
585  0)
586  G_set_key_value("name", name, *projinfo);
587  else
588  G_set_key_value("name", pszProj, *projinfo);
589  }
590  }
591 
592  /* -------------------------------------------------------------------- */
593  /* Find the GRASS datum name and choose parameters either */
594  /* interactively or not. */
595  /* -------------------------------------------------------------------- */
596 
597  {
598  const char *pszDatumNameConst = OSRGetAttrValue(hSRS, "DATUM", 0);
599  struct datum_list *list, *listhead;
600  char *dum1, *dum2, *pszDatumName;
601  int paramspresent =
602  GPJ__get_datum_params(temp_projinfo, &dum1, &dum2);
603 
604  if (pszDatumNameConst) {
605  /* Need to make a new copy of the string so we don't mess
606  * around with the memory inside the OGRSpatialReferenceH? */
607 
608  pszDatumName = G_store(pszDatumNameConst);
609  DatumNameMassage(&pszDatumName);
610  G_debug(3, "GPJ_osr_to_grass: pszDatumNameConst: <%s>", pszDatumName);
611 
612  list = listhead = read_datum_table();
613 
614  while (list != NULL) {
615  if (G_strcasecmp(pszDatumName, list->longname) == 0) {
616  datum = G_store(list->name);
617  break;
618  }
619  list = list->next;
620  }
621  free_datum_list(listhead);
622 
623  if (datum == NULL) {
624  if (paramspresent < 2)
625  /* Only give warning if no parameters present */
626  G_warning(_("Datum <%s> not recognised by GRASS and no parameters found"),
627  pszDatumName);
628  }
629  else {
630  G_set_key_value("datum", datum, *projinfo);
631 
632  if (paramspresent < 2) {
633  /* If no datum parameters were imported from the OSR
634  * object then we should use the set specified by datumtrans */
635  char *params, *chosenparams = NULL;
636  int paramsets;
637 
638  paramsets =
639  GPJ_get_default_datum_params_by_name(datum, &params);
640 
641  if (paramsets < 0)
642  G_warning(_("Datum <%s> apparently recognised by GRASS but no parameters found. "
643  "You may want to look into this."), datum);
644  else if (datumtrans > paramsets) {
645 
646  G_warning(_("Invalid transformation number %d; valid range is 1 to %d. "
647  "Leaving datum transform parameters unspecified."),
648  datumtrans, paramsets);
649  datumtrans = 0;
650  }
651 
652  if (paramsets > 0) {
653  struct gpj_datum_transform_list *tlist, *old;
654 
655  tlist = GPJ_get_datum_transform_by_name(datum);
656 
657  if (tlist != NULL) {
658  do {
659  if (tlist->count == datumtrans)
660  chosenparams = G_store(tlist->params);
661  old = tlist;
662  tlist = tlist->next;
664  } while (tlist != NULL);
665  }
666  }
667 
668  if (chosenparams != NULL) {
669  char *paramkey, *paramvalue;
670 
671  paramkey = strtok(chosenparams, "=");
672  paramvalue = chosenparams + strlen(paramkey) + 1;
673  G_set_key_value(paramkey, paramvalue, *projinfo);
674  G_free(chosenparams);
675  }
676 
677  if (paramsets > 0)
678  G_free(params);
679  }
680 
681  }
682  G_free(pszDatumName);
683  }
684  }
685 
686  /* -------------------------------------------------------------------- */
687  /* Determine an appropriate GRASS ellipsoid name if possible, or */
688  /* else just put a and es values into PROJ_INFO */
689  /* -------------------------------------------------------------------- */
690 
691  if ((datum != NULL) && (GPJ_get_datum_by_name(datum, &dstruct) > 0)) {
692  /* Use ellps name associated with datum */
693  G_set_key_value("ellps", dstruct.ellps, *projinfo);
694  GPJ_free_datum(&dstruct);
695  G_free(datum);
696  }
697  else {
698  /* If we can't determine the ellipsoid from the datum, derive it
699  * directly from "SPHEROID" parameters in WKT */
700  const char *pszSemiMajor = OSRGetAttrValue(hSRS, "SPHEROID", 1);
701  const char *pszInvFlat = OSRGetAttrValue(hSRS, "SPHEROID", 2);
702 
703  if (pszSemiMajor != NULL && pszInvFlat != NULL) {
704  char *ellps = NULL;
705  struct ellps_list *list, *listhead;
706  double a = atof(pszSemiMajor), invflat = atof(pszInvFlat), flat;
707  double es;
708 
709  /* Allow for incorrect WKT describing a sphere where InvFlat
710  * is given as 0 rather than inf */
711  if (invflat > 0)
712  flat = 1 / invflat;
713  else
714  flat = 0;
715 
716  es = flat * (2.0 - flat);
717 
718  list = listhead = read_ellipsoid_table(0);
719 
720  while (list != NULL) {
721  /* Try and match a and es against GRASS defined ellipsoids;
722  * accept first one that matches. These numbers were found
723  * by trial and error and could be fine-tuned, or possibly
724  * a direct comparison of IEEE floating point values used. */
725  if ((a == list->a || fabs(a - list->a) < 0.1 || fabs(1 - a / list->a) < 0.0000001) &&
726  ((es == 0 && list->es == 0) ||
727  /* Special case for sphere */
728  (invflat == list->rf || fabs(invflat - list->rf) < 0.0000001))) {
729  ellps = G_store(list->name);
730  break;
731  }
732  list = list->next;
733  }
734  if (listhead != NULL)
735  free_ellps_list(listhead);
736 
737  if (ellps == NULL) {
738  /* If we weren't able to find a matching ellps name, set
739  * a and es values directly from WKT-derived data */
740  char es_str[100];
741 
742  G_set_key_value("a", (char *)pszSemiMajor, *projinfo);
743 
744  sprintf(es_str, "%.16g", es);
745  G_set_key_value("es", es_str, *projinfo);
746  }
747  else {
748  /* else specify the GRASS ellps name for readability */
749  G_set_key_value("ellps", ellps, *projinfo);
750  G_free(ellps);
751  }
752 
753  }
754 
755  }
756 
757  /* -------------------------------------------------------------------- */
758  /* Finally append the detailed projection parameters to the end */
759  /* -------------------------------------------------------------------- */
760 
761  {
762  int i;
763 
764  for (i = 0; i < temp_projinfo->nitems; i++)
765  G_set_key_value(temp_projinfo->key[i],
766  temp_projinfo->value[i], *projinfo);
767 
768  G_free_key_value(temp_projinfo);
769  }
770 
771  G_free(pszProj4);
772 
773  /* -------------------------------------------------------------------- */
774  /* Set the linear units. */
775  /* -------------------------------------------------------------------- */
776  *projunits = G_create_key_value();
777 
778  if (OSRIsGeographic(hSRS)) {
779  /* We assume degrees ... someday we will be wrong! */
780  G_set_key_value("unit", "degree", *projunits);
781  G_set_key_value("units", "degrees", *projunits);
782  G_set_key_value("meters", "1.0", *projunits);
783  }
784  else {
785  char szFormatBuf[256];
786  char *pszUnitsName = NULL;
787  double dfToMeters;
788  char *pszUnitsPlural, *pszStringEnd;
789 
790  dfToMeters = OSRGetLinearUnits(hSRS, &pszUnitsName);
791 
792  /* Workaround for the most obvious case when unit name is unknown */
793  if ((G_strcasecmp(pszUnitsName, "unknown") == 0) &&
794  (dfToMeters == 1.))
795  G_asprintf(&pszUnitsName, "meter");
796 
797  if ((G_strcasecmp(pszUnitsName, "metre") == 0))
798  G_asprintf(&pszUnitsName, "meter");
799  if ((G_strcasecmp(pszUnitsName, "kilometre") == 0))
800  G_asprintf(&pszUnitsName, "kilometer");
801 
802  G_set_key_value("unit", pszUnitsName, *projunits);
803 
804  /* Attempt at plural formation (WKT format doesn't store plural
805  * form of unit name) */
806  pszUnitsPlural = G_malloc(strlen(pszUnitsName) + 3);
807  strcpy(pszUnitsPlural, pszUnitsName);
808  pszStringEnd = pszUnitsPlural + strlen(pszUnitsPlural) - 4;
809  if (G_strcasecmp(pszStringEnd, "foot") == 0) {
810  /* Special case for foot - change two o's to e's */
811  pszStringEnd[1] = 'e';
812  pszStringEnd[2] = 'e';
813  }
814  else if (G_strcasecmp(pszStringEnd, "inch") == 0) {
815  /* Special case for inch - add es */
816  pszStringEnd[4] = 'e';
817  pszStringEnd[5] = 's';
818  pszStringEnd[6] = '\0';
819  }
820  else {
821  /* For anything else add an s at the end */
822  pszStringEnd[4] = 's';
823  pszStringEnd[5] = '\0';
824  }
825 
826  G_set_key_value("units", pszUnitsPlural, *projunits);
827  G_free(pszUnitsPlural);
828 
829  sprintf(szFormatBuf, "%.16g", dfToMeters);
830  G_set_key_value("meters", szFormatBuf, *projunits);
831 
832  }
833 
834  if (hSRS != hSRS1)
835  OSRDestroySpatialReference(hSRS);
836 
837  return 2;
838 
839  /* -------------------------------------------------------------------- */
840  /* Fallback to returning an ungeoreferenced definition. */
841  /* -------------------------------------------------------------------- */
842  default_to_xy:
843  if (cellhd != NULL) {
844  cellhd->proj = PROJECTION_XY;
845  cellhd->zone = 0;
846  }
847  if (*projinfo)
848  G_free_key_value(*projinfo);
849 
850  *projinfo = NULL;
851  *projunits = NULL;
852 
853  if (hSRS != hSRS1)
854  OSRDestroySpatialReference(hSRS);
855 
856  return 1;
857 }
858 #endif
859 
860 /*!
861  * \brief Converts a WKT projection description to a GRASS co-ordinate system.
862  *
863  * \param cellhd Pointer to a GRASS Cell_head structure that will have its
864  * projection-related members populated with appropriate values
865  * \param projinfo Pointer to a pointer which will have a GRASS Key_Value
866  * structure allocated containing a set of GRASS PROJ_INFO values
867  * \param projunits Pointer to a pointer which will have a GRASS Key_Value
868  * structure allocated containing a set of GRASS PROJ_UNITS values
869  * \param wkt Well-known Text (WKT) description of the co-ordinate
870  * system to be converted
871  * \param datumtrans Index number of datum parameter set to use, 0 to leave
872  * unspecified
873  *
874  * \return 2 if a projected or lat/long co-ordinate system has been
875  * defined
876  * \return 1 if an unreferenced XY co-ordinate system has
877  * been defined
878  * \return -1 on error
879  */
880 int GPJ_wkt_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo,
881  struct Key_Value **projunits, const char *wkt,
882  int datumtrans)
883 {
884 #ifdef HAVE_OGR
885  int retval;
886 
887  if (wkt == NULL)
888  retval =
889  GPJ_osr_to_grass(cellhd, projinfo, projunits, NULL, datumtrans);
890  else {
891  OGRSpatialReferenceH hSRS;
892 
893  /* Set finder function for locating OGR csv co-ordinate system tables */
894  /* SetCSVFilenameHook(GPJ_set_csv_loc); */
895 
896  hSRS = OSRNewSpatialReference(wkt);
897  retval =
898  GPJ_osr_to_grass(cellhd, projinfo, projunits, hSRS, datumtrans);
899  OSRDestroySpatialReference(hSRS);
900  }
901 
902  return retval;
903 #else
904  return -1;
905 #endif
906 }
907 
908 #ifdef HAVE_OGR
909 /* GPJ_set_csv_loc()
910  * 'finder function' for use with OGR SetCSVFilenameHook() function */
911 
912 const char *GPJ_set_csv_loc(const char *name)
913 {
914  const char *gisbase = G_gisbase();
915  static char *buf = NULL;
916 
917  if (buf != NULL)
918  G_free(buf);
919 
920  G_asprintf(&buf, "%s%s/%s", gisbase, CSVDIR, name);
921 
922  return buf;
923 }
924 
925 
926 /* The list below is only for files that use a non-standard name for a
927  * datum that is already supported in GRASS. The number of entries must be even;
928  * they are all in pairs. The first one in the pair is the non-standard name;
929  * the second is the GRASS/GDAL name. If a name appears more than once (as for
930  * European_Terrestrial_Reference_System_1989) then it means there was more
931  * than one non-standard name for it that needs to be accounted for.
932  *
933  * N.B. The order of these pairs is different from that in
934  * ogr/ogrfromepsg.cpp in the GDAL source tree! GRASS uses the EPSG
935  * names in its WKT representation except WGS_1984 and WGS_1972 as
936  * these shortened versions seem to be standard.
937  * Below order:
938  * the equivalent name comes first in the pair, and
939  * the EPSG name (as used in the GRASS datum.table file) comes second.
940  *
941  * The datum parameters are stored in
942  * ../gis/datum.table # 3 parameters
943  * ../gis/datumtransform.table # 7 parameters (requires entry in datum.table)
944  *
945  * Hint: use GDAL's "testepsg" to identify the canonical name, e.g.
946  * testepsg epsg:4674
947  */
948 
949 static const char *papszDatumEquiv[] = {
950  "Militar_Geographische_Institute",
951  "Militar_Geographische_Institut",
952  "World_Geodetic_System_1984",
953  "WGS_1984",
954  "World_Geodetic_System_1972",
955  "WGS_1972",
956  "European_Terrestrial_Reference_System_89",
957  "European_Terrestrial_Reference_System_1989",
958  "European_Reference_System_1989",
959  "European_Terrestrial_Reference_System_1989",
960  "ETRS_1989",
961  "European_Terrestrial_Reference_System_1989",
962  "ETRS89",
963  "European_Terrestrial_Reference_System_1989",
964  "ETRF_1989",
965  "European_Terrestrial_Reference_System_1989",
966  "NZGD_2000",
967  "New_Zealand_Geodetic_Datum_2000",
968  "Monte_Mario_Rome",
969  "Monte_Mario",
970  "MONTROME",
971  "Monte_Mario",
972  "Campo_Inchauspe_1969",
973  "Campo_Inchauspe",
974  "S_JTSK",
975  "System_Jednotne_Trigonometricke_Site_Katastralni",
976  "S_JTSK_Ferro",
977  "Militar_Geographische_Institut",
978  "Potsdam_Datum_83",
979  "Deutsches_Hauptdreiecksnetz",
980  "South_American_1969",
981  "South_American_Datum_1969",
982  "ITRF_1992",
983  "ITRF92",
984  NULL
985 };
986 
987 /************************************************************************/
988 /* OGREPSGDatumNameMassage() */
989 /* */
990 /* Massage an EPSG datum name into WMT format. Also transform */
991 /* specific exception cases into WKT versions. */
992 
993 /************************************************************************/
994 
995 static void DatumNameMassage(char **ppszDatum)
996 {
997  int i, j;
998  char *pszDatum = *ppszDatum;
999 
1000  G_debug(3, "DatumNameMassage: Raw string found <%s>", (char *)pszDatum);
1001  /* -------------------------------------------------------------------- */
1002  /* Translate non-alphanumeric values to underscores. */
1003  /* -------------------------------------------------------------------- */
1004  for (i = 0; pszDatum[i] != '\0'; i++) {
1005  if (!(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z')
1006  && !(pszDatum[i] >= 'a' && pszDatum[i] <= 'z')
1007  && !(pszDatum[i] >= '0' && pszDatum[i] <= '9')) {
1008  pszDatum[i] = '_';
1009  }
1010  }
1011 
1012  /* -------------------------------------------------------------------- */
1013  /* Remove repeated and trailing underscores. */
1014  /* -------------------------------------------------------------------- */
1015  for (i = 1, j = 0; pszDatum[i] != '\0'; i++) {
1016  if (pszDatum[j] == '_' && pszDatum[i] == '_')
1017  continue;
1018 
1019  pszDatum[++j] = pszDatum[i];
1020  }
1021  if (pszDatum[j] == '_')
1022  pszDatum[j] = '\0';
1023  else
1024  pszDatum[j + 1] = '\0';
1025 
1026  /* -------------------------------------------------------------------- */
1027  /* Search for datum equivalences. Specific massaged names get */
1028  /* mapped to OpenGIS specified names. */
1029  /* -------------------------------------------------------------------- */
1030  G_debug(3, "DatumNameMassage: Search for datum equivalences of <%s>", (char *)pszDatum);
1031  for (i = 0; papszDatumEquiv[i] != NULL; i += 2) {
1032  if (EQUAL(*ppszDatum, papszDatumEquiv[i])) {
1033  G_free(*ppszDatum);
1034  *ppszDatum = G_store(papszDatumEquiv[i + 1]);
1035  break;
1036  }
1037  }
1038 }
1039 
1040 #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:880
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 G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
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
2D/3D raster map header (used also for region)
Definition: gis.h:390
struct ellps_list * read_ellipsoid_table(int fatal)
Definition: ellipse.c:224
char proj[100]
Definition: gprojects.h:38
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:87
int old
Definition: raster3d/cats.c:84
char ** G_tokenize(const char *buf, const char *delim)
Tokenize string.
Definition: gis/token.c:48
char * G_store(const char *s)
Copy string to allocated memory.
Definition: strings.c:86
double a
Definition: gprojects.h:66
int G_asprintf(char **out, const char *fmt,...)
Definition: asprintf.c:70
#define PROJECTION_XY
Projection code - XY coordinate system (unreferenced data)
Definition: gis.h:93
char * GPJ_grass_to_wkt2(const struct Key_Value *proj_info, const struct Key_Value *proj_units, const struct Key_Value *proj_epsg, int esri_style, int prettify)
Converts a GRASS co-ordinate system representation to WKT style. EPSG code is preferred if available...
Definition: convert.c:118
#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 PROJECTION_UTM
Projection code - UTM.
Definition: gis.h:95
int G_number_of_tokens(char **tokens)
Return number of tokens.
Definition: gis/token.c:185
projPJ pj
Definition: gprojects.h:35
#define PROJECTION_OTHER
Projection code - other projection (other then noted above)
Definition: gis.h:101
void free_datum_list(struct datum_list *dstruct)
Free the memory used by a datum_list linked list structure.
Definition: proj/datum.c:417
char ** value
Definition: gis.h:484
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:912
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
char * ellps
Definition: gprojects.h:43
int zone
Projection zone (UTM)
Definition: gis.h:424
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:136
OGRSpatialReferenceH GPJ_grass_to_osr2(const struct Key_Value *proj_info, const struct Key_Value *proj_units, const struct Key_Value *proj_epsg)
Converts a GRASS co-ordinate system to an OGRSpatialReferenceH object. EPSG code is preferred if avai...
Definition: convert.c:320
#define PROJECTION_LL
Projection code - Latitude-Longitude.
Definition: gis.h:99
char * longname
Definition: gprojects.h:43
char ** key
Definition: gis.h:483
int proj
Projection code.
Definition: gis.h:422
struct datum_list * read_datum_table(void)
Read the current GRASS datum.table from disk and store in memory.
Definition: proj/datum.c:345
Definition: gis.h:479
char * longname
Definition: gprojects.h:65
double es
Definition: gprojects.h:66
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)
&quot;Last resort&quot; function to retrieve a &quot;default&quot; set of datum parameters for a datum (N...
Definition: proj/datum.c:85
Definition: path.h:16
struct gpj_datum_transform_list * next
Definition: gprojects.h:59
void G_free_tokens(char **tokens)
Free memory allocated to tokens.
Definition: gis/token.c:204
#define _(str)
Definition: glocale.h:13
int nitems
Definition: gis.h:481
char buff[1024]
Definition: raster3d/cats.c:81
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
int GPJ_osr_to_grass(struct Cell_head *cellhd, struct Key_Value **projinfo, struct Key_Value **projunits, OGRSpatialReferenceH hSRS1, int datumtrans)
Converts an OGRSpatialReferenceH object to a GRASS co-ordinate system.
Definition: convert.c:381
const char * G_gisbase(void)
Get full path name of the top level module directory.
Definition: gisbase.c:41
double rf
Definition: gprojects.h:66
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
#define EQUAL
Definition: sqlp.tab.c:159
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