PostGIS  2.1.10dev-r@@SVN_REVISION@@
lwalgorithm.c
Go to the documentation of this file.
1 /**********************************************************************
2  * $Id: lwalgorithm.c 14427 2015-11-25 20:12:55Z pramsey $
3  *
4  * PostGIS - Spatial Types for PostgreSQL
5  * http://postgis.net
6  * Copyright 2008 Paul Ramsey
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include "liblwgeom_internal.h"
14 #include "lwgeom_log.h"
15 #include <ctype.h> /* for tolower */
16 
17 
21 int signum(double n)
22 {
23  if( n < 0 ) return -1;
24  if( n > 0 ) return 1;
25  return 0;
26 }
27 
28 int
29 p4d_same(const POINT4D *p1, const POINT4D *p2)
30 {
31  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) )
32  return LW_TRUE;
33  else
34  return LW_FALSE;
35 }
36 
37 int
38 p3d_same(const POINT3D *p1, const POINT3D *p2)
39 {
40  if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) && FP_EQUALS(p1->z,p2->z) )
41  return LW_TRUE;
42  else
43  return LW_FALSE;
44 }
45 
46 int
47 p2d_same(const POINT2D *p1, const POINT2D *p2)
48 {
49  if( FP_EQUALS(p1->x,p2->x) && FP_EQUALS(p1->y,p2->y) )
50  return LW_TRUE;
51  else
52  return LW_FALSE;
53 }
54 
62 int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
63 {
64  double side = ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
65  if ( side == 0.0 )
66  return 0;
67  else
68  return signum(side);
69 }
70 
74 double
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 
85 int
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 
95 int
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 
105 int
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 
118 double
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 
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 
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  int 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  int i = 0, j = 0;
464  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 
471  pa1 = (POINTARRAY*)l1->points;
472  pa2 = (POINTARRAY*)l2->points;
473 
474  /* One-point lines can't intersect (and shouldn't exist). */
475  if ( pa1->npoints < 2 || pa2->npoints < 2 )
476  return LINE_NO_CROSS;
477 
478  LWDEBUGF(4, "l1 = %s", lwgeom_to_ewkt((LWGEOM*)l1));
479  LWDEBUGF(4, "l2 = %s", lwgeom_to_ewkt((LWGEOM*)l2));
480 
481  /* Initialize first point of q */
482  getPoint2d_p(pa2, 0, &q1);
483 
484  for ( i = 1; i < pa2->npoints; i++ )
485  {
486 
487  /* Update second point of q to next value */
488  getPoint2d_p(pa2, i, &q2);
489 
490  /* Initialize first point of p */
491  getPoint2d_p(pa1, 0, &p1);
492 
493  for ( j = 1; j < pa1->npoints; j++ )
494  {
495 
496  /* Update second point of p to next value */
497  getPoint2d_p(pa1, j, &p2);
498 
499  this_cross = lw_segment_intersects(&p1, &p2, &q1, &q2);
500 
501  LWDEBUGF(4, "i=%d, j=%d (%.8g %.8g, %.8g %.8g)", this_cross, i, j, p1.x, p1.y, p2.x, p2.y);
502 
503  if ( this_cross == SEG_CROSS_LEFT )
504  {
505  LWDEBUG(4,"this_cross == SEG_CROSS_LEFT");
506  cross_left++;
507  if ( ! first_cross )
508  first_cross = SEG_CROSS_LEFT;
509  }
510 
511  if ( this_cross == SEG_CROSS_RIGHT )
512  {
513  LWDEBUG(4,"this_cross == SEG_CROSS_RIGHT");
514  cross_right++;
515  if ( ! first_cross )
516  first_cross = SEG_CROSS_LEFT;
517  }
518 
519  /*
520  ** Crossing at a co-linearity can be turned handled by extending
521  ** segment to next vertext and seeing if the end points straddle
522  ** the co-linear segment.
523  */
524  if ( this_cross == SEG_COLINEAR )
525  {
526  LWDEBUG(4,"this_cross == SEG_COLINEAR");
527  /* TODO: Add logic here and in segment_intersects()
528  continue;
529  */
530  }
531 
532  LWDEBUG(4,"this_cross == SEG_NO_INTERSECTION");
533 
534  /* Turn second point of p into first point */
535  p1 = p2;
536 
537  }
538 
539  /* Turn second point of q into first point */
540  q1 = q2;
541 
542  }
543 
544  LWDEBUGF(4, "first_cross=%d, cross_left=%d, cross_right=%d", first_cross, cross_left, cross_right);
545 
546  if ( !cross_left && !cross_right )
547  return LINE_NO_CROSS;
548 
549  if ( !cross_left && cross_right == 1 )
550  return LINE_CROSS_RIGHT;
551 
552  if ( !cross_right && cross_left == 1 )
553  return LINE_CROSS_LEFT;
554 
555  if ( cross_left - cross_right == 1 )
557 
558  if ( cross_left - cross_right == -1 )
560 
561  if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_LEFT )
563 
564  if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_RIGHT )
566 
567  return LINE_NO_CROSS;
568 
569 }
570 
571 
572 
573 
574 
575 static char *base32 = "0123456789bcdefghjkmnpqrstuvwxyz";
576 
577 /*
578 ** Calculate the geohash, iterating downwards and gaining precision.
579 ** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>
580 ** Released under the MIT License.
581 */
582 char *geohash_point(double longitude, double latitude, int precision)
583 {
584  int is_even=1, i=0;
585  double lat[2], lon[2], mid;
586  char bits[] = {16,8,4,2,1};
587  int bit=0, ch=0;
588  char *geohash = NULL;
589 
590  geohash = lwalloc(precision + 1);
591 
592  lat[0] = -90.0;
593  lat[1] = 90.0;
594  lon[0] = -180.0;
595  lon[1] = 180.0;
596 
597  while (i < precision)
598  {
599  if (is_even)
600  {
601  mid = (lon[0] + lon[1]) / 2;
602  if (longitude >= mid)
603  {
604  ch |= bits[bit];
605  lon[0] = mid;
606  }
607  else
608  {
609  lon[1] = mid;
610  }
611  }
612  else
613  {
614  mid = (lat[0] + lat[1]) / 2;
615  if (latitude >= mid)
616  {
617  ch |= bits[bit];
618  lat[0] = mid;
619  }
620  else
621  {
622  lat[1] = mid;
623  }
624  }
625 
626  is_even = !is_even;
627  if (bit < 4)
628  {
629  bit++;
630  }
631  else
632  {
633  geohash[i++] = base32[ch];
634  bit = 0;
635  ch = 0;
636  }
637  }
638  geohash[i] = 0;
639  return geohash;
640 }
641 
642 
643 /*
644 ** Calculate the geohash, iterating downwards and gaining precision.
645 ** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>
646 ** Released under the MIT License.
647 */
648 unsigned int geohash_point_as_int(POINT2D *pt)
649 {
650  int is_even=1;
651  double lat[2], lon[2], mid;
652  int bit=32;
653  unsigned int ch = 0;
654 
655  double longitude = pt->x;
656  double latitude = pt->y;
657 
658  lat[0] = -90.0;
659  lat[1] = 90.0;
660  lon[0] = -180.0;
661  lon[1] = 180.0;
662 
663  while (--bit >= 0)
664  {
665  if (is_even)
666  {
667  mid = (lon[0] + lon[1]) / 2;
668  if (longitude > mid)
669  {
670  ch |= 0x0001 << bit;
671  lon[0] = mid;
672  }
673  else
674  {
675  lon[1] = mid;
676  }
677  }
678  else
679  {
680  mid = (lat[0] + lat[1]) / 2;
681  if (latitude > mid)
682  {
683  ch |= 0x0001 << bit;
684  lat[0] = mid;
685  }
686  else
687  {
688  lat[1] = mid;
689  }
690  }
691 
692  is_even = !is_even;
693  }
694  return ch;
695 }
696 
697 /*
698 ** Decode a GeoHash into a bounding box. The lat and lon arguments should
699 ** both be passed as double arrays of length 2 at a minimum where the values
700 ** set in them will be the southwest and northeast coordinates of the bounding
701 ** box accordingly. A precision less than 0 indicates that the entire length
702 ** of the GeoHash should be used.
703 */
704 void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision)
705 {
706  int i, j, hashlen;
707  char c, cd, mask, is_even = 1;
708  static char bits[] = {16, 8, 4, 2, 1};
709 
710  lat[0] = -90.0;
711  lat[1] = 90.0;
712  lon[0] = -180.0;
713  lon[1] = 180.0;
714 
715  hashlen = strlen(geohash);
716 
717  if (precision < 0 || precision > hashlen)
718  {
719  precision = hashlen;
720  }
721 
722  for (i = 0; i < precision; i++)
723  {
724  c = tolower(geohash[i]);
725  cd = strchr(base32, c) - base32;
726 
727  for (j = 0; j < 5; j++)
728  {
729  mask = bits[j];
730  if (is_even)
731  {
732  lon[!(cd & mask)] = (lon[0] + lon[1]) / 2;
733  }
734  else
735  {
736  lat[!(cd & mask)] = (lat[0] + lat[1]) / 2;
737  }
738  is_even = !is_even;
739  }
740  }
741 }
742 
744 {
745  double minx, miny, maxx, maxy;
746  double latmax, latmin, lonmax, lonmin;
747  double lonwidth, latwidth;
748  double latmaxadjust, lonmaxadjust, latminadjust, lonminadjust;
749  int precision = 0;
750 
751  /* Get the bounding box, return error if things don't work out. */
752  minx = bbox.xmin;
753  miny = bbox.ymin;
754  maxx = bbox.xmax;
755  maxy = bbox.ymax;
756 
757  if ( minx == maxx && miny == maxy )
758  {
759  /* It's a point. Doubles have 51 bits of precision.
760  ** 2 * 51 / 5 == 20 */
761  return 20;
762  }
763 
764  lonmin = -180.0;
765  latmin = -90.0;
766  lonmax = 180.0;
767  latmax = 90.0;
768 
769  /* Shrink a world bounding box until one of the edges interferes with the
770  ** bounds of our rectangle. */
771  while ( 1 )
772  {
773  lonwidth = lonmax - lonmin;
774  latwidth = latmax - latmin;
775  latmaxadjust = lonmaxadjust = latminadjust = lonminadjust = 0.0;
776 
777  if ( minx > lonmin + lonwidth / 2.0 )
778  {
779  lonminadjust = lonwidth / 2.0;
780  }
781  else if ( maxx < lonmax - lonwidth / 2.0 )
782  {
783  lonmaxadjust = -1 * lonwidth / 2.0;
784  }
785  if ( miny > latmin + latwidth / 2.0 )
786  {
787  latminadjust = latwidth / 2.0;
788  }
789  else if (maxy < latmax - latwidth / 2.0 )
790  {
791  latmaxadjust = -1 * latwidth / 2.0;
792  }
793  /* Only adjust if adjustments are legal (we haven't crossed any edges). */
794  if ( (lonminadjust || lonmaxadjust) && (latminadjust || latmaxadjust ) )
795  {
796  latmin += latminadjust;
797  lonmin += lonminadjust;
798  latmax += latmaxadjust;
799  lonmax += lonmaxadjust;
800  /* Each adjustment cycle corresponds to 2 bits of storage in the
801  ** geohash. */
802  precision += 2;
803  }
804  else
805  {
806  break;
807  }
808  }
809 
810  /* Save the edges of our bounds, in case someone cares later. */
811  bounds->xmin = lonmin;
812  bounds->xmax = lonmax;
813  bounds->ymin = latmin;
814  bounds->ymax = latmax;
815 
816  /* Each geohash character (base32) can contain 5 bits of information.
817  ** We are returning the precision in characters, so here we divide. */
818  return precision / 5;
819 }
820 
821 
822 /*
823 ** Return a geohash string for the geometry. <http://geohash.org>
824 ** Where the precision is non-positive, calculate a precision based on the
825 ** bounds of the feature. Big features have loose precision.
826 ** Small features have tight precision.
827 */
828 char *lwgeom_geohash(const LWGEOM *lwgeom, int precision)
829 {
830  GBOX gbox;
831  GBOX gbox_bounds;
832  double lat, lon;
833  int result;
834 
835  gbox_init(&gbox);
836  gbox_init(&gbox_bounds);
837 
838  result = lwgeom_calculate_gbox_cartesian(lwgeom, &gbox);
839  if ( result == LW_FAILURE ) return NULL;
840 
841  /* Return error if we are being fed something outside our working bounds */
842  if ( gbox.xmin < -180 || gbox.ymin < -90 || gbox.xmax > 180 || gbox.ymax > 90 )
843  {
844  lwerror("Geohash requires inputs in decimal degrees.");
845  return NULL;
846  }
847 
848  /* What is the center of our geometry bounds? We'll use that to
849  ** approximate location. */
850  lon = gbox.xmin + (gbox.xmax - gbox.xmin) / 2;
851  lat = gbox.ymin + (gbox.ymax - gbox.ymin) / 2;
852 
853  if ( precision <= 0 )
854  {
855  precision = lwgeom_geohash_precision(gbox, &gbox_bounds);
856  }
857 
858  /*
859  ** Return the geohash of the center, with a precision determined by the
860  ** extent of the bounds.
861  ** Possible change: return the point at the center of the precision bounds?
862  */
863  return geohash_point(lon, lat, precision);
864 }
865 
866 
867 
868 
869 
870 
871 
872 
873 
874 
875 
876 
877 
878 
879 
880 
881 
882 
883 
884 
885 
886 
887 
double x
Definition: liblwgeom.h:308
double m
Definition: liblwgeom.h:308
double y
Definition: liblwgeom.h:296
int npoints
Definition: liblwgeom.h:327
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:425
double xmax
Definition: liblwgeom.h:249
int lwgeom_geohash_precision(GBOX bbox, GBOX *bounds)
Definition: lwalgorithm.c:743
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:106
#define LWDEBUG(level, msg)
Definition: lwgeom_log.h:50
double x
Definition: liblwgeom.h:296
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:86
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
The old function nessecary for ptarray_segmentize2d in ptarray.c.
Definition: measures.c:2123
#define FP_MIN(A, B)
char ** result
Definition: liblwgeom.h:218
#define LW_FAILURE
Definition: liblwgeom.h:54
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:296
char * geohash_point(double longitude, double latitude, int precision)
Definition: lwalgorithm.c:582
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:67
double x
Definition: liblwgeom.h:284
double ymin
Definition: liblwgeom.h:250
int pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring)
Definition: lwalgorithm.c:281
int signum(double n)
Returns -1 if n < 0.0 and 1 if n > 0.0.
Definition: lwalgorithm.c:21
void gbox_init(GBOX *gbox)
Zero out all the entries in the GBOX.
Definition: g_box.c:34
double xmin
Definition: liblwgeom.h:248
#define LW_FALSE
Definition: liblwgeom.h:52
const POINT2D * getPoint2d_cp(const POINTARRAY *pa, int n)
Returns a POINT2D pointer into the POINTARRAY serialized_ptlist, suitable for reading from...
Definition: lwgeom_api.c:458
#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:51
double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
Returns the length of a circular arc segment.
Definition: lwalgorithm.c:119
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:96
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:251
double y
Definition: liblwgeom.h:284
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:434
double z
Definition: liblwgeom.h:308
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:179
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:29
unsigned int geohash_point_as_int(POINT2D *pt)
Definition: lwalgorithm.c:648
int p2d_same(const POINT2D *p1, const POINT2D *p2)
Definition: lwalgorithm.c:47
char * lwgeom_geohash(const LWGEOM *lwgeom, int precision)
Calculate the GeoHash (http://geohash.org) string for a geometry.
Definition: lwalgorithm.c:828
#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:75
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:62
static char * base32
Definition: lwalgorithm.c:575
void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision)
Definition: lwalgorithm.c:704
void * lwalloc(size_t size)
Definition: lwutil.c:175
double y
Definition: liblwgeom.h:308
int p3d_same(const POINT3D *p1, const POINT3D *p2)
Definition: lwalgorithm.c:38
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:55
#define FP_MAX(A, B)
POINTARRAY * points
Definition: liblwgeom.h:378
int lwgeom_calculate_gbox_cartesian(const LWGEOM *lwgeom, GBOX *gbox)
Calculate the 2-4D bounding box of a geometry.
Definition: g_box.c:607