1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * PostGIS is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright 2008 Paul Ramsey
22 *
23 **********************************************************************/
24
25
26 #include "liblwgeom_internal.h"
27 #include "lwgeom_log.h"
28 #include <ctype.h> /* for tolower */
29 #include <stdbool.h>
30
31 int
p4d_same(const POINT4D * p1,const POINT4D * p2)32 p4d_same(const POINT4D *p1, const POINT4D *p2)
33 {
34 if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) && FP_EQUALS(p1->z,p2->z) && FP_EQUALS(p1->m,p2->m) )
35 return LW_TRUE;
36 else
37 return LW_FALSE;
38 }
39
40 int
p3d_same(const POINT3D * p1,const POINT3D * p2)41 p3d_same(const POINT3D *p1, const POINT3D *p2)
42 {
43 if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) && FP_EQUALS(p1->z,p2->z) )
44 return LW_TRUE;
45 else
46 return LW_FALSE;
47 }
48
49 int
p2d_same(const POINT2D * p1,const POINT2D * p2)50 p2d_same(const POINT2D *p1, const POINT2D *p2)
51 {
52 if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) )
53 return LW_TRUE;
54 else
55 return LW_FALSE;
56 }
57
58 /**
59 * lw_segment_side()
60 *
61 * Return -1 if point Q is left of segment P
62 * Return 1 if point Q is right of segment P
63 * Return 0 if point Q in on segment P
64 */
lw_segment_side(const POINT2D * p1,const POINT2D * p2,const POINT2D * q)65 int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
66 {
67 double side = ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
68 return SIGNUM(side);
69 }
70
71 /**
72 * Returns the length of a linear segment
73 */
74 double
lw_seg_length(const POINT2D * A1,const POINT2D * A2)75 lw_seg_length(const POINT2D *A1, const POINT2D *A2)
76 {
77 return sqrt((A1->x-A2->x)*(A1->x-A2->x)+(A1->y-A2->y)*(A1->y-A2->y));
78 }
79
80 /**
81 * Returns true if P is on the same side of the plane partition
82 * defined by A1/A3 as A2 is. Only makes sense if P has already been
83 * determined to be on the circle defined by A1/A2/A3.
84 */
85 int
lw_pt_in_arc(const POINT2D * P,const POINT2D * A1,const POINT2D * A2,const POINT2D * A3)86 lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
87 {
88 return lw_segment_side(A1, A3, A2) == lw_segment_side(A1, A3, P);
89 }
90
91 /**
92 * Returns true if P is between A1/A2. Only makes sense if P has already been
93 * deterined to be on the line defined by A1/A2.
94 */
95 int
lw_pt_in_seg(const POINT2D * P,const POINT2D * A1,const POINT2D * A2)96 lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2)
97 {
98 return ((A1->x <= P->x && P->x < A2->x) || (A1->x >= P->x && P->x > A2->x)) ||
99 ((A1->y <= P->y && P->y < A2->y) || (A1->y >= P->y && P->y > A2->y));
100 }
101
102 /**
103 * Returns true if arc A is actually a point (all vertices are the same) .
104 */
105 int
lw_arc_is_pt(const POINT2D * A1,const POINT2D * A2,const POINT2D * A3)106 lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
107 {
108 if ( A1->x == A2->x && A2->x == A3->x &&
109 A1->y == A2->y && A2->y == A3->y )
110 return LW_TRUE;
111 else
112 return LW_FALSE;
113 }
114
115 /**
116 * Returns the length of a circular arc segment
117 */
118 double
lw_arc_length(const POINT2D * A1,const POINT2D * A2,const POINT2D * A3)119 lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
120 {
121 POINT2D C;
122 double radius_A, circumference_A;
123 int a2_side, clockwise;
124 double a1, a3;
125 double angle;
126
127 if ( lw_arc_is_pt(A1, A2, A3) )
128 return 0.0;
129
130 radius_A = lw_arc_center(A1, A2, A3, &C);
131
132 /* Co-linear! Return linear distance! */
133 if ( radius_A < 0 )
134 {
135 double dx = A1->x - A3->x;
136 double dy = A1->y - A3->y;
137 return sqrt(dx*dx + dy*dy);
138 }
139
140 /* Closed circle! Return the circumference! */
141 circumference_A = M_PI * 2 * radius_A;
142 if ( p2d_same(A1, A3) )
143 return circumference_A;
144
145 /* Determine the orientation of the arc */
146 a2_side = lw_segment_side(A1, A3, A2);
147
148 /* The side of the A1/A3 line that A2 falls on dictates the sweep
149 direction from A1 to A3. */
150 if ( a2_side == -1 )
151 clockwise = LW_TRUE;
152 else
153 clockwise = LW_FALSE;
154
155 /* Angles of each point that defines the arc section */
156 a1 = atan2(A1->y - C.y, A1->x - C.x);
157 a3 = atan2(A3->y - C.y, A3->x - C.x);
158
159 /* What's the sweep from A1 to A3? */
160 if ( clockwise )
161 {
162 if ( a1 > a3 )
163 angle = a1 - a3;
164 else
165 angle = 2*M_PI + a1 - a3;
166 }
167 else
168 {
169 if ( a3 > a1 )
170 angle = a3 - a1;
171 else
172 angle = 2*M_PI + a3 - a1;
173 }
174
175 /* Length as proportion of circumference */
176 return circumference_A * (angle / (2*M_PI));
177 }
178
lw_arc_side(const POINT2D * A1,const POINT2D * A2,const POINT2D * A3,const POINT2D * Q)179 int lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q)
180 {
181 POINT2D C;
182 double radius_A;
183 double side_Q, side_A2;
184 double d;
185
186 side_Q = lw_segment_side(A1, A3, Q);
187 radius_A = lw_arc_center(A1, A2, A3, &C);
188 side_A2 = lw_segment_side(A1, A3, A2);
189
190 /* Linear case */
191 if ( radius_A < 0 )
192 return side_Q;
193
194 d = distance2d_pt_pt(Q, &C);
195
196 /* Q is on the arc boundary */
197 if ( d == radius_A && side_Q == side_A2 )
198 {
199 return 0;
200 }
201
202 /* Q on A1-A3 line, so its on opposite side to A2 */
203 if ( side_Q == 0 )
204 {
205 return -1 * side_A2;
206 }
207
208 /*
209 * Q is inside the arc boundary, so it's not on the side we
210 * might think from examining only the end points
211 */
212 if ( d < radius_A && side_Q == side_A2 )
213 {
214 side_Q *= -1;
215 }
216
217 return side_Q;
218 }
219
220 /**
221 * Determines the center of the circle defined by the three given points.
222 * In the event the circle is complete, the midpoint of the segment defined
223 * by the first and second points is returned. If the points are collinear,
224 * as determined by equal slopes, then -1.0 is returned. If the interior
225 * point is coincident with either end point, they are taken as collinear.
226 * For non-collinear cases, arc radious is returned.
227 */
228 double
lw_arc_center(const POINT2D * p1,const POINT2D * p2,const POINT2D * p3,POINT2D * result)229 lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
230 {
231 POINT2D c;
232 double cx, cy, cr;
233 double dx21, dy21, dx31, dy31, h21, h31, d;
234
235 c.x = c.y = 0.0;
236
237 LWDEBUGF(2, "lw_arc_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);
238
239 /* Closed circle */
240 if (fabs(p1->x - p3->x) < EPSILON_SQLMM &&
241 fabs(p1->y - p3->y) < EPSILON_SQLMM)
242 {
243 cx = p1->x + (p2->x - p1->x) / 2.0;
244 cy = p1->y + (p2->y - p1->y) / 2.0;
245 c.x = cx;
246 c.y = cy;
247 *result = c;
248 cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0));
249 return cr;
250 }
251
252 /* Using cartesian eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle */
253 dx21 = p2->x - p1->x;
254 dy21 = p2->y - p1->y;
255 dx31 = p3->x - p1->x;
256 dy31 = p3->y - p1->y;
257
258 h21 = pow(dx21, 2.0) + pow(dy21, 2.0);
259 h31 = pow(dx31, 2.0) + pow(dy31, 2.0);
260
261 /* 2 * |Cross product|, d<0 means clockwise and d>0 counterclockwise sweeping angle */
262 d = 2 * (dx21 * dy31 - dx31 * dy21);
263
264 /* Check colinearity, |Cross product| = 0 */
265 if (fabs(d) < EPSILON_SQLMM)
266 return -1.0;
267
268 /* Calculate centroid coordinates and radius */
269 cx = p1->x + (h21 * dy31 - h31 * dy21) / d;
270 cy = p1->y - (h21 * dx31 - h31 * dx21) / d;
271 c.x = cx;
272 c.y = cy;
273 *result = c;
274 cr = sqrt(pow(cx - p1->x, 2) + pow(cy - p1->y, 2));
275
276 LWDEBUGF(2, "lw_arc_center center is (%.16f,%.16f)", result->x, result->y);
277
278 return cr;
279 }
280
281 int
pt_in_ring_2d(const POINT2D * p,const POINTARRAY * ring)282 pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring)
283 {
284 int cn = 0; /* the crossing number counter */
285 uint32_t i;
286 const POINT2D *v1, *v2;
287 const POINT2D *first, *last;
288
289 first = getPoint2d_cp(ring, 0);
290 last = getPoint2d_cp(ring, ring->npoints-1);
291 if ( memcmp(first, last, sizeof(POINT2D)) )
292 {
293 lwerror("pt_in_ring_2d: V[n] != V[0] (%g %g != %g %g)",
294 first->x, first->y, last->x, last->y);
295 return LW_FALSE;
296
297 }
298
299 LWDEBUGF(2, "pt_in_ring_2d called with point: %g %g", p->x, p->y);
300 /* printPA(ring); */
301
302 /* loop through all edges of the polygon */
303 v1 = getPoint2d_cp(ring, 0);
304 for (i=0; i<ring->npoints-1; i++)
305 {
306 double vt;
307 v2 = getPoint2d_cp(ring, i+1);
308
309 /* edge from vertex i to vertex i+1 */
310 if
311 (
312 /* an upward crossing */
313 ((v1->y <= p->y) && (v2->y > p->y))
314 /* a downward crossing */
315 || ((v1->y > p->y) && (v2->y <= p->y))
316 )
317 {
318
319 vt = (double)(p->y - v1->y) / (v2->y - v1->y);
320
321 /* P->x <intersect */
322 if (p->x < v1->x + vt * (v2->x - v1->x))
323 {
324 /* a valid crossing of y=p->y right of p->x */
325 ++cn;
326 }
327 }
328 v1 = v2;
329 }
330
331 LWDEBUGF(3, "pt_in_ring_2d returning %d", cn&1);
332
333 return (cn&1); /* 0 if even (out), and 1 if odd (in) */
334 }
335
336
337 static int
lw_seg_interact(const POINT2D * p1,const POINT2D * p2,const POINT2D * q1,const POINT2D * q2)338 lw_seg_interact(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
339 {
340 double minq=FP_MIN(q1->x,q2->x);
341 double maxq=FP_MAX(q1->x,q2->x);
342 double minp=FP_MIN(p1->x,p2->x);
343 double maxp=FP_MAX(p1->x,p2->x);
344
345 if (FP_GT(minp,maxq) || FP_LT(maxp,minq))
346 return LW_FALSE;
347
348 minq=FP_MIN(q1->y,q2->y);
349 maxq=FP_MAX(q1->y,q2->y);
350 minp=FP_MIN(p1->y,p2->y);
351 maxp=FP_MAX(p1->y,p2->y);
352
353 if (FP_GT(minp,maxq) || FP_LT(maxp,minq))
354 return LW_FALSE;
355
356 return LW_TRUE;
357 }
358
359 /**
360 ** @brief returns the kind of #CG_SEGMENT_INTERSECTION_TYPE behavior of lineseg 1 (constructed from p1 and p2) and lineseg 2 (constructed from q1 and q2)
361 ** @param p1 start point of first straight linesegment
362 ** @param p2 end point of first straight linesegment
363 ** @param q1 start point of second line segment
364 ** @param q2 end point of second line segment
365 ** @return a #CG_SEGMENT_INTERSECTION_TYPE
366 ** Returns one of
367 ** SEG_ERROR = -1,
368 ** SEG_NO_INTERSECTION = 0,
369 ** SEG_COLINEAR = 1,
370 ** SEG_CROSS_LEFT = 2,
371 ** SEG_CROSS_RIGHT = 3,
372 */
lw_segment_intersects(const POINT2D * p1,const POINT2D * p2,const POINT2D * q1,const POINT2D * q2)373 int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
374 {
375
376 int pq1, pq2, qp1, qp2;
377
378 /* No envelope interaction => we are done. */
379 if (!lw_seg_interact(p1, p2, q1, p2))
380 {
381 return SEG_NO_INTERSECTION;
382 }
383
384 /* Are the start and end points of q on the same side of p? */
385 pq1=lw_segment_side(p1,p2,q1);
386 pq2=lw_segment_side(p1,p2,q2);
387 if ((pq1>0 && pq2>0) || (pq1<0 && pq2<0))
388 {
389 return SEG_NO_INTERSECTION;
390 }
391
392 /* Are the start and end points of p on the same side of q? */
393 qp1=lw_segment_side(q1,q2,p1);
394 qp2=lw_segment_side(q1,q2,p2);
395 if ( (qp1 > 0.0 && qp2 > 0.0) || (qp1 < 0.0 && qp2 < 0.0) )
396 {
397 return SEG_NO_INTERSECTION;
398 }
399
400 /* Nobody is on one side or another? Must be colinear. */
401 if ( pq1 == 0.0 && pq2 == 0.0 && qp1 == 0.0 && qp2 == 0.0 )
402 {
403 return SEG_COLINEAR;
404 }
405
406 /*
407 ** When one end-point touches, the sidedness is determined by the
408 ** location of the other end-point. Only touches by the first point
409 ** will be considered "real" to avoid double counting.
410 */
411 LWDEBUGF(4, "pq1=%.15g pq2=%.15g", pq1, pq2);
412 LWDEBUGF(4, "qp1=%.15g qp2=%.15g", qp1, qp2);
413
414 /* Second point of p or q touches, it's not a crossing. */
415 if ( pq2 == 0 || qp2 == 0 )
416 {
417 return SEG_NO_INTERSECTION;
418 }
419
420 /* First point of p touches, it's a "crossing". */
421 if ( pq1 == 0 )
422 {
423 if ( pq2 > 0 )
424 return SEG_CROSS_RIGHT;
425 else
426 return SEG_CROSS_LEFT;
427 }
428
429 /* First point of q touches, it's a crossing. */
430 if ( qp1 == 0 )
431 {
432 if ( pq1 < pq2 )
433 return SEG_CROSS_RIGHT;
434 else
435 return SEG_CROSS_LEFT;
436 }
437
438 /* The segments cross, what direction is the crossing? */
439 if ( pq1 < pq2 )
440 return SEG_CROSS_RIGHT;
441 else
442 return SEG_CROSS_LEFT;
443
444 /* This should never happen! */
445 return SEG_ERROR;
446 }
447
448 /**
449 ** @brief lwline_crossing_direction: returns the kind of #CG_LINE_CROSS_TYPE behavior of 2 linestrings
450 ** @param l1 first line string
451 ** @param l2 second line string
452 ** @return a #CG_LINE_CROSS_TYPE
453 ** LINE_NO_CROSS = 0
454 ** LINE_CROSS_LEFT = -1
455 ** LINE_CROSS_RIGHT = 1
456 ** LINE_MULTICROSS_END_LEFT = -2
457 ** LINE_MULTICROSS_END_RIGHT = 2
458 ** LINE_MULTICROSS_END_SAME_FIRST_LEFT = -3
459 ** LINE_MULTICROSS_END_SAME_FIRST_RIGHT = 3
460 **
461 */
lwline_crossing_direction(const LWLINE * l1,const LWLINE * l2)462 int lwline_crossing_direction(const LWLINE *l1, const LWLINE *l2)
463 {
464 uint32_t i = 0, j = 0;
465 const POINT2D *p1, *p2, *q1, *q2;
466 POINTARRAY *pa1 = NULL, *pa2 = NULL;
467 int cross_left = 0;
468 int cross_right = 0;
469 int first_cross = 0;
470 int this_cross = 0;
471 #if POSTGIS_DEBUG_LEVEL >= 4
472 char *geom_ewkt;
473 #endif
474
475 pa1 = (POINTARRAY*)l1->points;
476 pa2 = (POINTARRAY*)l2->points;
477
478 /* One-point lines can't intersect (and shouldn't exist). */
479 if ( pa1->npoints < 2 || pa2->npoints < 2 )
480 return LINE_NO_CROSS;
481
482 #if POSTGIS_DEBUG_LEVEL >= 4
483 geom_ewkt = lwgeom_to_ewkt((LWGEOM*)l1);
484 LWDEBUGF(4, "l1 = %s", geom_ewkt);
485 lwfree(geom_ewkt);
486 geom_ewkt = lwgeom_to_ewkt((LWGEOM*)l2);
487 LWDEBUGF(4, "l2 = %s", geom_ewkt);
488 lwfree(geom_ewkt);
489 #endif
490
491 /* Initialize first point of q */
492 q1 = getPoint2d_cp(pa2, 0);
493
494 for ( i = 1; i < pa2->npoints; i++ )
495 {
496
497 /* Update second point of q to next value */
498 q2 = getPoint2d_cp(pa2, i);
499
500 /* Initialize first point of p */
501 p1 = getPoint2d_cp(pa1, 0);
502
503 for ( j = 1; j < pa1->npoints; j++ )
504 {
505
506 /* Update second point of p to next value */
507 p2 = getPoint2d_cp(pa1, j);
508
509 this_cross = lw_segment_intersects(p1, p2, q1, q2);
510
511 LWDEBUGF(4, "i=%d, j=%d (%.8g %.8g, %.8g %.8g)", this_cross, i, j, p1->x, p1->y, p2->x, p2->y);
512
513 if ( this_cross == SEG_CROSS_LEFT )
514 {
515 LWDEBUG(4,"this_cross == SEG_CROSS_LEFT");
516 cross_left++;
517 if ( ! first_cross )
518 first_cross = SEG_CROSS_LEFT;
519 }
520
521 if ( this_cross == SEG_CROSS_RIGHT )
522 {
523 LWDEBUG(4,"this_cross == SEG_CROSS_RIGHT");
524 cross_right++;
525 if ( ! first_cross )
526 first_cross = SEG_CROSS_LEFT;
527 }
528
529 /*
530 ** Crossing at a co-linearity can be turned handled by extending
531 ** segment to next vertex and seeing if the end points straddle
532 ** the co-linear segment.
533 */
534 if ( this_cross == SEG_COLINEAR )
535 {
536 LWDEBUG(4,"this_cross == SEG_COLINEAR");
537 /* TODO: Add logic here and in segment_intersects()
538 continue;
539 */
540 }
541
542 LWDEBUG(4,"this_cross == SEG_NO_INTERSECTION");
543
544 /* Turn second point of p into first point */
545 p1 = p2;
546
547 }
548
549 /* Turn second point of q into first point */
550 q1 = q2;
551
552 }
553
554 LWDEBUGF(4, "first_cross=%d, cross_left=%d, cross_right=%d", first_cross, cross_left, cross_right);
555
556 if ( !cross_left && !cross_right )
557 return LINE_NO_CROSS;
558
559 if ( !cross_left && cross_right == 1 )
560 return LINE_CROSS_RIGHT;
561
562 if ( !cross_right && cross_left == 1 )
563 return LINE_CROSS_LEFT;
564
565 if ( cross_left - cross_right == 1 )
566 return LINE_MULTICROSS_END_LEFT;
567
568 if ( cross_left - cross_right == -1 )
569 return LINE_MULTICROSS_END_RIGHT;
570
571 if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_LEFT )
572 return LINE_MULTICROSS_END_SAME_FIRST_LEFT;
573
574 if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_RIGHT )
575 return LINE_MULTICROSS_END_SAME_FIRST_RIGHT;
576
577 return LINE_NO_CROSS;
578
579 }
580
581
582
583
584
585 static char *base32 = "0123456789bcdefghjkmnpqrstuvwxyz";
586
587 /*
588 ** Calculate the geohash, iterating downwards and gaining precision.
589 ** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>
590 ** Released under the MIT License.
591 */
geohash_point(double longitude,double latitude,int precision)592 char *geohash_point(double longitude, double latitude, int precision)
593 {
594 int is_even=1, i=0;
595 double lat[2], lon[2], mid;
596 char bits[] = {16,8,4,2,1};
597 int bit=0, ch=0;
598 char *geohash = NULL;
599
600 geohash = lwalloc(precision + 1);
601
602 lat[0] = -90.0;
603 lat[1] = 90.0;
604 lon[0] = -180.0;
605 lon[1] = 180.0;
606
607 while (i < precision)
608 {
609 if (is_even)
610 {
611 mid = (lon[0] + lon[1]) / 2;
612 if (longitude >= mid)
613 {
614 ch |= bits[bit];
615 lon[0] = mid;
616 }
617 else
618 {
619 lon[1] = mid;
620 }
621 }
622 else
623 {
624 mid = (lat[0] + lat[1]) / 2;
625 if (latitude >= mid)
626 {
627 ch |= bits[bit];
628 lat[0] = mid;
629 }
630 else
631 {
632 lat[1] = mid;
633 }
634 }
635
636 is_even = !is_even;
637 if (bit < 4)
638 {
639 bit++;
640 }
641 else
642 {
643 geohash[i++] = base32[ch];
644 bit = 0;
645 ch = 0;
646 }
647 }
648 geohash[i] = 0;
649 return geohash;
650 }
651
652
653 /*
654 ** Calculate the geohash, iterating downwards and gaining precision.
655 ** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>
656 ** Released under the MIT License.
657 */
geohash_point_as_int(POINT2D * pt)658 unsigned int geohash_point_as_int(POINT2D *pt)
659 {
660 int is_even=1;
661 double lat[2], lon[2], mid;
662 int bit=32;
663 unsigned int ch = 0;
664
665 double longitude = pt->x;
666 double latitude = pt->y;
667
668 lat[0] = -90.0;
669 lat[1] = 90.0;
670 lon[0] = -180.0;
671 lon[1] = 180.0;
672
673 while (--bit >= 0)
674 {
675 if (is_even)
676 {
677 mid = (lon[0] + lon[1]) / 2;
678 if (longitude > mid)
679 {
680 ch |= 0x0001u << bit;
681 lon[0] = mid;
682 }
683 else
684 {
685 lon[1] = mid;
686 }
687 }
688 else
689 {
690 mid = (lat[0] + lat[1]) / 2;
691 if (latitude > mid)
692 {
693 ch |= 0x0001 << bit;
694 lat[0] = mid;
695 }
696 else
697 {
698 lat[1] = mid;
699 }
700 }
701
702 is_even = !is_even;
703 }
704 return ch;
705 }
706
707 /*
708 ** Decode a GeoHash into a bounding box. The lat and lon arguments should
709 ** both be passed as double arrays of length 2 at a minimum where the values
710 ** set in them will be the southwest and northeast coordinates of the bounding
711 ** box accordingly. A precision less than 0 indicates that the entire length
712 ** of the GeoHash should be used.
713 ** It will call `lwerror` if an invalid character is found
714 */
decode_geohash_bbox(char * geohash,double * lat,double * lon,int precision)715 void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision)
716 {
717 bool is_even = 1;
718
719 lat[0] = -90.0;
720 lat[1] = 90.0;
721 lon[0] = -180.0;
722 lon[1] = 180.0;
723
724 size_t hashlen = strlen(geohash);
725 if (precision < 0 || (size_t)precision > hashlen)
726 {
727 precision = (int)hashlen;
728 }
729
730 for (int i = 0; i < precision; i++)
731 {
732 char c = tolower(geohash[i]);
733
734 /* Valid characters are all digits in base32 */
735 char *base32_pos = strchr(base32, c);
736 if (!base32_pos)
737 {
738 lwerror("%s: Invalid character '%c'", __func__, geohash[i]);
739 return;
740 }
741 char cd = base32_pos - base32;
742
743 for (size_t j = 0; j < 5; j++)
744 {
745 const char bits[] = {16, 8, 4, 2, 1};
746 char mask = bits[j];
747 if (is_even)
748 {
749 lon[!(cd & mask)] = (lon[0] + lon[1]) / 2;
750 }
751 else
752 {
753 lat[!(cd & mask)] = (lat[0] + lat[1]) / 2;
754 }
755 is_even = !is_even;
756 }
757 }
758 }
759
lwgeom_geohash_precision(GBOX bbox,GBOX * bounds)760 int lwgeom_geohash_precision(GBOX bbox, GBOX *bounds)
761 {
762 double minx, miny, maxx, maxy;
763 double latmax, latmin, lonmax, lonmin;
764 double lonwidth, latwidth;
765 double latmaxadjust, lonmaxadjust, latminadjust, lonminadjust;
766 int precision = 0;
767
768 /* Get the bounding box, return error if things don't work out. */
769 minx = bbox.xmin;
770 miny = bbox.ymin;
771 maxx = bbox.xmax;
772 maxy = bbox.ymax;
773
774 if ( minx == maxx && miny == maxy )
775 {
776 /* It's a point. Doubles have 51 bits of precision.
777 ** 2 * 51 / 5 == 20 */
778 return 20;
779 }
780
781 lonmin = -180.0;
782 latmin = -90.0;
783 lonmax = 180.0;
784 latmax = 90.0;
785
786 /* Shrink a world bounding box until one of the edges interferes with the
787 ** bounds of our rectangle. */
788 while ( 1 )
789 {
790 lonwidth = lonmax - lonmin;
791 latwidth = latmax - latmin;
792 latmaxadjust = lonmaxadjust = latminadjust = lonminadjust = 0.0;
793
794 if ( minx > lonmin + lonwidth / 2.0 )
795 {
796 lonminadjust = lonwidth / 2.0;
797 }
798 else if ( maxx < lonmax - lonwidth / 2.0 )
799 {
800 lonmaxadjust = -1 * lonwidth / 2.0;
801 }
802 if ( lonminadjust || lonmaxadjust )
803 {
804 lonmin += lonminadjust;
805 lonmax += lonmaxadjust;
806 /* Each adjustment cycle corresponds to 2 bits of storage in the
807 ** geohash. */
808 precision++;
809 }
810 else
811 {
812 break;
813 }
814
815 if ( miny > latmin + latwidth / 2.0 )
816 {
817 latminadjust = latwidth / 2.0;
818 }
819 else if (maxy < latmax - latwidth / 2.0 )
820 {
821 latmaxadjust = -1 * latwidth / 2.0;
822 }
823 /* Only adjust if adjustments are legal (we haven't crossed any edges). */
824 if ( latminadjust || latmaxadjust )
825 {
826 latmin += latminadjust;
827 latmax += latmaxadjust;
828 /* Each adjustment cycle corresponds to 2 bits of storage in the
829 ** geohash. */
830 precision++;
831 }
832 else
833 {
834 break;
835 }
836 }
837
838 /* Save the edges of our bounds, in case someone cares later. */
839 bounds->xmin = lonmin;
840 bounds->xmax = lonmax;
841 bounds->ymin = latmin;
842 bounds->ymax = latmax;
843
844 /* Each geohash character (base32) can contain 5 bits of information.
845 ** We are returning the precision in characters, so here we divide. */
846 return precision / 5;
847 }
848
849
850 /*
851 ** Return a geohash string for the geometry. <http://geohash.org>
852 ** Where the precision is non-positive, calculate a precision based on the
853 ** bounds of the feature. Big features have loose precision.
854 ** Small features have tight precision.
855 */
lwgeom_geohash(const LWGEOM * lwgeom,int precision)856 char *lwgeom_geohash(const LWGEOM *lwgeom, int precision)
857 {
858 GBOX gbox;
859 GBOX gbox_bounds;
860 double lat, lon;
861 int result;
862
863 gbox_init(&gbox);
864 gbox_init(&gbox_bounds);
865
866 result = lwgeom_calculate_gbox_cartesian(lwgeom, &gbox);
867 if ( result == LW_FAILURE ) return NULL;
868
869 /* Return error if we are being fed something outside our working bounds */
870 if ( gbox.xmin < -180 || gbox.ymin < -90 || gbox.xmax > 180 || gbox.ymax > 90 )
871 {
872 lwerror("Geohash requires inputs in decimal degrees, got (%g %g, %g %g).",
873 gbox.xmin, gbox.ymin,
874 gbox.xmax, gbox.ymax);
875 return NULL;
876 }
877
878 /* What is the center of our geometry bounds? We'll use that to
879 ** approximate location. */
880 lon = gbox.xmin + (gbox.xmax - gbox.xmin) / 2;
881 lat = gbox.ymin + (gbox.ymax - gbox.ymin) / 2;
882
883 if ( precision <= 0 )
884 {
885 precision = lwgeom_geohash_precision(gbox, &gbox_bounds);
886 }
887
888 /*
889 ** Return the geohash of the center, with a precision determined by the
890 ** extent of the bounds.
891 ** Possible change: return the point at the center of the precision bounds?
892 */
893 return geohash_point(lon, lat, precision);
894 }
895