From ac094e8dce31cc6fe6893186fae4e9a6e73dc9ef Mon Sep 17 00:00:00 2001 From: Oliver Tonnhofer Date: Mon, 21 Oct 2013 14:59:31 +0200 Subject: [PATCH] replace OGR dependency for limitto with simple GeoJSON parser --- geom/geojson/doc.go | 4 + geom/geojson/geojson.go | 309 +++++++++++++++++++++++++++++++++++++++ geom/geom.go | 2 +- geom/geos/geos.go | 13 +- geom/limit/limit.go | 37 +++-- geom/limit/limit_test.go | 39 +++-- geom/ogr/doc.go | 4 - geom/ogr/ogr.go | 149 ------------------- imposm3.go | 4 +- 9 files changed, 371 insertions(+), 190 deletions(-) create mode 100644 geom/geojson/doc.go create mode 100644 geom/geojson/geojson.go delete mode 100644 geom/ogr/doc.go delete mode 100644 geom/ogr/ogr.go diff --git a/geom/geojson/doc.go b/geom/geojson/doc.go new file mode 100644 index 0000000..82cd715 --- /dev/null +++ b/geom/geojson/doc.go @@ -0,0 +1,4 @@ +/* +Package geojson creates GEOS geometries from GeoJSON files. +*/ +package geojson diff --git a/geom/geojson/geojson.go b/geom/geojson/geojson.go new file mode 100644 index 0000000..1a5df64 --- /dev/null +++ b/geom/geojson/geojson.go @@ -0,0 +1,309 @@ +package geojson + +import ( + "encoding/json" + "errors" + "fmt" + "imposm3/geom/geos" + "io" + "log" + "os" +) + +type object struct { + Type string `json:"type"` + Features []object `json:"features"` + Geometry *object `json:"geometry"` + Coordinates []interface{} `json:"coordinates"` + Properties map[string]interface{} `json:"properties"` +} + +type geometry struct { + Type string `json:"type"` + Coordinates []interface{} `json:"coordinates"` +} + +type Point struct { + Long float64 + Lat float64 +} + +func newPointFromCoords(coords []interface{}) (Point, error) { + p := Point{} + if len(coords) != 2 { + return p, errors.New("point list length not 2") + } + var ok bool + p.Long, ok = coords[0].(float64) + if !ok { + return p, errors.New("invalid lon") + } + p.Lat, ok = coords[1].(float64) + if !ok { + return p, errors.New("invalid lat") + } + return p, nil +} + +type LineString struct { + Points []Point +} + +func newLineStringFromCoords(coords []interface{}) (LineString, error) { + ls := LineString{} + + for _, part := range coords { + coord, ok := part.([]interface{}) + if !ok { + return ls, errors.New("point not a list") + } + p, err := newPointFromCoords(coord) + if err != nil { + return ls, err + } + ls.Points = append(ls.Points, p) + } + return ls, nil +} + +type Polygon struct { + LineStrings []LineString +} + +func newPolygonFromCoords(coords []interface{}) (Polygon, error) { + poly := Polygon{} + + for _, part := range coords { + lsCoords, ok := part.([]interface{}) + if !ok { + return poly, errors.New("polygon linestring not a list") + } + ls, err := newLineStringFromCoords(lsCoords) + if err != nil { + return poly, err + } + poly.LineStrings = append(poly.LineStrings, ls) + } + return poly, nil +} + +func newMultiPolygonFromCoords(coords []interface{}) ([]Polygon, error) { + mp := []Polygon{} + + for _, part := range coords { + polyCoords, ok := part.([]interface{}) + if !ok { + return mp, errors.New("multipolygon polygon not a list") + } + poly, err := newPolygonFromCoords(polyCoords) + if err != nil { + return mp, err + } + mp = append(mp, poly) + } + return mp, nil +} + +func ParseGeoJson(geojson string) (interface{}, error) { + obj := &object{} + err := json.Unmarshal([]byte(geojson), obj) + + if err != nil { + return nil, err + } + + return constructPolygons(obj) +} + +func ParseGeoJsonReader(r io.Reader) (interface{}, error) { + decoder := json.NewDecoder(r) + + obj := &object{} + err := decoder.Decode(obj) + if err != nil { + return nil, err + } + + return constructPolygons(obj) +} + +func newFeatureFromObj(obj *object) (interface{}, error) { + switch obj.Geometry.Type { + case "Point": + p, err := newPointFromCoords(obj.Geometry.Coordinates) + return p, err + case "LineString": + ls, err := newLineStringFromCoords(obj.Geometry.Coordinates) + return ls, err + case "Polygon": + poly, err := newPolygonFromCoords(obj.Geometry.Coordinates) + return poly, err + case "MultiPolygon": + poly, err := newMultiPolygonFromCoords(obj.Geometry.Coordinates) + return poly, err + default: + return nil, errors.New("unknown geometry type: " + obj.Geometry.Type) + } + return nil, nil +} + +func constructPolygons(obj *object) ([]Polygon, error) { + switch obj.Type { + case "Point": + return nil, errors.New("only Polygon or MultiPolygon are supported") + case "LineString": + return nil, errors.New("only Polygon or MultiPolygon are supported") + case "Polygon": + poly, err := newPolygonFromCoords(obj.Coordinates) + return []Polygon{poly}, err + case "MultiPolygon": + poly, err := newMultiPolygonFromCoords(obj.Coordinates) + return poly, err + case "Feature": + geom, err := constructPolygons(obj.Geometry) + return geom, err + case "FeatureCollection": + features := make([]Polygon, 0) + + for _, obj := range obj.Features { + geom, err := constructPolygons(&obj) + if err != nil { + return nil, err + } + features = append(features, geom...) + } + return features, nil + default: + return nil, errors.New("unknown type: " + obj.Type) + } + return nil, nil +} + +type GeoJson struct { + object object +} + +func NewGeoJson(r io.Reader) (*GeoJson, error) { + result := &GeoJson{} + + decoder := json.NewDecoder(r) + + err := decoder.Decode(&result.object) + if err != nil { + return nil, err + } + + return result, nil +} + +func (gj *GeoJson) Geoms() ([]*geos.Geom, error) { + + polygons, err := constructPolygons(&gj.object) + if err != nil { + return nil, err + } + + g := geos.NewGeos() + defer g.Finish() + result := []*geos.Geom{} + + for _, p := range polygons { + geom, err := geosPolygon(g, p) + if err != nil { + return nil, err + } + result = append(result, geom) + + } + return result, err +} + +func geosRing(g *geos.Geos, ls LineString) (*geos.Geom, error) { + coordSeq, err := g.CreateCoordSeq(uint32(len(ls.Points)), 2) + if err != nil { + return nil, err + } + + // coordSeq inherited by LinearRing, no destroy + for i, p := range ls.Points { + err := coordSeq.SetXY(g, uint32(i), p.Long, p.Lat) + if err != nil { + return nil, err + } + } + ring, err := coordSeq.AsLinearRing(g) + if err != nil { + // coordSeq gets Destroy by GEOS + return nil, err + } + + return ring, nil +} + +func geosPolygon(g *geos.Geos, polygon Polygon) (*geos.Geom, error) { + shell, err := geosRing(g, polygon.LineStrings[0]) + if err != nil { + return nil, err + } + + holes := make([]*geos.Geom, len(polygon.LineStrings)-1) + + for i, ls := range polygon.LineStrings[1:] { + hole, err := geosRing(g, ls) + if err != nil { + return nil, err + } + holes[i] = hole + } + + geom := g.Polygon(shell, holes) + if geom == nil { + g.Destroy(shell) + for _, hole := range holes { + g.Destroy(hole) + } + return nil, errors.New("unable to create polygon") + } + g.DestroyLater(geom) + return geom, nil +} + +var tests = []string{ + // `{"type": "Point", "coordinates": [102.0, 0.5]}`, + // `{"type": "LineString", "coordinates": [[102.1, 0.0], [103.0, 1.0], [104.0, 0.0], [105.0, 1.0]]}`, + `{"type": "Polygon", "coordinates": [[[102.1, 0.0], [103.0, 1.0], [104.0, 0.0], [105.0, 1.0]], [[0, 0]]]}`, + `{"type": "MultiPolygon", "coordinates": [[[[102.1, 0.0], [103.0, 1.0], [104.0, 0.0], [105.0, 1.0]], [[0, 0]]]]}`, +} + +func main() { + log.SetFlags(log.LstdFlags | log.Lshortfile) + for _, test := range tests { + fmt.Println("parsing: ", test) + geo, err := ParseGeoJson(test) + if err != nil { + log.Fatal(err) + } + fmt.Println(geo) + } + + fmt.Println("parsing file") + f, err := os.Open("/Users/olt/dev/cust/server/mapnik/conf/imposm/polygons/germany_clip_boundary_3857.geojson") + // f, err := os.Open("/Users/olt/dev/cust/server/mapnik/conf/imposm/polygons/germany_buffer.geojson") + if err != nil { + log.Fatal(err) + } + + gj, err := NewGeoJson(f) + if err != nil { + log.Fatal(err) + } + polygons, err := gj.Geoms() + if err != nil { + log.Fatal(err) + } + for _, p := range polygons { + fmt.Println(p.Area()) + } + // fmt.Println(geo) + +} diff --git a/geom/geom.go b/geom/geom.go index 9e9cb29..def273c 100644 --- a/geom/geom.go +++ b/geom/geom.go @@ -121,7 +121,7 @@ func Polygon(g *geos.Geos, nodes []element.Node) (*geos.Geom, error) { } // ring inherited by Polygon, no destroy - geom := g.CreatePolygon(ring, nil) + geom := g.Polygon(ring, nil) if geom == nil { g.Destroy(ring) return nil, errors.New("unable to create polygon") diff --git a/geom/geos/geos.go b/geom/geos/geos.go index 8a7d9df..fd9139a 100644 --- a/geom/geos/geos.go +++ b/geom/geos/geos.go @@ -104,17 +104,6 @@ func (this *Geos) SetHandleSrid(srid int) { this.srid = srid } -func (this *Geos) CreatePolygon(shell *Geom, holes []*Geom) *Geom { - if len(holes) > 0 { - panic("holes not implemented") - } - polygon := C.GEOSGeom_createPolygon_r(this.v, shell.v, nil, 0) - if polygon == nil { - return nil - } - return &Geom{polygon} -} - func (this *Geos) NumGeoms(geom *Geom) int32 { count := int32(C.GEOSGetNumGeometries_r(this.v, geom.v)) return count @@ -170,7 +159,7 @@ func (this *Geos) BoundsPolygon(bounds Bounds) *Geom { } // geom inherited by Polygon, no destroy - geom = this.CreatePolygon(geom, nil) + geom = this.Polygon(geom, nil) this.DestroyLater(geom) return geom diff --git a/geom/limit/limit.go b/geom/limit/limit.go index daf7893..1514ec4 100644 --- a/geom/limit/limit.go +++ b/geom/limit/limit.go @@ -2,10 +2,11 @@ package limit import ( "errors" + "imposm3/geom/geojson" "imposm3/geom/geos" - "imposm3/geom/ogr" "imposm3/logging" "math" + "os" "strings" "sync" ) @@ -58,7 +59,12 @@ func splitParams(bounds geos.Bounds, maxGrids int, minGridWidth float64) (float6 } gridWidth := math.Max(gridWidthX, gridWidthY) - return gridWidth, gridWidth * 100 + currentWidth := gridWidth + + for currentWidth <= width/2 { + currentWidth *= 2 + } + return gridWidth, currentWidth } func SplitPolygonAtAutoGrid(g *geos.Geos, geom *geos.Geom) ([]*geos.Geom, error) { @@ -89,7 +95,7 @@ func SplitPolygonAtGrid(g *geos.Geos, geom *geos.Geom, gridWidth, currentGridWid if gridWidth >= currentGridWidth { result = append(result, part) } else { - moreParts, err := SplitPolygonAtGrid(g, part, gridWidth, currentGridWidth/10.0) + moreParts, err := SplitPolygonAtGrid(g, part, gridWidth, currentGridWidth/2.0) if err != nil { return nil, err } @@ -107,12 +113,22 @@ type Limiter struct { bufferedBbox geos.Bounds } -func NewFromOgrSource(source string) (*Limiter, error) { - return NewFromOgrSourceWithBuffered(source, 0.0) +func NewFromGeoJson(source string) (*Limiter, error) { + return NewFromGeoJsonWithBuffered(source, 0.0) } -func NewFromOgrSourceWithBuffered(source string, buffer float64) (*Limiter, error) { - ds, err := ogr.Open(source) +func NewFromGeoJsonWithBuffered(source string, buffer float64) (*Limiter, error) { + + f, err := os.Open(source) + if err != nil { + return nil, err + } + + gj, err := geojson.NewGeoJson(f) + if err != nil { + return nil, err + } + geoms, err := gj.Geoms() if err != nil { return nil, err } @@ -120,11 +136,6 @@ func NewFromOgrSourceWithBuffered(source string, buffer float64) (*Limiter, erro g := geos.NewGeos() defer g.Finish() - layer, err := ds.Layer() - if err != nil { - return nil, err - } - index := g.CreateIndex() var polygons []*geos.Geom @@ -134,7 +145,7 @@ func NewFromOgrSourceWithBuffered(source string, buffer float64) (*Limiter, erro withBuffer = true } - for geom := range layer.Geoms() { + for _, geom := range geoms { if withBuffer { simplified := g.SimplifyPreserveTopology(geom, 1000) if simplified == nil { diff --git a/geom/limit/limit_test.go b/geom/limit/limit_test.go index 01e41aa..2d84522 100644 --- a/geom/limit/limit_test.go +++ b/geom/limit/limit_test.go @@ -25,9 +25,9 @@ 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, 0.1, 0.05, 0.11}, {0.05, 0.1, 0.1, 0.11}, {0.1, 0, 0.15, 0.05}, {0.1, 0.05, 0.15, 0.1}, @@ -39,7 +39,10 @@ func TestSplitPolygonAtGrids(t *testing.T) { geom := g.BoundsPolygon(geos.Bounds{0, 0, 0.15, 0.11}) - geoms, _ := SplitPolygonAtGrid(g, geom, 0.05, 5.0) + geoms, _ := SplitPolygonAtGrid(g, geom, 0.05, 0.2) + for _, geom := range geoms { + t.Log(geom.Bounds()) + } for i, geom := range geoms { if expected[i] != geom.Bounds() { t.Fatalf("%v != %v\n", expected[i], geom.Bounds()) @@ -315,37 +318,55 @@ func TestClipperWithBuffer(t *testing.T) { } func TestSplitParams(t *testing.T) { - var gridWidth float64 + var gridWidth, startWidth float64 - gridWidth, _ = splitParams(geos.Bounds{0, 0, 10000, 10000}, 10, 2000) + gridWidth, startWidth = splitParams(geos.Bounds{0, 0, 10000, 10000}, 10, 2000) if gridWidth != 2000.0 { t.Fatal(gridWidth) } + if startWidth != 8000.0 { + t.Fatal(startWidth) + } - gridWidth, _ = splitParams(geos.Bounds{0, 0, 10000, 10000}, 10, 1000) + gridWidth, startWidth = splitParams(geos.Bounds{0, 0, 10000, 10000}, 10, 1000) if gridWidth != 1000.0 { t.Fatal(gridWidth) } + if startWidth != 8000.0 { + t.Fatal(startWidth) + } - gridWidth, _ = splitParams(geos.Bounds{0, 0, 10000, 10000}, 10, 500) + gridWidth, startWidth = splitParams(geos.Bounds{0, 0, 10000, 10000}, 10, 500) if gridWidth != 1000.0 { t.Fatal(gridWidth) } + if startWidth != 8000.0 { + t.Fatal(startWidth) + } - gridWidth, _ = splitParams(geos.Bounds{0, 0, 10000, 5000}, 10, 500) + gridWidth, startWidth = splitParams(geos.Bounds{0, 0, 10000, 5000}, 10, 500) if gridWidth != 1000.0 { t.Fatal(gridWidth) } + if startWidth != 8000.0 { + t.Fatal(startWidth) + } - gridWidth, _ = splitParams(geos.Bounds{0, 0, 10000, 20000}, 10, 500) + gridWidth, startWidth = splitParams(geos.Bounds{0, 0, 10000, 20000}, 10, 500) if gridWidth != 2000.0 { t.Fatal(gridWidth) } + if startWidth != 8000.0 { + t.Fatal(startWidth) + } - gridWidth, _ = splitParams(geos.Bounds{0, 0, 10000, 20000}, 50, 100) + gridWidth, startWidth = splitParams(geos.Bounds{0, 0, 10000, 20000}, 50, 100) if gridWidth != 400.0 { t.Fatal(gridWidth) } + if startWidth != 6400.0 { + t.Fatal(startWidth) + } } diff --git a/geom/ogr/doc.go b/geom/ogr/doc.go deleted file mode 100644 index f68871b..0000000 --- a/geom/ogr/doc.go +++ /dev/null @@ -1,4 +0,0 @@ -/* -Package ogr provides a wrapper to the GDAL/OGR library. -*/ -package ogr diff --git a/geom/ogr/ogr.go b/geom/ogr/ogr.go deleted file mode 100644 index 4711070..0000000 --- a/geom/ogr/ogr.go +++ /dev/null @@ -1,149 +0,0 @@ -package ogr - -/* -#cgo LDFLAGS: -lgdal -#include "ogr_api.h" -#include "cpl_error.h" -#include "cpl_conv.h" -*/ -import "C" -import ( - "fmt" - "imposm3/geom/geos" - "strings" - "unsafe" -) - -func init() { - C.OGRRegisterAll() -} - -type DataSource struct { - v C.OGRDataSourceH -} - -type Layer struct { - v C.OGRLayerH -} -type OgrError struct { - message string -} - -func (e *OgrError) Error() string { - return e.message -} - -func lastOgrError(fallback string) error { - msg := C.CPLGetLastErrorMsg() - if msg == nil { - return &OgrError{fallback} - } - str := C.GoString(msg) - if str == "" { - return &OgrError{fallback} - } - return &OgrError{str} -} - -func Open(name string) (*DataSource, error) { - namec := C.CString(name) - defer C.free(unsafe.Pointer(namec)) - ds := C.OGROpen(namec, 0, nil) - if ds == nil { - return nil, lastOgrError("failed to open") - } - return &DataSource{ds}, nil -} - -func (ds *DataSource) Layer() (*Layer, error) { - layer := C.OGR_DS_GetLayer(ds.v, 0) - if layer == nil { - return nil, lastOgrError("failed to get layer 0") - } - return &Layer{layer}, nil -} - -func (ds *DataSource) Query(query string) (*Layer, error) { - // create select query if it is only a where statement - if !strings.HasPrefix(strings.ToLower(query), "select") { - layer, err := ds.Layer() - if err != nil { - return nil, err - } - layerDef := C.OGR_L_GetLayerDefn(layer.v) - name := C.OGR_FD_GetName(layerDef) - query = fmt.Sprintf("SELECT * FROM %s WHERE %s", C.GoString(name), query) - } - queryc := C.CString(query) - defer C.free(unsafe.Pointer(queryc)) - layer := C.OGR_DS_ExecuteSQL(ds.v, queryc, nil, nil) - if layer == nil { - return nil, lastOgrError("unable to execute query '" + query + "'") - } - return &Layer{layer}, nil -} - -func (layer *Layer) Wkts() chan string { - wkts := make(chan string) - - go func() { - defer close(wkts) - - C.OGR_L_ResetReading(layer.v) - for { - feature := C.OGR_L_GetNextFeature(layer.v) - if feature == nil { - break - } - geom := C.OGR_F_GetGeometryRef(feature) - var res *C.char - C.OGR_G_ExportToWkt(geom, &res) - wkts <- C.GoString(res) - C.CPLFree(unsafe.Pointer(res)) - C.OGR_F_Destroy(feature) - - } - }() - - return wkts -} - -func (layer *Layer) Wkbs() chan []byte { - wkbs := make(chan []byte) - - go func() { - defer close(wkbs) - - C.OGR_L_ResetReading(layer.v) - for { - feature := C.OGR_L_GetNextFeature(layer.v) - if feature == nil { - break - } - geom := C.OGR_F_GetGeometryRef(feature) - size := C.OGR_G_WkbSize(geom) - buf := make([]byte, size) - C.OGR_G_ExportToWkb(geom, C.wkbNDR, (*C.uchar)(&buf[0])) - wkbs <- buf - C.OGR_F_Destroy(feature) - } - }() - return wkbs -} - -func (layer *Layer) Geoms() chan *geos.Geom { - geoms := make(chan *geos.Geom) - - go func() { - defer close(geoms) - g := geos.NewGeos() - defer g.Finish() - - wkbs := layer.Wkbs() - for wkb := range wkbs { - geom := g.FromWkb(wkb) - geoms <- geom - } - }() - return geoms -} diff --git a/imposm3.go b/imposm3.go index 46a3f2b..a260deb 100644 --- a/imposm3.go +++ b/imposm3.go @@ -61,7 +61,7 @@ func main() { if config.BaseOptions.LimitTo != "" { var err error step := log.StartStep("Reading limitto geometries") - geometryLimiter, err = limit.NewFromOgrSourceWithBuffered( + geometryLimiter, err = limit.NewFromGeoJsonWithBuffered( config.BaseOptions.LimitTo, config.BaseOptions.LimitToCacheBuffer, ) @@ -135,7 +135,7 @@ func mainimport() { if config.ImportOptions.Write && config.BaseOptions.LimitTo != "" { var err error step := log.StartStep("Reading limitto geometries") - geometryLimiter, err = limit.NewFromOgrSource(config.BaseOptions.LimitTo) + geometryLimiter, err = limit.NewFromGeoJson(config.BaseOptions.LimitTo) if err != nil { log.Fatal(err) }