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