1 /*
2 
3  gg_transform.c -- Gaia PROJ.4 wrapping
4 
5  version 5.0, 2020 August 1
6 
7  Author: Sandro Furieri a.furieri@lqt.it
8 
9  ------------------------------------------------------------------------------
10 
11  Version: MPL 1.1/GPL 2.0/LGPL 2.1
12 
13  The contents of this file are subject to the Mozilla Public License Version
14  1.1 (the "License"); you may not use this file except in compliance with
15  the License. You may obtain a copy of the License at
16  http://www.mozilla.org/MPL/
17 
18 Software distributed under the License is distributed on an "AS IS" basis,
19 WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
20 for the specific language governing rights and limitations under the
21 License.
22 
23 The Original Code is the SpatiaLite library
24 
25 The Initial Developer of the Original Code is Alessandro Furieri
26 
27 Portions created by the Initial Developer are Copyright (C) 2008-2021
28 the Initial Developer. All Rights Reserved.
29 
30 Contributor(s):
31 
32 Alternatively, the contents of this file may be used under the terms of
33 either the GNU General Public License Version 2 or later (the "GPL"), or
34 the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
35 in which case the provisions of the GPL or the LGPL are applicable instead
36 of those above. If you wish to allow use of your version of this file only
37 under the terms of either the GPL or the LGPL, and not to allow others to
38 use your version of this file under the terms of the MPL, indicate your
39 decision by deleting the provisions above and replace them with the notice
40 and other provisions required by the GPL or the LGPL. If you do not delete
41 the provisions above, a recipient may use your version of this file under
42 the terms of any one of the MPL, the GPL or the LGPL.
43 
44 */
45 
46 #include <sys/types.h>
47 #include <stdlib.h>
48 #include <stdio.h>
49 #include <string.h>
50 #include <math.h>
51 
52 #if defined(_WIN32)
53 #include <windows.h>
54 #include <libloaderapi.h>
55 #endif
56 
57 #if defined(_WIN32) && !defined(__MINGW32__)
58 #include "config-msvc.h"
59 #else
60 #include "config.h"
61 #endif
62 
63 #include <spatialite/sqlite.h>
64 #include <spatialite/debug.h>
65 #include <spatialite_private.h>
66 #include <spatialite.h>
67 
68 #include <spatialite/gaiageo.h>
69 
70 #ifndef OMIT_PROJ		/* including PROJ.4 */
71 #ifdef PROJ_NEW			/* supporting PROJ.6 */
72 #include <proj.h>
73 #else /* supporting old PROJ.4 */
74 #include <proj_api.h>
75 #endif
76 
77 #ifndef OMIT_PROJ		/* including PROJ.4 */
78 
79 #ifndef PROJ_NEW		/* supporting old PROJ.4 */
80 static int
gaiaIsLongLat(const char * str)81 gaiaIsLongLat (const char *str)
82 {
83 /* checks if we have to do with ANGLES if +proj=longlat is defined */
84     if (str == NULL)
85 	return 0;
86     if (strstr (str, "+proj=longlat") != NULL)
87 	return 1;
88     return 0;
89 }
90 #endif
91 
92 #ifdef PROJ_NEW			/* only if PROJ.6 is supported */
93 
94 GAIAGEO_DECLARE void
gaiaResetProjErrorMsg_r(const void * p_cache)95 gaiaResetProjErrorMsg_r (const void *p_cache)
96 {
97 /* resets the PROJ error message */
98     struct splite_internal_cache *cache =
99 	(struct splite_internal_cache *) p_cache;
100     if (cache != NULL)
101       {
102 	  if (cache->magic1 == SPATIALITE_CACHE_MAGIC1
103 	      && cache->magic2 == SPATIALITE_CACHE_MAGIC2)
104 	    {
105 		if (cache->gaia_proj_error_msg != NULL)
106 		    sqlite3_free (cache->gaia_proj_error_msg);
107 		cache->gaia_proj_error_msg = NULL;
108 	    }
109       }
110 }
111 
112 GAIAGEO_DECLARE const char *
gaiaSetProjDatabasePath(const void * p_cache,const char * path)113 gaiaSetProjDatabasePath (const void *p_cache, const char *path)
114 {
115 /* return the currently set PATH leading to the private PROJ.6 database */
116     struct splite_internal_cache *cache =
117 	(struct splite_internal_cache *) p_cache;
118     if (cache != NULL)
119       {
120 	  if (cache->magic1 == SPATIALITE_CACHE_MAGIC1
121 	      && cache->magic2 == SPATIALITE_CACHE_MAGIC2)
122 	    {
123 		if (!proj_context_set_database_path
124 		    (cache->PROJ_handle, path, NULL, NULL))
125 		    return NULL;
126 		return proj_context_get_database_path (cache->PROJ_handle);
127 	    }
128       }
129     return NULL;
130 }
131 
132 GAIAGEO_DECLARE const char *
gaiaGetProjDatabasePath(const void * p_cache)133 gaiaGetProjDatabasePath (const void *p_cache)
134 {
135 /* return the currently set PATH leading to the private PROJ.6 database */
136     struct splite_internal_cache *cache =
137 	(struct splite_internal_cache *) p_cache;
138     if (cache != NULL)
139       {
140 	  if (cache->magic1 == SPATIALITE_CACHE_MAGIC1
141 	      && cache->magic2 == SPATIALITE_CACHE_MAGIC2)
142 	      return proj_context_get_database_path (cache->PROJ_handle);
143       }
144     return NULL;
145 }
146 
147 GAIAGEO_DECLARE const char *
gaiaGetProjErrorMsg_r(const void * p_cache)148 gaiaGetProjErrorMsg_r (const void *p_cache)
149 {
150 /* return the latest PROJ error message */
151     struct splite_internal_cache *cache =
152 	(struct splite_internal_cache *) p_cache;
153     if (cache != NULL)
154       {
155 	  if (cache->magic1 == SPATIALITE_CACHE_MAGIC1
156 	      && cache->magic2 == SPATIALITE_CACHE_MAGIC2)
157 	    {
158 		if (cache->gaia_proj_error_msg != NULL)
159 		    return cache->gaia_proj_error_msg;
160 	    }
161       }
162     return NULL;
163 }
164 
165 GAIAGEO_DECLARE void
gaiaSetProjErrorMsg_r(const void * p_cache,const char * msg)166 gaiaSetProjErrorMsg_r (const void *p_cache, const char *msg)
167 {
168 /* setting the latest PROJ error message */
169     struct splite_internal_cache *cache =
170 	(struct splite_internal_cache *) p_cache;
171     if (cache != NULL)
172       {
173 	  if (cache->magic1 == SPATIALITE_CACHE_MAGIC1
174 	      && cache->magic2 == SPATIALITE_CACHE_MAGIC2)
175 	    {
176 		if (cache->gaia_proj_error_msg != NULL)
177 		    sqlite3_free (cache->gaia_proj_error_msg);
178 		cache->gaia_proj_error_msg = sqlite3_mprintf ("%s", msg);
179 	    }
180       }
181 }
182 #endif
183 
184 
185 GAIAGEO_DECLARE double
gaiaRadsToDegs(double rads)186 gaiaRadsToDegs (double rads)
187 {
188 /* converts an ANGLE from radians to degrees */
189 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
190     return proj_todeg (rads);
191 #else /* supporting old PROJ.4 */
192     return rads * RAD_TO_DEG;
193 #endif
194 }
195 
196 GAIAGEO_DECLARE double
gaiaDegsToRads(double degs)197 gaiaDegsToRads (double degs)
198 {
199 /* converts an ANGLE from degrees to radians */
200 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
201     return proj_torad (degs);
202 #else /* supporting old PROJ.4 */
203     return degs * DEG_TO_RAD;
204 #endif
205 }
206 
207 static int
do_transfom_proj(gaiaGeomCollPtr org,gaiaGeomCollPtr dst,int ignore_z,int ignore_m,int from_radians,int to_radians,void * xfrom,void * xto,void * x_from_to)208 do_transfom_proj (gaiaGeomCollPtr org, gaiaGeomCollPtr dst, int ignore_z,
209 		  int ignore_m, int from_radians, int to_radians, void *xfrom,
210 		  void *xto, void *x_from_to)
211 {
212 /* supporting either old PROJ.4 or new PROJ.6 */
213     int error = 0;
214     int ret;
215 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
216     int result_count;
217 #endif
218     int ib;
219     int cnt;
220     int i;
221     double *xx;
222     double *yy;
223     double *zz;
224     double *old_zz = NULL;
225     double *mm = NULL;
226     double *old_mm = NULL;
227     double x;
228     double y;
229     double z;
230     double m;
231     gaiaPointPtr pt;
232     gaiaLinestringPtr ln;
233     gaiaLinestringPtr dst_ln;
234     gaiaPolygonPtr pg;
235     gaiaPolygonPtr dst_pg;
236     gaiaRingPtr rng;
237     gaiaRingPtr dst_rng;
238 #ifdef PROJ_NEW			/* only if new PROJ.6 is supported */
239     PJ *from_to_cs = (PJ *) x_from_to;
240     if (xfrom == NULL)
241 	xfrom = NULL;		/* silencing stupid compiler warnings - unused arg */
242     if (xto == NULL)
243 	xto = NULL;		/* silencing stupid compiler warnings - unused arg */
244 #else /* only id old PROJ.4 is supported */
245     projPJ from_cs = (projPJ) xfrom;
246     projPJ to_cs = (projPJ) xto;
247     if (x_from_to == NULL)
248 	x_from_to = NULL;	/* silencing stupid compiler warnings - unused arg */
249 #endif
250 
251     cnt = 0;
252     pt = org->FirstPoint;
253     while (pt)
254       {
255 	  /* counting POINTs */
256 	  cnt++;
257 	  pt = pt->Next;
258       }
259     if (cnt)
260       {
261 	  /* reprojecting POINTs */
262 	  xx = malloc (sizeof (double) * cnt);
263 	  yy = malloc (sizeof (double) * cnt);
264 	  zz = malloc (sizeof (double) * cnt);
265 	  mm = malloc (sizeof (double) * cnt);
266 	  for (i = 0; i < cnt; i++)
267 	    {
268 		/* zeroing unused coordinates */
269 		if (ignore_z || org->DimensionModel == GAIA_XY
270 		    || org->DimensionModel == GAIA_XY_M)
271 		    zz[i] = 0.0;
272 		if (ignore_m || org->DimensionModel == GAIA_XY
273 		    || org->DimensionModel == GAIA_XY_Z)
274 		    mm[i] = 0.0;
275 	    }
276 	  if (ignore_z
277 	      && (org->DimensionModel == GAIA_XY_Z
278 		  || org->DimensionModel == GAIA_XY_Z_M))
279 	      old_zz = malloc (sizeof (double) * cnt);
280 	  if (ignore_m
281 	      && (org->DimensionModel == GAIA_XY_M
282 		  || org->DimensionModel == GAIA_XY_Z_M))
283 	      old_mm = malloc (sizeof (double) * cnt);
284 	  i = 0;
285 	  pt = org->FirstPoint;
286 	  while (pt)
287 	    {
288 		/* inserting points to be converted in temporary arrays */
289 		if (from_radians)
290 		  {
291 		      xx[i] = gaiaDegsToRads (pt->X);
292 		      yy[i] = gaiaDegsToRads (pt->Y);
293 		  }
294 		else
295 		  {
296 		      xx[i] = pt->X;
297 		      yy[i] = pt->Y;
298 		  }
299 		if (org->DimensionModel == GAIA_XY_Z
300 		    || org->DimensionModel == GAIA_XY_Z_M)
301 		  {
302 		      if (ignore_z)
303 			  old_zz[i] = pt->Z;
304 		      else
305 			  zz[i] = pt->Z;
306 		  }
307 		if (org->DimensionModel == GAIA_XY_M
308 		    || org->DimensionModel == GAIA_XY_Z_M)
309 		  {
310 		      if (ignore_m)
311 			  old_mm[i] = pt->M;
312 		      else
313 			  mm[i] = pt->M;
314 		  }
315 		i++;
316 		pt = pt->Next;
317 	    }
318 	  /* applying reprojection */
319 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
320 	  result_count =
321 	      proj_trans_generic (from_to_cs, PJ_FWD, xx, sizeof (double), cnt,
322 				  yy, sizeof (double), cnt, zz, sizeof (double),
323 				  cnt, mm, sizeof (double), cnt);
324 	  if (result_count == cnt)
325 	      ret = 0;
326 	  else
327 	      ret = 1;
328 #else /* supporting old PROJ.4 */
329 	  ret = pj_transform (from_cs, to_cs, cnt, 0, xx, yy, zz);
330 #endif
331 	  if (ret == 0)
332 	    {
333 		/* inserting the reprojected POINTs in the new GEOMETRY */
334 		for (i = 0; i < cnt; i++)
335 		  {
336 		      if (to_radians)
337 			{
338 			    x = gaiaRadsToDegs (xx[i]);
339 			    y = gaiaRadsToDegs (yy[i]);
340 			}
341 		      else
342 			{
343 			    x = xx[i];
344 			    y = yy[i];
345 			}
346 		      if (org->DimensionModel == GAIA_XY_Z
347 			  || org->DimensionModel == GAIA_XY_Z_M)
348 			{
349 			    if (ignore_z)
350 				z = old_zz[i];
351 			    else
352 				z = zz[i];
353 			}
354 		      if (org->DimensionModel == GAIA_XY_M
355 			  || org->DimensionModel == GAIA_XY_Z_M)
356 			{
357 			    if (ignore_m)
358 				m = old_mm[i];
359 			    else
360 				m = mm[i];
361 			}
362 		      if (dst->DimensionModel == GAIA_XY_Z)
363 			  gaiaAddPointToGeomCollXYZ (dst, x, y, z);
364 		      else if (dst->DimensionModel == GAIA_XY_M)
365 			  gaiaAddPointToGeomCollXYM (dst, x, y, m);
366 		      else if (dst->DimensionModel == GAIA_XY_Z_M)
367 			  gaiaAddPointToGeomCollXYZM (dst, x, y, z, m);
368 		      else
369 			  gaiaAddPointToGeomColl (dst, x, y);
370 		  }
371 	    }
372 	  else
373 	      error = 1;
374 	  free (xx);
375 	  free (yy);
376 	  free (zz);
377 	  free (mm);
378 	  if (old_zz != NULL)
379 	      free (old_zz);
380 	  if (old_mm != NULL)
381 	      free (old_mm);
382 	  xx = NULL;
383 	  yy = NULL;
384 	  zz = NULL;
385 	  mm = NULL;
386 	  old_zz = NULL;
387 	  old_mm = NULL;
388       }
389     if (error)
390 	goto stop;
391     ln = org->FirstLinestring;
392     while (ln)
393       {
394 	  /* reprojecting LINESTRINGs */
395 	  cnt = ln->Points;
396 	  xx = malloc (sizeof (double) * cnt);
397 	  yy = malloc (sizeof (double) * cnt);
398 	  zz = malloc (sizeof (double) * cnt);
399 	  mm = malloc (sizeof (double) * cnt);
400 	  for (i = 0; i < cnt; i++)
401 	    {
402 		/* zeroing unused coordinates */
403 		if (ignore_z || org->DimensionModel == GAIA_XY
404 		    || org->DimensionModel == GAIA_XY_M)
405 		    zz[i] = 0.0;
406 		if (ignore_m || org->DimensionModel == GAIA_XY
407 		    || org->DimensionModel == GAIA_XY_Z)
408 		    mm[i] = 0.0;
409 	    }
410 	  if (ignore_z
411 	      && (org->DimensionModel == GAIA_XY_Z
412 		  || org->DimensionModel == GAIA_XY_Z_M))
413 	      old_zz = malloc (sizeof (double) * cnt);
414 	  if (ignore_m
415 	      && (org->DimensionModel == GAIA_XY_M
416 		  || org->DimensionModel == GAIA_XY_Z_M))
417 	      old_mm = malloc (sizeof (double) * cnt);
418 	  for (i = 0; i < cnt; i++)
419 	    {
420 		/* inserting points to be converted in temporary arrays */
421 		z = 0.0;
422 		m = 0.0;
423 		if (ln->DimensionModel == GAIA_XY_Z)
424 		  {
425 		      gaiaGetPointXYZ (ln->Coords, i, &x, &y, &z);
426 		  }
427 		else if (ln->DimensionModel == GAIA_XY_M)
428 		  {
429 		      gaiaGetPointXYM (ln->Coords, i, &x, &y, &m);
430 		  }
431 		else if (ln->DimensionModel == GAIA_XY_Z_M)
432 		  {
433 		      gaiaGetPointXYZM (ln->Coords, i, &x, &y, &z, &m);
434 		  }
435 		else
436 		  {
437 		      gaiaGetPoint (ln->Coords, i, &x, &y);
438 		  }
439 		if (from_radians)
440 		  {
441 		      xx[i] = gaiaDegsToRads (x);
442 		      yy[i] = gaiaDegsToRads (y);
443 		  }
444 		else
445 		  {
446 		      xx[i] = x;
447 		      yy[i] = y;
448 		  }
449 		if (ln->DimensionModel == GAIA_XY_Z
450 		    || ln->DimensionModel == GAIA_XY_Z_M)
451 		  {
452 		      if (ignore_z)
453 			  old_zz[i] = z;
454 		      else
455 			  zz[i] = z;
456 		  }
457 		if (ln->DimensionModel == GAIA_XY_M
458 		    || ln->DimensionModel == GAIA_XY_Z_M)
459 		  {
460 		      if (ignore_m)
461 			  old_mm[i] = m;
462 		      else
463 			  mm[i] = m;
464 		  }
465 	    }
466 	  /* applying reprojection */
467 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
468 	  result_count =
469 	      proj_trans_generic (from_to_cs, PJ_FWD, xx, sizeof (double), cnt,
470 				  yy, sizeof (double), cnt, zz, sizeof (double),
471 				  cnt, mm, sizeof (double), cnt);
472 	  if (result_count == cnt)
473 	      ret = 0;
474 	  else
475 	      ret = 1;
476 #else /* supporting old PROJ.4 */
477 	  ret = pj_transform (from_cs, to_cs, cnt, 0, xx, yy, zz);
478 #endif
479 	  if (ret == 0)
480 	    {
481 		/* inserting the reprojected LINESTRING in the new GEOMETRY */
482 		dst_ln = gaiaAddLinestringToGeomColl (dst, cnt);
483 		for (i = 0; i < cnt; i++)
484 		  {
485 		      /* setting LINESTRING points */
486 		      if (to_radians)
487 			{
488 			    x = gaiaRadsToDegs (xx[i]);
489 			    y = gaiaRadsToDegs (yy[i]);
490 			}
491 		      else
492 			{
493 			    x = xx[i];
494 			    y = yy[i];
495 			}
496 		      if (ln->DimensionModel == GAIA_XY_Z
497 			  || ln->DimensionModel == GAIA_XY_Z_M)
498 			{
499 			    if (ignore_z)
500 				z = old_zz[i];
501 			    else
502 				z = zz[i];
503 			}
504 		      if (ln->DimensionModel == GAIA_XY_M
505 			  || ln->DimensionModel == GAIA_XY_Z_M)
506 			{
507 			    if (ignore_m)
508 				m = old_mm[i];
509 			    else
510 				m = mm[i];
511 			}
512 		      if (dst_ln->DimensionModel == GAIA_XY_Z)
513 			{
514 			    gaiaSetPointXYZ (dst_ln->Coords, i, x, y, z);
515 			}
516 		      else if (dst_ln->DimensionModel == GAIA_XY_M)
517 			{
518 			    gaiaSetPointXYM (dst_ln->Coords, i, x, y, m);
519 			}
520 		      else if (dst_ln->DimensionModel == GAIA_XY_Z_M)
521 			{
522 			    gaiaSetPointXYZM (dst_ln->Coords, i, x, y, z, m);
523 			}
524 		      else
525 			{
526 			    gaiaSetPoint (dst_ln->Coords, i, x, y);
527 			}
528 		  }
529 	    }
530 	  else
531 	      error = 1;
532 	  free (xx);
533 	  free (yy);
534 	  free (zz);
535 	  free (mm);
536 	  if (old_zz != NULL)
537 	      free (old_zz);
538 	  if (old_mm != NULL)
539 	      free (old_mm);
540 	  xx = NULL;
541 	  yy = NULL;
542 	  zz = NULL;
543 	  mm = NULL;
544 	  old_zz = NULL;
545 	  old_mm = NULL;
546 	  if (error)
547 	      goto stop;
548 	  ln = ln->Next;
549       }
550     pg = org->FirstPolygon;
551     while (pg)
552       {
553 	  /* reprojecting POLYGONs */
554 	  rng = pg->Exterior;
555 	  cnt = rng->Points;
556 	  dst_pg = gaiaAddPolygonToGeomColl (dst, cnt, pg->NumInteriors);
557 	  xx = malloc (sizeof (double) * cnt);
558 	  yy = malloc (sizeof (double) * cnt);
559 	  zz = malloc (sizeof (double) * cnt);
560 	  mm = malloc (sizeof (double) * cnt);
561 	  for (i = 0; i < cnt; i++)
562 	    {
563 		/* zeroing unused coordinates */
564 		if (ignore_z || org->DimensionModel == GAIA_XY
565 		    || org->DimensionModel == GAIA_XY_M)
566 		    zz[i] = 0.0;
567 		if (ignore_m || org->DimensionModel == GAIA_XY
568 		    || org->DimensionModel == GAIA_XY_Z)
569 		    mm[i] = 0.0;
570 	    }
571 	  if (ignore_z
572 	      && (org->DimensionModel == GAIA_XY_Z
573 		  || org->DimensionModel == GAIA_XY_Z_M))
574 	      old_zz = malloc (sizeof (double) * cnt);
575 	  if (ignore_m
576 	      && (org->DimensionModel == GAIA_XY_M
577 		  || org->DimensionModel == GAIA_XY_Z_M))
578 	      old_mm = malloc (sizeof (double) * cnt);
579 	  for (i = 0; i < cnt; i++)
580 	    {
581 		/* inserting points to be converted in temporary arrays [EXTERIOR RING] */
582 		z = 0.0;
583 		m = 0.0;
584 		if (rng->DimensionModel == GAIA_XY_Z)
585 		  {
586 		      gaiaGetPointXYZ (rng->Coords, i, &x, &y, &z);
587 		  }
588 		else if (rng->DimensionModel == GAIA_XY_M)
589 		  {
590 		      gaiaGetPointXYM (rng->Coords, i, &x, &y, &m);
591 		  }
592 		else if (rng->DimensionModel == GAIA_XY_Z_M)
593 		  {
594 		      gaiaGetPointXYZM (rng->Coords, i, &x, &y, &z, &m);
595 		  }
596 		else
597 		  {
598 		      gaiaGetPoint (rng->Coords, i, &x, &y);
599 		  }
600 		if (from_radians)
601 		  {
602 		      xx[i] = gaiaDegsToRads (x);
603 		      yy[i] = gaiaDegsToRads (y);
604 		  }
605 		else
606 		  {
607 		      xx[i] = x;
608 		      yy[i] = y;
609 		  }
610 		if (rng->DimensionModel == GAIA_XY_Z
611 		    || rng->DimensionModel == GAIA_XY_Z_M)
612 		  {
613 		      if (ignore_z)
614 			  old_zz[i] = z;
615 		      else
616 			  zz[i] = z;
617 		  }
618 		if (rng->DimensionModel == GAIA_XY_M
619 		    || rng->DimensionModel == GAIA_XY_Z_M)
620 		  {
621 		      if (ignore_m)
622 			  old_mm[i] = m;
623 		      else
624 			  mm[i] = m;
625 		  }
626 	    }
627 	  /* applying reprojection */
628 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
629 	  result_count =
630 	      proj_trans_generic (from_to_cs, PJ_FWD, xx, sizeof (double), cnt,
631 				  yy, sizeof (double), cnt, zz, sizeof (double),
632 				  cnt, mm, sizeof (double), cnt);
633 	  if (result_count == cnt)
634 	      ret = 0;
635 	  else
636 	      ret = 1;
637 #else /* supporting old PROJ.4 */
638 	  ret = pj_transform (from_cs, to_cs, cnt, 0, xx, yy, zz);
639 #endif
640 	  if (ret == 0)
641 	    {
642 		/* inserting the reprojected POLYGON in the new GEOMETRY */
643 		dst_rng = dst_pg->Exterior;
644 		for (i = 0; i < cnt; i++)
645 		  {
646 		      /* setting EXTERIOR RING points */
647 		      if (to_radians)
648 			{
649 			    x = gaiaRadsToDegs (xx[i]);
650 			    y = gaiaRadsToDegs (yy[i]);
651 			}
652 		      else
653 			{
654 			    x = xx[i];
655 			    y = yy[i];
656 			}
657 		      if (rng->DimensionModel == GAIA_XY_Z
658 			  || rng->DimensionModel == GAIA_XY_Z_M)
659 			{
660 			    if (ignore_z)
661 				z = old_zz[i];
662 			    else
663 				z = zz[i];
664 			}
665 		      if (rng->DimensionModel == GAIA_XY_M
666 			  || rng->DimensionModel == GAIA_XY_Z_M)
667 			{
668 			    if (ignore_m)
669 				m = old_mm[i];
670 			    else
671 				m = mm[i];
672 			}
673 		      if (dst_rng->DimensionModel == GAIA_XY_Z)
674 			{
675 			    gaiaSetPointXYZ (dst_rng->Coords, i, x, y, z);
676 			}
677 		      else if (dst_rng->DimensionModel == GAIA_XY_M)
678 			{
679 			    gaiaSetPointXYM (dst_rng->Coords, i, x, y, m);
680 			}
681 		      else if (dst_rng->DimensionModel == GAIA_XY_Z_M)
682 			{
683 			    gaiaSetPointXYZM (dst_rng->Coords, i, x, y, z, m);
684 			}
685 		      else
686 			{
687 			    gaiaSetPoint (dst_rng->Coords, i, x, y);
688 			}
689 		  }
690 	    }
691 	  else
692 	      error = 1;
693 	  free (xx);
694 	  free (yy);
695 	  free (zz);
696 	  free (mm);
697 	  if (old_zz != NULL)
698 	      free (old_zz);
699 	  if (old_mm != NULL)
700 	      free (old_mm);
701 	  xx = NULL;
702 	  yy = NULL;
703 	  zz = NULL;
704 	  mm = NULL;
705 	  old_zz = NULL;
706 	  old_mm = NULL;
707 	  if (error)
708 	      goto stop;
709 	  for (ib = 0; ib < pg->NumInteriors; ib++)
710 	    {
711 		/* processing INTERIOR RINGS */
712 		rng = pg->Interiors + ib;
713 		cnt = rng->Points;
714 		xx = malloc (sizeof (double) * cnt);
715 		yy = malloc (sizeof (double) * cnt);
716 		zz = malloc (sizeof (double) * cnt);
717 		mm = malloc (sizeof (double) * cnt);
718 		for (i = 0; i < cnt; i++)
719 		  {
720 		      /* zeroing unused coordinates */
721 		      if (ignore_z || org->DimensionModel == GAIA_XY
722 			  || org->DimensionModel == GAIA_XY_M)
723 			  zz[i] = 0.0;
724 		      if (ignore_m || org->DimensionModel == GAIA_XY
725 			  || org->DimensionModel == GAIA_XY_Z)
726 			  mm[i] = 0.0;
727 		  }
728 		if (ignore_z
729 		    && (org->DimensionModel == GAIA_XY_Z
730 			|| org->DimensionModel == GAIA_XY_Z_M))
731 		    old_zz = malloc (sizeof (double) * cnt);
732 		if (ignore_m
733 		    && (org->DimensionModel == GAIA_XY_M
734 			|| org->DimensionModel == GAIA_XY_Z_M))
735 		    old_mm = malloc (sizeof (double) * cnt);
736 		for (i = 0; i < cnt; i++)
737 		  {
738 		      /* inserting points to be converted in temporary arrays [INTERIOR RING] */
739 		      z = 0.0;
740 		      m = 0.0;
741 		      if (rng->DimensionModel == GAIA_XY_Z)
742 			{
743 			    gaiaGetPointXYZ (rng->Coords, i, &x, &y, &z);
744 			}
745 		      else if (rng->DimensionModel == GAIA_XY_M)
746 			{
747 			    gaiaGetPointXYM (rng->Coords, i, &x, &y, &m);
748 			}
749 		      else if (rng->DimensionModel == GAIA_XY_Z_M)
750 			{
751 			    gaiaGetPointXYZM (rng->Coords, i, &x, &y, &z, &m);
752 			}
753 		      else
754 			{
755 			    gaiaGetPoint (rng->Coords, i, &x, &y);
756 			}
757 		      if (from_radians)
758 			{
759 			    xx[i] = gaiaDegsToRads (x);
760 			    yy[i] = gaiaDegsToRads (y);
761 			}
762 		      else
763 			{
764 			    xx[i] = x;
765 			    yy[i] = y;
766 			}
767 		      if (rng->DimensionModel == GAIA_XY_Z
768 			  || rng->DimensionModel == GAIA_XY_Z_M)
769 			{
770 			    if (ignore_z)
771 				old_zz[i] = z;
772 			    else
773 				zz[i] = z;
774 			}
775 		      if (rng->DimensionModel == GAIA_XY_M
776 			  || rng->DimensionModel == GAIA_XY_Z_M)
777 			{
778 			    if (ignore_m)
779 				old_mm[i] = m;
780 			    else
781 				mm[i] = m;
782 			}
783 		  }
784 		/* applying reprojection */
785 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
786 		result_count =
787 		    proj_trans_generic (from_to_cs, PJ_FWD, xx, sizeof (double),
788 					cnt, yy, sizeof (double), cnt, zz,
789 					sizeof (double), cnt, mm,
790 					sizeof (double), cnt);
791 		if (result_count == cnt)
792 		    ret = 0;
793 		else
794 		    ret = 1;
795 #else /* supporting old PROJ.4 */
796 		ret = pj_transform (from_cs, to_cs, cnt, 0, xx, yy, zz);
797 #endif
798 		if (ret == 0)
799 		  {
800 		      /* inserting the reprojected POLYGON in the new GEOMETRY */
801 		      dst_rng = gaiaAddInteriorRing (dst_pg, ib, cnt);
802 		      for (i = 0; i < cnt; i++)
803 			{
804 			    /* setting INTERIOR RING points */
805 			    if (to_radians)
806 			      {
807 				  x = gaiaRadsToDegs (xx[i]);
808 				  y = gaiaRadsToDegs (yy[i]);
809 			      }
810 			    else
811 			      {
812 				  x = xx[i];
813 				  y = yy[i];
814 			      }
815 			    if (rng->DimensionModel == GAIA_XY_Z
816 				|| rng->DimensionModel == GAIA_XY_Z_M)
817 			      {
818 				  if (ignore_z)
819 				      z = old_zz[i];
820 				  else
821 				      z = zz[i];
822 			      }
823 			    if (rng->DimensionModel == GAIA_XY_M
824 				|| rng->DimensionModel == GAIA_XY_Z_M)
825 			      {
826 				  if (ignore_m)
827 				      m = old_mm[i];
828 				  else
829 				      m = mm[i];
830 			      }
831 			    if (dst_rng->DimensionModel == GAIA_XY_Z)
832 			      {
833 				  gaiaSetPointXYZ (dst_rng->Coords, i, x, y, z);
834 			      }
835 			    else if (dst_rng->DimensionModel == GAIA_XY_M)
836 			      {
837 				  gaiaSetPointXYM (dst_rng->Coords, i, x, y, m);
838 			      }
839 			    else if (dst_rng->DimensionModel == GAIA_XY_Z_M)
840 			      {
841 				  gaiaSetPointXYZM (dst_rng->Coords, i, x, y, z,
842 						    m);
843 			      }
844 			    else
845 			      {
846 				  gaiaSetPoint (dst_rng->Coords, i, x, y);
847 			      }
848 			}
849 		  }
850 		else
851 		    error = 1;
852 		free (xx);
853 		free (yy);
854 		free (zz);
855 		free (mm);
856 		if (old_zz != NULL)
857 		    free (old_zz);
858 		if (old_mm != NULL)
859 		    free (old_mm);
860 		xx = NULL;
861 		yy = NULL;
862 		zz = NULL;
863 		mm = NULL;
864 		old_zz = NULL;
865 		old_mm = NULL;
866 		if (error)
867 		    goto stop;
868 	    }
869 	  pg = pg->Next;
870       }
871   stop:
872 #endif
873     return error;
874 }
875 
876 static gaiaGeomCollPtr
gaiaTransformCommon(void * x_handle,const void * p_cache,gaiaGeomCollPtr org,const char * proj_string_1,const char * proj_string_2,gaiaProjAreaPtr proj_bbox,int ignore_z,int ignore_m)877 gaiaTransformCommon (void *x_handle, const void *p_cache, gaiaGeomCollPtr org,
878 		     const char *proj_string_1,
879 		     const char *proj_string_2, gaiaProjAreaPtr proj_bbox,
880 		     int ignore_z, int ignore_m)
881 {
882 /* creates a new GEOMETRY reprojecting coordinates from the original one */
883     int error = 0;
884 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
885     PJ_CONTEXT *handle = (PJ_CONTEXT *) x_handle;
886     PJ *from_to_pre;
887     PJ *from_to_cs;
888     int proj_is_cached = 0;
889 #else /* supporting old PROJ.4 */
890     if (p_cache == NULL)
891 	p_cache = NULL;		/* silencing stupid compiler warnings about unused args */
892     projCtx handle = (projCtx) x_handle;
893     projPJ from_cs;
894     projPJ to_cs;
895 #endif
896     int from_radians;
897     int to_radians;
898     gaiaGeomCollPtr dst;
899 
900 /* preliminary validity check */
901 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
902     gaiaResetProjErrorMsg_r (p_cache);
903 #endif
904     if (proj_bbox == NULL)
905 	proj_bbox = NULL;	/* silencing stupid compiler warnings about unused args */
906     if (proj_string_1 == NULL)
907 	return NULL;
908 
909 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
910     if (gaiaCurrentCachedProjMatches
911 	(p_cache, proj_string_1, proj_string_2, proj_bbox))
912       {
913 	  from_to_cs = gaiaGetCurrentCachedProj (p_cache);
914 	  if (from_to_cs != NULL)
915 	    {
916 		proj_is_cached = 1;
917 		goto skip_cached;
918 	    }
919       }
920     if (proj_string_2 != NULL)
921       {
922 	  PJ_AREA *area = NULL;
923 	  if (proj_bbox != NULL)
924 	    {
925 		area = proj_area_create ();
926 		proj_area_set_bbox (area, proj_bbox->WestLongitude,
927 				    proj_bbox->SouthLatitude,
928 				    proj_bbox->EastLongitude,
929 				    proj_bbox->NorthLatitude);
930 	    }
931 	  from_to_pre =
932 	      proj_create_crs_to_crs (handle, proj_string_1, proj_string_2,
933 				      area);
934 	  if (!from_to_pre)
935 	      return NULL;
936 	  from_to_cs = proj_normalize_for_visualization (handle, from_to_pre);
937 	  proj_destroy (from_to_pre);
938 	  if (area != NULL)
939 	      proj_area_destroy (area);
940 	  if (!from_to_cs)
941 	      return NULL;
942 	  proj_is_cached =
943 	      gaiaSetCurrentCachedProj (p_cache, from_to_cs, proj_string_1,
944 					proj_string_2, proj_bbox);
945       }
946     else			/* PROJ_STRING - PIPELINE */
947       {
948 	  from_to_cs = proj_create (handle, proj_string_1);
949 	  if (!from_to_cs)
950 	      return NULL;
951 	  proj_is_cached =
952 	      gaiaSetCurrentCachedProj (p_cache, from_to_cs, proj_string_1,
953 					NULL, NULL);
954       }
955   skip_cached:
956 #else /* supporting old PROJ.4 */
957     if (proj_string_2 == NULL)
958 	return NULL;
959     if (handle != NULL)
960       {
961 	  from_cs = pj_init_plus_ctx (handle, proj_string_1);
962 	  to_cs = pj_init_plus_ctx (handle, proj_string_2);
963       }
964     else
965       {
966 	  from_cs = pj_init_plus (proj_string_1);
967 	  to_cs = pj_init_plus (proj_string_2);
968       }
969     if (!from_cs)
970       {
971 	  if (to_cs)
972 	      pj_free (to_cs);
973 	  return NULL;
974       }
975     if (!to_cs)
976       {
977 	  pj_free (from_cs);
978 	  return NULL;
979       }
980 #endif
981     if (org->DimensionModel == GAIA_XY_Z)
982 	dst = gaiaAllocGeomCollXYZ ();
983     else if (org->DimensionModel == GAIA_XY_M)
984 	dst = gaiaAllocGeomCollXYM ();
985     else if (org->DimensionModel == GAIA_XY_Z_M)
986 	dst = gaiaAllocGeomCollXYZM ();
987     else
988 	dst = gaiaAllocGeomColl ();
989 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
990     from_radians = proj_angular_input (from_to_cs, PJ_FWD);
991     to_radians = proj_angular_output (from_to_cs, PJ_FWD);
992     error =
993 	do_transfom_proj (org, dst, ignore_z, ignore_m, from_radians,
994 			  to_radians, NULL, NULL, from_to_cs);
995 #else /* supporting old PROJ.4 */
996     from_radians = gaiaIsLongLat (proj_string_1);
997     to_radians = gaiaIsLongLat (proj_string_2);
998     error =
999 	do_transfom_proj (org, dst, ignore_z, ignore_m, from_radians,
1000 			  to_radians, from_cs, to_cs, NULL);
1001 #endif
1002 
1003 /* destroying the PROJ4 params */
1004 #ifdef PROJ_NEW			/* supporting new PROJ.6 */
1005     if (!proj_is_cached)
1006 	proj_destroy (from_to_cs);
1007 #else /* supporting old PROJ.4 */
1008     pj_free (from_cs);
1009     pj_free (to_cs);
1010 #endif
1011     if (error)
1012       {
1013 	  /* some error occurred */
1014 	  gaiaPointPtr pP;
1015 	  gaiaPointPtr pPn;
1016 	  gaiaLinestringPtr pL;
1017 	  gaiaLinestringPtr pLn;
1018 	  gaiaPolygonPtr pA;
1019 	  gaiaPolygonPtr pAn;
1020 	  pP = dst->FirstPoint;
1021 	  while (pP != NULL)
1022 	    {
1023 		pPn = pP->Next;
1024 		gaiaFreePoint (pP);
1025 		pP = pPn;
1026 	    }
1027 	  pL = dst->FirstLinestring;
1028 	  while (pL != NULL)
1029 	    {
1030 		pLn = pL->Next;
1031 		gaiaFreeLinestring (pL);
1032 		pL = pLn;
1033 	    }
1034 	  pA = dst->FirstPolygon;
1035 	  while (pA != NULL)
1036 	    {
1037 		pAn = pA->Next;
1038 		gaiaFreePolygon (pA);
1039 		pA = pAn;
1040 	    }
1041 	  dst->FirstPoint = NULL;
1042 	  dst->LastPoint = NULL;
1043 	  dst->FirstLinestring = NULL;
1044 	  dst->LastLinestring = NULL;
1045 	  dst->FirstPolygon = NULL;
1046 	  dst->LastPolygon = NULL;
1047       }
1048     if (dst)
1049       {
1050 	  gaiaMbrGeometry (dst);
1051 	  dst->DeclaredType = org->DeclaredType;
1052       }
1053     return dst;
1054 }
1055 
1056 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaTransform(gaiaGeomCollPtr org,const char * proj_from,const char * proj_to)1057 gaiaTransform (gaiaGeomCollPtr org, const char *proj_from, const char *proj_to)
1058 {
1059     return gaiaTransformEx (org, proj_from, proj_to, NULL);
1060 }
1061 
1062 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaTransform_r(const void * p_cache,gaiaGeomCollPtr org,const char * proj_from,const char * proj_to)1063 gaiaTransform_r (const void *p_cache, gaiaGeomCollPtr org,
1064 		 const char *proj_from, const char *proj_to)
1065 {
1066     return gaiaTransformEx_r (p_cache, org, proj_from, proj_to, NULL);
1067 }
1068 
1069 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaTransformEx(gaiaGeomCollPtr org,const char * proj_string_1,const char * proj_string_2,gaiaProjAreaPtr proj_bbox)1070 gaiaTransformEx (gaiaGeomCollPtr org, const char *proj_string_1,
1071 		 const char *proj_string_2, gaiaProjAreaPtr proj_bbox)
1072 {
1073     return gaiaTransformCommon (NULL, NULL, org, proj_string_1, proj_string_2,
1074 				proj_bbox, 0, 0);
1075 }
1076 
1077 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaTransformEx_r(const void * p_cache,gaiaGeomCollPtr org,const char * proj_string_1,const char * proj_string_2,gaiaProjAreaPtr proj_bbox)1078 gaiaTransformEx_r (const void *p_cache, gaiaGeomCollPtr org,
1079 		   const char *proj_string_1, const char *proj_string_2,
1080 		   gaiaProjAreaPtr proj_bbox)
1081 {
1082     struct splite_internal_cache *cache =
1083 	(struct splite_internal_cache *) p_cache;
1084     void *handle = NULL;
1085     if (cache == NULL)
1086 	return NULL;
1087     if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
1088 	|| cache->magic2 != SPATIALITE_CACHE_MAGIC2)
1089 	return NULL;
1090     handle = cache->PROJ_handle;
1091     if (handle == NULL)
1092 	return NULL;
1093     return gaiaTransformCommon (handle, cache, org, proj_string_1,
1094 				proj_string_2, proj_bbox, 0, 0);
1095 }
1096 
1097 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaTransformXY(gaiaGeomCollPtr org,const char * proj_from,const char * proj_to)1098 gaiaTransformXY (gaiaGeomCollPtr org, const char *proj_from,
1099 		 const char *proj_to)
1100 {
1101     return gaiaTransformCommon (NULL, NULL, org,
1102 				proj_from, proj_to, NULL, 1, 1);
1103 }
1104 
1105 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaTransformXY_r(const void * p_cache,gaiaGeomCollPtr org,const char * proj_from,const char * proj_to)1106 gaiaTransformXY_r (const void *p_cache, gaiaGeomCollPtr org,
1107 		   const char *proj_from, const char *proj_to)
1108 {
1109     struct splite_internal_cache *cache =
1110 	(struct splite_internal_cache *) p_cache;
1111     void *handle = NULL;
1112     if (cache == NULL)
1113 	return NULL;
1114     if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
1115 	|| cache->magic2 != SPATIALITE_CACHE_MAGIC2)
1116 	return NULL;
1117     handle = cache->PROJ_handle;
1118     if (handle == NULL)
1119 	return NULL;
1120     return gaiaTransformCommon (handle, cache, org,
1121 				proj_from, proj_to, NULL, 1, 1);
1122 }
1123 
1124 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaTransformXYZ(gaiaGeomCollPtr org,const char * proj_from,const char * proj_to)1125 gaiaTransformXYZ (gaiaGeomCollPtr org, const char *proj_from,
1126 		  const char *proj_to)
1127 {
1128     return gaiaTransformCommon (NULL, NULL, org,
1129 				proj_from, proj_to, NULL, 0, 1);
1130 }
1131 
1132 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaTransformXYZ_r(const void * p_cache,gaiaGeomCollPtr org,const char * proj_from,const char * proj_to)1133 gaiaTransformXYZ_r (const void *p_cache, gaiaGeomCollPtr org,
1134 		    const char *proj_from, const char *proj_to)
1135 {
1136     struct splite_internal_cache *cache =
1137 	(struct splite_internal_cache *) p_cache;
1138     void *handle = NULL;
1139     if (cache == NULL)
1140 	return NULL;
1141     if (cache->magic1 != SPATIALITE_CACHE_MAGIC1
1142 	|| cache->magic2 != SPATIALITE_CACHE_MAGIC2)
1143 	return NULL;
1144     handle = cache->PROJ_handle;
1145     if (handle == NULL)
1146 	return NULL;
1147     return gaiaTransformCommon (handle, cache, org,
1148 				proj_from, proj_to, NULL, 0, 1);
1149 }
1150 
1151 #ifdef PROJ_NEW			/* only if new PROJ.6 is supported */
1152 GAIAGEO_DECLARE char *
gaiaGetProjString(const void * p_cache,const char * auth_name,int auth_srid)1153 gaiaGetProjString (const void *p_cache, const char *auth_name, int auth_srid)
1154 {
1155 /* return the proj-string expression corresponding to some CRS */
1156     PJ *crs_def;
1157     struct splite_internal_cache *cache =
1158 	(struct splite_internal_cache *) p_cache;
1159     const char *proj_string;
1160     int len;
1161     char *result;
1162     char xsrid[64];
1163 
1164     sprintf (xsrid, "%d", auth_srid);
1165     crs_def =
1166 	proj_create_from_database (cache->PROJ_handle, auth_name, xsrid,
1167 				   PJ_CATEGORY_CRS, 0, NULL);
1168     if (crs_def == NULL)
1169 	return NULL;
1170     proj_string =
1171 	proj_as_proj_string (cache->PROJ_handle, crs_def, PJ_PROJ_5, NULL);
1172     if (proj_string == NULL)
1173       {
1174 	  proj_destroy (crs_def);
1175 	  return NULL;
1176       }
1177     len = strlen (proj_string);
1178     result = malloc (len + 1);
1179     strcpy (result, proj_string);
1180     proj_destroy (crs_def);
1181     return result;
1182 }
1183 
1184 GAIAGEO_DECLARE char *
gaiaGetProjWKT(const void * p_cache,const char * auth_name,int auth_srid,int style,int indented,int indentation)1185 gaiaGetProjWKT (const void *p_cache, const char *auth_name, int auth_srid,
1186 		int style, int indented, int indentation)
1187 {
1188 /* return the WKT expression corresponding to some CRS */
1189     PJ *crs_def;
1190     struct splite_internal_cache *cache =
1191 	(struct splite_internal_cache *) p_cache;
1192     int proj_style;
1193     const char *wkt;
1194     int len;
1195     char *result;
1196     char xsrid[64];
1197     char dummy[64];
1198     char *options[4];
1199     options[1] = dummy;
1200     options[2] = "OUTPUT_AXIS=AUTO";
1201     options[3] = NULL;
1202 
1203     sprintf (xsrid, "%d", auth_srid);
1204     crs_def =
1205 	proj_create_from_database (cache->PROJ_handle, auth_name, xsrid,
1206 				   PJ_CATEGORY_CRS, 0, NULL);
1207     if (crs_def == NULL)
1208 	return NULL;
1209 
1210     switch (style)
1211       {
1212       case GAIA_PROJ_WKT_GDAL:
1213 	  proj_style = PJ_WKT1_GDAL;
1214 	  break;
1215       case GAIA_PROJ_WKT_ESRI:
1216 	  proj_style = PJ_WKT1_ESRI;
1217 	  break;
1218       case GAIA_PROJ_WKT_ISO_2015:
1219 	  proj_style = PJ_WKT2_2015;
1220 	  break;
1221       default:
1222 	  proj_style = PJ_WKT2_2015;
1223 	  break;
1224       };
1225     if (indented)
1226 	options[0] = "MULTILINE=YES";
1227     else
1228 	options[0] = "MULTILINE=NO";
1229     if (indentation < 1)
1230 	indentation = 1;
1231     if (indentation > 8)
1232 	indentation = 8;
1233     sprintf (dummy, "INDENTATION_WIDTH=%d", indentation);
1234     wkt =
1235 	proj_as_wkt (cache->PROJ_handle, crs_def, proj_style,
1236 		     (const char *const *) options);
1237     if (wkt == NULL)
1238       {
1239 	  proj_destroy (crs_def);
1240 	  return NULL;
1241       }
1242     len = strlen (wkt);
1243     result = malloc (len + 1);
1244     strcpy (result, wkt);
1245     proj_destroy (crs_def);
1246     return result;
1247 }
1248 
1249 GAIAGEO_DECLARE int
gaiaGuessSridFromWKT(sqlite3 * sqlite,const void * p_cache,const char * wkt,int * srid)1250 gaiaGuessSridFromWKT (sqlite3 * sqlite, const void *p_cache, const char *wkt,
1251 		      int *srid)
1252 {
1253 /* return the SRID value corresponding to a given WKT expression */
1254     PJ *crs1 = NULL;
1255     PJ *crs2 = NULL;
1256     struct splite_internal_cache *cache =
1257 	(struct splite_internal_cache *) p_cache;
1258     sqlite3_stmt *stmt = NULL;
1259     int ret;
1260     const char *sql;
1261     int xsrid = -1;
1262 
1263 /* sanity check */
1264     if (cache == NULL)
1265 	goto error;
1266     if (cache->PROJ_handle == NULL)
1267 	goto error;
1268 
1269 /* attempting to parse the WKT expression */
1270     crs1 = proj_create_from_wkt (cache->PROJ_handle, wkt, NULL, NULL, NULL);
1271     if (crs1 == NULL)
1272       {
1273 	  spatialite_e
1274 	      ("gaiaGuessSridFromWKT: invalid/malformed WKT expression\n");
1275 	  goto error;
1276       }
1277 
1278     sql = "SELECT srid, Upper(auth_name), auth_srid FROM spatial_ref_sys";
1279     ret = sqlite3_prepare_v2 (sqlite, sql, strlen (sql), &stmt, NULL);
1280     if (ret != SQLITE_OK)
1281       {
1282 	  spatialite_e ("gaiaGuessSridFromWKT: %s\n", sqlite3_errmsg (sqlite));
1283 	  goto error;
1284       }
1285     while (1)
1286       {
1287 	  /* scrolling the result set rows */
1288 	  ret = sqlite3_step (stmt);
1289 	  if (ret == SQLITE_DONE)
1290 	      break;		/* end of result set */
1291 	  if (ret == SQLITE_ROW)
1292 	    {
1293 		char dummy[64];
1294 		int srid = sqlite3_column_int (stmt, 0);
1295 		const char *auth_name =
1296 		    (const char *) sqlite3_column_text (stmt, 1);
1297 		int auth_srid = sqlite3_column_int (stmt, 2);
1298 		sprintf (dummy, "%d", auth_srid);
1299 		/* parsing some CRS */
1300 		crs2 =
1301 		    proj_create_from_database (cache->PROJ_handle, auth_name,
1302 					       dummy, PJ_CATEGORY_CRS, 0, NULL);
1303 		if (crs2 != NULL)
1304 		  {
1305 		      /* ok, it's a valid CRS - comparing */
1306 		      if (proj_is_equivalent_to
1307 			  (crs1, crs2,
1308 			   PJ_COMP_EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS))
1309 			{
1310 			    xsrid = srid;
1311 			    proj_destroy (crs2);
1312 			    break;
1313 			}
1314 		      proj_destroy (crs2);
1315 		  }
1316 	    }
1317       }
1318     sqlite3_finalize (stmt);
1319     proj_destroy (crs1);
1320     *srid = xsrid;
1321     gaiaResetProjErrorMsg_r (cache);
1322     return 1;
1323 
1324   error:
1325     if (stmt != NULL)
1326 	sqlite3_finalize (stmt);
1327     if (crs1 != NULL)
1328 	proj_destroy (crs1);
1329     *srid = xsrid;
1330     return 0;
1331 }
1332 #endif
1333 
1334 #endif /* end including PROJ.4 */
1335 
1336 GAIAGEO_DECLARE void
gaiaShiftLongitude(gaiaGeomCollPtr geom)1337 gaiaShiftLongitude (gaiaGeomCollPtr geom)
1338 {
1339 /* returns a geometry that is the old geometry with negative longitudes shift by 360 */
1340     int ib;
1341     int iv;
1342     double x;
1343     double y;
1344     double z;
1345     double m;
1346     gaiaPointPtr point;
1347     gaiaPolygonPtr polyg;
1348     gaiaLinestringPtr line;
1349     gaiaRingPtr ring;
1350     if (!geom)
1351 	return;
1352     point = geom->FirstPoint;
1353     while (point)
1354       {
1355 	  /* shifting POINTs */
1356 	  if (point->X < 0)
1357 	    {
1358 		point->X += 360.0;
1359 	    }
1360 	  point = point->Next;
1361       }
1362     line = geom->FirstLinestring;
1363     while (line)
1364       {
1365 	  /* shifting LINESTRINGs */
1366 	  for (iv = 0; iv < line->Points; iv++)
1367 	    {
1368 		m = 0.0;
1369 		z = 0.0;
1370 		if (line->DimensionModel == GAIA_XY_Z)
1371 		  {
1372 		      gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
1373 		  }
1374 		else if (line->DimensionModel == GAIA_XY_M)
1375 		  {
1376 		      gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
1377 		  }
1378 		else if (line->DimensionModel == GAIA_XY_Z_M)
1379 		  {
1380 		      gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
1381 		  }
1382 		else
1383 		  {
1384 		      gaiaGetPoint (line->Coords, iv, &x, &y);
1385 		  }
1386 		if (x < 0)
1387 		  {
1388 		      x += 360.0;
1389 		  }
1390 		if (line->DimensionModel == GAIA_XY_Z)
1391 		  {
1392 		      gaiaSetPointXYZ (line->Coords, iv, x, y, z);
1393 		  }
1394 		else if (line->DimensionModel == GAIA_XY_M)
1395 		  {
1396 		      gaiaSetPointXYM (line->Coords, iv, x, y, m);
1397 		  }
1398 		else if (line->DimensionModel == GAIA_XY_Z_M)
1399 		  {
1400 		      gaiaSetPointXYZM (line->Coords, iv, x, y, z, m);
1401 		  }
1402 		else
1403 		  {
1404 		      gaiaSetPoint (line->Coords, iv, x, y);
1405 		  }
1406 	    }
1407 	  line = line->Next;
1408       }
1409     polyg = geom->FirstPolygon;
1410     while (polyg)
1411       {
1412 	  /* shifting POLYGONs */
1413 	  m = 0.0;
1414 	  z = 0.0;
1415 	  ring = polyg->Exterior;
1416 	  for (iv = 0; iv < ring->Points; iv++)
1417 	    {
1418 		/* shifting the EXTERIOR RING */
1419 		if (ring->DimensionModel == GAIA_XY_Z)
1420 		  {
1421 		      gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
1422 		  }
1423 		else if (ring->DimensionModel == GAIA_XY_M)
1424 		  {
1425 		      gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
1426 		  }
1427 		else if (ring->DimensionModel == GAIA_XY_Z_M)
1428 		  {
1429 		      gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
1430 		  }
1431 		else
1432 		  {
1433 		      gaiaGetPoint (ring->Coords, iv, &x, &y);
1434 		  }
1435 		if (x < 0)
1436 		  {
1437 		      x += 360.0;
1438 		  }
1439 		if (ring->DimensionModel == GAIA_XY_Z)
1440 		  {
1441 		      gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
1442 		  }
1443 		else if (ring->DimensionModel == GAIA_XY_M)
1444 		  {
1445 		      gaiaSetPointXYM (ring->Coords, iv, x, y, m);
1446 		  }
1447 		else if (ring->DimensionModel == GAIA_XY_Z_M)
1448 		  {
1449 		      gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
1450 		  }
1451 		else
1452 		  {
1453 		      gaiaSetPoint (ring->Coords, iv, x, y);
1454 		  }
1455 	    }
1456 	  for (ib = 0; ib < polyg->NumInteriors; ib++)
1457 	    {
1458 		/* shifting the INTERIOR RINGs */
1459 		ring = polyg->Interiors + ib;
1460 		for (iv = 0; iv < ring->Points; iv++)
1461 		  {
1462 		      if (ring->DimensionModel == GAIA_XY_Z)
1463 			{
1464 			    gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
1465 			}
1466 		      else if (ring->DimensionModel == GAIA_XY_M)
1467 			{
1468 			    gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
1469 			}
1470 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
1471 			{
1472 			    gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
1473 			}
1474 		      else
1475 			{
1476 			    gaiaGetPoint (ring->Coords, iv, &x, &y);
1477 			}
1478 		      if (x < 0)
1479 			{
1480 			    x += 360.0;
1481 			}
1482 		      if (ring->DimensionModel == GAIA_XY_Z)
1483 			{
1484 			    gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
1485 			}
1486 		      else if (ring->DimensionModel == GAIA_XY_M)
1487 			{
1488 			    gaiaSetPointXYM (ring->Coords, iv, x, y, m);
1489 			}
1490 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
1491 			{
1492 			    gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
1493 			}
1494 		      else
1495 			{
1496 			    gaiaSetPoint (ring->Coords, iv, x, y);
1497 			}
1498 		  }
1499 	    }
1500 	  polyg = polyg->Next;
1501       }
1502     gaiaMbrGeometry (geom);
1503 }
1504 
1505 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaMakeCircle(double cx,double cy,double radius,double step)1506 gaiaMakeCircle (double cx, double cy, double radius, double step)
1507 {
1508 /* return a Linestring approximating a Circle */
1509     return gaiaMakeEllipse (cx, cy, radius, radius, step);
1510 }
1511 
1512 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaMakeArc(double cx,double cy,double radius,double start,double stop,double step)1513 gaiaMakeArc (double cx,
1514 	     double cy, double radius, double start, double stop, double step)
1515 {
1516 /* return a Linestring approximating a Circular Arc */
1517     return gaiaMakeEllipticArc (cx, cy, radius, radius, start, stop, step);
1518 }
1519 
1520 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaMakeEllipse(double cx,double cy,double x_axis,double y_axis,double step)1521 gaiaMakeEllipse (double cx, double cy, double x_axis, double y_axis,
1522 		 double step)
1523 {
1524 /* return a Linestring approximating an Ellipse */
1525     gaiaDynamicLinePtr dyn;
1526     double x;
1527     double y;
1528     double angle = 0.0;
1529     int points = 0;
1530     gaiaPointPtr pt;
1531     gaiaGeomCollPtr geom;
1532     gaiaLinestringPtr ln;
1533     int iv = 0;
1534 
1535     if (step < 0.0)
1536 	step *= -1.0;
1537     if (step == 0.0)
1538 	step = 10.0;
1539     if (step < 0.1)
1540 	step = 0.1;
1541     if (step > 45.0)
1542 	step = 45.0;
1543     if (x_axis < 0.0)
1544 	x_axis *= -1.0;
1545     if (y_axis < 0.0)
1546 	y_axis *= -1.0;
1547     dyn = gaiaAllocDynamicLine ();
1548     while (angle < 360.0)
1549       {
1550 	  double rads = angle * .0174532925199432958;
1551 	  x = cx + (x_axis * cos (rads));
1552 	  y = cy + (y_axis * sin (rads));
1553 	  gaiaAppendPointToDynamicLine (dyn, x, y);
1554 	  angle += step;
1555       }
1556 /* closing the ellipse */
1557     gaiaAppendPointToDynamicLine (dyn, dyn->First->X, dyn->First->Y);
1558 
1559     pt = dyn->First;
1560     while (pt)
1561       {
1562 	  /* counting how many points */
1563 	  points++;
1564 	  pt = pt->Next;
1565       }
1566     if (points == 0)
1567       {
1568 	  gaiaFreeDynamicLine (dyn);
1569 	  return NULL;
1570       }
1571 
1572     geom = gaiaAllocGeomColl ();
1573     ln = gaiaAddLinestringToGeomColl (geom, points);
1574     pt = dyn->First;
1575     while (pt)
1576       {
1577 	  /* setting Vertices */
1578 	  gaiaSetPoint (ln->Coords, iv, pt->X, pt->Y);
1579 	  iv++;
1580 	  pt = pt->Next;
1581       }
1582     gaiaFreeDynamicLine (dyn);
1583     return geom;
1584 }
1585 
1586 GAIAGEO_DECLARE gaiaGeomCollPtr
gaiaMakeEllipticArc(double cx,double cy,double x_axis,double y_axis,double start,double stop,double step)1587 gaiaMakeEllipticArc (double cx,
1588 		     double cy, double x_axis, double y_axis, double start,
1589 		     double stop, double step)
1590 {
1591 /* return a Linestring approximating an Elliptic Arc */
1592     gaiaDynamicLinePtr dyn;
1593     double x;
1594     double y;
1595     double angle;
1596     double rads;
1597     int points = 0;
1598     gaiaPointPtr pt;
1599     gaiaGeomCollPtr geom;
1600     gaiaLinestringPtr ln;
1601     int iv = 0;
1602 
1603     if (step < 0.0)
1604 	step *= -1.0;
1605     if (step == 0.0)
1606 	step = 10.0;
1607     if (step < 0.1)
1608 	step = 0.1;
1609     if (step > 45.0)
1610 	step = 45.0;
1611     if (x_axis < 0.0)
1612 	x_axis *= -1.0;
1613     if (y_axis < 0.0)
1614 	y_axis *= -1.0;
1615 /* normalizing Start/Stop angles */
1616     while (start >= 360.0)
1617 	start -= 360.0;
1618     while (start <= -720.0)
1619 	start += 360;
1620     while (stop >= 360.0)
1621 	stop -= 360.0;
1622     while (stop <= -720.0)
1623 	stop += 360;
1624     if (start < 0.0)
1625 	start = 360.0 + start;
1626     if (stop < 0.0)
1627 	stop = 360.0 + stop;
1628     if (start > stop)
1629 	stop += 360.0;
1630 
1631     dyn = gaiaAllocDynamicLine ();
1632     angle = start;
1633     while (angle < stop)
1634       {
1635 	  rads = angle * .0174532925199432958;
1636 	  x = cx + (x_axis * cos (rads));
1637 	  y = cy + (y_axis * sin (rads));
1638 	  gaiaAppendPointToDynamicLine (dyn, x, y);
1639 	  angle += step;
1640 	  points++;
1641       }
1642     if (points == 0)
1643       {
1644 	  gaiaFreeDynamicLine (dyn);
1645 	  return NULL;
1646       }
1647 /* closing the arc */
1648     rads = stop * .0174532925199432958;
1649     x = cx + (x_axis * cos (rads));
1650     y = cy + (y_axis * sin (rads));
1651     if (x != dyn->Last->X || y != dyn->Last->Y)
1652 	gaiaAppendPointToDynamicLine (dyn, x, y);
1653 
1654     pt = dyn->First;
1655     points = 0;
1656     while (pt)
1657       {
1658 	  /* counting how many points */
1659 	  points++;
1660 	  pt = pt->Next;
1661       }
1662     if (points == 0)
1663       {
1664 	  gaiaFreeDynamicLine (dyn);
1665 	  return NULL;
1666       }
1667 
1668     geom = gaiaAllocGeomColl ();
1669     ln = gaiaAddLinestringToGeomColl (geom, points);
1670     pt = dyn->First;
1671     while (pt)
1672       {
1673 	  /* setting Vertices */
1674 	  gaiaSetPoint (ln->Coords, iv, pt->X, pt->Y);
1675 	  iv++;
1676 	  pt = pt->Next;
1677       }
1678     gaiaFreeDynamicLine (dyn);
1679     return geom;
1680 }
1681 
1682 GAIAGEO_DECLARE void
gaiaShiftCoords(gaiaGeomCollPtr geom,double shift_x,double shift_y)1683 gaiaShiftCoords (gaiaGeomCollPtr geom, double shift_x, double shift_y)
1684 {
1685 /* returns a geometry that is the old geometry with required shifting applied to coordinates */
1686     int ib;
1687     int iv;
1688     double x;
1689     double y;
1690     double z;
1691     double m;
1692     gaiaPointPtr point;
1693     gaiaPolygonPtr polyg;
1694     gaiaLinestringPtr line;
1695     gaiaRingPtr ring;
1696     if (!geom)
1697 	return;
1698     point = geom->FirstPoint;
1699     while (point)
1700       {
1701 	  /* shifting POINTs */
1702 	  point->X += shift_x;
1703 	  point->Y += shift_y;
1704 	  point = point->Next;
1705       }
1706     line = geom->FirstLinestring;
1707     while (line)
1708       {
1709 	  /* shifting LINESTRINGs */
1710 	  for (iv = 0; iv < line->Points; iv++)
1711 	    {
1712 		m = 0.0;
1713 		z = 0.0;
1714 		if (line->DimensionModel == GAIA_XY_Z)
1715 		  {
1716 		      gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
1717 		  }
1718 		else if (line->DimensionModel == GAIA_XY_M)
1719 		  {
1720 		      gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
1721 		  }
1722 		else if (line->DimensionModel == GAIA_XY_Z_M)
1723 		  {
1724 		      gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
1725 		  }
1726 		else
1727 		  {
1728 		      gaiaGetPoint (line->Coords, iv, &x, &y);
1729 		  }
1730 		x += shift_x;
1731 		y += shift_y;
1732 		if (line->DimensionModel == GAIA_XY_Z)
1733 		  {
1734 		      gaiaSetPointXYZ (line->Coords, iv, x, y, z);
1735 		  }
1736 		else if (line->DimensionModel == GAIA_XY_M)
1737 		  {
1738 		      gaiaSetPointXYM (line->Coords, iv, x, y, m);
1739 		  }
1740 		else if (line->DimensionModel == GAIA_XY_Z_M)
1741 		  {
1742 		      gaiaSetPointXYZM (line->Coords, iv, x, y, z, m);
1743 		  }
1744 		else
1745 		  {
1746 		      gaiaSetPoint (line->Coords, iv, x, y);
1747 		  }
1748 	    }
1749 	  line = line->Next;
1750       }
1751     polyg = geom->FirstPolygon;
1752     while (polyg)
1753       {
1754 	  /* shifting POLYGONs */
1755 	  m = 0.0;
1756 	  z = 0.0;
1757 	  ring = polyg->Exterior;
1758 	  for (iv = 0; iv < ring->Points; iv++)
1759 	    {
1760 		/* shifting the EXTERIOR RING */
1761 		if (ring->DimensionModel == GAIA_XY_Z)
1762 		  {
1763 		      gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
1764 		  }
1765 		else if (ring->DimensionModel == GAIA_XY_M)
1766 		  {
1767 		      gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
1768 		  }
1769 		else if (ring->DimensionModel == GAIA_XY_Z_M)
1770 		  {
1771 		      gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
1772 		  }
1773 		else
1774 		  {
1775 		      gaiaGetPoint (ring->Coords, iv, &x, &y);
1776 		  }
1777 		x += shift_x;
1778 		y += shift_y;
1779 		if (ring->DimensionModel == GAIA_XY_Z)
1780 		  {
1781 		      gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
1782 		  }
1783 		else if (ring->DimensionModel == GAIA_XY_M)
1784 		  {
1785 		      gaiaSetPointXYM (ring->Coords, iv, x, y, m);
1786 		  }
1787 		else if (ring->DimensionModel == GAIA_XY_Z_M)
1788 		  {
1789 		      gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
1790 		  }
1791 		else
1792 		  {
1793 		      gaiaSetPoint (ring->Coords, iv, x, y);
1794 		  }
1795 	    }
1796 	  for (ib = 0; ib < polyg->NumInteriors; ib++)
1797 	    {
1798 		/* shifting the INTERIOR RINGs */
1799 		ring = polyg->Interiors + ib;
1800 		for (iv = 0; iv < ring->Points; iv++)
1801 		  {
1802 		      if (ring->DimensionModel == GAIA_XY_Z)
1803 			{
1804 			    gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
1805 			}
1806 		      else if (ring->DimensionModel == GAIA_XY_M)
1807 			{
1808 			    gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
1809 			}
1810 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
1811 			{
1812 			    gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
1813 			}
1814 		      else
1815 			{
1816 			    gaiaGetPoint (ring->Coords, iv, &x, &y);
1817 			}
1818 		      x += shift_x;
1819 		      y += shift_y;
1820 		      if (ring->DimensionModel == GAIA_XY_Z)
1821 			{
1822 			    gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
1823 			}
1824 		      else if (ring->DimensionModel == GAIA_XY_M)
1825 			{
1826 			    gaiaSetPointXYM (ring->Coords, iv, x, y, m);
1827 			}
1828 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
1829 			{
1830 			    gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
1831 			}
1832 		      else
1833 			{
1834 			    gaiaSetPoint (ring->Coords, iv, x, y);
1835 			}
1836 		  }
1837 	    }
1838 	  polyg = polyg->Next;
1839       }
1840     gaiaMbrGeometry (geom);
1841 }
1842 
1843 GAIAGEO_DECLARE void
gaiaShiftCoords3D(gaiaGeomCollPtr geom,double shift_x,double shift_y,double shift_z)1844 gaiaShiftCoords3D (gaiaGeomCollPtr geom, double shift_x, double shift_y,
1845 		   double shift_z)
1846 {
1847 /* returns a geometry that is the old geometry with required shifting applied to coordinates */
1848     int ib;
1849     int iv;
1850     double x;
1851     double y;
1852     double z;
1853     double m;
1854     gaiaPointPtr point;
1855     gaiaPolygonPtr polyg;
1856     gaiaLinestringPtr line;
1857     gaiaRingPtr ring;
1858     if (!geom)
1859 	return;
1860     point = geom->FirstPoint;
1861     while (point)
1862       {
1863 	  /* shifting POINTs */
1864 	  point->X += shift_x;
1865 	  point->Y += shift_y;
1866 	  if (point->DimensionModel == GAIA_XY_Z
1867 	      || point->DimensionModel == GAIA_XY_Z_M)
1868 	      point->Z += shift_z;
1869 	  point = point->Next;
1870       }
1871     line = geom->FirstLinestring;
1872     while (line)
1873       {
1874 	  /* shifting LINESTRINGs */
1875 	  for (iv = 0; iv < line->Points; iv++)
1876 	    {
1877 		m = 0.0;
1878 		z = 0.0;
1879 		if (line->DimensionModel == GAIA_XY_Z)
1880 		  {
1881 		      gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
1882 		  }
1883 		else if (line->DimensionModel == GAIA_XY_M)
1884 		  {
1885 		      gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
1886 		  }
1887 		else if (line->DimensionModel == GAIA_XY_Z_M)
1888 		  {
1889 		      gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
1890 		  }
1891 		else
1892 		  {
1893 		      gaiaGetPoint (line->Coords, iv, &x, &y);
1894 		  }
1895 		x += shift_x;
1896 		y += shift_y;
1897 		z += shift_z;
1898 		if (line->DimensionModel == GAIA_XY_Z)
1899 		  {
1900 		      gaiaSetPointXYZ (line->Coords, iv, x, y, z);
1901 		  }
1902 		else if (line->DimensionModel == GAIA_XY_M)
1903 		  {
1904 		      gaiaSetPointXYM (line->Coords, iv, x, y, m);
1905 		  }
1906 		else if (line->DimensionModel == GAIA_XY_Z_M)
1907 		  {
1908 		      gaiaSetPointXYZM (line->Coords, iv, x, y, z, m);
1909 		  }
1910 		else
1911 		  {
1912 		      gaiaSetPoint (line->Coords, iv, x, y);
1913 		  }
1914 	    }
1915 	  line = line->Next;
1916       }
1917     polyg = geom->FirstPolygon;
1918     while (polyg)
1919       {
1920 	  /* shifting POLYGONs */
1921 	  m = 0.0;
1922 	  z = 0.0;
1923 	  ring = polyg->Exterior;
1924 	  for (iv = 0; iv < ring->Points; iv++)
1925 	    {
1926 		/* shifting the EXTERIOR RING */
1927 		if (ring->DimensionModel == GAIA_XY_Z)
1928 		  {
1929 		      gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
1930 		  }
1931 		else if (ring->DimensionModel == GAIA_XY_M)
1932 		  {
1933 		      gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
1934 		  }
1935 		else if (ring->DimensionModel == GAIA_XY_Z_M)
1936 		  {
1937 		      gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
1938 		  }
1939 		else
1940 		  {
1941 		      gaiaGetPoint (ring->Coords, iv, &x, &y);
1942 		  }
1943 		x += shift_x;
1944 		y += shift_y;
1945 		z += shift_z;
1946 		if (ring->DimensionModel == GAIA_XY_Z)
1947 		  {
1948 		      gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
1949 		  }
1950 		else if (ring->DimensionModel == GAIA_XY_M)
1951 		  {
1952 		      gaiaSetPointXYM (ring->Coords, iv, x, y, m);
1953 		  }
1954 		else if (ring->DimensionModel == GAIA_XY_Z_M)
1955 		  {
1956 		      gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
1957 		  }
1958 		else
1959 		  {
1960 		      gaiaSetPoint (ring->Coords, iv, x, y);
1961 		  }
1962 	    }
1963 	  for (ib = 0; ib < polyg->NumInteriors; ib++)
1964 	    {
1965 		/* shifting the INTERIOR RINGs */
1966 		ring = polyg->Interiors + ib;
1967 		for (iv = 0; iv < ring->Points; iv++)
1968 		  {
1969 		      if (ring->DimensionModel == GAIA_XY_Z)
1970 			{
1971 			    gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
1972 			}
1973 		      else if (ring->DimensionModel == GAIA_XY_M)
1974 			{
1975 			    gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
1976 			}
1977 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
1978 			{
1979 			    gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
1980 			}
1981 		      else
1982 			{
1983 			    gaiaGetPoint (ring->Coords, iv, &x, &y);
1984 			}
1985 		      x += shift_x;
1986 		      y += shift_y;
1987 		      z += shift_z;
1988 		      if (ring->DimensionModel == GAIA_XY_Z)
1989 			{
1990 			    gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
1991 			}
1992 		      else if (ring->DimensionModel == GAIA_XY_M)
1993 			{
1994 			    gaiaSetPointXYM (ring->Coords, iv, x, y, m);
1995 			}
1996 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
1997 			{
1998 			    gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
1999 			}
2000 		      else
2001 			{
2002 			    gaiaSetPoint (ring->Coords, iv, x, y);
2003 			}
2004 		  }
2005 	    }
2006 	  polyg = polyg->Next;
2007       }
2008     gaiaMbrGeometry (geom);
2009 }
2010 
2011 static void
normalizePoint(double * x,double * y)2012 normalizePoint (double *x, double *y)
2013 {
2014     if ((-180.0 <= *x) && (*x <= 180.0) && (-90.0 <= *y) && (*y <= 90.0))
2015       {
2016 	  /* then this point is already OK */
2017 	  return;
2018       }
2019     if ((*x > 180.0) || (*x < -180.0))
2020       {
2021 	  int numCycles = (int) (*x / 360.0);
2022 	  *x -= numCycles * 360.0;
2023       }
2024     if (*x > 180.0)
2025       {
2026 	  *x -= 360.0;
2027       }
2028     if (*x < -180.0)
2029       {
2030 	  *x += 360.0;
2031       }
2032     if ((*y > 90.0) || (*y < -90.0))
2033       {
2034 	  int numCycles = (int) (*y / 360.0);
2035 	  *y -= numCycles * 360.0;
2036       }
2037     if (*y > 180.0)
2038       {
2039 	  *y = -1.0 * (*y - 180.0);
2040       }
2041     if (*y < -180.0)
2042       {
2043 	  *y = -1.0 * (*y + 180.0);
2044       }
2045     if (*y > 90.0)
2046       {
2047 	  *y = 180 - *y;
2048       }
2049     if (*y < -90.0)
2050       {
2051 	  *y = -180.0 - *y;
2052       }
2053 }
2054 
2055 GAIAGEO_DECLARE void
gaiaNormalizeLonLat(gaiaGeomCollPtr geom)2056 gaiaNormalizeLonLat (gaiaGeomCollPtr geom)
2057 {
2058 /* returns a geometry that is the old geometry with all latitudes shifted
2059  into the range -90 to 90, and all longitudes shifted into the range -180 to
2060  180.
2061 */
2062     int ib;
2063     int iv;
2064     double x = 0.0;
2065     double y = 0.0;
2066     double z = 0.0;
2067     double m = 0.0;
2068     gaiaPointPtr point;
2069     gaiaPolygonPtr polyg;
2070     gaiaLinestringPtr line;
2071     gaiaRingPtr ring;
2072     if (!geom)
2073 	return;
2074     point = geom->FirstPoint;
2075     while (point)
2076       {
2077 	  normalizePoint (&(point->X), &(point->Y));
2078 	  point = point->Next;
2079       }
2080     line = geom->FirstLinestring;
2081     while (line)
2082       {
2083 	  /* shifting LINESTRINGs */
2084 	  for (iv = 0; iv < line->Points; iv++)
2085 	    {
2086 		if (line->DimensionModel == GAIA_XY_Z)
2087 		  {
2088 		      gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
2089 		  }
2090 		else if (line->DimensionModel == GAIA_XY_M)
2091 		  {
2092 		      gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
2093 		  }
2094 		else if (line->DimensionModel == GAIA_XY_Z_M)
2095 		  {
2096 		      gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
2097 		  }
2098 		else
2099 		  {
2100 		      gaiaGetPoint (line->Coords, iv, &x, &y);
2101 		  }
2102 		normalizePoint (&x, &y);
2103 		if (line->DimensionModel == GAIA_XY_Z)
2104 		  {
2105 		      gaiaSetPointXYZ (line->Coords, iv, x, y, z);
2106 		  }
2107 		else if (line->DimensionModel == GAIA_XY_M)
2108 		  {
2109 		      gaiaSetPointXYM (line->Coords, iv, x, y, m);
2110 		  }
2111 		else if (line->DimensionModel == GAIA_XY_Z_M)
2112 		  {
2113 		      gaiaSetPointXYZM (line->Coords, iv, x, y, z, m);
2114 		  }
2115 		else
2116 		  {
2117 		      gaiaSetPoint (line->Coords, iv, x, y);
2118 		  }
2119 	    }
2120 	  line = line->Next;
2121       }
2122     polyg = geom->FirstPolygon;
2123     while (polyg)
2124       {
2125 	  /* shifting POLYGONs */
2126 	  ring = polyg->Exterior;
2127 	  for (iv = 0; iv < ring->Points; iv++)
2128 	    {
2129 		/* shifting the EXTERIOR RING */
2130 		if (ring->DimensionModel == GAIA_XY_Z)
2131 		  {
2132 		      gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2133 		  }
2134 		else if (ring->DimensionModel == GAIA_XY_M)
2135 		  {
2136 		      gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2137 		  }
2138 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2139 		  {
2140 		      gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2141 		  }
2142 		else
2143 		  {
2144 		      gaiaGetPoint (ring->Coords, iv, &x, &y);
2145 		  }
2146 		normalizePoint (&x, &y);
2147 		if (ring->DimensionModel == GAIA_XY_Z)
2148 		  {
2149 		      gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
2150 		  }
2151 		else if (ring->DimensionModel == GAIA_XY_M)
2152 		  {
2153 		      gaiaSetPointXYM (ring->Coords, iv, x, y, m);
2154 		  }
2155 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2156 		  {
2157 		      gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
2158 		  }
2159 		else
2160 		  {
2161 		      gaiaSetPoint (ring->Coords, iv, x, y);
2162 		  }
2163 	    }
2164 	  for (ib = 0; ib < polyg->NumInteriors; ib++)
2165 	    {
2166 		/* shifting the INTERIOR RINGs */
2167 		ring = polyg->Interiors + ib;
2168 		for (iv = 0; iv < ring->Points; iv++)
2169 		  {
2170 		      if (ring->DimensionModel == GAIA_XY_Z)
2171 			{
2172 			    gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2173 			}
2174 		      else if (ring->DimensionModel == GAIA_XY_M)
2175 			{
2176 			    gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2177 			}
2178 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2179 			{
2180 			    gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2181 			}
2182 		      else
2183 			{
2184 			    gaiaGetPoint (ring->Coords, iv, &x, &y);
2185 			}
2186 		      normalizePoint (&x, &y);
2187 		      if (ring->DimensionModel == GAIA_XY_Z)
2188 			{
2189 			    gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
2190 			}
2191 		      else if (ring->DimensionModel == GAIA_XY_M)
2192 			{
2193 			    gaiaSetPointXYM (ring->Coords, iv, x, y, m);
2194 			}
2195 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2196 			{
2197 			    gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
2198 			}
2199 		      else
2200 			{
2201 			    gaiaSetPoint (ring->Coords, iv, x, y);
2202 			}
2203 		  }
2204 	    }
2205 	  polyg = polyg->Next;
2206       }
2207     gaiaMbrGeometry (geom);
2208 }
2209 
2210 GAIAGEO_DECLARE void
gaiaScaleCoords(gaiaGeomCollPtr geom,double scale_x,double scale_y)2211 gaiaScaleCoords (gaiaGeomCollPtr geom, double scale_x, double scale_y)
2212 {
2213 /* returns a geometry that is the old geometry with required scaling applied to coordinates */
2214     int ib;
2215     int iv;
2216     double x;
2217     double y;
2218     double z;
2219     double m;
2220     gaiaPointPtr point;
2221     gaiaPolygonPtr polyg;
2222     gaiaLinestringPtr line;
2223     gaiaRingPtr ring;
2224     if (!geom)
2225 	return;
2226     point = geom->FirstPoint;
2227     while (point)
2228       {
2229 	  /* scaling POINTs */
2230 	  point->X *= scale_x;
2231 	  point->Y *= scale_y;
2232 	  point = point->Next;
2233       }
2234     line = geom->FirstLinestring;
2235     while (line)
2236       {
2237 	  /* scaling LINESTRINGs */
2238 	  for (iv = 0; iv < line->Points; iv++)
2239 	    {
2240 		m = 0.0;
2241 		z = 0.0;
2242 		if (line->DimensionModel == GAIA_XY_Z)
2243 		  {
2244 		      gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
2245 		  }
2246 		else if (line->DimensionModel == GAIA_XY_M)
2247 		  {
2248 		      gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
2249 		  }
2250 		else if (line->DimensionModel == GAIA_XY_Z_M)
2251 		  {
2252 		      gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
2253 		  }
2254 		else
2255 		  {
2256 		      gaiaGetPoint (line->Coords, iv, &x, &y);
2257 		  }
2258 		x *= scale_x;
2259 		y *= scale_y;
2260 		if (line->DimensionModel == GAIA_XY_Z)
2261 		  {
2262 		      gaiaSetPointXYZ (line->Coords, iv, x, y, z);
2263 		  }
2264 		else if (line->DimensionModel == GAIA_XY_M)
2265 		  {
2266 		      gaiaSetPointXYM (line->Coords, iv, x, y, m);
2267 		  }
2268 		else if (line->DimensionModel == GAIA_XY_Z_M)
2269 		  {
2270 		      gaiaSetPointXYZM (line->Coords, iv, x, y, z, m);
2271 		  }
2272 		else
2273 		  {
2274 		      gaiaSetPoint (line->Coords, iv, x, y);
2275 		  }
2276 	    }
2277 	  line = line->Next;
2278       }
2279     polyg = geom->FirstPolygon;
2280     while (polyg)
2281       {
2282 	  /* scaling POLYGONs */
2283 	  m = 0.0;
2284 	  z = 0.0;
2285 	  ring = polyg->Exterior;
2286 	  for (iv = 0; iv < ring->Points; iv++)
2287 	    {
2288 		/* scaling the EXTERIOR RING */
2289 		if (ring->DimensionModel == GAIA_XY_Z)
2290 		  {
2291 		      gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2292 		  }
2293 		else if (ring->DimensionModel == GAIA_XY_M)
2294 		  {
2295 		      gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2296 		  }
2297 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2298 		  {
2299 		      gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2300 		  }
2301 		else
2302 		  {
2303 		      gaiaGetPoint (ring->Coords, iv, &x, &y);
2304 		  }
2305 		x *= scale_x;
2306 		y *= scale_y;
2307 		if (ring->DimensionModel == GAIA_XY_Z)
2308 		  {
2309 		      gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
2310 		  }
2311 		else if (ring->DimensionModel == GAIA_XY_M)
2312 		  {
2313 		      gaiaSetPointXYM (ring->Coords, iv, x, y, m);
2314 		  }
2315 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2316 		  {
2317 		      gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
2318 		  }
2319 		else
2320 		  {
2321 		      gaiaSetPoint (ring->Coords, iv, x, y);
2322 		  }
2323 	    }
2324 	  for (ib = 0; ib < polyg->NumInteriors; ib++)
2325 	    {
2326 		/* scaling the INTERIOR RINGs */
2327 		ring = polyg->Interiors + ib;
2328 		for (iv = 0; iv < ring->Points; iv++)
2329 		  {
2330 		      if (ring->DimensionModel == GAIA_XY_Z)
2331 			{
2332 			    gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2333 			}
2334 		      else if (ring->DimensionModel == GAIA_XY_M)
2335 			{
2336 			    gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2337 			}
2338 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2339 			{
2340 			    gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2341 			}
2342 		      else
2343 			{
2344 			    gaiaGetPoint (ring->Coords, iv, &x, &y);
2345 			}
2346 		      x *= scale_x;
2347 		      y *= scale_y;
2348 		      if (ring->DimensionModel == GAIA_XY_Z)
2349 			{
2350 			    gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
2351 			}
2352 		      else if (ring->DimensionModel == GAIA_XY_M)
2353 			{
2354 			    gaiaSetPointXYM (ring->Coords, iv, x, y, m);
2355 			}
2356 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2357 			{
2358 			    gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
2359 			}
2360 		      else
2361 			{
2362 			    gaiaSetPoint (ring->Coords, iv, x, y);
2363 			}
2364 		  }
2365 	    }
2366 	  polyg = polyg->Next;
2367       }
2368     gaiaMbrGeometry (geom);
2369 }
2370 
2371 GAIAGEO_DECLARE void
gaiaRotateCoords(gaiaGeomCollPtr geom,double angle)2372 gaiaRotateCoords (gaiaGeomCollPtr geom, double angle)
2373 {
2374 /* returns a geometry that is the old geometry with required rotation applied to coordinates */
2375     int ib;
2376     int iv;
2377     double x;
2378     double y;
2379     double z;
2380     double m;
2381     double nx;
2382     double ny;
2383     double rad = angle * 0.0174532925199432958;
2384     double cosine = cos (rad);
2385     double sine = sin (rad);
2386     gaiaPointPtr point;
2387     gaiaPolygonPtr polyg;
2388     gaiaLinestringPtr line;
2389     gaiaRingPtr ring;
2390     if (!geom)
2391 	return;
2392     point = geom->FirstPoint;
2393     while (point)
2394       {
2395 	  /* shifting POINTs */
2396 	  x = point->X;
2397 	  y = point->Y;
2398 	  point->X = (x * cosine) + (y * sine);
2399 	  point->Y = (y * cosine) - (x * sine);
2400 	  point = point->Next;
2401       }
2402     line = geom->FirstLinestring;
2403     while (line)
2404       {
2405 	  /* rotating LINESTRINGs */
2406 	  for (iv = 0; iv < line->Points; iv++)
2407 	    {
2408 		m = 0.0;
2409 		z = 0.0;
2410 		if (line->DimensionModel == GAIA_XY_Z)
2411 		  {
2412 		      gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
2413 		  }
2414 		else if (line->DimensionModel == GAIA_XY_M)
2415 		  {
2416 		      gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
2417 		  }
2418 		else if (line->DimensionModel == GAIA_XY_Z_M)
2419 		  {
2420 		      gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
2421 		  }
2422 		else
2423 		  {
2424 		      gaiaGetPoint (line->Coords, iv, &x, &y);
2425 		  }
2426 		nx = (x * cosine) + (y * sine);
2427 		ny = (y * cosine) - (x * sine);
2428 		if (line->DimensionModel == GAIA_XY_Z)
2429 		  {
2430 		      gaiaSetPointXYZ (line->Coords, iv, nx, ny, z);
2431 		  }
2432 		else if (line->DimensionModel == GAIA_XY_M)
2433 		  {
2434 		      gaiaSetPointXYM (line->Coords, iv, nx, ny, m);
2435 		  }
2436 		else if (line->DimensionModel == GAIA_XY_Z_M)
2437 		  {
2438 		      gaiaSetPointXYZM (line->Coords, iv, nx, ny, z, m);
2439 		  }
2440 		else
2441 		  {
2442 		      gaiaSetPoint (line->Coords, iv, nx, ny);
2443 		  }
2444 	    }
2445 	  line = line->Next;
2446       }
2447     polyg = geom->FirstPolygon;
2448     while (polyg)
2449       {
2450 	  /* rotating POLYGONs */
2451 	  m = 0.0;
2452 	  z = 0.0;
2453 	  ring = polyg->Exterior;
2454 	  for (iv = 0; iv < ring->Points; iv++)
2455 	    {
2456 		/* rotating the EXTERIOR RING */
2457 		if (ring->DimensionModel == GAIA_XY_Z)
2458 		  {
2459 		      gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2460 		  }
2461 		else if (ring->DimensionModel == GAIA_XY_M)
2462 		  {
2463 		      gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2464 		  }
2465 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2466 		  {
2467 		      gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2468 		  }
2469 		else
2470 		  {
2471 		      gaiaGetPoint (ring->Coords, iv, &x, &y);
2472 		  }
2473 		nx = (x * cosine) + (y * sine);
2474 		ny = (y * cosine) - (x * sine);
2475 		if (ring->DimensionModel == GAIA_XY_Z)
2476 		  {
2477 		      gaiaSetPointXYZ (ring->Coords, iv, nx, ny, z);
2478 		  }
2479 		else if (ring->DimensionModel == GAIA_XY_M)
2480 		  {
2481 		      gaiaSetPointXYM (ring->Coords, iv, nx, ny, m);
2482 		  }
2483 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2484 		  {
2485 		      gaiaSetPointXYZM (ring->Coords, iv, nx, ny, z, m);
2486 		  }
2487 		else
2488 		  {
2489 		      gaiaSetPoint (ring->Coords, iv, nx, ny);
2490 		  }
2491 	    }
2492 	  for (ib = 0; ib < polyg->NumInteriors; ib++)
2493 	    {
2494 		/* rotating the INTERIOR RINGs */
2495 		m = 0.0;
2496 		z = 0.0;
2497 		ring = polyg->Interiors + ib;
2498 		for (iv = 0; iv < ring->Points; iv++)
2499 		  {
2500 		      if (ring->DimensionModel == GAIA_XY_Z)
2501 			{
2502 			    gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2503 			}
2504 		      else if (ring->DimensionModel == GAIA_XY_M)
2505 			{
2506 			    gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2507 			}
2508 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2509 			{
2510 			    gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2511 			}
2512 		      else
2513 			{
2514 			    gaiaGetPoint (ring->Coords, iv, &x, &y);
2515 			}
2516 		      nx = (x * cosine) + (y * sine);
2517 		      ny = (y * cosine) - (x * sine);
2518 		      if (ring->DimensionModel == GAIA_XY_Z)
2519 			{
2520 			    gaiaSetPointXYZ (ring->Coords, iv, nx, ny, z);
2521 			}
2522 		      else if (ring->DimensionModel == GAIA_XY_M)
2523 			{
2524 			    gaiaSetPointXYM (ring->Coords, iv, nx, ny, m);
2525 			}
2526 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2527 			{
2528 			    gaiaSetPointXYZM (ring->Coords, iv, nx, ny, z, m);
2529 			}
2530 		      else
2531 			{
2532 			    gaiaSetPoint (ring->Coords, iv, nx, ny);
2533 			}
2534 		  }
2535 	    }
2536 	  polyg = polyg->Next;
2537       }
2538     gaiaMbrGeometry (geom);
2539 }
2540 
2541 GAIAGEO_DECLARE void
gaiaReflectCoords(gaiaGeomCollPtr geom,int x_axis,int y_axis)2542 gaiaReflectCoords (gaiaGeomCollPtr geom, int x_axis, int y_axis)
2543 {
2544 /* returns a geometry that is the old geometry with required reflection applied to coordinates */
2545     int ib;
2546     int iv;
2547     double x;
2548     double y;
2549     double z = 0.0;
2550     double m = 0.0;
2551     gaiaPointPtr point;
2552     gaiaPolygonPtr polyg;
2553     gaiaLinestringPtr line;
2554     gaiaRingPtr ring;
2555     if (!geom)
2556 	return;
2557     point = geom->FirstPoint;
2558     while (point)
2559       {
2560 	  /* reflecting POINTs */
2561 	  if (x_axis)
2562 	      point->X *= -1.0;
2563 	  if (y_axis)
2564 	      point->Y *= -1.0;
2565 	  point = point->Next;
2566       }
2567     line = geom->FirstLinestring;
2568     while (line)
2569       {
2570 	  /* reflecting LINESTRINGs */
2571 	  for (iv = 0; iv < line->Points; iv++)
2572 	    {
2573 		if (line->DimensionModel == GAIA_XY_Z)
2574 		  {
2575 		      gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
2576 		  }
2577 		else if (line->DimensionModel == GAIA_XY_M)
2578 		  {
2579 		      gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
2580 		  }
2581 		else if (line->DimensionModel == GAIA_XY_Z_M)
2582 		  {
2583 		      gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
2584 		  }
2585 		else
2586 		  {
2587 		      gaiaGetPoint (line->Coords, iv, &x, &y);
2588 		  }
2589 		if (x_axis)
2590 		    x *= -1.0;
2591 		if (y_axis)
2592 		    y *= -1.0;
2593 		if (line->DimensionModel == GAIA_XY_Z)
2594 		  {
2595 		      gaiaSetPointXYZ (line->Coords, iv, x, y, z);
2596 		  }
2597 		else if (line->DimensionModel == GAIA_XY_M)
2598 		  {
2599 		      gaiaSetPointXYM (line->Coords, iv, x, y, m);
2600 		  }
2601 		else if (line->DimensionModel == GAIA_XY_Z_M)
2602 		  {
2603 		      gaiaSetPointXYZM (line->Coords, iv, x, y, z, m);
2604 		  }
2605 		else
2606 		  {
2607 		      gaiaSetPoint (line->Coords, iv, x, y);
2608 		  }
2609 	    }
2610 	  line = line->Next;
2611       }
2612     polyg = geom->FirstPolygon;
2613     while (polyg)
2614       {
2615 	  /* reflecting POLYGONs */
2616 	  ring = polyg->Exterior;
2617 	  for (iv = 0; iv < ring->Points; iv++)
2618 	    {
2619 		/* reflecting the EXTERIOR RING */
2620 		if (ring->DimensionModel == GAIA_XY_Z)
2621 		  {
2622 		      gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2623 		  }
2624 		else if (ring->DimensionModel == GAIA_XY_M)
2625 		  {
2626 		      gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2627 		  }
2628 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2629 		  {
2630 		      gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2631 		  }
2632 		else
2633 		  {
2634 		      gaiaGetPoint (ring->Coords, iv, &x, &y);
2635 		  }
2636 		if (x_axis)
2637 		    x *= -1.0;
2638 		if (y_axis)
2639 		    y *= -1.0;
2640 		if (ring->DimensionModel == GAIA_XY_Z)
2641 		  {
2642 		      gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
2643 		  }
2644 		else if (ring->DimensionModel == GAIA_XY_M)
2645 		  {
2646 		      gaiaSetPointXYM (ring->Coords, iv, x, y, m);
2647 		  }
2648 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2649 		  {
2650 		      gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
2651 		  }
2652 		else
2653 		  {
2654 		      gaiaSetPoint (ring->Coords, iv, x, y);
2655 		  }
2656 	    }
2657 	  for (ib = 0; ib < polyg->NumInteriors; ib++)
2658 	    {
2659 		/* reflecting the INTERIOR RINGs */
2660 		ring = polyg->Interiors + ib;
2661 		for (iv = 0; iv < ring->Points; iv++)
2662 		  {
2663 		      if (ring->DimensionModel == GAIA_XY_Z)
2664 			{
2665 			    gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2666 			}
2667 		      else if (ring->DimensionModel == GAIA_XY_M)
2668 			{
2669 			    gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2670 			}
2671 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2672 			{
2673 			    gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2674 			}
2675 		      else
2676 			{
2677 			    gaiaGetPoint (ring->Coords, iv, &x, &y);
2678 			}
2679 		      if (x_axis)
2680 			  x *= -1.0;
2681 		      if (y_axis)
2682 			  y *= -1.0;
2683 		      if (ring->DimensionModel == GAIA_XY_Z)
2684 			{
2685 			    gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
2686 			}
2687 		      else if (ring->DimensionModel == GAIA_XY_M)
2688 			{
2689 			    gaiaSetPointXYM (ring->Coords, iv, x, y, m);
2690 			}
2691 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2692 			{
2693 			    gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
2694 			}
2695 		      else
2696 			{
2697 			    gaiaSetPoint (ring->Coords, iv, x, y);
2698 			}
2699 		  }
2700 	    }
2701 	  polyg = polyg->Next;
2702       }
2703     gaiaMbrGeometry (geom);
2704 }
2705 
2706 GAIAGEO_DECLARE void
gaiaSwapCoords(gaiaGeomCollPtr geom)2707 gaiaSwapCoords (gaiaGeomCollPtr geom)
2708 {
2709 /* returns a geometry that is the old geometry with swapped x- and y-coordinates */
2710     int ib;
2711     int iv;
2712     double x;
2713     double y;
2714     double z;
2715     double m;
2716     double sv;
2717     gaiaPointPtr point;
2718     gaiaPolygonPtr polyg;
2719     gaiaLinestringPtr line;
2720     gaiaRingPtr ring;
2721     if (!geom)
2722 	return;
2723     point = geom->FirstPoint;
2724     while (point)
2725       {
2726 	  /* swapping POINTs */
2727 	  sv = point->X;
2728 	  point->X = point->Y;
2729 	  point->Y = sv;
2730 	  point = point->Next;
2731       }
2732     line = geom->FirstLinestring;
2733     while (line)
2734       {
2735 	  /* swapping LINESTRINGs */
2736 	  for (iv = 0; iv < line->Points; iv++)
2737 	    {
2738 		m = 0.0;
2739 		z = 0.0;
2740 		if (line->DimensionModel == GAIA_XY_Z)
2741 		  {
2742 		      gaiaGetPointXYZ (line->Coords, iv, &x, &y, &z);
2743 		  }
2744 		else if (line->DimensionModel == GAIA_XY_M)
2745 		  {
2746 		      gaiaGetPointXYM (line->Coords, iv, &x, &y, &m);
2747 		  }
2748 		else if (line->DimensionModel == GAIA_XY_Z_M)
2749 		  {
2750 		      gaiaGetPointXYZM (line->Coords, iv, &x, &y, &z, &m);
2751 		  }
2752 		else
2753 		  {
2754 		      gaiaGetPoint (line->Coords, iv, &x, &y);
2755 		  }
2756 		sv = x;
2757 		x = y;
2758 		y = sv;
2759 		if (line->DimensionModel == GAIA_XY_Z)
2760 		  {
2761 		      gaiaSetPointXYZ (line->Coords, iv, x, y, z);
2762 		  }
2763 		else if (line->DimensionModel == GAIA_XY_M)
2764 		  {
2765 		      gaiaSetPointXYM (line->Coords, iv, x, y, m);
2766 		  }
2767 		else if (line->DimensionModel == GAIA_XY_Z_M)
2768 		  {
2769 		      gaiaSetPointXYZM (line->Coords, iv, x, y, z, m);
2770 		  }
2771 		else
2772 		  {
2773 		      gaiaSetPoint (line->Coords, iv, x, y);
2774 		  }
2775 	    }
2776 	  line = line->Next;
2777       }
2778     polyg = geom->FirstPolygon;
2779     while (polyg)
2780       {
2781 	  /* swapping POLYGONs */
2782 	  m = 0.0;
2783 	  z = 0.0;
2784 	  ring = polyg->Exterior;
2785 	  for (iv = 0; iv < ring->Points; iv++)
2786 	    {
2787 		/* shifting the EXTERIOR RING */
2788 		if (ring->DimensionModel == GAIA_XY_Z)
2789 		  {
2790 		      gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2791 		  }
2792 		else if (ring->DimensionModel == GAIA_XY_M)
2793 		  {
2794 		      gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2795 		  }
2796 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2797 		  {
2798 		      gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2799 		  }
2800 		else
2801 		  {
2802 		      gaiaGetPoint (ring->Coords, iv, &x, &y);
2803 		  }
2804 		sv = x;
2805 		x = y;
2806 		y = sv;
2807 		if (ring->DimensionModel == GAIA_XY_Z)
2808 		  {
2809 		      gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
2810 		  }
2811 		else if (ring->DimensionModel == GAIA_XY_M)
2812 		  {
2813 		      gaiaSetPointXYM (ring->Coords, iv, x, y, m);
2814 		  }
2815 		else if (ring->DimensionModel == GAIA_XY_Z_M)
2816 		  {
2817 		      gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
2818 		  }
2819 		else
2820 		  {
2821 		      gaiaSetPoint (ring->Coords, iv, x, y);
2822 		  }
2823 	    }
2824 	  for (ib = 0; ib < polyg->NumInteriors; ib++)
2825 	    {
2826 		/* swapping the INTERIOR RINGs */
2827 		m = 0.0;
2828 		z = 0.0;
2829 		ring = polyg->Interiors + ib;
2830 		for (iv = 0; iv < ring->Points; iv++)
2831 		  {
2832 		      if (ring->DimensionModel == GAIA_XY_Z)
2833 			{
2834 			    gaiaGetPointXYZ (ring->Coords, iv, &x, &y, &z);
2835 			}
2836 		      else if (ring->DimensionModel == GAIA_XY_M)
2837 			{
2838 			    gaiaGetPointXYM (ring->Coords, iv, &x, &y, &m);
2839 			}
2840 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2841 			{
2842 			    gaiaGetPointXYZM (ring->Coords, iv, &x, &y, &z, &m);
2843 			}
2844 		      else
2845 			{
2846 			    gaiaGetPoint (ring->Coords, iv, &x, &y);
2847 			}
2848 		      sv = x;
2849 		      x = y;
2850 		      y = sv;
2851 		      if (ring->DimensionModel == GAIA_XY_Z)
2852 			{
2853 			    gaiaSetPointXYZ (ring->Coords, iv, x, y, z);
2854 			}
2855 		      else if (ring->DimensionModel == GAIA_XY_M)
2856 			{
2857 			    gaiaSetPointXYM (ring->Coords, iv, x, y, m);
2858 			}
2859 		      else if (ring->DimensionModel == GAIA_XY_Z_M)
2860 			{
2861 			    gaiaSetPointXYZM (ring->Coords, iv, x, y, z, m);
2862 			}
2863 		      else
2864 			{
2865 			    gaiaSetPoint (ring->Coords, iv, x, y);
2866 			}
2867 		  }
2868 	    }
2869 	  polyg = polyg->Next;
2870       }
2871     gaiaMbrGeometry (geom);
2872 }
2873