PostGIS  2.5.0beta2dev-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 
227 double
228 lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
229 {
230  POINT2D c;
231  double cx, cy, cr;
232  double dx21, dy21, dx31, dy31, h21, h31, d;
233 
234  c.x = c.y = 0.0;
235 
236  LWDEBUGF(2, "lw_arc_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);
237 
238  /* Closed circle */
239  if (fabs(p1->x - p3->x) < EPSILON_SQLMM &&
240  fabs(p1->y - p3->y) < EPSILON_SQLMM)
241  {
242  cx = p1->x + (p2->x - p1->x) / 2.0;
243  cy = p1->y + (p2->y - p1->y) / 2.0;
244  c.x = cx;
245  c.y = cy;
246  *result = c;
247  cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0));
248  return cr;
249  }
250 
251  /* Using cartesian eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle */
252  dx21 = p2->x - p1->x;
253  dy21 = p2->y - p1->y;
254  dx31 = p3->x - p1->x;
255  dy31 = p3->y - p1->y;
256 
257  h21 = pow(dx21, 2.0) + pow(dy21, 2.0);
258  h31 = pow(dx31, 2.0) + pow(dy31, 2.0);
259 
260  /* 2 * |Cross product|, d<0 means clockwise and d>0 counterclockwise sweeping angle */
261  d = 2 * (dx21 * dy31 - dx31 * dy21);
262 
263  /* Check colinearity, |Cross product| = 0 */
264  if (fabs(d) < EPSILON_SQLMM)
265  return -1.0;
266 
267  /* Calculate centroid coordinates and radius */
268  cx = p1->x + (h21 * dy31 - h31 * dy21) / d;
269  cy = p1->y - (h21 * dx31 - h31 * dx21) / d;
270  c.x = cx;
271  c.y = cy;
272  *result = c;
273  cr = sqrt(pow(cx - p1->x, 2) + pow(cy - p1->y, 2));
274 
275  LWDEBUGF(2, "lw_arc_center center is (%.16f,%.16f)", result->x, result->y);
276 
277  return cr;
278 }
279 
280 int
281 pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring)
282 {
283  int cn = 0; /* the crossing number counter */
284  uint32_t i;
285  const POINT2D *v1, *v2;
286  const POINT2D *first, *last;
287 
288  first = getPoint2d_cp(ring, 0);
289  last = getPoint2d_cp(ring, ring->npoints-1);
290  if ( memcmp(first, last, sizeof(POINT2D)) )
291  {
292  lwerror("pt_in_ring_2d: V[n] != V[0] (%g %g != %g %g)",
293  first->x, first->y, last->x, last->y);
294  return LW_FALSE;
295 
296  }
297 
298  LWDEBUGF(2, "pt_in_ring_2d called with point: %g %g", p->x, p->y);
299  /* printPA(ring); */
300 
301  /* loop through all edges of the polygon */
302  v1 = getPoint2d_cp(ring, 0);
303  for (i=0; i<ring->npoints-1; i++)
304  {
305  double vt;
306  v2 = getPoint2d_cp(ring, i+1);
307 
308  /* edge from vertex i to vertex i+1 */
309  if
310  (
311  /* an upward crossing */
312  ((v1->y <= p->y) && (v2->y > p->y))
313  /* a downward crossing */
314  || ((v1->y > p->y) && (v2->y <= p->y))
315  )
316  {
317 
318  vt = (double)(p->y - v1->y) / (v2->y - v1->y);
319 
320  /* P->x <intersect */
321  if (p->x < v1->x + vt * (v2->x - v1->x))
322  {
323  /* a valid crossing of y=p->y right of p->x */
324  ++cn;
325  }
326  }
327  v1 = v2;
328  }
329 
330  LWDEBUGF(3, "pt_in_ring_2d returning %d", cn&1);
331 
332  return (cn&1); /* 0 if even (out), and 1 if odd (in) */
333 }
334 
335 
336 static int
337 lw_seg_interact(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
338 {
339  double minq=FP_MIN(q1->x,q2->x);
340  double maxq=FP_MAX(q1->x,q2->x);
341  double minp=FP_MIN(p1->x,p2->x);
342  double maxp=FP_MAX(p1->x,p2->x);
343 
344  if (FP_GT(minp,maxq) || FP_LT(maxp,minq))
345  return LW_FALSE;
346 
347  minq=FP_MIN(q1->y,q2->y);
348  maxq=FP_MAX(q1->y,q2->y);
349  minp=FP_MIN(p1->y,p2->y);
350  maxp=FP_MAX(p1->y,p2->y);
351 
352  if (FP_GT(minp,maxq) || FP_LT(maxp,minq))
353  return LW_FALSE;
354 
355  return LW_TRUE;
356 }
357 
372 int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
373 {
374 
375  int pq1, pq2, qp1, qp2;
376 
377  /* No envelope interaction => we are done. */
378  if (!lw_seg_interact(p1, p2, q1, p2))
379  {
380  return SEG_NO_INTERSECTION;
381  }
382 
383  /* Are the start and end points of q on the same side of p? */
384  pq1=lw_segment_side(p1,p2,q1);
385  pq2=lw_segment_side(p1,p2,q2);
386  if ((pq1>0 && pq2>0) || (pq1<0 && pq2<0))
387  {
388  return SEG_NO_INTERSECTION;
389  }
390 
391  /* Are the start and end points of p on the same side of q? */
392  qp1=lw_segment_side(q1,q2,p1);
393  qp2=lw_segment_side(q1,q2,p2);
394  if ( (qp1 > 0.0 && qp2 > 0.0) || (qp1 < 0.0 && qp2 < 0.0) )
395  {
396  return SEG_NO_INTERSECTION;
397  }
398 
399  /* Nobody is on one side or another? Must be colinear. */
400  if ( pq1 == 0.0 && pq2 == 0.0 && qp1 == 0.0 && qp2 == 0.0 )
401  {
402  return SEG_COLINEAR;
403  }
404 
405  /*
406  ** When one end-point touches, the sidedness is determined by the
407  ** location of the other end-point. Only touches by the first point
408  ** will be considered "real" to avoid double counting.
409  */
410  LWDEBUGF(4, "pq1=%.15g pq2=%.15g", pq1, pq2);
411  LWDEBUGF(4, "qp1=%.15g qp2=%.15g", qp1, qp2);
412 
413  /* Second point of p or q touches, it's not a crossing. */
414  if ( pq2 == 0 || qp2 == 0 )
415  {
416  return SEG_NO_INTERSECTION;
417  }
418 
419  /* First point of p touches, it's a "crossing". */
420  if ( pq1 == 0 )
421  {
422  if ( pq2 > 0 )
423  return SEG_CROSS_RIGHT;
424  else
425  return SEG_CROSS_LEFT;
426  }
427 
428  /* First point of q touches, it's a crossing. */
429  if ( qp1 == 0 )
430  {
431  if ( pq1 < pq2 )
432  return SEG_CROSS_RIGHT;
433  else
434  return SEG_CROSS_LEFT;
435  }
436 
437  /* The segments cross, what direction is the crossing? */
438  if ( pq1 < pq2 )
439  return SEG_CROSS_RIGHT;
440  else
441  return SEG_CROSS_LEFT;
442 
443  /* This should never happen! */
444  return SEG_ERROR;
445 }
446 
461 int lwline_crossing_direction(const LWLINE *l1, const LWLINE *l2)
462 {
463  uint32_t i = 0, j = 0;
464  const POINT2D *p1, *p2, *q1, *q2;
465  POINTARRAY *pa1 = NULL, *pa2 = NULL;
466  int cross_left = 0;
467  int cross_right = 0;
468  int first_cross = 0;
469  int this_cross = 0;
470 #if POSTGIS_DEBUG_LEVEL >= 4
471  char *geom_ewkt;
472 #endif
473 
474  pa1 = (POINTARRAY*)l1->points;
475  pa2 = (POINTARRAY*)l2->points;
476 
477  /* One-point lines can't intersect (and shouldn't exist). */
478  if ( pa1->npoints < 2 || pa2->npoints < 2 )
479  return LINE_NO_CROSS;
480 
481 #if POSTGIS_DEBUG_LEVEL >= 4
482  geom_ewkt = lwgeom_to_ewkt((LWGEOM*)l1);
483  LWDEBUGF(4, "l1 = %s", geom_ewkt);
484  lwfree(geom_ewkt);
485  geom_ewkt = lwgeom_to_ewkt((LWGEOM*)l2);
486  LWDEBUGF(4, "l2 = %s", geom_ewkt);
487  lwfree(geom_ewkt);
488 #endif
489 
490  /* Initialize first point of q */
491  q1 = getPoint2d_cp(pa2, 0);
492 
493  for ( i = 1; i < pa2->npoints; i++ )
494  {
495 
496  /* Update second point of q to next value */
497  q2 = getPoint2d_cp(pa2, i);
498 
499  /* Initialize first point of p */
500  p1 = getPoint2d_cp(pa1, 0);
501 
502  for ( j = 1; j < pa1->npoints; j++ )
503  {
504 
505  /* Update second point of p to next value */
506  p2 = getPoint2d_cp(pa1, j);
507 
508  this_cross = lw_segment_intersects(p1, p2, q1, q2);
509 
510  LWDEBUGF(4, "i=%d, j=%d (%.8g %.8g, %.8g %.8g)", this_cross, i, j, p1->x, p1->y, p2->x, p2->y);
511 
512  if ( this_cross == SEG_CROSS_LEFT )
513  {
514  LWDEBUG(4,"this_cross == SEG_CROSS_LEFT");
515  cross_left++;
516  if ( ! first_cross )
517  first_cross = SEG_CROSS_LEFT;
518  }
519 
520  if ( this_cross == SEG_CROSS_RIGHT )
521  {
522  LWDEBUG(4,"this_cross == SEG_CROSS_RIGHT");
523  cross_right++;
524  if ( ! first_cross )
525  first_cross = SEG_CROSS_LEFT;
526  }
527 
528  /*
529  ** Crossing at a co-linearity can be turned handled by extending
530  ** segment to next vertex and seeing if the end points straddle
531  ** the co-linear segment.
532  */
533  if ( this_cross == SEG_COLINEAR )
534  {
535  LWDEBUG(4,"this_cross == SEG_COLINEAR");
536  /* TODO: Add logic here and in segment_intersects()
537  continue;
538  */
539  }
540 
541  LWDEBUG(4,"this_cross == SEG_NO_INTERSECTION");
542 
543  /* Turn second point of p into first point */
544  p1 = p2;
545 
546  }
547 
548  /* Turn second point of q into first point */
549  q1 = q2;
550 
551  }
552 
553  LWDEBUGF(4, "first_cross=%d, cross_left=%d, cross_right=%d", first_cross, cross_left, cross_right);
554 
555  if ( !cross_left && !cross_right )
556  return LINE_NO_CROSS;
557 
558  if ( !cross_left && cross_right == 1 )
559  return LINE_CROSS_RIGHT;
560 
561  if ( !cross_right && cross_left == 1 )
562  return LINE_CROSS_LEFT;
563 
564  if ( cross_left - cross_right == 1 )
566 
567  if ( cross_left - cross_right == -1 )
569 
570  if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_LEFT )
572 
573  if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_RIGHT )
575 
576  return LINE_NO_CROSS;
577 
578 }
579 
580 
581 
582 
583 
584 static char *base32 = "0123456789bcdefghjkmnpqrstuvwxyz";
585 
586 /*
587 ** Calculate the geohash, iterating downwards and gaining precision.
588 ** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>
589 ** Released under the MIT License.
590 */
591 char *geohash_point(double longitude, double latitude, int precision)
592 {
593  int is_even=1, i=0;
594  double lat[2], lon[2], mid;
595  char bits[] = {16,8,4,2,1};
596  int bit=0, ch=0;
597  char *geohash = NULL;
598 
599  geohash = lwalloc(precision + 1);
600 
601  lat[0] = -90.0;
602  lat[1] = 90.0;
603  lon[0] = -180.0;
604  lon[1] = 180.0;
605 
606  while (i < precision)
607  {
608  if (is_even)
609  {
610  mid = (lon[0] + lon[1]) / 2;
611  if (longitude >= mid)
612  {
613  ch |= bits[bit];
614  lon[0] = mid;
615  }
616  else
617  {
618  lon[1] = mid;
619  }
620  }
621  else
622  {
623  mid = (lat[0] + lat[1]) / 2;
624  if (latitude >= mid)
625  {
626  ch |= bits[bit];
627  lat[0] = mid;
628  }
629  else
630  {
631  lat[1] = mid;
632  }
633  }
634 
635  is_even = !is_even;
636  if (bit < 4)
637  {
638  bit++;
639  }
640  else
641  {
642  geohash[i++] = base32[ch];
643  bit = 0;
644  ch = 0;
645  }
646  }
647  geohash[i] = 0;
648  return geohash;
649 }
650 
651 
652 /*
653 ** Calculate the geohash, iterating downwards and gaining precision.
654 ** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>
655 ** Released under the MIT License.
656 */
657 unsigned int geohash_point_as_int(POINT2D *pt)
658 {
659  int is_even=1;
660  double lat[2], lon[2], mid;
661  int bit=32;
662  unsigned int ch = 0;
663 
664  double longitude = pt->x;
665  double latitude = pt->y;
666 
667  lat[0] = -90.0;
668  lat[1] = 90.0;
669  lon[0] = -180.0;
670  lon[1] = 180.0;
671 
672  while (--bit >= 0)
673  {
674  if (is_even)
675  {
676  mid = (lon[0] + lon[1]) / 2;
677  if (longitude > mid)
678  {
679  ch |= 0x0001u << bit;
680  lon[0] = mid;
681  }
682  else
683  {
684  lon[1] = mid;
685  }
686  }
687  else
688  {
689  mid = (lat[0] + lat[1]) / 2;
690  if (latitude > mid)
691  {
692  ch |= 0x0001 << bit;
693  lat[0] = mid;
694  }
695  else
696  {
697  lat[1] = mid;
698  }
699  }
700 
701  is_even = !is_even;
702  }
703  return ch;
704 }
705 
706 /*
707 ** Decode a GeoHash into a bounding box. The lat and lon arguments should
708 ** both be passed as double arrays of length 2 at a minimum where the values
709 ** set in them will be the southwest and northeast coordinates of the bounding
710 ** box accordingly. A precision less than 0 indicates that the entire length
711 ** of the GeoHash should be used.
712 */
713 void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision)
714 {
715  int i, j, hashlen;
716  char c, cd, mask, is_even = 1;
717  static char bits[] = {16, 8, 4, 2, 1};
718 
719  lat[0] = -90.0;
720  lat[1] = 90.0;
721  lon[0] = -180.0;
722  lon[1] = 180.0;
723 
724  hashlen = strlen(geohash);
725 
726  if (precision < 0 || precision > hashlen)
727  {
728  precision = hashlen;
729  }
730 
731  for (i = 0; i < precision; i++)
732  {
733  c = tolower(geohash[i]);
734  cd = strchr(base32, c) - base32;
735 
736  for (j = 0; j < 5; j++)
737  {
738  mask = bits[j];
739  if (is_even)
740  {
741  lon[!(cd & mask)] = (lon[0] + lon[1]) / 2;
742  }
743  else
744  {
745  lat[!(cd & mask)] = (lat[0] + lat[1]) / 2;
746  }
747  is_even = !is_even;
748  }
749  }
750 }
751 
753 {
754  double minx, miny, maxx, maxy;
755  double latmax, latmin, lonmax, lonmin;
756  double lonwidth, latwidth;
757  double latmaxadjust, lonmaxadjust, latminadjust, lonminadjust;
758  int precision = 0;
759 
760  /* Get the bounding box, return error if things don't work out. */
761  minx = bbox.xmin;
762  miny = bbox.ymin;
763  maxx = bbox.xmax;
764  maxy = bbox.ymax;
765 
766  if ( minx == maxx && miny == maxy )
767  {
768  /* It's a point. Doubles have 51 bits of precision.
769  ** 2 * 51 / 5 == 20 */
770  return 20;
771  }
772 
773  lonmin = -180.0;
774  latmin = -90.0;
775  lonmax = 180.0;
776  latmax = 90.0;
777 
778  /* Shrink a world bounding box until one of the edges interferes with the
779  ** bounds of our rectangle. */
780  while ( 1 )
781  {
782  lonwidth = lonmax - lonmin;
783  latwidth = latmax - latmin;
784  latmaxadjust = lonmaxadjust = latminadjust = lonminadjust = 0.0;
785 
786  if ( minx > lonmin + lonwidth / 2.0 )
787  {
788  lonminadjust = lonwidth / 2.0;
789  }
790  else if ( maxx < lonmax - lonwidth / 2.0 )
791  {
792  lonmaxadjust = -1 * lonwidth / 2.0;
793  }
794  if ( lonminadjust || lonmaxadjust )
795  {
796  lonmin += lonminadjust;
797  lonmax += lonmaxadjust;
798  /* Each adjustment cycle corresponds to 2 bits of storage in the
799  ** geohash. */
800  precision++;
801  }
802  else
803  {
804  break;
805  }
806 
807  if ( miny > latmin + latwidth / 2.0 )
808  {
809  latminadjust = latwidth / 2.0;
810  }
811  else if (maxy < latmax - latwidth / 2.0 )
812  {
813  latmaxadjust = -1 * latwidth / 2.0;
814  }
815  /* Only adjust if adjustments are legal (we haven't crossed any edges). */
816  if ( latminadjust || latmaxadjust )
817  {
818  latmin += latminadjust;
819  latmax += latmaxadjust;
820  /* Each adjustment cycle corresponds to 2 bits of storage in the
821  ** geohash. */
822  precision++;
823  }
824  else
825  {
826  break;
827  }
828  }
829 
830  /* Save the edges of our bounds, in case someone cares later. */
831  bounds->xmin = lonmin;
832  bounds->xmax = lonmax;
833  bounds->ymin = latmin;
834  bounds->ymax = latmax;
835 
836  /* Each geohash character (base32) can contain 5 bits of information.
837  ** We are returning the precision in characters, so here we divide. */
838  return precision / 5;
839 }
840 
841 
842 /*
843 ** Return a geohash string for the geometry. <http://geohash.org>
844 ** Where the precision is non-positive, calculate a precision based on the
845 ** bounds of the feature. Big features have loose precision.
846 ** Small features have tight precision.
847 */
848 char *lwgeom_geohash(const LWGEOM *lwgeom, int precision)
849 {
850  GBOX gbox;
851  GBOX gbox_bounds;
852  double lat, lon;
853  int result;
854 
855  gbox_init(&gbox);
856  gbox_init(&gbox_bounds);
857 
858  result = lwgeom_calculate_gbox_cartesian(lwgeom, &gbox);
859  if ( result == LW_FAILURE ) return NULL;
860 
861  /* Return error if we are being fed something outside our working bounds */
862  if ( gbox.xmin < -180 || gbox.ymin < -90 || gbox.xmax > 180 || gbox.ymax > 90 )
863  {
864  lwerror("Geohash requires inputs in decimal degrees, got (%g %g, %g %g).",
865  gbox.xmin, gbox.ymin,
866  gbox.xmax, gbox.ymax);
867  return NULL;
868  }
869 
870  /* What is the center of our geometry bounds? We'll use that to
871  ** approximate location. */
872  lon = gbox.xmin + (gbox.xmax - gbox.xmin) / 2;
873  lat = gbox.ymin + (gbox.ymax - gbox.ymin) / 2;
874 
875  if ( precision <= 0 )
876  {
877  precision = lwgeom_geohash_precision(gbox, &gbox_bounds);
878  }
879 
880  /*
881  ** Return the geohash of the center, with a precision determined by the
882  ** extent of the bounds.
883  ** Possible change: return the point at the center of the precision bounds?
884  */
885  return geohash_point(lon, lat, precision);
886 }
887 
888 
889 
890 
891 
892 
893 
894 
895 
896 
897 
898 
899 
900 
901 
902 
903 
904 
905 
906 
907 
908 
909 
double x
Definition: liblwgeom.h:354
double m
Definition: liblwgeom.h:354
void lwfree(void *mem)
Definition: lwutil.c:244
double y
Definition: liblwgeom.h:342
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:556
double xmax
Definition: liblwgeom.h:295
int lwgeom_geohash_precision(GBOX bbox, GBOX *bounds)
Definition: lwalgorithm.c:752
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:342
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)
Definition: measures.c:2313
#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:372
double z
Definition: liblwgeom.h:342
char * geohash_point(double longitude, double latitude, int precision)
Definition: lwalgorithm.c:591
unsigned int uint32_t
Definition: uthash.h:78
double x
Definition: liblwgeom.h:330
double ymin
Definition: liblwgeom.h:296
int pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring)
Definition: lwalgorithm.c:281
void gbox_init(GBOX *gbox)
Zero out all the entries in the GBOX.
Definition: g_box.c:47
double xmin
Definition: liblwgeom.h:294
#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:461
double ymax
Definition: liblwgeom.h:297
double y
Definition: liblwgeom.h:330
double z
Definition: liblwgeom.h:354
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:228
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:337
#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:657
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:848
#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:584
void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision)
Definition: lwalgorithm.c:713
void * lwalloc(size_t size)
Definition: lwutil.c:229
double y
Definition: liblwgeom.h:354
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:424
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:373