add first geometry clipper

master
Oliver Tonnhofer 2013-05-27 15:26:42 +02:00
parent 40b4fffc4f
commit 35933e85ce
4 changed files with 723 additions and 11 deletions

265
geom/clipper/clipper.go Normal file
View File

@ -0,0 +1,265 @@
package clipper
import (
"errors"
"goposm/geom/geos"
"goposm/geom/ogr"
"math"
"strings"
)
// Tile bbox into multiple sub-boxes, each of `width` size.
// >>> list(tile_bbox((-1, 1, 0.49, 1.51), 0.5)) #doctest: +NORMALIZE_WHITESPACE
// [(-1.0, 1.0, -0.5, 1.5),
// (-1.0, 1.5, -0.5, 2.0),
// (-0.5, 1.0, 0.0, 1.5),
// (-0.5, 1.5, 0.0, 2.0),
// (0.0, 1.0, 0.5, 1.5),
// (0.0, 1.5, 0.5, 2.0)]
func tileBounds(bounds geos.Bounds, width float64) []geos.Bounds {
var results []geos.Bounds
minX := math.Floor(bounds.MinX/width) * width
minY := math.Floor(bounds.MinY/width) * width
maxX := math.Ceil(bounds.MaxX/width) * width
maxY := math.Ceil(bounds.MaxY/width) * width
xSteps := math.Ceil((maxX - minX) / width)
ySteps := math.Ceil((maxY - minY) / width)
for x := 0; x < int(xSteps); x++ {
for y := 0; y < int(ySteps); y++ {
results = append(results, geos.Bounds{
minX + float64(x)*width,
minY + float64(y)*width,
minX + float64(x+1)*width,
minY + float64(y+1)*width,
})
}
}
return results
}
func SplitPolygonAtGrid(g *geos.Geos, geom *geos.Geom, gridWidth, currentGridWidth float64) ([]*geos.Geom, error) {
// >>> p = list(split_polygon_at_grid(geometry.box(-0.5, 1, 0.2, 2), 1))
// >>> p[0].contains(geometry.box(-0.5, 1, 0, 2))
// True
// >>> p[0].area == geometry.box(-0.5, 1, 0, 2).area
// True
// >>> p[1].contains(geometry.box(0, 1, 0.2, 2))
// True
// >>> p[1].area == geometry.box(0, 1, 0.2, 2).area
// True
// if not geom.is_valid:
// geom = geom.buffer(0)
var result []*geos.Geom
geomBounds := geom.Bounds()
if geomBounds == geos.NilBounds {
return nil, errors.New("couldn't create bounds for geom")
}
for _, bounds := range tileBounds(geom.Bounds(), currentGridWidth) {
clipGeom := g.BoundsPolygon(bounds)
if clipGeom == nil {
return nil, errors.New("couldn't create bounds polygon")
}
part := g.Intersection(geom, clipGeom)
if part == nil {
return nil, errors.New("couldn't create intersection")
}
if !g.IsEmpty(part) && strings.HasPrefix(g.Type(part), "Polygon") {
if gridWidth >= currentGridWidth {
result = append(result, part)
} else {
moreParts, err := SplitPolygonAtGrid(g, part, gridWidth, currentGridWidth/10.0)
if err != nil {
return nil, err
}
result = append(result, moreParts...)
}
}
}
return result, nil
}
type Clipper struct {
index *geos.Index
}
func NewFromOgrSource(source string) (*Clipper, error) {
ds, err := ogr.Open(source)
if err != nil {
return nil, err
}
g := geos.NewGeos()
defer g.Finish()
layer, err := ds.Layer()
if err != nil {
return nil, err
}
index := g.CreateIndex()
for geom := range layer.Geoms() {
parts, err := SplitPolygonAtGrid(g, geom, 10000, 10000*100)
if err != nil {
return nil, err
}
for _, part := range parts {
g.IndexAdd(index, part)
}
}
return &Clipper{index}, nil
}
func filterGeometryByType(g *geos.Geos, geom *geos.Geom, targetType string) []*geos.Geom {
// Filter (multi)geometry for compatible `geom_type`,
// because we can't insert points into linestring tables for example
geomType := g.Type(geom)
if geomType == targetType {
// same type is fine
return []*geos.Geom{geom}
}
if geomType == "Polygon" && targetType == "MultiPolygon" {
// multipolygon mappings also support polygons
return []*geos.Geom{geom}
}
if geomType == "MultiPolygon" && targetType == "Polygon" {
// polygon mappings should also support multipolygons
return []*geos.Geom{geom}
}
if g.NumGeoms(geom) >= 1 {
// GeometryCollection or MultiLineString? return list of geometries
var geoms []*geos.Geom
for _, part := range g.Geoms(geom) {
// only parts with same type
if g.Type(part) == targetType {
geoms = append(geoms, part)
}
}
if len(geoms) != 0 {
return geoms
}
}
return []*geos.Geom{}
}
func (clipper *Clipper) clip(geom *geos.Geom) ([]*geos.Geom, error) {
g := geos.NewGeos()
defer g.Finish()
hits := g.IndexQuery(clipper.index, geom)
if len(hits) == 0 {
return nil, nil
}
geomType := g.Type(geom)
var intersections []*geos.Geom
for _, hit := range hits {
if g.Contains(hit, geom) {
return []*geos.Geom{geom}, nil
}
if g.Intersects(hit, geom) {
newPart := g.Intersection(hit, geom)
newParts := filterGeometryByType(g, newPart, geomType)
for _, p := range newParts {
intersections = append(intersections, p)
}
}
}
return mergeGeometries(g, intersections, geomType), nil
// if not intersections:
// raise EmtpyGeometryError('No intersection or empty geometry')
// # intersections from multiple sub-polygons
// # try to merge them back to a single geometry
// try:
// if geom.type.endswith('Polygon'):
// union = cascaded_union(list(flatten_polygons(intersections)))
// elif geom.type.endswith('LineString'):
// linestrings = flatten_linestrings(intersections)
// linestrings = list(filter_invalid_linestrings(linestrings))
// if not linestrings:
// raise EmtpyGeometryError()
// union = linemerge(linestrings)
// if union.type == 'MultiLineString':
// union = list(union.geoms)
// elif geom.type == 'Point':
// union = intersections[0]
// else:
// log.warn('unexpexted geometry type %s', geom.type)
// raise EmtpyGeometryError()
// except ValueError, ex:
// # likely an 'No Shapely geometry can be created from null value' error
// log.warn('could not create union: %s', ex)
// raise EmtpyGeometryError()
// return union
}
func flattenPolygons(g *geos.Geos, geoms []*geos.Geom) []*geos.Geom {
var result []*geos.Geom
for _, geom := range geoms {
if g.Type(geom) == "MultiPolygon" {
result = append(result, g.Geoms(geom)...)
} else {
result = append(result, geom)
}
}
return result
}
func flattenLineStrings(g *geos.Geos, geoms []*geos.Geom) []*geos.Geom {
var result []*geos.Geom
for _, geom := range geoms {
if g.Type(geom) == "MultiLineString" {
result = append(result, g.Geoms(geom)...)
} else {
result = append(result, geom)
}
}
return result
}
func filterInvalidLineStrings(g *geos.Geos, geoms []*geos.Geom) []*geos.Geom {
var result []*geos.Geom
for _, geom := range geoms {
if geom.Length() > 1e-9 {
result = append(result, geom)
}
}
return result
}
// mergeGeometries
func mergeGeometries(g *geos.Geos, geoms []*geos.Geom, geomType string) []*geos.Geom {
// intersections from multiple sub-polygons
// try to merge them back to a single geometry
if strings.HasSuffix(geomType, "Polygon") {
polygons := flattenPolygons(g, geoms)
polygon := g.UnionPolygons(polygons)
return []*geos.Geom{polygon}
} else if strings.HasSuffix(geomType, "LineString") {
linestrings := flattenLineStrings(g, geoms)
linestrings = filterInvalidLineStrings(g, linestrings)
if len(linestrings) == 0 {
return nil
}
union := g.LineMerge(linestrings)
return union
} else if geomType == "Point" {
return geoms[0:1]
} else {
panic("unexpected geometry type" + geomType)
}
}

