first mulitppolygon builder

cache/db.go

@@ -315,26 +315,25 @@ func (p *NodesCache) GetNode(id int64) (*element.Node, error) {
 }
 
 func (p *NodesCache) Iter() chan *element.Node {
-    node := make(chan *element.Node)
-    go func() {
-        ro := levigo.NewReadOptions()
-        ro.SetFillCache(false)
-        it := p.db.NewIterator(ro)
-        defer it.Close()
-        it.SeekToFirst()
-        for ; it.Valid(); it.Next() {
-            nodes, err := binary.UnmarshalNode(it.Value())
-            if err != nil {
-                panic(err)
-            }
-            node <- nodes
-        }
-        close(node)
-    }()
-    return node
+	node := make(chan *element.Node)
+	go func() {
+		ro := levigo.NewReadOptions()
+		ro.SetFillCache(false)
+		it := p.db.NewIterator(ro)
+		defer it.Close()
+		it.SeekToFirst()
+		for ; it.Valid(); it.Next() {
+			nodes, err := binary.UnmarshalNode(it.Value())
+			if err != nil {
+				panic(err)
+			}
+			node <- nodes
+		}
+		close(node)
+	}()
+	return node
 }
 
-
 func (p *WaysCache) PutWay(way *element.Way) error {
 	keyBuf := idToKeyBuf(way.Id)
 	data, err := binary.MarshalWay(way)
@@ -416,6 +415,23 @@ func (p *WaysCache) Iter() chan *element.Way {
 	return way
 }
 
+func (self *WaysCache) FillMembers(members []element.Member) bool {
+	if members == nil || len(members) == 0 {
+		return false
+	}
+	for i, member := range members {
+		if member.Type != element.WAY {
+			continue
+		}
+		way, err := self.GetWay(member.Id)
+		if err != nil || way == nil {
+			return false
+		}
+		members[i].Way = way
+	}
+	return true
+}
+
 func (p *RelationsCache) PutRelation(relation *element.Relation) error {
 	keyBuf := idToKeyBuf(relation.Id)
 	data, err := binary.MarshalRelation(relation)

element/element.go

@@ -45,6 +45,7 @@ type Member struct {
 	Id   int64
 	Type MemberType
 	Role string
+	Way  *Way
 }
 
 type Relation struct {

geom/geom.go

@@ -77,6 +77,21 @@ func LineStringWKB(g *geos.GEOS, nodes []element.Node) (*element.Geometry, error
 }
 
 func PolygonWKB(g *geos.GEOS, nodes []element.Node) (*element.Geometry, error) {
+	geom, err := Polygon(g, nodes)
+	if err != nil {
+		return nil, err
+	}
+	wkb := g.AsWKB(geom)
+	if wkb == nil {
+		return nil, errors.New("could not create wkb")
+	}
+	return &element.Geometry{
+		Wkb:  wkb,
+		Geom: geom,
+	}, nil
+}
+
+func Polygon(g *geos.GEOS, nodes []element.Node) (*geos.Geom, error) {
 	coordSeq, err := g.CreateCoordSeq(uint32(len(nodes)), 2)
 	if err != nil {
 		return nil, err
@@ -93,18 +108,11 @@ func PolygonWKB(g *geos.GEOS, nodes []element.Node) (*element.Geometry, error) {
 		return nil, err
 	}
 	// geom inherited by Polygon, no destroy
+
 	geom = g.CreatePolygon(geom, nil)
 	if err != nil {
 		return nil, err
 	}
-	wkb := g.AsWKB(geom)
-	if wkb == nil {
-		g.Destroy(geom)
-		return nil, errors.New("could not create wkb")
-	}
 	g.DestroyLater(geom)
-	return &element.Geometry{
-		Wkb:  wkb,
-		Geom: geom,
-	}, nil
+	return geom, nil
 }

geom/geom_test.go

@@ -13,12 +13,12 @@ func TestLineString(t *testing.T) {
 	nodes[1] = element.Node{Lat: 0, Long: 10}
 	g := geos.NewGEOS()
 	defer g.Finish()
-	wkt, err := LineStringWKB(g, nodes)
+	geom, err := LineStringWKB(g, nodes)
 	if err != nil {
 		t.Fatal(err)
 	}
-	if bytes.Compare(wkt[0:2], []byte{0x1, 0x2}) != 0 {
-		t.Errorf("%#v", wkt)
+	if bytes.Compare(geom.Wkb[0:2], []byte{0x1, 0x2}) != 0 {
+		t.Errorf("%#v", geom.Wkb)
 	}
 }
 
@@ -31,13 +31,13 @@ func TestPolygon(t *testing.T) {
 	}
 	g := geos.NewGEOS()
 	defer g.Finish()
-	wkt, err := PolygonWKB(g, nodes)
+	geom, err := PolygonWKB(g, nodes)
 	if err != nil {
 		t.Fatal(err)
 	}
 
-	if bytes.Compare(wkt[0:2], []byte{0x1, 0x3}) != 0 {
-		t.Errorf("%#v", wkt)
+	if bytes.Compare(geom.Wkb[0:2], []byte{0x1, 0x3}) != 0 {
+		t.Errorf("%#v", geom.Wkb)
 	}
 }
 

geom/geos/geos.go

@@ -134,6 +134,55 @@ 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) Contains(a, b *Geom) bool {
+	result := C.GEOSContains_r(this.v, a.v, b.v)
+	if result == 1 {
+		return true
+	}
+	// result == 2 -> exception (already logged to console)
+	return false
+}
+
+func (this *GEOS) ExteriorRing(geom *Geom) *Geom {
+	ring := C.GEOSGetExteriorRing_r(this.v, geom.v)
+	if ring == nil {
+		return nil
+	}
+	return &Geom{ring}
+}
+
+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))
+		if geom == nil {
+			return nil
+		}
+		return &Geom{geom}
+	}
+
+	interiorPtr := make([]*C.GEOSGeometry, len(interiors))
+	for i, geom := range interiors {
+		interiorPtr[i] = geom.v
+	}
+	geom := C.GEOSGeom_createPolygon_r(this.v, exterior.v, &interiorPtr[0], C.uint(len(interiors)))
+	if geom == nil {
+		return nil
+	}
+	return &Geom{geom}
+}
+
+func (this *GEOS) MultiPolygon(polygons []*Geom) *Geom {
+	polygonPtr := make([]*C.GEOSGeometry, len(polygons))
+	for i, geom := range polygons {
+		polygonPtr[i] = geom.v
+	}
+	geom := C.GEOSGeom_createCollection_r(this.v, C.GEOS_MULTIPOLYGON, &polygonPtr[0], C.uint(len(polygons)))
+	if geom == nil {
+		return nil
+	}
+	return &Geom{geom}
+}
+
 func (this *GEOS) AsWKT(geom *Geom) string {
 	str := C.GEOSGeomToWKT_r(this.v, geom.v)
 	result := C.GoString(str)
@@ -160,7 +209,7 @@ func (this *GEOS) IsValid(geom *Geom) bool {
 
 func (this *Geom) Area() float64 {
 	var area C.double
-	if ret := C.GEOSArea(this.v, &area); ret == 0 {
+	if ret := C.GEOSArea(this.v, &area); ret == 1 {
 		return float64(area)
 	} else {
 		return 0

geom/multipolygon.go

@@ -0,0 +1,169 @@
+package geom
+
+import (
+	"errors"
+	"fmt"
+	"goposm/element"
+	"goposm/geom/geos"
+	"sort"
+)
+
+func BuildRelation(rel *element.Relation) error {
+	rings, err := BuildRings(rel)
+	if err != nil {
+		return err
+	}
+	geom, err := BuildGeometry(rings)
+	if err != nil {
+		return err
+	}
+	rel.Geom = &element.Geometry{Geom: geom}
+	return nil
+}
+
+func BuildRings(rel *element.Relation) ([]*Ring, error) {
+	var rings []*Ring
+	var incompleteRings []*Ring
+	var completeRings []*Ring
+	var err error
+
+	// create rings for all WAY members
+	for _, member := range rel.Members {
+		if member.Way == nil {
+			continue
+		}
+		rings = append(rings, NewRing(member.Way))
+	}
+
+	g := geos.NewGEOS()
+	defer g.Finish()
+
+	// create geometries for closed rings, collect incomplete rings
+	for _, r := range rings {
+		if r.IsClosed() {
+			r.geom, err = Polygon(g, r.nodes)
+			if err != nil {
+				return nil, err
+			}
+			completeRings = append(completeRings, r)
+		} else {
+			incompleteRings = append(incompleteRings, r)
+		}
+	}
+	// merge incomplete rings
+	mergedRings := mergeRings(incompleteRings)
+	if len(completeRings)+len(mergedRings) == 0 {
+		return nil, errors.New(
+			fmt.Sprintf("linestring from relation %d has no rings", rel.Id),
+		)
+	}
+	// create geometries for merged rings
+	for _, ring := range mergedRings {
+		if !ring.IsClosed() {
+			return nil, errors.New(
+				fmt.Sprintf("linestrings from relation %d do not form a ring", rel.Id),
+			)
+		}
+		ring.geom, err = Polygon(g, ring.nodes)
+		if err != nil {
+			return nil, err
+		}
+	}
+
+	completeRings = append(completeRings, mergedRings...)
+	return completeRings, nil
+}
+
+type SortableRingsDesc []*Ring
+
+func (r SortableRingsDesc) Len() int           { return len(r) }
+func (r SortableRingsDesc) Less(i, j int) bool { return r[i].area > r[j].area }
+func (r SortableRingsDesc) Swap(i, j int)      { r[i], r[j] = r[j], r[i] }
+
+func BuildGeometry(rings []*Ring) (*geos.Geom, error) {
+	g := geos.NewGEOS()
+	defer g.Finish()
+
+	// sort by area (large to small)
+	for _, r := range rings {
+		r.area = r.geom.Area()
+	}
+	sort.Sort(SortableRingsDesc(rings))
+
+	totalRings := len(rings)
+	shells := map[*Ring]bool{rings[0]: true}
+	for i := 0; i < totalRings; i++ {
+		testGeom := rings[i].geom //TODO prepared
+		for j := i + 1; j < totalRings; j++ {
+			if g.Contains(testGeom, rings[j].geom) {
+				if rings[j].containedBy != -1 {
+					// j is inside a larger ring, remove that relationship
+					// e.g. j is hole inside a hole (i)
+					delete(rings[rings[j].containedBy].holes, rings[j])
+					delete(shells, rings[j])
+				}
+				// remember parent
+				rings[j].containedBy = i
+				// add ring as hole or shell
+				if ringIsHole(rings, j) {
+					rings[i].holes[rings[j]] = true
+				} else {
+					shells[rings[j]] = true
+				}
+			}
+		}
+		if rings[i].containedBy == -1 {
+			// add as shell if it is not a hole
+			shells[rings[i]] = true
+		}
+	}
+
+	var polygons []*geos.Geom
+	for shell, _ := range shells {
+		var interiors []*geos.Geom
+		for hole, _ := range shell.holes {
+			ring := g.ExteriorRing(hole.geom)
+			if ring == nil {
+				return nil, errors.New("Error while getting exterior ring.")
+			}
+			interiors = append(interiors, ring)
+		}
+		exterior := g.ExteriorRing(shell.geom)
+		if exterior == nil {
+			return nil, errors.New("Error while getting exterior ring.")
+		}
+		polygon := g.Polygon(exterior, interiors)
+		if polygon == nil {
+			return nil, errors.New("Error while building polygon.")
+		}
+		polygons = append(polygons, polygon)
+	}
+	var result *geos.Geom
+
+	if len(polygons) == 1 {
+		result = polygons[0]
+	} else {
+		result = g.MultiPolygon(polygons)
+		if result == nil {
+			return nil, errors.New("Error while building multi-polygon.")
+		}
+	}
+
+	return result, nil
+}
+
+// ringIsHole returns true if rings[idx] is a hole, False if it is a
+// shell (also if hole in a hole, etc)
+func ringIsHole(rings []*Ring, idx int) bool {
+
+	containedCounter := 0
+	for {
+		idx = rings[idx].containedBy
+		if idx == -1 {
+
+			break
+		}
+		containedCounter += 1
+	}
+	return containedCounter%2 == 1
+}

geom/multipolygon_test.go

@@ -0,0 +1,192 @@
+package geom
+
+import (
+	"goposm/element"
+	"goposm/geom/geos"
+	"testing"
+)
+
+type coord struct {
+	id   int64
+	long float64
+	lat  float64
+}
+
+func makeWay(id int64, tags element.Tags, coords []coord) element.Way {
+	way := element.Way{}
+	way.Id = id
+	for _, coord := range coords {
+		way.Refs = append(way.Refs, coord.id)
+		way.Nodes = append(way.Nodes,
+			element.Node{OSMElem: element.OSMElem{Id: coord.id}, Long: coord.long, Lat: coord.lat})
+	}
+	return way
+}
+
+func TestMultiPolygonWithHole(t *testing.T) {
+	w1 := makeWay(1, nil, []coord{
+		{1, 0, 0},
+		{2, 10, 0},
+		{3, 10, 10},
+		{4, 0, 10},
+		{1, 0, 0},
+	})
+	w2 := makeWay(1, nil, []coord{
+		{5, 2, 2},
+		{6, 8, 2},
+		{7, 8, 8},
+		{8, 2, 8},
+		{5, 2, 2},
+	})
+
+	rel := element.Relation{OSMElem: element.OSMElem{Id: 1}}
+	rel.Members = []element.Member{
+		{1, element.WAY, "outer", &w1},
+		{2, element.WAY, "inner", &w2},
+	}
+
+	BuildRelation(&rel)
+	g := geos.NewGEOS()
+	defer g.Finish()
+
+	if !g.IsValid(rel.Geom.Geom) {
+		t.Fatal("geometry not valid", g.AsWKT(rel.Geom.Geom))
+	}
+
+	if area := rel.Geom.Geom.Area(); area != 64 {
+		t.Fatal("aread not 64", area)
+	}
+}
+
+func TestMultiPolygonWithMultipleHoles(t *testing.T) {
+	w1 := makeWay(1, nil, []coord{
+		{1, 0, 0},
+		{2, 10, 0},
+		{3, 10, 10},
+		{4, 0, 10},
+		{1, 0, 0},
+	})
+	w2 := makeWay(1, nil, []coord{
+		{1, 1, 1},
+		{2, 2, 1},
+		{3, 2, 2},
+		{4, 1, 2},
+		{1, 1, 1},
+	})
+	w3 := makeWay(3, nil, []coord{
+		{1, 3, 3},
+		{2, 4, 3},
+		{3, 4, 4},
+		{4, 3, 4},
+		{1, 3, 3},
+	})
+
+	rel := element.Relation{OSMElem: element.OSMElem{Id: 1}}
+	rel.Members = []element.Member{
+		{1, element.WAY, "outer", &w1},
+		{2, element.WAY, "inner", &w2},
+		{3, element.WAY, "inner", &w3},
+	}
+
+	BuildRelation(&rel)
+	g := geos.NewGEOS()
+	defer g.Finish()
+
+	if !g.IsValid(rel.Geom.Geom) {
+		t.Fatal("geometry not valid", g.AsWKT(rel.Geom.Geom))
+	}
+
+	if area := rel.Geom.Geom.Area(); area != 100-1-1 {
+		t.Fatal("area invalid", area)
+	}
+}
+
+func TestMultiPolygonWithNeastedHoles(t *testing.T) {
+	w1 := makeWay(1, nil, []coord{
+		{1, 0, 0},
+		{2, 10, 0},
+		{3, 10, 10},
+		{4, 0, 10},
+		{1, 0, 0},
+	})
+	w2 := makeWay(2, nil, []coord{
+		{1, 1, 1},
+		{2, 9, 1},
+		{3, 9, 9},
+		{4, 1, 9},
+		{1, 1, 1},
+	})
+	w3 := makeWay(3, nil, []coord{
+		{1, 2, 2},
+		{2, 8, 2},
+		{3, 8, 8},
+		{4, 2, 8},
+		{1, 2, 2},
+	})
+	w4 := makeWay(4, nil, []coord{
+		{1, 3, 3},
+		{2, 7, 3},
+		{3, 7, 7},
+		{4, 3, 7},
+		{1, 3, 3},
+	})
+	w5 := makeWay(5, nil, []coord{
+		{1, 4, 4},
+		{2, 6, 4},
+		{3, 6, 6},
+		{4, 4, 6},
+		{1, 4, 4},
+	})
+
+	rel := element.Relation{OSMElem: element.OSMElem{Id: 1}}
+	rel.Members = []element.Member{
+		{1, element.WAY, "outer", &w1},
+		{2, element.WAY, "inner", &w2},
+		{3, element.WAY, "inner", &w3},
+		{4, element.WAY, "inner", &w4},
+		{5, element.WAY, "inner", &w5},
+	}
+
+	BuildRelation(&rel)
+	g := geos.NewGEOS()
+	defer g.Finish()
+
+	if !g.IsValid(rel.Geom.Geom) {
+		t.Fatal("geometry not valid", g.AsWKT(rel.Geom.Geom))
+	}
+
+	if area := rel.Geom.Geom.Area(); area != 100-64+36-16+4 {
+		t.Fatal("area invalid", area)
+	}
+}
+
+func TestPolygonFromTwoWays(t *testing.T) {
+	w1 := makeWay(1, nil, []coord{
+		{1, 0, 0},
+		{2, 10, 0},
+		{3, 10, 10},
+	})
+	w2 := makeWay(2, nil, []coord{
+		{3, 10, 10},
+		{4, 0, 10},
+		{1, 0, 0},
+	})
+
+	rel := element.Relation{OSMElem: element.OSMElem{Id: 1}}
+	rel.Members = []element.Member{
+		{1, element.WAY, "outer", &w1},
+		{2, element.WAY, "inner", &w2},
+	}
+
+	BuildRelation(&rel)
+	g := geos.NewGEOS()
+	defer g.Finish()
+
+	if !g.IsValid(rel.Geom.Geom) {
+		t.Fatal("geometry not valid", g.AsWKT(rel.Geom.Geom))
+	}
+
+	if area := rel.Geom.Geom.Area(); area != 100 {
+		t.Fatal("area invalid", area)
+	}
+}

geom/ring.go

@@ -2,12 +2,33 @@ package geom
 
 import (
 	"goposm/element"
+	"goposm/geom/geos"
 )
 
 type Ring struct {
-	ways  []*element.Way
-	refs  []int64
-	nodes []element.Node
+	ways        []*element.Way
+	refs        []int64
+	nodes       []element.Node
+	geom        *geos.Geom
+	holes       map[*Ring]bool
+	containedBy int
+	area        float64
+}
+
+func (r *Ring) IsClosed() bool {
+	return len(r.refs) >= 4 && r.refs[0] == r.refs[len(r.refs)-1]
+}
+
+func NewRing(way *element.Way) *Ring {
+	ring := Ring{}
+	ring.ways = []*element.Way{way}
+	ring.refs = make([]int64, len(way.Refs))
+	ring.nodes = make([]element.Node, len(way.Nodes))
+	ring.containedBy = -1
+	ring.holes = make(map[*Ring]bool)
+	copy(ring.refs, way.Refs)
+	copy(ring.nodes, way.Nodes)
+	return &ring
 }
 
 func reverseRefs(refs []int64) {

geom/ring_test.go

@@ -10,23 +10,23 @@ func TestRingMerge(t *testing.T) {
 	w1 := element.Way{}
 	w1.Id = 1
 	w1.Refs = []int64{1, 2, 3}
-	nodes := []element.Node{
+	w1.Nodes = []element.Node{
 		element.Node{},
 		element.Node{},
 		element.Node{},
 	}
-	r1 := Ring{[]*element.Way{&w1}, w1.Refs, nodes}
+	r1 := NewRing(&w1)
 
 	w2 := element.Way{}
 	w2.Id = 2
 	w2.Refs = []int64{3, 4, 1}
-	nodes = []element.Node{
+	w2.Nodes = []element.Node{
 		element.Node{},
 		element.Node{},
 		element.Node{},
 	}
-	r2 := Ring{[]*element.Way{&w2}, w2.Refs, nodes}
-	rings := []*Ring{&r1, &r2}
+	r2 := NewRing(&w2)
+	rings := []*Ring{r1, r2}
 
 	result := mergeRings(rings)
 	if len(result) != 1 {
@@ -45,35 +45,35 @@ func TestRingMergeReverseEndpoints(t *testing.T) {
 	w1 := element.Way{}
 	w1.Id = 1
 	w1.Refs = []int64{1, 2, 3, 4}
-	nodes := []element.Node{
+	w1.Nodes = []element.Node{
 		element.Node{},
 		element.Node{},
 		element.Node{},
 		element.Node{},
 	}
-	r1 := Ring{[]*element.Way{&w1}, w1.Refs, nodes}
+	r1 := NewRing(&w1)
 
 	w2 := element.Way{}
 	w2.Id = 2
 	w2.Refs = []int64{6, 5, 4}
-	nodes = []element.Node{
+	w2.Nodes = []element.Node{
 		element.Node{},
 		element.Node{},
 		element.Node{},
 	}
-	r2 := Ring{[]*element.Way{&w2}, w2.Refs, nodes}
+	r2 := NewRing(&w2)
 
 	w3 := element.Way{}
 	w3.Id = 3
 	w3.Refs = []int64{1, 7, 6}
-	nodes = []element.Node{
+	w3.Nodes = []element.Node{
 		element.Node{},
 		element.Node{},
 		element.Node{},
 	}
-	r3 := Ring{[]*element.Way{&w3}, w3.Refs, nodes}
+	r3 := NewRing(&w3)
 
-	rings := []*Ring{&r1, &r2, &r3}
+	rings := []*Ring{r1, r2, r3}
 
 	result := mergeRings(rings)
 	if len(result) != 1 {
@@ -140,10 +140,10 @@ func TestRingMergePermutations(t *testing.T) {
 			w4.Nodes = []element.Node{element.Node{}, element.Node{}, element.Node{}, element.Node{}}
 
 			rings := []*Ring{
-				&Ring{[]*element.Way{&w1}, w1.Refs, w1.Nodes},
-				&Ring{[]*element.Way{&w2}, w2.Refs, w2.Nodes},
-				&Ring{[]*element.Way{&w3}, w3.Refs, w3.Nodes},
-				&Ring{[]*element.Way{&w4}, w4.Refs, w4.Nodes},
+				&Ring{ways: []*element.Way{&w1}, refs: w1.Refs, nodes: w1.Nodes},
+				&Ring{ways: []*element.Way{&w2}, refs: w2.Refs, nodes: w2.Nodes},
+				&Ring{ways: []*element.Way{&w3}, refs: w3.Refs, nodes: w3.Nodes},
+				&Ring{ways: []*element.Way{&w4}, refs: w4.Refs, nodes: w4.Nodes},
 			}
 			result := mergeRings(rings)
 			if len(result) != 1 {

goposm.go

@@ -2,6 +2,7 @@ package main
 
 import (
 	"flag"
+	"fmt"
 	"goposm/cache"
 	"goposm/database"
 	_ "goposm/database/postgis"
@@ -232,8 +233,19 @@ func main() {
 	if *write {
 		progress.Reset()
 		rel := osmCache.Relations.Iter()
-		for _ = range rel {
+		for r := range rel {
 			progress.AddRelations(1)
+			if !osmCache.Ways.FillMembers(r.Members) {
+				fmt.Println("missing ways")
+			}
+			for _, m := range r.Members {
+				if m.Way == nil {
+					continue
+				}
+				if !osmCache.Coords.FillWay(m.Way) {
+					fmt.Println("missing nodes", m.Way)
+				}
+			}
 			// fmt.Println(r)
 		}