Skip to content

Commit

Permalink
Added basic support for 4326 Datasource.
Browse files Browse the repository at this point in the history
Add functions to convert a geometry in 4326 to WebMercator and back. Added support to the postgis provider to do the conversion.
I don't think this is the right place for the conversion to take place, but is good enought for now. It should live in the feature
but that requires the features to know the encoding of the geometries, instead of assuming WebMercator.

I'm using the more correct rather then the faster conversion algo.

Relates to issue go-spatial#47
  • Loading branch information
gdey committed Jul 29, 2016
1 parent f4dab5e commit 2de9474
Show file tree
Hide file tree
Showing 5 changed files with 290 additions and 5 deletions.
158 changes: 155 additions & 3 deletions basic/geometry_math.go
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,15 @@ import (
"fmt"

"github.com/terranodo/tegola"
"github.com/terranodo/tegola/maths/webmercator"
)

// RehomeGeometry will make sure to normalize all points to the ulx and uly coordinates.
//TODO I don't think this is needed anymore. Marking for Cleanup. — gdey
func RehomeGeometry(geometery tegola.Geometry, ulx, uly float64) (tegola.Geometry, error) {
switch geo := geometery.(type) {
func RehomeGeometry(geometry tegola.Geometry, ulx, uly float64) (tegola.Geometry, error) {
switch geo := geometry.(type) {
default:
return nil, fmt.Errorf("Unknown Geometry: %+v", geometery)
return nil, fmt.Errorf("Unknown Geometry: %+v", geometry)
case tegola.Point:
return &Point{geo.X() - ulx, geo.Y() - uly}, nil
case tegola.Point3:
Expand Down Expand Up @@ -63,3 +64,154 @@ func RehomeGeometry(geometery tegola.Geometry, ulx, uly float64) (tegola.Geometr
return &mpoly, nil
}
}

func ApplyToPoints(geometry tegola.Geometry, f func(coords ...float64) ([]float64, error)) (tegola.Geometry, error) {
switch geo := geometry.(type) {
default:
return nil, fmt.Errorf("Unknown Geometry: %+v", geometry)
case tegola.Point:
c, err := f(geo.X(), geo.Y())
if err != nil {
return nil, err
}
if len(c) < 2 {
return nil, fmt.Errorf("Function did not return minimum number of coordinates got %v expected 2", len(c))
}
return &Point{c[0], c[1]}, nil
case tegola.Point3:
c, err := f(geo.X(), geo.Y(), geo.Z())
if err != nil {
return nil, err
}
if len(c) < 3 {
return nil, fmt.Errorf("Function did not return minimum number of coordinates got %v expected 3", len(c))
}
return &Point3{c[0], c[1], c[2]}, nil
case tegola.MultiPoint:
var pts MultiPoint
for _, pt := range geo.Points() {
c, err := f(pt.X(), pt.Y())
if err != nil {
return nil, err
}
if len(c) < 2 {
return nil, fmt.Errorf("Function did not return minimum number of coordinates got %v expected 2", len(c))
}
pts = append(pts, Point{c[0], c[1]})
}
return pts, nil
case tegola.LineString:
var line Line
for _, ptGeo := range geo.Subpoints() {
c, err := f(ptGeo.X(), ptGeo.Y())
if err != nil {
return nil, err
}
if len(c) < 2 {
return nil, fmt.Errorf("Function did not return minimum number of coordinates got %v expected 2", len(c))
}
line = append(line, Point{c[0], c[1]})
}
return &line, nil
case tegola.MultiLine:
var line MultiLine
for i, lineGeo := range geo.Lines() {
geoLine, err := ApplyToPoints(lineGeo, f)
if err != nil {
return nil, fmt.Errorf("Got error converting line(%v) of multiline: %v", i, err)
}
ln, _ := geoLine.(Line)
line = append(line, ln)
}
return line, nil
case tegola.Polygon:
var poly Polygon
for i, line := range geo.Sublines() {
geoLine, err := ApplyToPoints(line, f)
if err != nil {
return nil, fmt.Errorf("Got error converting line(%v) of polygon: %v", i, err)
}
ln, _ := geoLine.(Line)
poly = append(poly, ln)
}
return &poly, nil
case tegola.MultiPolygon:
var mpoly MultiPolygon
for i, poly := range geo.Polygons() {
geoPoly, err := ApplyToPoints(poly, f)
if err != nil {
return nil, fmt.Errorf("Got error converting poly(%v) of multipolygon: %v", i, err)
}
p, _ := geoPoly.(Polygon)
mpoly = append(mpoly, p)
}
return &mpoly, nil
}
}

func CloneGeometry(geometry tegola.Geometry) (tegola.Geometry, error) {
switch geo := geometry.(type) {
default:
return nil, fmt.Errorf("Unknown Geometry: %+v", geometry)
case tegola.Point:
return &Point{geo.X(), geo.Y()}, nil
case tegola.Point3:
return &Point3{geo.X(), geo.Y(), geo.Z()}, nil
case tegola.MultiPoint:
var pts MultiPoint
for _, pt := range geo.Points() {
pts = append(pts, Point{pt.X(), pt.Y()})
}
return pts, nil
case tegola.LineString:
var line Line
for _, ptGeo := range geo.Subpoints() {
line = append(line, Point{ptGeo.X(), ptGeo.Y()})
}
return &line, nil
case tegola.MultiLine:
var line MultiLine
for i, lineGeo := range geo.Lines() {
geoLine, err := CloneGeometry(lineGeo)
if err != nil {
return nil, fmt.Errorf("Got error converting line(%v) of multiline: %v", i, err)
}
ln, _ := geoLine.(Line)
line = append(line, ln)
}
return line, nil
case tegola.Polygon:
var poly Polygon
for i, line := range geo.Sublines() {
geoLine, err := CloneGeometry(line)
if err != nil {
return nil, fmt.Errorf("Got error converting line(%v) of polygon: %v", i, err)
}
ln, _ := geoLine.(Line)
poly = append(poly, ln)
}
return &poly, nil
case tegola.MultiPolygon:
var mpoly MultiPolygon
for i, poly := range geo.Polygons() {
geoPoly, err := CloneGeometry(poly)
if err != nil {
return nil, fmt.Errorf("Got error converting poly(%v) of multipolygon: %v", i, err)
}
p, _ := geoPoly.(Polygon)
mpoly = append(mpoly, p)
}
return &mpoly, nil
}
}
func ToWebMercator(SRID int, geometry tegola.Geometry) (tegola.Geometry, error) {
// We are converting to the same set so, just return a clone.
switch SRID {
default:
return nil, fmt.Errorf("Don't know how to convert from %v to %v.", tegola.WebMercator)
case tegola.WebMercator:
return CloneGeometry(geometry)
case tegola.WGS84:
return ApplyToPoints(geometry, webmercator.ToXY)
}
}
7 changes: 6 additions & 1 deletion geometry_math.go
Original file line number Diff line number Diff line change
@@ -1,7 +1,11 @@
package tegola

import "math"
const (
WebMercator = 3857
WGS84 = 4326
)

/*
// AreaOfPolygon will calculate the Area of a polygon using the surveyor's formula
// (https://en.wikipedia.org/wiki/Shoelace_formula)
func AreaOfPolygon(p Polygon) (area float64) {
Expand All @@ -17,3 +21,4 @@ func AreaOfPolygon(p Polygon) (area float64) {
}
return math.Abs(area) / 2.0
}
*/
43 changes: 43 additions & 0 deletions maths/maths.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
/*
Package math contins generic math functions that we need for doing transforms.
this package will augment the go math library.
*/
package maths

import (
"math"

"github.com/terranodo/tegola"
)

const (
WebMercator = tegola.WebMercator
WGS84 = tegola.WGS84
Deg2Rad = math.Pi / 180
Rad2Deg = 180 / math.Pi
PiDiv2 = math.Pi / 2.0
)

// AreaOfPolygon will calculate the Area of a polygon using the surveyor's formula
// (https://en.wikipedia.org/wiki/Shoelace_formula)
func AreaOfPolygon(p tegola.Polygon) (area float64) {
var points []tegola.Point
for _, l := range p.Sublines() {
points = append(points, l.Subpoints()...)
}
n := len(points)
for i := range points {
j := (i + 1) % n
area += points[i].X() * points[j].Y()
area -= points[j].X() * points[i].Y()
}
return math.Abs(area) / 2.0
}

func RadToDeg(rad float64) float64 {
return rad * Rad2Deg
}

func DegToRad(deg float64) float64 {
return deg * Deg2Rad
}
76 changes: 76 additions & 0 deletions maths/webmercator/main.go
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
/*
Package webmercator does the translation to and from WebMercator and WGS84
Gotten from: http://wiki.openstreetmap.org/wiki/Mercator#C.23
*/
package webmercator

import (
"fmt"
"math"

"github.com/terranodo/tegola/maths"
)

const (
RMajor = 6378137.0
RMinor = 6356752.3142
Ratio = RMinor / RMajor
)

var Eccent float64
var Com float64

func init() {
Eccent = math.Sqrt(1.0 - (Ratio * Ratio))
Com = 0.5 * Eccent
}

func LonToX(lon float64) float64 {
return RMajor * maths.DegToRad(lon)
}

func con(phi float64) float64 {
v := Eccent * math.Sin(phi)
return math.Pow(((1.0 - v) / (1.0 + v)), Com)
}

func LatToY(lat float64) float64 {
lat = math.Min(89.5, math.Max(lat, -89.5))
phi := maths.DegToRad(lat)
ts := math.Tan(0.5*((math.Pi*0.5)-phi)) / con(phi)
return 0 - RMajor*math.Log(ts)
}

func XToLon(x float64) float64 {
return maths.RadToDeg(x) / RMajor
}

func YToLat(y float64) float64 {
ts := math.Exp(-y / RMajor)
phi := maths.PiDiv2 - 2*math.Atan(ts)
dphi := 1.0
i := 0
for (math.Abs(dphi) > 0.000000001) && (i < 15) {
dphi = maths.PiDiv2 - 2*math.Atan(ts*con(phi)) - phi
phi += dphi
i++
}
return maths.RadToDeg(phi)
}

func ToLonLat(c ...float64) ([]float64, error) {
if len(c) < 2 {
return c, fmt.Errorf("Coords should have at least 2 coords")
}
crds := []float64{XToLon(c[0]), YToLat(c[1])}
crds = append(crds, c[2:]...)
return crds, nil
}
func ToXY(c ...float64) ([]float64, error) {
if len(c) < 2 {
return c, fmt.Errorf("Coords should have at least 2 coords")
}
crds := []float64{LonToX(c[0]), LatToY(c[1])}
crds = append(crds, c[2:]...)
return crds, nil
}
11 changes: 10 additions & 1 deletion provider/postgis/postgis.go
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ import (
"github.com/jackc/pgx"

"github.com/terranodo/tegola"
"github.com/terranodo/tegola/basic"
"github.com/terranodo/tegola/mvt"
"github.com/terranodo/tegola/mvt/provider"
"github.com/terranodo/tegola/util/dict"
Expand Down Expand Up @@ -48,7 +49,7 @@ const fldsSQL = "SELECT * FROM %[1]v LIMIT 0;"

const Name = "postgis"
const DefaultPort = 5432
const DefaultSRID = 3857
const DefaultSRID = tegola.WebMercator
const DefaultMaxConn = 5

const (
Expand Down Expand Up @@ -324,6 +325,14 @@ func (p Provider) MVTLayer(layerName string, tile tegola.Tile, tags map[string]i
if geom, err = wkb.DecodeBytes(geobytes); err != nil {
return nil, fmt.Errorf("Was unable to decode geometry field(%v) into wkb for layer %v.", plyr.GeomFieldName, layerName)
}
// TODO: Need to move this from being the responsiblity of the provider to the responsibility of the feature. But that means a feature should know
// how the points are encoded.
if p.srid != DefaultSRID {
// We need to convert our points to Webmercator.
if geom, err = basic.ToWebMercator(p.srid, geom); err != nil {
return nil, fmt.Errorf("Was unable to transform geometry to webmercator from SRID(%v) for layer %v.", p.srid, layerName)
}
}
case plyr.IDFieldName:
switch aval := v.(type) {
case int64:
Expand Down

0 comments on commit 2de9474

Please sign in to comment.