Initial commit

LICENSE

+MIT License
+
+Copyright (c) 2016 Jonny Schäfer
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

+# radolan
+### Go package radolan parses the DWD RADOLAN / RADVOR radar composite format.
+This format is used by the [DWD](http://www.dwd.de/DE/leistungen/radolan/radolan.html)
+for weather radar data.
+
+The obtained results can be processed and visualized with additional functions.
+The example program `radolan2png` is included to quickly convert composite files to png images.
+
+Currently the national and the extended european grids are supported.
+Tested input products are PG, FZ, SF, RW, RX and EX. Those can be considered working with
+sufficient accuracy.
+
+### Documentation
+Documentation is included in the corresponding source files and also available at
+https://godoc.org/gitlab.cs.fau.de/since/radolan
+
+### Installation
+```
+mkdir -p ~/go/src ~/go/pkg ~/go/bin
+GOPATH="~/go" go get gitlab.cs.fau.de/since/radolan/radolan2png
+```
+
+### Sample image
+This image shows the radar reflectivity (dBZ) captured 31.07.2016 18:50 CEST
+![alt text](https://gitlab.cs.fau.de/since/radolan/raw/master/assets/31-07-2016-1850.png)
+```Datenbasis: Deutscher Wetterdienst, Radardaten bildlich wiedergegeben```

conversion.go

+package radolan
+
+import (
+	"math"
+)
+
+// Radar reflectivity factor Z in logarithmic representation dBZ: dBZ = 10 * log(Z)
+type DBZ float64
+
+// Raw radar video processor value.
+type RVP6 float64
+
+// PrecipitationRate returns the estimated precipitation rate in mm/h for the given
+// reflectivity factor. The used Z-R relation is descibed in [6].
+func (z DBZ) PrecipitationRate() float64 {
+	return math.Pow(math.Pow(10, float64(z)/10)/256, 1/1.42)
+}
+
+// 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)
+}
+
+// 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)
+}
+
+// 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))
+}

conversion_test.go

+package radolan
+
+import (
+	"math"
+	"testing"
+)
+
+func TestConversion(t *testing.T) {
+	testcases := []struct {
+		rvp  RVP6
+		dbz  DBZ
+		rate float64
+	}{
+		{0, -32.5, 0.0001},
+		{65, 0, 0.0201},
+		{100, 17.5, 0.3439},
+		{200, 67.5, 1141.7670},
+	}
+
+	for _, test := range testcases {
+		dbz := test.rvp.ToDBZ()
+		rate := dbz.PrecipitationRate()
+		rvp := dbz.ToRVP6()
+
+		if dbz != test.dbz {
+			t.Errorf("RVP6(%f).ToDBZ() = %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)
+		}
+		if math.Abs(test.rate-rate) > 0.0001 {
+			t.Errorf("RVP6(%f).ToDBZ().PrecipitationRate() = %f; expected: %f", test.rvp, rate, test.rate)
+		}
+	}
+}

data.go

+package radolan
+
+import (
+	"bufio"
+)
+
+// encoding types of the composite
+type encoding int
+
+const (
+	runlength encoding = iota
+	littleEndian
+	singleByte
+	unknown
+)
+
+// parsing methods
+var parse = [4]func(c *Composite, rd *bufio.Reader) error{}
+
+// init maps the parsing methods to the encoding type
+func init() {
+	parse[runlength] = (*Composite).parseRunlength
+	parse[littleEndian] = (*Composite).parseLittleEndian
+	parse[singleByte] = (*Composite).parseSingleByte
+	parse[unknown] = (*Composite).parseUnknown
+}
+
+// identifyEncoing identifies the encoding type of the data section by
+// only comparing header characteristics.
+// This method requires header data to be already written.
+func (c *Composite) identifyEncoding() encoding {
+	values := c.Dx * c.Dy
+
+	if c.level != nil {
+		return runlength
+	}
+	if c.dataLength == values*2 {
+		return littleEndian
+	}
+	if c.dataLength == values {
+		return singleByte
+	}
+
+	return unknown
+}
+
+// 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 {
+		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)
+	}
+
+	return parse[c.identifyEncoding()](c, reader)
+}
+
+// parseUnknown performs no action and always returns an error.
+func (c *Composite) parseUnknown(rd *bufio.Reader) error {
+	return newError("parseUnknown", "unknown encoding")
+}