View File

@ -0,0 +1,301 @@
package clipper
import (
"goposm/geom/geos"
"testing"
)
func TestTileBounds(t *testing.T) {
expected := []geos.Bounds{
{-1.0, 1.0, -0.5, 1.5},
{-1.0, 1.5, -0.5, 2.0},
{-0.5, 1.0, 0.0, 1.5},
{-0.5, 1.5, 0.0, 2.0},
{0.0, 1.0, 0.5, 1.5},
{0.0, 1.5, 0.5, 2.0},
}
for i, bounds := range tileBounds(geos.Bounds{-1, 1, 0.49, 1.51}, 0.5) {
if expected[i] != bounds {
t.Fatalf("%v != %v\n", expected[i], bounds)
}
}
}
func TestSplitPolygonAtGrids(t *testing.T) {
expected := []geos.Bounds{
{0, 0, 0.05, 0.05},
{0, 0.05, 0.05, 0.1},
{0, 0.1, 0.05, 0.11},
{0.05, 0, 0.1, 0.05},
{0.05, 0.05, 0.1, 0.1},
{0.05, 0.1, 0.1, 0.11},
{0.1, 0, 0.15, 0.05},
{0.1, 0.05, 0.15, 0.1},
{0.1, 0.1, 0.15, 0.11},
}
g := geos.NewGeos()
defer g.Finish()
geom := g.BoundsPolygon(geos.Bounds{0, 0, 0.15, 0.11})
geoms, _ := SplitPolygonAtGrid(g, geom, 0.05, 5.0)
for i, geom := range geoms {
if expected[i] != geom.Bounds() {
t.Fatalf("%v != %v\n", expected[i], geom.Bounds())
}
}
}
func TestMergePolygonGeometries(t *testing.T) {
g := geos.NewGeos()
defer g.Finish()
// check non intersecting polygons
// should return multipolygon
geoms := []*geos.Geom{
g.BoundsPolygon(geos.Bounds{0, 0, 10, 10}),
g.BoundsPolygon(geos.Bounds{20, 20, 30, 30}),
}
result := mergeGeometries(g, geoms, "Polygon")
if len(result) != 1 {
t.Fatal("not a single geometrie")
}
if g.Type(result[0]) != "MultiPolygon" {
t.Fatal("not a multipolygon")
}
if !g.IsValid(result[0]) {
t.Fatal("not valid")
}
// check intersecting polygons
// should return single polygon
geoms = []*geos.Geom{
g.BoundsPolygon(geos.Bounds{0, 0, 10, 10}),
g.BoundsPolygon(geos.Bounds{5, 5, 30, 30}),
}
result = mergeGeometries(g, geoms, "Polygon")
if len(result) != 1 {
t.Fatal("not a single geometrie")
}
if g.Type(result[0]) != "Polygon" {
t.Fatal("not a polygon")
}
if !g.IsValid(result[0]) {
t.Fatal("not valid")
}
// same with multipolygon type
geoms = []*geos.Geom{
g.BoundsPolygon(geos.Bounds{0, 0, 10, 10}),
g.BoundsPolygon(geos.Bounds{5, 5, 30, 30}),
}
result = mergeGeometries(g, geoms, "MultiPolygon")
if len(result) != 1 {
t.Fatal("not a single geometrie")
}
if g.Type(result[0]) != "Polygon" {
t.Fatal("not a polygon")
}
if !g.IsValid(result[0]) {
t.Fatal("not valid")
}
// strip non Polygons
geoms = []*geos.Geom{
g.FromWkt("POINT(0 0)"),
g.BoundsPolygon(geos.Bounds{0, 0, 10, 10}),
g.FromWkt("LINESTRING(0 0, 0 10)"),
g.BoundsPolygon(geos.Bounds{5, 5, 30, 30}),
}
result = mergeGeometries(g, geoms, "Polygon")
if len(result) != 1 {
t.Fatal("not a single geometrie")
}
if g.Type(result[0]) != "Polygon" {
t.Fatal("not a polygon")
}
if !g.IsValid(result[0]) {
t.Fatal("not valid")
}
}
func TestMergeLineStringGeometries(t *testing.T) {
g := geos.NewGeos()
defer g.Finish()
// check non intersecting linestrings
// should return slice of two linestrings
geoms := []*geos.Geom{
g.FromWkt("LINESTRING(0 0, 10 0)"),
g.FromWkt("LINESTRING(20 0, 30 0)"),
}
result := mergeGeometries(g, geoms, "LineString")
if len(result) != 2 {
t.Fatal("not two lines")
}
if g.Type(result[0]) != "LineString" || g.Type(result[1]) != "LineString" {
t.Fatal("not two lines")
}
if !g.IsValid(result[0]) || !g.IsValid(result[1]) {
t.Fatal("not valid")
}
// check intersecting linestrings
// should return slice of a single merged linestring
geoms = []*geos.Geom{
g.FromWkt("LINESTRING(0 0, 10 0)"),
g.FromWkt("LINESTRING(0 0, 0 10)"),
g.FromWkt("LINESTRING(10 0, 10 10)"),
}
result = mergeGeometries(g, geoms, "LineString")
if len(result) != 1 {
t.Fatal("not a single geometrie")
}
if g.Type(result[0]) != "LineString" {
t.Fatal("not a linestring")
}
if !g.IsValid(result[0]) {
t.Fatal("not valid")
}
if result[0].Length() != 30 {
t.Fatal("wrong length", result[0].Length())
}
// same but with multilinestring type
geoms = []*geos.Geom{
g.FromWkt("LINESTRING(0 0, 10 0)"),
g.FromWkt("MULTILINESTRING((0 0, 0 10), (10 0, 10 10))"),
}
result = mergeGeometries(g, geoms, "MultiLineString")
if len(result) != 1 {
t.Fatal("not a single geometrie")
}
if g.Type(result[0]) != "LineString" {
t.Fatal("not a linestring")
}
if !g.IsValid(result[0]) {
t.Fatal("not valid")
}
if result[0].Length() != 30 {
t.Fatal("wrong length", result[0].Length())
}
// strip non LineStrings and tiny LineStrings
geoms = []*geos.Geom{
g.FromWkt("POINT(0 0)"),
g.FromWkt("LINESTRING(0 0, 0 10)"),
g.FromWkt("LINESTRING(20 20, 20.00000000001 20)"), // tiny length
g.FromWkt("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))"),
}
result = mergeGeometries(g, geoms, "LineString")
if len(result) != 1 {
t.Fatal("not a single geometrie")
}
if g.Type(result[0]) != "LineString" {
t.Fatal("not a linestring")
}
if !g.IsValid(result[0]) {
t.Fatal("not valid")
}
}
func TestFilterGeometryByType(t *testing.T) {
g := geos.NewGeos()
defer g.Finish()
var result []*geos.Geom
// filtered out
result = filterGeometryByType(g, g.FromWkt("POINT(0 0)"), "Polygon")
if len(result) != 0 {
t.Fatal()
}
result = filterGeometryByType(g, g.FromWkt("POINT(0 0)"), "Point")
if len(result) != 1 {
t.Fatal()
}
// filtered out
result = filterGeometryByType(g, g.FromWkt("LINESTRING(0 0, 10 0)"), "Polygon")
if len(result) != 0 {
t.Fatal()
}
// polygon <-> multipolygon types are compatible in both directions
result = filterGeometryByType(g, g.FromWkt("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))"), "Polygon")
if len(result) != 1 {
t.Fatal()
}
result = filterGeometryByType(g, g.FromWkt("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))"), "MultiPolygon")
if len(result) != 1 {
t.Fatal()
}
result = filterGeometryByType(g, g.FromWkt("MULTIPOLYGON(((0 0, 10 0, 10 10, 0 10, 0 0)))"), "Polygon")
if len(result) != 1 {
t.Fatal()
}
result = filterGeometryByType(g, g.FromWkt("LINESTRING(0 0, 10 0)"), "LineString")
if len(result) != 1 {
t.Fatal()
}
// multilinestrings are split
result = filterGeometryByType(g, g.FromWkt("MULTILINESTRING((0 0, 10 0), (20 0, 30 0))"), "LineString")
if len(result) != 2 {
t.Fatal()
}
}
func TestClipper(t *testing.T) {
g := geos.NewGeos()
defer g.Finish()
clipper, err := NewFromOgrSource("./hamburg_clip.geojson")
if err != nil {
t.Fatal(err)
}
result, err := clipper.clip(g.FromWkt("POINT(0 0)"))
if err != nil || result != nil {
t.Fatal(err)
}
result, err = clipper.clip(g.FromWkt("POINT(1106543 7082055)"))
if err != nil {
t.Fatal(err)
}
if len(result) != 1 {
t.Fatal()
}
result, err = clipper.clip(g.FromWkt("LINESTRING(1106543 7082055, 1107105.2 7087540.0)"))
if err != nil {
t.Fatal(err)
}
if len(result) != 2 {
t.Fatal()
}
geom := g.FromWkt("POLYGON((1106543 7082055, 1107105.2 7087540.0, 1112184.9 7084424.5, 1106543 7082055))")
result, err = clipper.clip(geom)
if err != nil {
t.Fatal(err)
}
if len(result) != 1 {
t.Fatal()
}
if geom.Area() <= result[0].Area() {
t.Fatalf("%f <= %f", geom.Area(), result[0].Area())
}
}

