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