package crs
import (
"math"
)
type AlbersEqualArea struct {
Lonf float64
Latf float64
Sp1 float64
Sp2 float64
Eastf float64
Northf float64
}
func (cs AlbersEqualArea) String() string {
return build("albers_equal_area").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"sp1", cs.Sp1,
"sp2", cs.Sp2,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (aea AlbersEqualArea) consts(s Spheroid) (n, c, rf float64) {
phif := radian(aea.Latf)
phi1 := radian(aea.Sp1)
phi2 := radian(aea.Sp2)
e, e2 := s.E(), s.E2()
alphaf := authalicQ(math.Sin(phif), e, e2)
alpha1 := authalicQ(math.Sin(phi1), e, e2)
alpha2 := authalicQ(math.Sin(phi2), e, e2)
m1 := math.Cos(phi1) / math.Sqrt(1-e2*sin2(phi1))
m2 := math.Cos(phi2) / math.Sqrt(1-e2*sin2(phi2))
if math.Abs(phi1-phi2) < 1e-10 {
n = math.Sin(phi1)
} else {
n = (m1*m1 - m2*m2) / (alpha2 - alpha1)
}
c = m1*m1 + n*alpha1
rf = (s.A * math.Sqrt(c-n*alphaf)) / n
return n, c, rf
}
func (aea AlbersEqualArea) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
lambdaf := radian(aea.Lonf)
n, c, rf := aea.consts(s)
lambda := radian(lon)
phi := radian(lat)
alpha := authalicQ(math.Sin(phi), s.E(), s.E2())
theta := n * (lambda - lambdaf)
r := (s.A * math.Sqrt(c-n*alpha)) / n
east := aea.Eastf + r*math.Sin(theta)
north := aea.Northf + rf - r*math.Cos(theta)
return east, north, h
}
func (aea AlbersEqualArea) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
lambdaf := radian(aea.Lonf)
n, c, rf := aea.consts(s)
qp := authalicQ(1, s.E(), s.E2())
ri := math.Sqrt(intPow(east-aea.Eastf, 2) + intPow(rf-(north-aea.Northf), 2))
alphai := (c - (intPow(ri, 2) * intPow(n, 2) / s.A2())) / n
betai := math.Asin(clamp(alphai/qp, -1, 1))
var theta float64
if n > 0 {
theta = math.Atan2((east - aea.Eastf), (rf - (north - aea.Northf)))
} else {
theta = math.Atan2(-(east - aea.Eastf), -(rf - (north - aea.Northf)))
}
phi := authalicToGeodetic(betai, s)
lambda := lambdaf + (theta / n)
return degree(lambda), degree(phi), h
}
package crs
import (
"math"
)
type AzimuthalEquidistant struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs AzimuthalEquidistant) String() string {
return build("azimuthal_equidistant").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs AzimuthalEquidistant) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
a := s.A
phi0 := radian(cs.Latf)
lam0 := radian(cs.Lonf)
phi := radian(lat)
lam := radian(lon)
sinPhi0, cosPhi0 := math.Sincos(phi0)
sinPhi, cosPhi := math.Sincos(phi)
dLam := lam - lam0
cosC := sinPhi0*sinPhi + cosPhi0*cosPhi*math.Cos(dLam)
c := math.Acos(clamp(cosC, -1, 1))
if math.Abs(c) < 1e-14 {
return cs.Eastf, cs.Northf, h
}
k := c / math.Sin(c)
east := cs.Eastf + a*k*cosPhi*math.Sin(dLam)
north := cs.Northf + a*k*(cosPhi0*sinPhi-sinPhi0*cosPhi*math.Cos(dLam))
return east, north, h
}
func (cs AzimuthalEquidistant) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
a := s.A
phi0 := radian(cs.Latf)
lam0 := radian(cs.Lonf)
x := (east - cs.Eastf) / a
y := (north - cs.Northf) / a
rho := math.Hypot(x, y)
if rho < 1e-14 {
return cs.Lonf, cs.Latf, h
}
c := rho
sinC, cosC := math.Sincos(c)
sinPhi0, cosPhi0 := math.Sincos(phi0)
phi := math.Asin(clamp(cosC*sinPhi0+y*sinC*cosPhi0/rho, -1, 1))
lam := lam0 + math.Atan2(x*sinC, rho*cosPhi0*cosC-y*sinPhi0*sinC)
return degree(lam), degree(phi), h
}
package crs
import (
"fmt"
"math"
)
var World = BoundingBox{
MinLon: -180,
MinLat: -90,
MaxLon: 180,
MaxLat: 90,
}
type BoundingBox struct {
MinLon float64
MinLat float64
MaxLon float64
MaxLat float64
}
func (b BoundingBox) String() string {
return fmt.Sprintf("%s,%s,%s,%s", formatFloat(b.MinLon), formatFloat(b.MinLat), formatFloat(b.MaxLon), formatFloat(b.MaxLat))
}
// Contains reports whether (lon, lat) lies in the box. Longitude must be east of
// Greenwich (EPSG area-of-use convention). When the point is relative to a local
// prime meridian, pass lon + datum.primeMeridianLongitude().
func (b BoundingBox) Contains(lon, lat float64) bool {
if b.MaxLat < b.MinLat {
return false
}
if lat < b.MinLat || lat > b.MaxLat {
return false
}
if b.MinLon <= b.MaxLon {
return lon >= b.MinLon && lon <= b.MaxLon
}
return lon >= b.MinLon || lon <= b.MaxLon
}
func (b BoundingBox) Area() float64 {
lonSpan := b.MaxLon - b.MinLon
if lonSpan < 0 {
lonSpan = (180 - b.MinLon) + (b.MaxLon + 180)
}
latSpan := b.MaxLat - b.MinLat
if lonSpan <= 0 || latSpan <= 0 {
return math.MaxFloat64
}
return lonSpan * latSpan
}
func lonIntervals(minLon, maxLon float64) [][2]float64 {
if minLon <= maxLon {
return [][2]float64{{minLon, maxLon}}
}
return [][2]float64{{minLon, 180}, {-180, maxLon}}
}
func intervalsOverlap(a, b [2]float64) bool {
return a[0] <= b[1] && b[0] <= a[1]
}
func (b BoundingBox) Intersects(bbox BoundingBox) bool {
if b.MaxLat < b.MinLat || bbox.MaxLat < bbox.MinLat {
return false
}
if b.MaxLat < bbox.MinLat || bbox.MaxLat < b.MinLat {
return false
}
for _, bi := range lonIntervals(b.MinLon, b.MaxLon) {
for _, bj := range lonIntervals(bbox.MinLon, bbox.MaxLon) {
if intervalsOverlap(bi, bj) {
return true
}
}
}
return false
}
package crs
import (
"math"
)
type BonneSouthOrientated struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs BonneSouthOrientated) String() string {
return build("bonne_south_orientated").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs BonneSouthOrientated) am1(s Spheroid) float64 {
phi1 := radian(cs.Latf)
sin1 := math.Sin(phi1)
cos1 := math.Cos(phi1)
return cos1 / (math.Sqrt(1-s.E2()*sin1*sin1) * sin1)
}
func (cs BonneSouthOrientated) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
phi1 := radian(cs.Latf)
phi := radian(lat)
lam := radian(lon) - radian(cs.Lonf)
am1 := cs.am1(s)
m1 := meridianDistance(s, phi1) / s.A
m := meridianDistance(s, phi) / s.A
sinPhi, cosPhi := math.Sincos(phi)
rh := am1 + m1 - m
var x, y float64
if math.Abs(rh) > 1e-14 {
e := cosPhi * lam / (rh * math.Sqrt(1-s.E2()*sinPhi*sinPhi))
x = rh * math.Sin(e)
y = am1 - rh*math.Cos(e)
}
return cs.Eastf - s.A*x, cs.Northf - s.A*y, h
}
func (cs BonneSouthOrientated) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
am1 := cs.am1(s)
x := (cs.Eastf - east) / s.A
y := am1 - (cs.Northf-north)/s.A
phi1 := radian(cs.Latf)
rh := math.Copysign(math.Hypot(x, y), phi1)
m1 := meridianDistance(s, phi1) / s.A
phi := footpointLatitude(s, s.A*(am1+m1-rh))
sinPhi, cosPhi := math.Sincos(phi)
var lam float64
if math.Abs(math.Abs(phi)-math.Pi/2) <= 1e-14 {
lam = 0
} else {
lm := rh * math.Sqrt(1-s.E2()*sinPhi*sinPhi) / cosPhi
if phi1 > 0 {
lam = lm * math.Atan2(x, y)
} else {
lam = lm * math.Atan2(-x, -y)
}
}
return degree(lam + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"math"
)
type CassiniSoldner struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs CassiniSoldner) String() string {
return build("cassini_soldner").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func meridianDistance(s Spheroid, phi float64) float64 {
e2 := s.E2()
e4 := e2 * e2
e6 := e4 * e2
return s.A * ((1-e2/4-3*e4/64-5*e6/256)*phi -
(3*e2/8+3*e4/32+45*e6/1024)*math.Sin(2*phi) +
(15*e4/256+45*e6/1024)*math.Sin(4*phi) -
(35*e6/3072)*math.Sin(6*phi))
}
func footpointLatitude(s Spheroid, m float64) float64 {
e2 := s.E2()
e4 := e2 * e2
e6 := e4 * e2
e1 := (1 - math.Sqrt(1-e2)) / (1 + math.Sqrt(1-e2))
e12 := e1 * e1
e13 := e12 * e1
e14 := e12 * e12
mu := m / (s.A * (1 - e2/4 - 3*e4/64 - 5*e6/256))
return mu +
(3*e1/2-27*e13/32)*math.Sin(2*mu) +
(21*e12/16-55*e14/32)*math.Sin(4*mu) +
(151*e13/96)*math.Sin(6*mu) +
(1097*e14/512)*math.Sin(8*mu)
}
func (cs CassiniSoldner) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
phi := radian(lat)
lam := radian(lon)
phi0 := radian(cs.Latf)
lam0 := radian(cs.Lonf)
sinPhi, cosPhi := math.Sincos(phi)
e2 := s.E2()
nu := s.A / math.Sqrt(1-e2*sinPhi*sinPhi)
t := math.Tan(phi)
t2 := t * t
c := e2 * cosPhi * cosPhi / (1 - e2)
a := (lam - lam0) * cosPhi
a2 := a * a
a3 := a2 * a
a4 := a2 * a2
a5 := a4 * a
east := cs.Eastf + nu*(a-t2*a3/6-(8-t2+8*c)*t2*a5/120)
north := cs.Northf + (meridianDistance(s, phi) - meridianDistance(s, phi0)) +
nu*t*(a2/2+(5-t2+6*c)*a4/24)
return east, north, h
}
func (cs CassiniSoldner) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
phi0 := radian(cs.Latf)
lam0 := radian(cs.Lonf)
e2 := s.E2()
m1 := meridianDistance(s, phi0) + (north - cs.Northf)
phi1 := footpointLatitude(s, m1)
if math.Abs(math.Abs(phi1)-math.Pi/2) < 1e-12 {
return degree(lam0), degree(phi1), h
}
sinPhi1, cosPhi1 := math.Sincos(phi1)
nu1 := s.A / math.Sqrt(1-e2*sinPhi1*sinPhi1)
rho1 := s.A * (1 - e2) / math.Pow(1-e2*sinPhi1*sinPhi1, 1.5)
t1 := math.Tan(phi1)
t12 := t1 * t1
d := (east - cs.Eastf) / nu1
d2 := d * d
d3 := d2 * d
d4 := d2 * d2
d5 := d4 * d
phi := phi1 - (nu1*t1/rho1)*(d2/2-(1+3*t12)*d4/24)
lam := lam0 + (d-t12*d3/3+(1+3*t12)*t12*d5/15)/cosPhi1
return degree(lam), degree(phi), h
}
package crs
import (
"math"
)
type ColombiaUrban struct {
Lonf, Latf, Eastf, Northf, H0 float64
}
func (cs ColombiaUrban) String() string {
return build("colombia_urban").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
"h0", cs.H0,
).String()
}
func (cs ColombiaUrban) consts(s Spheroid) (A, B, C, D, rho0, h0a float64) {
h0a = cs.H0 / s.A
phi0 := radian(cs.Latf)
sin0 := math.Sin(phi0)
nu0 := 1 / math.Sqrt(1-s.E2()*sin0*sin0)
A = 1 + h0a/nu0
rho0 = (1 - s.E2()) / math.Pow(1-s.E2()*sin0*sin0, 1.5)
B = math.Tan(phi0) / (2 * rho0 * nu0)
C = 1 + h0a
D = rho0 * (1 + h0a/(1-s.E2()))
return A, B, C, D, rho0, h0a
}
func (cs ColombiaUrban) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
A, B, _, _, rho0, h0a := cs.consts(s)
phi := radian(lat)
dLam := radian(lon - cs.Lonf)
sinPhi := math.Sin(phi)
nu := 1 / math.Sqrt(1-s.E2()*sinPhi*sinPhi)
lamNuCos := dLam * nu * math.Cos(phi)
east := cs.Eastf + s.A*A*lamNuCos
sinPhiM := math.Sin(0.5 * (phi + radian(cs.Latf)))
rhoM := (1 - s.E2()) / math.Pow(1-s.E2()*sinPhiM*sinPhiM, 1.5)
G := 1 + h0a/rhoM
north := cs.Northf + s.A*G*rho0*((phi-radian(cs.Latf))+B*lamNuCos*lamNuCos)
return east, north, h
}
func (cs ColombiaUrban) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
_, B, C, D, _, _ := cs.consts(s)
de := (east - cs.Eastf) / s.A
dn := (north - cs.Northf) / s.A
phi := radian(cs.Latf) + dn/D - B*(de/C)*(de/C)
sinPhi := math.Sin(phi)
nu := 1 / math.Sqrt(1-s.E2()*sinPhi*sinPhi)
lon := radian(cs.Lonf) + de/(C*nu*math.Cos(phi))
return degree(lon), degree(phi), h
}
// Package crs provides coordinate reference systems, map projections, and
// datum transformations.
package crs
import (
"fmt"
)
type Func func(float64, float64, float64) (float64, float64, float64, error)
func (f Func) Round(dec int) Func {
return f.RoundAxis(dec, dec, dec)
}
func (f Func) RoundAxis(decA, decB, decC int) Func {
return func(a, b, c float64) (float64, float64, float64, error) {
a, b, c, err := f(a, b, c)
return round(a, decA), round(b, decB), round(c, decC), err
}
}
type Conversion interface {
ToGeographic(s Spheroid, a, b, c float64) (lon, lat, h float64)
FromGeographic(s Spheroid, lon, lat, h float64) (a, b, c float64)
}
type CoordinateReferenceSystem struct {
Conversion Conversion
Datum Datum
BoundingBox BoundingBox
}
func (crs CoordinateReferenceSystem) String() string {
return build("").addAll(
"conversion", crs.Conversion,
"bbox", crs.BoundingBox,
"datum", crs.Datum.Name,
).String()
}
func (crs CoordinateReferenceSystem) TransformTo(to CoordinateReferenceSystem) (Func, error) {
var err error
crs, err = crs.Intersects(to.BoundingBox)
if err != nil {
return nil, err
}
to, err = to.Intersects(crs.BoundingBox)
if err != nil {
return nil, err
}
return func(a, b, c float64) (float64, float64, float64, error) {
lon, lat, h := crs.Conversion.ToGeographic(crs.Datum.Spheroid, a, b, c)
// EPSG may publish different areas of use for geographic vs projected CRS
// that share a datum (e.g. Guam 1963 geographic vs Yap Islands). Same-datum
// hops are allowed when the point lies in either CRS extent.
inFrom := crs.BoundingBox.Contains(lon, lat)
inTo := to.BoundingBox.Contains(lon, lat)
if crs.Datum.Name == to.Datum.Name {
if !inFrom && !inTo {
return 0, 0, 0, fmt.Errorf("out of bounds: [%f,%f]", lon, lat)
}
} else if !inFrom || !inTo {
return 0, 0, 0, fmt.Errorf("out of bounds: [%f,%f]", lon, lat)
}
if crs.Datum.Name != to.Datum.Name {
var err error
lon, lat, h, err = crs.Datum.TransformTo(to.Datum, lon, lat, h)
if err != nil {
return 0, 0, 0, err
}
}
a1, b1, c1 := to.Conversion.FromGeographic(to.Datum.Spheroid, lon, lat, h)
return a1, b1, c1, nil
}, nil
}
func (crs CoordinateReferenceSystem) Intersects(bbox ...BoundingBox) (CoordinateReferenceSystem, error) {
datum, err := crs.Datum.Intersects(bbox...)
if err != nil {
return crs, err
}
crs.Datum = datum
return crs, nil
}
// AtEpoch evaluates time-dependent datum operations at the given coordinate epoch
// (decimal year), mirroring Intersects as a pre-Transform datum rewrite.
// Time-specific Helmerts are only kept when epoch matches their transformation epoch.
func (crs CoordinateReferenceSystem) AtEpoch(epoch float64) CoordinateReferenceSystem {
crs.Datum = crs.Datum.AtEpoch(epoch)
return crs
}
type CRS interface {
CoordinateReferenceSystem | int | string
}
func Load[C CRS](crs C, intersects ...BoundingBox) (CoordinateReferenceSystem, error) {
switch c := any(crs).(type) {
case int:
out, err := loadEPSG(c)
if err != nil {
return CoordinateReferenceSystem{}, err
}
return out.Intersects(intersects...)
case string:
out, err := parseCoordinateReferenceSystem(c)
if err != nil {
return CoordinateReferenceSystem{}, err
}
return out.Intersects(intersects...)
case CoordinateReferenceSystem:
return c.Intersects(intersects...)
}
return CoordinateReferenceSystem{}, fmt.Errorf("invalid crs")
}
func Transform[F, T CRS](from F, to T, intersects ...BoundingBox) (Func, error) {
fromCRS, err := Load(from, intersects...)
if err != nil {
return nil, err
}
toCRS, err := Load(to, intersects...)
if err != nil {
return nil, err
}
return fromCRS.TransformTo(toCRS)
}
// TransformAt is Transform with both datums evaluated at epoch (decimal year).
func TransformAt[F, T CRS](from F, to T, epoch float64, intersects ...BoundingBox) (Func, error) {
fromCRS, err := Load(from, intersects...)
if err != nil {
return nil, err
}
toCRS, err := Load(to, intersects...)
if err != nil {
return nil, err
}
return fromCRS.AtEpoch(epoch).TransformTo(toCRS.AtEpoch(epoch))
}
package crs
import (
"embed"
"errors"
"fmt"
"io"
"slices"
"strconv"
"strings"
"sync"
)
var (
WGS84 = Datum{
Name: "wgs84",
Spheroid: Spheroid{
Name: "wgs84",
A: 6378137,
Fi: 298.257223563,
},
}
)
//go:embed datum/*.txt
var datumDir embed.FS
var datumStore sync.Map
func RegisterDatum(d Datum) {
datumStore.Store(d.Name, d)
}
func loadDatum(name string) (d Datum, err error) {
name = strings.ToLower(name)
v, ok := datumStore.Load(name)
if ok {
return v.(Datum), nil
}
switch name {
case "wgs84":
return WGS84, nil
}
defer func() {
if err == nil {
datumStore.Store(name, d)
}
}()
file, err := datumDir.Open(fmt.Sprintf("datum/%s.txt", name))
if err != nil {
return d, fmt.Errorf("datum not found: %s", name)
}
data, err := io.ReadAll(file)
if err != nil {
return d, err
}
d, err = parseDatum(string(data))
if err != nil {
return d, err
}
d.Name = name
return d, nil
}
func stripComment(line string) string {
if i := strings.Index(line, "#"); i >= 0 {
line = line[:i]
}
return strings.TrimSpace(line)
}
func asParts(txt string) parts {
txt = flattenDSL(txt)
var pp []part
for f := range strings.FieldsSeq(txt) {
key, value, ok := strings.Cut(f, "=")
if ok {
pp = append(pp, part{
key: strings.Trim(key, "-+"),
value: value,
})
continue
}
key, value, ok = strings.Cut(f, ":")
if ok {
pp = append(pp, part{
key: strings.Trim(key, "-+"),
value: value,
})
continue
}
pp = append(pp, part{
key: strings.Trim(key, "-+"),
})
}
return pp
}
type parts []part
func (pp parts) getPart(find ...string) part {
for _, part := range pp {
if slices.Contains(find, part.key) {
return part
}
}
return part{}
}
func (pp parts) asTransformation() (Transformation, error) {
var t Transformation
t.Accuracy, _ = pp.getPart("accuracy").asFloat()
bbox, ok := pp.getPart("bbox").asBoundingBox()
if ok {
t.BoundingBox = bbox
} else {
t.BoundingBox = World
}
targetName, ok := pp.getPart("target").asString()
if !ok {
targetName = "wgs84"
}
target, err := loadDatum(targetName)
if err != nil {
if _, serr := loadSpheroid(targetName); serr == nil {
return Transformation{}, fmt.Errorf("target %q is a spheroid, not a datum", targetName)
}
return Transformation{}, err
}
if target.Name == "" {
return Transformation{}, UnsupportedError{
Err: errors.New("no target"),
}
}
t.Target = &target
op, err := pp.asOperation()
if err != nil {
if errors.As(err, &UnsupportedError{}) {
return t, nil
}
return Transformation{}, err
}
t.Operation = op
return t, nil
}
type part struct {
key string
value string
err error
}
func (p part) asString() (string, bool) {
if p.err != nil || p.value == "" {
return "", false
}
return p.value, true
}
func (p part) asFloat() (float64, bool) {
if p.err != nil || p.value == "" {
return 0, false
}
v, err := strconv.ParseFloat(p.value, 64)
if err != nil {
return 0, false
}
return v, true
}
func (pp parts) float(key string) (float64, bool) {
return pp.getPart(key).asFloat()
}
func (pp parts) floatOrZero(key string) float64 {
v, _ := pp.float(key)
return v
}
func (pp parts) hasInverse() bool {
for _, p := range pp {
if p.key == "inverse" {
return true
}
}
return false
}
func (pp parts) asOperation() (Operation, error) {
opName, ok := pp.getPart("operation").asString()
if !ok {
return nil, nil
}
var op Operation
switch opName {
case "position_vector":
op = PositionVector{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: pp.floatOrZero("rx"), Ry: pp.floatOrZero("ry"), Rz: pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
}
case "coordinate_frame":
op = PositionVector{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: -pp.floatOrZero("rx"), Ry: -pp.floatOrZero("ry"), Rz: -pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
}
case "horizontal_grid":
grid, ok := pp.getPart("grid").asString()
if !ok {
return nil, fmt.Errorf("horizontal_grid missing grid=")
}
op = HorizontalGrid(grid)
case "vertical_grid":
grid, ok := pp.getPart("grid").asString()
if !ok {
return nil, fmt.Errorf("vertical_grid missing grid=")
}
op = VerticalGrid(grid)
case "vertical_offset":
op = VerticalOffset{Dh: pp.floatOrZero("dh")}
case "vertical_offset_and_slope":
op = VerticalOffsetAndSlope{
Lat0: pp.floatOrZero("lat0"),
Lon0: pp.floatOrZero("lon0"),
Dh: pp.floatOrZero("dh"),
SlopeLat: pp.floatOrZero("slope_lat"),
SlopeLon: pp.floatOrZero("slope_lon"),
}
case "geocentric_translations":
op = PositionVector{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
}
case "identity":
op = Identity{}
case "longitude_rotation":
op = LongitudeRotation{Lon: pp.floatOrZero("lon")}
case "geographic_offset":
op = GeographicOffset{
Lat: pp.floatOrZero("lat"),
Lon: pp.floatOrZero("lon"),
}
case "geographic_3d_offset":
op = Geographic3dOffset{
Lat: pp.floatOrZero("lat"),
Lon: pp.floatOrZero("lon"),
H: pp.floatOrZero("h"),
}
case "molodensky_badekas":
op = MolodenskyBadekas{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: -pp.floatOrZero("rx"), Ry: -pp.floatOrZero("ry"), Rz: -pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
Px: pp.floatOrZero("px"), Py: pp.floatOrZero("py"), Pz: pp.floatOrZero("pz"),
}
case "molodensky_badekas_pv", "molodensky_badekas_pv_geocentric":
op = MolodenskyBadekas{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: pp.floatOrZero("rx"), Ry: pp.floatOrZero("ry"), Rz: pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
Px: pp.floatOrZero("px"), Py: pp.floatOrZero("py"), Pz: pp.floatOrZero("pz"),
}
case "time_specific_position_vector":
op = TimeSpecificPositionVector{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: pp.floatOrZero("rx"), Ry: pp.floatOrZero("ry"), Rz: pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
TransformationEpoch: pp.floatOrZero("epoch"),
}
case "time_specific_coordinate_frame":
op = TimeSpecificPositionVector{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: -pp.floatOrZero("rx"), Ry: -pp.floatOrZero("ry"), Rz: -pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
TransformationEpoch: pp.floatOrZero("epoch"),
}
case "time_dependent_position_vector":
op = TimeDependentPositionVector{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: pp.floatOrZero("rx"), Ry: pp.floatOrZero("ry"), Rz: pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
TxRate: pp.floatOrZero("dtx"),
TyRate: pp.floatOrZero("dty"),
TzRate: pp.floatOrZero("dtz"),
RxRate: pp.floatOrZero("drx"),
RyRate: pp.floatOrZero("dry"),
RzRate: pp.floatOrZero("drz"),
DsRate: pp.floatOrZero("dds"),
ReferenceEpoch: pp.floatOrZero("epoch"),
}
case "time_dependent_coordinate_frame":
op = TimeDependentPositionVector{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: -pp.floatOrZero("rx"), Ry: -pp.floatOrZero("ry"), Rz: -pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
TxRate: pp.floatOrZero("dtx"),
TyRate: pp.floatOrZero("dty"),
TzRate: pp.floatOrZero("dtz"),
RxRate: -pp.floatOrZero("drx"),
RyRate: -pp.floatOrZero("dry"),
RzRate: -pp.floatOrZero("drz"),
DsRate: pp.floatOrZero("dds"),
ReferenceEpoch: pp.floatOrZero("epoch"),
}
case "coordinate_frame_full_matrix":
op = CoordinateFrameFullMatrix{
Tx: pp.floatOrZero("tx"), Ty: pp.floatOrZero("ty"), Tz: pp.floatOrZero("tz"),
Rx: -pp.floatOrZero("rx"), Ry: -pp.floatOrZero("ry"), Rz: -pp.floatOrZero("rz"),
Ds: pp.floatOrZero("ds"),
}
case "velocity_grid":
grid, ok := pp.getPart("grid").asString()
if !ok {
return nil, fmt.Errorf("velocity_grid missing grid=")
}
op = VelocityGrid{
Grid: grid,
Dt: pp.floatOrZero("dt"),
Epoch: pp.floatOrZero("epoch"),
}
default:
return nil, UnsupportedError{
Err: fmt.Errorf("invalid operation: %s", opName),
}
}
if pp.hasInverse() {
return Inverse{
Operation: op,
}, nil
}
return op, nil
}
func (p part) asBoundingBox() (BoundingBox, bool) {
var (
b BoundingBox
err error
)
if p.err != nil || p.value == "" {
return b, false
}
for i, v := range strings.Split(p.value, ",") {
if p.value == "" {
continue
}
switch i {
case 0:
b.MinLon, err = strconv.ParseFloat(v, 64)
if err != nil {
return b, false
}
case 1:
b.MinLat, err = strconv.ParseFloat(v, 64)
if err != nil {
return b, false
}
case 2:
b.MaxLon, err = strconv.ParseFloat(v, 64)
if err != nil {
return b, false
}
case 3:
b.MaxLat, err = strconv.ParseFloat(v, 64)
if err != nil {
return b, false
}
}
}
return b, true
}
type Datum struct {
Name string
Spheroid Spheroid
Transformations []Transformation
// coordinateEpoch is set by AtEpoch / TransformAt. When hasCoordinateEpoch
// is false (plain Transform), time-specific Helmerts are excluded from paths.
coordinateEpoch float64
hasCoordinateEpoch bool
}
func (d Datum) String() string {
return build("").add("spheroid", d.Spheroid.Name).add("", d.Transformations).String()
}
func (d Datum) Intersects(bbox ...BoundingBox) (Datum, error) {
if len(bbox) == 0 {
return d, nil
}
if len(d.Transformations) == 0 {
return d, nil
}
// Clone so DeleteFunc cannot mutate the cached LoadDatum slice backing array.
t := slices.DeleteFunc(slices.Clone(d.Transformations), func(t Transformation) bool {
for _, b := range bbox {
if !t.BoundingBox.Intersects(b) {
return true
}
}
return false
})
if len(t) == 0 {
// No published hub ops cover this area (e.g. Guam1963 on Yap). Keep the
// spheroid so same-datum projection still works; cross-datum will fail later.
return Datum{
Name: d.Name,
Spheroid: d.Spheroid,
coordinateEpoch: d.coordinateEpoch,
hasCoordinateEpoch: d.hasCoordinateEpoch,
}, nil
}
return Datum{
Name: d.Name,
Spheroid: d.Spheroid,
Transformations: t,
coordinateEpoch: d.coordinateEpoch,
hasCoordinateEpoch: d.hasCoordinateEpoch,
}, nil
}
// AtEpoch returns a copy of d with time-dependent operations evaluated at epoch
// (decimal year). TimeDependentPositionVector becomes PositionVector; VelocityGrid
// with a non-zero Epoch gets Dt = epoch - Epoch. Time-specific Helmerts that do
// not match epoch are dropped from the in-memory list (pathfinding also gates them).
func (d Datum) AtEpoch(epoch float64) Datum {
out := Datum{
Name: d.Name,
Spheroid: d.Spheroid,
coordinateEpoch: epoch,
hasCoordinateEpoch: true,
}
if len(d.Transformations) == 0 {
return out
}
ts := make([]Transformation, 0, len(d.Transformations))
for _, t := range d.Transformations {
op := operationAtEpoch(t.Operation, epoch)
if !timeSpecificAllowed(op, true, epoch) {
continue
}
t.Operation = op
ts = append(ts, t)
}
out.Transformations = ts
return out
}
func operationAtEpoch(op Operation, epoch float64) Operation {
switch o := op.(type) {
case TimeDependentPositionVector:
return o.at(epoch)
case Inverse:
return Inverse{Operation: operationAtEpoch(o.Operation, epoch)}
case VelocityGrid:
if o.Epoch != 0 {
o.Dt = epoch - o.Epoch
}
return o
default:
return op
}
}
// primeMeridianLongitude is degrees east of Greenwich for this datum's prime
// meridian, taken from a longitude_rotation hop (0 when the datum is Greenwich-based).
// EPSG areas of use are Greenwich-oriented; use lon+primeMeridianLongitude() for Contains.
func (d Datum) primeMeridianLongitude() float64 {
for _, t := range d.Transformations {
if lr, ok := t.Operation.(LongitudeRotation); ok {
return lr.Lon
}
}
return 0
}
func (d Datum) ToWGS84(lon, lat, h float64) (float64, float64, float64, error) {
return d.toWGS84Visited(lon, lat, h, nil)
}
func (d Datum) toWGS84Visited(lon, lat, h float64, visited map[string]bool) (float64, float64, float64, error) {
if len(d.Transformations) == 0 {
return lon, lat, h, nil
}
if visited == nil {
visited = make(map[string]bool)
}
if visited[d.Name] {
return 0, 0, 0, UnsupportedError{
Err: errors.New("datum transformation cycle"),
}
}
visited[d.Name] = true
defer delete(visited, d.Name)
glon := lon + d.primeMeridianLongitude()
var lastErr error
for _, t := range d.Transformations {
if !timeSpecificAllowed(t.Operation, d.hasCoordinateEpoch, d.coordinateEpoch) {
continue
}
lon0, lat0, h0, err := t.toWGS84Visited(d.Spheroid, lon, lat, h, visited, glon)
if err != nil {
lastErr = err
continue
}
return lon0, lat0, h0, nil
}
if lastErr != nil {
return 0, 0, 0, lastErr
}
return 0, 0, 0, UnsupportedError{
Err: errors.New("no valid transformation found"),
}
}
func (d Datum) FromWGS84(lon0, lat0, h0 float64) (float64, float64, float64, error) {
return d.fromWGS84Visited(lon0, lat0, h0, nil)
}
func (d Datum) fromWGS84Visited(lon0, lat0, h0 float64, visited map[string]bool) (float64, float64, float64, error) {
if len(d.Transformations) == 0 {
return lon0, lat0, h0, nil
}
if visited == nil {
visited = make(map[string]bool)
}
if visited[d.Name] {
return 0, 0, 0, UnsupportedError{
Err: errors.New("datum transformation cycle"),
}
}
visited[d.Name] = true
defer delete(visited, d.Name)
var lastErr error
for _, t := range d.Transformations {
if !timeSpecificAllowed(t.Operation, d.hasCoordinateEpoch, d.coordinateEpoch) {
continue
}
lon, lat, h, err := t.fromWGS84Visited(d.Spheroid, lon0, lat0, h0, visited)
if err != nil {
lastErr = err
continue
}
return lon, lat, h, nil
}
if lastErr != nil {
return 0, 0, 0, lastErr
}
return 0, 0, 0, UnsupportedError{
Err: errors.New("no valid transformation found"),
}
}
func (d Datum) TransformTo(target Datum, lon, lat, h float64) (float64, float64, float64, error) {
if d.Name == target.Name {
return lon, lat, h, nil
}
hasEpoch := d.hasCoordinateEpoch || target.hasCoordinateEpoch
epoch := d.coordinateEpoch
if target.hasCoordinateEpoch {
epoch = target.coordinateEpoch
}
excluded := make(map[edgeKey]bool)
var lastErr error
for {
path, err := findBestPath(d, target, lon, lat, excluded, hasEpoch, epoch)
if err != nil {
if lastErr != nil {
return 0, 0, 0, lastErr
}
return 0, 0, 0, err
}
outLon, outLat, outH, failed, err := applyPath(path, lon, lat, h, excluded, hasEpoch, epoch)
if err == nil {
return outLon, outLat, outH, nil
}
lastErr = err
if len(failed) == 0 {
return 0, 0, 0, lastErr
}
grew := false
for _, k := range failed {
if !excluded[k] {
excluded[k] = true
grew = true
}
}
if !grew {
return 0, 0, 0, lastErr
}
}
}
package crs
import (
"container/heap"
"errors"
"fmt"
"io/fs"
"strings"
"sync"
)
const maxPathHops = 8
// isHubEquivalent reports whether d may be treated as a free identity hop to/from
// WGS84. Datums with no exported ops (typically PROJ-ballpark-only in the wild)
// stay hub-equivalent so CRS remain reachable without exporting ballpark edges.
// Published null hubs with accuracy/bbox use operation=identity instead of empty.
func isHubEquivalent(d Datum) bool {
return len(d.Transformations) == 0
}
type revEdge struct {
owner string
index int
}
type edgeKey struct {
owner string
index int
inverse bool
}
type graphEdge struct {
key edgeKey
to string
accuracy float64
inverse bool
}
var (
reverseOnce sync.Once
reverseIndex map[string][]revEdge
)
func ensureReverseIndex() {
reverseOnce.Do(func() {
reverseIndex = make(map[string][]revEdge)
entries, err := fs.ReadDir(datumDir, "datum")
if err != nil {
return
}
for _, e := range entries {
if e.IsDir() || !strings.HasSuffix(e.Name(), ".txt") {
continue
}
name := strings.TrimSuffix(e.Name(), ".txt")
d, err := loadDatum(name)
if err != nil {
continue
}
for i, t := range d.Transformations {
if t.Target == nil || t.Operation == nil {
continue
}
tgt := strings.ToLower(t.Target.Name)
reverseIndex[tgt] = append(reverseIndex[tgt], revEdge{owner: name, index: i})
}
}
})
}
func datumName(d Datum) string {
return strings.ToLower(d.Name)
}
// outgoingEdges returns graph edges from node that cover (lon, lat), excluding keys in excluded.
// Forward and inverse candidates are already ordered by the owner's transformation sort
// (accuracy, then area); we emit one best edge per neighbor for Dijkstra, but
// edgesBetween returns all candidates for apply.
// Time-specific Helmerts are included only when hasCoordEpoch and coordEpoch match.
func outgoingEdges(node string, lon, lat float64, excluded map[edgeKey]bool, hasCoordEpoch bool, coordEpoch float64) []graphEdge {
ensureReverseIndex()
d, err := loadDatum(node)
if err != nil {
return nil
}
// AoU bboxes are Greenwich-based; lon may be relative to this node's PM.
glon := lon + d.primeMeridianLongitude()
best := make(map[string]graphEdge) // neighbor -> best edge
consider := func(e graphEdge) {
if excluded[e.key] {
return
}
cur, ok := best[e.to]
if !ok || e.accuracy < cur.accuracy {
best[e.to] = e
}
}
for i, t := range d.Transformations {
if t.Target == nil || t.Operation == nil {
continue
}
if !timeSpecificAllowed(t.Operation, hasCoordEpoch, coordEpoch) {
continue
}
if !t.BoundingBox.Contains(glon, lat) {
continue
}
tgt := strings.ToLower(t.Target.Name)
if tgt == "" || tgt == node {
continue
}
consider(graphEdge{
key: edgeKey{owner: node, index: i, inverse: false},
to: tgt,
accuracy: t.Accuracy,
inverse: false,
})
}
for _, rev := range reverseIndex[node] {
owner, err := loadDatum(rev.owner)
if err != nil || rev.index < 0 || rev.index >= len(owner.Transformations) {
continue
}
t := owner.Transformations[rev.index]
if t.Target == nil || t.Operation == nil {
continue
}
if !timeSpecificAllowed(t.Operation, hasCoordEpoch, coordEpoch) {
continue
}
if !t.BoundingBox.Contains(glon, lat) {
continue
}
consider(graphEdge{
key: edgeKey{owner: rev.owner, index: rev.index, inverse: true},
to: rev.owner,
accuracy: t.Accuracy,
inverse: true,
})
}
out := make([]graphEdge, 0, len(best))
for _, e := range best {
out = append(out, e)
}
return out
}
// edgesBetween returns all covering edges from→to (forward and inverse), in
// owner transformation order (accuracy ascending within each owner list).
func edgesBetween(from, to string, lon, lat float64, excluded map[edgeKey]bool, hasCoordEpoch bool, coordEpoch float64) []graphEdge {
ensureReverseIndex()
from = strings.ToLower(from)
to = strings.ToLower(to)
var out []graphEdge
glon := lon
if fromDatum, err := loadDatum(from); err == nil {
glon = lon + fromDatum.primeMeridianLongitude()
for i, t := range fromDatum.Transformations {
if t.Target == nil || t.Operation == nil {
continue
}
if !timeSpecificAllowed(t.Operation, hasCoordEpoch, coordEpoch) {
continue
}
if !strings.EqualFold(t.Target.Name, to) {
continue
}
if !t.BoundingBox.Contains(glon, lat) {
continue
}
k := edgeKey{owner: from, index: i, inverse: false}
if excluded[k] {
continue
}
out = append(out, graphEdge{key: k, to: to, accuracy: t.Accuracy, inverse: false})
}
}
if toDatum, err := loadDatum(to); err == nil {
// Inverse edge: coordinates remain in `from`'s frame (Greenwich AoU check).
for i, t := range toDatum.Transformations {
if t.Target == nil || t.Operation == nil {
continue
}
if !timeSpecificAllowed(t.Operation, hasCoordEpoch, coordEpoch) {
continue
}
if !strings.EqualFold(t.Target.Name, from) {
continue
}
if !t.BoundingBox.Contains(glon, lat) {
continue
}
k := edgeKey{owner: to, index: i, inverse: true}
if excluded[k] {
continue
}
out = append(out, graphEdge{key: k, to: to, accuracy: t.Accuracy, inverse: true})
}
}
return out
}
type pathNode struct {
name string
cost float64
hops int
viaHub bool // true if wgs84 appears as a non-terminal intermediate so far
prev int // index in settled slice; -1 for start
edge graphEdge
}
type pqItem struct {
idx int
cost float64
hops int
hub bool
}
type pathPQ []pqItem
func (p pathPQ) Len() int {
return len(p)
}
func (p pathPQ) Less(i, j int) bool {
if p[i].cost != p[j].cost {
return p[i].cost < p[j].cost
}
if p[i].hops != p[j].hops {
return p[i].hops < p[j].hops
}
if p[i].hub != p[j].hub {
return !p[i].hub && p[j].hub
}
return p[i].idx < p[j].idx
}
func (p pathPQ) Swap(i, j int) {
p[i], p[j] = p[j], p[i]
}
func (p *pathPQ) Push(x any) {
*p = append(*p, x.(pqItem))
}
func (p *pathPQ) Pop() any {
old := *p
n := len(old)
item := old[n-1]
*p = old[:n-1]
return item
}
func betterPath(aCost float64, aHops int, aHub bool, bCost float64, bHops int, bHub bool) bool {
if aCost != bCost {
return aCost < bCost
}
if aHops != bHops {
return aHops < bHops
}
if aHub != bHub {
return !aHub && bHub
}
return false
}
// findBestPath finds a minimum accumulated-accuracy path from→to covering (lon,lat).
func findBestPath(from, to Datum, lon, lat float64, excluded map[edgeKey]bool, hasCoordEpoch bool, coordEpoch float64) ([]string, error) {
ensureReverseIndex()
start := datumName(from)
goal := datumName(to)
if start == "" || goal == "" {
return nil, UnsupportedError{Err: errors.New("missing datum name")}
}
if start == goal {
return []string{start}, nil
}
if excluded == nil {
excluded = map[edgeKey]bool{}
}
// Datums with no published ops are treated as WGS84-equivalent (identity hops).
searchFrom, searchTo := from, to
prefix, suffix := "", ""
if start != "wgs84" && isHubEquivalent(from) {
prefix = start
searchFrom = WGS84
}
if goal != "wgs84" && isHubEquivalent(to) {
suffix = goal
searchTo = WGS84
}
coreStart := datumName(searchFrom)
coreGoal := datumName(searchTo)
var core []string
if coreStart == coreGoal {
core = []string{coreStart}
} else {
var err error
core, err = findBestPathCore(coreStart, coreGoal, lon, lat, excluded, hasCoordEpoch, coordEpoch)
if err != nil {
return nil, err
}
}
path := make([]string, 0, len(core)+2)
if prefix != "" {
path = append(path, prefix)
}
path = append(path, core...)
if suffix != "" {
path = append(path, suffix)
}
return dedupeConsecutive(path), nil
}
func dedupeConsecutive(path []string) []string {
if len(path) == 0 {
return path
}
out := []string{path[0]}
for _, n := range path[1:] {
if n != out[len(out)-1] {
out = append(out, n)
}
}
return out
}
func findBestPathCore(start, goal string, lon, lat float64, excluded map[edgeKey]bool, hasCoordEpoch bool, coordEpoch float64) ([]string, error) {
nodes := []pathNode{{
name: start,
cost: 0,
hops: 0,
prev: -1,
}}
best := map[string]int{start: 0}
pq := &pathPQ{{idx: 0, cost: 0, hops: 0, hub: false}}
heap.Init(pq)
for pq.Len() > 0 {
item := heap.Pop(pq).(pqItem)
cur := nodes[item.idx]
if bestIdx, ok := best[cur.name]; ok && bestIdx != item.idx {
continue
}
if cur.name == goal {
return reconstructPath(nodes, item.idx), nil
}
if cur.hops >= maxPathHops {
continue
}
for _, e := range outgoingEdges(cur.name, lon, lat, excluded, hasCoordEpoch, coordEpoch) {
nextCost := cur.cost + e.accuracy
nextHops := cur.hops + 1
viaHub := cur.viaHub
if cur.name == "wgs84" && start != "wgs84" {
viaHub = true
}
if prevIdx, ok := best[e.to]; ok {
prev := nodes[prevIdx]
if !betterPath(nextCost, nextHops, viaHub, prev.cost, prev.hops, prev.viaHub) {
continue
}
}
idx := len(nodes)
nodes = append(nodes, pathNode{
name: e.to,
cost: nextCost,
hops: nextHops,
viaHub: viaHub,
prev: item.idx,
edge: e,
})
best[e.to] = idx
heap.Push(pq, pqItem{idx: idx, cost: nextCost, hops: nextHops, hub: viaHub})
}
}
return nil, UnsupportedError{Err: fmt.Errorf("no transformation path from %s to %s", start, goal)}
}
func reconstructPath(nodes []pathNode, idx int) []string {
var rev []string
for idx >= 0 {
rev = append(rev, nodes[idx].name)
idx = nodes[idx].prev
}
for i, j := 0, len(rev)-1; i < j; i, j = i+1, j-1 {
rev[i], rev[j] = rev[j], rev[i]
}
return rev
}
func applyEdge(from Datum, e graphEdge, lon, lat, h float64) (float64, float64, float64, error) {
owner, err := loadDatum(e.key.owner)
if err != nil {
return 0, 0, 0, err
}
if e.key.index < 0 || e.key.index >= len(owner.Transformations) {
return 0, 0, 0, UnsupportedError{Err: errors.New("invalid transformation index")}
}
t := owner.Transformations[e.key.index]
if e.inverse {
return t.FromDatum(owner.Spheroid, lon, lat, h)
}
// Forward: owner must be `from`.
_ = from
return t.ToDatum(owner.Spheroid, lon, lat, h)
}
func applyHop(fromName, toName string, lon, lat, h float64, excluded map[edgeKey]bool, hasCoordEpoch bool, coordEpoch float64) (float64, float64, float64, []edgeKey, error) {
from, err := loadDatum(fromName)
if err != nil {
return 0, 0, 0, nil, err
}
toDatum, toErr := loadDatum(toName)
cands := edgesBetween(fromName, toName, lon, lat, excluded, hasCoordEpoch, coordEpoch)
if len(cands) == 0 {
// Identity hop between WGS84 and a hub-equivalent datum (no published ops).
if toErr == nil && identityHopAllowed(from, toDatum, fromName, toName) {
return lon, lat, h, nil, nil
}
return 0, 0, 0, nil, UnsupportedError{
Err: fmt.Errorf("no transformation from %s to %s", fromName, toName),
}
}
var lastErr error
tried := make([]edgeKey, 0, len(cands))
for _, e := range cands {
tried = append(tried, e.key)
outLon, outLat, outH, err := applyEdge(from, e, lon, lat, h)
if err != nil {
lastErr = err
continue
}
return outLon, outLat, outH, tried, nil
}
if lastErr == nil {
lastErr = UnsupportedError{Err: errors.New("no valid transformation found")}
}
return 0, 0, 0, tried, lastErr
}
func identityHopAllowed(from, to Datum, fromName, toName string) bool {
fromEq := fromName == "wgs84" || isHubEquivalent(from)
toEq := toName == "wgs84" || isHubEquivalent(to)
return fromEq && toEq
}
func applyPath(nodes []string, lon, lat, h float64, excluded map[edgeKey]bool, hasCoordEpoch bool, coordEpoch float64) (float64, float64, float64, []edgeKey, error) {
if len(nodes) < 2 {
return lon, lat, h, nil, nil
}
var failed []edgeKey
for i := 0; i < len(nodes)-1; i++ {
var tried []edgeKey
var err error
lon, lat, h, tried, err = applyHop(nodes[i], nodes[i+1], lon, lat, h, excluded, hasCoordEpoch, coordEpoch)
failed = append(failed, tried...)
if err != nil {
return 0, 0, 0, failed, err
}
}
return lon, lat, h, nil, nil
}
package crs
import (
"embed"
"fmt"
"io"
"strconv"
"strings"
"sync"
)
//go:embed epsg/*.txt
var epsgDir embed.FS
var epsgStore sync.Map
func RegisterEPSG(code int, c CoordinateReferenceSystem) {
epsgStore.Store(code, c)
}
func loadEPSG(code int) (c CoordinateReferenceSystem, err error) {
v, ok := epsgStore.Load(code)
if ok {
return v.(CoordinateReferenceSystem), nil
}
defer func() {
if err == nil {
epsgStore.Store(code, c)
}
}()
file, err := epsgDir.Open(fmt.Sprintf("epsg/%d.txt", code))
if err != nil {
return c, fmt.Errorf("epsg not found: %d", code)
}
data, err := io.ReadAll(file)
if err != nil {
return c, err
}
return parseCoordinateReferenceSystem(string(data))
}
func parseBoundingBox(s string) (BoundingBox, bool) {
parts := strings.Split(s, ",")
if len(parts) != 4 {
return BoundingBox{}, false
}
vals := make([]float64, 4)
for i, p := range parts {
v, err := strconv.ParseFloat(strings.TrimSpace(p), 64)
if err != nil {
return BoundingBox{}, false
}
vals[i] = v
}
return BoundingBox{
MinLon: vals[0],
MinLat: vals[1],
MaxLon: vals[2],
MaxLat: vals[3],
}, true
}
type epsgParams struct {
str map[string]string
f map[string]float64
}
func newEPSGParams(fields map[string]string) epsgParams {
p := epsgParams{
str: fields,
f: make(map[string]float64),
}
for k, v := range fields {
if fv, err := strconv.ParseFloat(v, 64); err == nil {
p.f[k] = fv
}
}
return p
}
func (p epsgParams) ifStr(str string, then, els float64) float64 {
if _, ok := p.str[str]; ok {
return then
}
return els
}
func (p epsgParams) float(keys ...string) float64 {
return p.floatOr(0, keys...)
}
func (p epsgParams) floatOr(def float64, keys ...string) float64 {
for _, k := range keys {
if v, ok := p.f[k]; ok {
return v
}
}
return def
}
func buildConversion(name string, fields map[string]string) (Conversion, error) {
p := newEPSGParams(fields)
switch name {
case "geographic":
return Geographic{}, nil
case "geocentric":
return Geocentric{}, nil
case "transverse_mercator":
return TransverseMercator{
Lonf: p.float("lonf"), Latf: p.float("latf"), Scale: p.floatOr(1, "scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
ZoneWidth: p.float("zone_width"),
}, nil
case "utm":
return TransverseMercator{
Lonf: p.float("zone")*6 - 183, Scale: 0.9996,
Eastf: 500000, Northf: p.ifStr("southern", 1e7, 0),
}, nil
case "web_mercator":
return WebMercator{
Lonf: p.float("lonf"), Latf: p.float("latf"), Scale: p.float("scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_conformal_conic":
return LambertConformalConic{
Lonf: p.float("lonf"), Latf: p.float("latf"), Scale: p.floatOr(1, "scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_conformal_conic_1sp_variant_b":
return LambertConformalConic1SPVariantB{
Lonf: p.float("lonf"), Lat0: p.float("lat0"), Latf: p.float("latf"), Scale: p.floatOr(1, "scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_conformal_conic_2sp":
return LambertConformalConic2SP{
Lonf: p.float("lonf"), Latf: p.float("latf"), Sp1: p.float("sp1"), Sp2: p.float("sp2"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_conformal_conic_2sp_michigan":
return LambertConformalConic2SPMichigan{
Lonf: p.float("lonf"), Latf: p.float("latf"), Sp1: p.float("sp1"), Sp2: p.float("sp2"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_conformal_conic_2sp_belgium":
return LambertConformalConic2SPBelgium{
Lonf: p.float("lonf"), Latf: p.float("latf"), Sp1: p.float("sp1"), Sp2: p.float("sp2"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_conic_near_conformal":
return LambertConicNearConformal{
Lonf: p.float("lonf"), Latf: p.float("latf"), Scale: p.floatOr(1, "scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_azimuthal_equal_area":
return LambertAzimuthalEqualArea{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_azimuthal_equal_area_spherical":
return LambertAzimuthalEqualAreaSpherical{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_cylindrical_equal_area":
return LambertCylindricalEqualArea{
Lonf: p.float("lonf"), Sp1: p.float("sp1", "latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "lambert_cylindrical_equal_area_spherical":
return LambertCylindricalEqualAreaSpherical{
Lonf: p.float("lonf"), Sp1: p.float("sp1", "latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "albers_equal_area":
return AlbersEqualArea{
Lonf: p.float("lonf"), Latf: p.float("latf"), Sp1: p.float("sp1"), Sp2: p.float("sp2"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "krovak":
return Krovak{
Lonf: p.float("lonf"), Latf: p.float("latf"), Alpha: p.float("alpha"), Scale: p.float("scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "krovak_modified", "krovak_modified_north_orientated": // alias: always east,north
return KrovakModified{
Lonf: p.float("lonf"), Latf: p.float("latf"), Alpha: p.float("alpha"), Scale: p.float("scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "hotine_oblique_mercator_a":
return HotineObliqueMercatorA{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Alpha: p.float("azimuth", "alpha"), Gamma: p.float("gamma"),
Scale: p.float("scale"), Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "swiss_oblique_mercator":
return SwissObliqueMercator{
Lonf: p.float("lonf"), Latf: p.float("latf"), Alpha: p.float("azimuth", "alpha"), Gamma: p.float("gamma"),
Scale: p.float("scale"), Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "laborde_oblique_mercator":
return LabordeObliqueMercator{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Alpha: p.float("azimuth", "alpha"),
Scale: p.float("scale"), Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "cassini_soldner":
return CassiniSoldner{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "hyperbolic_cassini_soldner":
return HyperbolicCassiniSoldner{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "polar_stereographic_a":
return PolarStereographicA{
Lonf: p.float("lonf"), Latf: p.float("latf"), Scale: p.float("scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "polar_stereographic_b":
return PolarStereographicB{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "polar_stereographic_c":
return PolarStereographicC{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "oblique_stereographic":
return ObliqueStereographic{
Lonf: p.float("lonf"), Latf: p.float("latf"), Scale: p.float("scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "mercator_a":
return MercatorA{
Lonf: p.float("lonf"), Latf: p.float("latf"), Scale: p.float("scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "mercator_b":
return MercatorB{
Lonf: p.float("lonf"), Sp1: p.float("sp1"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "azimuthal_equidistant":
return AzimuthalEquidistant{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "modified_azimuthal_equidistant":
return ModifiedAzimuthalEquidistant{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "american_polyconic":
return AmericanPolyconic{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "bonne_south_orientated":
return BonneSouthOrientated{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "equal_earth":
return EqualEarth{
Lonf: p.float("lonf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "equidistant_cylindrical":
return EquidistantCylindrical{
Lonf: p.float("lonf"), Latf: p.float("latf"), Sp1: p.float("sp1", "latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "colombia_urban":
return ColombiaUrban{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"), H0: p.float("h0"),
}, nil
case "guam_projection":
return GuamProjection{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "local_orthographic":
return LocalOrthographic{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Azimuth: p.float("azimuth", "alpha"), Scale: p.floatOr(1, "scale"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "tunisia_mining_grid":
return TunisiaMiningGrid{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
case "new_zealand_map_grid":
return NewZealandMapGrid{
Lonf: p.float("lonf"), Latf: p.float("latf"),
Eastf: p.float("eastf"), Northf: p.float("northf"),
}, nil
}
return nil, fmt.Errorf("unknown conversion: %s", name)
}
package crs
import (
"math"
)
type EquidistantCylindrical struct {
Lonf, Latf, Sp1, Eastf, Northf float64
}
func (cs EquidistantCylindrical) String() string {
return build("equidistant_cylindrical").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"sp1", cs.Sp1,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs EquidistantCylindrical) rc(s Spheroid) float64 {
phi1 := radian(cs.Sp1)
sin1 := math.Sin(phi1)
nu1 := 1 / math.Sqrt(1-s.E2()*sin1*sin1)
return nu1 * math.Cos(phi1)
}
func (cs EquidistantCylindrical) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
rc := cs.rc(s)
m0 := meridianDistance(s, radian(cs.Latf))
east := cs.Eastf + s.A*rc*(radian(lon)-radian(cs.Lonf))
north := cs.Northf + meridianDistance(s, radian(lat)) - m0
return east, north, h
}
func (cs EquidistantCylindrical) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
rc := cs.rc(s)
m0 := meridianDistance(s, radian(cs.Latf))
lam := radian(cs.Lonf) + (east-cs.Eastf)/(s.A*rc)
phi := footpointLatitude(s, north-cs.Northf+m0)
return degree(lam), degree(phi), h
}
package crs
import (
"math"
)
// EqualEarth is EPSG method 1078 (PROJ +proj=eqearth).
// Formulas: IOGP Guidance Note 7-2 §3.4.4.
type EqualEarth struct {
Lonf, Eastf, Northf float64
}
func (cs EqualEarth) String() string {
return build("equal_earth").addAll(
"lonf", cs.Lonf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
var equalEarthA = [...]float64{1.340264, -0.081106, 0.000893, 0.003796}
func (cs EqualEarth) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
e := s.E()
e2 := s.E2()
qp := authalicQ(1, e, e2)
rqda := math.Sqrt(0.5 * qp)
sbeta := math.Sin(radian(lat))
if e >= 1e-12 {
sbeta = authalicQ(sbeta, e, e2) / qp
sbeta = clamp(sbeta, -1, 1)
}
m := math.Sqrt(3) / 2
psi := math.Asin(m * sbeta)
psi2 := psi * psi
psi6 := psi2 * psi2 * psi2
a1, a2, a3, a4 := equalEarthA[0], equalEarthA[1], equalEarthA[2], equalEarthA[3]
lam := radian(lon - cs.Lonf)
x := lam * math.Cos(psi) / (m * (a1 + 3*a2*psi2 + psi6*(7*a3+9*a4*psi2)))
y := psi * (a1 + a2*psi2 + psi6*(a3+a4*psi2))
return cs.Eastf + s.A*rqda*x, cs.Northf + s.A*rqda*y, h
}
func (cs EqualEarth) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
e := s.E()
e2 := s.E2()
qp := authalicQ(1, e, e2)
rqda := math.Sqrt(0.5 * qp)
a1, a2, a3, a4 := equalEarthA[0], equalEarthA[1], equalEarthA[2], equalEarthA[3]
m := math.Sqrt(3) / 2
x := (east - cs.Eastf) / (s.A * rqda)
y := (north - cs.Northf) / (s.A * rqda)
yc := clamp(y, -1.3173627591574, 1.3173627591574)
for range 12 {
y2 := yc * yc
y6 := y2 * y2 * y2
f := yc*(a1+a2*y2+y6*(a3+a4*y2)) - y
fder := a1 + 3*a2*y2 + y6*(7*a3+9*a4*y2)
dy := f / fder
yc -= dy
if math.Abs(dy) < 1e-11 {
break
}
}
y2 := yc * yc
y6 := y2 * y2 * y2
lam := m * x * (a1 + 3*a2*y2 + y6*(7*a3+9*a4*y2)) / math.Cos(yc)
sbeta := math.Sin(yc) / m
beta := math.Asin(clamp(sbeta, -1, 1))
phi := beta
if e >= 1e-12 {
phi = authalicToGeodetic(beta, s)
}
return degree(lam + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"fmt"
)
type UnsupportedError struct {
Err error
}
func (e UnsupportedError) Unwrap() error {
return e.Err
}
func (e UnsupportedError) Error() string {
return fmt.Sprintf("unsupported: %s", e.Err)
}
type OutOfBoundsError struct {
Err error
}
func (e OutOfBoundsError) Unwrap() error {
return e.Err
}
func (e OutOfBoundsError) Error() string {
return fmt.Sprintf("out of bounds: %s", e.Err)
}
type GridNotFoundError struct {
Err error
}
func (e GridNotFoundError) Unwrap() error {
return e.Err
}
func (e GridNotFoundError) Error() string {
return fmt.Sprintf("grid not found: %s", e.Err)
}
package crs
type Geocentric struct{}
func (g Geocentric) String() string {
return "geocentric"
}
func (g Geocentric) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
return s.GeographicToGeocentric(lon, lat, h)
}
func (g Geocentric) ToGeographic(s Spheroid, x, y, z float64) (float64, float64, float64) {
return s.GeocentricToGeographic(x, y, z)
}
package crs
type Geographic struct{}
func (g Geographic) String() string {
return "geographic"
}
func (g Geographic) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
return lon, lat, h
}
func (g Geographic) ToGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
return lon, lat, h
}
package crs
// LongitudeRotation shifts longitude between a local prime meridian and Greenwich
// (or another Greenwich-oriented sibling). Lon is degrees east of Greenwich for
// the local meridian: ToTarget adds Lon (local → Greenwich-oriented), FromTarget subtracts.
type LongitudeRotation struct {
Lon float64
}
func (l LongitudeRotation) String() string {
return build("longitude_rotation").addAll("lon", l.Lon).String()
}
func (l LongitudeRotation) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return lon + l.Lon, lat, h, nil
}
func (l LongitudeRotation) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return lon0 - l.Lon, lat0, h0, nil
}
// Identity is a published null datum link (e.g. EPSG:1149): coordinates unchanged,
// kept so accuracy and area-of-use still participate in path selection.
type Identity struct{}
func (Identity) String() string {
return "identity"
}
func (Identity) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return lon, lat, h, nil
}
func (Identity) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return lon0, lat0, h0, nil
}
type GeographicOffset struct {
Lat float64 // arc-seconds
Lon float64 // arc-seconds
}
func (o GeographicOffset) String() string {
return build("geographic_offset").addAll(
"lat", o.Lat,
"lon", o.Lon,
).String()
}
func (o GeographicOffset) deg() (dLon, dLat float64) {
return o.Lon / 3600, o.Lat / 3600
}
func (o GeographicOffset) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
dLon, dLat := o.deg()
return lon + dLon, lat + dLat, h, nil
}
func (o GeographicOffset) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
dLon, dLat := o.deg()
return lon0 - dLon, lat0 - dLat, h0, nil
}
type Geographic3dOffset struct {
Lat float64 // arc-seconds
Lon float64 // arc-seconds
H float64 // metres
}
func (o Geographic3dOffset) String() string {
return build("geographic_3d_offset").addAll(
"lat", o.Lat,
"lon", o.Lon,
"h", o.H,
).String()
}
func (o Geographic3dOffset) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return lon + o.Lon/3600, lat + o.Lat/3600, h + o.H, nil
}
func (o Geographic3dOffset) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return lon0 - o.Lon/3600, lat0 - o.Lat/3600, h0 - o.H, nil
}
package crs
import "math"
type CoordinateFrameFullMatrix struct {
Tx, Ty, Tz, Rx, Ry, Rz, Ds float64
}
func (m CoordinateFrameFullMatrix) String() string {
return build("coordinate_frame_full_matrix").addAll(
"tx", m.Tx,
"ty", m.Ty,
"tz", m.Tz,
"rx", m.Rx,
"ry", m.Ry,
"rz", m.Rz,
"ds", m.Ds,
).String()
}
func (m CoordinateFrameFullMatrix) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
x, y, z := source.GeographicToGeocentric(lon, lat, h)
x0, y0, z0 := calcHelmertExact(x, y, z, m.Tx, m.Ty, m.Tz, m.Rx, m.Ry, m.Rz, m.Ds)
lon0, lat0, h0 := target.GeocentricToGeographic(x0, y0, z0)
return lon0, lat0, h0, nil
}
func (m CoordinateFrameFullMatrix) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
x0, y0, z0 := source.GeographicToGeocentric(lon0, lat0, h0)
x, y, z := calcHelmertExactInverse(x0, y0, z0, m.Tx, m.Ty, m.Tz, m.Rx, m.Ry, m.Rz, m.Ds)
lon, lat, h := target.GeocentricToGeographic(x, y, z)
return lon, lat, h, nil
}
const (
asec = math.Pi / 648000
ppm = 0.000001
)
func calcHelmertExact(x, y, z, tx, ty, tz, rx, ry, rz, ds float64) (x1, y1, z1 float64) {
r00, r01, r02, r10, r11, r12, r20, r21, r22 := rotationMatrixExact(rx*asec, ry*asec, rz*asec)
s := 1 + ds*ppm
x1 = s*(r00*x+r01*y+r02*z) + tx
y1 = s*(r10*x+r11*y+r12*z) + ty
z1 = s*(r20*x+r21*y+r22*z) + tz
return x1, y1, z1
}
func calcHelmertExactInverse(x1, y1, z1, tx, ty, tz, rx, ry, rz, ds float64) (x, y, z float64) {
r00, r01, r02, r10, r11, r12, r20, r21, r22 := rotationMatrixExact(rx*asec, ry*asec, rz*asec)
s := 1 + ds*ppm
dx, dy, dz := x1-tx, y1-ty, z1-tz
x = (r00*dx + r10*dy + r20*dz) / s
y = (r01*dx + r11*dy + r21*dz) / s
z = (r02*dx + r12*dy + r22*dz) / s
return x, y, z
}
func rotationMatrixExact(rx, ry, rz float64) (r00, r01, r02, r10, r11, r12, r20, r21, r22 float64) {
sx, cx := math.Sincos(rx)
sy, cy := math.Sincos(ry)
sz, cz := math.Sincos(rz)
r00 = cy * cz
r01 = -cx*sz + sx*sy*cz
r02 = sx*sz + cx*sy*cz
r10 = cy * sz
r11 = cx*cz + sx*sy*sz
r12 = -sx*cz + cx*sy*sz
r20 = -sy
r21 = sx * cy
r22 = cx * cy
return
}
package crs
import (
"errors"
"fmt"
"math"
"reflect"
"strconv"
"strings"
)
func formatFloat(v float64) string {
return strconv.FormatFloat(v, 'f', -1, 64)
}
func degree(r float64) float64 {
return r * 180 / math.Pi
}
func radian(g float64) float64 {
return g * math.Pi / 180
}
func intPow(val float64, times int) float64 {
result := 1.0
for range times {
result *= val
}
return result
}
func round(val float64, dec int) float64 {
factor := math.Pow(10, float64(dec))
r := math.Round(val*factor) / factor
if r == -0 {
return 0
}
return r
}
func sin2(r float64) float64 {
return intPow(math.Sin(r), 2)
}
func sign(x float64) float64 {
if x < 0 {
return -1
}
return 1
}
func clamp(x, min, max float64) float64 {
if x < min {
return min
}
if x > max {
return max
}
return x
}
func authalicQ(sinPhi, e, e2 float64) float64 {
if e < 1e-12 {
return 2 * sinPhi
}
es := e * sinPhi
return (1 - e2) * (sinPhi/(1-e2*sinPhi*sinPhi) - (1/(2*e))*math.Log((1-es)/(1+es)))
}
func authalicToGeodetic(beta float64, s Spheroid) float64 {
e2, e4, e6 := s.E2(), s.E4(), s.E6()
return beta +
(e2/3+31*e4/180+517*e6/5040)*math.Sin(2*beta) +
(23*e4/360+251*e6/3780)*math.Sin(4*beta) +
(761*e6/45360)*math.Sin(6*beta)
}
func build(txt string) builder {
b := &strings.Builder{}
if len(txt) > 0 {
b.WriteString(txt)
}
return builder{
s: b,
}
}
type builder struct {
s *strings.Builder
err error
}
func (b builder) addKey(key string) builder {
b.s.WriteString(" " + key)
return b
}
func (b builder) addErr(err error) builder {
b.err = errors.Join(b.err, err)
return b
}
func (b builder) addAll(kv ...any) builder {
for i := 0; i+1 < len(kv); i += 2 {
k := kv[i]
kk, ok := k.(string)
if !ok {
return b.addErr(fmt.Errorf("invalid key: %v", k))
}
b = b.add(kk, kv[i+1])
}
return b
}
func (b builder) add(key string, val any) builder {
switch v := val.(type) {
case nil:
return b.addKey(key)
case string:
if v == "" {
return b
}
case fmt.Stringer:
val = v.String()
if val == "" {
return b
}
case float64:
if v == 0 {
return b
}
val = formatFloat(v)
}
if key == "" {
if val == nil {
return b
}
r := reflect.ValueOf(val)
if r.Kind() == reflect.Slice && r.Type().Elem().Implements(reflect.TypeFor[fmt.Stringer]()) {
for l := range r.Len() {
fmt.Fprintf(b.s, " %v", r.Index(l))
}
return b
}
b.s.WriteString(" " + r.String())
return b
}
fmt.Fprintf(b.s, " %s=%v", key, val)
return b
}
func (b builder) String() string {
if b.err != nil {
return b.err.Error()
}
return strings.TrimSpace(b.s.String())
}
package crs
import (
"fmt"
)
type HorizontalGrid string
func (hg HorizontalGrid) String() string {
return fmt.Sprintf("operation=horizontal_grid grid=%s", string(hg))
}
func (hg HorizontalGrid) FromTarget(source Spheroid, target Spheroid, lon0 float64, lat0 float64, h0 float64) (float64, float64, float64, error) {
grid, err := loadGrid(string(hg))
if err != nil {
return 0, 0, 0, err
}
return grid.FromTarget(source, target, lon0, lat0, h0)
}
func (hg HorizontalGrid) ToTarget(source Spheroid, target Spheroid, lon float64, lat float64, h float64) (float64, float64, float64, error) {
grid, err := loadGrid(string(hg))
if err != nil {
return 0, 0, 0, err
}
return grid.ToTarget(source, target, lon, lat, h)
}
type horizontalGridData struct {
File string
subGrids []subGrid
}
func (g *horizontalGridData) FromTarget(source Spheroid, target Spheroid, lon0 float64, lat0 float64, h0 float64) (float64, float64, float64, error) {
lon, lat, err := g.FromWGS84(lon0, lat0)
if err != nil {
return 0, 0, 0, err
}
return lon, lat, h0, nil
}
func (g *horizontalGridData) ToTarget(source Spheroid, target Spheroid, lon float64, lat float64, h float64) (float64, float64, float64, error) {
lon0, lat0, err := g.ToWGS84(lon, lat)
if err != nil {
return 0, 0, 0, err
}
return lon0, lat0, h, nil
}
type Inverse struct {
Operation Operation
}
func (i Inverse) String() string {
return fmt.Sprintf("%s inverse", i.Operation)
}
func (i Inverse) FromTarget(source Spheroid, target Spheroid, lon0 float64, lat0 float64, h0 float64) (float64, float64, float64, error) {
return i.Operation.ToTarget(target, source, lon0, lat0, h0)
}
func (i Inverse) ToTarget(source Spheroid, target Spheroid, lon float64, lat float64, h float64) (float64, float64, float64, error) {
return i.Operation.FromTarget(target, source, lon, lat, h)
}
package crs
import (
"math"
)
type HotineObliqueMercatorA struct {
Lonf, Latf, Alpha, Gamma, Scale, Eastf, Northf float64
}
func (cs HotineObliqueMercatorA) String() string {
return build("hotine_oblique_mercator_a").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"azimuth", cs.Alpha,
"gamma", cs.Gamma,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
type hotineConsts struct {
A, B, H, gamma0, lon0, gammaC float64
}
func (cs HotineObliqueMercatorA) consts(s Spheroid) hotineConsts {
a := s.A
e := s.E()
e2 := s.E2()
kC := cs.Scale
if kC == 0 {
kC = 1
}
phiC := radian(cs.Latf)
lambdaC := radian(cs.Lonf)
alphaC := radian(cs.Alpha)
gammaC := radian(cs.Gamma)
sinPhi := math.Sin(phiC)
cosPhi := math.Cos(phiC)
sLat := sign(sinPhi)
B := math.Sqrt(1 + (e2*math.Pow(cosPhi, 4))/(1-e2))
A := a * B * kC * math.Sqrt(1-e2) / (1 - e2*sinPhi*sinPhi)
t0 := math.Tan(math.Pi/4-phiC/2) /
math.Pow((1-e*sinPhi)/(1+e*sinPhi), e/2)
D := B * math.Sqrt(1-e2) / (cosPhi * math.Sqrt(1-e2*sinPhi*sinPhi))
D2 := D * D
if D < 1 {
D2 = 1
}
F := D + math.Sqrt(D2-1)*sLat
H := F * math.Pow(t0, B)
G := (F - 1/F) / 2
gamma0 := math.Asin(clamp(math.Sin(alphaC)/D, -1, 1))
lon0 := lambdaC - math.Asin(clamp(G*math.Tan(gamma0), -1, 1))/B
return hotineConsts{A: A, B: B, H: H, gamma0: gamma0, lon0: lon0, gammaC: gammaC}
}
func (cs HotineObliqueMercatorA) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
e := s.E()
c := cs.consts(s)
phi := radian(lat)
lambda := radian(lon)
t := math.Tan(math.Pi/4-phi/2) /
math.Pow((1-e*math.Sin(phi))/(1+e*math.Sin(phi)), e/2)
Q := c.H / math.Pow(t, c.B)
S := (Q - 1/Q) / 2
T := (Q + 1/Q) / 2
V := math.Sin(c.B * (lambda - c.lon0))
U := (-V*math.Cos(c.gamma0) + S*math.Sin(c.gamma0)) / T
v := c.A * math.Log((1-U)/(1+U)) / (2 * c.B)
u := c.A * math.Atan2(
S*math.Cos(c.gamma0)+V*math.Sin(c.gamma0),
math.Cos(c.B*(lambda-c.lon0)),
) / c.B
east := v*math.Cos(c.gammaC) + u*math.Sin(c.gammaC) + cs.Eastf
north := u*math.Cos(c.gammaC) - v*math.Sin(c.gammaC) + cs.Northf
return east, north, h
}
func (cs HotineObliqueMercatorA) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
e2 := s.E2()
e4 := e2 * e2
e6 := e4 * e2
e8 := e4 * e4
c := cs.consts(s)
v := (east-cs.Eastf)*math.Cos(c.gammaC) - (north-cs.Northf)*math.Sin(c.gammaC)
u := (north-cs.Northf)*math.Cos(c.gammaC) + (east-cs.Eastf)*math.Sin(c.gammaC)
Q := math.Exp(-(c.B * v / c.A))
S := (Q - 1/Q) / 2
T := (Q + 1/Q) / 2
V := math.Sin(c.B * u / c.A)
U := (V*math.Cos(c.gamma0) + S*math.Sin(c.gamma0)) / T
t := math.Pow(c.H/math.Sqrt((1+U)/(1-U)), 1/c.B)
chi := math.Pi/2 - 2*math.Atan(t)
phi := chi +
(e2/2+5*e4/24+e6/12+13*e8/360)*math.Sin(2*chi) +
(7*e4/48+29*e6/240+811*e8/11520)*math.Sin(4*chi) +
(7*e6/120+81*e8/1120)*math.Sin(6*chi) +
(4279*e8/161280)*math.Sin(8*chi)
lon := c.lon0 - math.Atan2(
S*math.Cos(c.gamma0)-V*math.Sin(c.gamma0),
math.Cos(c.B*u/c.A),
)/c.B
return degree(lon), degree(phi), h
}
package crs
import (
"math"
)
type HyperbolicCassiniSoldner struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs HyperbolicCassiniSoldner) String() string {
return build("hyperbolic_cassini_soldner").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs HyperbolicCassiniSoldner) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
base := CassiniSoldner(cs)
east, north, h := base.FromGeographic(s, lon, lat, h)
phi := radian(lat)
sinPhi := math.Sin(phi)
e2 := s.E2()
oneMe2Sin2 := 1 - e2*sinPhi*sinPhi
nu := s.A / math.Sqrt(oneMe2Sin2)
rho := s.A * (1 - e2) / math.Pow(oneMe2Sin2, 1.5)
x := north - cs.Northf
north = cs.Northf + x - x*x*x/(6*rho*nu)
return east, north, h
}
func (cs HyperbolicCassiniSoldner) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
phi0 := radian(cs.Latf)
target := north - cs.Northf
x := target
e2 := s.E2()
for range 10 {
m1 := meridianDistance(s, phi0) + x
phi1 := footpointLatitude(s, m1)
sinPhi1 := math.Sin(phi1)
oneMe2Sin2 := 1 - e2*sinPhi1*sinPhi1
nu := s.A / math.Sqrt(oneMe2Sin2)
rho := s.A * (1 - e2) / math.Pow(oneMe2Sin2, 1.5)
f := x - x*x*x/(6*rho*nu) - target
df := 1 - x*x/(2*rho*nu)
dx := f / df
x -= dx
if math.Abs(dx) < 1e-12 {
break
}
}
base := CassiniSoldner(cs)
return base.ToGeographic(s, east, cs.Northf+x, h)
}
package crs
import (
"math"
)
type Krovak struct {
Lonf float64
Latf float64
Alpha float64
Scale float64
Eastf float64
Northf float64
}
func (cs Krovak) String() string {
return build("krovak").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"alpha", cs.Alpha,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (k Krovak) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
phic := radian(k.Latf)
lambda0 := radian(k.Lonf)
phip := radian(78.5)
alphac := radian(k.Alpha)
A := s.A * math.Sqrt(1-s.E2()) / (1 - s.E2()*sin2(phic))
B := math.Sqrt(1 + (s.E2() * intPow(math.Cos(phic), 4) / (1 - s.E2())))
gamma0 := math.Asin(math.Sin(phic) / B)
t0 := math.Tan(math.Pi/4+gamma0/2) * math.Pow((1+s.E()*math.Sin(phic))/(1-s.E()*math.Sin(phic)), s.E()*B/2) / math.Pow(math.Tan(math.Pi/4+phic/2), B)
n := math.Sin(phip)
r0 := k.Scale * A / math.Tan(phip)
phi := radian(lat)
lambda := radian(lon)
U := 2 * (math.Atan(t0*math.Pow(math.Tan(phi/2+math.Pi/4), B)/math.Pow((1+s.E()*math.Sin(phi))/(1-s.E()*math.Sin(phi)), s.E()*B/2)) - math.Pi/4)
V := B * (lambda0 - lambda)
T := math.Asin(math.Cos(alphac)*math.Sin(U) + math.Sin(alphac)*math.Cos(U)*math.Cos(V))
D := math.Asin(math.Cos(U) * math.Sin(V) / math.Cos(T))
theta := n * D
r := r0 * math.Pow(math.Tan(math.Pi/4+phip/2), n) / math.Pow(math.Tan(T/2+math.Pi/4), n)
Xp := r * math.Cos(theta)
Yp := r * math.Sin(theta)
return -(Yp + k.Eastf), -(Xp + k.Northf), h
}
func (k Krovak) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
phic := radian(k.Latf)
lambda0 := radian(k.Lonf)
phip := radian(78.5)
alphac := radian(k.Alpha)
A := s.A * math.Sqrt(1-s.E2()) / (1 - s.E2()*sin2(phic))
B := math.Sqrt(1 + (s.E2() * intPow(math.Cos(phic), 4) / (1 - s.E2())))
gamma0 := math.Asin(math.Sin(phic) / B)
t0 := math.Tan(math.Pi/4+gamma0/2) * math.Pow((1+s.E()*math.Sin(phic))/(1-s.E()*math.Sin(phic)), s.E()*B/2) / math.Pow(math.Tan(math.Pi/4+phic/2), B)
n := math.Sin(phip)
r0 := k.Scale * A / math.Tan(phip)
Xpi := (-north) - k.Northf
Ypi := (-east) - k.Eastf
ri := math.Sqrt(intPow(Xpi, 2) + intPow(Ypi, 2))
thetai := math.Atan2(Ypi, Xpi)
di := thetai / math.Sin(phip)
ti := 2 * (math.Atan(math.Pow(r0/ri, 1/n)*math.Tan(math.Pi/4+phip/2)) - math.Pi/4)
ui := math.Asin(math.Cos(alphac)*math.Sin(ti) - math.Sin(alphac)*math.Cos(ti)*math.Cos(di))
vi := math.Asin(math.Cos(ti) * math.Sin(di) / math.Cos(ui))
phi := ui
for range 12 {
phiNext := 2 * (math.Atan(math.Pow(t0, -1/B)*math.Pow(math.Tan(ui/2+math.Pi/4), 1/B)*math.Pow((1+s.E()*math.Sin(phi))/(1-s.E()*math.Sin(phi)), s.E()/2)) - math.Pi/4)
if math.Abs(phiNext-phi) < 1e-12 {
phi = phiNext
break
}
phi = phiNext
}
lambda := lambda0 - vi/B
return degree(lambda), degree(phi), h
}
package crs
import (
"math"
)
const (
krovakModX0 = 1089000.0
krovakModY0 = 654000.0
)
const (
krovakModC1 = 2.946529277e-02
krovakModC2 = 2.515965696e-02
krovakModC3 = 1.193845912e-07
krovakModC4 = -4.668270147e-07
krovakModC5 = 9.233980362e-12
krovakModC6 = 1.523735715e-12
krovakModC7 = 1.696780024e-18
krovakModC8 = 4.408314235e-18
krovakModC9 = -8.331083518e-24
krovakModC10 = -3.689471323e-24
)
type KrovakModified struct {
Lonf, Latf, Alpha, Scale, Eastf, Northf float64
}
func (cs KrovakModified) String() string {
return build("krovak_modified").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"alpha", cs.Alpha,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs KrovakModified) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
xp, yp := krovakXY(s, lon, lat, cs.Lonf, cs.Latf, cs.Alpha, cs.Scale)
xp, yp = krovakModifiedApply(xp, yp)
return -(yp + cs.Eastf), -(xp + cs.Northf), h
}
func (cs KrovakModified) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
xp, yp := krovakModifiedInvert((-north)-cs.Northf, (-east)-cs.Eastf)
return krovakLonLat(s, xp, yp, cs.Lonf, cs.Latf, cs.Alpha, cs.Scale, h)
}
func krovakXY(s Spheroid, lon, lat, lonf, latf, alpha, scale float64) (xp, yp float64) {
phic := radian(latf)
lambda0 := radian(lonf)
phip := radian(78.5)
alphac := radian(alpha)
A := s.A * math.Sqrt(1-s.E2()) / (1 - s.E2()*sin2(phic))
B := math.Sqrt(1 + (s.E2() * intPow(math.Cos(phic), 4) / (1 - s.E2())))
gamma0 := math.Asin(math.Sin(phic) / B)
t0 := math.Tan(math.Pi/4+gamma0/2) * math.Pow((1+s.E()*math.Sin(phic))/(1-s.E()*math.Sin(phic)), s.E()*B/2) / math.Pow(math.Tan(math.Pi/4+phic/2), B)
n := math.Sin(phip)
r0 := scale * A / math.Tan(phip)
phi := radian(lat)
lambda := radian(lon)
U := 2 * (math.Atan(t0*math.Pow(math.Tan(phi/2+math.Pi/4), B)/math.Pow((1+s.E()*math.Sin(phi))/(1-s.E()*math.Sin(phi)), s.E()*B/2)) - math.Pi/4)
V := B * (lambda0 - lambda)
T := math.Asin(math.Cos(alphac)*math.Sin(U) + math.Sin(alphac)*math.Cos(U)*math.Cos(V))
D := math.Asin(math.Cos(U) * math.Sin(V) / math.Cos(T))
theta := n * D
r := r0 * math.Pow(math.Tan(math.Pi/4+phip/2), n) / math.Pow(math.Tan(T/2+math.Pi/4), n)
return r * math.Cos(theta), r * math.Sin(theta)
}
func krovakLonLat(s Spheroid, xp, yp, lonf, latf, alpha, scale, h float64) (float64, float64, float64) {
phic := radian(latf)
lambda0 := radian(lonf)
phip := radian(78.5)
alphac := radian(alpha)
A := s.A * math.Sqrt(1-s.E2()) / (1 - s.E2()*sin2(phic))
B := math.Sqrt(1 + (s.E2() * intPow(math.Cos(phic), 4) / (1 - s.E2())))
gamma0 := math.Asin(math.Sin(phic) / B)
t0 := math.Tan(math.Pi/4+gamma0/2) * math.Pow((1+s.E()*math.Sin(phic))/(1-s.E()*math.Sin(phic)), s.E()*B/2) / math.Pow(math.Tan(math.Pi/4+phic/2), B)
n := math.Sin(phip)
r0 := scale * A / math.Tan(phip)
ri := math.Sqrt(intPow(xp, 2) + intPow(yp, 2))
thetai := math.Atan2(yp, xp)
di := thetai / math.Sin(phip)
ti := 2 * (math.Atan(math.Pow(r0/ri, 1/n)*math.Tan(math.Pi/4+phip/2)) - math.Pi/4)
ui := math.Asin(math.Cos(alphac)*math.Sin(ti) - math.Sin(alphac)*math.Cos(ti)*math.Cos(di))
vi := math.Asin(math.Cos(ti) * math.Sin(di) / math.Cos(ui))
phi := ui
for range 12 {
phiNext := 2 * (math.Atan(math.Pow(t0, -1/B)*math.Pow(math.Tan(ui/2+math.Pi/4), 1/B)*math.Pow((1+s.E()*math.Sin(phi))/(1-s.E()*math.Sin(phi)), s.E()/2)) - math.Pi/4)
if math.Abs(phiNext-phi) < 1e-12 {
phi = phiNext
break
}
phi = phiNext
}
return degree(lambda0 - vi/B), degree(phi), h
}
func krovakModifiedDXDY(xr, yr float64) (dX, dY float64) {
xr2 := xr * xr
yr2 := yr * yr
xr4 := xr2 * xr2
yr4 := yr2 * yr2
dX = krovakModC1 + krovakModC3*xr - krovakModC4*yr - 2*krovakModC6*xr*yr + krovakModC5*(xr2-yr2) +
krovakModC7*xr*(xr2-3*yr2) - krovakModC8*yr*(3*xr2-yr2) +
4*krovakModC9*xr*yr*(xr2-yr2) + krovakModC10*(xr4+yr4-6*xr2*yr2)
dY = krovakModC2 + krovakModC3*yr + krovakModC4*xr + 2*krovakModC5*xr*yr + krovakModC6*(xr2-yr2) +
krovakModC8*xr*(xr2-3*yr2) + krovakModC7*yr*(3*xr2-yr2) -
4*krovakModC10*xr*yr*(xr2-yr2) + krovakModC9*(xr4+yr4-6*xr2*yr2)
return dX, dY
}
func krovakModifiedApply(xp, yp float64) (float64, float64) {
dX, dY := krovakModifiedDXDY(xp-krovakModX0, yp-krovakModY0)
return xp - dX, yp - dY
}
func krovakModifiedInvert(xp, yp float64) (float64, float64) {
u, v := xp, yp
for range 10 {
dX, dY := krovakModifiedDXDY(u-krovakModX0, v-krovakModY0)
nu, nv := xp+dX, yp+dY
if math.Abs(nu-u) < 1e-12 && math.Abs(nv-v) < 1e-12 {
return nu, nv
}
u, v = nu, nv
}
return u, v
}
package crs
import (
"math"
)
type LabordeObliqueMercator struct {
Lonf, Latf, Alpha, Scale, Eastf, Northf float64
}
func (cs LabordeObliqueMercator) String() string {
return build("laborde_oblique_mercator").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"azimuth", cs.Alpha,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
type labordeConsts struct {
kRg, p0s, A, C, Ca, Cb, Cc, Cd float64
}
func (cs LabordeObliqueMercator) consts(s Spheroid) labordeConsts {
k0 := cs.Scale
if k0 == 0 {
k0 = 1
}
phi0 := radian(cs.Latf)
az := radian(cs.Alpha)
e := s.E()
e2 := s.E2()
sinp := math.Sin(phi0)
t := 1 - e2*sinp*sinp
N := 1 / math.Sqrt(t)
R := (1 - e2) * N / t
kRg := k0 * math.Sqrt(N*R)
p0s := math.Atan(math.Sqrt(R/N) * math.Tan(phi0))
A := sinp / math.Sin(p0s)
te := e * sinp
C := 0.5*e*A*math.Log((1+te)/(1-te)) - A*math.Log(math.Tan(math.Pi/4+0.5*phi0)) +
math.Log(math.Tan(math.Pi/4+0.5*p0s))
tAz := az + az
Cb := 1 / (12 * kRg * kRg)
Ca := (1 - math.Cos(tAz)) * Cb
Cb *= math.Sin(tAz)
Cc := 3 * (Ca*Ca - Cb*Cb)
Cd := 6 * Ca * Cb
return labordeConsts{kRg: kRg, p0s: p0s, A: A, C: C, Ca: Ca, Cb: Cb, Cc: Cc, Cd: Cd}
}
func (cs LabordeObliqueMercator) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
q := cs.consts(s)
e := s.E()
phi := radian(lat)
lam := radian(lon - cs.Lonf)
V1 := q.A * math.Log(math.Tan(math.Pi/4+0.5*phi))
t := e * math.Sin(phi)
V2 := 0.5 * e * q.A * math.Log((1+t)/(1-t))
ps := 2 * (math.Atan(math.Exp(V1-V2+q.C)) - math.Pi/4)
I1 := ps - q.p0s
cosps := math.Cos(ps)
sinps := math.Sin(ps)
cosps2 := cosps * cosps
sinps2 := sinps * sinps
I4 := q.A * cosps
I2 := 0.5 * q.A * I4 * sinps
I3 := I2 * q.A * q.A * (5*cosps2 - sinps2) / 12
I6 := I4 * q.A * q.A
I5 := I6 * (cosps2 - sinps2) / 6
I6 *= q.A * q.A * (5*cosps2*cosps2 + sinps2*(sinps2-18*cosps2)) / 120
t2 := lam * lam
x := q.kRg * lam * (I4 + t2*(I5+t2*I6))
y := q.kRg * (I1 + t2*(I2+t2*I3))
x2 := x * x
y2 := y * y
V1 = 3*x*y2 - x*x2
V2 = y*y2 - 3*x2*y
x += q.Ca*V1 + q.Cb*V2
y += q.Ca*V2 - q.Cb*V1
return cs.Eastf + s.A*x, cs.Northf + s.A*y, h
}
func (cs LabordeObliqueMercator) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
q := cs.consts(s)
e := s.E()
e2 := s.E2()
phi0 := radian(cs.Latf)
x := (east - cs.Eastf) / s.A
y := (north - cs.Northf) / s.A
x2 := x * x
y2 := y * y
V1 := 3*x*y2 - x*x2
V2 := y*y2 - 3*x2*y
V3 := x * (5*y2*y2 + x2*(-10*y2+x2))
V4 := y * (5*x2*x2 + y2*(-10*x2+y2))
x += -q.Ca*V1 - q.Cb*V2 + q.Cc*V3 + q.Cd*V4
y += q.Cb*V1 - q.Ca*V2 - q.Cd*V3 + q.Cc*V4
ps := q.p0s + y/q.kRg
pe := ps + phi0 - q.p0s
for range 20 {
V1 = q.A * math.Log(math.Tan(math.Pi/4+0.5*pe))
tpe := e * math.Sin(pe)
V2 = 0.5 * e * q.A * math.Log((1+tpe)/(1-tpe))
t := ps - 2*(math.Atan(math.Exp(V1-V2+q.C))-math.Pi/4)
pe += t
if math.Abs(t) < 1e-10 {
break
}
}
t := e * math.Sin(pe)
t = 1 - t*t
Re := (1 - e2) / (t * math.Sqrt(t))
k0 := cs.Scale
if k0 == 0 {
k0 = 1
}
tanps := math.Tan(ps)
t2 := tanps * tanps
s2 := q.kRg * q.kRg
d := Re * k0 * q.kRg
I7 := tanps / (2 * d)
I8 := tanps * (5 + 3*t2) / (24 * d * s2)
d = math.Cos(ps) * q.kRg * q.A
I9 := 1 / d
d *= s2
I10 := (1 + 2*t2) / (6 * d)
I11 := (5 + t2*(28+24*t2)) / (120 * d * s2)
x2 = x * x
phi := pe + x2*(-I7+I8*x2)
lam := x * (I9 + x2*(-I10+x2*I11))
return degree(lam + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"math"
)
type LambertAzimuthalEqualArea struct {
Lonf float64
Latf float64
Eastf float64
Northf float64
}
func (cs LambertAzimuthalEqualArea) String() string {
return build("lambert_azimuthal_equal_area").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (laea LambertAzimuthalEqualArea) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
phi0 := radian(laea.Latf)
lambda0 := radian(laea.Lonf)
e, e2 := s.E(), s.E2()
q0 := authalicQ(math.Sin(phi0), e, e2)
qp := authalicQ(1, e, e2)
beta0 := math.Asin(q0 / qp)
rq := s.A * math.Sqrt(qp/2)
g := s.A * (math.Cos(phi0) / math.Sqrt(1-e2*sin2(phi0))) / (rq * math.Cos(beta0))
phi := radian(lat)
lambda := radian(lon)
q := authalicQ(math.Sin(phi), e, e2)
beta := math.Asin(q / qp)
b := rq * math.Sqrt(2/(1+math.Sin(beta0)*math.Sin(beta)+(math.Cos(beta0)*math.Cos(beta)*math.Cos(lambda-lambda0))))
east := laea.Eastf + ((b * g) * (math.Cos(beta) * math.Sin(lambda-lambda0)))
north := laea.Northf + (b/g)*((math.Cos(beta0)*math.Sin(beta))-(math.Sin(beta0)*math.Cos(beta)*math.Cos(lambda-lambda0)))
return east, north, h
}
func (laea LambertAzimuthalEqualArea) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
phi0 := radian(laea.Latf)
lambda0 := radian(laea.Lonf)
e, e2 := s.E(), s.E2()
q0 := authalicQ(math.Sin(phi0), e, e2)
qp := authalicQ(1, e, e2)
beta0 := math.Asin(q0 / qp)
rq := s.A * math.Sqrt(qp/2)
g := s.A * (math.Cos(phi0) / math.Sqrt(1-e2*sin2(phi0))) / (rq * math.Cos(beta0))
rho := math.Sqrt(intPow((east-laea.Eastf)/g, 2) + intPow(g*(north-laea.Northf), 2))
if rho < 1e-14 {
return laea.Lonf, laea.Latf, h
}
c := 2 * math.Asin(rho/(2*rq))
betai := math.Asin((math.Cos(c) * math.Sin(beta0)) + ((g * (north - laea.Northf) * math.Sin(c) * math.Cos(beta0)) / rho))
phi := authalicToGeodetic(betai, s)
lambda := lambda0 + math.Atan2((east-laea.Eastf)*math.Sin(c), (g*rho*math.Cos(beta0)*math.Cos(c)-intPow(g, 2)*(north-laea.Northf)*math.Sin(beta0)*math.Sin(c)))
return degree(lambda), degree(phi), h
}
package crs
import (
"math"
)
type LambertAzimuthalEqualAreaSpherical struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs LambertAzimuthalEqualAreaSpherical) String() string {
return build("lambert_azimuthal_equal_area_spherical").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs LambertAzimuthalEqualAreaSpherical) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
r := s.A
phi0 := radian(cs.Latf)
lam0 := radian(cs.Lonf)
phi := radian(lat)
lam := radian(lon)
sinPhi0, cosPhi0 := math.Sincos(phi0)
sinPhi, cosPhi := math.Sincos(phi)
dLam := lam - lam0
denom := 1 + sinPhi0*sinPhi + cosPhi0*cosPhi*math.Cos(dLam)
if denom < 1e-14 {
return cs.Eastf, cs.Northf, h
}
k := math.Sqrt(2 / denom)
east := cs.Eastf + r*k*cosPhi*math.Sin(dLam)
north := cs.Northf + r*k*(cosPhi0*sinPhi-sinPhi0*cosPhi*math.Cos(dLam))
return east, north, h
}
func (cs LambertAzimuthalEqualAreaSpherical) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
r := s.A
phi0 := radian(cs.Latf)
lam0 := radian(cs.Lonf)
x := (east - cs.Eastf) / r
y := (north - cs.Northf) / r
rho := math.Hypot(x, y)
if rho < 1e-14 {
return cs.Lonf, cs.Latf, h
}
c := 2 * math.Asin(clamp(rho/2, -1, 1))
sinC, cosC := math.Sincos(c)
sinPhi0, cosPhi0 := math.Sincos(phi0)
phi := math.Asin(clamp(cosC*sinPhi0+y*sinC*cosPhi0/rho, -1, 1))
lam := lam0 + math.Atan2(x*sinC, rho*cosPhi0*cosC-y*sinPhi0*sinC)
return degree(lam), degree(phi), h
}
package crs
import (
"math"
)
type LambertConformalConic struct {
Lonf float64
Latf float64
Scale float64
Eastf float64
Northf float64
}
func (cs LambertConformalConic) String() string {
return build("lambert_conformal_conic").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (lcc LambertConformalConic) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
phif := radian(lcc.Latf)
lambdaf := radian(lcc.Lonf)
tf := math.Tan(math.Pi/4-phif/2) / math.Pow((1-s.E()*math.Sin(phif))/(1+s.E()*math.Sin(phif)), s.E()/2)
m1 := math.Cos(phif) / math.Sqrt(1-s.E2()*sin2(phif))
n := math.Sin(phif)
f := lcc.Scale * m1 / (n * math.Pow(tf, n))
rf := s.A * f * math.Pow(tf, n)
phi := radian(lat)
lambda := radian(lon)
t := math.Tan(math.Pi/4-phi/2) / math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2)
r := s.A * f * math.Pow(t, n)
theta := n * (lambda - lambdaf)
east := lcc.Eastf + r*math.Sin(theta)
north := lcc.Northf + rf - r*math.Cos(theta)
return east, north, h
}
func (lcc LambertConformalConic) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
phif := radian(lcc.Latf)
lambdaf := radian(lcc.Lonf)
tf := math.Tan(math.Pi/4-phif/2) / math.Pow((1-s.E()*math.Sin(phif))/(1+s.E()*math.Sin(phif)), s.E()/2)
m1 := math.Cos(phif) / math.Sqrt(1-s.E2()*sin2(phif))
n := math.Sin(phif)
f := lcc.Scale * m1 / (n * math.Pow(tf, n))
rf := s.A * f * math.Pow(tf, n)
ri := math.Hypot(east-lcc.Eastf, rf-(north-lcc.Northf))
ti := math.Pow(ri/(s.A*f), 1/n)
theta := math.Atan2(east-lcc.Eastf, rf-(north-lcc.Northf))
phi := math.Pi/2 - 2*math.Atan(ti)
for range 6 {
next := math.Pi/2 - 2*math.Atan(ti*math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2))
if math.Abs(next-phi) < 1e-12 {
phi = next
break
}
phi = next
}
lambda := theta/n + lambdaf
return degree(lambda), degree(phi), h
}
package crs
import (
"math"
)
// LambertConformalConic1SPVariantB is EPSG method 1102.
// Same cone as LCC 1SP at Lat0 (natural origin); FE/FN apply at Latf (false origin).
type LambertConformalConic1SPVariantB struct {
Lonf, Lat0, Latf, Scale, Eastf, Northf float64
}
func (cs LambertConformalConic1SPVariantB) String() string {
return build("lambert_conformal_conic_1sp_variant_b").addAll(
"lonf", cs.Lonf,
"lat0", cs.Lat0,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs LambertConformalConic1SPVariantB) consts(s Spheroid) (n, f, rf float64) {
k0 := cs.Scale
if k0 == 0 {
k0 = 1
}
phi0 := radian(cs.Lat0)
phif := radian(cs.Latf)
t0 := math.Tan(math.Pi/4-phi0/2) / math.Pow((1-s.E()*math.Sin(phi0))/(1+s.E()*math.Sin(phi0)), s.E()/2)
tf := math.Tan(math.Pi/4-phif/2) / math.Pow((1-s.E()*math.Sin(phif))/(1+s.E()*math.Sin(phif)), s.E()/2)
m0 := math.Cos(phi0) / math.Sqrt(1-s.E2()*sin2(phi0))
n = math.Sin(phi0)
f = k0 * m0 / (n * math.Pow(t0, n))
rf = s.A * f * math.Pow(tf, n)
return n, f, rf
}
func (cs LambertConformalConic1SPVariantB) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
n, f, rf := cs.consts(s)
phi := radian(lat)
t := math.Tan(math.Pi/4-phi/2) / math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2)
r := s.A * f * math.Pow(t, n)
theta := n * (radian(lon) - radian(cs.Lonf))
east := cs.Eastf + r*math.Sin(theta)
north := cs.Northf + rf - r*math.Cos(theta)
return east, north, h
}
func (cs LambertConformalConic1SPVariantB) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
n, f, rf := cs.consts(s)
ri := math.Hypot(east-cs.Eastf, rf-(north-cs.Northf))
ti := math.Pow(ri/(s.A*f), 1/n)
theta := math.Atan2(east-cs.Eastf, rf-(north-cs.Northf))
phi := math.Pi/2 - 2*math.Atan(ti)
for range 6 {
next := math.Pi/2 - 2*math.Atan(ti*math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2))
if math.Abs(next-phi) < 1e-12 {
phi = next
break
}
phi = next
}
return degree(theta/n + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"math"
)
type LambertConformalConic2SP struct {
Lonf float64
Latf float64
Sp1 float64
Sp2 float64
Eastf float64
Northf float64
}
func (cs LambertConformalConic2SP) String() string {
return build("lambert_conformal_conic_2sp").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"sp1", cs.Sp1,
"sp2", cs.Sp2,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (lcc LambertConformalConic2SP) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
phif := radian(lcc.Latf)
phi1 := radian(lcc.Sp1)
phi2 := radian(lcc.Sp2)
lambdaf := radian(lcc.Lonf)
tf := math.Tan(math.Pi/4-phif/2) / math.Pow((1-s.E()*math.Sin(phif))/(1+s.E()*math.Sin(phif)), s.E()/2)
t1 := math.Tan(math.Pi/4-phi1/2) / math.Pow((1-s.E()*math.Sin(phi1))/(1+s.E()*math.Sin(phi1)), s.E()/2)
t2 := math.Tan(math.Pi/4-phi2/2) / math.Pow((1-s.E()*math.Sin(phi2))/(1+s.E()*math.Sin(phi2)), s.E()/2)
m1 := math.Cos(phi1) / math.Sqrt(1-s.E2()*sin2(phi1))
m2 := math.Cos(phi2) / math.Sqrt(1-s.E2()*sin2(phi2))
var n float64
if math.Abs(phi1-phi2) < 1e-14 {
n = math.Sin(phi1)
} else {
n = math.Log(m1/m2) / math.Log(t1/t2)
}
f := m1 / (n * math.Pow(t1, n))
rf := s.A * f * math.Pow(tf, n)
phi := radian(lat)
lambda := radian(lon)
t := math.Tan(math.Pi/4-phi/2) / math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2)
r := s.A * f * math.Pow(t, n)
theta := n * (lambda - lambdaf)
east := lcc.Eastf + r*math.Sin(theta)
north := lcc.Northf + rf - r*math.Cos(theta)
return east, north, h
}
func (lcc LambertConformalConic2SP) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
phif := radian(lcc.Latf)
phi1 := radian(lcc.Sp1)
phi2 := radian(lcc.Sp2)
lambdaf := radian(lcc.Lonf)
tf := math.Tan(math.Pi/4-phif/2) / math.Pow((1-s.E()*math.Sin(phif))/(1+s.E()*math.Sin(phif)), s.E()/2)
t1 := math.Tan(math.Pi/4-phi1/2) / math.Pow((1-s.E()*math.Sin(phi1))/(1+s.E()*math.Sin(phi1)), s.E()/2)
t2 := math.Tan(math.Pi/4-phi2/2) / math.Pow((1-s.E()*math.Sin(phi2))/(1+s.E()*math.Sin(phi2)), s.E()/2)
m1 := math.Cos(phi1) / math.Sqrt(1-s.E2()*sin2(phi1))
m2 := math.Cos(phi2) / math.Sqrt(1-s.E2()*sin2(phi2))
var n float64
if math.Abs(phi1-phi2) < 1e-14 {
n = math.Sin(phi1)
} else {
n = math.Log(m1/m2) / math.Log(t1/t2)
}
f := m1 / (n * math.Pow(t1, n))
rf := s.A * f * math.Pow(tf, n)
ri := math.Hypot(east-lcc.Eastf, rf-(north-lcc.Northf))
ti := math.Pow(ri/(s.A*f), 1/n)
theta := math.Atan2(east-lcc.Eastf, rf-(north-lcc.Northf))
phi := math.Pi/2 - 2*math.Atan(ti)
for range 6 {
next := math.Pi/2 - 2*math.Atan(ti*math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2))
if math.Abs(next-phi) < 1e-12 {
phi = next
break
}
phi = next
}
lambda := theta/n + lambdaf
return degree(lambda), degree(phi), h
}
package crs
import (
"math"
)
const belgiumLCCAngle = 29.2985 * math.Pi / (180 * 3600)
type LambertConformalConic2SPBelgium struct {
Lonf, Latf, Sp1, Sp2, Eastf, Northf float64
}
func (cs LambertConformalConic2SPBelgium) String() string {
return build("lambert_conformal_conic_2sp_belgium").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"sp1", cs.Sp1,
"sp2", cs.Sp2,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs LambertConformalConic2SPBelgium) consts(s Spheroid) (n, f, rf float64) {
phif := radian(cs.Latf)
phi1 := radian(cs.Sp1)
phi2 := radian(cs.Sp2)
tf := math.Tan(math.Pi/4-phif/2) / math.Pow((1-s.E()*math.Sin(phif))/(1+s.E()*math.Sin(phif)), s.E()/2)
t1 := math.Tan(math.Pi/4-phi1/2) / math.Pow((1-s.E()*math.Sin(phi1))/(1+s.E()*math.Sin(phi1)), s.E()/2)
t2 := math.Tan(math.Pi/4-phi2/2) / math.Pow((1-s.E()*math.Sin(phi2))/(1+s.E()*math.Sin(phi2)), s.E()/2)
m1 := math.Cos(phi1) / math.Sqrt(1-s.E2()*sin2(phi1))
m2 := math.Cos(phi2) / math.Sqrt(1-s.E2()*sin2(phi2))
if math.Abs(phi1-phi2) < 1e-14 {
n = math.Sin(phi1)
} else {
n = math.Log(m1/m2) / math.Log(t1/t2)
}
f = m1 / (n * math.Pow(t1, n))
rf = s.A * f * math.Pow(tf, n)
return n, f, rf
}
func (cs LambertConformalConic2SPBelgium) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
n, f, rf := cs.consts(s)
phi := radian(lat)
t := math.Tan(math.Pi/4-phi/2) / math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2)
r := s.A * f * math.Pow(t, n)
theta := n*(radian(lon)-radian(cs.Lonf)) - belgiumLCCAngle
east := cs.Eastf + r*math.Sin(theta)
north := cs.Northf + rf - r*math.Cos(theta)
return east, north, h
}
func (cs LambertConformalConic2SPBelgium) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
n, f, rf := cs.consts(s)
ri := math.Hypot(east-cs.Eastf, rf-(north-cs.Northf))
ti := math.Pow(ri/(s.A*f), 1/n)
theta := math.Atan2(east-cs.Eastf, rf-(north-cs.Northf)) + belgiumLCCAngle
phi := math.Pi/2 - 2*math.Atan(ti)
for range 6 {
next := math.Pi/2 - 2*math.Atan(ti*math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2))
if math.Abs(next-phi) < 1e-12 {
phi = next
break
}
phi = next
}
return degree(theta/n + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"math"
)
const michiganEllipsoidScale = 1.0000382
type LambertConformalConic2SPMichigan struct {
Lonf, Latf, Sp1, Sp2, Eastf, Northf float64
}
func (cs LambertConformalConic2SPMichigan) String() string {
return build("lambert_conformal_conic_2sp_michigan").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"sp1", cs.Sp1,
"sp2", cs.Sp2,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs LambertConformalConic2SPMichigan) consts(s Spheroid) (n, f, rf float64) {
phif := radian(cs.Latf)
phi1 := radian(cs.Sp1)
phi2 := radian(cs.Sp2)
tf := math.Tan(math.Pi/4-phif/2) / math.Pow((1-s.E()*math.Sin(phif))/(1+s.E()*math.Sin(phif)), s.E()/2)
t1 := math.Tan(math.Pi/4-phi1/2) / math.Pow((1-s.E()*math.Sin(phi1))/(1+s.E()*math.Sin(phi1)), s.E()/2)
t2 := math.Tan(math.Pi/4-phi2/2) / math.Pow((1-s.E()*math.Sin(phi2))/(1+s.E()*math.Sin(phi2)), s.E()/2)
m1 := math.Cos(phi1) / math.Sqrt(1-s.E2()*sin2(phi1))
m2 := math.Cos(phi2) / math.Sqrt(1-s.E2()*sin2(phi2))
if math.Abs(phi1-phi2) < 1e-14 {
n = math.Sin(phi1)
} else {
n = math.Log(m1/m2) / math.Log(t1/t2)
}
f = m1 / (n * math.Pow(t1, n))
rf = s.A * michiganEllipsoidScale * f * math.Pow(tf, n)
return n, f, rf
}
func (cs LambertConformalConic2SPMichigan) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
n, f, rf := cs.consts(s)
phi := radian(lat)
t := math.Tan(math.Pi/4-phi/2) / math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2)
r := s.A * michiganEllipsoidScale * f * math.Pow(t, n)
theta := n * (radian(lon) - radian(cs.Lonf))
east := cs.Eastf + r*math.Sin(theta)
north := cs.Northf + rf - r*math.Cos(theta)
return east, north, h
}
func (cs LambertConformalConic2SPMichigan) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
n, f, rf := cs.consts(s)
ri := math.Hypot(east-cs.Eastf, rf-(north-cs.Northf))
ti := math.Pow(ri/(s.A*michiganEllipsoidScale*f), 1/n)
theta := math.Atan2(east-cs.Eastf, rf-(north-cs.Northf))
phi := math.Pi/2 - 2*math.Atan(ti)
for range 6 {
next := math.Pi/2 - 2*math.Atan(ti*math.Pow((1-s.E()*math.Sin(phi))/(1+s.E()*math.Sin(phi)), s.E()/2))
if math.Abs(next-phi) < 1e-12 {
phi = next
break
}
phi = next
}
return degree(theta/n + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"math"
)
type LambertConicNearConformal struct {
Lonf, Latf, Scale, Eastf, Northf float64
}
func (cs LambertConicNearConformal) String() string {
return build("lambert_conic_near_conformal").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
type lccaConsts struct {
l, m0, r0, c, k0 float64
}
func (cs LambertConicNearConformal) consts(s Spheroid) lccaConsts {
k0 := cs.Scale
if k0 == 0 {
k0 = 1
}
phi0 := radian(cs.Latf)
l := math.Sin(phi0)
m0 := meridianDistance(s, phi0) / s.A
s2p0 := l * l
r0 := 1 / (1 - s.E2()*s2p0)
n0 := math.Sqrt(r0)
r0 *= (1 - s.E2()) * n0
tan0 := math.Tan(phi0)
return lccaConsts{
l: l, m0: m0, r0: n0 / tan0, c: 1 / (6 * r0 * n0), k0: k0,
}
}
func lccaFS(S, C float64) float64 { return S * (1 + S*S*C) }
func lccaFSp(S, C float64) float64 { return 1 + 3*S*S*C }
func (cs LambertConicNearConformal) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
q := cs.consts(s)
S := meridianDistance(s, radian(lat))/s.A - q.m0
dr := lccaFS(S, q.c)
r := q.r0 - dr
lamL := (radian(lon) - radian(cs.Lonf)) * q.l
east := cs.Eastf + s.A*q.k0*(r*math.Sin(lamL))
north := cs.Northf + s.A*q.k0*(q.r0-r*math.Cos(lamL))
return east, north, h
}
func (cs LambertConicNearConformal) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
q := cs.consts(s)
x := (east - cs.Eastf) / (s.A * q.k0)
y := (north - cs.Northf) / (s.A * q.k0)
theta := math.Atan2(x, q.r0-y)
dr := y - x*math.Tan(0.5*theta)
S := dr
for range 10 {
dif := (lccaFS(S, q.c) - dr) / lccaFSp(S, q.c)
S -= dif
if math.Abs(dif) < 1e-12 {
break
}
}
phi := footpointLatitude(s, s.A*(S+q.m0))
lam := radian(cs.Lonf) + theta/q.l
return degree(lam), degree(phi), h
}
package crs
import (
"math"
)
type LambertCylindricalEqualArea struct {
Lonf, Sp1, Eastf, Northf float64
}
func (cs LambertCylindricalEqualArea) String() string {
return build("lambert_cylindrical_equal_area").addAll(
"lonf", cs.Lonf,
"sp1", cs.Sp1,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs LambertCylindricalEqualArea) k0(s Spheroid) float64 {
phi1 := radian(cs.Sp1)
sin1 := math.Sin(phi1)
return math.Cos(phi1) / math.Sqrt(1-s.E2()*sin1*sin1)
}
func (cs LambertCylindricalEqualArea) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
k0 := cs.k0(s)
east := cs.Eastf + s.A*k0*(radian(lon)-radian(cs.Lonf))
north := cs.Northf + s.A*0.5*authalicQ(math.Sin(radian(lat)), s.E(), s.E2())/k0
return east, north, h
}
func (cs LambertCylindricalEqualArea) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
k0 := cs.k0(s)
qp := authalicQ(1, s.E(), s.E2())
lam := radian(cs.Lonf) + (east-cs.Eastf)/(s.A*k0)
sinBeta := clamp(2*(north-cs.Northf)*k0/(s.A*qp), -1, 1)
beta := math.Asin(sinBeta)
return degree(lam), degree(authalicToGeodetic(beta, s)), h
}
type LambertCylindricalEqualAreaSpherical struct {
Lonf, Sp1, Eastf, Northf float64
}
func (cs LambertCylindricalEqualAreaSpherical) String() string {
return build("lambert_cylindrical_equal_area_spherical").addAll(
"lonf", cs.Lonf,
"sp1", cs.Sp1,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs LambertCylindricalEqualAreaSpherical) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
r := s.A
k0 := math.Cos(radian(cs.Sp1))
east := cs.Eastf + r*k0*(radian(lon)-radian(cs.Lonf))
north := cs.Northf + r*math.Sin(radian(lat))/k0
return east, north, h
}
func (cs LambertCylindricalEqualAreaSpherical) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
r := s.A
k0 := math.Cos(radian(cs.Sp1))
lam := radian(cs.Lonf) + (east-cs.Eastf)/(r*k0)
phi := math.Asin(clamp((north-cs.Northf)*k0/r, -1, 1))
return degree(lam), degree(phi), h
}
package crs
import (
"math"
)
type LocalOrthographic struct {
Lonf, Latf, Azimuth, Scale, Eastf, Northf float64
}
func (cs LocalOrthographic) String() string {
return build("local_orthographic").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"azimuth", cs.Azimuth,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs LocalOrthographic) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
k0 := cs.Scale
if k0 == 0 {
k0 = 1
}
e2 := s.E2()
phi0 := radian(cs.Latf)
phi := radian(lat)
lam := radian(lon) - radian(cs.Lonf)
sinPhi0, cosPhi0 := math.Sincos(phi0)
sinPhi, cosPhi := math.Sincos(phi)
sinLam, cosLam := math.Sincos(lam)
nu0 := 1 / math.Sqrt(1-e2*sinPhi0*sinPhi0)
nu := 1 / math.Sqrt(1-e2*sinPhi*sinPhi)
xp := nu * cosPhi * sinLam
yp := nu*(sinPhi*cosPhi0-cosPhi*sinPhi0*cosLam) + e2*(nu0*sinPhi0-nu*sinPhi)*cosPhi0
sinA, cosA := math.Sincos(radian(cs.Azimuth))
east := cs.Eastf + s.A*k0*(cosA*xp-sinA*yp)
north := cs.Northf + s.A*k0*(sinA*xp+cosA*yp)
return east, north, h
}
func (cs LocalOrthographic) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
k0 := cs.Scale
if k0 == 0 {
k0 = 1
}
e2 := s.E2()
sinA, cosA := math.Sincos(radian(cs.Azimuth))
xf := (east - cs.Eastf) / (s.A * k0)
yf := (north - cs.Northf) / (s.A * k0)
x := cosA*xf + sinA*yf
y := -sinA*xf + cosA*yf
phi0 := radian(cs.Latf)
sinPhi0, cosPhi0 := math.Sincos(phi0)
nu0 := 1 / math.Sqrt(1-e2*sinPhi0*sinPhi0)
yShift := e2 * nu0 * sinPhi0 * cosPhi0
yScale := 1 / math.Sqrt(1-e2*cosPhi0*cosPhi0)
yRec := (y - yShift) / yScale
rh := math.Hypot(x, yRec)
sinc := rh
if sinc > 1 {
sinc = 1
}
cosc := math.Sqrt(math.Max(0, 1-sinc*sinc))
var phi, lam float64
if rh < 1e-14 {
phi, lam = phi0, 0
} else {
phi = math.Asin(clamp(cosc*sinPhi0+yRec*sinc*cosPhi0/rh, -1, 1))
lam = math.Atan2(x*sinc, rh*cosPhi0*cosc-yRec*sinPhi0*sinc)
}
for range 20 {
sinPhi, cosPhi := math.Sincos(phi)
sinLam, cosLam := math.Sincos(lam)
oneMinus := 1 - e2*sinPhi*sinPhi
nu := 1 / math.Sqrt(oneMinus)
xp := nu * cosPhi * sinLam
yp := nu*(sinPhi*cosPhi0-cosPhi*sinPhi0*cosLam) + e2*(nu0*sinPhi0-nu*sinPhi)*cosPhi0
rho := (1 - e2) * nu / oneMinus
j11 := -rho * sinPhi * sinLam
j12 := nu * cosPhi * cosLam
j21 := rho * (cosPhi*cosPhi0 + sinPhi*sinPhi0*cosLam)
j22 := nu * sinPhi0 * cosPhi * sinLam
d := j11*j22 - j12*j21
dx := x - xp
dy := y - yp
dPhi := (j22*dx - j12*dy) / d
dLam := (-j21*dx + j11*dy) / d
phi += dPhi
lam += dLam
if math.Abs(dPhi) < 1e-12 && math.Abs(dLam) < 1e-12 {
break
}
}
return degree(lam + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"math"
)
// MercatorA is EPSG method 9804 (PROJ +proj=merc with k / 1SP).
// Latf is unused in the formulas (natural origin is on the equator).
type MercatorA struct {
Lonf, Latf, Scale, Eastf, Northf float64
}
func (cs MercatorA) String() string {
return build("mercator_a").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs MercatorA) k0() float64 {
if cs.Scale == 0 {
return 1
}
return cs.Scale
}
func (cs MercatorA) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
return mercatorFromGeographic(s, lon, lat, h, cs.Lonf, cs.Eastf, cs.Northf, cs.k0())
}
func (cs MercatorA) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
return mercatorToGeographic(s, east, north, h, cs.Lonf, cs.Eastf, cs.Northf, cs.k0())
}
// MercatorB is EPSG method 9805 (PROJ +proj=merc with lat_ts / 2SP).
// Sp1 is the latitude of the first standard parallel; k0 is derived from it.
type MercatorB struct {
Lonf, Sp1, Eastf, Northf float64
}
func (cs MercatorB) String() string {
return build("mercator_b").addAll(
"lonf", cs.Lonf,
"sp1", cs.Sp1,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs MercatorB) k0(s Spheroid) float64 {
phi1 := radian(cs.Sp1)
sin1 := math.Sin(phi1)
return math.Cos(phi1) / math.Sqrt(1-s.E2()*sin1*sin1)
}
func (cs MercatorB) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
return mercatorFromGeographic(s, lon, lat, h, cs.Lonf, cs.Eastf, cs.Northf, cs.k0(s))
}
func (cs MercatorB) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
return mercatorToGeographic(s, east, north, h, cs.Lonf, cs.Eastf, cs.Northf, cs.k0(s))
}
func mercatorFromGeographic(s Spheroid, lon, lat, h, lonf, eastf, northf, k0 float64) (float64, float64, float64) {
e := s.E()
phi := radian(lat)
sinPhi := math.Sin(phi)
east := eastf + s.A*k0*(radian(lon)-radian(lonf))
north := northf + s.A*k0*math.Log(
math.Tan(math.Pi/4+phi/2)*math.Pow((1-e*sinPhi)/(1+e*sinPhi), e/2),
)
return east, north, h
}
func mercatorToGeographic(s Spheroid, east, north, h, lonf, eastf, northf, k0 float64) (float64, float64, float64) {
e2 := s.E2()
e4 := e2 * e2
e6 := e4 * e2
e8 := e4 * e4
t := math.Exp((northf - north) / (s.A * k0))
chi := math.Pi/2 - 2*math.Atan(t)
phi := chi +
(e2/2+5*e4/24+e6/12+13*e8/360)*math.Sin(2*chi) +
(7*e4/48+29*e6/240+811*e8/11520)*math.Sin(4*chi) +
(7*e6/120+81*e8/1120)*math.Sin(6*chi) +
(4279*e8/161280)*math.Sin(8*chi)
lambda := radian(lonf) + (east-eastf)/(s.A*k0)
return degree(lambda), degree(phi), h
}
package crs
import (
"math"
)
type ModifiedAzimuthalEquidistant struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs ModifiedAzimuthalEquidistant) String() string {
return build("modified_azimuthal_equidistant").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs ModifiedAzimuthalEquidistant) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
a, e, e2 := s.A, s.E(), s.E2()
phi0 := radian(cs.Latf)
lam0 := radian(cs.Lonf)
phi := radian(lat)
dLam := radian(lon) - lam0
sinPhi0, cosPhi0 := math.Sincos(phi0)
sinPhi, cosPhi := math.Sincos(phi)
nu0 := a / math.Sqrt(1-e2*sinPhi0*sinPhi0)
nu := a / math.Sqrt(1-e2*sinPhi*sinPhi)
psi := math.Atan((1-e2)*math.Tan(phi) + e2*nu0*sinPhi0/(nu*cosPhi))
alpha := math.Atan2(math.Sin(dLam), cosPhi0*math.Tan(psi)-sinPhi0*math.Cos(dLam))
sinAlpha, cosAlpha := math.Sincos(alpha)
g := e * sinPhi0 / math.Sqrt(1-e2)
hCoeff := e * cosPhi0 * cosAlpha / math.Sqrt(1-e2)
var sAng float64
if math.Abs(sinAlpha) < 1e-14 {
sAng = math.Asin(clamp(cosPhi0*math.Sin(psi)-sinPhi0*math.Cos(psi), -1, 1))
if cosAlpha < 0 {
sAng = -sAng
}
} else {
sAng = math.Asin(clamp(math.Sin(dLam)*math.Cos(psi)/sinAlpha, -1, 1))
}
s2 := sAng * sAng
s3 := s2 * sAng
s4 := s2 * s2
s5 := s4 * sAng
h2 := hCoeff * hCoeff
c := nu0 * sAng * (1 - s2*h2*(1-h2)/6 + s3/8*g*hCoeff*(1-2*h2) +
s4/120*(h2*(4-7*h2)-3*g*g*(1-7*h2)) - s5/48*g*hCoeff)
return cs.Eastf + c*sinAlpha, cs.Northf + c*cosAlpha, h
}
func (cs ModifiedAzimuthalEquidistant) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
a, e2 := s.A, s.E2()
phi0 := radian(cs.Latf)
sinPhi0, cosPhi0 := math.Sincos(phi0)
nu0 := a / math.Sqrt(1-e2*sinPhi0*sinPhi0)
x := east - cs.Eastf
y := north - cs.Northf
c := math.Hypot(x, y)
if c < 1e-14 {
return cs.Lonf, cs.Latf, h
}
alpha := math.Atan2(x, y)
sinAlpha, cosAlpha := math.Sincos(alpha)
aa := -e2 * cosPhi0 * cosPhi0 * cosAlpha * cosAlpha / (1 - e2)
b := 3 * e2 * (1 - aa) * sinPhi0 * cosPhi0 * cosAlpha / (1 - e2)
d := c / nu0
j := d - aa*(1+aa)*d*d*d/6 - b*(1+3*aa)*d*d*d*d/24
k := 1 - aa*j*j/2 - b*j*j*j/6
psi := math.Asin(clamp(sinPhi0*math.Cos(j)+cosPhi0*math.Sin(j)*cosAlpha, -1, 1))
phi := math.Atan((1 - e2*k*sinPhi0/math.Sin(psi)) * math.Tan(psi) / (1 - e2))
lam := radian(cs.Lonf) + math.Asin(clamp(sinAlpha*math.Sin(j)/math.Cos(psi), -1, 1))
return degree(lam), degree(phi), h
}
type GuamProjection struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs GuamProjection) String() string {
return build("guam_projection").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs GuamProjection) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
phi := radian(lat)
dLam := radian(lon) - radian(cs.Lonf)
sinPhi, cosPhi := math.Sincos(phi)
t := math.Sqrt(1 - s.E2()*sinPhi*sinPhi)
x := s.A * dLam * cosPhi / t
m0 := meridianDistance(s, radian(cs.Latf))
m := meridianDistance(s, phi)
east := cs.Eastf + x
north := cs.Northf + m - m0 + 0.5*x*x*math.Tan(phi)*t/s.A
return east, north, h
}
func (cs GuamProjection) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
m0 := meridianDistance(s, radian(cs.Latf))
x := east - cs.Eastf
x2 := 0.5 * x * x
phi := radian(cs.Latf)
for range 3 {
sinPhi := math.Sin(phi)
t := math.Sqrt(1 - s.E2()*sinPhi*sinPhi)
phi = footpointLatitude(s, m0+(north-cs.Northf)-x2*math.Tan(phi)*t/s.A)
}
sinPhi, cosPhi := math.Sincos(phi)
t := math.Sqrt(1 - s.E2()*sinPhi*sinPhi)
lam := radian(cs.Lonf) + x*t/(s.A*cosPhi)
return degree(lam), degree(phi), h
}
package crs
type MolodenskyBadekas struct {
Tx, Ty, Tz, Rx, Ry, Rz, Ds float64
Px, Py, Pz float64
}
func (m MolodenskyBadekas) String() string {
return build("molodensky_badekas").addAll(
"tx", m.Tx,
"ty", m.Ty,
"tz", m.Tz,
"rx", m.Rx,
"ry", m.Ry,
"rz", m.Rz,
"ds", m.Ds,
"px", m.Px,
"py", m.Py,
"pz", m.Pz,
).String()
}
func (m MolodenskyBadekas) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
x, y, z := source.GeographicToGeocentric(lon, lat, h)
x0, y0, z0 := calcMolodenskyBadekas(x, y, z, m.Tx, m.Ty, m.Tz, m.Rx, m.Ry, m.Rz, m.Ds, m.Px, m.Py, m.Pz)
lon0, lat0, h0 := target.GeocentricToGeographic(x0, y0, z0)
return lon0, lat0, h0, nil
}
func (m MolodenskyBadekas) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
x0, y0, z0 := source.GeographicToGeocentric(lon0, lat0, h0)
x, y, z := calcMolodenskyBadekasInverse(x0, y0, z0, m.Tx, m.Ty, m.Tz, m.Rx, m.Ry, m.Rz, m.Ds, m.Px, m.Py, m.Pz)
lon, lat, h := target.GeocentricToGeographic(x, y, z)
return lon, lat, h, nil
}
func calcMolodenskyBadekas(x, y, z, tx, ty, tz, rx, ry, rz, ds, px, py, pz float64) (x1, y1, z1 float64) {
x1, y1, z1 = calcHelmert(x-px, y-py, z-pz, tx, ty, tz, rx, ry, rz, ds)
return x1 + px, y1 + py, z1 + pz
}
func calcMolodenskyBadekasInverse(x1, y1, z1, tx, ty, tz, rx, ry, rz, ds, px, py, pz float64) (x, y, z float64) {
x, y, z = calcHelmert(x1-tx-px, y1-ty-py, z1-tz-pz, 0, 0, 0, -rx, -ry, -rz, -ds)
return x + px, y + py, z + pz
}
package crs
import (
"bytes"
"fmt"
"io"
"io/fs"
"math"
"net/http"
"os"
"path"
"path/filepath"
"slices"
"sort"
"strings"
"sync"
"time"
)
type gridFilesystem struct {
prefix string
fsys fs.FS
}
var (
gridMu sync.RWMutex
gridCDNBase string
gridCDNCache string
gridSearchDirs []string
gridFS []gridFilesystem // sorted by prefix length, longest first
gridHTTPClient = &http.Client{Timeout: 10 * time.Minute}
)
func SetGridCDN(uri string, cacheDir ...string) {
gridMu.Lock()
defer gridMu.Unlock()
uri = strings.TrimSpace(uri)
if uri == "" {
gridCDNBase = ""
gridCDNCache = ""
return
}
if !strings.HasSuffix(uri, "/") {
uri += "/"
}
gridCDNBase = uri
gridCDNCache = ""
if len(cacheDir) > 0 {
gridCDNCache = strings.TrimSpace(cacheDir[0])
}
}
func RegisterGridDir(dir string) {
dir = strings.TrimSpace(dir)
if dir == "" {
return
}
gridMu.Lock()
gridSearchDirs = append(gridSearchDirs, dir)
gridMu.Unlock()
}
func RegisterGridFS(prefix string, fsys fs.FS) {
if fsys == nil {
panic("crs: RegisterGridFS: fsys is nil")
}
if prefix == "/" {
prefix = ""
}
gridMu.Lock()
gridFS = append(gridFS, gridFilesystem{prefix: prefix, fsys: fsys})
sort.Slice(gridFS, func(i, j int) bool {
return len(gridFS[i].prefix) > len(gridFS[j].prefix)
})
gridMu.Unlock()
}
func gridFilename(p string) string {
return path.Base(strings.ReplaceAll(p, `\`, `/`))
}
func openGridReader(name string) (io.ReadCloser, error) {
gridMu.RLock()
cdnBase := gridCDNBase
cacheDir := gridCDNCache
searchDirs := append([]string(nil), gridSearchDirs...)
filesystems := append([]gridFilesystem(nil), gridFS...)
gridMu.RUnlock()
if cacheDir != "" {
found := slices.Contains(searchDirs, cacheDir)
if !found {
searchDirs = append(searchDirs, cacheDir)
}
}
for _, dir := range searchDirs {
if dir == "" {
continue
}
if f, err := os.Open(filepath.Join(dir, name)); err == nil {
return f, nil
}
}
for _, reg := range filesystems {
p := name
if reg.prefix != "" {
p = strings.TrimSuffix(reg.prefix, "/") + "/" + name
}
if f, err := reg.fsys.Open(p); err == nil {
return f, nil
}
}
if cdnBase == "" {
return nil, fmt.Errorf("crs: grid %q not found", name)
}
resp, err := gridHTTPClient.Get(cdnBase + name)
if err != nil {
return nil, fmt.Errorf("crs: grid %q not found", name)
}
if resp.StatusCode == http.StatusNotFound {
_ = resp.Body.Close()
return nil, fmt.Errorf("crs: grid %q not found", name)
}
if resp.StatusCode != http.StatusOK {
_ = resp.Body.Close()
return nil, fmt.Errorf("cdn: grid %q: HTTP %s", name, resp.Status)
}
if cacheDir == "" {
return resp.Body, nil
}
data, err := io.ReadAll(resp.Body)
_ = resp.Body.Close()
if err != nil {
return nil, err
}
final := filepath.Join(cacheDir, name)
tmp := final + ".tmp"
if err = os.MkdirAll(cacheDir, 0o755); err != nil {
return io.NopCloser(bytes.NewReader(data)), nil
}
if err = os.WriteFile(tmp, data, 0o644); err != nil {
return io.NopCloser(bytes.NewReader(data)), nil
}
if err = os.Rename(tmp, final); err != nil {
_ = os.Remove(tmp)
return io.NopCloser(bytes.NewReader(data)), nil
}
if f, err := os.Open(final); err == nil {
return f, nil
}
return io.NopCloser(bytes.NewReader(data)), nil
}
// gridStore caches parsed grids in memory (local embed, disk, and CDN sources).
var gridStore sync.Map
func loadGrid(name string) (grid *horizontalGridData, err error) {
name = gridFilename(name)
if v, ok := gridStore.Load(name); ok {
return v.(*horizontalGridData), nil
}
defer func() {
if err == nil {
grid.File = name
gridStore.Store(name, grid)
}
}()
r, err := openGridReader(name)
if err != nil {
return nil, GridNotFoundError{
Err: fmt.Errorf("grid: %s", name),
}
}
defer r.Close() //nolint:errcheck
switch strings.ToLower(path.Ext(name)) {
case ".tif":
return parseTaggedImageFileFormat(r)
case ".gsb":
return parseGridShiftBinary(r)
default:
return nil, fmt.Errorf("crs: unknown grid extension %q", path.Ext(name))
}
}
type subGrid struct {
Name string
Parent string
Columns int
Rows int
SLat float64
NLat float64
ELong float64
WLong float64
LatInc float64
LongInc float64
Values [][2]float32
}
func (g *horizontalGridData) ToWGS84(lon, lat float64) (float64, float64, error) {
dlon, dlat, err := g.Shift(lon, lat)
if err != nil {
return 0, 0, err
}
return lon + dlon, lat + dlat, nil
}
func (g *horizontalGridData) FromWGS84(lon, lat float64) (float64, float64, error) {
qlon, qlat := lon, lat
for range 10 {
dlon, dlat, err := g.Shift(qlon, qlat)
if err != nil {
return 0, 0, err
}
newLon := lon - dlon
newLat := lat - dlat
if math.Abs(newLon-qlon) < 1e-12 && math.Abs(newLat-qlat) < 1e-12 {
break
}
qlon, qlat = newLon, newLat
}
return qlon, qlat, nil
}
func (g *horizontalGridData) Shift(lon, lat float64) (dlon, dlat float64, err error) {
idx := g.selectSubgrid(lon, lat)
if idx < 0 {
return 0, 0, OutOfBoundsError{
Err: fmt.Errorf("coordinate: [%f, %f]: %s", lon, lat, g.File),
}
}
return g.subGrids[idx].shift(lon, lat)
}
func (g *horizontalGridData) selectSubgrid(lon, lat float64) int {
lam := -lon * 3600
phi := lat * 3600
for i, sg := range g.subGrids {
if sg.Parent != "" && sg.Parent != "NONE" {
continue
}
if sg.contains(phi, lam) {
return g.deepestSubgridRec(i, phi, lam, make(map[int]bool))
}
}
return -1
}
func (g *horizontalGridData) deepestSubgridRec(idx int, phi, lam float64, seen map[int]bool) int {
if seen[idx] {
return idx
}
seen[idx] = true
best := idx
name := g.subGrids[idx].Name
for i, sg := range g.subGrids {
if sg.Parent != name || !sg.contains(phi, lam) {
continue
}
if deeper := g.deepestSubgridRec(i, phi, lam, seen); deeper >= 0 {
best = deeper
}
}
return best
}
const relToleranceHGridShift = 1e-5
func (sg subGrid) gridEpsilon() float64 {
return (sg.LongInc + sg.LatInc) * relToleranceHGridShift
}
func (sg subGrid) contains(phi, lam float64) bool {
eps := sg.gridEpsilon()
return phi >= sg.SLat-eps && phi <= sg.NLat+eps &&
lam >= sg.ELong-eps && lam <= sg.WLong+eps
}
func interpolateIndex(f float64, size int) (idx int, frac float64, ok bool) {
if math.IsNaN(f) {
return 0, 0, false
}
idx = int(math.Round(math.Floor(f)))
frac = f - float64(idx)
if idx < 0 {
if idx == -1 && frac > 1-10*relToleranceHGridShift {
idx++
frac = 0
} else {
return 0, 0, false
}
} else if idx+1 >= size {
if idx+1 == size && frac < 10*relToleranceHGridShift {
idx--
frac = 1
} else {
return 0, 0, false
}
}
return idx, frac, true
}
func (sg subGrid) shift(lon, lat float64) (dlon, dlat float64, err error) {
if sg.Columns < 2 || sg.Rows < 2 || len(sg.Values) == 0 {
return 0, 0, OutOfBoundsError{
Err: fmt.Errorf("coordinate: [%f, %f]: %s", lon, lat, sg.Name),
}
}
phi := lat * 3600
lam := -lon * 3600
if !sg.contains(phi, lam) {
return 0, 0, OutOfBoundsError{
Err: fmt.Errorf("coordinate: [%f, %f]: %s", lon, lat, sg.Name),
}
}
fcol := (lam - sg.ELong) / sg.LongInc
frow := (phi - sg.SLat) / sg.LatInc
ppr := sg.Columns
col, dx, ok := interpolateIndex(fcol, ppr)
if !ok {
return 0, 0, OutOfBoundsError{
Err: fmt.Errorf("coordinate: [%f, %f]: %s", lon, lat, sg.Name),
}
}
row, dy, ok := interpolateIndex(frow, sg.Rows)
if !ok {
return 0, 0, OutOfBoundsError{
Err: fmt.Errorf("coordinate: [%f, %f]: %s", lon, lat, sg.Name),
}
}
se := row*ppr + col
sw := se + 1
ne := se + ppr
nw := ne + 1
sse := sg.Values[se]
ssw := sg.Values[sw]
sne := sg.Values[ne]
snw := sg.Values[nw]
latsv := (1-dx)*(1-dy)*float64(sse[0]) + dx*(1-dy)*float64(ssw[0]) +
(1-dx)*dy*float64(sne[0]) + dx*dy*float64(snw[0])
lonsv := (1-dx)*(1-dy)*float64(sse[1]) + dx*(1-dy)*float64(ssw[1]) +
(1-dx)*dy*float64(sne[1]) + dx*dy*float64(snw[1])
return -lonsv / 3600, latsv / 3600, nil
}
func (sg subGrid) validate() error {
want := sg.Columns * sg.Rows
if want == 0 {
return fmt.Errorf("subgrid %q: zero grid dimensions", sg.Name)
}
if len(sg.Values) != want {
return fmt.Errorf("subgrid %q: got %d values, want %d", sg.Name, len(sg.Values), want)
}
return nil
}
package crs
import (
"encoding/binary"
"fmt"
"io"
"math"
"strings"
)
func parseGridShiftBinary(r io.Reader) (*horizontalGridData, error) {
first, err := readGridRecord(r)
if err != nil {
return nil, fmt.Errorf("read file header: %w", err)
}
order := detectGSBEndian(first)
numSrec, numFile, err := parseGridFileHeader(first, r, order)
if err != nil {
return nil, fmt.Errorf("read file header: %w", err)
}
n := []subGrid{}
for range numFile {
sg, count, err := readGridSubgridHeader(r, numSrec, order)
if err != nil {
return nil, fmt.Errorf("read subgrid header: %w", err)
}
sg.Values = make([][2]float32, count)
for i := range sg.Values {
rec, err := readGridRecord(r)
if err != nil {
return nil, fmt.Errorf("read grid value %d: %w", i, err)
}
sg.Values[i][0] = recordFloat32(rec, 0, order)
sg.Values[i][1] = recordFloat32(rec, 4, order)
}
if err := sg.validate(); err != nil {
return nil, err
}
n = append(n, sg)
}
return &horizontalGridData{
subGrids: n,
}, nil
}
type gridRecord [16]byte
func detectGSBEndian(rec gridRecord) binary.ByteOrder {
if recordName(rec) != "NUM_OREC" {
return binary.LittleEndian
}
le := int32(binary.LittleEndian.Uint32(rec[8:12]))
if le > 0 && le < 100 {
return binary.LittleEndian
}
be := int32(binary.BigEndian.Uint32(rec[8:12]))
if be > 0 && be < 100 {
return binary.BigEndian
}
return binary.LittleEndian
}
func parseGridFileHeader(first gridRecord, r io.Reader, order binary.ByteOrder) (numSrec, numFile int, err error) {
rec := first
numOrec := int(recordValueInt32(rec, order))
if recordName(rec) == "NUM_OREC" && numOrec > 0 {
for i := 1; i < numOrec; i++ {
rec, err = readGridRecord(r)
if err != nil {
return 0, 0, err
}
switch recordName(rec) {
case "NUM_SREC":
numSrec = int(recordValueInt32(rec, order))
case "NUM_FILE":
numFile = int(recordValueInt32(rec, order))
}
}
}
if numSrec <= 0 {
numSrec = 11
}
return numSrec, numFile, nil
}
func readGridSubgridHeader(r io.Reader, numSrec int, order binary.ByteOrder) (subGrid, int, error) {
var (
sg subGrid
count int
cn [4]string
)
for range numSrec {
rec, err := readGridRecord(r)
if err != nil {
return subGrid{}, 0, err
}
switch recordName(rec) {
case "SUB_NAME":
sg.Name = recordValueString(rec)
case "CN_1":
cn[0] = recordValueString(rec)
case "CN_2":
cn[1] = recordValueString(rec)
case "CN_3":
cn[2] = recordValueString(rec)
case "CN_4":
cn[3] = recordValueString(rec)
case "PARENT":
sg.Parent = recordValueString(rec)
case "S_LAT":
sg.SLat = recordValueFloat64(rec, order)
case "N_LAT":
sg.NLat = recordValueFloat64(rec, order)
case "E_LONG":
sg.ELong = recordValueFloat64(rec, order)
case "W_LONG":
sg.WLong = recordValueFloat64(rec, order)
case "LAT_INC":
sg.LatInc = recordValueFloat64(rec, order)
case "LONG_INC":
sg.LongInc = recordValueFloat64(rec, order)
case "COLUMNS":
sg.Columns = int(recordValueInt32(rec, order))
case "ROWS":
sg.Rows = int(recordValueInt32(rec, order))
case "GS_COUNT":
count = int(recordValueInt32(rec, order))
}
}
if sg.Name == "" {
sg.Name = trimGrid(strings.Join(cn[:], ""))
}
if count <= 0 {
return subGrid{}, 0, fmt.Errorf("subgrid %q: missing GS_COUNT", sg.Name)
}
if sg.Columns == 0 || sg.Rows == 0 {
if sg.LatInc == 0 || sg.LongInc == 0 {
return subGrid{}, 0, fmt.Errorf("subgrid %q: missing grid dimensions", sg.Name)
}
sg.Rows = int(math.Floor((sg.NLat-sg.SLat)/sg.LatInc+0.5)) + 1
sg.Columns = int(math.Floor((sg.WLong-sg.ELong)/sg.LongInc+0.5)) + 1
}
if sg.Columns*sg.Rows != count {
return subGrid{}, 0, fmt.Errorf("subgrid %q: GS_COUNT %d != %d×%d", sg.Name, count, sg.Columns, sg.Rows)
}
return sg, count, nil
}
func trimGrid(s string) string {
return strings.TrimFunc(s, func(r rune) bool {
return r == 0 || r == ' '
})
}
func recordValueString(rec gridRecord) string {
return trimGrid(string(rec[8:16]))
}
func recordValueInt32(rec gridRecord, order binary.ByteOrder) int32 {
return int32(order.Uint32(rec[8:12]))
}
func recordValueFloat64(rec gridRecord, order binary.ByteOrder) float64 {
return math.Float64frombits(order.Uint64(rec[8:16]))
}
func readGridRecord(r io.Reader) (gridRecord, error) {
var rec gridRecord
_, err := io.ReadFull(r, rec[:])
return rec, err
}
func recordName(rec gridRecord) string {
return trimGrid(string(rec[:8]))
}
func recordFloat32(rec gridRecord, off int, order binary.ByteOrder) float32 {
return math.Float32frombits(order.Uint32(rec[off : off+4]))
}
package crs
import (
"bytes"
"compress/zlib"
"encoding/binary"
"fmt"
"io"
"math"
"regexp"
"strings"
)
func parseTaggedImageFileFormat(r io.Reader) (*horizontalGridData, error) {
data, err := io.ReadAll(r)
if err != nil {
return nil, err
}
p := &tifParser{data: data}
if err := p.readHeader(); err != nil {
return nil, err
}
var subgrids []subGrid
for ifdOff := p.firstIFD; ifdOff != 0; {
sg, next, err := p.readIFD(ifdOff)
if err != nil {
return nil, err
}
if err := sg.validate(); err != nil {
return nil, err
}
subgrids = append(subgrids, sg)
ifdOff = next
}
if len(subgrids) == 0 {
return nil, fmt.Errorf("tiff: no image directories")
}
return &horizontalGridData{subGrids: subgrids}, nil
}
type tifParser struct {
data []byte
order binary.ByteOrder
firstIFD uint32
}
func (p *tifParser) readHeader() error {
if len(p.data) < 8 {
return fmt.Errorf("tiff: file too short")
}
switch string(p.data[:2]) {
case "II":
p.order = binary.LittleEndian
case "MM":
p.order = binary.BigEndian
default:
return fmt.Errorf("tiff: invalid byte order")
}
if p.u16(2) != 42 {
return fmt.Errorf("tiff: invalid magic")
}
p.firstIFD = p.u32(4)
return nil
}
func (p *tifParser) readIFD(off uint32) (subGrid, uint32, error) {
if int(off)+2 > len(p.data) {
return subGrid{}, 0, fmt.Errorf("tiff: IFD out of range")
}
n := int(p.u16(int(off)))
base := int(off) + 2
need := base + n*12 + 4
if need > len(p.data) {
return subGrid{}, 0, fmt.Errorf("tiff: IFD truncated")
}
tags := map[uint16]ifdEntry{}
for i := range n {
eo := base + i*12
e := ifdEntry{
tag: p.u16(eo),
typ: p.u16(eo + 2),
count: p.u32(eo + 4),
val: p.u32(eo + 8),
}
tags[e.tag] = e
}
next := p.u32(base + n*12)
width := int(p.tagLong(tags, 256))
height := int(p.tagLong(tags, 257))
if width <= 0 || height <= 0 {
return subGrid{}, next, fmt.Errorf("tiff: invalid dimensions %dx%d", width, height)
}
meta := parseGDALMetadata(p.tagString(tags, 42112))
if typ := meta["TYPE"]; typ != "" && typ != "HORIZONTAL_OFFSET" {
return subGrid{}, next, fmt.Errorf("tiff: unsupported TYPE %q", typ)
}
samples := int(p.tagLong(tags, 277))
if samples < 2 {
return subGrid{}, next, fmt.Errorf("tiff: need at least 2 samples, got %d", samples)
}
bands, err := p.readBands(tags, width, height, samples)
if err != nil {
return subGrid{}, next, err
}
latBand := bands[0]
lonBand := bands[1]
posLon := meta["positive_value"]
if posLon == "" {
posLon = "east"
}
tie := p.tagDoubles(tags, 33922, 6)
scale := p.tagDoubles(tags, 33550, 3)
if len(tie) < 6 || len(scale) < 2 {
return subGrid{}, next, fmt.Errorf("tiff: missing ModelTiepointTag/ModelPixelScaleTag")
}
lonUL := tie[3] // degrees east
latUL := tie[4] // degrees north
dLon := scale[0]
dLat := math.Abs(scale[1])
if dLon == 0 || dLat == 0 {
return subGrid{}, next, fmt.Errorf("tiff: zero pixel scale")
}
lonWest := lonUL
lonEast := lonUL + float64(width-1)*dLon
latNorth := latUL
latSouth := latUL - float64(height-1)*dLat
sg := subGrid{
Name: meta["grid_name"],
Parent: meta["parent_grid_name"],
Columns: width,
Rows: height,
SLat: latSouth * 3600,
NLat: latNorth * 3600,
ELong: (-lonEast) * 3600, // eastern edge (min lam = -lon×3600)
WLong: (-lonWest) * 3600, // western edge (max lam)
LatInc: math.Abs(dLat) * 3600,
LongInc: math.Abs(dLon) * 3600,
Values: make([][2]float32, width*height),
}
for row := range height {
tiffRow := height - 1 - row
for col := range width {
iTiff := tiffRow*width + col
iNtv2 := row*width + (width - 1 - col)
latOff := latBand[iTiff]
lonOff := lonBand[iTiff]
switch posLon {
case "east":
lonOff = -lonOff // match Grid west-positive storage
case "west":
// already west-positive
default:
return subGrid{}, next, fmt.Errorf("tiff: unknown positive_value %q", posLon)
}
sg.Values[iNtv2] = [2]float32{latOff, lonOff}
}
}
return sg, next, nil
}
type ifdEntry struct {
tag, typ uint16
count uint32
val uint32
}
func (p *tifParser) u16(off int) uint16 {
return p.order.Uint16(p.data[off:])
}
func (p *tifParser) u32(off int) uint32 {
return p.order.Uint32(p.data[off:])
}
func (p *tifParser) ifdValueSize(typ uint16) int {
switch typ {
case 3, 8: // SHORT, SBYTE
return 2
case 4, 9, 11: // LONG, SLONG, FLOAT
return 4
default:
return 0
}
}
func (p *tifParser) ifdInline(e ifdEntry) bool {
sz := p.ifdValueSize(e.typ)
return sz > 0 && int(e.count)*sz <= 4
}
func (p *tifParser) ifdInlineU16(e ifdEntry, i int) uint16 {
var b [4]byte
p.order.PutUint32(b[:], e.val)
return p.order.Uint16(b[i*2 : i*2+2])
}
func (p *tifParser) tagLong(tags map[uint16]ifdEntry, tag uint16) uint32 {
e, ok := tags[tag]
if !ok {
return 0
}
switch e.typ {
case 3: // SHORT
if e.count == 1 {
return uint32(uint16(e.val))
}
if p.ifdInline(e) {
return uint32(p.ifdInlineU16(e, 0))
}
return uint32(p.u16(int(e.val)))
case 4: // LONG
if e.count == 1 {
return e.val
}
if p.ifdInline(e) {
return e.val
}
return p.u32(int(e.val))
default:
return 0
}
}
func (p *tifParser) tagString(tags map[uint16]ifdEntry, tag uint16) string {
e, ok := tags[tag]
if !ok || e.typ != 2 {
return ""
}
off := int(e.val)
end := off + int(e.count)
if off < 0 || end > len(p.data) {
return ""
}
return strings.TrimRight(string(p.data[off:end]), "\x00")
}
func (p *tifParser) tagDoubles(tags map[uint16]ifdEntry, tag uint16, want int) []float64 {
e, ok := tags[tag]
if !ok {
return nil
}
off := int(e.val)
out := make([]float64, 0, want)
switch e.typ {
case 5: // RATIONAL
for i := 0; i < int(e.count) && i < want; i++ {
o := off + i*8
num := float64(p.order.Uint32(p.data[o:]))
den := float64(p.order.Uint32(p.data[o+4:]))
if den == 0 {
return nil
}
out = append(out, num/den)
}
case 12: // DOUBLE
for i := 0; i < int(e.count) && i < want; i++ {
o := off + i*8
out = append(out, math.Float64frombits(p.order.Uint64(p.data[o:])))
}
default:
return nil
}
return out
}
func (p *tifParser) readBands(tags map[uint16]ifdEntry, width, height, samples int) ([][]float32, error) {
planar := int(p.tagLong(tags, 284)) // 1=contig, 2=separate
if planar == 0 {
planar = 1
}
comp := int(p.tagLong(tags, 259))
predictor := int(p.tagLong(tags, 317))
bps := int(p.tagLong(tags, 258))
if bps == 0 {
bps = 32
}
sfmt := int(p.tagLong(tags, 339)) // SampleFormat (338 is ExtraSamples)
if sfmt == 0 && bps == 32 {
sfmt = 3 // IEEE float
}
if bps != 32 || sfmt != 3 {
return nil, fmt.Errorf("tiff: only Float32 supported (bps=%d fmt=%d)", bps, sfmt)
}
n := width * height
bands := make([][]float32, samples)
for i := range bands {
bands[i] = make([]float32, n)
}
if p.isTiled(tags) {
return p.readTiledBands(tags, width, height, samples, comp, predictor, planar, bands)
}
return p.readStripBands(tags, width, height, samples, comp, predictor, planar, bands)
}
func (p *tifParser) isTiled(tags map[uint16]ifdEntry) bool {
if len(p.tagLongArray(tags, 273)) > 0 {
return false
}
if len(p.tagLongArray(tags, 324)) > 1 {
return true
}
if len(p.tagLongArray(tags, 322)) > 1 {
return true
}
return false
}
func (p *tifParser) tileLayout(tags map[uint16]ifdEntry) (offsets, counts []uint32, tileW, tileH int) {
// Standard TIFF: 322/323 offsets, 324/325 dimensions.
offsets = p.tagLongArray(tags, 322)
counts = p.tagLongArray(tags, 323)
tileW = int(p.tagLong(tags, 324))
tileH = int(p.tagLong(tags, 325))
// PROJ GTG: 322/323 hold tile size, 324/325 hold offset tables.
if len(offsets) <= 1 && len(p.tagLongArray(tags, 324)) > 1 {
offsets = p.tagLongArray(tags, 324)
counts = p.tagLongArray(tags, 325)
tileW = int(p.tagLong(tags, 322))
tileH = int(p.tagLong(tags, 323))
}
return offsets, counts, tileW, tileH
}
func (p *tifParser) readTiledBands(
tags map[uint16]ifdEntry,
width, height, samples, comp, predictor, planar int,
bands [][]float32,
) ([][]float32, error) {
offsets, counts, tileW, tileH := p.tileLayout(tags)
if len(offsets) == 0 || len(offsets) != len(counts) || tileW <= 0 || tileH <= 0 {
return nil, fmt.Errorf("tiff: invalid tile tables")
}
tilesX := (width + tileW - 1) / tileW
tilesY := (height + tileH - 1) / tileH
tilesPerSample := tilesX * tilesY
if planar == 2 {
if len(offsets)%samples != 0 || len(offsets)/samples != tilesPerSample {
return nil, fmt.Errorf("tiff: tile count %d, want %d×%d", len(offsets), samples, tilesPerSample)
}
} else if len(offsets) != tilesPerSample {
return nil, fmt.Errorf("tiff: tile count %d, want %d", len(offsets), tilesPerSample)
}
readSampleTiles := func(sample int) error {
for ty := range tilesY {
for tx := range tilesX {
tw := min(tileW, width-tx*tileW)
th := min(tileH, height-ty*tileH)
ti := ty*tilesX + tx
idx := ti
if planar == 2 {
idx = sample*tilesPerSample + ti
}
dec, err := p.readChunk(offsets[idx], counts[idx], comp, tileW*tileH*4, predictor, tileW)
if err != nil {
return fmt.Errorf("tile (%d,%d) sample %d: %w", tx, ty, sample, err)
}
for row := range th {
for col := range tw {
src := (row*tileW + col) * 4
dst := (ty*tileH+row)*width + (tx*tileW + col)
bands[sample][dst] = math.Float32frombits(p.order.Uint32(dec[src:]))
}
}
}
}
return nil
}
if planar == 2 {
for s := range samples {
if err := readSampleTiles(s); err != nil {
return nil, err
}
}
return bands, nil
}
if err := readSampleTiles(0); err != nil {
return nil, err
}
if samples > 1 {
return nil, fmt.Errorf("tiff: contig tiled multi-sample layout not supported")
}
return bands, nil
}
func (p *tifParser) readStripBands(
tags map[uint16]ifdEntry,
width, height, samples, comp, predictor, planar int,
bands [][]float32,
) ([][]float32, error) {
offsets := p.tagLongArray(tags, 273)
counts := p.tagLongArray(tags, 279)
if len(offsets) == 0 || len(offsets) != len(counts) {
return nil, fmt.Errorf("tiff: invalid strip tables")
}
n := width * height
rowsPerStrip := int(p.tagLong(tags, 278))
if rowsPerStrip == 0 {
rowsPerStrip = height
}
stripsPerSample := len(offsets)
if planar == 2 {
if len(offsets)%samples != 0 {
return nil, fmt.Errorf("tiff: strip count %d not divisible by samples %d", len(offsets), samples)
}
stripsPerSample = len(offsets) / samples
}
readSample := func(sample int) error {
stripBase := 0
if planar == 2 {
stripBase = sample * stripsPerSample
}
row := 0
for si := range stripsPerSample {
idx := stripBase + si
stripRows := rowsPerStrip
if row+stripRows > height {
stripRows = height - row
}
want := width * stripRows * 4
dec, err := p.readChunk(offsets[idx], counts[idx], comp, want, predictor, width)
if err != nil {
return fmt.Errorf("strip %d sample %d: %w", si, sample, err)
}
for r := range stripRows {
for col := range width {
src := (r*width + col) * 4
dst := (row+r)*width + col
bands[sample][dst] = math.Float32frombits(p.order.Uint32(dec[src:]))
}
}
row += stripRows
}
return nil
}
if planar == 2 {
for s := range samples {
if err := readSample(s); err != nil {
return nil, fmt.Errorf("sample %d: %w", s, err)
}
}
return bands, nil
}
dec, err := p.decompressStrips(offsets, counts, comp, n*samples*4, predictor, width*samples)
if err != nil {
return nil, err
}
for i := range n {
for s := range samples {
off := (i*samples + s) * 4
bands[s][i] = math.Float32frombits(p.order.Uint32(dec[off:]))
}
}
return bands, nil
}
func (p *tifParser) readChunk(offset, count uint32, comp, want, predictor, rowWidth int) ([]byte, error) {
off := int(offset)
n := int(count)
if off < 0 || off+n > len(p.data) {
return nil, fmt.Errorf("tiff: chunk out of range")
}
dec, err := decompressTIFF(p.data[off:off+n], comp)
if err != nil {
return nil, err
}
if predictor == 3 {
rowBytes := rowWidth * 4
for rowOff := 0; rowOff+rowBytes <= len(dec); rowOff += rowBytes {
undoFloatPredictorRow(dec[rowOff : rowOff+rowBytes])
}
}
if want >= 0 && len(dec) != want {
return nil, fmt.Errorf("tiff: chunk size %d, want %d", len(dec), want)
}
return dec, nil
}
func (p *tifParser) decompressStrips(offsets, counts []uint32, comp, want, predictor, rowWidth int) ([]byte, error) {
var out []byte
for i := range offsets {
dec, err := p.readChunk(offsets[i], counts[i], comp, -1, predictor, rowWidth)
if err != nil {
return nil, fmt.Errorf("strip %d: %w", i, err)
}
out = append(out, dec...)
}
if want >= 0 && len(out) != want {
return nil, fmt.Errorf("tiff: strip total size %d, want %d", len(out), want)
}
return out, nil
}
// undoFloatPredictorRow reverses libtiff predictor 3 for one scanline/tile row.
func undoFloatPredictorRow(row []byte) {
if len(row) < 8 {
return
}
const bps = 4
wc := len(row) / bps
for i := 1; i < len(row); i++ {
row[i] = row[i] + row[i-1]
}
tmp := make([]byte, len(row))
copy(tmp, row)
for i := range wc {
for b := range bps {
row[i*bps+b] = tmp[(bps-b-1)*wc+i]
}
}
}
func (p *tifParser) tagLongArray(tags map[uint16]ifdEntry, tag uint16) []uint32 {
e, ok := tags[tag]
if !ok {
return nil
}
if e.count == 1 {
switch e.typ {
case 3:
return []uint32{uint32(uint16(e.val))}
case 4:
return []uint32{e.val}
}
}
if p.ifdInline(e) {
out := make([]uint32, e.count)
for i := range e.count {
switch e.typ {
case 3:
out[i] = uint32(p.ifdInlineU16(e, int(i)))
case 4:
out[i] = e.val
default:
return nil
}
}
return out
}
off := int(e.val)
out := make([]uint32, e.count)
for i := range e.count {
switch e.typ {
case 4:
out[i] = p.u32(off + int(i)*4)
case 3:
out[i] = uint32(p.u16(off + int(i)*2))
default:
return nil
}
}
return out
}
func decompressTIFF(raw []byte, compression int) ([]byte, error) {
switch compression {
case 1: // none
return raw, nil
case 8, 32946: // deflate / adobe deflate
zr, err := zlib.NewReader(bytes.NewReader(raw))
if err != nil {
return nil, fmt.Errorf("tiff: deflate: %w", err)
}
defer zr.Close() //nolint:errcheck
return io.ReadAll(zr)
default:
return nil, fmt.Errorf("tiff: unsupported compression %d", compression)
}
}
var gdalItemRE = regexp.MustCompile(`<Item\s+name="([^"]+)"[^>]*>([^<]*)</Item>`)
func parseGDALMetadata(s string) map[string]string {
out := map[string]string{}
for _, m := range gdalItemRE.FindAllStringSubmatch(s, -1) {
out[m[1]] = strings.TrimSpace(m[2])
}
return out
}
package crs
import (
"math"
)
type NewZealandMapGrid struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs NewZealandMapGrid) String() string {
return build("new_zealand_map_grid").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
const (
nzmgSec5ToRad = 0.4848136811095359935899141023
nzmgRadToSec5 = 2.062648062470963551564733573
)
var nzmgBf = [...]complex128{
complex(0.7557853228, 0),
complex(0.249204646, 0.003371507),
complex(-0.001541739, 0.041058560),
complex(-0.10162907, 0.01727609),
complex(-0.26623489, -0.36249218),
complex(-0.6870983, -1.1651967),
}
var nzmgTpsi = [...]float64{
0.6399175073, -0.1358797613, 0.063294409, -0.02526853, 0.0117879,
-0.0055161, 0.0026906, -0.001333, 0.00067, -0.00034,
}
var nzmgTphi = [...]float64{
1.5627014243, 0.5185406398, -0.03333098, -0.1052906, -0.0368594,
0.007317, 0.01220, 0.00394, -0.0013,
}
func nzmgZpoly1(z complex128) complex128 {
a := nzmgBf[5]
for i := 4; i >= 0; i-- {
a = nzmgBf[i] + z*a
}
return z * a
}
func nzmgZpolyd1(z complex128) (f, fp complex128) {
a := nzmgBf[5]
b := a
for i := 4; i >= 0; i-- {
b = a + z*b
a = nzmgBf[i] + z*a
}
b = a + z*b
return z * a, b
}
func (cs NewZealandMapGrid) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
dPhi := (radian(lat) - radian(cs.Latf)) * nzmgRadToSec5
p := nzmgTpsi[9]
for i := 8; i >= 0; i-- {
p = nzmgTpsi[i] + dPhi*p
}
p *= dPhi
z := nzmgZpoly1(complex(p, radian(lon)-radian(cs.Lonf)))
return cs.Eastf + s.A*imag(z), cs.Northf + s.A*real(z), h
}
func (cs NewZealandMapGrid) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
z := complex((north-cs.Northf)/s.A, (east-cs.Eastf)/s.A)
p := z
for range 20 {
f, fp := nzmgZpolyd1(p)
f -= z
den := real(fp)*real(fp) + imag(fp)*imag(fp)
dp := complex(-(real(f)*real(fp)+imag(f)*imag(fp))/den, -(imag(f)*real(fp)-real(f)*imag(fp))/den)
p += dp
if math.Abs(real(dp))+math.Abs(imag(dp)) <= 1e-10 {
break
}
}
lam := imag(p)
phiNom := nzmgTphi[8]
for i := 7; i >= 0; i-- {
phiNom = nzmgTphi[i] + real(p)*phiNom
}
phi := radian(cs.Latf) + real(p)*phiNom*nzmgSec5ToRad
return degree(lam + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"math"
)
type ObliqueStereographic struct {
Lonf, Latf, Scale, Eastf, Northf float64
}
func (cs ObliqueStereographic) String() string {
return build("oblique_stereographic").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
type obliqueStereoConsts struct {
R, n, c, chi0, k0, lon0 float64
}
func (cs ObliqueStereographic) consts(s Spheroid) obliqueStereoConsts {
e := s.E()
e2 := s.E2()
k0 := cs.Scale
if k0 == 0 {
k0 = 1
}
phi0 := radian(cs.Latf)
lon0 := radian(cs.Lonf)
sin0 := math.Sin(phi0)
cos0 := math.Cos(phi0)
oneMe2Sin2 := 1 - e2*sin0*sin0
rho0 := s.A * (1 - e2) / math.Pow(oneMe2Sin2, 1.5)
nu0 := s.A / math.Sqrt(oneMe2Sin2)
R := math.Sqrt(rho0 * nu0)
n := math.Sqrt(1 + (e2*math.Pow(cos0, 4))/(1-e2))
S1 := (1 + sin0) / (1 - sin0)
S2 := (1 - e*sin0) / (1 + e*sin0)
w1 := math.Pow(S1*math.Pow(S2, e), n)
sinChi00 := (w1 - 1) / (w1 + 1)
c := ((n + sin0) * (1 - sinChi00)) / ((n - sin0) * (1 + sinChi00))
w2 := c * w1
chi0 := math.Asin((w2 - 1) / (w2 + 1))
return obliqueStereoConsts{R: R, n: n, c: c, chi0: chi0, k0: k0, lon0: lon0}
}
func (cs ObliqueStereographic) conformalLat(s Spheroid, phi, c, n float64) float64 {
e := s.E()
sinPhi := math.Sin(phi)
Sa := (1 + sinPhi) / (1 - sinPhi)
Sb := (1 - e*sinPhi) / (1 + e*sinPhi)
w := c * math.Pow(Sa*math.Pow(Sb, e), n)
return math.Asin((w - 1) / (w + 1))
}
func (cs ObliqueStereographic) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
k := cs.consts(s)
chi := cs.conformalLat(s, radian(lat), k.c, k.n)
dLam := k.n * (radian(lon) - k.lon0)
B := 1 + math.Sin(chi)*math.Sin(k.chi0) + math.Cos(chi)*math.Cos(k.chi0)*math.Cos(dLam)
east := cs.Eastf + 2*k.R*k.k0*math.Cos(chi)*math.Sin(dLam)/B
north := cs.Northf + 2*k.R*k.k0*(math.Sin(chi)*math.Cos(k.chi0)-math.Cos(chi)*math.Sin(k.chi0)*math.Cos(dLam))/B
return east, north, h
}
func (cs ObliqueStereographic) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
e := s.E()
e2 := s.E2()
k := cs.consts(s)
de := east - cs.Eastf
dn := north - cs.Northf
g := 2 * k.R * k.k0 * math.Tan(math.Pi/4-k.chi0/2)
hh := 4*k.R*k.k0*math.Tan(k.chi0) + g
i := math.Atan2(de, hh+dn)
j := math.Atan2(de, g-dn) - i
chi := k.chi0 + 2*math.Atan((dn-de*math.Tan(j/2))/(2*k.R*k.k0))
lon := (j+2*i)/k.n + k.lon0
psi := math.Log((1+math.Sin(chi))/(k.c*(1-math.Sin(chi)))) / (2 * k.n)
phi := 2*math.Atan(math.Exp(psi)) - math.Pi/2
for range 10 {
sinPhi := math.Sin(phi)
psiI := math.Log(math.Tan(phi/2+math.Pi/4) * math.Pow((1-e*sinPhi)/(1+e*sinPhi), e/2))
dPhi := (psiI - psi) * math.Cos(phi) * (1 - e2*sinPhi*sinPhi) / (1 - e2)
phi -= dPhi
if math.Abs(dPhi) < 1e-14 {
break
}
}
return degree(lon), degree(phi), h
}
package crs
import (
"fmt"
"slices"
"strconv"
"strings"
)
var (
datumStartKeys = []string{"accuracy"}
crsEmbedKeys = []string{"datum", "spheroid", "accuracy"}
)
func flattenDSL(txt string) string {
var b strings.Builder
for line := range strings.SplitSeq(txt, "\n") {
line = stripComment(line)
if line == "" {
continue
}
if b.Len() > 0 {
b.WriteByte(' ')
}
b.WriteString(line)
}
return b.String()
}
func (pp parts) indexOfKey(keys ...string) int {
for i, p := range pp {
if slices.Contains(keys, p.key) {
return i
}
}
return -1
}
func (pp parts) hasKey(key string) bool {
return pp.indexOfKey(key) >= 0
}
func (pp parts) segmentByStartKeys(keys ...string) []parts {
if len(pp) == 0 {
return nil
}
starts := []int{0}
for i := 1; i < len(pp); i++ {
if slices.Contains(keys, pp[i].key) {
starts = append(starts, i)
}
}
out := make([]parts, 0, len(starts))
for i, start := range starts {
end := len(pp)
if i+1 < len(starts) {
end = starts[i+1]
}
out = append(out, pp[start:end])
}
return out
}
func parseSpheroid(txt string) (Spheroid, error) {
return asParts(txt).asSpheroid()
}
func (pp parts) asSpheroid() (Spheroid, error) {
var s Spheroid
if name, ok := pp.getPart("spheroid").asString(); ok {
loaded, err := loadSpheroid(name)
if err != nil {
return s, err
}
return loaded, nil
}
a, hasA := pp.float("a")
fi, hasFi := pp.float("fi")
if hasA {
s.A = a
}
if hasFi {
s.Fi = fi // fi=0 is a sphere (authalic / planetary)
}
if !hasA || s.A == 0 || !hasFi {
return s, fmt.Errorf("spheroid requires spheroid=<name> or a= and fi=")
}
return s, nil
}
func parseDatum(txt string) (Datum, error) {
return asParts(txt).asDatum()
}
func (pp parts) asDatum() (Datum, error) {
var d Datum
if len(pp) == 0 {
return d, fmt.Errorf("empty datum")
}
headerEnd := len(pp)
if i := pp.indexOfKey(datumStartKeys...); i >= 0 {
headerEnd = i
}
header := pp[:headerEnd]
if name, ok := header.getPart("datum").asString(); ok {
d.Name = strings.ToLower(name)
}
if name, ok := header.getPart("spheroid").asString(); ok {
s, err := loadSpheroid(name)
if err != nil {
return d, err
}
d.Spheroid = s
}
if a, ok := header.float("a"); ok {
d.Spheroid.A = a
}
if fi, ok := header.float("fi"); ok {
d.Spheroid.Fi = fi
}
if d.Spheroid.A == 0 {
if d.Name == "" {
return d, fmt.Errorf("datum missing spheroid (need spheroid=<name> or a= and fi=)")
}
loaded, err := loadDatum(d.Name)
if err != nil {
return d, err
}
if headerEnd == len(pp) {
return loaded, nil
}
d.Spheroid = loaded.Spheroid
}
for _, seg := range pp[headerEnd:].segmentByStartKeys(datumStartKeys...) {
if !seg.hasKey("accuracy") {
continue
}
t, err := seg.asTransformation()
if err != nil {
return d, err
}
if t.Operation != nil {
d.Transformations = append(d.Transformations, t)
}
}
sortTransformations(d.Transformations)
return d, nil
}
func sortTransformations(ts []Transformation) {
slices.SortFunc(ts, func(a, b Transformation) int {
if a.Accuracy < b.Accuracy {
return -1
}
if a.Accuracy > b.Accuracy {
return 1
}
if a.BoundingBox.Area() < b.BoundingBox.Area() {
return -1
}
if a.BoundingBox.Area() > b.BoundingBox.Area() {
return 1
}
return 0
})
}
func parseCoordinateReferenceSystem(txt string) (CoordinateReferenceSystem, error) {
return asParts(txt).asCRS()
}
func (pp parts) asCRS() (CoordinateReferenceSystem, error) {
var c CoordinateReferenceSystem
if len(pp) == 0 {
return c, fmt.Errorf("empty crs")
}
if !pp.hasKey("conversion") {
return c, fmt.Errorf("missing conversion")
}
expanded := pp.hasKey("accuracy")
if !expanded {
return pp.asCRSCompact()
}
return pp.asCRSExpanded()
}
func (pp parts) asCRSCompact() (CoordinateReferenceSystem, error) {
var c CoordinateReferenceSystem
fields := make(map[string]string, len(pp))
for _, p := range pp {
if p.key == "" {
continue
}
fields[p.key] = p.value
}
csName, ok := fields["conversion"]
if !ok {
return c, fmt.Errorf("missing conversion")
}
if bboxStr, ok := fields["bbox"]; ok {
bbox, ok := parseBoundingBox(bboxStr)
if !ok {
return c, fmt.Errorf("invalid bbox: %s", bboxStr)
}
c.BoundingBox = bbox
} else {
c.BoundingBox = World
}
cs, err := buildConversion(csName, fields)
if err != nil {
return c, err
}
c.Conversion = cs
switch {
case fields["datum"] != "":
c.Datum, err = loadDatum(fields["datum"])
if err != nil {
return c, fmt.Errorf("datum %s: %w", fields["datum"], err)
}
case fields["spheroid"] != "":
s, err := loadSpheroid(fields["spheroid"])
if err != nil {
return c, fmt.Errorf("spheroid %s: %w", fields["spheroid"], err)
}
c.Datum = Datum{Spheroid: s}
default:
a, okA := parseFloatField(fields, "a")
fi, okFi := parseFloatField(fields, "fi")
if !okA || !okFi {
return c, fmt.Errorf("missing datum or spheroid (need datum=, spheroid=, or a= and fi=)")
}
c.Datum = Datum{Spheroid: Spheroid{A: a, Fi: fi}}
}
c.Datum, err = c.Datum.Intersects(c.BoundingBox)
if err != nil {
return c, fmt.Errorf("datum: %w", err)
}
return c, nil
}
func parseFloatField(fields map[string]string, key string) (float64, bool) {
v, ok := fields[key]
if !ok || v == "" {
return 0, false
}
f, err := strconv.ParseFloat(v, 64)
if err != nil {
return 0, false
}
return f, true
}
func (pp parts) asCRSExpanded() (CoordinateReferenceSystem, error) {
var c CoordinateReferenceSystem
split := pp.indexOfKey(crsEmbedKeys...)
if split < 0 {
return c, fmt.Errorf("expanded crs missing datum/spheroid/accuracy")
}
header := pp[:split]
rest := pp[split:]
fields := make(map[string]string, len(header))
for _, p := range header {
if p.key == "" {
continue
}
fields[p.key] = p.value
}
csName, ok := fields["conversion"]
if !ok {
return c, fmt.Errorf("missing conversion")
}
if bboxStr, ok := fields["bbox"]; ok {
bbox, ok := parseBoundingBox(bboxStr)
if !ok {
return c, fmt.Errorf("invalid bbox: %s", bboxStr)
}
c.BoundingBox = bbox
} else {
c.BoundingBox = World
}
cs, err := buildConversion(csName, fields)
if err != nil {
return c, err
}
c.Conversion = cs
c.Datum, err = rest.asDatum()
if err != nil {
return c, err
}
c.Datum, err = c.Datum.Intersects(c.BoundingBox)
if err != nil {
return c, fmt.Errorf("datum: %w", err)
}
return c, nil
}
package crs
import (
"math"
)
type PolarStereographicA struct {
Lonf, Latf, Scale, Eastf, Northf float64
}
func (cs PolarStereographicA) String() string {
return build("polar_stereographic_a").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs PolarStereographicA) south() bool {
return cs.Latf < 0
}
func (cs PolarStereographicA) k0() float64 {
if cs.Scale == 0 {
return 1
}
return cs.Scale
}
func (cs PolarStereographicA) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
return polarStereoFromGeographic(s, lon, lat, h, cs.Lonf, cs.Eastf, cs.Northf, cs.k0(), cs.south())
}
func (cs PolarStereographicA) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
return polarStereoToGeographic(s, east, north, h, cs.Lonf, cs.Eastf, cs.Northf, cs.k0(), cs.south())
}
type PolarStereographicB struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs PolarStereographicB) String() string {
return build("polar_stereographic_b").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs PolarStereographicB) south() bool {
return cs.Latf < 0
}
func (cs PolarStereographicB) k0(s Spheroid) float64 {
e := polarStereoE(s)
mF, tF := polarStereoMF(s, radian(cs.Latf), cs.south())
return mF * polarStereoC(e) / (2 * tF)
}
func (cs PolarStereographicB) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
return polarStereoFromGeographic(s, lon, lat, h, cs.Lonf, cs.Eastf, cs.Northf, cs.k0(s), cs.south())
}
func (cs PolarStereographicB) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
return polarStereoToGeographic(s, east, north, h, cs.Lonf, cs.Eastf, cs.Northf, cs.k0(s), cs.south())
}
// PolarStereographicC is EPSG method 9830.
// Latf is the standard parallel (latitude of false origin); FE/FN apply there, not at the pole.
type PolarStereographicC struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs PolarStereographicC) String() string {
return build("polar_stereographic_c").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs PolarStereographicC) south() bool {
return cs.Latf < 0
}
func (cs PolarStereographicC) k0(s Spheroid) float64 {
e := polarStereoE(s)
mF, tF := polarStereoMF(s, radian(cs.Latf), cs.south())
return mF * polarStereoC(e) / (2 * tF)
}
func (cs PolarStereographicC) rhoF(s Spheroid) float64 {
mF, _ := polarStereoMF(s, radian(cs.Latf), cs.south())
return s.A * mF
}
// northfAtPole shifts NF so shared A/B helpers (origin at pole) match variant C.
func (cs PolarStereographicC) northfAtPole(s Spheroid) float64 {
rhoF := cs.rhoF(s)
if cs.south() {
return cs.Northf - rhoF
}
return cs.Northf + rhoF
}
func (cs PolarStereographicC) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
return polarStereoFromGeographic(s, lon, lat, h, cs.Lonf, cs.Eastf, cs.northfAtPole(s), cs.k0(s), cs.south())
}
func (cs PolarStereographicC) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
return polarStereoToGeographic(s, east, north, h, cs.Lonf, cs.Eastf, cs.northfAtPole(s), cs.k0(s), cs.south())
}
func polarStereoE(s Spheroid) float64 {
return math.Sqrt(s.E2())
}
func polarStereoC(e float64) float64 {
return math.Exp(0.5 * ((1+e)*math.Log(1+e) + (1-e)*math.Log(1-e)))
}
func polarStereoT(phi, e float64, south bool) float64 {
sinPhi := math.Sin(phi)
es := e * sinPhi
f := math.Pow((1+es)/(1-es), e/2)
if south {
return math.Tan(math.Pi/4+phi/2) / f
}
return math.Tan(math.Pi/4-phi/2) * f
}
func polarStereoMF(s Spheroid, phiF float64, south bool) (mF, tF float64) {
e := polarStereoE(s)
sinF := math.Sin(phiF)
mF = math.Cos(phiF) / math.Sqrt(1-s.E2()*sinF*sinF)
tF = polarStereoT(phiF, e, south)
return mF, tF
}
func polarStereoFromGeographic(s Spheroid, lon, lat, h, lonf, eastf, northf, k0 float64, south bool) (float64, float64, float64) {
e := polarStereoE(s)
c := polarStereoC(e)
phi := radian(lat)
dLam := radian(lon - lonf)
t := polarStereoT(phi, e, south)
rho := 2 * s.A * k0 * t / c
sinD := math.Sin(dLam)
cosD := math.Cos(dLam)
east := eastf + rho*sinD
north := northf + rho*cosD
if !south {
north = northf - rho*cosD
}
return east, north, h
}
func polarStereoToGeographic(s Spheroid, east, north, h, lonf, eastf, northf, k0 float64, south bool) (float64, float64, float64) {
e2 := s.E2()
e4 := e2 * e2
e6 := e4 * e2
e8 := e4 * e4
c := polarStereoC(polarStereoE(s))
de := east - eastf
dn := north - northf
rho := math.Hypot(de, dn)
t := rho * c / (2 * s.A * k0)
var chi float64
if south {
chi = 2*math.Atan(t) - math.Pi/2
} else {
chi = math.Pi/2 - 2*math.Atan(t)
}
phi := chi +
(e2/2+5*e4/24+e6/12+13*e8/360)*math.Sin(2*chi) +
(7*e4/48+29*e6/240+811*e8/11520)*math.Sin(4*chi) +
(7*e6/120+81*e8/1120)*math.Sin(6*chi) +
(4279*e8/161280)*math.Sin(8*chi)
var lon float64
if de == 0 {
lon = lonf
} else if south {
lon = lonf + degree(math.Atan2(de, dn))
} else {
lon = lonf + degree(math.Atan2(de, -dn))
}
return lon, degree(phi), h
}
package crs
import (
"math"
)
type AmericanPolyconic struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs AmericanPolyconic) String() string {
return build("american_polyconic").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs AmericanPolyconic) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
phi := radian(lat)
phi0 := radian(cs.Latf)
dLam := radian(lon - cs.Lonf)
m0 := meridianDistance(s, phi0)
if math.Abs(phi) <= 1e-10 {
return cs.Eastf + s.A*dLam, cs.Northf - m0, h
}
sinPhi, cosPhi := math.Sincos(phi)
ms := cosPhi / (math.Sqrt(1-s.E2()*sinPhi*sinPhi) * sinPhi)
e := dLam * sinPhi
east := cs.Eastf + s.A*ms*math.Sin(e)
north := cs.Northf + (meridianDistance(s, phi) - m0) + s.A*ms*(1-math.Cos(e))
return east, north, h
}
func (cs AmericanPolyconic) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
a := s.A
e2 := s.E2()
oneEs := 1 - e2
m0 := meridianDistance(s, radian(cs.Latf))
x := (east - cs.Eastf) / a
y := (north - cs.Northf) / a
y += m0 / a
if math.Abs(y) <= 1e-10 {
return degree(radian(cs.Lonf) + x), 0, h
}
r := y*y + x*x
phi := y
for range 20 {
sinPhi, cosPhi := math.Sincos(phi)
if math.Abs(cosPhi) < 1e-12 {
return cs.Lonf, degree(math.Copysign(math.Pi/2, phi)), h
}
s2ph := sinPhi * cosPhi
mlp := math.Sqrt(1 - e2*sinPhi*sinPhi)
c := sinPhi * mlp / cosPhi
ml := meridianDistance(s, phi) / a
mlb := ml*ml + r
mlp = oneEs / (mlp * mlp * mlp)
dPhi := (ml + ml + c*mlb - 2*y*(c*ml+1)) /
(e2*s2ph*(mlb-2*y*ml)/c + 2*(y-ml)*(c*mlp-1/s2ph) - mlp - mlp)
phi += dPhi
if math.Abs(dPhi) <= 1e-12 {
break
}
}
sinPhi := math.Sin(phi)
if math.Abs(sinPhi) < 1e-12 {
return degree(radian(cs.Lonf) + x), 0, h
}
lam := math.Asin(clamp(x*math.Tan(phi)*math.Sqrt(1-e2*sinPhi*sinPhi), -1, 1)) / sinPhi
return degree(lam + radian(cs.Lonf)), degree(phi), h
}
package crs
import (
"math"
)
type PositionVector struct {
Tx float64
Ty float64
Tz float64
Rx float64
Ry float64
Rz float64
Ds float64
}
func (pc PositionVector) String() string {
return build("position_vector").addAll(
"tx", pc.Tx,
"ty", pc.Ty,
"tz", pc.Tz,
"rx", pc.Rx,
"ry", pc.Ry,
"rz", pc.Rz,
"ds", pc.Ds,
).String()
}
func (pv PositionVector) FromTarget(source Spheroid, target Spheroid, lon0 float64, lat0 float64, h0 float64) (float64, float64, float64, error) {
x0, y0, z0 := source.GeographicToGeocentric(lon0, lat0, h0)
x, y, z := calcHelmert(x0, y0, z0, -pv.Tx, -pv.Ty, -pv.Tz, -pv.Rx, -pv.Ry, -pv.Rz, -pv.Ds)
lon, lat, h := target.GeocentricToGeographic(x, y, z)
return lon, lat, h, nil
}
func (pv PositionVector) ToTarget(source Spheroid, target Spheroid, lon float64, lat float64, h float64) (float64, float64, float64, error) {
x, y, z := source.GeographicToGeocentric(lon, lat, h)
x0, y0, z0 := calcHelmert(x, y, z, pv.Tx, pv.Ty, pv.Tz, pv.Rx, pv.Ry, pv.Rz, pv.Ds)
lon0, lat0, h0 := target.GeocentricToGeographic(x0, y0, z0)
return lon0, lat0, h0, nil
}
func calcHelmert(x, y, z, tx, ty, tz, rx, ry, rz, ds float64) (x1, y1, z1 float64) {
const (
asec = math.Pi / 648000
ppm = 0.000001
)
x1 = (1+ds*ppm)*(x+z*ry*asec-y*rz*asec) + tx
y1 = (1+ds*ppm)*(y+x*rz*asec-z*rx*asec) + ty
z1 = (1+ds*ppm)*(z+y*rx*asec-x*ry*asec) + tz
return x1, y1, z1
}
package crs
import (
"math"
)
// SwissObliqueMercator is EPSG method 9814 (PROJ +proj=somerc).
// Alpha/Gamma are retained for EPSG parameter round-trip; the classic Swiss
// formulas assume azimuth = rectified grid angle = 90°.
type SwissObliqueMercator struct {
Lonf float64
Latf float64
Scale float64
Eastf float64
Northf float64
Alpha float64
Gamma float64
}
func (cs SwissObliqueMercator) String() string {
return build("swiss_oblique_mercator").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
"alpha", cs.Alpha,
"gamma", cs.Gamma,
).String()
}
type somercConsts struct {
c, K, kR, cosp0, sinp0, hlfE, lon0 float64
}
func (som SwissObliqueMercator) consts(s Spheroid) somercConsts {
e := s.E()
e2 := s.E2()
k0 := som.Scale
if k0 == 0 {
k0 = 1
}
phi0 := radian(som.Latf)
cos0 := math.Cos(phi0)
cos02 := cos0 * cos0
c := math.Sqrt(1 + e2*cos02*cos02/(1-e2))
sp := math.Sin(phi0)
sinp0 := sp / c
phip0 := math.Asin(clamp(sinp0, -1, 1))
cosp0 := math.Cos(phip0)
esp := e * sp
K := math.Log(math.Tan(math.Pi/4+phip0/2)) -
c*(math.Log(math.Tan(math.Pi/4+phi0/2))-0.5*e*math.Log((1+esp)/(1-esp)))
kR := s.A * k0 * math.Sqrt(1-e2) / (1 - esp*esp)
return somercConsts{
c: c, K: K, kR: kR, cosp0: cosp0, sinp0: sinp0, hlfE: 0.5 * e, lon0: radian(som.Lonf),
}
}
func (som SwissObliqueMercator) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
q := som.consts(s)
phi := radian(lat)
lam := radian(lon) - q.lon0
sp := s.E() * math.Sin(phi)
phip := 2*math.Atan(math.Exp(
q.c*(math.Log(math.Tan(math.Pi/4+phi/2))-q.hlfE*math.Log((1+sp)/(1-sp)))+q.K,
)) - math.Pi/2
lamp := q.c * lam
cp := math.Cos(phip)
phipp := math.Asin(clamp(q.cosp0*math.Sin(phip)-q.sinp0*cp*math.Cos(lamp), -1, 1))
lampp := math.Asin(clamp(cp*math.Sin(lamp)/math.Cos(phipp), -1, 1))
east := q.kR*lampp + som.Eastf
north := q.kR*math.Log(math.Tan(math.Pi/4+phipp/2)) + som.Northf
return east, north, h
}
func (som SwissObliqueMercator) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
q := som.consts(s)
e := s.E()
roneEs := 1 / (1 - s.E2())
phipp := 2 * (math.Atan(math.Exp((north-som.Northf)/q.kR)) - math.Pi/4)
lampp := (east - som.Eastf) / q.kR
cp := math.Cos(phipp)
phip := math.Asin(clamp(q.cosp0*math.Sin(phipp)+q.sinp0*cp*math.Cos(lampp), -1, 1))
lamp := math.Asin(clamp(cp*math.Sin(lampp)/math.Cos(phip), -1, 1))
con := (q.K - math.Log(math.Tan(math.Pi/4+phip/2))) / q.c
phi := phip
for range 6 {
esp := e * math.Sin(phi)
delp := (con + math.Log(math.Tan(math.Pi/4+phi/2)) -
q.hlfE*math.Log((1+esp)/(1-esp))) *
(1 - esp*esp) * math.Cos(phi) * roneEs
phi -= delp
if math.Abs(delp) < 1e-10 {
break
}
}
return degree(q.lon0 + lamp/q.c), degree(phi), h
}
package crs
import (
"fmt"
"math"
)
type Spheroid struct {
Name string
A float64
Fi float64
}
func (s Spheroid) String() string {
return fmt.Sprintf("a=%s fi=%s", formatFloat(s.A), formatFloat(s.Fi))
}
func (s Spheroid) F() float64 {
if s.Fi == 0 {
return 0
}
return 1 / s.Fi
}
func (s Spheroid) A2() float64 {
return s.A * s.A
}
func (s Spheroid) F2() float64 {
f := s.F()
return f * f
}
func (s Spheroid) B() float64 {
return s.A * (1 - s.F())
}
func (s Spheroid) E2() float64 {
if s.Fi == 0 {
return 0
}
return 2/s.Fi - s.F2()
}
func (s Spheroid) E() float64 {
return math.Sqrt(s.E2())
}
func (s Spheroid) E4() float64 {
e2 := s.E2()
return e2 * e2
}
func (s Spheroid) E6() float64 {
e2 := s.E2()
return e2 * e2 * e2
}
func (s Spheroid) Ei() float64 {
e2 := s.E2()
t := math.Sqrt(1 - e2)
return (1 - t) / (1 + t)
}
func (s Spheroid) Ei2() float64 {
ei := s.Ei()
return ei * ei
}
func (s Spheroid) Ei3() float64 {
ei := s.Ei()
return ei * ei * ei
}
func (s Spheroid) Ei4() float64 {
ei := s.Ei()
return ei * ei * ei * ei
}
func (s Spheroid) GeographicToGeocentric(lon, lat, h float64) (x, y, z float64) {
n := s.A / math.Sqrt(1-s.E2()*intPow(math.Sin(radian(lat)), 2))
x = (n + h) * math.Cos(radian(lon)) * math.Cos(radian(lat))
y = (n + h) * math.Cos(radian(lat)) * math.Sin(radian(lon))
z = (n*intPow(s.A*(1-s.F()), 2)/(s.A2()) + h) * math.Sin(radian(lat))
return x, y, z
}
func (s Spheroid) GeocentricToGeographic(x, y, z float64) (lon, lat, h float64) {
sd := math.Sqrt(x*x + y*y)
T := math.Atan(z * s.A / (sd * s.B()))
B := math.Atan((z + s.E2()*(s.A2())/s.B()*
intPow(math.Sin(T), 3)) / (sd - s.E2()*s.A*intPow(math.Cos(T), 3)))
n := s.A / math.Sqrt(1-s.E2()*intPow(math.Sin(B), 2))
h = sd/math.Cos(B) - n
lon = degree(math.Atan2(y, x))
lat = degree(B)
return lon, lat, h
}
package crs
import (
"embed"
"fmt"
"io"
"strings"
"sync"
)
//go:embed spheroid/*.txt
var spheroidDir embed.FS
var spheroidStore sync.Map
func RegisterSpheroid(s Spheroid) {
spheroidStore.Store(s.Name, s)
}
func loadSpheroid(name string) (s Spheroid, err error) {
name = strings.ToLower(strings.TrimSpace(name))
v, ok := spheroidStore.Load(name)
if ok {
return v.(Spheroid), nil
}
defer func() {
if err == nil {
spheroidStore.Store(name, s)
}
}()
file, err := spheroidDir.Open(fmt.Sprintf("spheroid/%s.txt", name))
if err != nil {
return s, err
}
data, err := io.ReadAll(file)
if err != nil {
return s, err
}
s, err = parseSpheroid(string(data))
if err != nil {
return s, err
}
s.Name = name
return s, nil
}
package crs
import "math"
// timeSpecificEpochTol is the max |coordEpoch − TransformationEpoch| (years)
// for a time-specific Helmert to participate in path selection.
const timeSpecificEpochTol = 1e-6
type TimeSpecificPositionVector struct {
Tx, Ty, Tz, Rx, Ry, Rz, Ds float64
TransformationEpoch float64
}
// timeSpecificEpoch returns the transformation epoch when op is (or wraps)
// TimeSpecificPositionVector.
func timeSpecificEpoch(op Operation) (epoch float64, ok bool) {
switch o := op.(type) {
case TimeSpecificPositionVector:
return o.TransformationEpoch, true
case Inverse:
return timeSpecificEpoch(o.Operation)
default:
return 0, false
}
}
// timeSpecificAllowed reports whether a time-specific op may be used for the
// given coordinate epoch. Non-time-specific ops always return true.
// Without a coordinate epoch (Transform), time-specific edges are excluded.
func timeSpecificAllowed(op Operation, hasCoordEpoch bool, coordEpoch float64) bool {
te, isTS := timeSpecificEpoch(op)
if !isTS {
return true
}
if !hasCoordEpoch {
return false
}
return math.Abs(coordEpoch-te) <= timeSpecificEpochTol
}
func (t TimeSpecificPositionVector) String() string {
return build("time_specific_position_vector").addAll(
"tx", t.Tx,
"ty", t.Ty,
"tz", t.Tz,
"rx", t.Rx,
"ry", t.Ry,
"rz", t.Rz,
"ds", t.Ds,
"epoch", t.TransformationEpoch,
).String()
}
func (t TimeSpecificPositionVector) asPositionVector() PositionVector {
return PositionVector{
Tx: t.Tx, Ty: t.Ty, Tz: t.Tz,
Rx: t.Rx, Ry: t.Ry, Rz: t.Rz,
Ds: t.Ds,
}
}
func (t TimeSpecificPositionVector) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return t.asPositionVector().ToTarget(source, target, lon, lat, h)
}
func (t TimeSpecificPositionVector) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return t.asPositionVector().FromTarget(source, target, lon0, lat0, h0)
}
type TimeDependentPositionVector struct {
Tx, Ty, Tz, Rx, Ry, Rz, Ds float64
TxRate, TyRate, TzRate, RxRate, RyRate, RzRate, DsRate float64
ReferenceEpoch float64
}
func (t TimeDependentPositionVector) String() string {
return build("time_dependent_position_vector").addAll(
"tx", t.Tx,
"ty", t.Ty,
"tz", t.Tz,
"rx", t.Rx,
"ry", t.Ry,
"rz", t.Rz,
"ds", t.Ds,
"dtx", t.TxRate,
"dty", t.TyRate,
"dtz", t.TzRate,
"drx", t.RxRate,
"dry", t.RyRate,
"drz", t.RzRate,
"dds", t.DsRate,
"epoch", t.ReferenceEpoch,
).String()
}
func (t TimeDependentPositionVector) at(epoch float64) PositionVector {
dt := epoch - t.ReferenceEpoch
return PositionVector{
Tx: t.Tx + t.TxRate*dt,
Ty: t.Ty + t.TyRate*dt,
Tz: t.Tz + t.TzRate*dt,
Rx: t.Rx + t.RxRate*dt,
Ry: t.Ry + t.RyRate*dt,
Rz: t.Rz + t.RzRate*dt,
Ds: t.Ds + t.DsRate*dt,
}
}
// ToTarget evaluates Helmert parameters at ReferenceEpoch (dt=0). Use
// Datum.AtEpoch / TransformAt so rates apply at a coordinate epoch.
func (t TimeDependentPositionVector) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return t.at(t.ReferenceEpoch).ToTarget(source, target, lon, lat, h)
}
// FromTarget evaluates Helmert parameters at ReferenceEpoch (dt=0). Use
// Datum.AtEpoch / TransformAt so rates apply at a coordinate epoch.
func (t TimeDependentPositionVector) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return t.at(t.ReferenceEpoch).FromTarget(source, target, lon0, lat0, h0)
}
package crs
import (
"math"
)
type TransverseMercator struct {
Lonf float64
Latf float64
Scale float64
Eastf float64
Northf float64
ZoneWidth float64
}
func (cs TransverseMercator) String() string {
return build("transverse_mercator").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
"zone_width", cs.ZoneWidth,
).String()
}
func (cs TransverseMercator) zone(lon float64) int {
w := cs.ZoneWidth
d := lon - cs.Lonf
for d < 0 {
d += 360
}
for d >= 360 {
d -= 360
}
z := 1 + int(math.Floor(d/w))
nZones := int(math.Round(360 / w))
if z > nZones {
z = nZones
}
if z < 1 {
z = 1
}
return z
}
func (cs TransverseMercator) centralMeridian(zone int) float64 {
return cs.Lonf + (float64(zone)-0.5)*cs.ZoneWidth
}
// fixedLonf returns a single-zone TM with the given central meridian.
func (cs TransverseMercator) fixedLonf(lonf float64) TransverseMercator {
cs.Lonf = lonf
cs.ZoneWidth = 0
return cs
}
func (cs TransverseMercator) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
if cs.ZoneWidth > 0 {
return cs.fixedLonf(cs.centralMeridian(cs.zone(lon))).FromGeographic(s, lon, lat, h)
}
phiO := radian(cs.Latf)
lambdaO := radian(cs.Lonf)
n := s.F() / (2 - s.F())
n2 := n * n
n3 := n * n * n
n4 := n * n * n * n
B := (s.A / (1 + n)) * (1 + n2/4 + n4/64)
h1 := n/2.0 - (2/3.0)*n2 + (5/16.0)*n3 + (41/180.0)*n4
h2 := (13/48.0)*n2 - (3/5.0)*n3 + (557/1440.0)*n4
h3 := (61/240.0)*n3 - (103/140.0)*n4
h4 := (49561 / 161280.0) * n4
var MO float64
switch phiO {
case 0:
MO = 0
case math.Pi / 2:
MO = B * (math.Pi / 2)
case -math.Pi / 2:
MO = B * (-math.Pi / 2)
default:
Q0 := math.Asinh(math.Tan(phiO)) - (s.E() * math.Atanh(s.E()*math.Sin(phiO)))
xi00 := math.Atan(math.Sinh(Q0))
xi01 := h1 * math.Sin(2*xi00)
xi02 := h2 * math.Sin(4*xi00)
xi03 := h3 * math.Sin(6*xi00)
xi04 := h4 * math.Sin(8*xi00)
xi0 := xi00 + xi01 + xi02 + xi03 + xi04
MO = B * xi0
}
phi := radian(lat)
lambda := radian(lon)
Q := math.Asinh(math.Tan(phi)) - s.E()*math.Atanh(s.E()*math.Sin(phi))
beta := math.Atan(math.Sinh(Q))
eta0 := math.Atanh(math.Cos(beta) * math.Sin(lambda-lambdaO))
xi0 := math.Asin(math.Sin(beta) * math.Cosh(eta0))
xi1 := h1 * math.Sin(2*xi0) * math.Cosh(2*eta0)
xi2 := h2 * math.Sin(4*xi0) * math.Cosh(4*eta0)
xi3 := h3 * math.Sin(6*xi0) * math.Cosh(6*eta0)
xi4 := h4 * math.Sin(8*xi0) * math.Cosh(8*eta0)
xi := xi0 + xi1 + xi2 + xi3 + xi4
eta1 := h1 * math.Cos(2*xi0) * math.Sinh(2*eta0)
eta2 := h2 * math.Cos(4*xi0) * math.Sinh(4*eta0)
eta3 := h3 * math.Cos(6*xi0) * math.Sinh(6*eta0)
eta4 := h4 * math.Cos(8*xi0) * math.Sinh(8*eta0)
eta := eta0 + eta1 + eta2 + eta3 + eta4
east := cs.Eastf + cs.Scale*B*eta
north := cs.Northf + cs.Scale*(B*xi-MO)
return east, north, h
}
func (cs TransverseMercator) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
if cs.ZoneWidth > 0 {
nZones := int(math.Round(360 / cs.ZoneWidth))
bestDist := math.Inf(1)
bestLon, bestLat, bestH := 0.0, 0.0, h
found := false
for z := 1; z <= nZones; z++ {
lon, lat, hOut := cs.fixedLonf(cs.centralMeridian(z)).ToGeographic(s, east, north, h)
minLon := cs.Lonf + float64(z-1)*cs.ZoneWidth
maxLon := minLon + cs.ZoneWidth
lonN := lon
for lonN < cs.Lonf {
lonN += 360
}
for lonN >= cs.Lonf+360 {
lonN -= 360
}
if lonN < minLon || lonN >= maxLon {
continue
}
dist := math.Abs(lon - cs.centralMeridian(z))
if dist < bestDist {
bestDist = dist
bestLon, bestLat, bestH = lon, lat, hOut
found = true
}
}
if found {
return bestLon, bestLat, bestH
}
return cs.fixedLonf(cs.centralMeridian(1)).ToGeographic(s, east, north, h)
}
phiO := radian(cs.Latf)
lambdaO := radian(cs.Lonf)
n := s.F() / (2 - s.F())
n2 := n * n
n3 := n * n * n
n4 := n * n * n * n
B := (s.A / (1 + n)) * (1 + n2/4 + n4/64)
h1 := n/2.0 - (2/3.0)*n2 + (5/16.0)*n3 + (41/180.0)*n4
h2 := (13/48.0)*n2 - (3/5.0)*n3 + (557/1440.0)*n4
h3 := (61/240.0)*n3 - (103/140.0)*n4
h4 := (49561 / 161280.0) * n4
e := math.Sqrt(s.E2())
var MO float64
switch phiO {
case 0:
MO = 0
case math.Pi / 2:
MO = B * (math.Pi / 2)
case -math.Pi / 2:
MO = B * (-math.Pi / 2)
default:
Q0 := math.Asinh(math.Tan(phiO)) - (e * math.Atanh(e*math.Sin(phiO)))
xi00 := math.Atan(math.Sinh(Q0))
xi01 := h1 * math.Sin(2*xi00)
xi02 := h2 * math.Sin(4*xi00)
xi03 := h3 * math.Sin(6*xi00)
xi04 := h4 * math.Sin(8*xi00)
xi0 := xi00 + xi01 + xi02 + xi03 + xi04
MO = B * xi0
}
h1i := n/2.0 - (2/3.0)*n2 + (37/96.0)*n3 + (1/360.0)*n4
h2i := (1/48.0)*n2 - (1/15.0)*n3 + (437/1440.0)*n4
h3i := (17/480.0)*n3 - (37/840.0)*n4
h4i := (4397 / 161280.0) * n4
etai := (east - cs.Eastf) / (B * cs.Scale)
xii := ((north - cs.Northf) + cs.Scale*MO) / (B * cs.Scale)
xi1i := h1i * math.Sin(2*xii) * math.Cosh(2*etai)
xi2i := h2i * math.Sin(4*xii) * math.Cosh(4*etai)
xi3i := h3i * math.Sin(6*xii) * math.Cosh(6*etai)
xi4i := h4i * math.Sin(8*xii) * math.Cosh(8*etai)
xi0i := xii - (xi1i + xi2i + xi3i + xi4i)
eta1i := h1i * math.Cos(2*xii) * math.Sinh(2*etai)
eta2i := h2i * math.Cos(4*xii) * math.Sinh(4*etai)
eta3i := h3i * math.Cos(6*xii) * math.Sinh(6*etai)
eta4i := h4i * math.Cos(8*xii) * math.Sinh(8*etai)
eta0i := etai - (eta1i + eta2i + eta3i + eta4i)
betai := math.Asin(math.Sin(xi0i) / math.Cosh(eta0i))
Qi := math.Asinh(math.Tan(betai))
Qii := Qi + (s.E() * math.Atanh(s.E()*math.Tanh(Qi)))
for range 15 {
newQ := Qi + (s.E() * math.Atanh(s.E()*math.Tanh(Qii)))
if math.Abs(newQ-Qii) < 1e-14 {
break
}
Qii = newQ
}
phi := math.Atan(math.Sinh(Qii))
lambda := lambdaO + math.Asin(math.Tanh(eta0i)/math.Cos(betai))
return degree(lambda), degree(phi), h
}
package crs
import "fmt"
type Operation interface {
ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error)
FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error)
}
type Transformation struct {
Target *Datum
Accuracy float64
BoundingBox BoundingBox
Operation Operation
}
func (t Transformation) requireTarget() error {
if t.Target == nil || t.Target.Name == "" {
return UnsupportedError{
Err: fmt.Errorf("no target"),
}
}
return nil
}
func (t Transformation) String() string {
return build("").addAll(
"operation", t.Operation,
"target", t.Target.Name,
"accuracy", t.Accuracy,
"bbox", t.BoundingBox,
).String()
}
func (t Transformation) ToWGS84(source Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return t.toWGS84Visited(source, lon, lat, h, nil, 0)
}
// toWGS84Visited applies this hop toward WGS 84. greenwichLon is lon expressed east
// of Greenwich for area-of-use checks; lon itself may be relative to a local PM.
func (t Transformation) toWGS84Visited(source Spheroid, lon, lat, h float64, visited map[string]bool, greenwichLon float64) (float64, float64, float64, error) {
if err := t.requireTarget(); err != nil {
return 0, 0, 0, err
}
if !t.BoundingBox.Contains(greenwichLon, lat) {
return 0, 0, 0, OutOfBoundsError{
Err: fmt.Errorf("[%f,%f]: %s", greenwichLon, lat, t.BoundingBox),
}
}
lon, lat, h, err := t.Operation.ToTarget(source, t.Target.Spheroid, lon, lat, h)
if err != nil {
return 0, 0, 0, err
}
return t.Target.toWGS84Visited(lon, lat, h, visited)
}
func (t Transformation) ToDatum(source Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
if err := t.requireTarget(); err != nil {
return 0, 0, 0, err
}
return t.Operation.ToTarget(source, t.Target.Spheroid, lon, lat, h)
}
func (t Transformation) FromDatum(owner Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
if err := t.requireTarget(); err != nil {
return 0, 0, 0, err
}
return t.Operation.FromTarget(t.Target.Spheroid, owner, lon, lat, h)
}
func (t Transformation) FromWGS84(source Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return t.fromWGS84Visited(source, lon0, lat0, h0, nil)
}
func (t Transformation) fromWGS84Visited(source Spheroid, lon0, lat0, h0 float64, visited map[string]bool) (float64, float64, float64, error) {
if err := t.requireTarget(); err != nil {
return 0, 0, 0, err
}
if !t.BoundingBox.Contains(lon0, lat0) {
return 0, 0, 0, OutOfBoundsError{
Err: fmt.Errorf("[%f,%f]: %s", lon0, lat0, t.BoundingBox),
}
}
lon0, lat0, h0, err := t.Target.fromWGS84Visited(lon0, lat0, h0, visited)
if err != nil {
return 0, 0, 0, err
}
return t.Operation.FromTarget(t.Target.Spheroid, source, lon0, lat0, h0)
}
package crs
type TunisiaMiningGrid struct {
Lonf, Latf, Eastf, Northf float64
}
func (cs TunisiaMiningGrid) String() string {
return build("tunisia_mining_grid").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
const tunisiaLonPerKm = 0.012185 // grads per kilometre easting
func tunisiaLatPerKm(aboveOrigin bool) float64 {
if aboveOrigin {
return 0.010015
}
return 0.01002
}
func (cs TunisiaMiningGrid) FromGeographic(_ Spheroid, lon, lat, h float64) (float64, float64, float64) {
latG := lat / 0.9
lonG := lon / 0.9
lat0G := cs.Latf / 0.9
lon0G := cs.Lonf / 0.9
a := tunisiaLatPerKm(latG > lat0G)
east := cs.Eastf + (lonG-lon0G)/tunisiaLonPerKm*1000
north := cs.Northf + (latG-lat0G)/a*1000
return east, north, h
}
func (cs TunisiaMiningGrid) ToGeographic(_ Spheroid, east, north, h float64) (float64, float64, float64) {
lat0G := cs.Latf / 0.9
lon0G := cs.Lonf / 0.9
eKm := (east - cs.Eastf) / 1000
nKm := (north - cs.Northf) / 1000
a := tunisiaLatPerKm(nKm > 0)
lonG := lon0G + eKm*tunisiaLonPerKm
latG := lat0G + nKm*a
return lonG * 0.9, latG * 0.9, h
}
package crs
import (
"errors"
"fmt"
"io"
"math"
"path"
"strings"
"sync"
)
// VelocityGrid is PROJ +proj=deformation with a fixed timespan dt (years).
// Grid bands are east/north/up velocities in mm/year (GTG TYPE=VELOCITY).
// When Epoch != 0, Datum.AtEpoch sets Dt = epoch - Epoch.
type VelocityGrid struct {
Grid string
Dt float64
Epoch float64
}
func (v VelocityGrid) String() string {
return build("velocity_grid").addAll(
"grid", v.Grid,
"dt", v.Dt,
"epoch", v.Epoch,
).String()
}
// ToTarget implements [Operation]: apply +dt · V in geocentric space.
func (v VelocityGrid) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return v.apply(source, target, lon, lat, h, v.Dt)
}
// FromTarget implements [Operation]: apply −dt · V.
func (v VelocityGrid) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return v.apply(source, target, lon0, lat0, h0, -v.Dt)
}
func (v VelocityGrid) apply(source, target Spheroid, lon, lat, h, dt float64) (float64, float64, float64, error) {
if dt == 0 {
return lon, lat, h, nil
}
grid, err := loadVelocityGrid(v.Grid)
if err != nil {
return 0, 0, 0, err
}
e, n, u, err := grid.Velocity(lon, lat) // mm/year
if err != nil {
return 0, 0, 0, err
}
// mm/year → m/year
e *= 0.001
n *= 0.001
u *= 0.001
phi := radian(lat)
lam := radian(lon)
sinPhi, cosPhi := math.Sincos(phi)
sinLam, cosLam := math.Sincos(lam)
vx := (-sinPhi*cosLam*n - sinLam*e + cosPhi*cosLam*u) * dt
vy := (-sinPhi*sinLam*n + cosLam*e + cosPhi*sinLam*u) * dt
vz := (cosPhi*n + sinPhi*u) * dt
x, y, z := source.GeographicToGeocentric(lon, lat, h)
lon1, lat1, h1 := target.GeocentricToGeographic(x+vx, y+vy, z+vz)
return lon1, lat1, h1, nil
}
type velocityGridData struct {
File string
LonUL float64
LatUL float64
DLon float64
DLat float64
Columns int
Rows int
East []float32
North []float32
Up []float32
}
var velocityGridStore sync.Map
func loadVelocityGrid(name string) (grid *velocityGridData, err error) {
name = gridFilename(name)
if v, ok := velocityGridStore.Load(name); ok {
return v.(*velocityGridData), nil
}
defer func() {
if err == nil {
grid.File = name
velocityGridStore.Store(name, grid)
}
}()
r, err := openGridReader(name)
if err != nil {
return nil, UnsupportedError{
Err: fmt.Errorf("velocity grid: %s", name),
}
}
defer r.Close() //nolint:errcheck
switch strings.ToLower(path.Ext(name)) {
case ".tif":
return parseVelocityGridTIFF(r)
default:
return nil, fmt.Errorf("crs: unsupported velocity grid extension %q", path.Ext(name))
}
}
// Velocity returns east/north/up in mm/year at (lon, lat) degrees.
func (g *velocityGridData) Velocity(lon, lat float64) (east, north, up float64, err error) {
if g.Columns < 2 || g.Rows < 2 || len(g.East) == 0 {
return 0, 0, 0, UnsupportedError{
Err: errors.New("invalid grid"),
}
}
colF := (lon - g.LonUL) / g.DLon
rowF := (g.LatUL - lat) / g.DLat
col, dx, ok := interpolateIndex(colF, g.Columns)
if !ok {
return 0, 0, 0, UnsupportedError{
Err: errors.New("cannot interpolate grid"),
}
}
row, dy, ok := interpolateIndex(rowF, g.Rows)
if !ok {
return 0, 0, 0, UnsupportedError{
Err: errors.New("cannot interpolate grid"),
}
}
i00 := row*g.Columns + col
i10 := i00 + 1
i01 := i00 + g.Columns
i11 := i01 + 1
bilinear := func(band []float32) (float64, error) {
v00 := float64(band[i00])
v10 := float64(band[i10])
v01 := float64(band[i01])
v11 := float64(band[i11])
if math.IsNaN(v00) || math.IsNaN(v10) || math.IsNaN(v01) || math.IsNaN(v11) {
return 0, UnsupportedError{
Err: errors.New("cannot interpolate grid"),
}
}
return (1-dx)*(1-dy)*v00 + dx*(1-dy)*v10 + (1-dx)*dy*v01 + dx*dy*v11, nil
}
east, err = bilinear(g.East)
if err != nil {
return 0, 0, 0, err
}
north, err = bilinear(g.North)
if err != nil {
return 0, 0, 0, err
}
up, err = bilinear(g.Up)
if err != nil {
return 0, 0, 0, err
}
return east, north, up, nil
}
func parseVelocityGridTIFF(r io.Reader) (*velocityGridData, error) {
data, err := io.ReadAll(r)
if err != nil {
return nil, err
}
p := &tifParser{data: data}
if err := p.readHeader(); err != nil {
return nil, err
}
var grids []*velocityGridData
for ifdOff := p.firstIFD; ifdOff != 0; {
g, next, err := p.readVelocityIFD(ifdOff)
if err != nil {
return nil, err
}
if g != nil {
grids = append(grids, g)
}
ifdOff = next
}
if len(grids) == 0 {
return nil, UnsupportedError{
Err: errors.New("no velocity grid directories"),
}
}
return grids[0], nil
}
func (p *tifParser) readVelocityIFD(off uint32) (*velocityGridData, uint32, error) {
if int(off)+2 > len(p.data) {
return nil, 0, fmt.Errorf("tiff: IFD out of range")
}
n := int(p.u16(int(off)))
base := int(off) + 2
need := base + n*12 + 4
if need > len(p.data) {
return nil, 0, fmt.Errorf("tiff: IFD truncated")
}
tags := map[uint16]ifdEntry{}
for i := range n {
eo := base + i*12
e := ifdEntry{
tag: p.u16(eo),
typ: p.u16(eo + 2),
count: p.u32(eo + 4),
val: p.u32(eo + 8),
}
tags[e.tag] = e
}
next := p.u32(base + n*12)
meta := parseGDALMetadata(p.tagString(tags, 42112))
typ := meta["TYPE"]
switch typ {
case "VELOCITY", "VELOCITY_CARTOGRAPHIC":
// ok
case "", "HORIZONTAL_OFFSET", "VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL", "VERTICAL_OFFSET_VERTICAL_TO_VERTICAL":
return nil, next, nil // skip non-velocity IFDs
default:
return nil, next, fmt.Errorf("tiff: unsupported velocity TYPE %q", typ)
}
width := int(p.tagLong(tags, 256))
height := int(p.tagLong(tags, 257))
if width <= 0 || height <= 0 {
return nil, next, fmt.Errorf("tiff: invalid dimensions %dx%d", width, height)
}
samples := int(p.tagLong(tags, 277))
if samples < 3 {
return nil, next, fmt.Errorf("tiff: velocity grid needs 3 samples, got %d", samples)
}
bands, err := p.readBands(tags, width, height, samples)
if err != nil {
return nil, next, err
}
tie := p.tagDoubles(tags, 33922, 6)
scale := p.tagDoubles(tags, 33550, 3)
if len(tie) < 6 || len(scale) < 2 {
return nil, next, fmt.Errorf("tiff: missing ModelTiepointTag/ModelPixelScaleTag")
}
dLon := scale[0]
dLat := math.Abs(scale[1])
if dLon == 0 || dLat == 0 {
return nil, next, fmt.Errorf("tiff: zero pixel scale")
}
return &velocityGridData{
LonUL: tie[3],
LatUL: tie[4],
DLon: dLon,
DLat: dLat,
Columns: width,
Rows: height,
East: bands[0],
North: bands[1],
Up: bands[2],
}, next, nil
}
package crs
import (
"errors"
"fmt"
"io"
"math"
"path"
"strings"
"sync"
)
// VerticalGrid is a geoid / vertical undulation grid (PROJ vgridshift / GTG).
// Grid samples are the geoid undulation N (metres), matching PROJ egm08-style
// grids used with +inv +proj=vgridshift: H = h − N, h = H + N.
type VerticalGrid string
func (vg VerticalGrid) String() string {
return fmt.Sprintf("operation=vertical_grid grid=%s", string(vg))
}
// ToTarget implements [Operation]: orthometric → ellipsoidal (h = H + N).
func (vg VerticalGrid) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
grid, err := loadVerticalGrid(string(vg))
if err != nil {
return 0, 0, 0, err
}
n, err := grid.Undulation(lon, lat)
if err != nil {
return 0, 0, 0, err
}
return lon, lat, h + n, nil
}
// FromTarget implements [Operation]: ellipsoidal → orthometric (H = h − N).
func (vg VerticalGrid) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
grid, err := loadVerticalGrid(string(vg))
if err != nil {
return 0, 0, 0, err
}
n, err := grid.Undulation(lon0, lat0)
if err != nil {
return 0, 0, 0, err
}
return lon0, lat0, h0 - n, nil
}
type verticalGridData struct {
File string
LonUL float64 // degrees east at pixel (0,0) center
LatUL float64 // degrees north at pixel (0,0) center
DLon float64 // degrees per column (east positive)
DLat float64 // degrees per row magnitude (southward in TIFF)
Columns int
Rows int
Values []float32 // row-major, TIFF order (row 0 = north)
}
var verticalGridStore sync.Map
func loadVerticalGrid(name string) (grid *verticalGridData, err error) {
name = gridFilename(name)
if v, ok := verticalGridStore.Load(name); ok {
return v.(*verticalGridData), nil
}
defer func() {
if err == nil {
grid.File = name
verticalGridStore.Store(name, grid)
}
}()
r, err := openGridReader(name)
if err != nil {
return nil, UnsupportedError{
Err: fmt.Errorf("grid: %s", name),
}
}
defer r.Close() //nolint:errcheck
switch strings.ToLower(path.Ext(name)) {
case ".tif":
return parseVerticalGridTIFF(r)
default:
return nil, fmt.Errorf("crs: unsupported vertical grid extension %q", path.Ext(name))
}
}
// Undulation returns geoid undulation N in metres at (lon, lat) degrees.
func (g *verticalGridData) Undulation(lon, lat float64) (float64, error) {
if g.Columns < 2 || g.Rows < 2 || len(g.Values) == 0 {
return 0, UnsupportedError{
Err: errors.New("invalid grid"),
}
}
colF := (lon - g.LonUL) / g.DLon
rowF := (g.LatUL - lat) / g.DLat
col, dx, ok := interpolateIndex(colF, g.Columns)
if !ok {
return 0, UnsupportedError{
Err: errors.New("cannot interpolate grid"),
}
}
row, dy, ok := interpolateIndex(rowF, g.Rows)
if !ok {
return 0, UnsupportedError{
Err: errors.New("cannot interpolate grid"),
}
}
i00 := row*g.Columns + col
i10 := i00 + 1
i01 := i00 + g.Columns
i11 := i01 + 1
v00 := float64(g.Values[i00])
v10 := float64(g.Values[i10])
v01 := float64(g.Values[i01])
v11 := float64(g.Values[i11])
if math.IsNaN(v00) || math.IsNaN(v10) || math.IsNaN(v01) || math.IsNaN(v11) {
return 0, UnsupportedError{
Err: errors.New("cannot interpolate grid"),
}
}
n := (1-dx)*(1-dy)*v00 + dx*(1-dy)*v10 + (1-dx)*dy*v01 + dx*dy*v11
return n, nil
}
// parseVerticalGridTIFF reads a PROJ vertical/geoid GeoTIFF (GTG).
func parseVerticalGridTIFF(r io.Reader) (*verticalGridData, error) {
data, err := io.ReadAll(r)
if err != nil {
return nil, err
}
p := &tifParser{data: data}
if err := p.readHeader(); err != nil {
return nil, err
}
var grids []*verticalGridData
for ifdOff := p.firstIFD; ifdOff != 0; {
g, next, err := p.readVerticalIFD(ifdOff)
if err != nil {
return nil, err
}
if g != nil {
grids = append(grids, g)
}
ifdOff = next
}
if len(grids) == 0 {
return nil, fmt.Errorf("tiff: no vertical grid directories")
}
return grids[0], nil
}
func (p *tifParser) readVerticalIFD(off uint32) (*verticalGridData, uint32, error) {
if int(off)+2 > len(p.data) {
return nil, 0, fmt.Errorf("tiff: IFD out of range")
}
n := int(p.u16(int(off)))
base := int(off) + 2
need := base + n*12 + 4
if need > len(p.data) {
return nil, 0, fmt.Errorf("tiff: IFD truncated")
}
tags := map[uint16]ifdEntry{}
for i := range n {
eo := base + i*12
e := ifdEntry{
tag: p.u16(eo),
typ: p.u16(eo + 2),
count: p.u32(eo + 4),
val: p.u32(eo + 8),
}
tags[e.tag] = e
}
next := p.u32(base + n*12)
meta := parseGDALMetadata(p.tagString(tags, 42112))
typ := meta["TYPE"]
switch typ {
case "", "VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL", "VERTICAL_OFFSET_VERTICAL_TO_VERTICAL":
// ok
case "HORIZONTAL_OFFSET":
return nil, next, nil // skip horizontal IFDs in mixed files
default:
return nil, next, fmt.Errorf("tiff: unsupported vertical TYPE %q", typ)
}
width := int(p.tagLong(tags, 256))
height := int(p.tagLong(tags, 257))
if width <= 0 || height <= 0 {
return nil, next, fmt.Errorf("tiff: invalid dimensions %dx%d", width, height)
}
samples := max(int(p.tagLong(tags, 277)), 1)
bands, err := p.readBands(tags, width, height, samples)
if err != nil {
return nil, next, err
}
tie := p.tagDoubles(tags, 33922, 6)
scale := p.tagDoubles(tags, 33550, 3)
if len(tie) < 6 || len(scale) < 2 {
return nil, next, fmt.Errorf("tiff: missing ModelTiepointTag/ModelPixelScaleTag")
}
dLon := scale[0]
dLat := math.Abs(scale[1])
if dLon == 0 || dLat == 0 {
return nil, next, fmt.Errorf("tiff: zero pixel scale")
}
return &verticalGridData{
LonUL: tie[3],
LatUL: tie[4],
DLon: dLon,
DLat: dLat,
Columns: width,
Rows: height,
Values: bands[0],
}, next, nil
}
package crs
import (
"math"
)
type VerticalOffset struct {
Dh float64
}
func (o VerticalOffset) String() string {
return build("").addAll(
"operation", "vertical_offset",
"dh", o.Dh,
).String()
}
func (o VerticalOffset) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return lon, lat, h + o.Dh, nil
}
func (o VerticalOffset) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return lon0, lat0, h0 - o.Dh, nil
}
type VerticalOffsetAndSlope struct {
Lat0, Lon0 float64
Dh float64
SlopeLat, SlopeLon float64 // arc-seconds
}
func (o VerticalOffsetAndSlope) String() string {
return build("").addAll(
"operation", "vertical_offset_and_slope",
"lat0", o.Lat0,
"lon0", o.Lon0,
"dh", o.Dh,
"slope_lat", o.SlopeLat,
"slope_lon", o.SlopeLon,
).String()
}
func (o VerticalOffsetAndSlope) delta(s Spheroid, lon, lat float64) float64 {
phi0 := radian(o.Lat0)
lam0 := radian(o.Lon0)
phi := radian(lat)
lam := radian(lon)
sinPhi0 := math.Sin(phi0)
e2 := s.E2()
oneMe2Sin2 := 1 - e2*sinPhi0*sinPhi0
rho0 := s.A * (1 - e2) / (oneMe2Sin2 * math.Sqrt(oneMe2Sin2))
nu0 := s.A / math.Sqrt(oneMe2Sin2)
slopeLat := radian(o.SlopeLat / 3600)
slopeLon := radian(o.SlopeLon / 3600)
return o.Dh +
slopeLat*rho0*(phi-phi0) +
slopeLon*nu0*(lam-lam0)*math.Cos(phi)
}
// ToTarget implements [Operation].
func (o VerticalOffsetAndSlope) ToTarget(source, target Spheroid, lon, lat, h float64) (float64, float64, float64, error) {
return lon, lat, h + o.delta(source, lon, lat), nil
}
// FromTarget implements [Operation].
func (o VerticalOffsetAndSlope) FromTarget(source, target Spheroid, lon0, lat0, h0 float64) (float64, float64, float64, error) {
return lon0, lat0, h0 - o.delta(source, lon0, lat0), nil
}
package crs
import (
"math"
)
type WebMercator struct {
Lonf float64
Latf float64
Scale float64
Eastf float64
Northf float64
}
func (cs WebMercator) String() string {
return build("web_mercator").addAll(
"lonf", cs.Lonf,
"latf", cs.Latf,
"scale", cs.Scale,
"eastf", cs.Eastf,
"northf", cs.Northf,
).String()
}
func (cs WebMercator) k() float64 {
if cs.Scale == 0 {
return 1
}
return cs.Scale
}
func mercatorM(phi float64) float64 {
return math.Log(math.Tan(math.Pi/4 + phi/2))
}
func (cs WebMercator) FromGeographic(s Spheroid, lon, lat, h float64) (float64, float64, float64) {
lambda := radian(lon)
phi := radian(lat)
lambda0 := radian(cs.Lonf)
phi0 := radian(cs.Latf)
k := cs.k()
rk := k * s.A
east := cs.Eastf + rk*(lambda-lambda0)
north := cs.Northf + rk*(mercatorM(phi)-mercatorM(phi0))
return east, north, h
}
func (cs WebMercator) ToGeographic(s Spheroid, east, north, h float64) (float64, float64, float64) {
lambda0 := radian(cs.Lonf)
phi0 := radian(cs.Latf)
k := cs.k()
rk := k * s.A
psi := (north-cs.Northf)/rk + mercatorM(phi0)
phi := 2*math.Atan(math.Exp(psi)) - math.Pi/2
lambda := lambda0 + (east-cs.Eastf)/rk
return degree(lambda), degree(phi), h
}