From 1fa92e18c81ca954fa151d9ffba57b2aecc9efb1 Mon Sep 17 00:00:00 2001
From: Since <ax20yhum@cip.cs.fau.de>
Date: Mon, 13 Nov 2017 18:59:40 +0100
Subject: [PATCH] Fix grid identification
An identificated grid is now indicated by the HasProjection flag.
---
radolan.go | 3 +++
translate.go | 53 ++++++++++++++++++++++++++++++++++++++++++++--------
2 files changed, 48 insertions(+), 8 deletions(-)
diff --git a/radolan.go b/radolan.go
index 44244cb..068c802 100644
--- a/radolan.go
+++ b/radolan.go
@@ -52,6 +52,7 @@ import (
// dBZ = 10 * log(Z)
// Real world geographical coordinates (latitude, longitude) can be projected into the
// coordinate system of the composite by using the translation method:
+// // if c.HasProjection
// x, y := c.Translate(52.51861, 13.40833) // Berlin (lat, lon)
//
// rvp := c.At(int(x), int(y)) // Raw value (rvp-6)
@@ -75,6 +76,8 @@ type Composite struct {
Rx float64 // horizontal resolution in km/px
Ry float64 // vertical resolution in km/px
+ HasProjection bool // coordinate translation available
+
dataLength int // length of binary section in bytes
precision int // multiplicator 10^precision for each raw value
diff --git a/translate.go b/translate.go
index 9c9e864..cfe635e 100644
--- a/translate.go
+++ b/translate.go
@@ -12,10 +12,21 @@ const (
junctionEast = 10.0 // E
)
-// cornerPoints returns corner coordinates of the national or extended european grid
+// 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
// based on the product label or resolution of the composite. The used values are
// described in [1], [4] and [5].
-func (c *Composite) cornerPoints() (originTop, originLeft, edgeBottom, edgeRight float64) {
+// If an error is returned, translation methods will not work.
+func (c *Composite) cornerPoints() (originTop, originLeft, edgeBottom, edgeRight float64, err error) {
// national grid (pg) values described in [4]
if c.Product == "PG" {
originTop, originLeft = 54.6547, 01.9178 // N, E
@@ -24,22 +35,42 @@ func (c *Composite) cornerPoints() (originTop, originLeft, edgeBottom, edgeRight
}
// national grid values described in [1]
- if c.Dx == c.Dy {
+ if c.isScaled(900, 900) {
originTop, originLeft = 54.5877, 02.0715 // N, E
edgeBottom, edgeRight = 47.0705, 14.6209 // N, E
return
}
- // extended european grid described in [5]
- originTop, originLeft = 56.5423, -0.8654 // N, E
- edgeBottom, edgeRight = 43.8736, 18.2536 // N, E
+ // 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
return
}
// calibrateProjection initializes fields that are necessary for coordinate translation
func (c *Composite) calibrateProjection() {
// get corner points
- originTop, originLeft, edgeBottom, edgeRight := c.cornerPoints()
+ 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
// set resolution to 1 km for calibration
c.Rx, c.Ry = 1.0, 1.0
@@ -54,8 +85,14 @@ func (c *Composite) calibrateProjection() {
}
// Translate translates geographical coordinates (latitude north, longitude east) to the
-// according data indices in the coordinate system of the composite. Procedures adapted from [1].
+// according data indices in the coordinate system of the composite.
+// NaN is returned when no projection is available. Procedures adapted from [1].
func (c *Composite) Translate(north, east float64) (x, y float64) {
+ if !c.HasProjection {
+ x, y = math.NaN(), math.NaN()
+ return
+ }
+
rad := func(deg float64) float64 {
return deg * (math.Pi / 180.0)
}
--
GitLab