radolan.go 7.74 KB
Newer Older
Jonny Schäfer's avatar
Jonny Schäfer committed
1
// Package radolan parses the DWD RADOLAN / RADVOR radar composite format. This data
Jonny Schäfer's avatar
Jonny Schäfer committed
2
// is available at the Open Data Portal (https://www.dwd.de/DE/leistungen/opendata/opendata.html).
Jonny Schäfer's avatar
Jonny Schäfer committed
3
4
// The obtained results can be processed and visualized with additional functions.
//
Jonny Schäfer's avatar
Jonny Schäfer committed
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// Tested input products and grids:
//
//	| Product | Grid              | Description             |
//	| ------- | ----------------- | ----------------------- |
//	| EX      | middle-european   | reflectivity            |
//	| FX      | national          | nowcast reflectivity    |
//	| FZ      | national          | nowcast reflectivity    |
//	| PE      | local             | echo top                |
//	| PF      | local             | reflectivity            |
//	| PG      | national picture  | reflectivity            |
//	| PR      | local             | doppler radial velocity |
//	| PX      | local             | reflectivity            |
//	| PZ      | local             | 3D reflectivity CAPPI   |
//	| RW      | national          | hourly accumulated      |
//	| RX      | national          | reflectivity            |
//	| SF      | national          | daily accumulated       |
//	| WX      | extended national | reflectivity            |
//
// Those can be considered working with sufficient accuracy.
Jonny Schäfer's avatar
Jonny Schäfer committed
24
25
26
27
28
29
30
//
// 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
Jonny Schäfer's avatar
Jonny Schäfer committed
31
32
33
34
//	[3] https://www.dwd.de/DE/leistungen/radarprodukte/formatbeschreibung_fxdaten.pdf
//	[4] https://www.dwd.de/DE/leistungen/opendata/help/radar/radar_pg_coordinates_pdf.pdf
//	[5] https://www.dwd.de/DE/leistungen/radarniederschlag/rn_info/download_niederschlagsbestimmung.pdf
//	[6] hex editor and much reverse engineering
Jonny Schäfer's avatar
Jonny Schäfer committed
35
36
37
package radolan

import (
Jonny Schäfer's avatar
Jonny Schäfer committed
38
	"archive/tar"
Jonny Schäfer's avatar
Jonny Schäfer committed
39
	"bufio"
40
	"compress/bzip2"
Jonny Schäfer's avatar
Jonny Schäfer committed
41
42
	"fmt"
	"io"
Jonny Schäfer's avatar
Jonny Schäfer committed
43
	"sort"
Jonny Schäfer's avatar
Jonny Schäfer committed
44
45
46
	"time"
)

47
48
49
50
51
52
53
// 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
54
55
56
// 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 ].
Jonny Schäfer's avatar
Jonny Schäfer committed
57
//
58
// The data value is used differently depending on the product type:
Jonny Schäfer's avatar
Jonny Schäfer committed
59
// (also consult the DataUnit field of the Composite)
60
//
61
//	Product label            | values represent         | unit
62
63
//	-------------------------+--------------------------+------------------------
//	 PG, PC, PX*, ...        | cloud reflectivity       | dBZ
64
65
//	 RX, WX, EX, FZ, FX, ... | cloud reflectivity	    | dBZ
//	 RW, SF,  ...            | aggregated precipitation | mm/interval
66
67
68
69
//	 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
70
71
72
73
74
//
// 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
75
//	// if c.HasProjection
Jonny Schäfer's avatar
Jonny Schäfer committed
76
77
//	x, y := c.Translate(52.51861, 13.40833)	// Berlin (lat, lon)
//
78
79
//	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
Jonny Schäfer's avatar
Jonny Schäfer committed
80
81
82
83
84
85
//
//	fmt.Println("Rainfall in Berlin [mm/h]:", rat)
//
type Composite struct {
	Product string // composite product label

86
87
88
	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
89

Jonny Schäfer's avatar
Jonny Schäfer committed
90
	DataUnit Unit
91

92
93
94
95
96
	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)
Jonny Schäfer's avatar
Jonny Schäfer committed
97
	Data  [][]float32   // data for each pixel [y][x] at layer 0 (alias for DataZ[0][x][y])
Jonny Schäfer's avatar
Jonny Schäfer committed
98
99
100

	Dx int // data width
	Dy int // data height
101
	Dz int // data layer
Jonny Schäfer's avatar
Jonny Schäfer committed
102
103
104
105

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

Jonny Schäfer's avatar
Jonny Schäfer committed
106
107
	HasProjection bool // coordinate translation available

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

110
111
	precision int       // multiplicator 10^precision for each raw value
	level     []float32 // maps data value to corresponding index value in runlength based formats
Jonny Schäfer's avatar
Jonny Schäfer committed
112
113
114
115
116

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

117
118
119
120
121
122
123
124
125
126
// 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.
Jonny Schäfer's avatar
Jonny Schäfer committed
127
128
129
130
131
132
133
134
135
136
137
138
139
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
	}
140
	comp.arrangeData()
Jonny Schäfer's avatar
Jonny Schäfer committed
141
142
143

	comp.calibrateProjection()

144
145
146
147
	if comp.DataUnit == Unit_unknown {
		err = ErrUnknownUnit
	}

Jonny Schäfer's avatar
Jonny Schäfer committed
148
149
150
	return
}

151
// NewComposites reads .tar.bz2 data from rd and returns the parsed composites sorted by
Jonny Schäfer's avatar
Jonny Schäfer committed
152
153
// ForecastTime in ascending order.
func NewComposites(rd io.Reader) ([]*Composite, error) {
154
	bzipReader := bzip2.NewReader(rd)
Jonny Schäfer's avatar
Jonny Schäfer committed
155

156
	tarReader := tar.NewReader(bzipReader)
Jonny Schäfer's avatar
Jonny Schäfer committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179

	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
180
181
182
183
184
185
186
187
// 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
}

188
189
190
191
// 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 {
192
193
194
	return c.AtZ(x, y, 0)
}

195
196
197
198
// 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 {
199
	if x < 0 || y < 0 || z < 0 || x >= c.Dx || y >= c.Dy || z >= c.Dz {
200
		return NaN
Jonny Schäfer's avatar
Jonny Schäfer committed
201
	}
202
	return c.DataZ[z][y][x]
Jonny Schäfer's avatar
Jonny Schäfer committed
203
204
205
206
207
208
}

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