View File

@ -0,0 +1,8 @@
{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "id": 1 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 1108158.325407607015222, 7089488.17479423712939 ], [ 1109408.259032631991431, 7091932.489438731223345 ], [ 1111213.718713223934174, 7093432.409788761287928 ], [ 1113241.388816042570397, 7094043.488449884578586 ], [ 1113741.362266052514315, 7094626.790808229707181 ], [ 1114907.966982742538676, 7094932.33013879135251 ], [ 1116574.545149442739785, 7094960.106441569514573 ], [ 1117046.74229667452164, 7092793.554824859835207 ], [ 1117352.281627236166969, 7090321.463877587579191 ], [ 1118324.452224477892742, 7088793.767224779352546 ], [ 1117380.057930014561862, 7087877.149233094416559 ], [ 1117907.807682802900672, 7086099.465855280868709 ], [ 1116574.545149442739785, 7085571.716102492064238 ], [ 1117324.505324458004907, 7083182.95406355522573 ], [ 1116241.229516102699563, 7080877.520932952873409 ], [ 1116741.202966112876311, 7079683.139913484454155 ], [ 1110602.640052100410685, 7079405.376885700970888 ], [ 1107408.3652325917501, 7078960.956041247583926 ], [ 1107408.3652325917501, 7078960.956041247583926 ], [ 1103158.590907506411895, 7079044.284949583001435 ], [ 1104464.07713808817789, 7082405.217585762031376 ], [ 1105241.813615881605074, 7083766.256421900354326 ], [ 1107297.26002147840336, 7083932.914238570258021 ], [ 1107491.694140926934779, 7084960.63744136877358 ], [ 1106741.733965911669657, 7087710.491416423581541 ], [ 1106741.733965911669657, 7087710.491416423581541 ], [ 1107241.707415921846405, 7088571.556802552193403 ], [ 1108158.325407607015222, 7089488.17479423712939 ] ] ] } }
]
}

