replace OGR dependency for limitto with simple GeoJSON parser

master
Oliver Tonnhofer 2013-10-21 14:59:31 +02:00
parent f0d5972908
commit ac094e8dce
9 changed files with 371 additions and 190 deletions

4
geom/geojson/doc.go Normal file
View File

@ -0,0 +1,4 @@
/*
Package geojson creates GEOS geometries from GeoJSON files.
*/
package geojson

309
geom/geojson/geojson.go Normal file
View File

@ -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)
}

View File

@ -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")

View File

@ -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

View File

@ -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 {

View File

@ -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)
}
}

View File

@ -1,4 +0,0 @@
/*
Package ogr provides a wrapper to the GDAL/OGR library.
*/
package ogr

View File

@ -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
}

View File

@ -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)
}