1 /******************************************************************************
2 * $Id$
3 *
4 * Project: MapServer
5 * Purpose: projectionObj / PROJ.4 interface.
6 * Author: Steve Lime and the MapServer team.
7 *
8 ******************************************************************************
9 * Copyright (c) 1996-2005 Regents of the University of Minnesota.
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included in
19 * all copies of this Software or works derived from this Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 ****************************************************************************/
29
30 #include "mapserver.h"
31 #include "mapproject.h"
32 #include "mapthread.h"
33 #include <assert.h>
34 #include <float.h>
35 #include <sys/types.h>
36 #include <sys/stat.h>
37 #include "mapaxisorder.h"
38
39 #include "ogr_srs_api.h"
40
41 static char *ms_proj_lib = NULL;
42 #if PROJ_VERSION_MAJOR >= 6
43 static unsigned ms_proj_lib_change_counter = 0;
44 #endif
45
46 typedef struct LinkedListOfProjContext LinkedListOfProjContext;
47 struct LinkedListOfProjContext
48 {
49 LinkedListOfProjContext* next;
50 projectionContext* context;
51 };
52
53 static LinkedListOfProjContext* headOfLinkedListOfProjContext = NULL;
54
55 static int msTestNeedWrap( pointObj pt1, pointObj pt2, pointObj pt2_geo,
56 reprojectionObj* reprojector );
57
58
59 static projectionContext* msProjectionContextCreate(void);
60 static void msProjectionContextUnref(projectionContext* ctx);
61
62 #if PROJ_VERSION_MAJOR >= 6
63
64 #include "proj_experimental.h"
65
66 struct reprojectionObj
67 {
68 projectionObj* in;
69 projectionObj* out;
70 PJ* pj;
71 int should_do_line_cutting;
72 shapeObj splitShape;
73 int bFreePJ;
74 };
75
76 /* Helps considerably for use cases like msautotest/wxs/wms_inspire.map */
77 /* which involve a number of layers with same SRS, and a number of exposed */
78 /* output SRS */
79 #define PJ_CACHE_ENTRY_SIZE 32
80
81 typedef struct
82 {
83 char* inStr;
84 char* outStr;
85 PJ* pj;
86 } pjCacheEntry;
87
88 struct projectionContext
89 {
90 PJ_CONTEXT* proj_ctx;
91 unsigned ms_proj_lib_change_counter;
92 int ref_count;
93 pjCacheEntry pj_cache[PJ_CACHE_ENTRY_SIZE];
94 int pj_cache_size;
95 };
96
97 /************************************************************************/
98 /* msProjectHasLonWrapOrOver() */
99 /************************************************************************/
100
msProjectHasLonWrapOrOver(projectionObj * in)101 static int msProjectHasLonWrapOrOver(projectionObj *in) {
102 int i;
103 for( i = 0; i < in->numargs; i++ )
104 {
105 if( strncmp(in->args[i], "lon_wrap=", strlen("lon_wrap=")) == 0 ||
106 strcmp(in->args[i], "+over") == 0 ||
107 strcmp(in->args[i], "over") == 0 )
108 {
109 return MS_TRUE;
110 }
111 }
112 return MS_FALSE;
113 }
114
115 /************************************************************************/
116 /* createNormalizedPJ() */
117 /************************************************************************/
118
119 /* Return to be freed with proj_destroy() if *pbFreePJ = TRUE */
createNormalizedPJ(projectionObj * in,projectionObj * out,int * pbFreePJ)120 static PJ* createNormalizedPJ(projectionObj *in, projectionObj *out, int* pbFreePJ)
121 {
122 if( in->proj == out->proj )
123 {
124 /* Special case to avoid out_str below to cause in_str to become invalid */
125 *pbFreePJ = TRUE;
126 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR == 0
127 /* 6.0 didn't support proj=noop */
128 return proj_create(in->proj_ctx->proj_ctx, "+proj=affine");
129 #else
130 return proj_create(in->proj_ctx->proj_ctx, "+proj=noop");
131 #endif
132 }
133
134 const char* const wkt_options[] = { "MULTILINE=NO", NULL };
135 const char* in_str = msProjectHasLonWrapOrOver(in) ?
136 proj_as_proj_string(in->proj_ctx->proj_ctx, in->proj, PJ_PROJ_4, NULL) :
137 proj_as_wkt(in->proj_ctx->proj_ctx, in->proj, PJ_WKT2_2018, wkt_options);
138 const char* out_str = msProjectHasLonWrapOrOver(out) ?
139 proj_as_proj_string(out->proj_ctx->proj_ctx, out->proj, PJ_PROJ_4, NULL) :
140 proj_as_wkt(out->proj_ctx->proj_ctx, out->proj, PJ_WKT2_2018, wkt_options);
141 PJ* pj_raw;
142 PJ* pj_normalized;
143 if( !in_str || !out_str )
144 return NULL;
145
146 if( in->proj_ctx->proj_ctx == out->proj_ctx->proj_ctx )
147 {
148 int i;
149 pjCacheEntry* pj_cache = in->proj_ctx->pj_cache;
150 for( i = 0; i < in->proj_ctx->pj_cache_size; i++ )
151 {
152 if (strcmp(pj_cache[i].inStr, in_str) == 0 &&
153 strcmp(pj_cache[i].outStr, out_str) == 0 )
154 {
155 PJ* ret = pj_cache[i].pj;
156 if( i != 0 )
157 {
158 /* Move entry to top */
159 pjCacheEntry tmp;
160 memcpy(&tmp, &pj_cache[i], sizeof(pjCacheEntry));
161 memmove(&pj_cache[1], &pj_cache[0], i * sizeof(pjCacheEntry));
162 memcpy(&pj_cache[0], &tmp, sizeof(pjCacheEntry));
163 }
164 #ifdef notdef
165 fprintf(stderr, "cache hit!\n");
166 #endif
167 *pbFreePJ = FALSE;
168 return ret;
169 }
170 }
171 }
172
173 #ifdef notdef
174 fprintf(stderr, "%s -> %s\n", in_str, out_str);
175 fprintf(stderr, "%p -> %p\n", in->proj_ctx->proj_ctx, out->proj_ctx->proj_ctx);
176 fprintf(stderr, "cache miss!\n");
177 #endif
178
179 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 2
180 if( strstr(in_str, "+proj=") && strstr(in_str, "+over") &&
181 strstr(out_str, "+proj=") && strstr(out_str, "+over") &&
182 strlen(in_str) < 400 && strlen(out_str) < 400 )
183 {
184 // Fixed per PROJ commit
185 // https://github.com/OSGeo/PROJ/commit/78302efb70eb4b49610cda6a60bf9ce39b82264f
186 // and
187 // https://github.com/OSGeo/PROJ/commit/ae70b26b9cbae85a38d5b26533ba06da0ea13940
188 // Fix for wcs_get_capabilities_tileindexmixedsrs_26711.xml and wcs_20_getcov_bands_name_new_reproject.dat
189 char szPipeline[1024];
190 strcpy(szPipeline, "+proj=pipeline");
191 if( msProjIsGeographicCRS(in) )
192 {
193 strcat(szPipeline, " +step +proj=unitconvert +xy_in=deg +xy_out=rad");
194 }
195 strcat(szPipeline, " +step +inv ");
196 strcat(szPipeline, in_str);
197 strcat(szPipeline, " +step ");
198 strcat(szPipeline, out_str);
199 if( msProjIsGeographicCRS(out) )
200 {
201 strcat(szPipeline, " +step +proj=unitconvert +xy_in=rad +xy_out=deg");
202 }
203 /* We do not want datum=NAD83 to imply a transformation with towgs84=0,0,0 */
204 {
205 char* ptr = szPipeline;
206 while(1)
207 {
208 ptr = strstr(ptr, " +datum=NAD83");
209 if( !ptr )
210 break;
211 memcpy(ptr, " +ellps=GRS80", 13);
212 }
213 }
214
215 /* Remove +nadgrids=@null as it doesn't work if going outside of [-180,180] */
216 /* Fixed per https://github.com/OSGeo/PROJ/commit/10a30bb539be1afb25952b19af8bbe72e1b13b56 */
217 {
218 char* ptr = strstr(szPipeline, " +nadgrids=@null");
219 if( ptr )
220 memcpy(ptr, " ", 16);
221 }
222
223 pj_normalized = proj_create(in->proj_ctx->proj_ctx, szPipeline);
224 }
225 else
226 #endif
227 {
228 pj_raw = proj_create_crs_to_crs(in->proj_ctx->proj_ctx, in_str, out_str, NULL);
229 if( !pj_raw )
230 return NULL;
231 pj_normalized = proj_normalize_for_visualization(in->proj_ctx->proj_ctx, pj_raw);
232 proj_destroy(pj_raw);
233 }
234 if( !pj_normalized )
235 return NULL;
236
237 if( in->proj_ctx->proj_ctx == out->proj_ctx->proj_ctx )
238 {
239 /* Insert entry into cache */
240 int i;
241 pjCacheEntry* pj_cache = in->proj_ctx->pj_cache;
242 if( in->proj_ctx->pj_cache_size < PJ_CACHE_ENTRY_SIZE )
243 {
244 i = in->proj_ctx->pj_cache_size;
245 assert ( pj_cache[i].inStr == NULL );
246 in->proj_ctx->pj_cache_size ++;
247 }
248 else
249 {
250 i = 0;
251 /* Evict oldest entry */
252 msFree(pj_cache[PJ_CACHE_ENTRY_SIZE - 1].inStr);
253 msFree(pj_cache[PJ_CACHE_ENTRY_SIZE - 1].outStr);
254 proj_destroy(pj_cache[PJ_CACHE_ENTRY_SIZE - 1].pj);
255 memmove(&pj_cache[1], &pj_cache[0],
256 (PJ_CACHE_ENTRY_SIZE - 1) * sizeof(pjCacheEntry));
257 }
258 pj_cache[i].inStr = msStrdup(in_str);
259 pj_cache[i].outStr = msStrdup(out_str);
260 pj_cache[i].pj = pj_normalized;
261 *pbFreePJ = FALSE;
262 }
263 else
264 {
265 *pbFreePJ = TRUE;
266 }
267
268 return pj_normalized;
269 }
270
271 /************************************************************************/
272 /* getBaseGeographicCRS() */
273 /************************************************************************/
274
275 /* Return to be freed with proj_destroy() */
getBaseGeographicCRS(projectionObj * in)276 static PJ* getBaseGeographicCRS(projectionObj* in)
277 {
278 PJ_CONTEXT* ctxt;
279 PJ_TYPE type;
280 assert(in && in->proj);
281 ctxt = in->proj_ctx->proj_ctx;
282 type = proj_get_type(in->proj);
283 if( type == PJ_TYPE_PROJECTED_CRS )
284 {
285 return proj_get_source_crs(ctxt, in->proj);
286 }
287 if( type == PJ_TYPE_BOUND_CRS )
288 {
289 /* If it is a boundCRS of a projectedCRS, extract the geographicCRS
290 * from the projectedCRS, and rewrap it in a boundCRS */
291 PJ* source_crs = proj_get_source_crs(ctxt, in->proj);
292 PJ* geog_source_crs;
293 PJ* hub_crs;
294 PJ* transf;
295 PJ* ret;
296 if( proj_get_type(source_crs) != PJ_TYPE_PROJECTED_CRS )
297 {
298 proj_destroy(source_crs);
299 return NULL;
300 }
301 geog_source_crs = proj_get_source_crs(ctxt, source_crs);
302 proj_destroy(source_crs);
303 hub_crs = proj_get_target_crs(ctxt, in->proj);
304 transf = proj_crs_get_coordoperation(ctxt, in->proj);
305 ret = proj_crs_create_bound_crs(ctxt, geog_source_crs,
306 hub_crs, transf);
307 proj_destroy(geog_source_crs);
308 proj_destroy(hub_crs);
309 proj_destroy(transf);
310 return ret;
311
312 }
313 return NULL;
314 }
315
316 /************************************************************************/
317 /* msProjErrorLogger() */
318 /************************************************************************/
319
msProjErrorLogger(void * user_data,int level,const char * message)320 static void msProjErrorLogger(void * user_data,
321 int level, const char * message)
322 {
323 (void)user_data;
324 if( level == PJ_LOG_ERROR )
325 {
326 msDebug( "PROJ: Error: %s\n", message );
327 }
328 else if( level == PJ_LOG_DEBUG )
329 {
330 msDebug( "PROJ: Debug: %s\n", message );
331 }
332 }
333
334 /************************************************************************/
335 /* msProjectionContextCreate() */
336 /************************************************************************/
337
msProjectionContextCreate(void)338 projectionContext* msProjectionContextCreate(void)
339 {
340 projectionContext* ctx = (projectionContext*)msSmallCalloc(1, sizeof(projectionContext));
341 ctx->proj_ctx = proj_context_create();
342 if( ctx->proj_ctx == NULL )
343 {
344 msFree(ctx);
345 return NULL;
346 }
347 ctx->ref_count = 1;
348 proj_context_use_proj4_init_rules(ctx->proj_ctx, TRUE);
349 proj_log_func (ctx->proj_ctx, NULL, msProjErrorLogger);
350 return ctx;
351 }
352
353 /************************************************************************/
354 /* msProjectionContextUnref() */
355 /************************************************************************/
356
msProjectionContextUnref(projectionContext * ctx)357 void msProjectionContextUnref(projectionContext* ctx)
358 {
359 if( !ctx )
360 return;
361 --ctx->ref_count;
362 if( ctx->ref_count == 0 )
363 {
364 int i;
365 for( i = 0; i < ctx->pj_cache_size; i++ )
366 {
367 msFree(ctx->pj_cache[i].inStr);
368 msFree(ctx->pj_cache[i].outStr);
369 proj_destroy(ctx->pj_cache[i].pj);
370 }
371 proj_context_destroy(ctx->proj_ctx);
372 msFree(ctx);
373 }
374 }
375
376 /************************************************************************/
377 /* msProjectCreateReprojector() */
378 /************************************************************************/
379
msProjectCreateReprojector(projectionObj * in,projectionObj * out)380 reprojectionObj* msProjectCreateReprojector(projectionObj* in, projectionObj* out)
381 {
382 reprojectionObj* obj = (reprojectionObj*)msSmallCalloc(1, sizeof(reprojectionObj));
383 obj->in = in;
384 obj->out = out;
385 obj->should_do_line_cutting = -1;
386
387 /* -------------------------------------------------------------------- */
388 /* If the source and destination are simple and equal, then do */
389 /* nothing. */
390 /* -------------------------------------------------------------------- */
391 if( in && in->numargs == 1 && out && out->numargs == 1
392 && strcmp(in->args[0],out->args[0]) == 0 ) {
393 /* do nothing, no transformation required */
394 }
395 /* -------------------------------------------------------------------- */
396 /* If we have a fully defined input coordinate system and */
397 /* output coordinate system, then we will use createNormalizedPJ */
398 /* -------------------------------------------------------------------- */
399 else if( in && in->proj && out && out->proj ) {
400 PJ* pj = createNormalizedPJ(in, out, &(obj->bFreePJ));
401 if( !pj )
402 {
403 msFree(obj);
404 return NULL;
405 }
406 obj->pj = pj;
407 }
408
409 /* nothing to do if the other coordinate system is also lat/long */
410 else if( (in == NULL || in->proj == NULL) &&
411 (out == NULL || out->proj == NULL || msProjIsGeographicCRS(out) ))
412 {
413 /* do nothing */
414 }
415 else if( out == NULL && in != NULL && msProjIsGeographicCRS(in) )
416 {
417 /* do nothing */
418 }
419
420 /* -------------------------------------------------------------------- */
421 /* Otherwise we assume that the NULL projectionObj is supposed to be */
422 /* lat/long in the same datum as the other projectionObj. This */
423 /* is essentially a backwards compatibility mode. */
424 /* -------------------------------------------------------------------- */
425 else {
426 PJ* pj = NULL;
427
428 if( (in==NULL || in->proj==NULL) && out && out->proj) { /* input coordinates are lat/lon */
429 PJ* source_crs = getBaseGeographicCRS(out);
430 projectionObj in_modified;
431 memset(&in_modified, 0, sizeof(in_modified));
432
433 in_modified.proj_ctx = out->proj_ctx;
434 in_modified.proj = source_crs;
435 pj = createNormalizedPJ(&in_modified, out, &(obj->bFreePJ));
436 proj_destroy(source_crs);
437 } else if( /* (out==NULL || out->proj==NULL) && */ in && in->proj ) {
438 PJ* target_crs = getBaseGeographicCRS(in);
439 projectionObj out_modified;
440 memset(&out_modified, 0, sizeof(out_modified));
441
442 out_modified.proj_ctx = in->proj_ctx;
443 out_modified.proj = target_crs;
444 pj = createNormalizedPJ(in, &out_modified, &(obj->bFreePJ));
445 proj_destroy(target_crs);
446 }
447 if( !pj )
448 {
449 msFree(obj);
450 return NULL;
451 }
452 obj->pj = pj;
453 }
454
455 return obj;
456 }
457
458 /************************************************************************/
459 /* msProjectDestroyReprojector() */
460 /************************************************************************/
461
msProjectDestroyReprojector(reprojectionObj * reprojector)462 void msProjectDestroyReprojector(reprojectionObj* reprojector)
463 {
464 if( !reprojector )
465 return;
466 if( reprojector->bFreePJ )
467 proj_destroy(reprojector->pj);
468 msFreeShape(&(reprojector->splitShape));
469 msFree(reprojector);
470 }
471
472 /************************************************************************/
473 /* msProjectTransformPoints() */
474 /************************************************************************/
475
msProjectTransformPoints(reprojectionObj * reprojector,int npoints,double * x,double * y)476 int msProjectTransformPoints( reprojectionObj* reprojector,
477 int npoints, double* x, double* y )
478 {
479 proj_trans_generic (reprojector->pj, PJ_FWD,
480 x, sizeof(double), npoints,
481 y, sizeof(double), npoints,
482 NULL, 0, 0,
483 NULL, 0, 0 );
484 return MS_SUCCESS;
485 }
486
487 #else
488
489 /************************************************************************/
490 /* msProjectionContextCreate() */
491 /************************************************************************/
492
msProjectionContextCreate(void)493 projectionContext* msProjectionContextCreate(void)
494 {
495 return NULL;
496 }
497
498 /************************************************************************/
499 /* msProjectionContextUnref() */
500 /************************************************************************/
501
msProjectionContextUnref(projectionContext * ctx)502 void msProjectionContextUnref(projectionContext* ctx)
503 {
504 (void)ctx;
505 }
506
507 struct reprojectionObj
508 {
509 projectionObj* in;
510 projectionObj* out;
511 int should_do_line_cutting;
512 shapeObj splitShape;
513 int no_op;
514 };
515
516 /************************************************************************/
517 /* msProjectCreateReprojector() */
518 /************************************************************************/
519
msProjectCreateReprojector(projectionObj * in,projectionObj * out)520 reprojectionObj* msProjectCreateReprojector(projectionObj* in, projectionObj* out)
521 {
522 reprojectionObj* obj;
523 obj = (reprojectionObj*)msSmallCalloc(1, sizeof(reprojectionObj));
524 obj->in = in;
525 obj->out = out;
526 obj->should_do_line_cutting = -1;
527
528 /* -------------------------------------------------------------------- */
529 /* If the source and destination are equal, then do nothing. */
530 /* -------------------------------------------------------------------- */
531 if( in && out && in->numargs == out->numargs && in->numargs > 0
532 && strcmp(in->args[0],out->args[0]) == 0 )
533 {
534 int i;
535 obj->no_op = MS_TRUE;
536 for( i = 1; i < in->numargs; i++ )
537 {
538 if( strcmp(in->args[i],out->args[i]) != 0 ) {
539 obj->no_op = MS_FALSE;
540 break;
541 }
542 }
543 }
544
545 /* -------------------------------------------------------------------- */
546 /* If we have a fully defined input coordinate system and */
547 /* output coordinate system, then we will use pj_transform */
548 /* -------------------------------------------------------------------- */
549 else if( in && in->proj && out && out->proj ) {
550 /* do nothing for now */
551 }
552 /* nothing to do if the other coordinate system is also lat/long */
553 else if( in == NULL && (out == NULL || msProjIsGeographicCRS(out) ))
554 {
555 obj->no_op = MS_TRUE;
556 }
557 else if( out == NULL && in != NULL && msProjIsGeographicCRS(in) )
558 {
559 obj->no_op = MS_TRUE;
560 }
561 return obj;
562 }
563
564 /************************************************************************/
565 /* msProjectDestroyReprojector() */
566 /************************************************************************/
567
msProjectDestroyReprojector(reprojectionObj * reprojector)568 void msProjectDestroyReprojector(reprojectionObj* reprojector)
569 {
570 if( !reprojector )
571 return;
572 msFreeShape(&(reprojector->splitShape));
573 msFree(reprojector);
574 }
575
576 #endif
577
578 /*
579 ** Initialize, load and free a projectionObj structure
580 */
msInitProjection(projectionObj * p)581 int msInitProjection(projectionObj *p)
582 {
583 p->gt.need_geotransform = MS_FALSE;
584 p->numargs = 0;
585 p->args = NULL;
586 p->wellknownprojection = wkp_none;
587 p->proj = NULL;
588 p->args = (char **)malloc(MS_MAXPROJARGS*sizeof(char *));
589 MS_CHECK_ALLOC(p->args, MS_MAXPROJARGS*sizeof(char *), -1);
590 #if PROJ_VERSION_MAJOR >= 6
591 p->proj_ctx = NULL;
592 #elif PJ_VERSION >= 480
593 p->proj_ctx = NULL;
594 #endif
595 return(0);
596 }
597
msFreeProjection(projectionObj * p)598 void msFreeProjection(projectionObj *p)
599 {
600 #if PROJ_VERSION_MAJOR >= 6
601 proj_destroy(p->proj);
602 p->proj = NULL;
603 msProjectionContextUnref(p->proj_ctx);
604 p->proj_ctx = NULL;
605 #else
606 if(p->proj) {
607 pj_free(p->proj);
608 p->proj = NULL;
609 }
610 #if PJ_VERSION >= 480
611 if(p->proj_ctx) {
612 pj_ctx_free(p->proj_ctx);
613 p->proj_ctx = NULL;
614 }
615 #endif
616 #endif
617
618 p->gt.need_geotransform = MS_FALSE;
619 p->wellknownprojection = wkp_none;
620 msFreeCharArray(p->args, p->numargs);
621 p->args = NULL;
622 p->numargs = 0;
623 }
624
msFreeProjectionExceptContext(projectionObj * p)625 void msFreeProjectionExceptContext(projectionObj *p)
626 {
627 #if PROJ_VERSION_MAJOR >= 6
628 projectionContext* ctx = p->proj_ctx;
629 p->proj_ctx = NULL;
630 msFreeProjection(p);
631 p->proj_ctx = ctx;
632 #else
633 msFreeProjection(p);
634 #endif
635 }
636
637 /************************************************************************/
638 /* msProjectionInheritContextFrom() */
639 /************************************************************************/
640
msProjectionInheritContextFrom(projectionObj * pDst,projectionObj * pSrc)641 void msProjectionInheritContextFrom(projectionObj *pDst, projectionObj* pSrc)
642 {
643 #if PROJ_VERSION_MAJOR >= 6
644 if( pDst->proj_ctx == NULL && pSrc->proj_ctx != NULL)
645 {
646 pDst->proj_ctx = pSrc->proj_ctx;
647 pDst->proj_ctx->ref_count ++;
648 }
649 #endif
650 }
651
652 /************************************************************************/
653 /* msProjectionSetContext() */
654 /************************************************************************/
655
msProjectionSetContext(projectionObj * p,projectionContext * ctx)656 void msProjectionSetContext(projectionObj *p, projectionContext* ctx)
657 {
658 #if PROJ_VERSION_MAJOR >= 6
659 if( p->proj_ctx == NULL && ctx != NULL)
660 {
661 p->proj_ctx = ctx;
662 p->proj_ctx->ref_count ++;
663 }
664 #endif
665 }
666
667 /*
668 ** Handle OGC WMS/WFS AUTO projection in the format:
669 ** "AUTO:proj_id,units_id,lon0,lat0"
670 */
671
_msProcessAutoProjection(projectionObj * p)672 static int _msProcessAutoProjection(projectionObj *p)
673 {
674 char **args;
675 int numargs, nProjId, nUnitsId, nZone;
676 double dLat0, dLon0;
677 const char *pszUnits = "m";
678 char szProjBuf[512]="";
679
680 /* WMS/WFS AUTO projection: "AUTO:proj_id,units_id,lon0,lat0" */
681 args = msStringSplit(p->args[0], ',', &numargs);
682 if (numargs != 4 ||
683 (strncasecmp(args[0], "AUTO:", 5) != 0 &&
684 strncasecmp(args[0], "AUTO2:", 6) != 0)) {
685 msSetError(MS_PROJERR,
686 "WMS/WFS AUTO/AUTO2 PROJECTION must be in the format "
687 "'AUTO:proj_id,units_id,lon0,lat0' or 'AUTO2:crs_id,factor,lon0,lat0'(got '%s').\n",
688 "_msProcessAutoProjection()", p->args[0]);
689 return -1;
690 }
691
692 if (strncasecmp(args[0], "AUTO:", 5)==0)
693 nProjId = atoi(args[0]+5);
694 else
695 nProjId = atoi(args[0]+6);
696
697 nUnitsId = atoi(args[1]);
698 dLon0 = atof(args[2]);
699 dLat0 = atof(args[3]);
700
701
702 /*There is no unit parameter for AUTO2. The 2nd parameter is
703 factor. Set the units to always be meter*/
704 if (strncasecmp(args[0], "AUTO2:", 6) == 0)
705 nUnitsId = 9001;
706
707 msFreeCharArray(args, numargs);
708
709 /* Handle EPSG Units. Only meters for now. */
710 switch(nUnitsId) {
711 case 9001: /* Meters */
712 pszUnits = "m";
713 break;
714 default:
715 msSetError(MS_PROJERR,
716 "WMS/WFS AUTO PROJECTION: EPSG Units %d not supported.\n",
717 "_msProcessAutoProjection()", nUnitsId);
718 return -1;
719 }
720
721 /* Build PROJ4 definition.
722 * This is based on the definitions found in annex E of the WMS 1.1.1
723 * spec and online at http://www.digitalearth.org/wmt/auto.html
724 * The conversion from the original WKT definitions to PROJ4 format was
725 * done using the MapScript setWKTProjection() function (based on OGR).
726 */
727 switch(nProjId) {
728 case 42001: /** WGS 84 / Auto UTM **/
729 nZone = (int) floor( (dLon0 + 180.0) / 6.0 ) + 1;
730 sprintf( szProjBuf,
731 "+proj=tmerc+lat_0=0+lon_0=%.16g+k=0.999600+x_0=500000"
732 "+y_0=%.16g+ellps=WGS84+datum=WGS84+units=%s+type=crs",
733 -183.0 + nZone * 6.0,
734 (dLat0 >= 0.0) ? 0.0 : 10000000.0,
735 pszUnits);
736 break;
737 case 42002: /** WGS 84 / Auto Tr. Mercator **/
738 sprintf( szProjBuf,
739 "+proj=tmerc+lat_0=0+lon_0=%.16g+k=0.999600+x_0=500000"
740 "+y_0=%.16g+ellps=WGS84+datum=WGS84+units=%s+type=crs",
741 dLon0,
742 (dLat0 >= 0.0) ? 0.0 : 10000000.0,
743 pszUnits);
744 break;
745 case 42003: /** WGS 84 / Auto Orthographic **/
746 sprintf( szProjBuf,
747 "+proj=ortho+lon_0=%.16g+lat_0=%.16g+x_0=0+y_0=0"
748 "+ellps=WGS84+datum=WGS84+units=%s+type=crs",
749 dLon0, dLat0, pszUnits );
750 break;
751 case 42004: /** WGS 84 / Auto Equirectangular **/
752 /* Note that we have to pass lon_0 as lon_ts for this one to */
753 /* work. Either a PROJ4 bug or a PROJ4 documentation issue. */
754 sprintf( szProjBuf,
755 "+proj=eqc+lon_ts=%.16g+lat_ts=%.16g+x_0=0+y_0=0"
756 "+ellps=WGS84+datum=WGS84+units=%s+type=crs",
757 dLon0, dLat0, pszUnits);
758 break;
759 case 42005: /** WGS 84 / Auto Mollweide **/
760 sprintf( szProjBuf,
761 "+proj=moll+lon_0=%.16g+x_0=0+y_0=0+ellps=WGS84"
762 "+datum=WGS84+units=%s+type=crs",
763 dLon0, pszUnits);
764 break;
765 default:
766 msSetError(MS_PROJERR,
767 "WMS/WFS AUTO PROJECTION %d not supported.\n",
768 "_msProcessAutoProjection()", nProjId);
769 return -1;
770 }
771
772 /* msDebug("%s = %s\n", p->args[0], szProjBuf); */
773
774 /* OK, pass the definition to pj_init() */
775 args = msStringSplit(szProjBuf, '+', &numargs);
776
777 #if PROJ_VERSION_MAJOR >= 6
778 if( !(p->proj = proj_create_argv(p->proj_ctx->proj_ctx, numargs, args)) ) {
779 int l_pj_errno = proj_context_errno (p->proj_ctx->proj_ctx);
780 msSetError(MS_PROJERR, "proj error \"%s\" for \"%s\"",
781 "msProcessProjection()", proj_errno_string(l_pj_errno), szProjBuf) ;
782 return(-1);
783 }
784 #else
785 msAcquireLock( TLOCK_PROJ );
786 if( !(p->proj = pj_init(numargs, args)) ) {
787 int *pj_errno_ref = pj_get_errno_ref();
788 msReleaseLock( TLOCK_PROJ );
789 msSetError(MS_PROJERR, "proj error \"%s\" for \"%s\"",
790 "msProcessProjection()", pj_strerrno(*pj_errno_ref), szProjBuf) ;
791 return(-1);
792 }
793 msReleaseLock( TLOCK_PROJ );
794 #endif
795
796 msFreeCharArray(args, numargs);
797
798 return(0);
799 }
800
msProcessProjection(projectionObj * p)801 int msProcessProjection(projectionObj *p)
802 {
803 assert( p->proj == NULL );
804
805 if( strcasecmp(p->args[0],"GEOGRAPHIC") == 0 ) {
806 msSetError(MS_PROJERR,
807 "PROJECTION 'GEOGRAPHIC' no longer supported.\n"
808 "Provide explicit definition.\n"
809 "ie. proj=latlong\n"
810 " ellps=clrk66\n",
811 "msProcessProjection()");
812 return(-1);
813 }
814
815 if (strcasecmp(p->args[0], "AUTO") == 0) {
816 p->proj = NULL;
817 return 0;
818 }
819
820 #if PROJ_VERSION_MAJOR >= 6
821 if( p->proj_ctx == NULL )
822 {
823 p->proj_ctx = msProjectionContextCreate();
824 if( p->proj_ctx == NULL )
825 {
826 return -1;
827 }
828 }
829 if( p->proj_ctx->ms_proj_lib_change_counter != ms_proj_lib_change_counter )
830 {
831 msAcquireLock( TLOCK_PROJ );
832 p->proj_ctx->ms_proj_lib_change_counter = ms_proj_lib_change_counter;
833 {
834 const char* const paths[1] = { ms_proj_lib };
835 proj_context_set_search_paths(p->proj_ctx->proj_ctx, 1, ms_proj_lib ? paths : NULL);
836 }
837 msReleaseLock( TLOCK_PROJ );
838 }
839 #endif
840
841 if (strncasecmp(p->args[0], "AUTO:", 5) == 0 ||
842 strncasecmp(p->args[0], "AUTO2:", 6) == 0) {
843 /* WMS/WFS AUTO projection: "AUTO:proj_id,units_id,lon0,lat0" */
844 /*WMS 1.3.0: AUTO2:auto_crs_id,factor,lon0,lat0*/
845 return _msProcessAutoProjection(p);
846 }
847
848 #if PROJ_VERSION_MAJOR >= 6
849 {
850 char** args = (char**)msSmallMalloc(sizeof(char*) * (p->numargs+1));
851 memcpy(args, p->args, sizeof(char*) * p->numargs);
852
853 #if PROJ_VERSION_MAJOR == 6 && PROJ_VERSION_MINOR < 2
854 /* PROJ lookups are faster with EPSG in uppercase. Fixed in PROJ 6.2 */
855 /* Do that only for those versions, as it can create confusion if using */
856 /* a real old-style 'epsg' file... */
857 char szTemp[24];
858 if( p->numargs && strncmp(args[0], "init=epsg:", strlen("init=epsg:")) == 0 &&
859 strlen(args[0]) < 24)
860 {
861 strcpy(szTemp, "init=EPSG:");
862 strcat(szTemp, args[0] + strlen("init=epsg:"));
863 args[0] = szTemp;
864 }
865 #endif
866
867 args[p->numargs] = (char*) "type=crs";
868 if( !(p->proj = proj_create_argv(p->proj_ctx->proj_ctx, p->numargs + 1, args)) ) {
869 int l_pj_errno = proj_context_errno (p->proj_ctx->proj_ctx);
870 if(p->numargs>1) {
871 msSetError(MS_PROJERR, "proj error \"%s\" for \"%s:%s\"",
872 "msProcessProjection()", proj_errno_string(l_pj_errno), p->args[0],p->args[1]) ;
873 } else {
874 msSetError(MS_PROJERR, "proj error \"%s\" for \"%s\"",
875 "msProcessProjection()", proj_errno_string(l_pj_errno), p->args[0]) ;
876 }
877 free(args);
878 return(-1);
879 }
880 free(args);
881 }
882 #else
883 msAcquireLock( TLOCK_PROJ );
884 #if PJ_VERSION < 480
885 if( !(p->proj = pj_init(p->numargs, p->args)) ) {
886 #else
887 p->proj_ctx = pj_ctx_alloc();
888 if( !(p->proj=pj_init_ctx(p->proj_ctx, p->numargs, p->args)) ) {
889 #endif
890
891 int *pj_errno_ref = pj_get_errno_ref();
892 msReleaseLock( TLOCK_PROJ );
893 if(p->numargs>1) {
894 msSetError(MS_PROJERR, "proj error \"%s\" for \"%s:%s\"",
895 "msProcessProjection()", pj_strerrno(*pj_errno_ref), p->args[0],p->args[1]) ;
896 } else {
897 msSetError(MS_PROJERR, "proj error \"%s\" for \"%s\"",
898 "msProcessProjection()", pj_strerrno(*pj_errno_ref), p->args[0]) ;
899 }
900 return(-1);
901 }
902
903 msReleaseLock( TLOCK_PROJ );
904 #endif
905
906 #ifdef USE_PROJ_FASTPATHS
907 if(strcasestr(p->args[0],"epsg:4326")) {
908 p->wellknownprojection = wkp_lonlat;
909 } else if(strcasestr(p->args[0],"epsg:3857")) {
910 p->wellknownprojection = wkp_gmerc;
911 } else {
912 p->wellknownprojection = wkp_none;
913 }
914 #endif
915
916
917 return(0);
918 }
919
920 /************************************************************************/
921 /* int msIsAxisInverted */
922 /* Check to see if we should invert the axis. */
923 /* */
924 /************************************************************************/
925 int msIsAxisInverted(int epsg_code)
926 {
927 const unsigned int row = epsg_code / 8;
928 const unsigned char index = epsg_code % 8;
929
930 /*check the static table*/
931 if ((row < sizeof(axisOrientationEpsgCodes)) && (axisOrientationEpsgCodes[row] & (1 << index)))
932 return MS_TRUE;
933 else
934 return MS_FALSE;
935 }
936
937 /************************************************************************/
938 /* msProjectPoint() */
939 /************************************************************************/
940 int msProjectPoint(projectionObj *in, projectionObj *out, pointObj *point)
941 {
942 int ret;
943 reprojectionObj* reprojector = msProjectCreateReprojector(in, out);
944 if( !reprojector )
945 return MS_FAILURE;
946 ret = msProjectPointEx(reprojector, point);
947 msProjectDestroyReprojector(reprojector);
948 return ret;
949 }
950
951 /************************************************************************/
952 /* msProjectPointEx() */
953 /************************************************************************/
954 int msProjectPointEx(reprojectionObj* reprojector, pointObj *point)
955 {
956 projectionObj* in = reprojector->in;
957 projectionObj* out = reprojector->out;
958
959 if( in && in->gt.need_geotransform ) {
960 double x_out, y_out;
961
962 x_out = in->gt.geotransform[0]
963 + in->gt.geotransform[1] * point->x
964 + in->gt.geotransform[2] * point->y;
965 y_out = in->gt.geotransform[3]
966 + in->gt.geotransform[4] * point->x
967 + in->gt.geotransform[5] * point->y;
968
969 point->x = x_out;
970 point->y = y_out;
971 }
972
973 #if PROJ_VERSION_MAJOR >= 6
974 if( reprojector->pj ) {
975 PJ_COORD c;
976 c.xyzt.x = point->x;
977 c.xyzt.y = point->y;
978 c.xyzt.z = 0;
979 c.xyzt.t = 0;
980 c = proj_trans (reprojector->pj, PJ_FWD, c);
981 if( c.xyzt.x == HUGE_VAL || c.xyzt.y == HUGE_VAL ) {
982 return MS_FAILURE;
983 }
984 point->x = c.xyzt.x;
985 point->y = c.xyzt.y;
986 }
987 #else
988 if( reprojector->no_op ) {
989 /* do nothing, no transformation required */
990 }
991
992 /* -------------------------------------------------------------------- */
993 /* If we have a fully defined input coordinate system and */
994 /* output coordinate system, then we will use pj_transform. */
995 /* -------------------------------------------------------------------- */
996 else if( in && in->proj && out && out->proj ) {
997 int error;
998 double z = 0.0;
999
1000 if( pj_is_latlong(in->proj) ) {
1001 point->x *= DEG_TO_RAD;
1002 point->y *= DEG_TO_RAD;
1003 }
1004
1005 #if PJ_VERSION < 480
1006 msAcquireLock( TLOCK_PROJ );
1007 #endif
1008 error = pj_transform( in->proj, out->proj, 1, 0,
1009 &(point->x), &(point->y), &z );
1010 #if PJ_VERSION < 480
1011 msReleaseLock( TLOCK_PROJ );
1012 #endif
1013
1014 if( error || point->x == HUGE_VAL || point->y == HUGE_VAL ) {
1015 // msSetError(MS_PROJERR,"proj says: %s","msProjectPoint()",pj_strerrno(error));
1016 return MS_FAILURE;
1017 }
1018
1019 if( pj_is_latlong(out->proj) ) {
1020 point->x *= RAD_TO_DEG;
1021 point->y *= RAD_TO_DEG;
1022 }
1023 }
1024
1025 /* -------------------------------------------------------------------- */
1026 /* Otherwise we fallback to using pj_fwd() or pj_inv() and */
1027 /* assuming that the NULL projectionObj is supposed to be */
1028 /* lat/long in the same datum as the other projectionObj. This */
1029 /* is essentially a backwards compatibility mode. */
1030 /* -------------------------------------------------------------------- */
1031 else {
1032 projUV p;
1033
1034 p.u = point->x;
1035 p.v = point->y;
1036
1037 if(in==NULL || in->proj==NULL) { /* input coordinates are lat/lon */
1038 p.u *= DEG_TO_RAD; /* convert to radians */
1039 p.v *= DEG_TO_RAD;
1040 p = pj_fwd(p, out->proj);
1041 } else /* if(out==NULL || out->proj==NULL) */ { /* output coordinates are lat/lon */
1042 p = pj_inv(p, in->proj);
1043 p.u *= RAD_TO_DEG; /* convert to decimal degrees */
1044 p.v *= RAD_TO_DEG;
1045 }
1046
1047 point->x = p.u;
1048 point->y = p.v;
1049 if( point->x == HUGE_VAL || point->y == HUGE_VAL ) {
1050 return MS_FAILURE;
1051 }
1052 }
1053 #endif
1054
1055 if( out && out->gt.need_geotransform ) {
1056 double x_out, y_out;
1057
1058 x_out = out->gt.invgeotransform[0]
1059 + out->gt.invgeotransform[1] * point->x
1060 + out->gt.invgeotransform[2] * point->y;
1061 y_out = out->gt.invgeotransform[3]
1062 + out->gt.invgeotransform[4] * point->x
1063 + out->gt.invgeotransform[5] * point->y;
1064
1065 point->x = x_out;
1066 point->y = y_out;
1067 }
1068
1069 return(MS_SUCCESS);
1070 }
1071
1072 /************************************************************************/
1073 /* msProjectGrowRect() */
1074 /************************************************************************/
1075 static void msProjectGrowRect(reprojectionObj* reprojector,
1076 rectObj *prj_rect,
1077 pointObj *prj_point, int *failure )
1078
1079 {
1080 if( msProjectPointEx(reprojector, prj_point) == MS_SUCCESS ) {
1081 prj_rect->miny = MS_MIN(prj_rect->miny, prj_point->y);
1082 prj_rect->maxy = MS_MAX(prj_rect->maxy, prj_point->y);
1083 prj_rect->minx = MS_MIN(prj_rect->minx, prj_point->x);
1084 prj_rect->maxx = MS_MAX(prj_rect->maxx, prj_point->x);
1085 } else
1086 (*failure)++;
1087 }
1088
1089 /************************************************************************/
1090 /* msProjectSegment() */
1091 /* */
1092 /* Interpolate along a line segment for which one end */
1093 /* reprojects and the other end does not. Finds the horizon. */
1094 /************************************************************************/
1095 static int msProjectSegment( reprojectionObj* reprojector,
1096 pointObj *start, pointObj *end )
1097
1098 {
1099 pointObj testPoint, subStart, subEnd;
1100
1101 /* -------------------------------------------------------------------- */
1102 /* Without loss of generality we assume the first point */
1103 /* reprojects, and the second doesn't. If that is not the case */
1104 /* then re-call with the points reversed. */
1105 /* -------------------------------------------------------------------- */
1106 testPoint = *start;
1107 if( msProjectPointEx( reprojector, &testPoint ) == MS_FAILURE ) {
1108 testPoint = *end;
1109 if( msProjectPointEx( reprojector, &testPoint ) == MS_FAILURE )
1110 return MS_FAILURE;
1111 else
1112 return msProjectSegment( reprojector, end, start );
1113 }
1114
1115 /* -------------------------------------------------------------------- */
1116 /* We will apply a binary search till we are within out */
1117 /* tolerance. */
1118 /* -------------------------------------------------------------------- */
1119 subStart = *start;
1120 subEnd = *end;
1121
1122 #define TOLERANCE 0.01
1123
1124 while( fabs(subStart.x - subEnd.x)
1125 + fabs(subStart.y - subEnd.y) > TOLERANCE ) {
1126 pointObj midPoint;
1127
1128 midPoint.x = (subStart.x + subEnd.x) * 0.5;
1129 midPoint.y = (subStart.y + subEnd.y) * 0.5;
1130
1131 testPoint = midPoint;
1132
1133 if( msProjectPointEx( reprojector, &testPoint ) == MS_FAILURE ) {
1134 if (midPoint.x == subEnd.x && midPoint.y == subEnd.y)
1135 return MS_FAILURE; /* break infinite loop */
1136
1137 subEnd = midPoint;
1138 }
1139 else {
1140 if (midPoint.x == subStart.x && midPoint.y == subStart.y)
1141 return MS_FAILURE; /* break infinite loop */
1142
1143 subStart = midPoint;
1144 }
1145 }
1146
1147 /* -------------------------------------------------------------------- */
1148 /* Now reproject the end points and return. */
1149 /* -------------------------------------------------------------------- */
1150 *end = subStart;
1151
1152 if( msProjectPointEx( reprojector, end ) == MS_FAILURE
1153 || msProjectPointEx( reprojector, start ) == MS_FAILURE )
1154 return MS_FAILURE;
1155 else
1156 return MS_SUCCESS;
1157 }
1158
1159 /************************************************************************/
1160 /* msProjectShapeShouldDoLineCutting() */
1161 /************************************************************************/
1162
1163 /* Detect projecting from north polar stereographic to longlat or EPSG:3857 */
1164 #ifdef USE_GEOS
1165 static int msProjectShapeShouldDoLineCutting(reprojectionObj* reprojector)
1166 {
1167 if( reprojector->should_do_line_cutting >= 0)
1168 return reprojector->should_do_line_cutting;
1169
1170 projectionObj *in = reprojector->in;
1171 projectionObj *out = reprojector->out;
1172 if( !(!in->gt.need_geotransform && !msProjIsGeographicCRS(in) &&
1173 (msProjIsGeographicCRS(out) ||
1174 (out->numargs == 1 && strcmp(out->args[0], "init=epsg:3857") == 0))) )
1175 {
1176 reprojector->should_do_line_cutting = MS_FALSE;
1177 return MS_FALSE;
1178 }
1179
1180 int srcIsPolar;
1181 double extremeLongEasting;
1182 if( msProjIsGeographicCRS(out) )
1183 {
1184 pointObj p;
1185 double gt3 = out->gt.need_geotransform ? out->gt.geotransform[3] : 0.0;
1186 double gt4 = out->gt.need_geotransform ? out->gt.geotransform[4] : 0.0;
1187 double gt5 = out->gt.need_geotransform ? out->gt.geotransform[5] : 1.0;
1188 p.x = 0.0;
1189 p.y = 0.0;
1190 srcIsPolar = msProjectPointEx(reprojector, &p) == MS_SUCCESS &&
1191 fabs(gt3 + p.x * gt4 + p.y * gt5 - 90) < 1e-8;
1192 extremeLongEasting = 180;
1193 }
1194 else
1195 {
1196 pointObj p1;
1197 pointObj p2;
1198 double gt1 = out->gt.need_geotransform ? out->gt.geotransform[1] : 1.0;
1199 p1.x = 0.0;
1200 p1.y = -0.1;
1201 p2.x = 0.0;
1202 p2.y = 0.1;
1203 srcIsPolar = msProjectPointEx(reprojector, &p1) == MS_SUCCESS &&
1204 msProjectPointEx(reprojector, &p2) == MS_SUCCESS &&
1205 fabs((p1.x - p2.x) * gt1) > 20e6;
1206 extremeLongEasting = 20037508.3427892;
1207 }
1208 if( !srcIsPolar )
1209 {
1210 reprojector->should_do_line_cutting = MS_FALSE;
1211 return MS_FALSE;
1212 }
1213
1214 pointObj p = {0,0,0,0}; // initialize
1215 double invgt0 = out->gt.need_geotransform ? out->gt.invgeotransform[0] : 0.0;
1216 double invgt1 = out->gt.need_geotransform ? out->gt.invgeotransform[1] : 1.0;
1217 double invgt3 = out->gt.need_geotransform ? out->gt.invgeotransform[3] : 0.0;
1218 double invgt4 = out->gt.need_geotransform ? out->gt.invgeotransform[4] : 0.0;
1219
1220 lineObj newLine = {0,NULL};
1221 const double EPS = 1e-10;
1222
1223 p.x = invgt0 + -extremeLongEasting * (1 - EPS) * invgt1;
1224 p.y = invgt3 + -extremeLongEasting * (1 - EPS) * invgt4;
1225 msProjectPoint(out, in, &p);
1226 pointObj first = p;
1227 msAddPointToLine(&newLine, &p );
1228
1229 p.x = invgt0 + extremeLongEasting * (1 - EPS) * invgt1;
1230 p.y = invgt3 + extremeLongEasting * (1 - EPS) * invgt4;
1231 msProjectPoint(out, in, &p);
1232 msAddPointToLine(&newLine, &p );
1233
1234 p.x = 0;
1235 p.y = 0;
1236 msAddPointToLine(&newLine, &p );
1237
1238 msAddPointToLine(&newLine, &first );
1239
1240 msInitShape(&(reprojector->splitShape));
1241 reprojector->splitShape.type = MS_SHAPE_POLYGON;
1242 msAddLineDirectly(&(reprojector->splitShape), &newLine);
1243
1244 reprojector->should_do_line_cutting = MS_TRUE;
1245 return MS_TRUE;
1246 }
1247 #endif
1248
1249 /************************************************************************/
1250 /* msProjectShapeLine() */
1251 /* */
1252 /* Reproject a single linestring, potentially splitting into */
1253 /* more than one line string if there are over the horizon */
1254 /* portions. */
1255 /* */
1256 /* For polygons, no splitting takes place, but over the horizon */
1257 /* points are clipped, and one segment is run from the fall */
1258 /* over the horizon point to the come back over the horizon point. */
1259 /************************************************************************/
1260
1261 static int
1262 msProjectShapeLine(reprojectionObj* reprojector,
1263 shapeObj *shape, int line_index)
1264
1265 {
1266 int i;
1267 pointObj lastPoint, thisPoint, wrkPoint;
1268 lineObj *line = shape->line + line_index;
1269 lineObj *line_out = line;
1270 int valid_flag = 0; /* 1=true, -1=false, 0=unknown */
1271 int numpoints_in = line->numpoints;
1272 int line_alloc = numpoints_in;
1273 int wrap_test;
1274 projectionObj *in = reprojector->in;
1275 projectionObj *out = reprojector->out;
1276
1277 #ifdef USE_PROJ_FASTPATHS
1278 #define MAXEXTENT 20037508.34
1279 #define M_PIby360 .0087266462599716479
1280 #define MAXEXTENTby180 111319.4907777777777777777
1281 #define p_x line->point[i].x
1282 #define p_y line->point[i].y
1283 if(in->wellknownprojection == wkp_lonlat && out->wellknownprojection == wkp_gmerc) {
1284 for( i = line->numpoints-1; i >= 0; i-- ) {
1285 p_x *= MAXEXTENTby180;
1286 p_y = log(tan((90 + p_y) * M_PIby360)) * MS_RAD_TO_DEG;
1287 p_y *= MAXEXTENTby180;
1288 if (p_x > MAXEXTENT) p_x = MAXEXTENT;
1289 if (p_x < -MAXEXTENT) p_x = -MAXEXTENT;
1290 if (p_y > MAXEXTENT) p_y = MAXEXTENT;
1291 if (p_y < -MAXEXTENT) p_y = -MAXEXTENT;
1292 }
1293 return MS_SUCCESS;
1294 }
1295 if(in->wellknownprojection == wkp_gmerc && out->wellknownprojection == wkp_lonlat) {
1296 for( i = line->numpoints-1; i >= 0; i-- ) {
1297 if (p_x > MAXEXTENT) p_x = MAXEXTENT;
1298 else if (p_x < -MAXEXTENT) p_x = -MAXEXTENT;
1299 if (p_y > MAXEXTENT) p_y = MAXEXTENT;
1300 else if (p_y < -MAXEXTENT) p_y = -MAXEXTENT;
1301 p_x = (p_x / MAXEXTENT) * 180;
1302 p_y = (p_y / MAXEXTENT) * 180;
1303 p_y = MS_RAD_TO_DEG * (2 * atan(exp(p_y * MS_DEG_TO_RAD)) - MS_PI2);
1304 }
1305 msComputeBounds( shape ); /* fixes bug 1586 */
1306 return MS_SUCCESS;
1307 }
1308 #undef p_x
1309 #undef p_y
1310 #endif
1311
1312 #ifdef USE_GEOS
1313 if( shape->type == MS_SHAPE_LINE &&
1314 msProjectShapeShouldDoLineCutting(reprojector) )
1315 {
1316 shapeObj tmpShapeInputLine;
1317 msInitShape(&tmpShapeInputLine);
1318 tmpShapeInputLine.type = MS_SHAPE_LINE;
1319 tmpShapeInputLine.numlines = 1;
1320 tmpShapeInputLine.line = line;
1321
1322 shapeObj* diff = NULL;
1323 if( msGEOSIntersects(&tmpShapeInputLine, &(reprojector->splitShape)) )
1324 {
1325 diff = msGEOSDifference(&tmpShapeInputLine, &(reprojector->splitShape));
1326 }
1327
1328 tmpShapeInputLine.numlines = 0;
1329 tmpShapeInputLine.line = NULL;
1330 msFreeShape(&tmpShapeInputLine);
1331
1332 if( diff )
1333 {
1334 for(int j = 0; j < diff->numlines; j++ )
1335 {
1336 for( i=0; i < diff->line[j].numpoints; i++ ) {
1337 msProjectPointEx(reprojector, &(diff->line[j].point[i]));
1338 }
1339 if( j == 0 )
1340 {
1341 line_out->numpoints = diff->line[j].numpoints;
1342 memcpy(line_out->point, diff->line[0].point, sizeof(pointObj) * line_out->numpoints);
1343 }
1344 else
1345 {
1346 msAddLineDirectly(shape, &(diff->line[j]));
1347 }
1348 }
1349 msFreeShape(diff);
1350 msFree(diff);
1351 return MS_SUCCESS;
1352 }
1353 }
1354 #endif
1355
1356 wrap_test = out != NULL && out->proj != NULL && msProjIsGeographicCRS(out)
1357 && !msProjIsGeographicCRS(in);
1358
1359 line->numpoints = 0;
1360
1361 memset( &lastPoint, 0, sizeof(lastPoint) );
1362
1363 /* -------------------------------------------------------------------- */
1364 /* Loop over all input points in linestring. */
1365 /* -------------------------------------------------------------------- */
1366 for( i=0; i < numpoints_in; i++ ) {
1367 int ms_err;
1368 wrkPoint = thisPoint = line->point[i];
1369
1370 ms_err = msProjectPointEx(reprojector, &wrkPoint );
1371
1372 /* -------------------------------------------------------------------- */
1373 /* Apply wrap logic. */
1374 /* -------------------------------------------------------------------- */
1375 if( wrap_test && i > 0 && ms_err != MS_FAILURE ) {
1376 double dist;
1377 pointObj pt1Geo;
1378
1379 if( line_out->numpoints > 0 )
1380 pt1Geo = line_out->point[line_out->numpoints-1];
1381 else
1382 pt1Geo = wrkPoint; /* this is a cop out */
1383
1384 if( out->gt.need_geotransform && out->gt.geotransform[2] == 0 ) {
1385 dist = out->gt.geotransform[1] * (wrkPoint.x - pt1Geo.x);
1386 } else {
1387 dist = wrkPoint.x - pt1Geo.x;
1388 }
1389
1390 if( fabs(dist) > 180.0
1391 && msTestNeedWrap( thisPoint, lastPoint,
1392 pt1Geo, reprojector ) ) {
1393 if( out->gt.need_geotransform && out->gt.geotransform[2] == 0 ) {
1394 if( dist > 0.0 )
1395 wrkPoint.x -= 360.0 * out->gt.invgeotransform[1];
1396 else if( dist < 0.0 )
1397 wrkPoint.x += 360.0 * out->gt.invgeotransform[1];
1398 }
1399 else {
1400 if( dist > 0.0 )
1401 wrkPoint.x -= 360.0;
1402 else if( dist < 0.0 )
1403 wrkPoint.x += 360.0;
1404 }
1405 }
1406 }
1407
1408 /* -------------------------------------------------------------------- */
1409 /* Put result into output line with appropriate logic for */
1410 /* failure breaking lines, etc. */
1411 /* -------------------------------------------------------------------- */
1412 if( ms_err == MS_FAILURE ) {
1413 /* We have started out invalid */
1414 if( i == 0 ) {
1415 valid_flag = -1;
1416 }
1417
1418 /* valid data has ended, we need to work out the horizon */
1419 else if( valid_flag == 1 ) {
1420 pointObj startPoint, endPoint;
1421
1422 startPoint = lastPoint;
1423 endPoint = thisPoint;
1424 if( msProjectSegment( reprojector, &startPoint, &endPoint )
1425 == MS_SUCCESS ) {
1426 line_out->point[line_out->numpoints++] = endPoint;
1427 }
1428 valid_flag = -1;
1429 }
1430
1431 /* Still invalid ... */
1432 else if( valid_flag == -1 ) {
1433 /* do nothing */
1434 }
1435 }
1436
1437 else {
1438 /* starting out valid. */
1439 if( i == 0 ) {
1440 line_out->point[line_out->numpoints++] = wrkPoint;
1441 valid_flag = 1;
1442 }
1443
1444 /* Still valid, nothing special */
1445 else if( valid_flag == 1 ) {
1446 line_out->point[line_out->numpoints++] = wrkPoint;
1447 }
1448
1449 /* we have come over the horizon, figure out where, start newline*/
1450 else {
1451 pointObj startPoint, endPoint;
1452
1453 startPoint = lastPoint;
1454 endPoint = thisPoint;
1455 if( msProjectSegment( reprojector, &endPoint, &startPoint )
1456 == MS_SUCCESS ) {
1457 lineObj newLine;
1458
1459 /* force pre-allocation of lots of points room */
1460 if( line_out->numpoints > 0
1461 && shape->type == MS_SHAPE_LINE ) {
1462 newLine.numpoints = numpoints_in - i + 1;
1463 newLine.point = line->point;
1464 msAddLine( shape, &newLine );
1465
1466 /* new line is now lineout, but start without any points */
1467 line_out = shape->line + shape->numlines-1;
1468
1469 line_out->numpoints = 0;
1470
1471 /* the shape->line array is realloc, refetch "line" */
1472 line = shape->line + line_index;
1473 } else if( line_out == line
1474 && line->numpoints >= i-2 ) {
1475 newLine.numpoints = numpoints_in;
1476 newLine.point = line->point;
1477 msAddLine( shape, &newLine );
1478
1479 line = shape->line + line_index;
1480
1481 line_out = shape->line + shape->numlines-1;
1482 line_out->numpoints = line->numpoints;
1483 line->numpoints = 0;
1484
1485 /*
1486 * Now realloc this array large enough to hold all
1487 * the points we could possibly need to add.
1488 */
1489 line_alloc = line_alloc * 2;
1490
1491 line_out->point = (pointObj *)
1492 realloc(line_out->point,
1493 sizeof(pointObj) * line_alloc);
1494 }
1495
1496 line_out->point[line_out->numpoints++] = startPoint;
1497 }
1498 line_out->point[line_out->numpoints++] = wrkPoint;
1499 valid_flag = 1;
1500 }
1501 }
1502
1503 lastPoint = thisPoint;
1504 }
1505
1506 /* -------------------------------------------------------------------- */
1507 /* Make sure that polygons are closed, even if the trip over */
1508 /* the horizon left them unclosed. */
1509 /* -------------------------------------------------------------------- */
1510 if( shape->type == MS_SHAPE_POLYGON
1511 && line_out->numpoints > 2
1512 && (line_out->point[0].x != line_out->point[line_out->numpoints-1].x
1513 || line_out->point[0].y != line_out->point[line_out->numpoints-1].y) ) {
1514 /* make a copy because msAddPointToLine can realloc the array */
1515 pointObj sFirstPoint = line_out->point[0];
1516 msAddPointToLine( line_out, &sFirstPoint );
1517 }
1518
1519 return(MS_SUCCESS);
1520 }
1521
1522 /************************************************************************/
1523 /* msProjectShape() */
1524 /************************************************************************/
1525 int msProjectShape(projectionObj *in, projectionObj *out, shapeObj *shape)
1526 {
1527 int ret;
1528 reprojectionObj* reprojector = msProjectCreateReprojector(in, out);
1529 if( !reprojector )
1530 return MS_FAILURE;
1531 ret = msProjectShapeEx(reprojector, shape);
1532 msProjectDestroyReprojector(reprojector);
1533 return ret;
1534 }
1535
1536 /************************************************************************/
1537 /* msProjectShapeEx() */
1538 /************************************************************************/
1539 int msProjectShapeEx(reprojectionObj* reprojector, shapeObj *shape)
1540 {
1541 int i;
1542 #ifdef USE_PROJ_FASTPATHS
1543 int j;
1544
1545 #define p_x shape->line[i].point[j].x
1546 #define p_y shape->line[i].point[j].y
1547 if(in->wellknownprojection == wkp_lonlat && out->wellknownprojection == wkp_gmerc) {
1548 for( i = shape->numlines-1; i >= 0; i-- ) {
1549 for( j = shape->line[i].numpoints-1; j >= 0; j-- ) {
1550 p_x *= MAXEXTENTby180;
1551 p_y = log(tan((90 + p_y) * M_PIby360)) * MS_RAD_TO_DEG;
1552 p_y *= MAXEXTENTby180;
1553 if (p_x > MAXEXTENT) p_x = MAXEXTENT;
1554 else if (p_x < -MAXEXTENT) p_x = -MAXEXTENT;
1555 if (p_y > MAXEXTENT) p_y = MAXEXTENT;
1556 else if (p_y < -MAXEXTENT) p_y = -MAXEXTENT;
1557 }
1558 }
1559 msComputeBounds( shape ); /* fixes bug 1586 */
1560 return MS_SUCCESS;
1561 }
1562 if(in->wellknownprojection == wkp_gmerc && out->wellknownprojection == wkp_lonlat) {
1563 for( i = shape->numlines-1; i >= 0; i-- ) {
1564 for( j = shape->line[i].numpoints-1; j >= 0; j-- ) {
1565 if (p_x > MAXEXTENT) p_x = MAXEXTENT;
1566 else if (p_x < -MAXEXTENT) p_x = -MAXEXTENT;
1567 if (p_y > MAXEXTENT) p_y = MAXEXTENT;
1568 else if (p_y < -MAXEXTENT) p_y = -MAXEXTENT;
1569 p_x = (p_x / MAXEXTENT) * 180;
1570 p_y = (p_y / MAXEXTENT) * 180;
1571 p_y = MS_RAD_TO_DEG * (2 * atan(exp(p_y * MS_DEG_TO_RAD)) - MS_PI2);
1572 }
1573 }
1574 msComputeBounds( shape ); /* fixes bug 1586 */
1575 return MS_SUCCESS;
1576 }
1577 #undef p_x
1578 #undef p_y
1579 #endif
1580
1581 for( i = shape->numlines-1; i >= 0; i-- ) {
1582 if( shape->type == MS_SHAPE_LINE || shape->type == MS_SHAPE_POLYGON ) {
1583 if( msProjectShapeLine( reprojector, shape, i ) == MS_FAILURE )
1584 msShapeDeleteLine( shape, i );
1585 } else if( msProjectLineEx(reprojector, shape->line+i ) == MS_FAILURE ) {
1586 msShapeDeleteLine( shape, i );
1587 }
1588 }
1589
1590 if( shape->numlines == 0 ) {
1591 msFreeShape( shape );
1592 return MS_FAILURE;
1593 } else {
1594 msComputeBounds( shape ); /* fixes bug 1586 */
1595 return(MS_SUCCESS);
1596 }
1597 }
1598
1599 /************************************************************************/
1600 /* msProjectLine() */
1601 /************************************************************************/
1602 int msProjectLine(projectionObj *in, projectionObj *out, lineObj *line)
1603 {
1604 int ret;
1605 reprojectionObj* reprojector = msProjectCreateReprojector(in, out);
1606 if( !reprojector )
1607 return MS_FAILURE;
1608 ret = msProjectLineEx(reprojector, line);
1609 msProjectDestroyReprojector(reprojector);
1610 return ret;
1611 }
1612
1613 /************************************************************************/
1614 /* msProjectLineEx() */
1615 /* */
1616 /* This function is now normally only used for point data. */
1617 /* msProjectShapeLine() is used for lines and polygons and has */
1618 /* lots of logic to handle horizon crossing. */
1619 /************************************************************************/
1620
1621 int msProjectLineEx(reprojectionObj* reprojector, lineObj *line)
1622 {
1623 int i, be_careful = 1;
1624
1625 if( be_careful )
1626 be_careful = reprojector->out->proj != NULL &&
1627 msProjIsGeographicCRS(reprojector->out)
1628 && !msProjIsGeographicCRS(reprojector->in);
1629
1630 if( be_careful ) {
1631 pointObj startPoint, thisPoint; /* locations in projected space */
1632
1633 startPoint = line->point[0];
1634
1635 for(i=0; i<line->numpoints; i++) {
1636 double dist;
1637
1638 thisPoint = line->point[i];
1639
1640 /*
1641 ** Read comments before msTestNeedWrap() to better understand
1642 ** this dateline wrapping logic.
1643 */
1644 msProjectPointEx(reprojector, &(line->point[i]));
1645 if( i > 0 ) {
1646 dist = line->point[i].x - line->point[0].x;
1647 if( fabs(dist) > 180.0 ) {
1648 if( msTestNeedWrap( thisPoint, startPoint,
1649 line->point[0], reprojector ) ) {
1650 if( dist > 0.0 ) {
1651 line->point[i].x -= 360.0;
1652 } else if( dist < 0.0 ) {
1653 line->point[i].x += 360.0;
1654 }
1655 }
1656 }
1657
1658 }
1659 }
1660 } else {
1661 for(i=0; i<line->numpoints; i++) {
1662 if( msProjectPointEx(reprojector, &(line->point[i])) == MS_FAILURE )
1663 return MS_FAILURE;
1664 }
1665 }
1666
1667 return(MS_SUCCESS);
1668 }
1669
1670 /************************************************************************/
1671 /* msProjectRectGrid() */
1672 /************************************************************************/
1673
1674 #define NUMBER_OF_SAMPLE_POINTS 100
1675
1676 static
1677 int msProjectRectGrid(reprojectionObj* reprojector, rectObj *rect)
1678 {
1679 pointObj prj_point;
1680 rectObj prj_rect;
1681 int failure=0;
1682 int ix, iy;
1683
1684 double dx, dy;
1685 double x, y;
1686
1687 prj_rect.minx = DBL_MAX;
1688 prj_rect.miny = DBL_MAX;
1689 prj_rect.maxx = -DBL_MAX;
1690 prj_rect.maxy = -DBL_MAX;
1691
1692 dx = (rect->maxx - rect->minx)/NUMBER_OF_SAMPLE_POINTS;
1693 dy = (rect->maxy - rect->miny)/NUMBER_OF_SAMPLE_POINTS;
1694
1695 /* first ensure the top left corner is processed, even if the rect
1696 turns out to be degenerate. */
1697
1698 prj_point.x = rect->minx;
1699 prj_point.y = rect->miny;
1700 #ifdef USE_POINT_Z_M
1701 prj_point.z = 0.0;
1702 prj_point.m = 0.0;
1703 #endif /* USE_POINT_Z_M */
1704
1705 msProjectGrowRect(reprojector,&prj_rect,&prj_point,
1706 &failure);
1707
1708 failure = 0;
1709 for(ix = 0; ix <= NUMBER_OF_SAMPLE_POINTS; ix++ ) {
1710 x = rect->minx + ix * dx;
1711
1712 for(iy = 0; iy <= NUMBER_OF_SAMPLE_POINTS; iy++ ) {
1713 y = rect->miny + iy * dy;
1714
1715 prj_point.x = x;
1716 prj_point.y = y;
1717 msProjectGrowRect(reprojector,&prj_rect,&prj_point,
1718 &failure);
1719 }
1720 }
1721
1722 if( prj_rect.minx > prj_rect.maxx ) {
1723 rect->minx = 0;
1724 rect->maxx = 0;
1725 rect->miny = 0;
1726 rect->maxy = 0;
1727
1728 msSetError(MS_PROJERR, "All points failed to reproject.", "msProjectRect()");
1729 return MS_FAILURE;
1730 }
1731
1732 if( failure ) {
1733 msDebug( "msProjectRect(): some points failed to reproject, doing internal sampling.\n" );
1734 }
1735
1736 rect->minx = prj_rect.minx;
1737 rect->miny = prj_rect.miny;
1738 rect->maxx = prj_rect.maxx;
1739 rect->maxy = prj_rect.maxy;
1740
1741 return(MS_SUCCESS);
1742 }
1743
1744 /************************************************************************/
1745 /* msProjectRectAsPolygon() */
1746 /************************************************************************/
1747
1748 static int
1749 msProjectRectAsPolygon(reprojectionObj* reprojector, rectObj *rect)
1750 {
1751 shapeObj polygonObj;
1752 lineObj ring;
1753 /* pointObj ringPoints[NUMBER_OF_SAMPLE_POINTS*4+4]; */
1754 pointObj *ringPoints;
1755 int ix, iy, ixy;
1756
1757 double dx, dy;
1758
1759 /* If projecting from longlat to projected */
1760 /* This hack was introduced for WFS 2.0 compliance testing, but is far */
1761 /* from being perfect */
1762 if( reprojector->out && !msProjIsGeographicCRS(reprojector->out) &&
1763 reprojector->in && msProjIsGeographicCRS(reprojector->in) &&
1764 fabs(rect->minx - -180.0) < 1e-5 && fabs(rect->miny - -90.0) < 1e-5 &&
1765 fabs(rect->maxx - 180.0) < 1e-5 && fabs(rect->maxy - 90.0) < 1e-5) {
1766 pointObj pointTest;
1767 pointTest.x = -180;
1768 pointTest.y = 85;
1769 msProjectPointEx(reprojector, &pointTest);
1770 /* Detect if we are reprojecting from EPSG:4326 to EPSG:3857 */
1771 /* and if so use more plausible bounds to avoid issues with computed */
1772 /* resolution for WCS */
1773 if (fabs(pointTest.x - -20037508.3427892) < 1e-5 && fabs(pointTest.y - 19971868.8804086) )
1774 {
1775 rect->minx = -20037508.3427892;
1776 rect->miny = -20037508.3427892;
1777 rect->maxx = 20037508.3427892;
1778 rect->maxy = 20037508.3427892;
1779 }
1780 else
1781 {
1782 rect->minx = -1e15;
1783 rect->miny = -1e15;
1784 rect->maxx = 1e15;
1785 rect->maxy = 1e15;
1786 }
1787 return MS_SUCCESS;
1788 }
1789
1790 /* -------------------------------------------------------------------- */
1791 /* Build polygon as steps around the source rectangle */
1792 /* and possibly its diagonal. */
1793 /* -------------------------------------------------------------------- */
1794 dx = (rect->maxx - rect->minx)/NUMBER_OF_SAMPLE_POINTS;
1795 dy = (rect->maxy - rect->miny)/NUMBER_OF_SAMPLE_POINTS;
1796
1797 if(dx==0 && dy==0) {
1798 pointObj foo;
1799 msDebug( "msProjectRect(): Warning: degenerate rect {%f,%f,%f,%f}\n",rect->minx,rect->miny,rect->minx,rect->miny );
1800 foo.x = rect->minx;
1801 foo.y = rect->miny;
1802 msProjectPointEx(reprojector,&foo);
1803 rect->minx=rect->maxx=foo.x;
1804 rect->miny=rect->maxy=foo.y;
1805 return MS_SUCCESS;
1806 }
1807
1808 /* If there is more than two sample points we will also get samples from the diagonal line */
1809 ringPoints = (pointObj*) calloc(sizeof(pointObj),NUMBER_OF_SAMPLE_POINTS*5+3);
1810 ring.point = ringPoints;
1811 ring.numpoints = 0;
1812
1813 msInitShape( &polygonObj );
1814 polygonObj.type = MS_SHAPE_POLYGON;
1815
1816 /* sample along top */
1817 if(dx != 0) {
1818 for(ix = 0; ix <= NUMBER_OF_SAMPLE_POINTS; ix++ ) {
1819 ringPoints[ring.numpoints].x = rect->minx + ix * dx;
1820 ringPoints[ring.numpoints++].y = rect->miny;
1821 }
1822 }
1823
1824 /* sample along right side */
1825 if(dy != 0) {
1826 for(iy = 1; iy <= NUMBER_OF_SAMPLE_POINTS; iy++ ) {
1827 ringPoints[ring.numpoints].x = rect->maxx;
1828 ringPoints[ring.numpoints++].y = rect->miny + iy * dy;
1829 }
1830 }
1831
1832 /* sample along bottom */
1833 if(dx != 0) {
1834 for(ix = NUMBER_OF_SAMPLE_POINTS-1; ix >= 0; ix-- ) {
1835 ringPoints[ring.numpoints].x = rect->minx + ix * dx;
1836 ringPoints[ring.numpoints++].y = rect->maxy;
1837 }
1838 }
1839
1840 /* sample along left side */
1841 if(dy != 0) {
1842 for(iy = NUMBER_OF_SAMPLE_POINTS-1; iy >= 0; iy-- ) {
1843 ringPoints[ring.numpoints].x = rect->minx;
1844 ringPoints[ring.numpoints++].y = rect->miny + iy * dy;
1845 }
1846 }
1847
1848 /* sample along diagonal line */
1849 /* This is done to handle cases where reprojection from world covering projection to one */
1850 /* which isn't could cause min and max values of the projected rectangle to be invalid */
1851 if(dy != 0 && dx != 0) {
1852 /* No need to compute corners as they've already been computed */
1853 for(ixy = NUMBER_OF_SAMPLE_POINTS-2; ixy >= 1; ixy-- ) {
1854 ringPoints[ring.numpoints].x = rect->minx + ixy * dx;
1855 ringPoints[ring.numpoints++].y = rect->miny + ixy * dy;
1856 }
1857 }
1858
1859 msAddLineDirectly( &polygonObj, &ring );
1860
1861 #ifdef notdef
1862 FILE *wkt = fopen("/tmp/www-before.wkt","w");
1863 char *tmp = msShapeToWKT(&polygonObj);
1864 fprintf(wkt,"%s\n", tmp);
1865 free(tmp);
1866 fclose(wkt);
1867 #endif
1868
1869 /* -------------------------------------------------------------------- */
1870 /* Attempt to reproject. */
1871 /* -------------------------------------------------------------------- */
1872 msProjectShapeLine( reprojector, &polygonObj, 0 );
1873
1874 /* If no points reprojected, try a grid sampling */
1875 if( polygonObj.numlines == 0 || polygonObj.line[0].numpoints == 0 ) {
1876 msFreeShape( &polygonObj );
1877 return msProjectRectGrid( reprojector, rect );
1878 }
1879
1880 #ifdef notdef
1881 wkt = fopen("/tmp/www-after.wkt","w");
1882 tmp = msShapeToWKT(&polygonObj);
1883 fprintf(wkt,"%s\n", tmp);
1884 free(tmp);
1885 fclose(wkt);
1886 #endif
1887
1888 /* -------------------------------------------------------------------- */
1889 /* Collect bounds. */
1890 /* -------------------------------------------------------------------- */
1891 rect->minx = rect->maxx = polygonObj.line[0].point[0].x;
1892 rect->miny = rect->maxy = polygonObj.line[0].point[0].y;
1893
1894 for( ix = 1; ix < polygonObj.line[0].numpoints; ix++ ) {
1895 pointObj *pnt = polygonObj.line[0].point + ix;
1896
1897 rect->minx = MS_MIN(rect->minx,pnt->x);
1898 rect->maxx = MS_MAX(rect->maxx,pnt->x);
1899 rect->miny = MS_MIN(rect->miny,pnt->y);
1900 rect->maxy = MS_MAX(rect->maxy,pnt->y);
1901 }
1902
1903 msFreeShape( &polygonObj );
1904
1905 /* -------------------------------------------------------------------- */
1906 /* Special case to handle reprojection from "more than the */
1907 /* whole world" projected coordinates that sometimes produce a */
1908 /* region greater than 360 degrees wide due to various wrapping */
1909 /* logic. */
1910 /* -------------------------------------------------------------------- */
1911 if( reprojector->out && msProjIsGeographicCRS(reprojector->out) &&
1912 reprojector->in && !msProjIsGeographicCRS(reprojector->in)
1913 && rect->maxx - rect->minx > 360.0 &&
1914 !reprojector->out->gt.need_geotransform ) {
1915 rect->maxx = 180;
1916 rect->minx = -180;
1917 }
1918
1919 return MS_SUCCESS;
1920 }
1921
1922 /************************************************************************/
1923 /* msProjectHasLonWrap() */
1924 /************************************************************************/
1925
1926 int msProjectHasLonWrap(projectionObj *in, double* pdfLonWrap)
1927 {
1928 int i;
1929 if( pdfLonWrap )
1930 *pdfLonWrap = 0;
1931
1932 if( !msProjIsGeographicCRS(in) )
1933 return MS_FALSE;
1934
1935 for( i = 0; i < in->numargs; i++ )
1936 {
1937 if( strncmp(in->args[i], "lon_wrap=",
1938 strlen("lon_wrap=")) == 0 )
1939 {
1940 if( pdfLonWrap )
1941 *pdfLonWrap = atof(in->args[i] + strlen("lon_wrap="));
1942 return MS_TRUE;
1943 }
1944 }
1945 return MS_FALSE;
1946 }
1947
1948 /************************************************************************/
1949 /* msProjectRect() */
1950 /************************************************************************/
1951
1952 int msProjectRect(projectionObj *in, projectionObj *out, rectObj *rect)
1953 {
1954 char *over = "+over";
1955 int ret;
1956 int bFreeInOver = MS_FALSE;
1957 int bFreeOutOver = MS_FALSE;
1958 projectionObj in_over,out_over,*inp,*outp;
1959 double dfLonWrap = 0.0;
1960 reprojectionObj* reprojector = NULL;
1961
1962 /* Detect projecting from north polar stereographic to longlat */
1963 if( in && !in->gt.need_geotransform &&
1964 out && !out->gt.need_geotransform &&
1965 !msProjIsGeographicCRS(in) && msProjIsGeographicCRS(out) )
1966 {
1967 pointObj p;
1968 p.x = 0.0;
1969 p.y = 0.0;
1970 reprojector = msProjectCreateReprojector(in, out);
1971 if( reprojector &&
1972 msProjectPointEx(reprojector, &p) == MS_SUCCESS &&
1973 fabs(p.y - 90) < 1e-8 )
1974 {
1975 /* Is the pole in the rectangle ? */
1976 if( 0 >= rect->minx && 0 >= rect->miny &&
1977 0 <= rect->maxx && 0 <= rect->maxy )
1978 {
1979 if( msProjectRectAsPolygon(reprojector, rect ) == MS_SUCCESS )
1980 {
1981 rect->minx = -180.0;
1982 rect->maxx = 180.0;
1983 rect->maxy = 90.0;
1984 msProjectDestroyReprojector(reprojector);
1985 return MS_SUCCESS;
1986 }
1987 }
1988 /* Are we sure the dateline is not enclosed ? */
1989 else if( rect->maxy < 0 || rect->maxx < 0 || rect->minx > 0 )
1990 {
1991 ret = msProjectRectAsPolygon(reprojector, rect );
1992 msProjectDestroyReprojector(reprojector);
1993 return ret;
1994 }
1995 }
1996 }
1997
1998 if(in && msProjectHasLonWrap(in, &dfLonWrap) && dfLonWrap == 180.0) {
1999 inp = in;
2000 outp = out;
2001 if( rect->maxx > 180.0 ) {
2002 rect->minx = -180.0;
2003 rect->maxx = 180.0;
2004 }
2005 }
2006 /*
2007 * Issue #4892: When projecting a rectangle we do not want proj to wrap resulting
2008 * coordinates around the dateline, as in practice a requested bounding box of
2009 * -22.000.000, -YYY, 22.000.000, YYY should be projected as
2010 * -190,-YYY,190,YYY rather than 170,-YYY,-170,YYY as would happen when wrapping (and
2011 * vice-versa when projecting from lonlat to metric).
2012 * To enforce this, we clone the input projections and add the "+over" proj
2013 * parameter in order to disable dateline wrapping.
2014 */
2015 else {
2016 if(out) {
2017 bFreeOutOver = MS_TRUE;
2018 msInitProjection(&out_over);
2019 msCopyProjectionExtended(&out_over,out,&over,1);
2020 outp = &out_over;
2021 if( reprojector ) {
2022 msProjectDestroyReprojector(reprojector);
2023 reprojector = NULL;
2024 }
2025 } else {
2026 outp = out;
2027 }
2028 if(in) {
2029 bFreeInOver = MS_TRUE;
2030 msInitProjection(&in_over);
2031 msCopyProjectionExtended(&in_over,in,&over,1);
2032 inp = &in_over;
2033 if( reprojector ) {
2034 msProjectDestroyReprojector(reprojector);
2035 reprojector = NULL;
2036 }
2037 } else {
2038 inp = in;
2039 }
2040 }
2041 if( reprojector == NULL )
2042 {
2043 reprojector = msProjectCreateReprojector(inp, outp);
2044 }
2045 ret = reprojector ? msProjectRectAsPolygon(reprojector, rect ) : MS_FAILURE;
2046 msProjectDestroyReprojector(reprojector);
2047 if(bFreeInOver)
2048 msFreeProjection(&in_over);
2049 if(bFreeOutOver)
2050 msFreeProjection(&out_over);
2051 return ret;
2052 }
2053
2054 #if PROJ_VERSION_MAJOR < 6
2055
2056 static int msProjectSortString(const void* firstelt, const void* secondelt)
2057 {
2058 char* firststr = *(char**)firstelt;
2059 char* secondstr = *(char**)secondelt;
2060 return strcmp(firststr, secondstr);
2061 }
2062
2063 /************************************************************************/
2064 /* msGetProjectNormalized() */
2065 /************************************************************************/
2066
2067 static projectionObj* msGetProjectNormalized( const projectionObj* p )
2068 {
2069 int i;
2070 #if PROJ_VERSION_MAJOR >= 6
2071 const char* pszNewProj4Def;
2072 #else
2073 char* pszNewProj4Def;
2074 #endif
2075 projectionObj* pnew;
2076
2077 pnew = (projectionObj*)msSmallMalloc(sizeof(projectionObj));
2078 msInitProjection(pnew);
2079 msCopyProjection(pnew, (projectionObj*)p);
2080
2081 if(p->proj == NULL )
2082 return pnew;
2083
2084 /* Normalize definition so that msProjectDiffers() works better */
2085 #if PROJ_VERSION_MAJOR >= 6
2086 pszNewProj4Def = proj_as_proj_string(p->proj_ctx->proj_ctx, p->proj, PJ_PROJ_4, NULL);
2087 #else
2088 pszNewProj4Def = pj_get_def( p->proj, 0 );
2089 #endif
2090 msFreeCharArray(pnew->args, pnew->numargs);
2091 pnew->args = msStringSplit(pszNewProj4Def,'+', &pnew->numargs);
2092 for(i = 0; i < pnew->numargs; i++)
2093 {
2094 /* Remove trailing space */
2095 if( strlen(pnew->args[i]) > 0 && pnew->args[i][strlen(pnew->args[i])-1] == ' ' )
2096 pnew->args[i][strlen(pnew->args[i])-1] = '\0';
2097 /* Remove spurious no_defs or init= */
2098 if( strcmp(pnew->args[i], "no_defs") == 0 ||
2099 strncmp(pnew->args[i], "init=", 5) == 0 )
2100 {
2101 if( i < pnew->numargs - 1 )
2102 {
2103 msFree(pnew->args[i]);
2104 memmove(pnew->args + i, pnew->args + i + 1,
2105 sizeof(char*) * (pnew->numargs - 1 -i ));
2106 }
2107 else
2108 {
2109 msFree(pnew->args[i]);
2110 }
2111 pnew->numargs --;
2112 i --;
2113 continue;
2114 }
2115 }
2116 /* Sort the strings so they can be compared */
2117 qsort(pnew->args, pnew->numargs, sizeof(char*), msProjectSortString);
2118 /*{
2119 fprintf(stderr, "'%s' =\n", pszNewProj4Def);
2120 for(i = 0; i < p->numargs; i++)
2121 fprintf(stderr, "'%s' ", p->args[i]);
2122 fprintf(stderr, "\n");
2123 }*/
2124 #if PROJ_VERSION_MAJOR < 6
2125 pj_dalloc(pszNewProj4Def);
2126 #endif
2127
2128 return pnew;
2129 }
2130 #endif
2131
2132 /************************************************************************/
2133 /* msProjectionsDiffer() */
2134 /************************************************************************/
2135
2136 /*
2137 ** Compare two projections, and return MS_TRUE if they differ.
2138 **
2139 ** For now this is implemented by exact comparison of the projection
2140 ** arguments, but eventually this should call a PROJ.4 function with
2141 ** more awareness of the issues.
2142 **
2143 ** NOTE: MS_FALSE is returned if either of the projection objects
2144 ** has no arguments, since reprojection won't work anyways.
2145 */
2146
2147 static int msProjectionsDifferInternal( projectionObj *proj1, projectionObj *proj2 )
2148
2149 {
2150 int i;
2151
2152 if( proj1->numargs == 0 || proj2->numargs == 0 )
2153 return MS_FALSE;
2154
2155 if( proj1->numargs != proj2->numargs )
2156 return MS_TRUE;
2157
2158 /* This test should be more rigerous. */
2159 if( proj1->gt.need_geotransform
2160 || proj2->gt.need_geotransform )
2161 return MS_TRUE;
2162
2163 for( i = 0; i < proj1->numargs; i++ ) {
2164 if( strcmp(proj1->args[i],proj2->args[i]) != 0 )
2165 return MS_TRUE;
2166 }
2167
2168 return MS_FALSE;
2169 }
2170
2171 int msProjectionsDiffer( projectionObj *proj1, projectionObj *proj2 )
2172 {
2173 int ret;
2174
2175 ret = msProjectionsDifferInternal(proj1, proj2);
2176 #if PROJ_VERSION_MAJOR < 6
2177 if( ret &&
2178 /* to speed up things, do normalization only if one proj is */
2179 /* likely of the form init=epsg:XXX and the other proj=XXX datum=YYY... */
2180 ( (proj1->numargs == 1 && proj2->numargs > 1) ||
2181 (proj1->numargs > 1 && proj2->numargs == 1) ) )
2182 {
2183 projectionObj* p1normalized;
2184 projectionObj* p2normalized;
2185
2186 p1normalized = msGetProjectNormalized( proj1 );
2187 p2normalized = msGetProjectNormalized( proj2 );
2188 ret = msProjectionsDifferInternal(p1normalized, p2normalized);
2189 msFreeProjection(p1normalized);
2190 msFree(p1normalized);
2191 msFreeProjection(p2normalized);
2192 msFree(p2normalized);
2193 }
2194 #endif
2195 return ret;
2196 }
2197
2198 /************************************************************************/
2199 /* msTestNeedWrap() */
2200 /************************************************************************/
2201 /*
2202
2203 Frank Warmerdam, Nov, 2001.
2204
2205 See Also:
2206
2207 http://mapserver.gis.umn.edu/bugs/show_bug.cgi?id=15
2208
2209 Proposal:
2210
2211 Modify msProjectLine() so that it "dateline wraps" objects when necessary
2212 in order to preserve their shape when reprojecting to lat/long. This
2213 will be accomplished by:
2214
2215 1) As each vertex is reprojected, compare the X distance between that
2216 vertex and the previous vertex. If it is less than 180 then proceed to
2217 the next vertex without any special logic, otherwise:
2218
2219 2) Reproject the center point of the line segment from the last vertex to
2220 the current vertex into lat/long. Does it's longitude lie between the
2221 longitudes of the start and end point. If yes, return to step 1) for
2222 the next vertex ... everything is fine.
2223
2224 3) We have determined that this line segment is suffering from 360 degree
2225 wrap to keep in the -180 to +180 range. Now add or subtract 360 degrees
2226 as determined from the original sign of the distances.
2227
2228 This is similar to the code there now (though disabled in CVS); however,
2229 it will ensure that big boxes will remain big, and not get dateline wrapped
2230 because of the extra test in step 2). However step 2 is invoked only very
2231 rarely so this process takes little more than the normal process. In fact,
2232 if we were sufficiently concerned about performance we could do a test on the
2233 shape MBR in lat/long space, and if the width is less than 180 we know we never
2234 need to even do test 1).
2235
2236 What doesn't this resolve:
2237
2238 This ensures that individual lines are kept in the proper shape when
2239 reprojected to geographic space. However, it does not:
2240
2241 o Ensure that all rings of a polygon will get transformed to the same "side"
2242 of the world. Depending on starting points of the different rings it is
2243 entirely possible for one ring to end up in the -180 area and another ring
2244 from the same polygon to end up in the +180 area. We might possibly be
2245 able to achieve this though, by treating the multi-ring polygon as a whole
2246 and testing the first point of each additional ring against the last
2247 vertex of the previous ring (or any previous vertex for that matter).
2248
2249 o It does not address the need to cut up lines and polygons into distinct
2250 chunks to preserve the correct semantics. Really a polygon that
2251 spaces the dateline in a -180 to 180 view should get split into two
2252 polygons. We haven't addressed that, though if it were to be addressed,
2253 it could be done as a followon and distinct step from what we are doing
2254 above. In fact this sort of improvement (split polygons based on dateline
2255 or view window) should be done for all lat/long shapes regardless of
2256 whether they are being reprojected from another projection.
2257
2258 o It does not address issues related to viewing rectangles that go outside
2259 the -180 to 180 longitude range. For instance, it is entirely reasonable
2260 to want a 160 to 200 longitude view to see an area on the dateline clearly.
2261 Currently shapes in the -180 to -160 range which should be displayed in the
2262 180 to 200 portion of that view will not be because there is no recogition
2263 that they belong there.
2264
2265
2266 */
2267
2268 static int msTestNeedWrap( pointObj pt1, pointObj pt2, pointObj pt2_geo,
2269 reprojectionObj* reprojector )
2270
2271 {
2272 pointObj middle;
2273 projectionObj* out = reprojector->out;
2274
2275 middle.x = (pt1.x + pt2.x) * 0.5;
2276 middle.y = (pt1.y + pt2.y) * 0.5;
2277
2278 if( msProjectPointEx( reprojector, &pt1 ) == MS_FAILURE
2279 || msProjectPointEx( reprojector, &pt2 ) == MS_FAILURE
2280 || msProjectPointEx( reprojector, &middle ) == MS_FAILURE )
2281 return 0;
2282
2283 /*
2284 * If the last point was moved, then we are considered due for a
2285 * move to.
2286 */
2287 if( out->gt.need_geotransform && out->gt.geotransform[2] == 0 ) {
2288 if( fabs( (pt2_geo.x-pt2.x) * out->gt.geotransform[1] ) > 180.0 )
2289 return 1;
2290 } else {
2291 if( fabs(pt2_geo.x-pt2.x) > 180.0 )
2292 return 1;
2293 }
2294
2295 /*
2296 * Otherwise, test to see if the middle point transforms
2297 * to be between the end points. If yes, no wrapping is needed.
2298 * Otherwise wrapping is needed.
2299 */
2300 if( (middle.x < pt1.x && middle.x < pt2_geo.x)
2301 || (middle.x > pt1.x && middle.x > pt2_geo.x) )
2302 return 1;
2303 else
2304 return 0;
2305 }
2306
2307 /************************************************************************/
2308 /* msProjFinder() */
2309 /************************************************************************/
2310
2311 #if PROJ_VERSION_MAJOR < 6
2312 static char *last_filename = NULL;
2313
2314 static const char *msProjFinder( const char *filename)
2315
2316 {
2317 if( last_filename != NULL )
2318 free( last_filename );
2319
2320 if( filename == NULL )
2321 return NULL;
2322
2323 if( ms_proj_lib == NULL )
2324 return filename;
2325
2326 last_filename = (char *) malloc(strlen(filename)+strlen(ms_proj_lib)+2);
2327 sprintf( last_filename, "%s/%s", ms_proj_lib, filename );
2328
2329 return last_filename;
2330 }
2331 #endif
2332
2333 /************************************************************************/
2334 /* msProjLibInitFromEnv() */
2335 /************************************************************************/
2336 void msProjLibInitFromEnv()
2337 {
2338 const char *val;
2339
2340 if( (val=getenv( "PROJ_LIB" )) != NULL ) {
2341 msSetPROJ_LIB(val, NULL);
2342 }
2343 }
2344
2345 /************************************************************************/
2346 /* msSetPROJ_LIB() */
2347 /************************************************************************/
2348 void msSetPROJ_LIB( const char *proj_lib, const char *pszRelToPath )
2349
2350 {
2351 char *extended_path = NULL;
2352
2353 /* Handle relative path if applicable */
2354 if( proj_lib && pszRelToPath
2355 && proj_lib[0] != '/'
2356 && proj_lib[0] != '\\'
2357 && !(proj_lib[0] != '\0' && proj_lib[1] == ':') ) {
2358 struct stat stat_buf;
2359 extended_path = (char*) msSmallMalloc(strlen(pszRelToPath)
2360 + strlen(proj_lib) + 10);
2361 sprintf( extended_path, "%s/%s", pszRelToPath, proj_lib );
2362
2363 #ifndef S_ISDIR
2364 # define S_ISDIR(x) ((x) & S_IFDIR)
2365 #endif
2366
2367 if( stat( extended_path, &stat_buf ) == 0
2368 && S_ISDIR(stat_buf.st_mode) )
2369 proj_lib = extended_path;
2370 }
2371
2372
2373 msAcquireLock( TLOCK_PROJ );
2374 #if PROJ_VERSION_MAJOR >= 6
2375 if( proj_lib == NULL && ms_proj_lib == NULL )
2376 {
2377 /* do nothing */
2378 }
2379 else if( proj_lib != NULL && ms_proj_lib != NULL &&
2380 strcmp(proj_lib, ms_proj_lib) == 0 )
2381 {
2382 /* do nothing */
2383 }
2384 else
2385 {
2386 ms_proj_lib_change_counter++;
2387 free( ms_proj_lib );
2388 ms_proj_lib = proj_lib ? msStrdup(proj_lib) : NULL;
2389 }
2390 #else
2391 {
2392 static int finder_installed = 0;
2393 if( finder_installed == 0 && proj_lib != NULL) {
2394 finder_installed = 1;
2395 pj_set_finder( msProjFinder );
2396 }
2397 }
2398
2399 if (proj_lib == NULL) pj_set_finder(NULL);
2400
2401 if( ms_proj_lib != NULL ) {
2402 free( ms_proj_lib );
2403 ms_proj_lib = NULL;
2404 }
2405
2406 if( last_filename != NULL ) {
2407 free( last_filename );
2408 last_filename = NULL;
2409 }
2410
2411 if( proj_lib != NULL )
2412 ms_proj_lib = msStrdup( proj_lib );
2413 #endif
2414 msReleaseLock( TLOCK_PROJ );
2415
2416 #if GDAL_VERSION_MAJOR >= 3
2417 if( ms_proj_lib != NULL )
2418 {
2419 const char* const apszPaths[] = { ms_proj_lib, NULL };
2420 OSRSetPROJSearchPaths(apszPaths);
2421 }
2422 #endif
2423
2424 if ( extended_path )
2425 msFree( extended_path );
2426 }
2427
2428 /************************************************************************/
2429 /* msGetProjectionString() */
2430 /* */
2431 /* Return the projection string. */
2432 /************************************************************************/
2433
2434 char *msGetProjectionString(projectionObj *proj)
2435 {
2436 char *pszProjString = NULL;
2437 int i = 0, nLen = 0;
2438
2439 if (proj) {
2440 /* -------------------------------------------------------------------- */
2441 /* Alloc buffer large enough to hold the whole projection defn */
2442 /* -------------------------------------------------------------------- */
2443 for (i=0; i<proj->numargs; i++) {
2444 if (proj->args[i])
2445 nLen += (strlen(proj->args[i]) + 2);
2446 }
2447
2448 pszProjString = (char*)malloc(sizeof(char) * nLen+1);
2449 pszProjString[0] = '\0';
2450
2451 /* -------------------------------------------------------------------- */
2452 /* Plug each arg into the string with a '+' prefix */
2453 /* -------------------------------------------------------------------- */
2454 for (i=0; i<proj->numargs; i++) {
2455 if (!proj->args[i] || strlen(proj->args[i]) == 0)
2456 continue;
2457 if (pszProjString[0] == '\0') {
2458 /* no space at beginning of line */
2459 if (proj->args[i][0] != '+')
2460 strcat(pszProjString, "+");
2461 } else {
2462 if (proj->args[i][0] != '+')
2463 strcat(pszProjString, " +");
2464 else
2465 strcat(pszProjString, " ");
2466 }
2467 strcat(pszProjString, proj->args[i]);
2468 }
2469 }
2470
2471 return pszProjString;
2472 }
2473
2474 /************************************************************************/
2475 /* msIsAxisInvertedProj() */
2476 /* */
2477 /* Return MS_TRUE is the proj object has epgsaxis=ne */
2478 /************************************************************************/
2479
2480 int msIsAxisInvertedProj( projectionObj *proj )
2481 {
2482 int i;
2483 const char *axis = NULL;
2484
2485 for( i = 0; i < proj->numargs; i++ ) {
2486 if( strstr(proj->args[i],"epsgaxis=") != NULL ) {
2487 axis = strstr(proj->args[i],"=") + 1;
2488 break;
2489 }
2490 }
2491
2492 if( axis == NULL )
2493 return MS_FALSE;
2494
2495 if( strcasecmp(axis,"en") == 0 )
2496 return MS_FALSE;
2497
2498 if( strcasecmp(axis,"ne") != 0 ) {
2499 msDebug( "msIsAxisInvertedProj(): odd +epsgaxis= value: '%s'.",
2500 axis );
2501 return MS_FALSE;
2502 }
2503
2504 return MS_TRUE;
2505 }
2506
2507 /************************************************************************/
2508 /* msAxisNormalizePoints() */
2509 /* */
2510 /* Convert the passed points to "easting, northing" axis */
2511 /* orientation if they are not already. */
2512 /************************************************************************/
2513
2514 void msAxisNormalizePoints( projectionObj *proj, int count,
2515 double *x, double *y )
2516
2517 {
2518 int i;
2519 if( !msIsAxisInvertedProj(proj ) )
2520 return;
2521
2522 /* Switch axes */
2523 for( i = 0; i < count; i++ ) {
2524 double tmp;
2525
2526 tmp = x[i];
2527 x[i] = y[i];
2528 y[i] = tmp;
2529 }
2530 }
2531
2532
2533
2534 /************************************************************************/
2535 /* msAxisSwapShape */
2536 /* */
2537 /* Utility function to swap x and y coordiatesn Use for now for */
2538 /* WFS 1.1.x */
2539 /************************************************************************/
2540 void msAxisSwapShape(shapeObj *shape)
2541 {
2542 double tmp;
2543 int i,j;
2544
2545 if (shape) {
2546 for(i=0; i<shape->numlines; i++) {
2547 for( j=0; j<shape->line[i].numpoints; j++ ) {
2548 tmp = shape->line[i].point[j].x;
2549 shape->line[i].point[j].x = shape->line[i].point[j].y;
2550 shape->line[i].point[j].y = tmp;
2551 }
2552 }
2553
2554 /*swap bounds*/
2555 tmp = shape->bounds.minx;
2556 shape->bounds.minx = shape->bounds.miny;
2557 shape->bounds.miny = tmp;
2558
2559 tmp = shape->bounds.maxx;
2560 shape->bounds.maxx = shape->bounds.maxy;
2561 shape->bounds.maxy = tmp;
2562 }
2563 }
2564
2565 /************************************************************************/
2566 /* msAxisDenormalizePoints() */
2567 /* */
2568 /* Convert points from easting,northing orientation to the */
2569 /* preferred epsg orientation of this projectionObj. */
2570 /************************************************************************/
2571
2572 void msAxisDenormalizePoints( projectionObj *proj, int count,
2573 double *x, double *y )
2574
2575 {
2576 /* For how this is essentially identical to normalizing */
2577 msAxisNormalizePoints( proj, count, x, y );
2578 }
2579
2580 /************************************************************************/
2581 /* msProjIsGeographicCRS() */
2582 /* */
2583 /* Returns whether a CRS is a geographic one. */
2584 /************************************************************************/
2585
2586 int msProjIsGeographicCRS(projectionObj* proj)
2587 {
2588 #if PROJ_VERSION_MAJOR >= 6
2589 PJ_TYPE type;
2590 if( !proj->proj )
2591 return FALSE;
2592 type = proj_get_type(proj->proj);
2593 if( type == PJ_TYPE_GEOGRAPHIC_2D_CRS || type == PJ_TYPE_GEOGRAPHIC_3D_CRS )
2594 return TRUE;
2595 if( type == PJ_TYPE_BOUND_CRS )
2596 {
2597 PJ* base_crs = proj_get_source_crs(proj->proj_ctx->proj_ctx, proj->proj);
2598 type = proj_get_type(base_crs);
2599 proj_destroy(base_crs);
2600 return type == PJ_TYPE_GEOGRAPHIC_2D_CRS || type == PJ_TYPE_GEOGRAPHIC_3D_CRS;
2601 }
2602 return FALSE;
2603 #else
2604 return pj_is_latlong(proj->proj);
2605 #endif
2606 }
2607
2608 /************************************************************************/
2609 /* ConvertProjUnitStringToMS */
2610 /* */
2611 /* Returns mapserver's unit code corresponding to the proj */
2612 /* unit passed as argument. */
2613 /* Please refer to ./src/pj_units.c file in the Proj.4 module. */
2614 /************************************************************************/
2615 static int ConvertProjUnitStringToMS(const char *pszProjUnit)
2616 {
2617 if (strcmp(pszProjUnit, "m") ==0) {
2618 return MS_METERS;
2619 } else if (strcmp(pszProjUnit, "km") ==0) {
2620 return MS_KILOMETERS;
2621 } else if (strcmp(pszProjUnit, "mi") ==0 || strcmp(pszProjUnit, "us-mi") ==0) {
2622 return MS_MILES;
2623 } else if (strcmp(pszProjUnit, "in") ==0 || strcmp(pszProjUnit, "us-in") ==0 ) {
2624 return MS_INCHES;
2625 } else if (strcmp(pszProjUnit, "ft") ==0 || strcmp(pszProjUnit, "us-ft") ==0) {
2626 return MS_FEET;
2627 } else if (strcmp(pszProjUnit, "kmi") == 0) {
2628 return MS_NAUTICALMILES;
2629 }
2630
2631 return -1;
2632 }
2633
2634 /************************************************************************/
2635 /* int GetMapserverUnitUsingProj(projectionObj *psProj) */
2636 /* */
2637 /* Return a mapserver unit corresponding to the projection */
2638 /* passed. Retunr -1 on failure */
2639 /************************************************************************/
2640 int GetMapserverUnitUsingProj(projectionObj *psProj)
2641 {
2642 #if PROJ_VERSION_MAJOR >= 6
2643 const char *proj_str;
2644 #else
2645 char *proj_str;
2646 #endif
2647
2648 if( msProjIsGeographicCRS( psProj ) )
2649 return MS_DD;
2650
2651 #if PROJ_VERSION_MAJOR >= 6
2652 proj_str = proj_as_proj_string(psProj->proj_ctx->proj_ctx, psProj->proj, PJ_PROJ_4, NULL);
2653 #else
2654 proj_str = pj_get_def( psProj->proj, 0 );
2655 #endif
2656
2657 /* -------------------------------------------------------------------- */
2658 /* Handle case of named units. */
2659 /* -------------------------------------------------------------------- */
2660 if( strstr(proj_str,"units=") != NULL ) {
2661 char units[32];
2662 char *blank;
2663
2664 strlcpy( units, (strstr(proj_str,"units=")+6), sizeof(units) );
2665 #if PROJ_VERSION_MAJOR < 6
2666 pj_dalloc( proj_str );
2667 #endif
2668
2669 blank = strchr(units, ' ');
2670 if( blank != NULL )
2671 *blank = '\0';
2672
2673 return ConvertProjUnitStringToMS( units );
2674 }
2675
2676 /* -------------------------------------------------------------------- */
2677 /* Handle case of to_meter value. */
2678 /* -------------------------------------------------------------------- */
2679 if( strstr(proj_str,"to_meter=") != NULL ) {
2680 char to_meter_str[32];
2681 char *blank;
2682 double to_meter;
2683
2684 strlcpy(to_meter_str,(strstr(proj_str,"to_meter=")+9),
2685 sizeof(to_meter_str));
2686 #if PROJ_VERSION_MAJOR < 6
2687 pj_dalloc( proj_str );
2688 #endif
2689
2690 blank = strchr(to_meter_str, ' ');
2691 if( blank != NULL )
2692 *blank = '\0';
2693
2694 to_meter = atof(to_meter_str);
2695
2696 if( fabs(to_meter-1.0) < 0.0000001 )
2697 return MS_METERS;
2698 else if( fabs(to_meter-1000.0) < 0.00001 )
2699 return MS_KILOMETERS;
2700 else if( fabs(to_meter-0.3048) < 0.0001 )
2701 return MS_FEET;
2702 else if( fabs(to_meter-0.0254) < 0.0001 )
2703 return MS_INCHES;
2704 else if( fabs(to_meter-1609.344) < 0.001 )
2705 return MS_MILES;
2706 else if( fabs(to_meter-1852.0) < 0.1 )
2707 return MS_NAUTICALMILES;
2708 else
2709 return -1;
2710 }
2711
2712 #if PROJ_VERSION_MAJOR < 6
2713 pj_dalloc( proj_str );
2714 #endif
2715 return -1;
2716 }
2717
2718 /************************************************************************/
2719 /* msProjectionContextGetFromPool() */
2720 /* */
2721 /* Returns a projection context from the pool, or create a new */
2722 /* one if the pool is empty. */
2723 /* After use, it should normally be returned with */
2724 /* msProjectionContextReleaseToPool() */
2725 /************************************************************************/
2726
2727 projectionContext* msProjectionContextGetFromPool()
2728 {
2729 projectionContext* context;
2730 msAcquireLock( TLOCK_PROJ );
2731
2732 if( headOfLinkedListOfProjContext )
2733 {
2734 LinkedListOfProjContext* next = headOfLinkedListOfProjContext->next;
2735 context = headOfLinkedListOfProjContext->context;
2736 msFree(headOfLinkedListOfProjContext);
2737 headOfLinkedListOfProjContext = next;
2738 }
2739 else
2740 {
2741 context = msProjectionContextCreate();
2742 }
2743
2744 msReleaseLock( TLOCK_PROJ );
2745 return context;
2746 }
2747
2748 /************************************************************************/
2749 /* msProjectionContextReleaseToPool() */
2750 /************************************************************************/
2751
2752 void msProjectionContextReleaseToPool(projectionContext* ctx)
2753 {
2754 LinkedListOfProjContext* link =
2755 (LinkedListOfProjContext*)msSmallMalloc(sizeof(LinkedListOfProjContext));
2756 link->context = ctx;
2757 msAcquireLock( TLOCK_PROJ );
2758 link->next = headOfLinkedListOfProjContext;
2759 headOfLinkedListOfProjContext = link;
2760 msReleaseLock( TLOCK_PROJ );
2761 }
2762
2763 /************************************************************************/
2764 /* msProjectionContextPoolCleanup() */
2765 /************************************************************************/
2766
2767 void msProjectionContextPoolCleanup()
2768 {
2769 LinkedListOfProjContext* link;
2770 msAcquireLock( TLOCK_PROJ );
2771 link = headOfLinkedListOfProjContext;
2772 while( link )
2773 {
2774 LinkedListOfProjContext* next = link->next;
2775 msProjectionContextUnref(link->context);
2776 msFree(link);
2777 link = next;
2778 }
2779 headOfLinkedListOfProjContext = NULL;
2780 msReleaseLock( TLOCK_PROJ );
2781 }
2782