[postgis-commits] svn - r3685 - in trunk/liblwgeom: . cunit
postgis-commits at postgis.refractions.net
postgis-commits at postgis.refractions.net
Tue Feb 10 18:11:25 PST 2009
Author: pramsey
Date: 2009-02-10 18:11:24 -0800 (Tue, 10 Feb 2009)
New Revision: 3685
Modified:
trunk/liblwgeom/cunit/cu_algorithm.c
trunk/liblwgeom/cunit/cu_algorithm.h
trunk/liblwgeom/liblwgeom.h
trunk/liblwgeom/lwalgorithm.c
trunk/liblwgeom/lwalgorithm.h
trunk/liblwgeom/lwcollection.c
trunk/liblwgeom/lwgeom.c
Log:
GeoHash implementation first cut.
Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c 2009-02-10 21:40:33 UTC (rev 3684)
+++ trunk/liblwgeom/cunit/cu_algorithm.c 2009-02-11 02:11:24 UTC (rev 3685)
@@ -35,7 +35,9 @@
(NULL == CU_add_test(pSuite, "test_lwpoint_interpolate()", test_lwpoint_interpolate)) ||
(NULL == CU_add_test(pSuite, "test_lwline_clip()", test_lwline_clip)) ||
(NULL == CU_add_test(pSuite, "test_lwline_clip_big()", test_lwline_clip_big)) ||
- (NULL == CU_add_test(pSuite, "test_lwmline_clip()", test_lwmline_clip))
+ (NULL == CU_add_test(pSuite, "test_lwmline_clip()", test_lwmline_clip)) ||
+ (NULL == CU_add_test(pSuite, "test_geohash_precision()", test_geohash_precision)) ||
+ (NULL == CU_add_test(pSuite, "test_geohash()", test_geohash))
)
{
CU_cleanup_registry();
@@ -701,3 +703,71 @@
lwline_free(line);
}
+void test_geohash_precision(void)
+{
+ BOX3D bbox;
+ int precision = 0;
+
+ bbox.xmin = 23.0;
+ bbox.xmax = 23.0;
+ bbox.ymin = 25.2;
+ bbox.ymax = 25.2;
+ precision = lwgeom_geohash_precision(bbox);
+ //printf("precision %d\n",precision);
+ CU_ASSERT_EQUAL(precision, 20);
+
+ bbox.xmin = 23.0;
+ bbox.ymin = 23.0;
+ bbox.xmax = 23.1;
+ bbox.ymax = 23.1;
+ precision = lwgeom_geohash_precision(bbox);
+ //printf("precision %d\n",precision);
+ CU_ASSERT_EQUAL(precision, 2);
+
+ bbox.xmin = 23.0;
+ bbox.ymin = 23.0;
+ bbox.xmax = 23.0001;
+ bbox.ymax = 23.0001;
+ precision = lwgeom_geohash_precision(bbox);
+ //printf("precision %d\n",precision);
+ CU_ASSERT_EQUAL(precision, 6);
+
+}
+
+void test_geohash(void)
+{
+ LWPOINT *lwpoint = NULL;
+ LWLINE *lwline = NULL;
+ LWMLINE *lwmline = NULL;
+ char *geohash = NULL;
+
+ lwpoint = (LWPOINT*)lwgeom_from_ewkt("POINT(23.0 25.2)", PARSER_CHECK_NONE);
+ geohash = lwgeom_geohash((LWGEOM*)lwpoint);
+ printf("geohash %s\n",geohash);
+ CU_ASSERT_STRING_EQUAL(geohash, "20");
+ lwfree(lwpoint);
+ lwfree(geohash);
+
+ lwline = (LWLINE*)lwgeom_from_ewkt("LINESTRING(23.0 23.0,23.1 23.1)", PARSER_CHECK_NONE);
+ geohash = lwgeom_geohash((LWGEOM*)lwline);
+ printf("geohash %s\n",geohash);
+ CU_ASSERT_STRING_EQUAL(geohash, "20");
+ lwfree(lwline);
+ lwfree(geohash);
+
+ lwline = (LWLINE*)lwgeom_from_ewkt("LINESTRING(23.0 23.0,23.001 23.001)", PARSER_CHECK_NONE);
+ geohash = lwgeom_geohash((LWGEOM*)lwline);
+ printf("geohash %s\n",geohash);
+ CU_ASSERT_STRING_EQUAL(geohash, "20");
+ lwfree(lwline);
+ lwfree(geohash);
+
+ lwmline = (LWMLINE*)lwgeom_from_ewkt("MULTILINESTRING((23.0 23.0,23.1 23.1),(23.0 23.0,23.1 23.1))", PARSER_CHECK_NONE);
+ geohash = lwgeom_geohash((LWGEOM*)lwmline);
+ printf("geohash %s\n",geohash);
+ CU_ASSERT_STRING_EQUAL(geohash, "20");
+ lwfree(lwmline);
+ lwfree(geohash);
+}
+
+
Modified: trunk/liblwgeom/cunit/cu_algorithm.h
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.h 2009-02-10 21:40:33 UTC (rev 3684)
+++ trunk/liblwgeom/cunit/cu_algorithm.h 2009-02-11 02:11:24 UTC (rev 3685)
@@ -37,4 +37,5 @@
void test_lwline_clip(void);
void test_lwline_clip_big(void);
void test_lwmline_clip(void);
-
+void test_geohash_precision(void);
+void test_geohash(void);
Modified: trunk/liblwgeom/liblwgeom.h
===================================================================
--- trunk/liblwgeom/liblwgeom.h 2009-02-10 21:40:33 UTC (rev 3684)
+++ trunk/liblwgeom/liblwgeom.h 2009-02-11 02:11:24 UTC (rev 3685)
@@ -843,9 +843,9 @@
* pointers to the serialized form (POINTARRAYs).
*/
LWGEOM *lwgeom_deserialize(uchar *serializedform);
+BOX3D *lwgeom_compute_box3d(const LWGEOM *geom);
-
/******************************************************************
* LWMULTIx and LWCOLLECTION functions
******************************************************************/
@@ -859,10 +859,10 @@
LWMCURVE *lwmcurve_deserialize(uchar *serialized_form);
LWMSURFACE *lwmsurface_deserialize(uchar *serialized_form);
-LWGEOM *lwcollection_getsubgeom(LWCOLLECTION *, int);
+LWGEOM *lwcollection_getsubgeom(LWCOLLECTION *col, int gnum);
+BOX3D *lwcollection_compute_box3d(LWCOLLECTION *col);
-
/******************************************************************
* SERIALIZED FORM functions
******************************************************************/
Modified: trunk/liblwgeom/lwalgorithm.c
===================================================================
--- trunk/liblwgeom/lwalgorithm.c 2009-02-10 21:40:33 UTC (rev 3684)
+++ trunk/liblwgeom/lwalgorithm.c 2009-02-11 02:11:24 UTC (rev 3685)
@@ -734,3 +734,185 @@
return NULL;
}
+
+int lwgeom_geohash_precision(BOX3D bbox)
+{
+ double minx, miny, maxx, maxy;
+ double latmax, latmin, lonmax, lonmin;
+ double lonwidth, latwidth;
+ double latmaxadjust, lonmaxadjust, latminadjust, lonminadjust;
+ int precision = 0;
+
+ /* Get the bounding box, return error if things don't work out. */
+ minx = bbox.xmin;
+ miny = bbox.ymin;
+ maxx = bbox.xmax;
+ maxy = bbox.ymax;
+
+ if( minx == maxx && miny == maxy )
+ {
+ /* It's a point. Doubles have 51 bits of precision.
+ ** 2 * 51 / 5 == 20 */
+ return 20;
+ }
+
+ lonmin = -180.0;
+ latmin = -90.0;
+ lonmax = 180.0;
+ latmax = 90.0;
+
+ /* Shrink a world bounding box until one of the edges interferes with the
+ ** bounds of our rectangle. */
+ while( 1 )
+ {
+ lonwidth = lonmax - lonmin;
+ latwidth = latmax - latmin;
+ latmaxadjust = lonmaxadjust = latminadjust = lonminadjust = 0.0;
+
+ if( minx > lonmin + lonwidth / 2.0 )
+ {
+ lonminadjust = lonwidth / 2.0;
+ }
+ else if ( maxx < lonmax - lonwidth / 2.0 )
+ {
+ lonmaxadjust = -1 * lonwidth / 2.0;
+ }
+ if( miny > latmin + latwidth / 2.0 )
+ {
+ latminadjust = latwidth / 2.0;
+ }
+ else if (maxy < latmax - latwidth / 2.0 )
+ {
+ latmaxadjust = -1 * latwidth / 2.0;
+ }
+ if ( (lonminadjust || lonmaxadjust) && (latminadjust || latmaxadjust ) )
+ {
+ latmin += latminadjust;
+ lonmin += lonminadjust;
+ latmax += latmaxadjust;
+ lonmax += lonmaxadjust;
+ precision++;
+ }
+ else
+ {
+ break;
+ }
+ }
+
+ /* Each adjustment above corresponds to 2 bits of storage in the
+ ** geohash. Each geohash character can contain 5 bits of information.
+ ** So geohashes have even lengths. 2 characters corresponds to 10 bits
+ ** of storage, which is 5 adjustments. */
+ return 2 * ( precision / 5 );
+}
+
+char *lwgeom_geohash(const LWGEOM *lwgeom)
+{
+ BOX3D *bbox = NULL;
+ int precision = 0;
+ double lat, lon;
+
+ bbox = lwgeom_compute_box3d(lwgeom);
+ if( ! bbox ) return NULL;
+
+ if ( bbox->xmin < -180 || bbox->ymin < -90 || bbox->xmax > 180 || bbox->ymax > 90 )
+ {
+ lwerror("Geohash requires inputs in decimal degrees.");
+ lwfree(bbox);
+ return NULL;
+ }
+
+ /* What is the center of our geometry? */
+ lon = bbox->xmin + (bbox->xmax - bbox->xmin) / 2;
+ lat = bbox->ymin + (bbox->ymax - bbox->ymin) / 2;
+
+ precision = lwgeom_geohash_precision(*bbox);
+
+ lwfree(bbox);
+
+ /* Return the geohash of the center, with a precision determined by the
+ ** extent of the bounds. */
+ return geohash_point(lat, lon, precision);
+}
+
+static char *base32 = "0123456789bcdefghjkmnpqrstuvwxyz";
+
+char *geohash_point(double latitude, double longitude, int precision)
+{
+ int is_even=1, i=0;
+ double lat[2], lon[2], mid;
+ char bits[] = {16,8,4,2,1};
+ int bit=0, ch=0;
+ char *geohash = NULL;
+
+ geohash = lwalloc(precision + 1);
+
+ lat[0] = -90.0; lat[1] = 90.0;
+ lon[0] = -180.0; lon[1] = 180.0;
+
+ while (i < precision)
+ {
+ if (is_even)
+ {
+ mid = (lon[0] + lon[1]) / 2;
+ if (longitude > mid)
+ {
+ ch |= bits[bit];
+ lon[0] = mid;
+ }
+ else
+ {
+ lon[1] = mid;
+ }
+ }
+ else
+ {
+ mid = (lat[0] + lat[1]) / 2;
+ if (latitude > mid)
+ {
+ ch |= bits[bit];
+ lat[0] = mid;
+ }
+ else
+ {
+ lat[1] = mid;
+ }
+ }
+
+ is_even = !is_even;
+ if (bit < 4)
+ {
+ bit++;
+ }
+ else
+ {
+ geohash[i++] = base32[ch];
+ bit = 0;
+ ch = 0;
+ }
+ }
+ geohash[i] = 0;
+ return geohash;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Modified: trunk/liblwgeom/lwalgorithm.h
===================================================================
--- trunk/liblwgeom/lwalgorithm.h 2009-02-10 21:40:33 UTC (rev 3684)
+++ trunk/liblwgeom/lwalgorithm.h 2009-02-11 02:11:24 UTC (rev 3685)
@@ -45,3 +45,8 @@
int lwpoint_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int ndims, int ordinate, double interpolation_value);
LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double from, double to);
LWCOLLECTION *lwmline_clip_to_ordinate_range(LWMLINE *mline, int ordinate, double from, double to);
+
+int lwgeom_geohash_precision(BOX3D bbox);
+char *lwgeom_geohash(const LWGEOM *lwgeom);
+char *geohash_point(double latitude, double longitude, int precision);
+
Modified: trunk/liblwgeom/lwcollection.c
===================================================================
--- trunk/liblwgeom/lwcollection.c 2009-02-10 21:40:33 UTC (rev 3684)
+++ trunk/liblwgeom/lwcollection.c 2009-02-11 02:11:24 UTC (rev 3685)
@@ -487,3 +487,54 @@
};
+BOX3D *lwcollection_compute_box3d(LWCOLLECTION *col)
+{
+ int i;
+ BOX3D *boxfinal = NULL;
+ BOX3D *boxtmp1 = NULL;
+ BOX3D *boxtmp2 = NULL;
+ for ( i = 0; i < col->ngeoms; i++ )
+ {
+ if( col->geoms[i] ) {
+ switch( TYPE_GETTYPE(col->geoms[i]->type) )
+ {
+ case POINTTYPE:
+ boxtmp1 = lwpoint_compute_box3d((LWPOINT*)(col->geoms[i]));
+ break;
+ case LINETYPE:
+ boxtmp1 = lwline_compute_box3d((LWLINE*)(col->geoms[i]));
+ break;
+ case POLYGONTYPE:
+ boxtmp1 = lwpoly_compute_box3d((LWPOLY*)(col->geoms[i]));
+ break;
+ case CIRCSTRINGTYPE:
+ boxtmp1 = lwcircstring_compute_box3d((LWCIRCSTRING *)(col->geoms[i]));
+ break;
+ case COMPOUNDTYPE:
+ case CURVEPOLYTYPE:
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case MULTICURVETYPE:
+ case MULTISURFACETYPE:
+ case COLLECTIONTYPE:
+ boxtmp1 = lwcollection_compute_box3d((LWCOLLECTION*)(col->geoms[i]));
+ boxfinal = box3d_union(boxtmp1, boxtmp2);
+ break;
+ }
+ boxtmp2 = boxfinal;
+ boxfinal = box3d_union(boxtmp1, boxtmp2);
+ if( boxtmp1 && boxtmp1 != boxfinal )
+ {
+ lwfree(boxtmp1);
+ boxtmp1 = NULL;
+ }
+ if( boxtmp2 && boxtmp2 != boxfinal )
+ {
+ lwfree(boxtmp2);
+ boxtmp2 = NULL;
+ }
+ }
+ }
+ return boxfinal;
+}
Modified: trunk/liblwgeom/lwgeom.c
===================================================================
--- trunk/liblwgeom/lwgeom.c 2009-02-10 21:40:33 UTC (rev 3684)
+++ trunk/liblwgeom/lwgeom.c 2009-02-11 02:11:24 UTC (rev 3685)
@@ -199,6 +199,34 @@
}
}
+BOX3D *lwgeom_compute_box3d(const LWGEOM *lwgeom)
+{
+ if( ! lwgeom ) return NULL;
+
+ switch(TYPE_GETTYPE(lwgeom->type))
+ {
+ case POINTTYPE:
+ return lwpoint_compute_box3d((LWPOINT *)lwgeom);
+ case LINETYPE:
+ return lwline_compute_box3d((LWLINE *)lwgeom);
+ case CIRCSTRINGTYPE:
+ return lwcircstring_compute_box3d((LWCIRCSTRING *)lwgeom);
+ case POLYGONTYPE:
+ return lwpoly_compute_box3d((LWPOLY *)lwgeom);
+ case COMPOUNDTYPE:
+ case CURVEPOLYTYPE:
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTICURVETYPE:
+ case MULTIPOLYGONTYPE:
+ case MULTISURFACETYPE:
+ case COLLECTIONTYPE:
+ return lwcollection_compute_box3d((LWCOLLECTION *)lwgeom);
+ }
+ /* Never get here, please. */
+ return NULL;
+}
+
int
lwgeom_compute_box2d_p(LWGEOM *lwgeom, BOX2DFLOAT4 *buf)
{
More information about the postgis-commits
mailing list