From e533068a6126717f45e5fe4e53f2d7a8d02d6d69 Mon Sep 17 00:00:00 2001 From: Peter Kleiweg Date: Sat, 31 Aug 2019 14:23:14 +0200 Subject: [PATCH] more work done, docs added --- README.md | 15 ++- v5/doc.go | 11 ++ v5/proj.c | 12 +- v5/proj.go | 333 ++++++++++++++++++++++++++++++++++++++++++++++++ v5/proj_go.h | 6 + v5/proj_test.go | 95 ++++++++++++++ v5/utm.go | 45 +++++++ v5/v5.go | 140 -------------------- 8 files changed, 514 insertions(+), 143 deletions(-) create mode 100644 v5/doc.go create mode 100644 v5/proj.go create mode 100644 v5/proj_test.go create mode 100644 v5/utm.go delete mode 100644 v5/v5.go diff --git a/README.md b/README.md index c904102..bdd1ab8 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,16 @@ -The [Go](http://golang.org/) package _proj_ provides an interface to the Cartographic Projections Library [PROJ](https://proj.org/). +The [Go](http://golang.org/) package _proj_ provides a limited interface to the Cartographic Projections Library [PROJ](https://proj.org/). + +This package supports PROJ version 5 and above. + +For PROJ.4, see: https://github.com/pebbe/go-proj-4 Keywords: cartography, cartographic projection -*Work in progress* +## Install + + go get github.com/pebbe/proj/v5 + +## Docs + + * [package help v5](http://godoc.org/github.com/pebbe/proj/v5) + diff --git a/v5/doc.go b/v5/doc.go new file mode 100644 index 0000000..b34ee8d --- /dev/null +++ b/v5/doc.go @@ -0,0 +1,11 @@ +/* +Package proj provides an interface to the Cartographic Projections Library PROJ [cartography]. + +See: https://proj.org/ + +This package supports PROJ version 5 and above. + +For PROJ.4, see: https://github.com/pebbe/go-proj-4 + +*/ +package proj diff --git a/v5/proj.c b/v5/proj.c index ef6f137..8313c75 100644 --- a/v5/proj.c +++ b/v5/proj.c @@ -1,7 +1,7 @@ #include "proj_go.h" int pjnull(PJ *pj) { - return pj == 0 ? 1 : 0; + return pj == 0 ? 0 : 1; } void trans(PJ *pj, PJ_DIRECTION direction, double u1, double v1, double w1, double t1, double *u2, double *v2, double *w2, double *t2) { @@ -18,3 +18,13 @@ void trans(PJ *pj, PJ_DIRECTION direction, double u1, double v1, double w1, doub *w2 = co2.uvwt.w; *t2 = co2.uvwt.t; } + +PJ_COORD uvwt(double u, double v, double w, double t) { + PJ_COORD + c; + c.uvwt.u = u; + c.uvwt.v = v; + c.uvwt.w = w; + c.uvwt.t = t; + return c; +} diff --git a/v5/proj.go b/v5/proj.go new file mode 100644 index 0000000..6c79603 --- /dev/null +++ b/v5/proj.go @@ -0,0 +1,333 @@ +package proj + +/* +#cgo darwin pkg-config: proj +#cgo !darwin LDFLAGS: -lproj +#include "proj_go.h" +*/ +import "C" + +import ( + "errors" + "math" + "runtime" + "unsafe" +) + +type Context struct { + pj_context *C.PJ_CONTEXT + opened bool + counter uint64 + projections map[uint64]*PJ +} + +// A projection object +type PJ struct { + pj *C.PJ + context *Context + index uint64 + opened bool +} + +type LibInfo struct { + Major int // Major version number. + Minor int // Minor version number. + Patch int // Patch level of release. + Release string // Release info. Version number and release date, e.g. “Rel. 4.9.3, 15 August 2016”. + Version string // Text representation of the full version number, e.g. “4.9.3”. + Searchpath string // Search path for PROJ. List of directories separated by semicolons (Windows) or colons (non-Windows). +} + +type ProjInfo struct { + ID string // Short ID of the operation the PJ object is based on, that is, what comes afther the +proj= in a proj-string, e.g. “merc”. + Description string // Long describes of the operation the PJ object is based on, e.g. “Mercator Cyl, Sph&Ell lat_ts=”. + Definition string // The proj-string that was used to create the PJ object with, e.g. “+proj=merc +lat_0=24 +lon_0=53 +ellps=WGS84”. + HasInverse bool // True if an inverse mapping of the defined operation exists, + Accuracy float64 // Expected accuracy of the transformation. -1 if unknown. +} + +// The direction of a transformation +type Direction C.PJ_DIRECTION + +const ( + Fwd = Direction(C.PJ_FWD) // Forward transformation + Ident = Direction(C.PJ_IDENT) // Do nothing + Inv = Direction(C.PJ_INV) // Inverse transformation +) + +var ( + errContextClosed = errors.New("Context is closed") + errDataSizeMismatch = errors.New("Data size mismatch") + errMissingData = errors.New("Missing data") + errProjectionClosed = errors.New("Projection is closed") +) + +// Create a context +func NewContext() *Context { + ctx := Context{ + pj_context: C.proj_context_create(), + counter: 0, + projections: make(map[uint64]*PJ), + opened: true, + } + runtime.SetFinalizer(&ctx, (*Context).Close) + return &ctx +} + +// Close a context +func (ctx *Context) Close() { + if ctx.opened { + indexen := make([]uint64, 0, len(ctx.projections)) + for i := range ctx.projections { + indexen = append(indexen, i) + } + for _, i := range indexen { + p := ctx.projections[i] + if p.opened { + C.proj_destroy(p.pj) + p.context = nil + p.opened = false + } + delete(ctx.projections, i) + } + + C.proj_context_destroy(ctx.pj_context) + ctx.pj_context = nil + ctx.opened = false + } +} + +// Create a transformation object +func (ctx *Context) Create(definition string) (*PJ, error) { + if !ctx.opened { + return &PJ{}, errContextClosed + } + + cs := C.CString(definition) + defer C.free(unsafe.Pointer(cs)) + pj := C.proj_create(ctx.pj_context, cs) + if C.pjnull(pj) == 0 { + errno := C.proj_context_errno(ctx.pj_context) + err := C.GoString(C.proj_errno_string(errno)) + return &PJ{}, errors.New(err) + } + + p := PJ{ + opened: true, + context: ctx, + index: ctx.counter, + pj: pj, + } + ctx.projections[ctx.counter] = &p + ctx.counter++ + + runtime.SetFinalizer(&p, (*PJ).Close) + return &p, nil +} + +// Close a transformation object +func (p *PJ) Close() { + if p.opened { + C.proj_destroy(p.pj) + if p.context.opened { + delete(p.context.projections, p.index) + } + p.context = nil + p.opened = false + } +} + +// Get information about the transformation object +func (p *PJ) Info() (ProjInfo, error) { + if !p.opened { + return ProjInfo{}, errProjectionClosed + } + info := C.proj_pj_info(p.pj) + hasInv := false + if info.has_inverse != 0 { + hasInv = true + } + return ProjInfo{ + ID: C.GoString(info.id), + Description: C.GoString(info.description), + Definition: C.GoString(info.definition), + HasInverse: hasInv, + Accuracy: float64(info.accuracy), + }, nil +} + +// Transform a single transformation +func (p *PJ) Trans(direction Direction, u1, v1, w1, t1 float64) (u2, v2, w2, t2 float64, err error) { + if !p.opened { + return 0, 0, 0, 0, errProjectionClosed + } + + var u, v, w, t C.double + C.trans(p.pj, C.PJ_DIRECTION(direction), C.double(u1), C.double(v1), C.double(w1), C.double(t1), &u, &v, &w, &t) + + e := C.proj_errno(p.pj) + if e != 0 { + return 0, 0, 0, 0, errors.New(C.GoString(C.proj_errno_string(e))) + } + + return float64(u), float64(v), float64(w), float64(t), nil +} + +func (p *PJ) TransSlice(direction Direction, u1, v1, w1, t1 []float64) (u2, v2, w2, t2 []float64, err error) { + if !p.opened { + return nil, nil, nil, nil, errProjectionClosed + } + if u1 == nil || v1 == nil { + return nil, nil, nil, nil, errMissingData + } + + var un, vn, wn, tn int + var u, v, w, t []C.double + var up, vp, wp, tp *C.double + var unc, vnc, wnc, tnc C.size_t + + un = len(u1) + unc = C.size_t(un) + vn = len(v1) + vnc = C.size_t(vn) + if w1 != nil { + wn = len(w1) + wnc = C.size_t(wn) + if t1 != nil { + tn = len(t1) + tnc = C.size_t(tn) + } + } + r := []int{un, vn, tn, wn} + var n int + for _, i := range r { + if i > n { + n = i + } + } + for _, i := range r { + if i > 1 && i < n { + return nil, nil, nil, nil, errDataSizeMismatch + } + } + + u = make([]C.double, un) + up = &u[0] + for i := 0; i < un; i++ { + u[i] = C.double(u1[i]) + } + v = make([]C.double, vn) + vp = &v[0] + for i := 0; i < vn; i++ { + v[i] = C.double(v1[i]) + } + if w1 != nil { + w = make([]C.double, wn) + wp = &w[0] + for i := 0; i < wn; i++ { + w[i] = C.double(w1[i]) + } + if t1 != nil { + t = make([]C.double, tn) + tp = &t[0] + for i := 0; i < tn; i++ { + t[i] = C.double(t1[i]) + } + } + } + + st := C.size_t(C.sizeof_double) + + C.proj_trans_generic( + p.pj, + C.PJ_DIRECTION(direction), + up, st, unc, + vp, st, vnc, + wp, st, wnc, + tp, st, tnc) + + e := C.proj_errno(p.pj) + if e != 0 { + return nil, nil, nil, nil, errors.New(C.GoString(C.proj_errno_string(e))) + } + + u2 = make([]float64, un) + for i := 0; i < un; i++ { + u2[i] = float64(u[i]) + } + v2 = make([]float64, vn) + for i := 0; i < vn; i++ { + v2[i] = float64(v[i]) + } + if w1 != nil { + w2 = make([]float64, wn) + for i := 0; i < wn; i++ { + w2[i] = float64(w[i]) + } + if t1 != nil { + t2 = make([]float64, tn) + for i := 0; i < tn; i++ { + t2[i] = float64(t[i]) + } + } + } + + return +} + +// Calculate geodesic distance between two points in geodetic coordinates +// +// The calculated distance is between the two points located on the ellipsoid +func (p *PJ) Dist(u1, v1, u2, v2 float64) (float64, error) { + if !p.opened { + return 0, errProjectionClosed + } + a := C.uvwt(C.double(u1), C.double(v1), 0, 0) + b := C.uvwt(C.double(u2), C.double(v2), 0, 0) + d := C.proj_lp_dist(p.pj, a, b) + e := C.proj_errno(p.pj) + if e != 0 { + return 0, errors.New(C.GoString(C.proj_errno_string(e))) + } + return float64(d), nil +} + +// Calculate geodesic distance between two points in geodetic coordinates +// +// Similar to Dist() but also takes the height above the ellipsoid into account +func (p *PJ) Dist3(u1, v1, w1, u2, v2, w2 float64) (float64, error) { + if !p.opened { + return 0, errProjectionClosed + } + a := C.uvwt(C.double(u1), C.double(v1), C.double(w1), 0) + b := C.uvwt(C.double(u2), C.double(v2), C.double(w2), 0) + d := C.proj_lpz_dist(p.pj, a, b) + e := C.proj_errno(p.pj) + if e != 0 { + return 0, errors.New(C.GoString(C.proj_errno_string(e))) + } + return float64(d), nil +} + +// Get information about the current instance of the PROJ library +func Info() LibInfo { + info := C.proj_info() + return LibInfo{ + Major: int(info.major), + Minor: int(info.minor), + Patch: int(info.patch), + Release: C.GoString(info.release), + Version: C.GoString(info.version), + Searchpath: C.GoString(info.searchpath), + } +} + +// Convert degrees to radians +func DegToRad(deg float64) float64 { + return deg / 180 * math.Pi +} + +// Convert radians to degrees +func RadToDeg(rad float64) float64 { + return rad / math.Pi * 180 +} diff --git a/v5/proj_go.h b/v5/proj_go.h index 3f89a3d..d2b1ff0 100644 --- a/v5/proj_go.h +++ b/v5/proj_go.h @@ -1,5 +1,11 @@ +#ifndef _PROJ_GO_H_ +#define _PROJ_GO_H_ + #include #include int pjnull(PJ *pj); void trans(PJ *pj, PJ_DIRECTION direction, double u1, double v1, double w1, double t1, double *u2, double *v2, double *w2, double *t2); +PJ_COORD uvwt(double u, double v, double w, double t); + +#endif diff --git a/v5/proj_test.go b/v5/proj_test.go new file mode 100644 index 0000000..a440a52 --- /dev/null +++ b/v5/proj_test.go @@ -0,0 +1,95 @@ +package proj_test + +import ( + "github.com/pebbe/proj/v5" + + "fmt" + "testing" +) + +func TestLatlongToMerc(t *testing.T) { + ctx := proj.NewContext() + defer ctx.Close() + + ll, err := ctx.Create("+proj=latlong +datum=WGS84") + if err != nil { + t.Fatal(err) + } + defer ll.Close() + + merc, err := ctx.Create("+proj=merc +ellps=clrk66 +lat_ts=33") + defer merc.Close() + if err != nil { + t.Fatal(err) + } + + u, v, _, _, err := ll.Trans(proj.Inv, proj.DegToRad(-16), proj.DegToRad(20.25), 0, 0) + if err != nil { + t.Fatal(err) + } + u, v, _, _, err = merc.Trans(proj.Fwd, u, v, 0, 0) + if err != nil { + t.Fatal(err) + } + s := fmt.Sprintf("%.2f %.2f", u, v) + s1 := "-1495284.21 1920596.79" + if s != s1 { + t.Fatalf("LatlongToMerc = %v, want %v", s, s1) + } + + pj, err := ctx.Create(` ++proj=pipeline ++step +proj=unitconvert +xy_in=deg +xy_out=rad ++step +inv +proj=latlong +datum=WGS84 ++step +proj=merc +ellps=clrk66 +lat_ts=33 +`) + defer pj.Close() + if err != nil { + t.Fatal(err) + } + x1 := []float64{-16, -10, 0, 30.4} + y1 := []float64{20.25, 25, 0, 40.8} + x2, y2, _, _, err := pj.TransSlice(proj.Fwd, x1, y1, nil, nil) + if err != nil { + t.Error(err) + } else { + s := fmt.Sprintf("[%.2f %.2f] [%.2f %.2f] [%.2f %.2f] [%.2f %.2f]", x2[0], y2[0], x2[1], y2[1], x2[2], y2[2], x2[3], y2[3]) + s1 := "[-1495284.21 1920596.79] [-934552.63 2398930.20] [0.00 0.00] [2841040.00 4159542.20]" + if s != s1 { + t.Errorf("LatlongToMerc = %v, want %v", s, s1) + } + } +} + +func TestInvalidErrorProblem(t *testing.T) { + ctx := proj.NewContext() + defer ctx.Close() + + ll, err := ctx.Create("+proj=latlong +datum=WGS84") + if err != nil { + t.Fatal(err) + } + defer ll.Close() + + merc, err := ctx.Create("+proj=merc +ellps=clrk66 +lat_ts=33") + defer merc.Close() + if err != nil { + t.Fatal(err) + } + + u, v, _, _, err := ll.Trans(proj.Inv, proj.DegToRad(3000), proj.DegToRad(500), 0, 0) + if err != nil { + t.Fatal(err) + } + _, _, _, _, err = merc.Trans(proj.Fwd, u, v, 0, 0) + if err == nil { + t.Error("err should not be nil") + } + + // Try create a new projection after an error + merc2, err := ctx.Create("+proj=merc +ellps=clrk66 +lat_ts=33") + defer merc2.Close() + if err != nil { + t.Error(err) + } +} diff --git a/v5/utm.go b/v5/utm.go new file mode 100644 index 0000000..443ce84 --- /dev/null +++ b/v5/utm.go @@ -0,0 +1,45 @@ +package proj + +import ( + "errors" +) + +// Get UTM zone for longitude and latitude in degrees. +// +// Reference: +// UTM Grid Zones of the World compiled by Alan Morton +// http://www.dmap.co.uk/utmworld.htm +func UTMzone(lng, lat float64) (xzone int, yzone string, err error) { + + if lat < -80 || lat > 84 { + err = errors.New("Arctic and antarctic region are not in UTM") + return + } + + for lng < -180 { + lng += 360 + } + for lng > 180 { + lng -= 360 + } + + xzone = 1 + int((lng + 180) / 6) + if lat > 72 && lng > 0 && lng < 42 { + if lng < 9 { + xzone = 31 + } else if lng < 21 { + xzone = 33 + } else if lng < 33 { + xzone = 35 + } else { + xzone = 37 + } + } + if lat > 56 && lat < 64 && lng > 3 && lng < 12 { + xzone = 32 + } + + yzone = string("CDEFGHJKLMNPQRSTUVWXX"[int((lat + 80) / 8)]) + + return +} diff --git a/v5/v5.go b/v5/v5.go deleted file mode 100644 index e0e5c8b..0000000 --- a/v5/v5.go +++ /dev/null @@ -1,140 +0,0 @@ -package proj - -/* -#cgo darwin pkg-config: proj -#cgo !darwin LDFLAGS: -lproj -#include "proj_go.h" -*/ -import "C" - -import ( - "errors" - "runtime" - "unsafe" -) - -type Context struct { - pj_context *C.PJ_CONTEXT - opened bool - counter uint64 - projections map[uint64]*Proj -} - -type Proj struct { - pj *C.PJ - context *Context - index uint64 - opened bool -} - -type Coord struct { - U, V, W, T float64 -} - -var ( - errContextClosed = errors.New("Context is closed") - errProjectionClosed = errors.New("Projection is closed") -) - -func NewContext() *Context { - ctx := Context{ - pj_context: C.proj_context_create(), - counter: 0, - projections: make(map[uint64]*Proj), - opened: true, - } - runtime.SetFinalizer(&ctx, (*Context).Close) - return &ctx -} - -func (ctx *Context) Close() { - if ctx.opened { - indexen := make([]uint64, 0, len(ctx.projections)) - for i := range ctx.projections { - indexen = append(indexen, i) - } - for _, i := range indexen { - p := ctx.projections[i] - if p.opened { - C.proj_destroy(p.pj) - p.context = nil - p.opened = false - } - delete(ctx.projections, i) - } - - C.proj_context_destroy(ctx.pj_context) - ctx.pj_context = nil - ctx.opened = false - } -} - -func (ctx *Context) Create(definition string) (*Proj, error) { - if !ctx.opened { - return nil, errContextClosed - } - - cs := C.CString(definition) - defer C.free(unsafe.Pointer(cs)) - pj := C.proj_create(ctx.pj_context, cs) - if C.pjnull(pj) != 0 { - errno := C.proj_context_errno(ctx.pj_context) - err := C.GoString(C.proj_errno_string(errno)) - return nil, errors.New(err) - } - - p := Proj{ - opened: true, - context: ctx, - index: ctx.counter, - pj: pj, - } - ctx.projections[ctx.counter] = &p - ctx.counter++ - - runtime.SetFinalizer(&p, (*Proj).Close) - return &p, nil -} - -func (p *Proj) Close() { - if p.opened { - C.proj_destroy(p.pj) - if p.context.opened { - delete(p.context.projections, p.index) - } - p.context = nil - p.opened = false - } -} - -func (p *Proj) Fwd(coord Coord) (Coord, error) { - return p.trans(coord, false) -} - -func (p *Proj) Inv(coord Coord) (Coord, error) { - return p.trans(coord, true) -} - -func (p *Proj) trans(coord Coord, inverse bool) (Coord, error) { - if !p.opened { - return Coord{}, errProjectionClosed - } - - var direction C.PJ_DIRECTION - if inverse { - direction = C.PJ_INV - } else { - direction = C.PJ_FWD - } - - var u, v, w, t C.double - C.trans(p.pj, direction, C.double(coord.U), C.double(coord.V), C.double(coord.W), C.double(coord.T), &u, &v, &w, &t) - - coord2 := Coord{ - U: float64(u), - V: float64(v), - W: float64(w), - T: float64(t), - } - return coord2, nil -}