header.go

+package radolan
+
+import (
+	"bufio"
+	"fmt"
+	"strings"
+	"time"
+	"unicode"
+)
+
+// splitHeader splits the given header string into its fields. The returned
+// map is using the field name as key and the field content as value.
+func splitHeader(header string) (m map[string]string) {
+	m = make(map[string]string)
+	var beginKey, endKey, beginValue, endValue int
+	var dispatch bool
+
+	for i, c := range header {
+		if unicode.IsUpper(c) {
+			if dispatch {
+				m[header[beginKey:endKey]] = header[beginValue:endValue]
+				beginKey = i
+				dispatch = false
+			}
+			endKey = i + 1
+		} else {
+			if i == 0 {
+				return // no key prefixing value
+			}
+			if !dispatch {
+				beginValue = i
+				dispatch = true
+			}
+			endValue = i + 1
+		}
+	}
+	m[header[beginKey:endKey]] = header[beginValue:endValue]
+
+	return
+}
+
+// parseHeader parses and the composite header and writes the related fields as
+// described in [1] and [3].
+func (c *Composite) parseHeader(reader *bufio.Reader) error {
+	header, err := reader.ReadString('\x03')
+	if err != nil || len(header) < 22 { // smaller length makes no sense
+		return newError("parseHeader", "header corrupted: too short")
+	}
+
+	// Split header segments
+	section := splitHeader(header[:len(header)-1]) // without delimiter
+
+	// Parse Product - Example: "PG" or "FZ"
+	c.Product = header[:2]
+
+	// Parse DataLength - Example: "BY 405160"
+	if _, err := fmt.Sscanf(section["BY"], "%d", &c.dataLength); err != nil {
+		return newError("parseHeader", "could not parse data length: "+err.Error())
+	}
+	c.dataLength -= len(header) // remove header length including delimiter
+
+	// Parse CaptureTime - Example: "PG262115100000616" or "FZ211615100000716"
+	date := header[2:8] + header[13:17] // cut WMO number
+	c.CaptureTime, err = time.Parse("0215040106", date)
+	if err != nil {
+		return newError("parseHeader", "could not parse capture time: "+err.Error())
+	}
+
+	// Parse ForecastTime - Example: "VV 005"
+	c.ForecastTime = c.CaptureTime
+	if vv, ok := section["VV"]; ok {
+		min := 0
+		if _, err := fmt.Sscanf(vv, "%d", &min); err != nil {
+			return newError("parseHeader", "could not parse forecast time: "+err.Error())
+		}
+		c.ForecastTime = c.CaptureTime.Add(time.Duration(min) * time.Minute)
+	}
+
+	// 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:]
+	}
+
+	if _, err := fmt.Sscanf(dim, "%dx%d", &c.Dy, &c.Dx); err != nil {
+		return newError("parseHeader", "could not parse dimensions: "+err.Error())
+	}
+
+	// Parse Precision - Example: "PR E-01" or "PR E+00"
+	if prec, ok := section["E"]; ok { // not that nice
+		if _, err := fmt.Sscanf(prec, "%d", &c.precision); err != nil {
+			return newError("parseHeader", "could not parse precision: "+err.Error())
+		}
+	}
+
+	// Parse Level - Example "LV 6  1.0 19.0 28.0 37.0 46.0 55.0"
+	if lv, ok := section["LV"]; ok {
+		level := strings.Fields(lv)
+		if len(level) < 2 {
+			return newError("parseHeader", "invalid level count")
+		}
+
+		c.level = make([]DBZ, len(level)-1)
+		for i, f := range level[1:] {
+			if _, err = fmt.Sscanf(f, "%f", &c.level[i]); err != nil {
+				return newError("parseHeader", "invalid level value: "+err.Error())
+			}
+		}
+	}
+
+	return nil
+}

