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