translate.go 3.16 KB
Newer Older
Jonny Schäfer's avatar
Jonny Schäfer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
package radolan

import (
	"math"
)

// values described in [1]
const (
	earthRadius = 6370.04 // km

	junctionNorth = 60.0 // N
	junctionEast  = 10.0 // E
)

Jonny Schäfer's avatar
Jonny Schäfer committed
15
16
17
18
19
20
21
22
23
24
25
// isScaled returns true if the composite can be scaled to the given dimensions
// while also maintaining the aspect ratio.
func (c *Composite) isScaled(dx, dy int) bool {
	epsilon := 0.00000001
	return math.Abs(1-float64(c.Dy*dx)/float64(c.Dx*dy)) < epsilon
}

// errNoProjection means that the projection grid could not be identified.
var errNoProjection = newError("cornerPoints", "warning: unable to identify grid")

// cornerPoints returns corner coordinates of the national, extended or middle-european grid
Jonny Schäfer's avatar
Jonny Schäfer committed
26
27
// based on the product label or resolution of the composite. The used values are
// described in [1], [4] and [5].
Jonny Schäfer's avatar
Jonny Schäfer committed
28
29
// If an error is returned, translation methods will not work.
func (c *Composite) cornerPoints() (originTop, originLeft, edgeBottom, edgeRight float64, err error) {
Jonny Schäfer's avatar
Jonny Schäfer committed
30
31
	// national grid (pg) values described in [4]
	if c.Product == "PG" {
32
33
		originTop, originLeft = 54.66218275, 1.900684377 // N, E
		edgeBottom, edgeRight = 46.98044293, 14.73300934 // N, E
Jonny Schäfer's avatar
Jonny Schäfer committed
34
35
36
37
		return
	}

	// national grid values described in [1]
Jonny Schäfer's avatar
Jonny Schäfer committed
38
	if c.isScaled(900, 900) {
Jonny Schäfer's avatar
Jonny Schäfer committed
39
40
41
42
43
		originTop, originLeft = 54.5877, 02.0715 // N, E
		edgeBottom, edgeRight = 47.0705, 14.6209 // N, E
		return
	}

Jonny Schäfer's avatar
Jonny Schäfer committed
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
	// extended national grid described in [5]
	if c.isScaled(900, 1100) {
		originTop, originLeft = 55.5482, 03.0889 // N, E
		edgeBottom, edgeRight = 46.1827, 15.4801 // N, E
		return
	}

	// middle european grid described in [5]
	if c.isScaled(1400, 1500) {
		originTop, originLeft = 56.5423, -0.8654 // N, E
		edgeBottom, edgeRight = 43.8736, 18.2536 // N, E
		return
	}

	err = errNoProjection
Jonny Schäfer's avatar
Jonny Schäfer committed
59
60
61
62
63
64
	return
}

// calibrateProjection initializes fields that are necessary for coordinate translation
func (c *Composite) calibrateProjection() {
	// get corner points
Jonny Schäfer's avatar
Jonny Schäfer committed
65
66
67
68
69
70
71
72
73
	originTop, originLeft, edgeBottom, edgeRight, err := c.cornerPoints()
	if err != nil {
		nan := math.NaN()
		c.Rx, c.Ry, c.offx, c.offy = nan, nan, nan, nan
		return
	}

	// found matching projection rule
	c.HasProjection = true
Jonny Schäfer's avatar
Jonny Schäfer committed
74
75
76
77
78
79
80
81
82

	// set resolution to 1 km for calibration
	c.Rx, c.Ry = 1.0, 1.0

	// calibrate offset correction
	c.offx, c.offy = c.Translate(originTop, originLeft)

	// calibrate scaling
	resx, resy := c.Translate(edgeBottom, edgeRight)
Jonny Schäfer's avatar
Jonny Schäfer committed
83
84
	c.Rx = (resx) / float64(c.Dx)
	c.Ry = (resy) / float64(c.Dy)
Jonny Schäfer's avatar
Jonny Schäfer committed
85
86
87
}

// Translate translates geographical coordinates (latitude north, longitude east) to the
Jonny Schäfer's avatar
Jonny Schäfer committed
88
89
// according data indices in the coordinate system of the composite.
// NaN is returned when no projection is available. Procedures adapted from [1].
Jonny Schäfer's avatar
Jonny Schäfer committed
90
func (c *Composite) Translate(north, east float64) (x, y float64) {
Jonny Schäfer's avatar
Jonny Schäfer committed
91
92
93
94
95
	if !c.HasProjection {
		x, y = math.NaN(), math.NaN()
		return
	}

Jonny Schäfer's avatar
Jonny Schäfer committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
	rad := func(deg float64) float64 {
		return deg * (math.Pi / 180.0)
	}

	lamda0, phi0 := rad(junctionEast), rad(junctionNorth)
	lamda, phi := rad(east), rad(north)

	m := (1.0 + math.Sin(phi0)) / (1.0 + math.Sin(phi))
	x = (earthRadius * m * math.Cos(phi) * math.Sin(lamda-lamda0))
	y = (earthRadius * m * math.Cos(phi) * math.Cos(lamda-lamda0))

	// offset correction
	x -= c.offx
	y -= c.offy

	// scaling
	x /= c.Rx
	y /= c.Ry

	return
}