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

Automatically convert values based on DataUnit field

It is no longer needed to keep track of conversions between RVP6 (which might
actually represent dBZ data) and DBZ. Those behaved differently depending on
the product type.
parent 5f567473
......@@ -12,7 +12,7 @@ type spec struct {
}
// local picture products do not provide dimensions in header
var catalog = map[string]spec{
var dimensionCatalog = 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
......
......@@ -4,11 +4,11 @@ import (
"math"
)
// Radar reflectivity factor Z in logarithmic representation dBZ: dBZ = 10 * log(Z)
type DBZ float64
var NaN = float32(math.NaN())
// Raw radar video processor value.
type RVP6 float64
func IsNaN(f float32) (is bool) {
return f != f
}
// Z-R relationship mathematically expressed as Z = a * R^b
type ZR struct {
......@@ -26,30 +26,30 @@ var (
// PrecipitationRate returns the estimated precipitation rate in mm/h for the given
// reflectivity factor and Z-R relationship.
func (z DBZ) PrecipitationRate(relation ZR) float64 {
return math.Pow(math.Pow(10, float64(z)/10)/relation.A, 1/relation.B)
func PrecipitationRate(relation ZR, dBZ float32) (rate float64) {
return math.Pow(math.Pow(10, float64(dBZ)/10)/relation.A, 1/relation.B)
}
// Reflectivity returns the estimated reflectivity factor for the given precipitation
// rate (mm/h) and Z-R relationship.
func Reflectivity(rate float64, relation ZR) DBZ {
return DBZ(10 * math.Log10(relation.A*math.Pow(rate, relation.B)))
func Reflectivity(relation ZR, rate float64) (dBZ float32) {
return float32(10 * math.Log10(relation.A*math.Pow(float64(rate), relation.B)))
}
// ToDBZ converts the given radar video processor values (rvp-6) to radar reflectivity
// toDBZ converts the given radar video processor values (rvp-6) to radar reflectivity
// factors in decibel relative to Z (dBZ).
func (r RVP6) ToDBZ() DBZ {
return DBZ(r/2.0 - 32.5)
func toDBZ(rvp6 float32) (dBZ float32) {
return rvp6/2.0 - 32.5
}
// ToRVP6 converts the given radar reflectivity factors (dBZ) to radar video processor
// toRVP6 converts the given radar reflectivity factors (dBZ) to radar video processor
// values (rvp-6).
func (z DBZ) ToRVP6() RVP6 {
return RVP6((z + 32.5) * 2)
func toRVP6(dBZ float32) float32 {
return (dBZ + 32.5) * 2
}
// rvp6Raw converts the raw value to radar video processor values (rvp-6) by applying the
// products precision field.
func (c *Composite) rvp6Raw(value int) RVP6 {
return RVP6(float64(value) * math.Pow10(c.precision))
func (c *Composite) rvp6Raw(value int) float32 {
return float32(value) * float32(math.Pow10(c.precision))
}
......@@ -7,8 +7,8 @@ import (
func TestConversion(t *testing.T) {
testcases := []struct {
rvp RVP6
dbz DBZ
rvp float32
dbz float32
zr float64
}{
{0, -32.5, 0.0001},
......@@ -18,22 +18,22 @@ func TestConversion(t *testing.T) {
}
for _, test := range testcases {
dbz := test.rvp.ToDBZ()
zr := dbz.PrecipitationRate(Aniol80)
rz := Reflectivity(zr, Aniol80)
rvp := dbz.ToRVP6()
dbz := toDBZ(test.rvp)
zr := PrecipitationRate(Aniol80, dbz)
rz := Reflectivity(Aniol80, zr)
rvp := toRVP6(dbz)
if dbz != test.dbz {
t.Errorf("RVP6(%f).ToDBZ() = %f; expected: %f", test.rvp, dbz, test.dbz)
t.Errorf("toDBZ(%f) = %f; expected: %f", test.rvp, dbz, test.dbz)
}
if rvp != test.rvp {
t.Errorf("RVP6(%f).ToDBZ().ToRVP6() = %f; expected: %f", test.rvp, rvp, test.rvp)
t.Errorf("toRVP6(toDBZ(%f)) = %f; expected: %f", test.rvp, rvp, test.rvp)
}
if math.Abs(test.zr-zr) > 0.0001 {
t.Errorf("RVP6(%f).ToDBZ().PrecipitationRate() = %f; expected: %f", test.rvp, zr, test.zr)
t.Errorf("PrecipitationRate(Aniol80, toDBZ(%f)) = %f; expected: %f", test.rvp, zr, test.zr)
}
if math.Abs(float64(test.dbz-rz)) > 0.0000001 {
t.Errorf("Reflectivity(RVP6(%f).ToDBZ().PrecipitationRate()) = %f; expected: %f",
t.Errorf("Reflectivity(PrecipitationRate(toDBZ(%f))) = %f; expected: %f",
test.rvp, rz, test.dbz)
}
}
......
......@@ -52,9 +52,9 @@ func (c *Composite) parseData(reader *bufio.Reader) error {
}
// create Data fields
c.PlainData = make([][]RVP6, c.Py)
c.PlainData = make([][]float32, c.Py)
for i := range c.PlainData {
c.PlainData[i] = make([]RVP6, c.Px)
c.PlainData[i] = make([]float32, c.Px)
}
return parse[c.identifyEncoding()](c, reader)
......@@ -64,12 +64,12 @@ func (c *Composite) parseData(reader *bufio.Reader) error {
// 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)
c.DataZ = make([][][]float32, 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.DataZ = [][][]float32{c.PlainData[c.Py-c.Dy:]} // strip elevation
}
c.Dz = len(c.DataZ)
......
......@@ -109,7 +109,7 @@ func (c *Composite) parseHeader(reader *bufio.Reader) 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
v, ok := dimensionCatalog[c.Product] // lookup in catalog
if !ok {
return newError("parseHeader", "no dimension information available")
}
......@@ -142,7 +142,7 @@ func (c *Composite) parseHeader(reader *bufio.Reader) error {
return newError("parseHeader", "invalid level format: "+lv)
}
c.level = make([]RVP6, cnt)
c.level = make([]float32, cnt)
for i := range c.level {
n := i * 5
if _, err = fmt.Sscanf(lv[n+2:n+7], "%f", &c.level[i]); err != nil {
......
......@@ -21,7 +21,7 @@ type headerTestcase struct {
expDy int
expDataLength int
expPrecision int
expLevel []RVP6
expLevel []float32
}
func TestParseHeaderPG(t *testing.T) {
......@@ -42,7 +42,7 @@ func TestParseHeaderPG(t *testing.T) {
ht.expDy = 460
ht.expDataLength = 22205 - 159 // BY - header_etx_length
ht.expPrecision = 0
ht.expLevel = []RVP6{1.0, 19.0, 28.0, 37.0, 46.0, 55.0}
ht.expLevel = []float32{1.0, 19.0, 28.0, 37.0, 46.0, 55.0}
if err1 != nil || err2 != nil {
t.Errorf("%s.parseHeader(): wrong testcase time.Parse", ht.expProduct)
......@@ -71,7 +71,7 @@ func TestParseHeaderFZ(t *testing.T) {
ht.expDy = 450
ht.expDataLength = 405160 - 154 // BY - header_etx_length
ht.expPrecision = -1
ht.expLevel = []RVP6(nil)
ht.expLevel = []float32(nil)
if err1 != nil || err2 != nil {
t.Errorf("%s.parseHeader(): wrong testcase time.Parse", ht.expProduct)
......
......@@ -3,7 +3,6 @@ package radolan
import (
"bufio"
"io"
"math"
)
// parseLittleEndian parses the little endian encoded composite as described in [1] and [3].
......@@ -37,7 +36,7 @@ func (c *Composite) readLineLittleEndian(rd *bufio.Reader) (line []byte, err err
}
// decodeLittleEndian decodes the source line and writes to the given destination.
func (c *Composite) decodeLittleEndian(dst []RVP6, line []byte) error {
func (c *Composite) decodeLittleEndian(dst []float32, line []byte) error {
if len(line)%2 != 0 || len(dst)*2 != len(line) {
return newError("decodeLittleEndian", "wrong destination or source size")
}
......@@ -52,18 +51,27 @@ func (c *Composite) decodeLittleEndian(dst []RVP6, line []byte) error {
// rvp6LittleEndian converts the raw two byte tuple of little endian encoded composite products
// to radar video processor values (rvp-6). NaN may be returned when the no-data flag is set.
func (c *Composite) rvp6LittleEndian(tuple [2]byte) RVP6 {
func (c *Composite) rvp6LittleEndian(tuple [2]byte) float32 {
var value int = 0x0F & int(tuple[1])
value = (value << 8) + int(tuple[0])
if tuple[1]&(1<<5) != 0 { // error code: no-data
return RVP6(math.NaN())
return NaN
}
if tuple[1]&(1<<6) != 0 { // flag: negative value
value *= -1
}
// set decimal point
return c.rvp6Raw(value)
conv := c.rvp6Raw(value) // set decimal point
// little endian encoded formats are also used for mm/h
if c.DataUnit != Unit_dBZ {
return conv
}
// Even though this format supports negative values and custom
// precision they do not make use of this and we still have to subtract
// the bias and scale it (RADVOR FX, dBZ)
return toDBZ(conv)
}
......@@ -24,7 +24,6 @@ import (
"compress/gzip"
"fmt"
"io"
"math"
"sort"
"time"
)
......@@ -36,18 +35,18 @@ import (
// (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 ].
// c.Data[ y ][ x ] and is stored as raw float value (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:
// The data value is used differently depending on the product type:
// (also consult the DataUnit field of the Composite)
//
// Product label | rvp-6 value represents | unit
// Product label | values represent | unit
// -------------------------+--------------------------+------------------------
// PG, PC, PX*, ... | cloud reflectivity | dBZ
// RX, WX, EX, FZ, FX, ... | raw value | convert to dBZ with ToDBZ()
// RW, SF, ... | aggregated precipitation | mm/h or mm/d
// RX, WX, EX, FZ, FX, ... | cloud reflectivity | dBZ
// RW, SF, ... | aggregated precipitation | mm/interval
// PR*, ... | doppler radial velocity | m/s
//
// The cloud reflectivity (in dBZ) can be converted to rainfall rate (in mm/h)
......@@ -60,9 +59,8 @@ import (
// // 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)
// dbz := rvp.ToDBZ() // Cloud reflectivity (dBZ)
// rat := dbz.PrecipitationRate(radolan.Doelling98) // Rainfall rate (mm/h) using Doelling98 as Z-R relationship
// dbz := c.At(int(x), int(y)) // Raw value is Cloud reflectivity (dBZ)
// rat := radolan.PrecipitationRate(radolan.Doelling98, dbz) // Rainfall rate (mm/h) using Doelling98 as Z-R relationship
//
// fmt.Println("Rainfall in Berlin [mm/h]:", rat)
//
......@@ -73,13 +71,14 @@ type Composite struct {
ForecastTime time.Time // data represents conditions predicted for this time
Interval time.Duration // time duration until next forecast
PlainData [][]RVP6 // rvp-6 data for parsed plain data element [y][x]
Px int // plain data width
Py int // plain data height
DataUnit Unit
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])
PlainData [][]float32 // data for parsed plain data element [y][x]
Px int // plain data width
Py int // plain data height
DataZ [][][]float32 // data for each voxel [z][y][x] (composites use only one z-layer)
Data [][]float32 // data for each pixel [y][x] at layer 0 (alias for DataZ[0][x][y])
Dx int // data width
Dy int // data height
......@@ -92,14 +91,23 @@ type Composite struct {
dataLength int // length of binary section in bytes
precision int // multiplicator 10^precision for each raw value
level []RVP6 // maps data value to corresponding rvp-6 value in runlength based formats
precision int // multiplicator 10^precision for each raw value
level []float32 // maps data value to corresponding index value in runlength based formats
offx float64 // horizontal projection offset
offy float64 // vertical projection offset
}
// NewComposite reads binary data from rd and parses the composite.
// ErrUnknownUnit indicates that the unit of the radar data is not defined in
// the catalog. The data values can be incorrect due to unit dependent
// conversions during parsing. Be careful when further processing the
// composite.
var ErrUnknownUnit = newError("NewComposite", "data unit not defined in catalog. data values can be incorrect")
// NewComposite reads binary data from rd and parses the composite. An error
// is returned on failure. When ErrUnknownUnit is returned, the data values can
// be incorrect due to unit dependent conversions during parsing. In this case
// be careful when further processing the composite.
func NewComposite(rd io.Reader) (comp *Composite, err error) {
reader := bufio.NewReader(rd)
comp = &Composite{}
......@@ -117,6 +125,10 @@ func NewComposite(rd io.Reader) (comp *Composite, err error) {
comp.calibrateProjection()
if comp.DataUnit == Unit_unknown {
err = ErrUnknownUnit
}
return
}
......@@ -161,19 +173,19 @@ func NewDummy(product string, dx, dy int) (comp *Composite) {
return
}
// At is shorthand for c.Data[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 area.
func (c *Composite) At(x, y int) RVP6 {
// At is shorthand for c.Data[y][x] and returns the radar video processor value
// at 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) float32 {
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 {
// AtZ is shorthand for c.DataZ[z][y][x] and returns the radar video processor
// value 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) float32 {
if x < 0 || y < 0 || z < 0 || x >= c.Dx || y >= c.Dy || z >= c.Dz {
return RVP6(math.NaN())
return NaN
}
return c.DataZ[z][y][x]
}
......
......@@ -2,7 +2,6 @@ package radolan
import (
"bufio"
"math"
)
// parseRunlength parses the runlength encoded composite and writes into the
......@@ -38,11 +37,10 @@ func (c *Composite) readLineRunlength(rd *bufio.Reader) (line []byte, err error)
}
// decodeRunlength decodes the source line and writes to the given destination.
func (c *Composite) decodeRunlength(dst []RVP6, line []byte) error {
func (c *Composite) decodeRunlength(dst []float32, line []byte) error {
// fill destination as runlength encoding will induce gaps
nan := RVP6(math.NaN())
for i := range dst {
dst[i] = nan
dst[i] = NaN
}
dstpos := 0
......@@ -78,14 +76,14 @@ func (c *Composite) decodeRunlength(dst []RVP6, line []byte) error {
// rvp6Runlength sets the value of level based composite products to radar
// video processor values (rvp-6).
func (c *Composite) rvp6Runlength(value byte) RVP6 {
func (c *Composite) rvp6Runlength(value byte) float32 {
if value == 0 {
return RVP6(math.NaN())
return NaN
}
value--
if int(value) >= len(c.level) { // border markings
return RVP6(math.NaN())
return NaN
}
return c.level[value]
}
......@@ -3,7 +3,6 @@ package radolan
import (
"bufio"
"io"
"math"
)
// parseSingleByte parses the single byte encoded composite as described in [1] and writes
......@@ -37,7 +36,7 @@ func (c *Composite) readLineSingleByte(rd *bufio.Reader) (line []byte, err error
}
// decodeSingleByte decodes the source line and writes to the given destination.
func (c *Composite) decodeSingleByte(dst []RVP6, line []byte) error {
func (c *Composite) decodeSingleByte(dst []float32, line []byte) error {
if len(dst) != len(line) {
return newError("decodeSingleByte", "wrong destination or source size")
}
......@@ -52,11 +51,17 @@ func (c *Composite) decodeSingleByte(dst []RVP6, line []byte) error {
// rvp6SingleByte converts the raw byte of single byte encoded
// composite products to radar video processor values (rvp-6). NaN may be returned
// when the no-data flag is set.
func (c *Composite) rvp6SingleByte(value byte) RVP6 {
func (c *Composite) rvp6SingleByte(value byte) float32 {
if value == 250 { // error code: no-data
return RVP6(math.NaN())
return NaN
}
// set decimal point
return c.rvp6Raw(int(value))
conv := c.rvp6Raw(int(value)) // set decimal point
// not sure if single byte formats are even used for other things than dBZ (RX, dBZ)
if c.DataUnit != Unit_dBZ {
return conv
}
return toDBZ(conv)
}
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