From cb1c930e79f5a16f13fffbe131653fcf6bdeb3ff Mon Sep 17 00:00:00 2001
From: Since <ax20yhum@cip.cs.fau.de>
Date: Wed, 20 Dec 2017 13:59:21 +0100
Subject: [PATCH] Add support for other runlength encoded files
This includes local radar images. DataZ and AtZ() is introduced for 3D CAPPI
products. Data and DataZ are both slices of PlainData, stripping elevation
views and splitting concatenated layers.
A lookup catalog is neccessary, as runlength decoding requires dimension
information, that is not specified in the local picture product headers.
---
catalog.go | 30 ++++++++++++++++++++++++++++++
data.go | 26 +++++++++++++++++++++-----
header.go | 28 +++++++++++++++++++++-------
littleendian.go | 8 ++++----
radolan.go | 39 ++++++++++++++++++++++++++++-----------
runlength.go | 6 +++---
singlebyte.go | 8 ++++----
7 files changed, 111 insertions(+), 34 deletions(-)
create mode 100644 catalog.go
diff --git a/catalog.go b/catalog.go
new file mode 100644
index 0000000..a4aad6c
--- /dev/null
+++ b/catalog.go
@@ -0,0 +1,30 @@
+package radolan
+
+type spec struct {
+ px int // plain data dimensions
+ py int
+
+ dx int // data (layer) dimensions
+ dy int
+
+ rx float64 // resolution
+ ry float64
+}
+
+// local picture products do not provide dimensions in header
+var catalog = map[string]spec{
+ "OL": {200, 224, 200, 200, 2, 2}, // reflectivity (no clutter detection)
+ "OX": {200, 224, 200, 200, 1, 1}, // reflectivity (no clutter detection)
+ "PD": {200, 224, 200, 200, 1, 1}, // radial velocity
+ "PE": {200, 224, 200, 200, 2, 2}, // echotop
+ "PF": {200, 224, 200, 200, 1, 1}, // reflectivity (15 classes)
+ "PH": {200, 224, 200, 200, 1, 1}, // accumulated rainfall
+ "PL": {200, 224, 200, 200, 2, 2}, // reflectivity
+ "PM": {200, 224, 200, 200, 2, 2}, // max. reflectivity
+ "PR": {200, 224, 200, 200, 1, 1}, // radial velocity
+ "PU": {200, 2400, 200, 200, 1, 1}, // 3D radial velocity
+ "PV": {200, 224, 200, 200, 1, 1}, // radial velocity
+ "PX": {200, 224, 200, 200, 1, 1}, // reflectivity (6 classes)
+ "PY": {200, 224, 200, 200, 1, 1}, // accumulated rainfall
+ "PZ": {200, 2400, 200, 200, 2, 2}, // 3D reflectivity CAPPI
+}
diff --git a/data.go b/data.go
index eddc7d6..28e6c31 100644
--- a/data.go
+++ b/data.go
@@ -29,7 +29,7 @@ func init() {
// only comparing header characteristics.
// This method requires header data to be already written.
func (c *Composite) identifyEncoding() encoding {
- values := c.Dx * c.Dy
+ values := c.Px * c.Py
if c.level != nil {
return runlength
@@ -47,19 +47,35 @@ func (c *Composite) identifyEncoding() encoding {
// parseData parses the composite data and writes the related fields.
// This method requires header data to be already written.
func (c *Composite) parseData(reader *bufio.Reader) error {
- if c.Dx == 0 || c.Dy == 0 {
+ if c.Px == 0 || c.Py == 0 {
return newError("parseData", "parsed header data required")
}
// create Data fields
- c.Data = make([][]RVP6, c.Dy)
- for i := range c.Data {
- c.Data[i] = make([]RVP6, c.Dx)
+ c.PlainData = make([][]RVP6, c.Py)
+ for i := range c.PlainData {
+ c.PlainData[i] = make([]RVP6, c.Px)
}
return parse[c.identifyEncoding()](c, reader)
}
+// arrangeData slices plain data into its data layers or strips preceeding
+// vertical projection
+func (c *Composite) arrangeData() {
+ if c.Py%c.Dy == 0 { // multiple layers are linked downwards
+ c.DataZ = make([][][]RVP6, c.Py/c.Dy)
+ for i := range c.DataZ {
+ c.DataZ[i] = c.PlainData[c.Dy*i : c.Dy*(i+1)] // split layers
+ }
+ } else { // only use bottom most part of plain data
+ c.DataZ = [][][]RVP6{c.PlainData[c.Py-c.Dy:]} // strip elevation
+ }
+
+ c.Dz = len(c.DataZ)
+ c.Data = c.DataZ[0] // alias
+}
+
// parseUnknown performs no action and always returns an error.
func (c *Composite) parseUnknown(rd *bufio.Reader) error {
return newError("parseUnknown", "unknown encoding")
diff --git a/header.go b/header.go
index 8206e73..918681a 100644
--- a/header.go
+++ b/header.go
@@ -89,14 +89,28 @@ func (c *Composite) parseHeader(reader *bufio.Reader) error {
}
}
- // Parse Dimensions - Example: "GP 450x 450" or "BG460460" or "GP 1500x1400"
- dim := section["GP"]
- if bg, ok := section["BG"]; ok {
- dim = bg[:len(bg)/2] + "x" + bg[len(bg)/2:]
- }
+ // Parse Dimensions - Example: "GP 450x 450" or "BG460460" or "GP 1500x1400" (if defined)
+ if dim, ok := section["GP"]; ok {
+ if _, err := fmt.Sscanf(dim, "%dx%d", &c.Dy, &c.Dx); err != nil {
+ return newError("parseHeader", "could not parse dimensions (GP): "+err.Error())
+ }
+ c.Px, c.Py = c.Dx, c.Dy // composite formats do not show elevation
+
+ } else if dim, ok := section["BG"]; ok {
+ if _, err := fmt.Sscanf(dim, "%3d%3d", &c.Dy, &c.Dx); err != nil {
+ return newError("parseHeader", "could not parse dimensions (BG): "+err.Error())
+ }
+ c.Px, c.Py = c.Dx, c.Dy // composite formats do not show elevation
+
+ } else { // dimensions of local picture products not defined in header
+ v, ok := catalog[c.Product] // lookup in catalog
+ if !ok {
+ return newError("parseHeader", "no dimension information available")
+ }
- if _, err := fmt.Sscanf(dim, "%dx%d", &c.Dy, &c.Dx); err != nil {
- return newError("parseHeader", "could not parse dimensions: "+err.Error())
+ c.Px, c.Py = v.px, v.py // plain data dimensions
+ c.Dx, c.Dy = v.dx, v.dy // data layer dimensions
+ c.Rx, c.Ry = v.rx, v.ry // data resolution
}
// Parse Precision - Example: "PR E-01" or "PR E+00"
diff --git a/littleendian.go b/littleendian.go
index c03e595..7b884e9 100644
--- a/littleendian.go
+++ b/littleendian.go
@@ -7,16 +7,16 @@ import (
)
// parseLittleEndian parses the little endian encoded composite as described in [1] and [3].
-// Result are written into the previously created Data field of the composite.
+// Result are written into the previously created PlainData field of the composite.
func (c *Composite) parseLittleEndian(reader *bufio.Reader) error {
- last := len(c.Data) - 1
- for i := range c.Data {
+ last := len(c.PlainData) - 1
+ for i := range c.PlainData {
line, err := c.readLineLittleEndian(reader)
if err != nil {
return err
}
- err = c.decodeLittleEndian(c.Data[last-i], line) // write vertically flipped
+ err = c.decodeLittleEndian(c.PlainData[last-i], line) // write vertically flipped
if err != nil {
return err
}
diff --git a/radolan.go b/radolan.go
index 5262b5f..712ed4c 100644
--- a/radolan.go
+++ b/radolan.go
@@ -29,13 +29,16 @@ import (
"time"
)
-// Radolan radar data is provided in so called composite formats. Each composite is a combined
-// image consisting of mulitiple radar sweeps spread over the composite area.
-// The composite c has a an internal resolution of c.Dx (horizontal) * c.Dy (vertical) records
-// covering a real surface of c.Dx * c.Rx * c.Dy * c.Dy square kilometers.
-// The pixel value at the position (x, y) is represented by c.Data[ y ][ x ] and is stored as
-// raw rvp-6 value (NaN if the no-data flag is set). This rvp-6 value is used differently
-// depending on the product type:
+// Radolan radar data is provided as single local sweeps or so called composite
+// formats. Each composite is a combined image consisting of mulitiple radar
+// sweeps spread over the composite area.
+// The 2D composite c has a an internal resolution of c.Dx (horizontal) * c.Dy
+// (vertical) records covering a real surface of c.Dx * c.Rx * c.Dy * c.Dy
+// square kilometers.
+// The pixel value at the position (x, y) is represented by
+// c.Data[ y ][ x ] and is stored as raw float value (called rvp-6) (NaN if the
+// no-data flag is set). Some 3D radar products feature multiple layers in
+// which the voxel at position (x, y, z) is accessible by c.DataZ[ z ][ y ][ x ].
//
// The rvp-6 value is used differently depending on the product type:
//
@@ -69,10 +72,16 @@ type Composite struct {
ForecastTime time.Time // data represents conditions predicted for this time
Interval time.Duration // time duration until next forecast
- Data [][]RVP6 // rvp-6 data for each point [y][x]
+ PlainData [][]RVP6 // rvp-6 data for parsed plain data element [y][x]
+ Px int // plain data width
+ Py int // plain data height
+
+ DataZ [][][]RVP6 // rvp-6 data for each voxel [z][y][x] (composites use only one z-layer)
+ Data [][]RVP6 // rvp-6 data for each pixel [y][x] at layer 0 (alias for DataZ[0][x][y])
Dx int // data width
Dy int // data height
+ Dz int // data layer
Rx float64 // horizontal resolution in km/px
Ry float64 // vertical resolution in km/px
@@ -81,7 +90,7 @@ type Composite struct {
dataLength int // length of binary section in bytes
- precision int // multiplicator 10^precision for each raw value
+ precision int // multiplicator 10^precision for each raw value
level []RVP6 // maps data value to corresponding rvp-6 value in runlength based formats
offx float64 // horizontal projection offset
@@ -102,6 +111,7 @@ func NewComposite(rd io.Reader) (comp *Composite, err error) {
if err != nil {
return
}
+ comp.arrangeData()
comp.calibrateProjection()
@@ -153,10 +163,17 @@ func NewDummy(product string, dx, dy int) (comp *Composite) {
// the given point. NaN is returned, if no data is available or the requested point is located
// outside the scanned area.
func (c *Composite) At(x, y int) RVP6 {
- if x < 0 || y < 0 || x >= c.Dx || y >= c.Dy {
+ return c.AtZ(x, y, 0)
+}
+
+// AtZ is shorthand for c.DataZ[z][y][x] and returns the radar video processor value (rvp-6) at
+// the given point. NaN is returned, if no data is available or the requested point is located
+// outside the scanned volume.
+func (c *Composite) AtZ(x, y, z int) RVP6 {
+ if x < 0 || y < 0 || z < 0 || x >= c.Dx || y >= c.Dy || z >= c.Dz {
return RVP6(math.NaN())
}
- return c.Data[y][x]
+ return c.DataZ[z][y][x]
}
// newError returns an error indicating the failed function and reason
diff --git a/runlength.go b/runlength.go
index 6d6e654..59f1dcc 100644
--- a/runlength.go
+++ b/runlength.go
@@ -6,15 +6,15 @@ import (
)
// parseRunlength parses the runlength encoded composite and writes into the
-// previously created Data field of the composite.
+// previously created PlainData field of the composite.
func (c *Composite) parseRunlength(reader *bufio.Reader) error {
- for i := range c.Data {
+ for i := range c.PlainData {
line, err := c.readLineRunlength(reader)
if err != nil {
return err
}
- err = c.decodeRunlength(c.Data[i], line)
+ err = c.decodeRunlength(c.PlainData[i], line)
if err != nil {
return err
}
diff --git a/singlebyte.go b/singlebyte.go
index 567b06d..1570ca2 100644
--- a/singlebyte.go
+++ b/singlebyte.go
@@ -7,16 +7,16 @@ import (
)
// parseSingleByte parses the single byte encoded composite as described in [1] and writes
-// into the previously created Data field of the composite.
+// into the previously created PlainData field of the composite.
func (c *Composite) parseSingleByte(reader *bufio.Reader) error {
- last := len(c.Data) - 1
- for i := range c.Data {
+ last := len(c.PlainData) - 1
+ for i := range c.PlainData {
line, err := c.readLineSingleByte(reader)
if err != nil {
return err
}
- err = c.decodeSingleByte(c.Data[last-i], line) // write vertically flipped
+ err = c.decodeSingleByte(c.PlainData[last-i], line) // write vertically flipped
if err != nil {
return err
}
--
GitLab