support -limitto with -srid 4326

- refactored geojson/limit package
- geojson files are now required to be in EPSG:4326
Oliver Tonnhofer 2015-03-02 09:22:08 +01:00
parent 60ff020b57
commit 112d889639
10 changed files with 206 additions and 195 deletions

View File

@ -62,9 +62,10 @@ func Main(usage func()) {
if config.BaseOptions.LimitTo != "" { if config.BaseOptions.LimitTo != "" {
var err error var err error
step := log.StartStep("Reading limitto geometries") step := log.StartStep("Reading limitto geometries")
geometryLimiter, err = limit.NewFromGeoJsonWithBuffered( geometryLimiter, err = limit.NewFromGeoJSON(
config.BaseOptions.LimitTo, config.BaseOptions.LimitTo,
config.BaseOptions.LimitToCacheBuffer, config.BaseOptions.LimitToCacheBuffer,
config.BaseOptions.Srid,
) )
if err != nil { if err != nil {
log.Fatal(err) log.Fatal(err)

View File

@ -4,11 +4,13 @@ import (
"encoding/json" "encoding/json"
"errors" "errors"
"fmt" "fmt"
"github.com/omniscale/imposm3/geom/geos"
"github.com/omniscale/imposm3/proj"
"io" "io"
"github.com/omniscale/imposm3/logging"
) )
var log = logging.NewLogger("geojson")
type object struct { type object struct {
Type string `json:"type"` Type string `json:"type"`
Features []object `json:"features"` Features []object `json:"features"`
@ -22,36 +24,37 @@ type geometry struct {
Coordinates []interface{} `json:"coordinates"` Coordinates []interface{} `json:"coordinates"`
} }
type point struct { type Point struct {
long float64 Long float64
lat float64 Lat float64
} }
func newPointFromCoords(coords []interface{}) (point, error) { func newPointFromCoords(coords []interface{}) (Point, error) {
p := point{} p := Point{}
if len(coords) != 2 && len(coords) != 3 { if len(coords) != 2 && len(coords) != 3 {
return p, errors.New("point list length not 2 or 3") return p, errors.New("point list length not 2 or 3")
} }
var ok bool var ok bool
p.long, ok = coords[0].(float64) p.Long, ok = coords[0].(float64)
if !ok { if !ok {
return p, errors.New("invalid lon") return p, errors.New("invalid lon")
} }
p.lat, ok = coords[1].(float64) p.Lat, ok = coords[1].(float64)
if !ok { if !ok {
return p, errors.New("invalid lat") return p, errors.New("invalid lat")
} }
if p.long >= -180.0 && p.long <= 180.0 && p.lat >= -90.0 && p.lat <= 90.0 { if p.Long >= 180.0 || p.Long <= -180.0 || p.Lat >= 90.0 || p.Lat <= -90.0 {
p.long, p.lat = proj.WgsToMerc(p.long, p.lat) log.Warn("coordinates outside of world boundary. non-4326?")
} }
return p, nil return p, nil
} }
type lineString []point type LineString []Point
func newLineStringFromCoords(coords []interface{}) (lineString, error) { func newLineStringFromCoords(coords []interface{}) (LineString, error) {
ls := lineString{} ls := LineString{}
for _, part := range coords { for _, part := range coords {
coord, ok := part.([]interface{}) coord, ok := part.([]interface{})
@ -67,15 +70,15 @@ func newLineStringFromCoords(coords []interface{}) (lineString, error) {
return ls, nil return ls, nil
} }
type polygon []lineString type Polygon []LineString
type polygonFeature struct { type polygonFeature struct {
polygon polygon polygon Polygon
properties map[string]string properties map[string]string
} }
type Feature struct { type Feature struct {
Geom *geos.Geom Polygon Polygon
Properties map[string]string Properties map[string]string
} }
@ -87,13 +90,13 @@ func stringProperties(properties map[string]interface{}) map[string]string {
return result return result
} }
func newPolygonFromCoords(coords []interface{}) (polygon, error) { func newPolygonFromCoords(coords []interface{}) (Polygon, error) {
poly := polygon{} poly := Polygon{}
for _, part := range coords { for _, part := range coords {
lsCoords, ok := part.([]interface{}) lsCoords, ok := part.([]interface{})
if !ok { if !ok {
return poly, errors.New("polygon lineString not a list") return poly, errors.New("polygon LineString not a list")
} }
ls, err := newLineStringFromCoords(lsCoords) ls, err := newLineStringFromCoords(lsCoords)
if err != nil { if err != nil {
@ -104,8 +107,8 @@ func newPolygonFromCoords(coords []interface{}) (polygon, error) {
return poly, nil return poly, nil
} }
func newMultiPolygonFeaturesFromCoords(coords []interface{}) ([]polygonFeature, error) { func newMultiPolygonFeaturesFromCoords(coords []interface{}) ([]Feature, error) {
features := []polygonFeature{} features := []Feature{}
for _, part := range coords { for _, part := range coords {
polyCoords, ok := part.([]interface{}) polyCoords, ok := part.([]interface{})
@ -116,12 +119,14 @@ func newMultiPolygonFeaturesFromCoords(coords []interface{}) ([]polygonFeature,
if err != nil { if err != nil {
return features, err return features, err
} }
features = append(features, polygonFeature{poly, nil}) features = append(features, Feature{poly, nil})
} }
return features, nil return features, nil
} }
func ParseGeoJson(r io.Reader) ([]Feature, error) { // ParseGeoJSON parses geojson from reader and returns []Feature in WGS84 and
// another []Feature tranformed in targetSRID.
func ParseGeoJSON(r io.Reader) ([]Feature, error) {
decoder := json.NewDecoder(r) decoder := json.NewDecoder(r)
obj := &object{} obj := &object{}
@ -136,22 +141,10 @@ func ParseGeoJson(r io.Reader) ([]Feature, error) {
return nil, err return nil, err
} }
g := geos.NewGeos() return polygons, nil
defer g.Finish()
result := []Feature{}
for _, p := range polygons {
geom, err := geosPolygon(g, p.polygon)
if err != nil {
return nil, err
}
result = append(result, Feature{geom, p.properties})
}
return result, err
} }
func constructPolygonFeatures(obj *object) ([]polygonFeature, error) { func constructPolygonFeatures(obj *object) ([]Feature, error) {
switch obj.Type { switch obj.Type {
case "Point": case "Point":
return nil, errors.New("only polygon or MultiPolygon are supported") return nil, errors.New("only polygon or MultiPolygon are supported")
@ -159,7 +152,7 @@ func constructPolygonFeatures(obj *object) ([]polygonFeature, error) {
return nil, errors.New("only polygon or MultiPolygon are supported") return nil, errors.New("only polygon or MultiPolygon are supported")
case "Polygon": case "Polygon":
poly, err := newPolygonFromCoords(obj.Coordinates) poly, err := newPolygonFromCoords(obj.Coordinates)
return []polygonFeature{{poly, nil}}, err return []Feature{{poly, nil}}, err
case "MultiPolygon": case "MultiPolygon":
poly, err := newMultiPolygonFeaturesFromCoords(obj.Coordinates) poly, err := newMultiPolygonFeaturesFromCoords(obj.Coordinates)
return poly, err return poly, err
@ -170,11 +163,11 @@ func constructPolygonFeatures(obj *object) ([]polygonFeature, error) {
} }
properties := stringProperties(obj.Properties) properties := stringProperties(obj.Properties)
for i, _ := range features { for i, _ := range features {
features[i].properties = properties features[i].Properties = properties
} }
return features, err return features, err
case "FeatureCollection": case "FeatureCollection":
features := make([]polygonFeature, 0) features := make([]Feature, 0)
for _, obj := range obj.Features { for _, obj := range obj.Features {
f, err := constructPolygonFeatures(&obj) f, err := constructPolygonFeatures(&obj)
@ -188,56 +181,3 @@ func constructPolygonFeatures(obj *object) ([]polygonFeature, error) {
return nil, errors.New("unknown type: " + obj.Type) return nil, errors.New("unknown type: " + obj.Type)
} }
} }
func geosRing(g *geos.Geos, ls lineString) (*geos.Geom, error) {
coordSeq, err := g.CreateCoordSeq(uint32(len(ls)), 2)
if err != nil {
return nil, err
}
// coordSeq inherited by LinearRing, no destroy
for i, p := range ls {
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) {
if len(polygon) == 0 {
return nil, errors.New("empty polygon")
}
shell, err := geosRing(g, polygon[0])
if err != nil {
return nil, err
}
holes := make([]*geos.Geom, len(polygon)-1)
for i, ls := range polygon[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")
}
return geom, nil
}

View File

@ -2,13 +2,12 @@ package geojson
import ( import (
"bytes" "bytes"
"math"
"testing" "testing"
) )
func TestParsePolygon(t *testing.T) { func TestParsePolygon(t *testing.T) {
r := bytes.NewBufferString(`{"type": "Polygon", "coordinates": [[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 2000], [1000, 1000]]]}`) r := bytes.NewBufferString(`{"type": "Polygon", "coordinates": [[[8, 50], [11, 50], [11, 53], [8, 53], [8, 50]]]}`)
features, err := ParseGeoJson(r) features, err := ParseGeoJSON(r)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
@ -18,13 +17,13 @@ func TestParsePolygon(t *testing.T) {
t.Fatal(features) t.Fatal(features)
} }
if math.Abs(features[0].Geom.Area()-1000000) > 0.00001 { if len(features[0].Polygon[0]) != 5 {
t.Fatal(features[0].Geom.Area()) t.Fatal(features)
} }
// ignore z values // ignore z values
r = bytes.NewBufferString(`{"type": "Polygon", "coordinates": [[[1000, 1000, 1000], [2000, 1000, 1000], [2000, 2000, 1000], [1000, 2000, 1000], [1000, 1000, 1000]]]}`) r = bytes.NewBufferString(`{"type": "Polygon", "coordinates": [[[8, 50, 0], [11, 50, 0], [11, 53, 0], [8, 53, 0], [8, 50, 0]]]}`)
features, err = ParseGeoJson(r) features, err = ParseGeoJSON(r)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
@ -34,12 +33,17 @@ func TestParsePolygon(t *testing.T) {
t.Fatal(features) t.Fatal(features)
} }
if math.Abs(features[0].Geom.Area()-1000000) > 0.00001 { if len(features[0].Polygon[0]) != 5 {
t.Fatal(features[0].Geom.Area()) t.Fatal(features)
} }
r = bytes.NewBufferString(`{"type": "Polygon", "coordinates": [[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 2000], [1000, 1000]], [[500, 500], [600, 500], [600, 600], [500, 600], [500, 500]]]}`) if p := features[0].Polygon[0][0]; p.Long != 8.0 || p.Lat != 50 {
features, err = ParseGeoJson(r) t.Fatal(features)
}
// with hole
r = bytes.NewBufferString(`{"type": "Polygon", "coordinates": [[[8, 50], [11, 50], [11, 53], [8, 53], [8, 50]], [[9, 51], [10, 51], [10, 52], [9, 52], [9, 51]]]}`)
features, err = ParseGeoJSON(r)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
@ -49,18 +53,18 @@ func TestParsePolygon(t *testing.T) {
t.Fatal(features) t.Fatal(features)
} }
if math.Abs(features[0].Geom.Area()-990000) > 0.00001 { if len(features[0].Polygon) != 2 {
t.Fatal(features[0].Geom.Area()) t.Fatal(features)
} }
} }
func TestParseMultiPolygon(t *testing.T) { func TestParseMultiPolygon(t *testing.T) {
r := bytes.NewBufferString(`{"type": "MultiPolygon", "coordinates": r := bytes.NewBufferString(`{"type": "MultiPolygon", "coordinates":
[[[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 1000]]], [[[[8, 50], [11, 50], [11, 53], [8, 50]]],
[[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 1000]]]] [[[8, 50], [11, 50], [11, 53], [8, 50]]]]
}`) }`)
features, err := ParseGeoJson(r) features, err := ParseGeoJSON(r)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
@ -73,9 +77,9 @@ func TestParseMultiPolygon(t *testing.T) {
func TestParseFeature(t *testing.T) { func TestParseFeature(t *testing.T) {
r := bytes.NewBufferString(`{"type": "Feature", "geometry": { r := bytes.NewBufferString(`{"type": "Feature", "geometry": {
"type": "Polygon", "coordinates": [[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 2000], [1000, 1000]]] "type": "Polygon", "coordinates": [[[8, 50], [11, 50], [11, 53], [8, 53], [8, 50]]]
}}`) }}`)
features, err := ParseGeoJson(r) features, err := ParseGeoJSON(r)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
@ -84,21 +88,21 @@ func TestParseFeature(t *testing.T) {
if len(features) != 1 { if len(features) != 1 {
t.Fatal(features) t.Fatal(features)
} }
if math.Abs(features[0].Geom.Area()-1000000) > 0.00001 { if len(features[0].Polygon[0]) != 5 {
t.Fatal(features[0].Geom.Area()) t.Fatal(features)
} }
} }
func TestParseFeatureCollection(t *testing.T) { func TestParseFeatureCollection(t *testing.T) {
r := bytes.NewBufferString(`{"type": "FeatureCollection", "features": [ r := bytes.NewBufferString(`{"type": "FeatureCollection", "features": [
{"type": "Feature", "geometry": {"type": "Feature", "geometry":
{"type": "Polygon", "coordinates": [[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 2000], [1000, 1000]]]} {"type": "Polygon", "coordinates": [[[8, 50], [11, 50], [11, 53], [8, 53], [8, 50]]]}
}, },
{"type": "Feature", "geometry": {"type": "Feature", "geometry":
{"type": "Polygon", "coordinates": [[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 2000], [1000, 1000]]]} {"type": "Polygon", "coordinates": [[[8, 50], [11, 50], [11, 53], [8, 53], [8, 50]]]}
} }
]}`) ]}`)
features, err := ParseGeoJson(r) features, err := ParseGeoJSON(r)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
@ -107,21 +111,24 @@ func TestParseFeatureCollection(t *testing.T) {
if len(features) != 2 { if len(features) != 2 {
t.Fatal(features) t.Fatal(features)
} }
if math.Abs(features[0].Geom.Area()-1000000) > 0.00001 { if len(features[0].Polygon[0]) != 5 {
t.Fatal(features[0].Geom.Area()) t.Fatal(features)
}
if len(features[1].Polygon[0]) != 5 {
t.Fatal(features)
} }
} }
func TestParseGeoJson(t *testing.T) { func TestParseProperties(t *testing.T) {
r := bytes.NewBufferString(`{"type": "FeatureCollection", "features": [ r := bytes.NewBufferString(`{"type": "FeatureCollection", "features": [
{"type": "Feature", "properties": {"foo": "bar", "baz": 42}, "geometry": {"type": "Feature", "properties": {"foo": "bar", "baz": 42}, "geometry":
{"type": "Polygon", "coordinates": [[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 2000], [1000, 1000]]]} {"type": "Polygon", "coordinates": [[[8, 50], [11, 50], [11, 53], [8, 53], [8, 50]]]}
}, },
{"type": "Feature", "geometry": {"type": "Feature", "geometry":
{"type": "Polygon", "coordinates": [[[1000, 1000], [2000, 1000], [2000, 2000], [1000, 2000], [1000, 1000]]]} {"type": "Polygon", "coordinates": [[[8, 50], [11, 50], [11, 53], [8, 53], [8, 50]]]}
} }
]}`) ]}`)
features, err := ParseGeoJson(r) features, err := ParseGeoJSON(r)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
@ -137,31 +144,10 @@ func TestParseGeoJson(t *testing.T) {
t.Errorf("baz != 42, but '%v'", v) t.Errorf("baz != 42, but '%v'", v)
} }
if math.Abs(features[0].Geom.Area()-1000000) > 0.00001 { if len(features[0].Polygon[0]) != 5 {
t.Fatal(features[0].Geom.Area())
}
}
func TestParseGeoJsonTransform(t *testing.T) {
// automatically transforms WGS84 to webmercator
r := bytes.NewBufferString(`{"type": "FeatureCollection", "features": [
{"type": "Feature", "geometry":
{"type": "Polygon", "coordinates": [[[8, 53], [9, 53], [9, 54], [8, 54], [8, 53]]]}
},
{"type": "Feature", "geometry":
{"type": "Polygon", "coordinates": [[[9, 53], [10, 53], [10, 54], [9, 54], [9, 53]]]}
}
]}`)
features, err := ParseGeoJson(r)
if err != nil {
t.Fatal(err)
}
if len(features) != 2 {
t.Fatal(features) t.Fatal(features)
} }
if math.Abs(features[0].Geom.Area()-20834374847.98027) > 0.01 { if len(features[1].Polygon[0]) != 5 {
t.Fatal(features[0].Geom.Area()) t.Fatal(features)
} }
} }

View File

@ -2,13 +2,15 @@ package limit
import ( import (
"errors" "errors"
"github.com/omniscale/imposm3/geom/geojson"
"github.com/omniscale/imposm3/geom/geos"
"github.com/omniscale/imposm3/logging"
"math" "math"
"os" "os"
"strings" "strings"
"sync" "sync"
"github.com/omniscale/imposm3/geom/geojson"
"github.com/omniscale/imposm3/geom/geos"
"github.com/omniscale/imposm3/logging"
"github.com/omniscale/imposm3/proj"
) )
var log = logging.NewLogger("limiter") var log = logging.NewLogger("limiter")
@ -67,16 +69,16 @@ func splitParams(bounds geos.Bounds, maxGrids int, minGridWidth float64) (float6
return gridWidth, currentWidth return gridWidth, currentWidth
} }
func SplitPolygonAtAutoGrid(g *geos.Geos, geom *geos.Geom) ([]*geos.Geom, error) { func splitPolygonAtAutoGrid(g *geos.Geos, geom *geos.Geom, minGridWidth float64) ([]*geos.Geom, error) {
geomBounds := geom.Bounds() geomBounds := geom.Bounds()
if geomBounds == geos.NilBounds { if geomBounds == geos.NilBounds {
return nil, errors.New("couldn't create bounds for geom") return nil, errors.New("couldn't create bounds for geom")
} }
gridWidth, currentGridWidth := splitParams(geomBounds, 32, 50000.0) gridWidth, currentGridWidth := splitParams(geomBounds, 32, minGridWidth)
return SplitPolygonAtGrid(g, geom, gridWidth, currentGridWidth) return splitPolygonAtGrid(g, geom, gridWidth, currentGridWidth)
} }
func SplitPolygonAtGrid(g *geos.Geos, geom *geos.Geom, gridWidth, currentGridWidth float64) ([]*geos.Geom, error) { func splitPolygonAtGrid(g *geos.Geos, geom *geos.Geom, gridWidth, currentGridWidth float64) ([]*geos.Geom, error) {
var result []*geos.Geom var result []*geos.Geom
geomBounds := geom.Bounds() geomBounds := geom.Bounds()
if geomBounds == geos.NilBounds { if geomBounds == geos.NilBounds {
@ -96,7 +98,7 @@ func SplitPolygonAtGrid(g *geos.Geos, geom *geos.Geom, gridWidth, currentGridWid
if gridWidth >= currentGridWidth { if gridWidth >= currentGridWidth {
result = append(result, part) result = append(result, part)
} else { } else {
moreParts, err := SplitPolygonAtGrid(g, part, gridWidth, currentGridWidth/2.0) moreParts, err := splitPolygonAtGrid(g, part, gridWidth, currentGridWidth/2.0)
g.Destroy(part) g.Destroy(part)
if err != nil { if err != nil {
return nil, err return nil, err
@ -128,19 +130,14 @@ type Limiter struct {
bufferedPrepMu *sync.Mutex bufferedPrepMu *sync.Mutex
} }
func NewFromGeoJson(source string) (*Limiter, error) { func NewFromGeoJSON(source string, buffer float64, targetSRID int) (*Limiter, error) {
return NewFromGeoJsonWithBuffered(source, 0.0)
}
func NewFromGeoJsonWithBuffered(source string, buffer float64) (*Limiter, error) {
f, err := os.Open(source) f, err := os.Open(source)
if err != nil { if err != nil {
return nil, err return nil, err
} }
defer f.Close() defer f.Close()
features, err := geojson.ParseGeoJson(f) features, err := geojson.ParseGeoJSON(f)
if err != nil { if err != nil {
return nil, err return nil, err
} }
@ -158,9 +155,13 @@ func NewFromGeoJsonWithBuffered(source string, buffer float64) (*Limiter, error)
withBuffer = true withBuffer = true
} }
for _, feature := range features { if withBuffer {
if withBuffer { for _, feature := range features {
simplified := g.SimplifyPreserveTopology(feature.Geom, 1000) geom, err := geosPolygon(g, feature.Polygon)
if err != nil {
return nil, err
}
simplified := g.SimplifyPreserveTopology(geom, 0.01)
if simplified == nil { if simplified == nil {
return nil, errors.New("couldn't simplify limitto") return nil, errors.New("couldn't simplify limitto")
} }
@ -172,9 +173,24 @@ func NewFromGeoJsonWithBuffered(source string, buffer float64) (*Limiter, error)
// buffered gets destroyed in UnionPolygons // buffered gets destroyed in UnionPolygons
bufferedPolygons = append(bufferedPolygons, buffered) bufferedPolygons = append(bufferedPolygons, buffered)
} }
polygons = append(polygons, feature.Geom) }
for _, feature := range features {
if targetSRID != 4326 {
// transforms polygon in-place
transformPolygon(feature.Polygon, targetSRID)
}
geom, err := geosPolygon(g, feature.Polygon)
if err != nil {
return nil, err
}
parts, err := SplitPolygonAtAutoGrid(g, feature.Geom) polygons = append(polygons, geom)
minGridWidth := 50000.0
if targetSRID == 4326 {
minGridWidth = 0.5
}
parts, err := splitPolygonAtAutoGrid(g, geom, minGridWidth)
if err != nil { if err != nil {
return nil, err return nil, err
@ -191,7 +207,7 @@ func NewFromGeoJsonWithBuffered(source string, buffer float64) (*Limiter, error)
if union == nil { if union == nil {
return nil, errors.New("unable to union limitto polygons") return nil, errors.New("unable to union limitto polygons")
} }
simplified := g.SimplifyPreserveTopology(union, 5000) simplified := g.SimplifyPreserveTopology(union, 0.05)
if simplified == nil { if simplified == nil {
return nil, errors.New("unable to simplify limitto polygons") return nil, errors.New("unable to simplify limitto polygons")
} }
@ -254,6 +270,10 @@ func filterGeometryByType(g *geos.Geos, geom *geos.Geom, targetType string) []*g
return []*geos.Geom{} return []*geos.Geom{}
} }
// Clip returns geom (in targetSRID) clipped to the LimitTo geometry.
// Returns nil if geom is outside of the LimitTo geometry.
// Returns only similar geometry types (e.g. clipped Polygon will return
// one or more Polygons, but no LineString or Point, etc.)
func (l *Limiter) Clip(geom *geos.Geom) ([]*geos.Geom, error) { func (l *Limiter) Clip(geom *geos.Geom) ([]*geos.Geom, error) {
g := geos.NewGeos() g := geos.NewGeos()
defer g.Finish() defer g.Finish()
@ -298,6 +318,8 @@ func (l *Limiter) Clip(geom *geos.Geom) ([]*geos.Geom, error) {
return mergeGeometries(g, intersections, geomType), nil return mergeGeometries(g, intersections, geomType), nil
} }
// IntersectsBuffer returns true if the point (EPSG:4326) intersects the buffered
// LimitTo geometry.
func (c *Limiter) IntersectsBuffer(g *geos.Geos, x, y float64) bool { func (c *Limiter) IntersectsBuffer(g *geos.Geos, x, y float64) bool {
if c.bufferedPrep == nil { if c.bufferedPrep == nil {
return true return true
@ -330,7 +352,7 @@ func flattenPolygons(g *geos.Geos, geoms []*geos.Geom) []*geos.Geom {
} else if g.Type(geom) == "Polygon" { } else if g.Type(geom) == "Polygon" {
result = append(result, geom) result = append(result, geom)
} else { } else {
log.Printf("unexpected geometry type in flattenPolygons") log.Printf("unexpected geometry type in flattenPolygons: %s", g.Type(geom))
g.Destroy(geom) g.Destroy(geom)
} }
} }
@ -348,7 +370,7 @@ func flattenLineStrings(g *geos.Geos, geoms []*geos.Geom) []*geos.Geom {
} else if g.Type(geom) == "LineString" { } else if g.Type(geom) == "LineString" {
result = append(result, geom) result = append(result, geom)
} else { } else {
log.Printf("unexpected geometry type in flattenPolygons") log.Printf("unexpected geometry type in flattenLineStrings: %s", g.Type(geom))
g.Destroy(geom) g.Destroy(geom)
} }
} }
@ -402,3 +424,67 @@ func mergeGeometries(g *geos.Geos, geoms []*geos.Geom, geomType string) []*geos.
panic("unexpected geometry type" + geomType) panic("unexpected geometry type" + geomType)
} }
} }
func geosRing(g *geos.Geos, ls geojson.LineString) (*geos.Geom, error) {
coordSeq, err := g.CreateCoordSeq(uint32(len(ls)), 2)
if err != nil {
return nil, err
}
// coordSeq inherited by LinearRing, no destroy
for i, p := range ls {
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 geojson.Polygon) (*geos.Geom, error) {
if len(polygon) == 0 {
return nil, errors.New("empty polygon")
}
shell, err := geosRing(g, polygon[0])
if err != nil {
return nil, err
}
holes := make([]*geos.Geom, len(polygon)-1)
for i, ls := range polygon[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")
}
return geom, nil
}
func transformPolygon(p geojson.Polygon, targetSRID int) {
if targetSRID != 3857 {
panic("transformation to non-4326/3856 not implemented")
}
for _, ls := range p {
for i := range ls {
ls[i].Long, ls[i].Lat = proj.WgsToMerc(ls[i].Long, ls[i].Lat)
}
}
}

View File

@ -1,8 +1,9 @@
package limit package limit
import ( import (
"github.com/omniscale/imposm3/geom/geos"
"testing" "testing"
"github.com/omniscale/imposm3/geom/geos"
) )
func TestTileBounds(t *testing.T) { func TestTileBounds(t *testing.T) {
@ -39,7 +40,7 @@ func TestSplitPolygonAtGrids(t *testing.T) {
geom := g.BoundsPolygon(geos.Bounds{0, 0, 0.15, 0.11}) geom := g.BoundsPolygon(geos.Bounds{0, 0, 0.15, 0.11})
geoms, _ := SplitPolygonAtGrid(g, geom, 0.05, 0.2) geoms, _ := splitPolygonAtGrid(g, geom, 0.05, 0.2)
for _, geom := range geoms { for _, geom := range geoms {
t.Log(geom.Bounds()) t.Log(geom.Bounds())
} }
@ -263,7 +264,7 @@ func TestFilterGeometryByType(t *testing.T) {
func TestClipper(t *testing.T) { func TestClipper(t *testing.T) {
g := geos.NewGeos() g := geos.NewGeos()
defer g.Finish() defer g.Finish()
limiter, err := NewFromGeoJson("./hamburg_clip.geojson") limiter, err := NewFromGeoJSON("./clipping.geojson", 0.0, 3857)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
} }
@ -305,16 +306,17 @@ func TestClipper(t *testing.T) {
func TestClipperWithBuffer(t *testing.T) { func TestClipperWithBuffer(t *testing.T) {
g := geos.NewGeos() g := geos.NewGeos()
defer g.Finish() defer g.Finish()
limiter, err := NewFromGeoJsonWithBuffered("./hamburg_clip.geojson", 10000.0) limiter, err := NewFromGeoJSON("./clipping.geojson", 0.1, 3857)
if err != nil { if err != nil {
t.Fatal(err) t.Fatal(err)
} }
if limiter.IntersectsBuffer(g, 1106543, 7082055) != true { if limiter.IntersectsBuffer(g, 9.94, 53.53) != true {
t.Fatal() t.Fatal()
} }
if limiter.IntersectsBuffer(g, 1006543, 7082055) != false { if limiter.IntersectsBuffer(g, 9.04, 53.53) != false {
t.Fatal() t.Fatal()
} }
} }
func TestSplitParams(t *testing.T) { func TestSplitParams(t *testing.T) {
@ -373,7 +375,7 @@ func TestSplitParams(t *testing.T) {
func BenchmarkClipper(b *testing.B) { func BenchmarkClipper(b *testing.B) {
g := geos.NewGeos() g := geos.NewGeos()
defer g.Finish() defer g.Finish()
limiter, err := NewFromGeoJson("./hamburg_clip.geojson") limiter, err := NewFromGeoJSON("./clipping.geojson", 1.0, 3857)
if err != nil { if err != nil {
b.Fatal(err) b.Fatal(err)
} }

View File

@ -40,9 +40,10 @@ func Import() {
if (config.ImportOptions.Write || config.ImportOptions.Read != "") && config.BaseOptions.LimitTo != "" { if (config.ImportOptions.Write || config.ImportOptions.Read != "") && config.BaseOptions.LimitTo != "" {
var err error var err error
step := log.StartStep("Reading limitto geometries") step := log.StartStep("Reading limitto geometries")
geometryLimiter, err = limit.NewFromGeoJsonWithBuffered( geometryLimiter, err = limit.NewFromGeoJSON(
config.BaseOptions.LimitTo, config.BaseOptions.LimitTo,
config.BaseOptions.LimitToCacheBuffer, config.BaseOptions.LimitToCacheBuffer,
config.BaseOptions.Srid,
) )
if err != nil { if err != nil {
log.Fatal(err) log.Fatal(err)

View File

@ -15,7 +15,6 @@ import (
"github.com/omniscale/imposm3/logging" "github.com/omniscale/imposm3/logging"
"github.com/omniscale/imposm3/mapping" "github.com/omniscale/imposm3/mapping"
"github.com/omniscale/imposm3/parser/pbf" "github.com/omniscale/imposm3/parser/pbf"
"github.com/omniscale/imposm3/proj"
"github.com/omniscale/imposm3/stats" "github.com/omniscale/imposm3/stats"
"github.com/omniscale/imposm3/util" "github.com/omniscale/imposm3/util"
) )
@ -193,9 +192,7 @@ func ReadPbf(cache *osmcache.OSMCache, progress *stats.Statistics,
} }
if withLimiter { if withLimiter {
for i, _ := range nds { for i, _ := range nds {
nd := element.Node{Long: nds[i].Long, Lat: nds[i].Lat} if !limiter.IntersectsBuffer(g, nds[i].Long, nds[i].Lat) {
proj.NodeToMerc(&nd)
if !limiter.IntersectsBuffer(g, nd.Long, nd.Lat) {
skip += 1 skip += 1
nds[i].Id = osmcache.SKIP nds[i].Id = osmcache.SKIP
} else { } else {
@ -228,9 +225,7 @@ func ReadPbf(cache *osmcache.OSMCache, progress *stats.Statistics,
numWithTags += 1 numWithTags += 1
} }
if withLimiter { if withLimiter {
nd := element.Node{Long: nds[i].Long, Lat: nds[i].Lat} if !limiter.IntersectsBuffer(g, nds[i].Long, nds[i].Lat) {
proj.NodeToMerc(&nd)
if !limiter.IntersectsBuffer(g, nd.Long, nd.Lat) {
nds[i].Id = osmcache.SKIP nds[i].Id = osmcache.SKIP
} }
} }

View File

@ -1,8 +0,0 @@
{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3857" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -19926188.851995971053839, -19971868.880408588796854 ], [ -19926188.851995971053839, 19971868.880408562719822 ], [ 19926188.851995974779129, 19971868.880408562719822 ], [ 19926188.851995974779129, -19971868.880408588796854 ], [ -19926188.851995971053839, -19971868.880408588796854 ] ] ] } }
]
}

8
test/clipping.geojson Normal file
View File

@ -0,0 +1,8 @@
{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -179.0, -85.0 ], [ -179.0, 85.0 ], [ 179.0, 85.0 ], [ 179.0, -85.0 ], [ -179.0, -85.0 ] ] ] } }
]
}

View File

@ -170,7 +170,7 @@ def imposm3_update(db_conf, osc, mapping_file):
print subprocess.check_output(( print subprocess.check_output((
"../imposm3 diff -connection %s" "../imposm3 diff -connection %s"
" -cachedir %s" " -cachedir %s"
" -limitto clipping-3857.geojson" " -limitto clipping.geojson"
" -dbschema-production " + TEST_SCHEMA_PRODUCTION + " -dbschema-production " + TEST_SCHEMA_PRODUCTION +
" -mapping %s %s") % ( " -mapping %s %s") % (
conn, tmpdir, mapping_file, osc, conn, tmpdir, mapping_file, osc,