Commit cb1c930e authored by Jonny Schäfer's avatar Jonny Schäfer
Browse files

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.
parent 6d2cdacc
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
}
......@@ -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")
......
......@@ -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"
......
......@@ -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
}
......
......@@ -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
......
......@@ -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
}
......
......@@ -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
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment