1
/**********************************************************************
2
* $Id: lwalgorithm.c 4168 2009-06-11 16:44:03Z pramsey $
4
* PostGIS - Spatial Types for PostgreSQL
5
* http://postgis.refractions.net
6
* Copyright 2008 Paul Ramsey
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.
11
**********************************************************************/
13
#include "lwalgorithm.h"
19
** Return < 0.0 if point Q is left of segment P
20
** Return > 0.0 if point Q is right of segment P
21
** Return = 0.0 if point Q in on segment P
23
double lw_segment_side(POINT2D *p1, POINT2D *p2, POINT2D *q)
25
return ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
28
int lw_segment_envelope_intersects(POINT2D p1, POINT2D p2, POINT2D q1, POINT2D q2)
30
double minq=LW_MIN(q1.x,q2.x);
31
double maxq=LW_MAX(q1.x,q2.x);
32
double minp=LW_MIN(p1.x,p2.x);
33
double maxp=LW_MAX(p1.x,p2.x);
35
if (minp>maxq || maxp<minq)
38
minq=LW_MIN(q1.y,q2.y);
39
maxq=LW_MAX(q1.y,q2.y);
40
minp=LW_MIN(p1.y,p2.y);
41
maxp=LW_MAX(p1.y,p2.y);
43
if (minp>maxq || maxp<minq)
50
** @brief returns the kind of #CG_SEGMENT_INTERSECTION_TYPE behavior of lineseg 1 (constructed from p1 and p2) and lineseg 2 (constructed from q1 and q2)
51
** @param p1 start point of first straight linesegment
52
** @param p2 end point of first straight linesegment
53
** @param q1 start point of second line segment
54
** @param q2 end point of second line segment
55
** @return a #CG_SEGMENT_INTERSECTION_TYPE
58
** SEG_NO_INTERSECTION = 0,
60
** SEG_CROSS_LEFT = 2,
61
** SEG_CROSS_RIGHT = 3,
62
** SEG_TOUCH_LEFT = 4,
63
** SEG_TOUCH_RIGHT = 5
65
int lw_segment_intersects(POINT2D *p1, POINT2D *p2, POINT2D *q1, POINT2D *q2)
68
double pq1, pq2, qp1, qp2;
70
/* No envelope interaction => we are done. */
71
if (!lw_segment_envelope_intersects(*p1, *p2, *q1, *p2))
73
return SEG_NO_INTERSECTION;
76
/* Are the start and end points of q on the same side of p? */
77
pq1=lw_segment_side(p1,p2,q1);
78
pq2=lw_segment_side(p1,p2,q2);
79
if ((pq1>0 && pq2>0) || (pq1<0 && pq2<0))
81
return SEG_NO_INTERSECTION;
84
/* Are the start and end points of p on the same side of q? */
85
qp1=lw_segment_side(q1,q2,p1);
86
qp2=lw_segment_side(q1,q2,p2);
87
if ((qp1>0 && qp2>0) || (qp1<0 && qp2<0))
89
return SEG_NO_INTERSECTION;
92
/* Nobody is on one side or another? Must be colinear. */
93
if (pq1 == 0.0 && pq2 == 0.0 && qp1 == 0.0 && qp2 == 0.0)
99
** When one end-point touches, the sidedness is determined by the
100
** location of the other end-point.
105
return SEG_TOUCH_LEFT;
107
return SEG_TOUCH_RIGHT;
112
return SEG_TOUCH_LEFT;
114
return SEG_TOUCH_RIGHT;
117
/* The segments cross, what direction is the crossing? */
119
return SEG_CROSS_RIGHT;
121
return SEG_CROSS_LEFT;
123
/* This should never happen! */
129
** @brief lwline_crossing_direction: returns the kind of #CG_LINE_CROSS_TYPE behavior of 2 linestrings
130
** @param l1 first line string
131
** @param l2 second line string
132
** @return a #CG_LINE_CROSS_TYPE
134
** LINE_CROSS_LEFT = -1
135
** LINE_CROSS_RIGHT = 1
136
** LINE_MULTICROSS_END_LEFT = -2
137
** LINE_MULTICROSS_END_RIGHT = 2
138
** LINE_MULTICROSS_END_SAME_FIRST_LEFT = -3
139
** LINE_MULTICROSS_END_SAME_FIRST_RIGHT = 3
142
int lwline_crossing_direction(LWLINE *l1, LWLINE *l2)
145
int i = 0, j = 0, rv = 0;
157
int vertex_touch = -1;
158
int vertex_touch_type = 0;
160
pa1 = (POINTARRAY*)l1->points;
161
pa2 = (POINTARRAY*)l2->points;
163
p1 = lwalloc(sizeof(POINT2D));
164
p2 = lwalloc(sizeof(POINT2D));
165
q1 = lwalloc(sizeof(POINT2D));
166
q2 = lwalloc(sizeof(POINT2D));
168
/* One-point lines can't intersect (and shouldn't exist). */
169
if ( pa1->npoints < 2 || pa2->npoints < 2 )
170
return LINE_NO_CROSS;
172
LWDEBUGF(4, "lineCrossingDirection: l1 = %s", lwgeom_to_ewkt((LWGEOM*)l1,0));
173
LWDEBUGF(4, "lineCrossingDirection: l2 = %s", lwgeom_to_ewkt((LWGEOM*)l2,0));
175
for ( i = 1; i < pa2->npoints; i++ )
178
rv = getPoint2d_p(pa2, i-1, q1);
179
rv = getPoint2d_p(pa2, i, q2);
181
for ( j = 1; j < pa1->npoints; j++ )
184
rv = getPoint2d_p(pa1, j-1, p1);
185
rv = getPoint2d_p(pa1, j, p2);
187
LWDEBUGF(4, "lineCrossingDirection: i=%d, j=%d", i, j);
189
this_cross = lw_segment_intersects(p1, p2, q1, q2);
191
if ( ! first_cross && this_cross )
192
first_cross = this_cross;
194
final_cross = this_cross;
196
if ( this_cross == SEG_CROSS_LEFT )
202
if ( this_cross == SEG_CROSS_RIGHT )
209
** Crossing at a co-linearity can be turned into crossing at
210
** a vertex by pulling the original touch point forward along
213
if ( this_cross == SEG_COLINEAR && vertex_touch == (i-1) )
220
** Crossing-at-vertex will cause four segment touch interactions, two at
221
** j-1 and two at j. We avoid incrementing our cross count by ignoring the
224
if ( this_cross == SEG_TOUCH_LEFT )
226
if ( vertex_touch == (i-1) && vertex_touch_type == SEG_TOUCH_RIGHT )
230
vertex_touch_type = 0;
235
vertex_touch_type = this_cross;
239
if ( this_cross == SEG_TOUCH_RIGHT )
241
if ( vertex_touch == (i-1) && vertex_touch_type == SEG_TOUCH_LEFT )
245
vertex_touch_type = 0;
250
vertex_touch_type = this_cross;
256
** TODO Handle co-linear cases.
259
LWDEBUGF(4, "lineCrossingDirection: this_cross=%d, vertex_touch=%d, vertex_touch_type=%d", this_cross, vertex_touch, vertex_touch_type);
265
LWDEBUGF(4, "first_cross=%d, final_cross=%d, cross_left=%d, cross_right=%d", first_cross, final_cross, cross_left, cross_right);
272
if ( !cross_left && !cross_right )
273
return LINE_NO_CROSS;
275
if ( !cross_left && cross_right == 1 )
276
return LINE_CROSS_RIGHT;
278
if ( !cross_right && cross_left == 1 )
279
return LINE_CROSS_LEFT;
281
if ( cross_left - cross_right == 1 )
282
return LINE_MULTICROSS_END_LEFT;
284
if ( cross_left - cross_right == -1 )
285
return LINE_MULTICROSS_END_RIGHT;
287
if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_LEFT )
288
return LINE_MULTICROSS_END_SAME_FIRST_LEFT;
290
if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_RIGHT )
291
return LINE_MULTICROSS_END_SAME_FIRST_RIGHT;
293
if ( cross_left - cross_right == 0 && first_cross == SEG_TOUCH_LEFT )
294
return LINE_MULTICROSS_END_SAME_FIRST_RIGHT;
296
if ( cross_left - cross_right == 0 && first_cross == SEG_TOUCH_RIGHT )
297
return LINE_MULTICROSS_END_SAME_FIRST_LEFT;
299
return LINE_NO_CROSS;
304
** lwpoint_get_ordinate(point, ordinate) => double
306
double lwpoint_get_ordinate(const POINT4D *p, int ordinate)
310
lwerror("Null input geometry.");
314
if ( ordinate > 3 || ordinate < 0 )
316
lwerror("Cannot extract ordinate %d.", ordinate);
330
void lwpoint_set_ordinate(POINT4D *p, int ordinate, double value)
334
lwerror("Null input geometry.");
338
if ( ordinate > 3 || ordinate < 0 )
340
lwerror("Cannot extract ordinate %d.", ordinate);
344
LWDEBUGF(4, " setting ordinate %d to %g", ordinate, value);
364
int lwpoint_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int ndims, int ordinate, double interpolation_value)
366
double p1_value = lwpoint_get_ordinate(p1, ordinate);
367
double p2_value = lwpoint_get_ordinate(p2, ordinate);
371
if ( ordinate < 0 || ordinate >= ndims )
373
lwerror("Ordinate (%d) is not within ndims (%d).", ordinate, ndims);
377
if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
378
FP_MAX(p1_value, p2_value) < interpolation_value )
380
lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
384
proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
386
for ( i = 0; i < ndims; i++ )
388
double newordinate = 0.0;
389
p1_value = lwpoint_get_ordinate(p1, i);
390
p2_value = lwpoint_get_ordinate(p2, i);
391
newordinate = p1_value + proportion * (p2_value - p1_value);
392
lwpoint_set_ordinate(p, i, newordinate);
393
LWDEBUGF(4, " clip ordinate(%d) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", i, p1_value, p2_value, proportion, newordinate );
399
LWCOLLECTION *lwmline_clip_to_ordinate_range(LWMLINE *mline, int ordinate, double from, double to)
401
LWCOLLECTION *lwgeom_out = NULL;
405
lwerror("Null input geometry.");
409
if ( mline->ngeoms == 1)
411
lwgeom_out = lwline_clip_to_ordinate_range(mline->geoms[0], ordinate, from, to);
416
char hasz = TYPE_HASZ(mline->type);
417
char hasm = TYPE_HASM(mline->type);
418
char hassrid = TYPE_HASSRID(mline->type);
420
char homogeneous = 1;
421
size_t geoms_size = 0;
422
lwgeom_out = lwcollection_construct_empty(mline->SRID, hasz, hasm);
423
lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, MULTILINETYPE);
424
for ( i = 0; i < mline->ngeoms; i ++ )
426
col = lwline_clip_to_ordinate_range(mline->geoms[i], ordinate, from, to);
428
{ /* Something was left after the clip. */
429
if ( lwgeom_out->ngeoms + col->ngeoms > geoms_size )
432
if ( lwgeom_out->geoms )
434
lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, geoms_size * sizeof(LWGEOM*));
438
lwgeom_out->geoms = lwalloc(geoms_size * sizeof(LWGEOM*));
441
for ( j = 0; j < col->ngeoms; j++ )
443
lwgeom_out->geoms[lwgeom_out->ngeoms] = col->geoms[j];
444
lwgeom_out->ngeoms++;
446
if ( TYPE_GETTYPE(col->type) != TYPE_GETTYPE(mline->type) )
450
/* Shallow free the struct, leaving the geoms behind. */
451
if ( col->bbox ) lwfree(col->bbox);
455
lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
456
lwgeom_add_bbox((LWGEOM*)lwgeom_out);
459
lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, COLLECTIONTYPE);
463
if ( ! lwgeom_out || lwgeom_out->ngeoms == 0 ) /* Nothing left after clip. */
474
** lwline_clip_to_ordinate_range(line, ordinate, from, to) => lwmline
476
** Take in a LINESTRING and return a MULTILINESTRING of those portions of the
477
** LINESTRING between the from/to range for the specified ordinate (XYZM)
479
LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double from, double to)
482
POINTARRAY *pa_in = NULL;
483
LWCOLLECTION *lwgeom_out = NULL;
484
POINTARRAY *pa_out = NULL;
485
DYNPTARRAY *dp = NULL;
487
int added_last_point = 0;
488
POINT4D *p = NULL, *q = NULL, *r = NULL;
489
double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
490
char hasz = TYPE_HASZ(line->type);
491
char hasm = TYPE_HASM(line->type);
492
char dims = TYPE_NDIMS(line->type);
493
char hassrid = TYPE_HASSRID(line->type);
495
LWDEBUGF(4, "hassrid = %d", hassrid);
497
/* Null input, nothing we can do. */
500
lwerror("Null input geometry.");
504
/* Ensure 'from' is less than 'to'. */
512
LWDEBUGF(4, "from = %g, to = %g, ordinate = %d", from, to, ordinate);
513
LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line, PARSER_CHECK_NONE));
515
/* Asking for an ordinate we don't have. Error. */
516
if ( ordinate >= dims )
518
lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
522
/* Prepare our working point objects. */
523
p = lwalloc(sizeof(POINT4D));
524
q = lwalloc(sizeof(POINT4D));
525
r = lwalloc(sizeof(POINT4D));
527
/* Construct a collection to hold our outputs. */
528
lwgeom_out = lwalloc(sizeof(LWCOLLECTION));
529
lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, MULTILINETYPE);
531
lwgeom_out->SRID = line->SRID;
533
lwgeom_out->SRID = -1;
534
lwgeom_out->bbox = NULL;
535
lwgeom_out->ngeoms = 0;
536
lwgeom_out->geoms = NULL;
538
pa_in = (POINTARRAY*)line->points;
540
for ( i = 0; i < pa_in->npoints; i++ )
542
LWDEBUGF(4, "Point #%d", i);
549
ordinate_value_q = ordinate_value_p;
551
rv = getPoint4d_p(pa_in, i, p);
552
ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
553
LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
554
LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
556
/* Is this point inside the ordinate range? Yes. */
557
if ( ordinate_value_p >= from && ordinate_value_p <= to )
559
LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
561
if ( ! added_last_point )
563
LWDEBUG(4," new ptarray required");
564
/* We didn't add the previous point, so this is a new segment.
565
* Make a new point array. */
566
if ( dp ) lwfree(dp);
567
dp = dynptarray_create(64, line->type);
569
/* We're transiting into the range so add an interpolated
570
* point at the range boundary.
571
* If we're on a boundary and crossing from the far side,
572
* we also need an interpolated point. */
573
if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
574
( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
575
( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
576
( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
578
double interpolation_value;
579
(ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
580
rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
581
rv = dynptarray_addPoint4d(dp, r, 0);
582
LWDEBUGF(4, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
585
/* Add the current vertex to the point array. */
586
rv = dynptarray_addPoint4d(dp, p, 0);
587
if ( ordinate_value_p == from || ordinate_value_p == to )
589
added_last_point = 2; /* Added on boundary. */
593
added_last_point = 1; /* Added inside range. */
596
/* Is this point inside the ordinate range? No. */
599
LWDEBUGF(4, " added_last_point (%d)", added_last_point);
600
if ( added_last_point == 1 )
602
/* We're transiting out of the range, so add an interpolated point
603
* to the point array at the range boundary. */
604
double interpolation_value;
605
(ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
606
rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
607
rv = dynptarray_addPoint4d(dp, r, 0);
608
LWDEBUGF(4, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
610
else if ( added_last_point == 2 )
612
/* We're out and the last point was on the boundary.
613
* If the last point was the near boundary, nothing to do.
614
* If it was the far boundary, we need an interpolated point. */
616
(ordinate_value_q == from && ordinate_value_p > from) ||
617
(ordinate_value_q == to && ordinate_value_p < to) ) )
619
double interpolation_value;
620
(ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
621
rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
622
rv = dynptarray_addPoint4d(dp, r, 0);
623
LWDEBUGF(4, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
626
else if ( i && ordinate_value_q < from && ordinate_value_p > to )
628
/* We just hopped over the whole range, from bottom to top,
629
* so we need to add *two* interpolated points! */
630
pa_out = ptarray_construct(hasz, hasm, 2);
631
/* Interpolate lower point. */
632
rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
633
setPoint4d(pa_out, 0, r);
634
/* Interpolate upper point. */
635
rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
636
setPoint4d(pa_out, 1, r);
638
else if ( i && ordinate_value_q > to && ordinate_value_p < from )
640
/* We just hopped over the whole range, from top to bottom,
641
* so we need to add *two* interpolated points! */
642
pa_out = ptarray_construct(hasz, hasm, 2);
643
/* Interpolate upper point. */
644
rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
645
setPoint4d(pa_out, 0, r);
646
/* Interpolate lower point. */
647
rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
648
setPoint4d(pa_out, 1, r);
650
/* We have an extant point-array, save it out to a multi-line. */
654
LWDEBUG(4, "saving pointarray to multi-line (1)");
657
/* Only one point, so we have to make an lwpoint to hold this
658
* and set the overall output type to a generic collection. */
659
if ( dp->pa->npoints == 1 )
661
oline = (LWGEOM*)lwpoint_construct(line->SRID, NULL, dp->pa);
662
oline->type = lwgeom_makeType(hasz, hasm, hassrid, POINTTYPE);
663
lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, COLLECTIONTYPE);
667
oline = (LWGEOM*)lwline_construct(line->SRID, NULL, dp->pa);
668
oline->type = lwgeom_makeType(hasz, hasm, hassrid, LINETYPE);
673
oline = (LWGEOM*)lwline_construct(line->SRID, NULL, pa_out);
675
lwgeom_out->ngeoms++;
676
if ( lwgeom_out->geoms ) /* We can't just realloc, since repalloc chokes on a starting null ptr. */
678
lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, sizeof(LWGEOM*) * lwgeom_out->ngeoms);
682
lwgeom_out->geoms = lwalloc(sizeof(LWGEOM*) * lwgeom_out->ngeoms);
684
lwgeom_out->geoms[lwgeom_out->ngeoms - 1] = oline;
685
lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
686
lwgeom_add_bbox((LWGEOM*)lwgeom_out);
687
if ( dp ) lwfree(dp);
689
if ( pa_out ) pa_out = NULL;
691
added_last_point = 0;
696
/* Still some points left to be saved out. */
697
if ( dp && dp->pa->npoints > 0 )
700
LWDEBUG(4, "saving pointarray to multi-line (2)");
701
LWDEBUGF(4, "dp->pa->npoints == %d", dp->pa->npoints);
702
LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
704
if ( dp->pa->npoints == 1 )
706
oline = (LWGEOM*)lwpoint_construct(line->SRID, NULL, dp->pa);
707
oline->type = lwgeom_makeType(hasz, hasm, hassrid, POINTTYPE);
708
lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, COLLECTIONTYPE);
712
oline = (LWGEOM*)lwline_construct(line->SRID, NULL, dp->pa);
713
oline->type = lwgeom_makeType(hasz, hasm, hassrid, LINETYPE);
716
lwgeom_out->ngeoms++;
717
if ( lwgeom_out->geoms ) /* We can't just realloc, since repalloc chokes on a starting null ptr. */
719
lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, sizeof(LWGEOM*) * lwgeom_out->ngeoms);
723
lwgeom_out->geoms = lwalloc(sizeof(LWGEOM*) * lwgeom_out->ngeoms);
725
lwgeom_out->geoms[lwgeom_out->ngeoms - 1] = oline;
726
lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
727
lwgeom_add_bbox((LWGEOM*)lwgeom_out);
728
if ( dp ) lwfree(dp);
736
if ( lwgeom_out->ngeoms > 0 )
743
static char *base32 = "0123456789bcdefghjkmnpqrstuvwxyz";
746
** Calculate the geohash, iterating downwards and gaining precision.
747
** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>
748
** Released under the MIT License.
750
char *geohash_point(double longitude, double latitude, int precision)
753
double lat[2], lon[2], mid;
754
char bits[] = {16,8,4,2,1};
756
char *geohash = NULL;
758
geohash = lwalloc(precision + 1);
765
while (i < precision)
769
mid = (lon[0] + lon[1]) / 2;
782
mid = (lat[0] + lat[1]) / 2;
801
geohash[i++] = base32[ch];
810
int lwgeom_geohash_precision(BOX3D bbox, BOX3D *bounds)
812
double minx, miny, maxx, maxy;
813
double latmax, latmin, lonmax, lonmin;
814
double lonwidth, latwidth;
815
double latmaxadjust, lonmaxadjust, latminadjust, lonminadjust;
818
/* Get the bounding box, return error if things don't work out. */
824
if ( minx == maxx && miny == maxy )
826
/* It's a point. Doubles have 51 bits of precision.
827
** 2 * 51 / 5 == 20 */
836
/* Shrink a world bounding box until one of the edges interferes with the
837
** bounds of our rectangle. */
840
lonwidth = lonmax - lonmin;
841
latwidth = latmax - latmin;
842
latmaxadjust = lonmaxadjust = latminadjust = lonminadjust = 0.0;
844
if ( minx > lonmin + lonwidth / 2.0 )
846
lonminadjust = lonwidth / 2.0;
848
else if ( maxx < lonmax - lonwidth / 2.0 )
850
lonmaxadjust = -1 * lonwidth / 2.0;
852
if ( miny > latmin + latwidth / 2.0 )
854
latminadjust = latwidth / 2.0;
856
else if (maxy < latmax - latwidth / 2.0 )
858
latmaxadjust = -1 * latwidth / 2.0;
860
/* Only adjust if adjustments are legal (we haven't crossed any edges). */
861
if ( (lonminadjust || lonmaxadjust) && (latminadjust || latmaxadjust ) )
863
latmin += latminadjust;
864
lonmin += lonminadjust;
865
latmax += latmaxadjust;
866
lonmax += lonmaxadjust;
867
/* Each adjustment cycle corresponds to 2 bits of storage in the
877
/* Save the edges of our bounds, in case someone cares later. */
878
bounds->xmin = lonmin;
879
bounds->xmax = lonmax;
880
bounds->ymin = latmin;
881
bounds->ymax = latmax;
883
/* Each geohash character (base32) can contain 5 bits of information.
884
** We are returning the precision in characters, so here we divide. */
885
return precision / 5;
890
** Return a geohash string for the geometry. <http://geohash.org>
891
** Where the precision is non-positive, calculate a precision based on the
892
** bounds of the feature. Big features have loose precision.
893
** Small features have tight precision.
895
char *lwgeom_geohash(const LWGEOM *lwgeom, int precision)
898
BOX3D precision_bounds;
901
bbox = lwgeom_compute_box3d(lwgeom);
902
if ( ! bbox ) return NULL;
904
/* Return error if we are being fed something outside our working bounds */
905
if ( bbox->xmin < -180 || bbox->ymin < -90 || bbox->xmax > 180 || bbox->ymax > 90 )
907
lwerror("Geohash requires inputs in decimal degrees.");
912
/* What is the center of our geometry bounds? We'll use that to
913
** approximate location. */
914
lon = bbox->xmin + (bbox->xmax - bbox->xmin) / 2;
915
lat = bbox->ymin + (bbox->ymax - bbox->ymin) / 2;
917
if ( precision <= 0 )
919
precision = lwgeom_geohash_precision(*bbox, &precision_bounds);
925
** Return the geohash of the center, with a precision determined by the
926
** extent of the bounds.
927
** Possible change: return the point at the center of the precision bounds?
929
return geohash_point(lon, lat, precision);
added
removed
liblwgeom/lwalgorithm.c
