1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * PostGIS is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright (C) 2001-2003 Refractions Research Inc.
22 *
23 **********************************************************************/
24
25
26 #include "../postgis_config.h"
27 #include "liblwgeom_internal.h"
28 #include "lwgeom_log.h"
29 #include <string.h>
30
31 /** convert decimal degress to radians */
32 static void
to_rad(POINT4D * pt)33 to_rad(POINT4D *pt)
34 {
35 pt->x *= M_PI/180.0;
36 pt->y *= M_PI/180.0;
37 }
38
39 /** convert radians to decimal degress */
40 static void
to_dec(POINT4D * pt)41 to_dec(POINT4D *pt)
42 {
43 pt->x *= 180.0/M_PI;
44 pt->y *= 180.0/M_PI;
45 }
46
47 /***************************************************************************/
48
49 #if POSTGIS_PROJ_VERSION < 60
50
51 static int
point4d_transform(POINT4D * pt,LWPROJ * pj)52 point4d_transform(POINT4D *pt, LWPROJ *pj)
53 {
54 POINT3D orig_pt = {pt->x, pt->y, pt->z}; /* Copy for error report*/
55
56 if (pj_is_latlong(pj->pj_from)) to_rad(pt) ;
57
58 LWDEBUGF(4, "transforming POINT(%f %f) from '%s' to '%s'",
59 orig_pt.x, orig_pt.y, pj_get_def(pj->pj_from,0), pj_get_def(pj->pj_to,0));
60
61 if (pj_transform(pj->pj_from, pj->pj_to, 1, 0, &(pt->x), &(pt->y), &(pt->z)) != 0)
62 {
63 int pj_errno_val = *pj_get_errno_ref();
64 if (pj_errno_val == -38)
65 {
66 lwnotice("PostGIS was unable to transform the point because either no grid"
67 " shift files were found, or the point does not lie within the "
68 "range for which the grid shift is defined. Refer to the "
69 "ST_Transform() section of the PostGIS manual for details on how "
70 "to configure PostGIS to alter this behaviour.");
71 lwerror("transform: couldn't project point (%g %g %g): %s (%d)",
72 orig_pt.x, orig_pt.y, orig_pt.z,
73 pj_strerrno(pj_errno_val), pj_errno_val);
74 }
75 else
76 {
77 lwerror("transform: %s (%d)",
78 pj_strerrno(pj_errno_val), pj_errno_val);
79 }
80 return LW_FAILURE;
81 }
82
83 if (pj_is_latlong(pj->pj_to)) to_dec(pt);
84 return LW_SUCCESS;
85 }
86
87 /**
88 * Transform given POINTARRAY
89 * from inpj projection to outpj projection
90 */
91 int
ptarray_transform(POINTARRAY * pa,LWPROJ * pj)92 ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
93 {
94 uint32_t i;
95 POINT4D p;
96
97 for ( i = 0; i < pa->npoints; i++ )
98 {
99 getPoint4d_p(pa, i, &p);
100 if ( ! point4d_transform(&p, pj) ) return LW_FAILURE;
101 ptarray_set_point4d(pa, i, &p);
102 }
103
104 return LW_SUCCESS;
105 }
106
107 int
lwgeom_transform_from_str(LWGEOM * geom,const char * instr,const char * outstr)108 lwgeom_transform_from_str(LWGEOM *geom, const char* instr, const char* outstr)
109 {
110 char *pj_errstr;
111 int rv;
112 LWPROJ pj;
113
114 pj.pj_from = projpj_from_string(instr);
115 if (!pj.pj_from)
116 {
117 pj_errstr = pj_strerrno(*pj_get_errno_ref());
118 if (!pj_errstr) pj_errstr = "";
119 lwerror("could not parse proj string '%s'", instr);
120 return LW_FAILURE;
121 }
122
123 pj.pj_to = projpj_from_string(outstr);
124 if (!pj.pj_to)
125 {
126 pj_free(pj.pj_from);
127 pj_errstr = pj_strerrno(*pj_get_errno_ref());
128 if (!pj_errstr) pj_errstr = "";
129 lwerror("could not parse proj string '%s'", outstr);
130 return LW_FAILURE;
131 }
132
133 rv = lwgeom_transform(geom, &pj);
134 pj_free(pj.pj_from);
135 pj_free(pj.pj_to);
136 return rv;
137 }
138
139 /**
140 * Transform given LWGEOM geometry
141 * from inpj projection to outpj projection
142 */
143 int
lwgeom_transform(LWGEOM * geom,LWPROJ * pj)144 lwgeom_transform(LWGEOM *geom, LWPROJ *pj)
145 {
146 uint32_t i;
147
148 /* No points to transform in an empty! */
149 if ( lwgeom_is_empty(geom) )
150 return LW_SUCCESS;
151
152 switch(geom->type)
153 {
154 case POINTTYPE:
155 case LINETYPE:
156 case CIRCSTRINGTYPE:
157 case TRIANGLETYPE:
158 {
159 LWLINE *g = (LWLINE*)geom;
160 if ( ! ptarray_transform(g->points, pj) ) return LW_FAILURE;
161 break;
162 }
163 case POLYGONTYPE:
164 {
165 LWPOLY *g = (LWPOLY*)geom;
166 for ( i = 0; i < g->nrings; i++ )
167 {
168 if ( ! ptarray_transform(g->rings[i], pj) ) return LW_FAILURE;
169 }
170 break;
171 }
172 case MULTIPOINTTYPE:
173 case MULTILINETYPE:
174 case MULTIPOLYGONTYPE:
175 case COLLECTIONTYPE:
176 case COMPOUNDTYPE:
177 case CURVEPOLYTYPE:
178 case MULTICURVETYPE:
179 case MULTISURFACETYPE:
180 case POLYHEDRALSURFACETYPE:
181 case TINTYPE:
182 {
183 LWCOLLECTION *g = (LWCOLLECTION*)geom;
184 for ( i = 0; i < g->ngeoms; i++ )
185 {
186 if ( ! lwgeom_transform(g->geoms[i], pj) ) return LW_FAILURE;
187 }
188 break;
189 }
190 default:
191 {
192 lwerror("lwgeom_transform: Cannot handle type '%s'",
193 lwtype_name(geom->type));
194 return LW_FAILURE;
195 }
196 }
197 return LW_SUCCESS;
198 }
199
200 projPJ
projpj_from_string(const char * str1)201 projpj_from_string(const char *str1)
202 {
203 if (!str1 || str1[0] == '\0')
204 {
205 return NULL;
206 }
207 return pj_init_plus(str1);
208 }
209
210 /***************************************************************************/
211
212 #else /* POSTGIS_PROJ_VERION >= 60 */
213
214 static PJ *
proj_cs_get_simplecs(const PJ * pj_crs)215 proj_cs_get_simplecs(const PJ *pj_crs)
216 {
217 PJ *pj_sub = NULL;
218 if (proj_get_type(pj_crs) == PJ_TYPE_COMPOUND_CRS)
219 {
220 /* Sub-CRS[0] is the horizontal component */
221 pj_sub = proj_crs_get_sub_crs(NULL, pj_crs, 0);
222 if (!pj_sub)
223 lwerror("%s: proj_crs_get_sub_crs(0) returned NULL", __func__);
224 }
225 else if (proj_get_type(pj_crs) == PJ_TYPE_BOUND_CRS)
226 {
227 pj_sub = proj_get_source_crs(NULL, pj_crs);
228 if (!pj_sub)
229 lwerror("%s: proj_get_source_crs returned NULL", __func__);
230 }
231 else
232 {
233 /* If this works, we have a CS so we can return */
234 pj_sub = proj_crs_get_coordinate_system(NULL, pj_crs);
235 if (pj_sub)
236 return pj_sub;
237 }
238
239 /* Only sub-components of the Compound or Bound CRS's get here */
240 /* If we failed to get sub-components, or we failed to extract */
241 /* a CS from a generic CRS, then this is another case we don't */
242 /* handle */
243 if (!pj_sub)
244 lwerror("%s: %s", __func__, proj_errno_string(proj_context_errno(NULL)));
245
246 /* If the components are usable, we can extract the CS and return */
247 int pj_type = proj_get_type(pj_sub);
248 if (pj_type == PJ_TYPE_GEOGRAPHIC_2D_CRS || pj_type == PJ_TYPE_PROJECTED_CRS)
249 {
250 PJ *pj_2d = proj_crs_get_coordinate_system(NULL, pj_sub);
251 proj_destroy(pj_sub);
252 return pj_2d;
253 }
254
255 /* If the components are *themselves* Bound/Compound, we can recurse */
256 if (pj_type == PJ_TYPE_COMPOUND_CRS || pj_type == PJ_TYPE_BOUND_CRS)
257 return proj_cs_get_simplecs(pj_sub);
258
259 /* This is a case we don't know how to handle */
260 lwerror("%s: un-handled CRS sub-type: %s", __func__, pj_type);
261 return NULL;
262 }
263
264
265 #define STR_EQUALS(A, B) strcmp((A), (B)) == 0
266 #define STR_IEQUALS(A, B) (strcasecmp((A), (B)) == 0)
267 #define STR_ISTARTS(A, B) (strncasecmp((A), (B), strlen((B))) == 0)
268
269 static uint8_t
proj_crs_is_swapped(const PJ * pj_crs)270 proj_crs_is_swapped(const PJ *pj_crs)
271 {
272 int axis_count;
273 PJ *pj_cs = proj_cs_get_simplecs(pj_crs);
274 if (!pj_cs)
275 lwerror("%s: proj_cs_get_simplecs returned NULL", __func__);
276
277 axis_count = proj_cs_get_axis_count(NULL, pj_cs);
278 if (axis_count >= 2)
279 {
280 const char *out_name1, *out_abbrev1, *out_direction1;
281 const char *out_name2, *out_abbrev2, *out_direction2;
282 /* Read first axis */
283 proj_cs_get_axis_info(NULL,
284 pj_cs, 0,
285 &out_name1, &out_abbrev1, &out_direction1,
286 NULL, NULL, NULL, NULL);
287 /* Read second axis */
288 proj_cs_get_axis_info(NULL,
289 pj_cs, 1,
290 &out_name2, &out_abbrev2, &out_direction2,
291 NULL, NULL, NULL, NULL);
292
293 proj_destroy(pj_cs);
294
295 /* Directions agree, this is a northing/easting CRS, so reverse it */
296 if(out_direction1 && STR_IEQUALS(out_direction1, "north") &&
297 out_direction2 && STR_IEQUALS(out_direction2, "east") )
298 {
299 return LW_TRUE;
300 }
301
302 /* Oddball case? Both axes north / both axes south, swap */
303 if(out_direction1 && out_direction2 &&
304 ((STR_IEQUALS(out_direction1, "north") && STR_IEQUALS(out_direction2, "north")) ||
305 (STR_IEQUALS(out_direction1, "south") && STR_IEQUALS(out_direction2, "south"))) &&
306 out_name1 && STR_ISTARTS(out_name1, "northing") &&
307 out_name2 && STR_ISTARTS(out_name2, "easting"))
308 {
309 return LW_TRUE;
310 }
311
312 /* Any lat/lon system with Lat in first axis gets swapped */
313 if (STR_ISTARTS(out_abbrev1, "Lat"))
314 return LW_TRUE;
315
316 return LW_FALSE;
317 }
318
319 /* Failed the axis count test, leave quietly */
320 proj_destroy(pj_cs);
321 return LW_FALSE;
322 }
323
324 LWPROJ *
lwproj_from_PJ(PJ * pj,int8_t extra_geography_data)325 lwproj_from_PJ(PJ *pj, int8_t extra_geography_data)
326 {
327 PJ *pj_source_crs = proj_get_source_crs(NULL, pj);
328 uint8_t source_is_latlong = LW_FALSE;
329 double out_semi_major_metre = DBL_MAX, out_semi_minor_metre = DBL_MAX;
330
331 if (!pj_source_crs)
332 {
333 lwerror("%s: unable to access source crs", __func__);
334 return NULL;
335 }
336 uint8_t source_swapped = proj_crs_is_swapped(pj_source_crs);
337
338 /* We only care about the extra values if there is no transformation */
339 if (!extra_geography_data)
340 {
341 proj_destroy(pj_source_crs);
342 }
343 else
344 {
345 PJ *pj_ellps;
346 double out_inv_flattening;
347 int out_is_semi_minor_computed;
348
349 PJ_TYPE pj_type = proj_get_type(pj_source_crs);
350 if (pj_type == PJ_TYPE_UNKNOWN)
351 {
352 proj_destroy(pj_source_crs);
353 lwerror("%s: unable to access source crs type", __func__);
354 return NULL;
355 }
356 source_is_latlong = (pj_type == PJ_TYPE_GEOGRAPHIC_2D_CRS) || (pj_type == PJ_TYPE_GEOGRAPHIC_3D_CRS);
357
358 pj_ellps = proj_get_ellipsoid(NULL, pj_source_crs);
359 proj_destroy(pj_source_crs);
360 if (!pj_ellps)
361 {
362 lwerror("%s: unable to access source crs ellipsoid", __func__);
363 return NULL;
364 }
365 if (!proj_ellipsoid_get_parameters(NULL,
366 pj_ellps,
367 &out_semi_major_metre,
368 &out_semi_minor_metre,
369 &out_is_semi_minor_computed,
370 &out_inv_flattening))
371 {
372 proj_destroy(pj_ellps);
373 lwerror("%s: unable to access source crs ellipsoid parameters", __func__);
374 return NULL;
375 }
376 proj_destroy(pj_ellps);
377 }
378
379 PJ *pj_target_crs = proj_get_target_crs(NULL, pj);
380 if (!pj_target_crs)
381 {
382 lwerror("%s: unable to access target crs", __func__);
383 return NULL;
384 }
385 uint8_t target_swapped = proj_crs_is_swapped(pj_target_crs);
386 proj_destroy(pj_target_crs);
387
388 LWPROJ *lp = lwalloc(sizeof(LWPROJ));
389 lp->pj = pj;
390 lp->source_swapped = source_swapped;
391 lp->target_swapped = target_swapped;
392 lp->source_is_latlong = source_is_latlong;
393 lp->source_semi_major_metre = out_semi_major_metre;
394 lp->source_semi_minor_metre = out_semi_minor_metre;
395
396 return lp;
397 }
398
399 int
lwgeom_transform_from_str(LWGEOM * geom,const char * instr,const char * outstr)400 lwgeom_transform_from_str(LWGEOM *geom, const char* instr, const char* outstr)
401 {
402 PJ *pj = proj_create_crs_to_crs(NULL, instr, outstr, NULL);
403 if (!pj)
404 {
405 PJ *pj_in = proj_create(NULL, instr);
406 if (!pj_in)
407 {
408 lwerror("could not parse proj string '%s'", instr);
409 }
410 proj_destroy(pj_in);
411
412 PJ *pj_out = proj_create(NULL, outstr);
413 if (!pj_out)
414 {
415 lwerror("could not parse proj string '%s'", outstr);
416 }
417 proj_destroy(pj_out);
418 lwerror("%s: Failed to transform", __func__);
419 return LW_FAILURE;
420 }
421
422 LWPROJ *lp = lwproj_from_PJ(pj, LW_FALSE);
423
424 int ret = lwgeom_transform(geom, lp);
425
426 proj_destroy(pj);
427 lwfree(lp);
428
429 return ret;
430 }
431
432 int
lwgeom_transform(LWGEOM * geom,LWPROJ * pj)433 lwgeom_transform(LWGEOM *geom, LWPROJ *pj)
434 {
435 uint32_t i;
436
437 /* No points to transform in an empty! */
438 if (lwgeom_is_empty(geom))
439 return LW_SUCCESS;
440
441 switch(geom->type)
442 {
443 case POINTTYPE:
444 case LINETYPE:
445 case CIRCSTRINGTYPE:
446 case TRIANGLETYPE:
447 {
448 LWLINE *g = (LWLINE*)geom;
449 if (!ptarray_transform(g->points, pj))
450 return LW_FAILURE;
451 break;
452 }
453 case POLYGONTYPE:
454 {
455 LWPOLY *g = (LWPOLY*)geom;
456 for (i = 0; i < g->nrings; i++)
457 {
458 if (!ptarray_transform(g->rings[i], pj))
459 return LW_FAILURE;
460 }
461 break;
462 }
463 case MULTIPOINTTYPE:
464 case MULTILINETYPE:
465 case MULTIPOLYGONTYPE:
466 case COLLECTIONTYPE:
467 case COMPOUNDTYPE:
468 case CURVEPOLYTYPE:
469 case MULTICURVETYPE:
470 case MULTISURFACETYPE:
471 case POLYHEDRALSURFACETYPE:
472 case TINTYPE:
473 {
474 LWCOLLECTION *g = (LWCOLLECTION*)geom;
475 for (i = 0; i < g->ngeoms; i++)
476 {
477 if (!lwgeom_transform(g->geoms[i], pj))
478 return LW_FAILURE;
479 }
480 break;
481 }
482 default:
483 {
484 lwerror("lwgeom_transform: Cannot handle type '%s'",
485 lwtype_name(geom->type));
486 return LW_FAILURE;
487 }
488 }
489 return LW_SUCCESS;
490 }
491
492 int
ptarray_transform(POINTARRAY * pa,LWPROJ * pj)493 ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
494 {
495 uint32_t i;
496 POINT4D p;
497 size_t n_converted;
498 size_t n_points = pa->npoints;
499 size_t point_size = ptarray_point_size(pa);
500 int has_z = ptarray_has_z(pa);
501 double *pa_double = (double*)(pa->serialized_pointlist);
502
503 /* Convert to radians if necessary */
504 if (proj_angular_input(pj->pj, PJ_FWD))
505 {
506 for (i = 0; i < pa->npoints; i++)
507 {
508 getPoint4d_p(pa, i, &p);
509 to_rad(&p);
510 }
511 }
512
513 if (pj->source_swapped)
514 ptarray_swap_ordinates(pa, LWORD_X, LWORD_Y);
515
516 if (n_points == 1)
517 {
518 /* For single points it's faster to call proj_trans */
519 PJ_XYZT v = {pa_double[0], pa_double[1], has_z ? pa_double[2] : 0.0, 0.0};
520 PJ_COORD t = proj_trans(pj->pj, PJ_FWD, (PJ_COORD)v);
521
522 int pj_errno_val = proj_errno(pj->pj);
523 if (pj_errno_val)
524 {
525 lwerror("transform: %s (%d)", proj_errno_string(pj_errno_val), pj_errno_val);
526 return LW_FAILURE;
527 }
528 pa_double[0] = ((PJ_XYZT)t.xyzt).x;
529 pa_double[1] = ((PJ_XYZT)t.xyzt).y;
530 if (has_z)
531 pa_double[2] = ((PJ_XYZT)t.xyzt).z;
532 }
533 else
534 {
535 /*
536 * size_t proj_trans_generic(PJ *P, PJ_DIRECTION direction,
537 * double *x, size_t sx, size_t nx,
538 * double *y, size_t sy, size_t ny,
539 * double *z, size_t sz, size_t nz,
540 * double *t, size_t st, size_t nt)
541 */
542
543 n_converted = proj_trans_generic(pj->pj,
544 PJ_FWD,
545 pa_double,
546 point_size,
547 n_points, /* X */
548 pa_double + 1,
549 point_size,
550 n_points, /* Y */
551 has_z ? pa_double + 2 : NULL,
552 has_z ? point_size : 0,
553 has_z ? n_points : 0, /* Z */
554 NULL,
555 0,
556 0 /* M */
557 );
558
559 if (n_converted != n_points)
560 {
561 lwerror("ptarray_transform: converted (%d) != input (%d)", n_converted, n_points);
562 return LW_FAILURE;
563 }
564
565 int pj_errno_val = proj_errno(pj->pj);
566 if (pj_errno_val)
567 {
568 lwerror("transform: %s (%d)", proj_errno_string(pj_errno_val), pj_errno_val);
569 return LW_FAILURE;
570 }
571 }
572
573 if (pj->target_swapped)
574 ptarray_swap_ordinates(pa, LWORD_X, LWORD_Y);
575
576 /* Convert radians to degrees if necessary */
577 if (proj_angular_output(pj->pj, PJ_FWD))
578 {
579 for (i = 0; i < pa->npoints; i++)
580 {
581 getPoint4d_p(pa, i, &p);
582 to_dec(&p);
583 }
584 }
585
586 return LW_SUCCESS;
587 }
588
589 #endif
590