radolan.go 6.4 KB
Newer Older
Jonny Schäfer's avatar
Jonny Schäfer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// Package radolan parses the DWD RADOLAN / RADVOR radar composite format. This data
// is available at the Global Basic Dataset (http://www.dwd.de/DE/leistungen/gds/gds.html).
// The obtained results can be processed and visualized with additional functions.
//
// Currently the national grid [1][4] and the extended european grid [5] are supported.
// Tested input products are PG, FZ, SF, RW, RX and EX. Those can be considered working with
// sufficient accuracy.
//
// In cases, where the publicly available format specification is unprecise or contradictory,
// reverse engineering was used to obtain reasonable approaches.
//	Used references:
//
//	[1] https://www.dwd.de/DE/leistungen/radolan/radolan_info/radolan_radvor_op_komposit_format_pdf.pdf
//	[2] https://www.dwd.de/DE/leistungen/gds/weiterfuehrende_informationen.zip
//	[3]  - legend_radar_products_fz_forecast.pdf
//	[4]  - legend_radar_products_pg_coordinates.pdf
//	[5]  - legend_radar_products_radolan_rw_sf.pdf
//	[6] https://www.dwd.de/DE/leistungen/radarniederschlag/rn_info/download_niederschlagsbestimmung.pdf
package radolan

import (
Jonny Schäfer's avatar
Jonny Schäfer committed
22
	"archive/tar"
Jonny Schäfer's avatar
Jonny Schäfer committed
23
	"bufio"
Jonny Schäfer's avatar
Jonny Schäfer committed
24
	"compress/gzip"
Jonny Schäfer's avatar
Jonny Schäfer committed
25
26
27
	"fmt"
	"io"
	"math"
Jonny Schäfer's avatar
Jonny Schäfer committed
28
	"sort"
Jonny Schäfer's avatar
Jonny Schäfer committed
29
30
31
	"time"
)

32
33
34
35
36
37
38
39
40
41
// 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 ].
Jonny Schäfer's avatar
Jonny Schäfer committed
42
//
43
// The rvp-6 value is used differently depending on the product type:
Jonny Schäfer's avatar
Jonny Schäfer committed
44
// (also consult the DataUnit field of the Composite)
45
46
47
48
49
50
51
52
53
54
//
//	Product label            | rvp-6 value represents   | 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
//	 PR*, ...                | doppler radial velocity  | m/s
//
// The cloud reflectivity (in dBZ) can be converted to rainfall rate (in mm/h)
// via PrecipitationRate().
Jonny Schäfer's avatar
Jonny Schäfer committed
55
56
57
58
59
//
// The cloud reflectivity factor Z is stored in its logarithmic representation dBZ:
//	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:
Jonny Schäfer's avatar
Jonny Schäfer committed
60
//	// if c.HasProjection
Jonny Schäfer's avatar
Jonny Schäfer committed
61
62
//	x, y := c.Translate(52.51861, 13.40833)	// Berlin (lat, lon)
//
Jonny Schäfer's avatar
Jonny Schäfer committed
63
64
65
//	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
Jonny Schäfer's avatar
Jonny Schäfer committed
66
67
68
69
70
71
//
//	fmt.Println("Rainfall in Berlin [mm/h]:", rat)
//
type Composite struct {
	Product string // composite product label

72
73
74
	CaptureTime  time.Time     // time of source data capture used for forcasting
	ForecastTime time.Time     // data represents conditions predicted for this time
	Interval     time.Duration // time duration until next forecast
Jonny Schäfer's avatar
Jonny Schäfer committed
75

76
77
78
	PlainData [][]RVP6 // rvp-6 data for parsed plain data element [y][x]
	Px        int      // plain data width
	Py        int      // plain data height
Jonny Schäfer's avatar
Jonny Schäfer committed
79
	DataUnit Unit
80
81
82

	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])
Jonny Schäfer's avatar
Jonny Schäfer committed
83
84
85

	Dx int // data width
	Dy int // data height
86
	Dz int // data layer
Jonny Schäfer's avatar
Jonny Schäfer committed
87
88
89
90

	Rx float64 // horizontal resolution in km/px
	Ry float64 // vertical resolution in km/px

Jonny Schäfer's avatar
Jonny Schäfer committed
91
92
	HasProjection bool // coordinate translation available

Jonny Schäfer's avatar
Jonny Schäfer committed
93
94
	dataLength int // length of binary section in bytes

95
	precision int    // multiplicator 10^precision for each raw value
96
	level     []RVP6 // maps data value to corresponding rvp-6 value in runlength based formats
Jonny Schäfer's avatar
Jonny Schäfer committed
97
98
99
100
101

	offx float64 // horizontal projection offset
	offy float64 // vertical projection offset
}

Jonny Schäfer's avatar
Jonny Schäfer committed
102
// NewComposite reads binary data from rd and parses the composite.
Jonny Schäfer's avatar
Jonny Schäfer committed
103
104
105
106
107
108
109
110
111
112
113
114
115
func NewComposite(rd io.Reader) (comp *Composite, err error) {
	reader := bufio.NewReader(rd)
	comp = &Composite{}

	err = comp.parseHeader(reader)
	if err != nil {
		return
	}

	err = comp.parseData(reader)
	if err != nil {
		return
	}
116
	comp.arrangeData()
Jonny Schäfer's avatar
Jonny Schäfer committed
117
118
119
120
121
122

	comp.calibrateProjection()

	return
}

Jonny Schäfer's avatar
Jonny Schäfer committed
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
// NewComposites reads tar gz data from rd and returns the parsed composites sorted by
// ForecastTime in ascending order.
func NewComposites(rd io.Reader) ([]*Composite, error) {
	gzipReader, err := gzip.NewReader(rd)
	if err != nil {
		return nil, err
	}
	defer gzipReader.Close()

	tarReader := tar.NewReader(gzipReader)

	var cs []*Composite
	for {
		_, err := tarReader.Next()
		if err == io.EOF {
			break
		}
		if err != nil {
			return nil, err
		}

		c, err := NewComposite(tarReader)
		if err != nil {
			return nil, err
		}
		cs = append(cs, c)
	}

	// sort composites in chronological order
	sort.Slice(cs, func(i, j int) bool { return cs[i].ForecastTime.Before(cs[j].ForecastTime) })
	return cs, nil
}

Jonny Schäfer's avatar
Jonny Schäfer committed
156
157
158
159
160
161
162
163
// NewDummy creates a blank dummy composite with the given product label and dimensions. It can
// be used for generic coordinate translation.
func NewDummy(product string, dx, dy int) (comp *Composite) {
	comp = &Composite{Product: product, Dx: dx, Dy: dy}
	comp.calibrateProjection()
	return
}

Jonny Schäfer's avatar
Jonny Schäfer committed
164
165
166
167
// 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 {
168
169
170
171
172
173
174
175
	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 {
Jonny Schäfer's avatar
Jonny Schäfer committed
176
177
		return RVP6(math.NaN())
	}
178
	return c.DataZ[z][y][x]
Jonny Schäfer's avatar
Jonny Schäfer committed
179
180
181
182
183
184
}

// newError returns an error indicating the failed function and reason
func newError(function, reason string) error {
	return fmt.Errorf("radolan.%s: %s", function, reason)
}