translate.go 3.95 KB
Newer Older
Jonny Schäfer's avatar
Jonny Schäfer committed
1
2
3
4
5
6
package radolan

import (
	"math"
)

7
8
9
10
11
12
13
14
15
16
17
// corner coordinates for projection are defined by the grid type
type grid int

const (
	unknownGrid          grid = iota
	nationalGrid              // resolution: 900km * 900km
	nationalPictureGrid       // resolution: 920km * 920km
	extendedNationalGrid      // resolution: 900km * 1100km
	middleEuropeanGrid        // resolution: 1400km * 1500km
)

Jonny Schäfer's avatar
Jonny Schäfer committed
18
19
20
21
22
23
24
25
// values described in [1]
const (
	earthRadius = 6370.04 // km

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

26
27
28
29
30
31
32
33
34
35
36
37
38
// minRes repeatedly bisects the given edges until no further step is possible
// for at least one edge. The resulting dimensions are the returned.
func minRes(dx, dy int) (rdx int, rdy int) {
	rdx, rdy = dx, dy

	if rdx == 0 || rdy == 0 {
		return
	}
	for rdx&1 == 0 && rdy&1 == 0 {
		rdx >>= 1
		rdy >>= 1
	}
	return
Jonny Schäfer's avatar
Jonny Schäfer committed
39
40
41
42
43
}

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

44
45
46
47
48
49
50
51
52
// detectGrid identifies the used projection grid based on the composite dimensions
func (c *Composite) detectGrid() grid {
	dx, dy := minRes(c.Dx, c.Dy)

	if mx, my := minRes(900, 900); dx == mx && dy == my {
		return nationalGrid
	}
	if mx, my := minRes(920, 920); dx == mx && dy == my {
		return nationalPictureGrid
Jonny Schäfer's avatar
Jonny Schäfer committed
53
	}
54
55
56
57
58
59
60
61
	if mx, my := minRes(900, 1100); dx == mx && dy == my {
		return extendedNationalGrid
	}
	if mx, my := minRes(1400, 1500); dx == mx && dy == my {
		return middleEuropeanGrid
	}
	return unknownGrid
}
Jonny Schäfer's avatar
Jonny Schäfer committed
62

63
64
// cornerPoints returns corner coordinates of the national, extended or
// middle-european grid based on the dimensions of the composite. The used
Jonny Schäfer's avatar
Jonny Schäfer committed
65
// values are described in [1], [4].  If an error is returned,
66
67
68
69
// translation methods will not work.
func (c *Composite) cornerPoints() (originTop, originLeft, edgeBottom, edgeRight float64, err error) {
	switch c.detectGrid() {
	case nationalGrid: // described in [1]
Jonny Schäfer's avatar
Jonny Schäfer committed
70
71
		originTop, originLeft = 54.5877, 02.0715 // N, E
		edgeBottom, edgeRight = 47.0705, 14.6209 // N, E
72
73
74
	case nationalPictureGrid: // (pg) described in [4]
		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
75
	case extendedNationalGrid: // described in [1]
Jonny Schäfer's avatar
Jonny Schäfer committed
76
77
		originTop, originLeft = 55.5482, 03.0889 // N, E
		edgeBottom, edgeRight = 46.1827, 15.4801 // N, E
Jonny Schäfer's avatar
Jonny Schäfer committed
78
	case middleEuropeanGrid: // described in [1]
Jonny Schäfer's avatar
Jonny Schäfer committed
79
80
		originTop, originLeft = 56.5423, -0.8654 // N, E
		edgeBottom, edgeRight = 43.8736, 18.2536 // N, E
81
82
	default:
		err = errNoProjection
Jonny Schäfer's avatar
Jonny Schäfer committed
83
	}
Jonny Schäfer's avatar
Jonny Schäfer committed
84
85
86
87
88
89
	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
90
91
92
93
94
95
96
97
98
	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
99
100
101
102
103
104
105
106
107

	// 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
108
109
	c.Rx = (resx) / float64(c.Dx)
	c.Ry = (resy) / float64(c.Dy)
Jonny Schäfer's avatar
Jonny Schäfer committed
110
111
112
}

// Translate translates geographical coordinates (latitude north, longitude east) to the
Jonny Schäfer's avatar
Jonny Schäfer committed
113
114
// 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
115
func (c *Composite) Translate(north, east float64) (x, y float64) {
Jonny Schäfer's avatar
Jonny Schäfer committed
116
117
118
119
120
	if !c.HasProjection {
		x, y = math.NaN(), math.NaN()
		return
	}

Jonny Schäfer's avatar
Jonny Schäfer committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
	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
}