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