header_test.go

+package radolan
+
+import (
+	"bufio"
+	"strings"
+	"testing"
+	"time"
+)
+
+type headerTestcase struct {
+	// head of file
+	test string
+
+	// expected
+	expBinary       string
+	expProduct      string
+	expCaptureTime  time.Time
+	expForecastTime time.Time
+	expDx           int
+	expDy           int
+	expDataLength   int
+	expPrecision    int
+	expLevel        []DBZ
+}
+
+func TestParseHeaderPG(t *testing.T) {
+	ht := &headerTestcase{}
+	var err1, err2 error
+
+	// head of file
+	ht.test = "PG262115100000616BY22205LV 6  1.0 19.0 28.0 37.0 46.0 55.0CS0MX 0MS " +
+		"88<boo,ros,emd,hnr,umd,pro,ess,fld,drs,neu,nhb,oft,eis,tur,isn,fbg,mem> " +
+		"are used, BG460460\x03binarycontent"
+
+	// expected
+	ht.expBinary = "binarycontent"
+	ht.expProduct = "PG"
+	ht.expCaptureTime, err1 = time.Parse(time.RFC1123, "Sun, 26 Jun 2016 23:15:00 CEST")
+	ht.expForecastTime, err2 = time.Parse(time.RFC1123, "Sun, 26 Jun 2016 23:15:00 CEST")
+	ht.expDx = 460
+	ht.expDy = 460
+	ht.expDataLength = 22205 - 159 // BY - header_etx_length
+	ht.expPrecision = 0
+	ht.expLevel = []DBZ{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)
+	}
+
+	testParseHeader(t, ht)
+}
+
+func TestParseHeaderFZ(t *testing.T) {
+	ht := &headerTestcase{}
+	var err1, err2 error
+
+	// head of file
+	ht.test = "FZ282105100000716BY 405160VS 3SW   2.13.1PR E-01INT   5GP 450x 450VV 100MF " +
+		"00000002MS 66<boo,ros,emd,hnr,umd,pro,ess,drs,neu,nhb,oft,eis,tur,isn,fbg,mem>" +
+		"\x03binarycontent"
+
+	// ht.expected values
+	ht.expBinary = "binarycontent"
+
+	ht.expProduct = "FZ"
+	ht.expCaptureTime, err1 = time.Parse(time.RFC1123, "Thu, 28 Jul 2016 23:05:00 CEST")
+	ht.expForecastTime, err2 = time.Parse(time.RFC1123, "Fri, 29 Jul 2016 00:45:00 CEST")
+	ht.expDx = 450
+	ht.expDy = 450
+	ht.expDataLength = 405160 - 154 // BY - header_etx_length
+	ht.expPrecision = -1
+	ht.expLevel = []DBZ(nil)
+
+	if err1 != nil || err2 != nil {
+		t.Errorf("%s.parseHeader(): wrong testcase time.Parse", ht.expProduct)
+	}
+
+	testParseHeader(t, ht)
+}
+
+func testParseHeader(t *testing.T, ht *headerTestcase) {
+	dummy := &Composite{}
+	reader := bufio.NewReader(strings.NewReader(ht.test))
+
+	// run
+	if err := dummy.parseHeader(reader); err != nil {
+		t.Errorf("%s.parseHeader(): returned error: %#v", err.Error())
+	}
+
+	// test results
+	// Product
+	if dummy.Product != ht.expProduct {
+		t.Errorf("%s.parseHeader(): Product: %#v; expected: %#v", ht.expProduct,
+			dummy.Product, ht.expProduct)
+	}
+
+	// CaptureTime
+	if !dummy.CaptureTime.Equal(ht.expCaptureTime) {
+		t.Errorf("%s.parseHeader(): CaptureTime: %#v; expected: %#v", ht.expProduct,
+			dummy.CaptureTime.String(), ht.expCaptureTime.String())
+	}
+
+	// ForecastTime
+	if !dummy.ForecastTime.Equal(ht.expForecastTime) {
+		t.Errorf("%s.parseHeader(): ForecastTime: %#v; expected: %#v", ht.expProduct,
+			dummy.ForecastTime.String(), ht.expForecastTime.String())
+	}
+
+	// Dx Dy
+	if dummy.Dx != ht.expDx || dummy.Dy != ht.expDy {
+		t.Errorf("%s.parseHeader(): Dx: %d Dy: %d; expected Dx: %d Dy: %d", ht.expProduct,
+			dummy.Dx, dummy.Dy, ht.expDx, ht.expDy)
+	}
+
+	// dataLength
+	if dummy.dataLength != ht.expDataLength {
+		t.Errorf("%s.parseHeader(): dataLength: %#v; expected: %#v", ht.expProduct,
+			dummy.dataLength, ht.expDataLength)
+	}
+
+	// precision
+	if dummy.precision != ht.expPrecision {
+		t.Errorf("%s.parseHeader(): precision: %#v; expected: %#v", ht.expProduct,
+			dummy.precision, ht.expPrecision)
+	}
+
+	// level
+	for i := range ht.expLevel {
+		if len(dummy.level) != len(ht.expLevel) || dummy.level[i] != ht.expLevel[i] {
+			t.Errorf("%s.parseHeader(): level: %#v; expected: %#v", ht.expProduct,
+				dummy.level, ht.expLevel)
+		}
+	}
+
+	// check consistency
+	if line, _ := reader.ReadString('\n'); line != ht.expBinary {
+		t.Errorf("%s.parseHeader(): binary data corrupted", ht.expProduct)
+	}
+}

