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 < 61
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
141 projPJ
projpj_from_string(const char * str1)142 projpj_from_string(const char *str1)
143 {
144 if (!str1 || str1[0] == '\0')
145 {
146 return NULL;
147 }
148 return pj_init_plus(str1);
149 }
150
151 /***************************************************************************/
152
153 #else /* POSTGIS_PROJ_VERION >= 61 */
154
155 LWPROJ *
lwproj_from_str(const char * str_in,const char * str_out)156 lwproj_from_str(const char* str_in, const char* str_out)
157 {
158 uint8_t source_is_latlong = LW_FALSE;
159 double semi_major_metre = DBL_MAX, semi_minor_metre = DBL_MAX;
160
161 /* Usable inputs? */
162 if (! (str_in && str_out))
163 return NULL;
164
165 PJ* pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX, str_in, str_out, NULL);
166 if (!pj)
167 return NULL;
168
169 /* Fill in geodetic parameter information when a null-transform */
170 /* is passed, because that's how we signal we want to store */
171 /* that info in the cache */
172 if (strcmp(str_in, str_out) == 0)
173 {
174 PJ *pj_source_crs = proj_get_source_crs(PJ_DEFAULT_CTX, pj);
175 PJ *pj_ellps;
176 PJ_TYPE pj_type = proj_get_type(pj_source_crs);
177 if (pj_type == PJ_TYPE_UNKNOWN)
178 {
179 proj_destroy(pj);
180 lwerror("%s: unable to access source crs type", __func__);
181 return NULL;
182 }
183 source_is_latlong = (pj_type == PJ_TYPE_GEOGRAPHIC_2D_CRS) || (pj_type == PJ_TYPE_GEOGRAPHIC_3D_CRS);
184
185 pj_ellps = proj_get_ellipsoid(PJ_DEFAULT_CTX, pj_source_crs);
186 proj_destroy(pj_source_crs);
187 if (!pj_ellps)
188 {
189 proj_destroy(pj);
190 lwerror("%s: unable to access source crs ellipsoid", __func__);
191 return NULL;
192 }
193 if (!proj_ellipsoid_get_parameters(PJ_DEFAULT_CTX,
194 pj_ellps,
195 &semi_major_metre,
196 &semi_minor_metre,
197 NULL,
198 NULL))
199 {
200 proj_destroy(pj_ellps);
201 proj_destroy(pj);
202 lwerror("%s: unable to access source crs ellipsoid parameters", __func__);
203 return NULL;
204 }
205 proj_destroy(pj_ellps);
206 }
207
208 /* Add in an axis swap if necessary */
209 PJ* pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, pj);
210 /* Swap failed for some reason? Fall back to coordinate operation */
211 if (!pj_norm)
212 pj_norm = pj;
213 /* Swap is not a copy of input? Clean up input */
214 else if (pj != pj_norm)
215 proj_destroy(pj);
216
217 /* Allocate and populate return value */
218 LWPROJ *lp = lwalloc(sizeof(LWPROJ));
219 lp->pj = pj_norm; /* Caller is going to have to explicitly proj_destroy this */
220 lp->source_is_latlong = source_is_latlong;
221 lp->source_semi_major_metre = semi_major_metre;
222 lp->source_semi_minor_metre = semi_minor_metre;
223 return lp;
224 }
225
226 int
lwgeom_transform_from_str(LWGEOM * geom,const char * instr,const char * outstr)227 lwgeom_transform_from_str(LWGEOM *geom, const char* instr, const char* outstr)
228 {
229 LWPROJ *lp = lwproj_from_str(instr, outstr);
230 if (!lp)
231 {
232 PJ *pj_in = proj_create(PJ_DEFAULT_CTX, instr);
233 if (!pj_in)
234 {
235 lwerror("could not parse proj string '%s'", instr);
236 }
237 proj_destroy(pj_in);
238
239 PJ *pj_out = proj_create(PJ_DEFAULT_CTX, outstr);
240 if (!pj_out)
241 {
242 lwerror("could not parse proj string '%s'", outstr);
243 }
244 proj_destroy(pj_out);
245 lwerror("%s: Failed to transform", __func__);
246 return LW_FAILURE;
247 }
248 int ret = lwgeom_transform(geom, lp);
249 proj_destroy(lp->pj);
250 lwfree(lp);
251 return ret;
252 }
253
254 int
ptarray_transform(POINTARRAY * pa,LWPROJ * pj)255 ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
256 {
257 uint32_t i;
258 POINT4D p;
259 size_t n_converted;
260 size_t n_points = pa->npoints;
261 size_t point_size = ptarray_point_size(pa);
262 int has_z = ptarray_has_z(pa);
263 double *pa_double = (double*)(pa->serialized_pointlist);
264
265 /* Convert to radians if necessary */
266 if (proj_angular_input(pj->pj, PJ_FWD))
267 {
268 for (i = 0; i < pa->npoints; i++)
269 {
270 getPoint4d_p(pa, i, &p);
271 to_rad(&p);
272 }
273 }
274
275 if (n_points == 1)
276 {
277 /* For single points it's faster to call proj_trans */
278 PJ_XYZT v = {pa_double[0], pa_double[1], has_z ? pa_double[2] : 0.0, 0.0};
279 PJ_COORD c;
280 c.xyzt = v;
281 PJ_COORD t = proj_trans(pj->pj, PJ_FWD, c);
282
283 int pj_errno_val = proj_errno_reset(pj->pj);
284 if (pj_errno_val)
285 {
286 lwerror("transform: %s (%d)", proj_errno_string(pj_errno_val), pj_errno_val);
287 return LW_FAILURE;
288 }
289 pa_double[0] = (t.xyzt).x;
290 pa_double[1] = (t.xyzt).y;
291 if (has_z)
292 pa_double[2] = (t.xyzt).z;
293 }
294 else
295 {
296 /*
297 * size_t proj_trans_generic(PJ *P, PJ_DIRECTION direction,
298 * double *x, size_t sx, size_t nx,
299 * double *y, size_t sy, size_t ny,
300 * double *z, size_t sz, size_t nz,
301 * double *t, size_t st, size_t nt)
302 */
303
304 n_converted = proj_trans_generic(pj->pj,
305 PJ_FWD,
306 pa_double,
307 point_size,
308 n_points, /* X */
309 pa_double + 1,
310 point_size,
311 n_points, /* Y */
312 has_z ? pa_double + 2 : NULL,
313 has_z ? point_size : 0,
314 has_z ? n_points : 0, /* Z */
315 NULL,
316 0,
317 0 /* M */
318 );
319
320 if (n_converted != n_points)
321 {
322 lwerror("ptarray_transform: converted (%d) != input (%d)", n_converted, n_points);
323 return LW_FAILURE;
324 }
325
326 int pj_errno_val = proj_errno_reset(pj->pj);
327 if (pj_errno_val)
328 {
329 lwerror("transform: %s (%d)", proj_errno_string(pj_errno_val), pj_errno_val);
330 return LW_FAILURE;
331 }
332 }
333
334 /* Convert radians to degrees if necessary */
335 if (proj_angular_output(pj->pj, PJ_FWD))
336 {
337 for (i = 0; i < pa->npoints; i++)
338 {
339 getPoint4d_p(pa, i, &p);
340 to_dec(&p);
341 }
342 }
343
344 return LW_SUCCESS;
345 }
346
347 #endif
348
349 /**
350 * Transform given LWGEOM geometry
351 * from inpj projection to outpj projection
352 */
353 int
lwgeom_transform(LWGEOM * geom,LWPROJ * pj)354 lwgeom_transform(LWGEOM *geom, LWPROJ *pj)
355 {
356 uint32_t i;
357
358 /* No points to transform in an empty! */
359 if ( lwgeom_is_empty(geom) )
360 return LW_SUCCESS;
361
362 switch(geom->type)
363 {
364 case POINTTYPE:
365 case LINETYPE:
366 case CIRCSTRINGTYPE:
367 case TRIANGLETYPE:
368 {
369 LWLINE *g = (LWLINE*)geom;
370 if ( ! ptarray_transform(g->points, pj) ) return LW_FAILURE;
371 break;
372 }
373 case POLYGONTYPE:
374 {
375 LWPOLY *g = (LWPOLY*)geom;
376 for ( i = 0; i < g->nrings; i++ )
377 {
378 if ( ! ptarray_transform(g->rings[i], pj) ) return LW_FAILURE;
379 }
380 break;
381 }
382 case MULTIPOINTTYPE:
383 case MULTILINETYPE:
384 case MULTIPOLYGONTYPE:
385 case COLLECTIONTYPE:
386 case COMPOUNDTYPE:
387 case CURVEPOLYTYPE:
388 case MULTICURVETYPE:
389 case MULTISURFACETYPE:
390 case POLYHEDRALSURFACETYPE:
391 case TINTYPE:
392 {
393 LWCOLLECTION *g = (LWCOLLECTION*)geom;
394 for ( i = 0; i < g->ngeoms; i++ )
395 {
396 if ( ! lwgeom_transform(g->geoms[i], pj) ) return LW_FAILURE;
397 }
398 break;
399 }
400 default:
401 {
402 lwerror("lwgeom_transform: Cannot handle type '%s'",
403 lwtype_name(geom->type));
404 return LW_FAILURE;
405 }
406 }
407 return LW_SUCCESS;
408 }
409