new buffered clipper

geom/clipper/clipper.go

@@ -87,10 +87,16 @@ func SplitPolygonAtGrid(g *geos.Geos, geom *geos.Geom, gridWidth, currentGridWid
 }
 
 type Clipper struct {
-	index *geos.Index
+	index        *geos.Index
+	bufferedPrep *geos.PreparedGeom
+	bufferedBbox geos.Bounds
 }
 
 func NewFromOgrSource(source string) (*Clipper, error) {
+	return NewFromOgrSourceWithBuffered(source, 0.0)
+}
+
+func NewFromOgrSourceWithBuffered(source string, buffer float64) (*Clipper, error) {
 	ds, err := ogr.Open(source)
 	if err != nil {
 		return nil, err
@@ -106,7 +112,27 @@ func NewFromOgrSource(source string) (*Clipper, error) {
 
 	index := g.CreateIndex()
 
+	var polygons []*geos.Geom
+
+	withBuffer := false
+	if buffer != 0.0 {
+		withBuffer = true
+	}
+
 	for geom := range layer.Geoms() {
+		if withBuffer {
+			simplified := g.SimplifyPreserveTopology(geom, 1000)
+			if simplified == nil {
+				return nil, errors.New("couldn't simplify limitto")
+			}
+			buffered := g.Buffer(simplified, buffer)
+			g.Destroy(simplified)
+			if buffered == nil {
+				return nil, errors.New("couldn't buffer limitto")
+			}
+			g.DestroyLater(buffered)
+			polygons = append(polygons, buffered)
+		}
 		parts, err := SplitPolygonAtGrid(g, geom, 20000, 20000*100)
 		if err != nil {
 			return nil, err
@@ -115,7 +141,27 @@ func NewFromOgrSource(source string) (*Clipper, error) {
 			g.IndexAdd(index, part)
 		}
 	}
-	return &Clipper{index}, nil
+
+	var bbox geos.Bounds
+	var prep *geos.PreparedGeom
+	if len(polygons) > 0 {
+		union := g.UnionPolygons(polygons)
+		if union == nil {
+			return nil, errors.New("unable to union limitto polygons")
+		}
+		simplified := g.SimplifyPreserveTopology(union, 5000)
+		if simplified == nil {
+			return nil, errors.New("unable to simplify limitto polygons")
+		}
+		g.Destroy(union)
+		bbox = simplified.Bounds()
+		prep = g.Prepare(simplified)
+		// keep simplified around for prepared geometry
+		if prep == nil {
+			return nil, errors.New("unable to prepare limitto polygons")
+		}
+	}
+	return &Clipper{index, prep, bbox}, nil
 }
 
 func filterGeometryByType(g *geos.Geos, geom *geos.Geom, targetType string) []*geos.Geom {
@@ -193,12 +239,23 @@ func (clipper *Clipper) Clip(geom *geos.Geom) ([]*geos.Geom, error) {
 	return mergeGeometries(g, intersections, geomType), nil
 }
 
-func (clipper *Clipper) Intersects(g *geos.Geos, geom *geos.Geom) bool {
-	hits := g.IndexQuery(clipper.index, geom)
-	if len(hits) == 0 {
+func (c *Clipper) IntersectsBuffer(g *geos.Geos, x, y float64) bool {
+	if c.bufferedPrep == nil {
+		return true
+	}
+	if x < c.bufferedBbox.MinX ||
+		y < c.bufferedBbox.MinY ||
+		x > c.bufferedBbox.MaxX ||
+		y > c.bufferedBbox.MaxY {
 		return false
 	}
-	return true
+	p := g.Point(x, y)
+	if p == nil {
+		return false
+	}
+	defer g.Destroy(p)
+
+	return g.PreparedIntersects(c.bufferedPrep, p)
 }
 
 func flattenPolygons(g *geos.Geos, geoms []*geos.Geom) []*geos.Geom {

geom/clipper/clipper_test.go

@@ -299,6 +299,21 @@ func TestClipper(t *testing.T) {
 	}
 }
 
+func TestClipperWithBuffer(t *testing.T) {
+	g := geos.NewGeos()
+	defer g.Finish()
+	clipper, err := NewFromOgrSourceWithBuffered("./hamburg_clip.geojson", 10000.0)
+	if err != nil {
+		t.Fatal(err)
+	}
+	if clipper.IntersectsBuffer(g, 1106543, 7082055) != true {
+		t.Fatal()
+	}
+	if clipper.IntersectsBuffer(g, 1006543, 7082055) != false {
+		t.Fatal()
+	}
+}
+
 func BenchmarkClipper(b *testing.B) {
 	g := geos.NewGeos()
 	defer g.Finish()

geom/geom.go

@@ -29,15 +29,9 @@ var (
 )
 
 func Point(g *geos.Geos, node element.Node) (*geos.Geom, error) {
-	coordSeq, err := g.CreateCoordSeq(1, 2)
-	if err != nil {
-		return nil, err
-	}
-	// coordSeq inherited by LineString
-	coordSeq.SetXY(g, 0, node.Long, node.Lat)
-	geom, err := coordSeq.AsPoint(g)
-	if err != nil {
-		return nil, err
+	geom := g.Point(node.Long, node.Lat)
+	if geom == nil {
+		return nil, NewGeomError("couldn't create point")
 	}
 	g.DestroyLater(geom)
 	return geom, nil

geom/geos/geos.go

@@ -176,6 +176,20 @@ func (this *Geos) BoundsPolygon(bounds Bounds) *Geom {
 
 }
 
+func (this *Geos) Point(x, y float64) *Geom {
+	coordSeq, err := this.CreateCoordSeq(1, 2)
+	if err != nil {
+		return nil
+	}
+	// coordSeq inherited by LineString
+	coordSeq.SetXY(this, 0, x, y)
+	geom, err := coordSeq.AsPoint(this)
+	if err != nil {
+		return nil
+	}
+	return geom
+}
+
 func (this *Geos) Polygon(exterior *Geom, interiors []*Geom) *Geom {
 	if len(interiors) == 0 {
 		geom := C.GEOSGeom_createPolygon_r(this.v, exterior.v, nil, C.uint(0))

geom/geos/ops.go

@@ -42,6 +42,14 @@ func (this *Geos) Buffer(geom *Geom, size float64) *Geom {
 	return &Geom{buffered}
 }
 
+func (this *Geos) SimplifyPreserveTopology(geom *Geom, tolerance float64) *Geom {
+	simplified := C.GEOSTopologyPreserveSimplify_r(this.v, geom.v, C.double(tolerance))
+	if simplified == nil {
+		return nil
+	}
+	return &Geom{simplified}
+}
+
 // UnionPolygons tries to merge polygons.
 // Returns a single (Multi)Polygon.
 // Destroys polygons and returns new allocated (Multi)Polygon as necessary.