image.go

+package radolan
+
+import (
+	"image"
+	"image/color"
+	"math"
+)
+
+// A ColorFunc can be used to assign colors to rvp-6 values for image creation.
+type ColorFunc func(rvp RVP6) color.RGBA
+
+// Sample color and grayscale gradients for visualization with the image method.
+var (
+	// HeatmapReflectivity is a color gradient for cloud reflectivity composites
+	// between 5dBZ and 75dBZ.
+	HeatmapReflectivity ColorFunc
+
+	// HeatmapReflectivityWide is a color gradient for cloud reflectivity composites
+	// between -32.5dBZ and 75dBZ.
+	HeatmapReflectivityWide ColorFunc
+
+	// HeatmapAccumulatedHour is a color gradient for accumulated rainfall composites
+	// (e.g RW) between 0.1mm/h and 100 mm/h using logarithmic compression.
+	HeatmapAccumulatedHour ColorFunc
+
+	// HeatmapAccumulatedDay is a color gradient for accumulated rainfall composites
+	// (e.g. SF) between 0.1mm and 200mm using logarithmic compression.
+	HeatmapAccumulatedDay ColorFunc
+
+	// GraymapLinear is a linear grayscale gradient between the (raw) rvp-6
+	// values 0 and 409.5.
+	GraymapLinear ColorFunc
+
+	// GraymapLinearWide is a linear grayscale gradient between the (raw) rvp-6
+	// values 0 and 4095.
+	GraymapLinearWide ColorFunc
+)
+
+func init() {
+	// return identity (no compression)
+	id := func(x RVP6) RVP6 {
+		return x
+	}
+	// return logarithm
+	log := func(x RVP6) RVP6 {
+		return RVP6(math.Log(float64(x)))
+	}
+
+	HeatmapReflectivity = heatmap(DBZ(5.0).ToRVP6(), DBZ(75.0).ToRVP6(), id)
+	HeatmapReflectivityWide = heatmap(DBZ(-32.5).ToRVP6(), DBZ(75.0).ToRVP6(), id)
+	HeatmapAccumulatedHour = heatmap(0.1, 100, log)
+	HeatmapAccumulatedDay = heatmap(0.1, 200, log)
+
+	GraymapLinear = graymap(0, 409.5, id)
+	GraymapLinearWide = graymap(0, 4095, id)
+}
+
+// Image creates an image by evaluating the color function fn for each raw rvp-6 value.
+func (c *Composite) Image(fn ColorFunc) *image.RGBA {
+	rec := image.Rect(0, 0, c.Dx, c.Dy)
+	img := image.NewRGBA(rec)
+
+	for y := 0; y < c.Dy; y++ {
+		for x := 0; x < c.Dx; x++ {
+			img.Set(x, y, fn(c.Data[y][x]))
+		}
+	}
+
+	return img
+}
+
+// graymap returns a grayscale gradient between min and max. A compression function is used to
+// make logarithmic scales possible.
+func graymap(min, max RVP6, compression func(RVP6) RVP6) ColorFunc {
+	min = compression(min)
+	max = compression(max)
+
+	return func(rvp RVP6) color.RGBA {
+		rvp = compression(rvp)
+
+		if rvp < min {
+			return color.RGBA{0x00, 0x00, 0x00, 0xFF} // black
+		}
+
+		p := float64((rvp - min) / (max - min))
+		if p > 1 {
+			p = 1
+		}
+
+		l := uint8(0xFF * p)
+		return color.RGBA{l, l, l, 0xFF}
+	}
+}
+
+// heatmap returns a colour gradient between min and max. A compression function is used to
+// make logarithmic scales possible.
+func heatmap(min, max RVP6, compression func(RVP6) RVP6) ColorFunc {
+	min = compression(min)
+	max = compression(max)
+
+	return func(rvp RVP6) color.RGBA {
+		rvp = compression(rvp)
+		if rvp < min {
+			return color.RGBA{0x00, 0x00, 0x00, 0xFF} // black
+		}
+
+		p := float64((rvp - min) / (max - min))
+		if p > 1 { // limit
+			p = 1
+		}
+		h := math.Mod(360-(330*p)+240, 360)
+
+		s := 1.0          // saturation
+		l := 0.5*p + 0.25 // lightness
+
+		// adapted from https://en.wikipedia.org/wiki/HSL_and_HSV#From_HSL
+		c := (1 - math.Abs(2*l-1)) * s // calculate chroma
+
+		hh := h / 60
+		x := c * (1 - math.Abs(math.Mod(hh, 2)-1))
+
+		if math.IsNaN(hh) {
+			hh = -1
+		}
+
+		var rr, gg, bb float64
+		switch int(hh) {
+		case 0:
+			rr, gg, bb = c, x, 0
+		case 1:
+			rr, gg, bb = x, c, 0
+		case 2:
+			rr, gg, bb = 0, c, x
+		case 3:
+			rr, gg, bb = 0, x, c
+		case 4:
+			rr, gg, bb = x, 0, c
+		case 5:
+			rr, gg, bb = c, 0, x
+		}
+
+		m := l - c/2
+		r, g, b := uint8(0xFF*(rr+m)), uint8(0xFF*(gg+m)), uint8(0xFF*(bb+m))
+
+		return color.RGBA{r, g, b, 0xFF}
+	}
+}

littleendian.go

+package radolan
+
+import (
+	"bufio"
+	"io"
+	"math"
+)
+
+// 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.
+func (c *Composite) parseLittleEndian(reader *bufio.Reader) error {
+	last := len(c.Data) - 1
+	for i := range c.Data {
+		line, err := c.readLineLittleEndian(reader)
+		if err != nil {
+			return err
+		}
+
+		err = c.decodeLittleEndian(c.Data[last-i], line) // write vertically flipped
+		if err != nil {
+			return err
+		}
+	}
+
+	return nil
+}
+
+// readLineLittleEndian reads a line until horizontal limit from the given reader
+// This method is used to get a line of little endian encoded data.
+func (c *Composite) readLineLittleEndian(rd *bufio.Reader) (line []byte, err error) {
+	line = make([]byte, c.Dx*2)
+	_, err = io.ReadFull(rd, line)