View File

@ -141,6 +141,24 @@ func (this *Geos) Buffer(geom *Geom, size float64) *Geom {
return &Geom{C.GEOSBuffer_r(this.v, geom.v, C.double(size), 50)}
}
func (this *Geos) NumGeoms(geom *Geom) int32 {
count := int32(C.GEOSGetNumGeometries_r(this.v, geom.v))
return count
}
func (this *Geos) Geoms(geom *Geom) []*Geom {
count := this.NumGeoms(geom)
var result []*Geom
for i := 0; int32(i) < count; i++ {
part := C.GEOSGetGeometryN_r(this.v, geom.v, C.int(i))
if part == nil {
return nil
}
result = append(result, &Geom{part})
}
return result
}
func (this *Geos) Contains(a, b *Geom) bool {
result := C.GEOSContains_r(this.v, a.v, b.v)
if result == 1 {
@ -150,6 +168,47 @@ func (this *Geos) Contains(a, b *Geom) bool {
return false
}
func (this *Geos) Intersects(a, b *Geom) bool {
result := C.GEOSIntersects_r(this.v, a.v, b.v)
if result == 1 {
return true
}
// result == 2 -> exception (already logged to console)
return false
}
func (this *Geos) Intersection(a, b *Geom) *Geom {
result := C.GEOSIntersection_r(this.v, a.v, b.v)
if result == nil {
return nil
}
geom := &Geom{result}
this.DestroyLater(geom)
return geom
}
func (this *Geos) UnionPolygons(polygons []*Geom) *Geom {
multiPolygon := this.MultiPolygon(polygons)
result := C.GEOSUnaryUnion_r(this.v, multiPolygon.v)
if result == nil {
return nil
}
return &Geom{result}
}
func (this *Geos) LineMerge(lines []*Geom) []*Geom {
multiLineString := this.MultiLineString(lines)
result := C.GEOSLineMerge_r(this.v, multiLineString.v)
if result == nil {
return nil
}
geom := &Geom{result}
if this.Type(geom) == "LineString" {
return []*Geom{geom}
}
return this.Geoms(geom)
}
func (this *Geos) ExteriorRing(geom *Geom) *Geom {
ring := C.GEOSGetExteriorRing_r(this.v, geom.v)
if ring == nil {
@ -158,6 +217,41 @@ func (this *Geos) ExteriorRing(geom *Geom) *Geom {
return &Geom{ring}
}
func (this *Geos) BoundsPolygon(bounds Bounds) *Geom {
coordSeq, err := this.CreateCoordSeq(5, 2)
if err != nil {
return nil
}
// coordSeq inherited by LineString, no destroy
if err := coordSeq.SetXY(this, 0, bounds.MinX, bounds.MinY); err != nil {
return nil
}
if err := coordSeq.SetXY(this, 1, bounds.MaxX, bounds.MinY); err != nil {
return nil
}
if err := coordSeq.SetXY(this, 2, bounds.MaxX, bounds.MaxY); err != nil {
return nil
}
if err := coordSeq.SetXY(this, 3, bounds.MinX, bounds.MaxY); err != nil {
return nil
}
if err := coordSeq.SetXY(this, 4, bounds.MinX, bounds.MinY); err != nil {
return nil
}
geom, err := coordSeq.AsLinearRing(this)
if err != nil {
return nil
}
// geom inherited by Polygon, no destroy
geom = this.CreatePolygon(geom, nil)
this.DestroyLater(geom)
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))
@ -199,6 +293,17 @@ func (this *Geos) MultiPolygon(polygons []*Geom) *Geom {
}
return &Geom{geom}
}
func (this *Geos) MultiLineString(lines []*Geom) *Geom {
linePtr := make([]*C.GEOSGeometry, len(lines))
for i, geom := range lines {
linePtr[i] = geom.v
}
geom := C.GEOSGeom_createCollection_r(this.v, C.GEOS_MULTILINESTRING, &linePtr[0], C.uint(len(lines)))
if geom == nil {
return nil
}
return &Geom{geom}
}
func (this *Geos) AsWkt(geom *Geom) string {
str := C.GEOSGeomToWKT_r(this.v, geom.v)
@ -232,6 +337,22 @@ func (this *Geos) IsValid(geom *Geom) bool {
return false
}
func (this *Geos) IsEmpty(geom *Geom) bool {
if C.GEOSisEmpty_r(this.v, geom.v) == 1 {
return true
}
return false
}
func (this *Geos) Type(geom *Geom) string {
geomType := C.GEOSGeomType_r(this.v, geom.v)
if geomType == nil {
return "Unknown"
}
defer C.free(unsafe.Pointer(geomType))
return C.GoString(geomType)
}
func (this *Geom) Area() float64 {
var area C.double
if ret := C.GEOSArea(this.v, &area); ret == 1 {
@ -241,14 +362,25 @@ func (this *Geom) Area() float64 {
}
}
func (this *Geom) Bounds() *Bounds {
func (this *Geom) Length() float64 {
var length C.double
if ret := C.GEOSLength(this.v, &length); ret == 1 {
return float64(length)
} else {
return 0
}
}
var NilBounds = Bounds{1e20, 1e20, -1e20, -1e20}
func (this *Geom) Bounds() Bounds {
geom := C.GEOSEnvelope(this.v)
if geom == nil {
return nil
return NilBounds
}
extRing := C.GEOSGetExteriorRing(geom)
if extRing == nil {
return nil
return NilBounds
}
cs := C.GEOSGeom_getCoordSeq(extRing)
var csLen C.uint
@ -277,7 +409,7 @@ func (this *Geom) Bounds() *Bounds {
}
}
return &Bounds{minx, miny, maxx, maxy}
return Bounds{minx, miny, maxx, maxy}
}
type Bounds struct {
@ -314,7 +446,8 @@ func (this *Geos) DestroyCoordSeq(coordSeq *CoordSeq) {
}
type Index struct {
v *C.GEOSSTRtree
v *C.GEOSSTRtree
geoms []*Geom
}
func (this *Geos) CreateIndex() *Index {
@ -322,17 +455,18 @@ func (this *Geos) CreateIndex() *Index {
if tree == nil {
panic("unable to create tree")
}
return &Index{tree}
return &Index{tree, []*Geom{}}
}
// IndexQuery adds a geom to the index with the id.
func (this *Geos) IndexAdd(index *Index, id int, geom *Geom) {
func (this *Geos) IndexAdd(index *Index, geom *Geom) {
id := len(index.geoms)
C.IndexAdd(this.v, index.v, geom.v, C.size_t(id))
index.geoms = append(index.geoms, geom)
}
// IndexQuery queries the index for intersections with geom and
// returns a channel with ids of each intersected geometries.
func (this *Geos) IndexQuery(index *Index, geom *Geom) <-chan int {
// IndexQuery queries the index for intersections with geom.
func (this *Geos) IndexQuery(index *Index, geom *Geom) []*Geom {
hits := make(chan int)
go func() {
//
@ -342,7 +476,11 @@ func (this *Geos) IndexQuery(index *Index, geom *Geom) <-chan int {
C.IndexQuery(this.v, index.v, geom.v, unsafe.Pointer(&hits))
close(hits)
}()
return hits
var geoms []*Geom
for idx := range hits {
geoms = append(geoms, index.geoms[idx])
}
return geoms
}
//export goIndexSendQueryResult