Changeset: f8a1ea4c5709 for MonetDB
URL: https://dev.monetdb.org/hg/MonetDB/rev/f8a1ea4c5709
Modified Files:
geom/monetdb5/geod.c
geom/monetdb5/geom.c
geom/monetdb5/geom.h
geom/monetdb5/geomBulk.c
geom/monetdb5/geom_atoms.c
geom/sql/40_geom.sql
Branch: geo-update-dev
Log Message:
Implemented the rtree build in the bulk mbr function. Implemented wkbIntersects
filter function using the rtree index to filter out candidates before doing the
actual intersects calculation. Implemented the bulk verison of st_transform,
with substancial improvements. Implemented mbr intersects.
diffs (truncated from 726 to 300 lines):
diff --git a/geom/monetdb5/geod.c b/geom/monetdb5/geod.c
--- a/geom/monetdb5/geod.c
+++ b/geom/monetdb5/geod.c
@@ -1166,7 +1166,6 @@ free:
return msg;
}
-//select st_distancegeographic(q1.geom,q2.geom) from brittany_ports as q1 join
wpi_ports as q2 on [q1.geom] st_dwithingeographic [q2.geom,5000];
str
wkbDWithinGeographicJoin(bat *lres_id, bat *rres_id, const bat *l_id, const
bat *r_id, const bat *d_id, const bat *ls_id, const bat *rs_id, bit
*nil_matches, lng *estimate, bit *anti) {
double distance_within = 0;
@@ -1184,19 +1183,16 @@ wkbDWithinGeographicJoin(bat *lres_id, b
return
filterJoinGeomGeomDoubleToBit(lres_id,rres_id,l_id,r_id,distance_within,ls_id,rs_id,*nil_matches,estimate,*anti,geosDistanceWithin,"geom.wkbDWithinGeographicJoin");
}
-//select
gml_id,st_distancegeographic(st_setsrid(st_makepoint(-4.50,48.32),4326),q1.geom)
from brittany_ports as q1 where [q1.geom] st_dwithingeographic
[st_setsrid(st_makepoint(-4.50,48.32),4326),5000];
str
wkbDWithinGeographicSelect(bat* outid, const bat *bid , const bat *sid, wkb
**wkb_const, dbl *distance_within, bit *anti) {
return
filterSelectGeomGeomDoubleToBit(outid,bid,sid,*wkb_const,*distance_within,*anti,geosDistanceWithin,"geom.wkbDWithinGeographicSelect");
}
-//select q1.gml_id,q2.gml_id,st_distancegeographic(q1.geom,q2.geom) from
brittany_ports as q1 join brittany_ports as q2 on [q1.geom]
ST_IntersectsGeographic [q2.geom];
str
wkbIntersectsGeographicJoin(bat *lres_id, bat *rres_id, const bat *l_id, const
bat *r_id, const bat *ls_id, const bat *rs_id, bit *nil_matches, lng *estimate,
bit *anti) {
return
filterJoinGeomGeomDoubleToBit(lres_id,rres_id,l_id,r_id,0,ls_id,rs_id,*nil_matches,estimate,*anti,geosDistanceWithin,"geom.wkbIntersectsGeographicJoin");
}
-//select
gml_id,st_distancegeographic(st_setsrid(st_makepoint(-1.7717988686592332,48.602742696725414),4326),q1.geom)
from brittany_ports as q1 where [q1.geom] ST_IntersectsGeographic
[st_setsrid(st_makepoint(-1.7717988686592332,48.602742696725414),4326)];
str
wkbIntersectsGeographicSelect(bat* outid, const bat *bid , const bat *sid, wkb
**wkb_const, bit *anti) {
return
filterSelectGeomGeomDoubleToBit(outid,bid,sid,*wkb_const,0,*anti,geosDistanceWithin,"geom.wkbIntersectsGeographicSelect");
diff --git a/geom/monetdb5/geom.c b/geom/monetdb5/geom.c
--- a/geom/monetdb5/geom.c
+++ b/geom/monetdb5/geom.c
@@ -16,6 +16,7 @@
#include "geom_atoms.h"
#include "gdk_logger.h"
#include "mal_exception.h"
+#include "gdk_rtree.h"
mbr mbrNIL = {0}; // will be initialized properly by geom prelude
@@ -369,7 +370,7 @@ transformCoordSeq(int idx, int coordinat
return MAL_SUCCEED;
}
-static str
+str
transformPoint(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P)
{
int coordinatesNum = 0;
@@ -409,7 +410,7 @@ transformPoint(GEOSGeometry **transforme
return MAL_SUCCEED;
}
-static str
+str
transformLine(GEOSCoordSeq *gcs_new, const GEOSGeometry *geosGeometry, PJ *P)
{
int coordinatesNum = 0;
@@ -447,7 +448,7 @@ transformLine(GEOSCoordSeq *gcs_new, con
return MAL_SUCCEED;
}
-static str
+str
transformLineString(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P)
{
GEOSCoordSeq coordSeq;
@@ -469,7 +470,7 @@ transformLineString(GEOSGeometry **trans
return ret;
}
-static str
+str
transformLinearRing(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P)
{
GEOSCoordSeq coordSeq = NULL;
@@ -491,7 +492,7 @@ transformLinearRing(GEOSGeometry **trans
return ret;
}
-static str
+str
transformPolygon(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P, int srid)
{
const GEOSGeometry *exteriorRingGeometry;
@@ -556,7 +557,7 @@ transformPolygon(GEOSGeometry **transfor
return ret;
}
-static str
+str
transformMultiGeometry(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P, int srid, int geometryType)
{
int geometriesNum, subGeometryType, i;
@@ -5513,11 +5514,16 @@ static mel_func geom_init_funcs[] = {
command("geom", "IntersectsGeographicselect", wkbIntersectsGeographicSelect,
false, "TODO", args(1, 5, batarg("", oid), batarg("b", wkb), batarg("s", oid),
arg("c", wkb), arg("anti",bit))),
command("geom", "IntersectsGeographicjoin", wkbIntersectsGeographicJoin,
false, "TODO", args(2, 9, batarg("lr",oid),batarg("rr",oid), batarg("a", wkb),
batarg("b", wkb),
batarg("sl",oid),batarg("sr",oid),arg("nil_matches",bit),arg("estimate",lng),arg("anti",bit))),
+ command("geom", "Intersects", wkbIntersects, false, "Returns true if these
Geometries 'spatially intersect in 2D'", args(1,3,
arg("",bit),arg("a",wkb),arg("b",wkb))),
+ command("geom", "Intersectsselect", wkbIntersectsSelect, false, "TODO",
args(1, 5, batarg("", oid), batarg("b", wkb), batarg("s", oid), arg("c", wkb),
arg("anti",bit))),
+ command("geom", "Intersectsjoin", wkbIntersectsJoin, false, "TODO", args(2,
8, batarg("lr",oid),batarg("rr",oid), batarg("a", wkb), batarg("b", wkb),
batarg("sl",oid),batarg("sr",oid),arg("nil_matches",bit),arg("estimate",lng))),
+
+ command("geom", "IntersectsMBR", mbrIntersects, false, "TODO", args(1,3,
arg("",bit),arg("a",mbr),arg("b",mbr))),
+
command("aggr", "Collect", wkbCollectAggr, false, "TODO", args(1, 2, arg("",
wkb), batarg("val", wkb))),
command("aggr", "subCollect", wkbCollectAggrSubGrouped, false, "TODO",
args(1, 5, batarg("", wkb), batarg("val", wkb), batarg("g", oid), batarg("e",
oid), arg("skip_nils", bit))),
command("aggr", "subCollect", wkbCollectAggrSubGroupedCand, false, "TODO",
args(1, 6, batarg("", wkb), batarg("val", wkb), batarg("g", oid),
batargany("e", 1), batarg("g", oid), arg("skip_nils", bit))),
- //TODO: See if we can remove (used in SQL level on 39_spatial_ref_sys.sql)
command("geom", "hasZ", geoHasZ, false, "returns 1 if the geometry has z
coordinate", args(1,2, arg("",int),arg("flags",int))),
command("geom", "hasM", geoHasM, false, "returns 1 if the geometry has m
coordinate", args(1,2, arg("",int),arg("flags",int))),
command("geom", "getType", geoGetType, false, "returns the str representation
of the geometry type", args(1,3,
arg("",str),arg("flags",int),arg("format",int))),
@@ -5567,7 +5573,6 @@ static mel_func geom_init_funcs[] = {
command("geom", "Crosses", wkbCrosses, false, "Returns TRUE if the supplied
geometries have some, but not all, interior points in common.", args(1,3,
arg("",bit),arg("a",wkb),arg("b",wkb))),
command("geom", "Disjoint", wkbDisjoint, false, "Returns true if these
Geometries are 'spatially disjoint'", args(1,3,
arg("",bit),arg("a",wkb),arg("b",wkb))),
command("geom", "Equals", wkbEquals, false, "Returns true if the given
geometries represent the same geometry. Directionality is ignored.", args(1,3,
arg("",bit),arg("a",wkb),arg("b",wkb))),
- command("geom", "Intersects", wkbIntersects, false, "Returns true if these
Geometries 'spatially intersect in 2D'", args(1,3,
arg("",bit),arg("a",wkb),arg("b",wkb))),
command("geom", "Overlaps", wkbOverlaps, false, "Returns TRUE if the
Geometries intersect but are not completely contained by each other.",
args(1,3, arg("",bit),arg("a",wkb),arg("b",wkb))),
command("geom", "Relate", wkbRelate, false, "Returns true if the Geometry a
'spatially related' to Geometry b, by testing for intersection between the
Interior, Boundary and Exterior of the two geometries as specified by the
values in the intersectionPatternMatrix.", args(1,4,
arg("",bit),arg("a",wkb),arg("b",wkb),arg("intersection_matrix_pattern",str))),
command("geom", "Touches", wkbTouches, false, "Returns TRUE if the geometries
have at least one point in common, but their interiors do not intersect.",
args(1,3, arg("",bit),arg("a",wkb),arg("b",wkb))),
@@ -5655,6 +5660,7 @@ static mel_func geom_init_funcs[] = {
command("batgeom", "mbr", wkbMBR_bat, false, "Creates the mbr for the given
wkb.", args(1,2, batarg("",mbr),batarg("",wkb))),
command("batgeom", "coordinateFromWKB", wkbCoordinateFromWKB_bat, false,
"returns xmin (=1), ymin (=2), xmax (=3) or ymax(=4) of the provided geometry",
args(1,3, batarg("",dbl),batarg("",wkb),arg("",int))),
command("batgeom", "coordinateFromMBR", wkbCoordinateFromMBR_bat, false,
"returns xmin (=1), ymin (=2), xmax (=3) or ymax(=4) of the provided mbr",
args(1,3, batarg("",dbl),batarg("",mbr),arg("",int))),
+ command("batgeom", "Transform", wkbTransform_bat, false, "Transforms a bat of
geometries from one srid to another", args(1,6,
batarg("",wkb),batarg("g",wkb),arg("srid_src",int),arg("srid_dst",int),arg("proj_src",str),arg("proj_dest",str))),
command("calc", "mbr", mbrFromString, false, "", args(1,2,
arg("",mbr),arg("v",str))),
command("calc", "mbr", mbrFromMBR, false, "", args(1,2,
arg("",mbr),arg("v",mbr))),
command("calc", "wkb", wkbFromWKB, false, "It is called when adding a new
geometry column to an existing table", args(1,2, arg("",wkb),arg("v",wkb))),
diff --git a/geom/monetdb5/geom.h b/geom/monetdb5/geom.h
--- a/geom/monetdb5/geom.h
+++ b/geom/monetdb5/geom.h
@@ -169,7 +169,15 @@ geom_export str wkbGeometryN_bat(bat *ou
geom_export str wkbNumGeometries(int* out, wkb** geom);
geom_export str wkbNumGeometries_bat(bat *outBAT_id, bat *inBAT_id);
+str transformPoint(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P);
+str transformLine(GEOSCoordSeq *gcs_new, const GEOSGeometry *geosGeometry, PJ
*P);
+str transformLineString(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P);
+str transformLinearRing(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P);
+str transformPolygon(GEOSGeometry **transformedGeometry, const GEOSGeometry
*geosGeometry, PJ *P, int srid);
+str transformMultiGeometry(GEOSGeometry **transformedGeometry, const
GEOSGeometry *geosGeometry, PJ *P, int srid, int geometryType);
geom_export str wkbTransform(wkb**, wkb**, int*, int*, char**, char**);
+geom_export str wkbTransform_bat(bat *outBAT_id, bat *inBAT_id, int *srid_src,
int *srid_dst, char **proj4_src_str, char **proj4_dst_str);
+geom_export str wkbTransform_bat_cand(bat *outBAT_id, bat *inBAT_id, bat
*s_id, int *srid_src, int *srid_dst, char **proj4_src_str, char
**proj4_dst_str);
geom_export str wkbTranslate(wkb**, wkb**, dbl*, dbl*, dbl*);
geom_export str wkbDelaunayTriangles(wkb**, wkb**, dbl*, int*);
geom_export str wkbPointOnSurface(wkb**, wkb**);
@@ -208,3 +216,7 @@ geom_export str wkbCoordinateFromWKB_bat
geom_export str wkbCoordinateFromMBR_bat(bat *outBAT_id, bat *inBAT_id, int*
coordinateIdx);
geom_export str geom_sql_upgrade(int);
+
+geom_export str wkbIntersectsJoin(bat *lres_id, bat *rres_id, const bat *l_id,
const bat *r_id, const bat *ls_id, const bat *rs_id, bit *nil_matches, lng
*estimate);
+geom_export str wkbIntersectsSelect(bat* outid, const bat *bid , const bat
*sid, wkb **wkb_const, bit *anti);
+geom_export str mbrIntersects(bit* out, mbr** mbr1, mbr** mbr2);
diff --git a/geom/monetdb5/geomBulk.c b/geom/monetdb5/geomBulk.c
--- a/geom/monetdb5/geomBulk.c
+++ b/geom/monetdb5/geomBulk.c
@@ -13,10 +13,460 @@
#include "geom.h"
#include "geod.h"
#include "geom_atoms.h"
+#include "gdk_rtree.h"
/********** Geo Update **********/
+static str
+filterSelectGeomGeomToBitIndex(bat* outid, const bat *bid , const bat *sid,
wkb *wkb_const, bit anti, char (*func) (const GEOSGeometry *, const
GEOSGeometry *), const char *name)
+{
+ BAT *out = NULL, *b = NULL, *s = NULL;
+ BATiter b_iter;
+ struct canditer ci;
+ GEOSGeom col_geom, const_geom;
-//TODO: Rename these functions with Stefanos
+ //Check if the geometry is null and convert to GEOS
+ if ((const_geom = wkb2geos(wkb_const)) == NULL) {
+ if ((out = BATdense(0, 0, 0)) == NULL)
+ throw(MAL, name, GDK_EXCEPTION);
+ *outid = out->batCacheid;
+ BBPkeepref(out);
+ return MAL_SUCCEED;
+ }
+
+ if ((b = BATdescriptor(*bid)) == NULL)
+ throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
+ if (sid && !is_bat_nil(*sid) && !(s = BATdescriptor(*sid))) {
+ BBPunfix(b->batCacheid);
+ throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
+ }
+ if ((out = COLnew(0, ATOMindex("oid"), ci.ncand, TRANSIENT)) == NULL) {
+ BBPunfix(b->batCacheid);
+ if (s)
+ BBPunfix(s->batCacheid);
+ throw(MAL, name, SQLSTATE(HY013) MAL_MALLOC_FAIL);
+ }
+
+ //Calculate the MBR for the constant geometry
+ mbr *const_mbr = NULL;
+ wkbMBR(&const_mbr,&wkb_const);
+
+ //Get a candidate list from searching on the rtree with the constant mbr
+ BUN* results_rtree = RTREEsearch(b,(mbr_t*)const_mbr, b->batCount);
+
+ canditer_init(&ci, b, s);
+ b_iter = bat_iterator(b);
+
+ //Intersect prev_cands with rtree_cands
+ //Cycle through rtree_cands
+ // if there is prev_cands -> bin search
+ // if there is not -> just cycle through rtree_cands
+
+ for (BUN i = 0; i < ci.ncand; i++) {
+ oid c_oid = canditer_next(&ci);
+
+ int i = 0;
+ while (results_rtree[i] != BUN_NONE && results_rtree[i] !=
c_oid) {
+ i++;
+ }
+ if (results_rtree[i] == BUN_NONE)
+ continue;
+
+ const wkb *col_wkb = BUNtvar(b_iter, c_oid - b->hseqbase);
+ if ((col_geom = wkb2geos(col_wkb)) == NULL)
+ continue;
+ if (GEOSGetSRID(col_geom) != GEOSGetSRID(const_geom)) {
+ GEOSGeom_destroy(col_geom);
+ GEOSGeom_destroy(const_geom);
+ bat_iterator_end(&b_iter);
+ BBPunfix(b->batCacheid);
+ if (s)
+ BBPunfix(s->batCacheid);
+ BBPreclaim(out);
+ throw(MAL, name, SQLSTATE(38000) "Geometries of
different SRID");
+ }
+
+ //GEOS functino returns 1 on true, 0 on false and 2 on exception
+ char ret = ((*func)(col_geom, const_geom));
+ bit cond = (ret == 1);
+ if (cond != anti) {
+ if (BUNappend(out, &c_oid, false) != GDK_SUCCEED) {
+ if (col_geom)
+ GEOSGeom_destroy(col_geom);
+ if (const_geom)
+ GEOSGeom_destroy(const_geom);
+ bat_iterator_end(&b_iter);
+ BBPunfix(b->batCacheid);
+ if (s)
+ BBPunfix(s->batCacheid);
+ BBPreclaim(out);
+ throw(MAL, name, SQLSTATE(HY013)
MAL_MALLOC_FAIL);
+ }
+ }
+ //TODO Deal with exception?
+ GEOSGeom_destroy(col_geom);
+ }
+ GEOSGeom_destroy(const_geom);
+ bat_iterator_end(&b_iter);
+ BBPunfix(b->batCacheid);
+ if (s)
+ BBPunfix(s->batCacheid);
+ *outid = out->batCacheid;
+ BBPkeepref(out);
+ return MAL_SUCCEED;
+}
+
+str
+wkbIntersectsSelect(bat* outid, const bat *bid , const bat *sid, wkb
**wkb_const, bit *anti) {
+ return
filterSelectGeomGeomToBitIndex(outid,bid,sid,*wkb_const,*anti,GEOSIntersects,"geom.wkbIntersectsSelect");
+}
+
+static str
+filterJoinGeomGeomDoubleToBit(bat *lres_id, bat *rres_id, const bat *l_id,
const bat *r_id, double double_flag, const bat *ls_id, const bat *rs_id, bit
nil_matches, lng *estimate, char (*func) (const GEOSGeometry *, const
GEOSGeometry *, double), const char *name)
+{
+ BAT *lres = NULL, *rres = NULL, *l = NULL, *r = NULL, *ls = NULL, *rs =
NULL;
+ BATiter l_iter, r_iter;
+ str msg = MAL_SUCCEED;
+ struct canditer l_ci, r_ci;
+ GEOSGeom l_geom, r_geom;
+ GEOSGeom *l_geoms = NULL, *r_geoms = NULL;
+ bool anti = false;
+
+ //get the input BATs
+ if ((l = BATdescriptor(*l_id)) == NULL || (r = BATdescriptor(*r_id)) ==
NULL) {
+ if (l)
+ BBPunfix(l->batCacheid);
+ if (r)
+ BBPunfix(r->batCacheid);
+ throw(MAL, name, SQLSTATE(HY002) RUNTIME_OBJECT_MISSING);
+ }
+ //get the candidate lists
+ if (ls_id && !is_bat_nil(*ls_id) && !(ls = BATdescriptor(*ls_id)) &&
rs_id && !is_bat_nil(*rs_id) && !(rs = BATdescriptor(*rs_id))) {
+ msg = createException(MAL, name, SQLSTATE(HY002)
RUNTIME_OBJECT_MISSING);
+ goto free;
+ }
+ canditer_init(&l_ci, l, ls);
+ canditer_init(&r_ci, r, rs);
+ //create new BATs for the output
+ if (is_lng_nil(*estimate) || *estimate == 0)
_______________________________________________
checkin-list mailing list -- [email protected]
To unsubscribe send an email to [email protected]
