Generate Geotiff from Rainbow file with wradlib

[d-siden]

Hi,

I’m trying to follow this tutorial in order to get data from a Rainbow (.vol) file and export it to a Geotiff file. However the tutorial refers to RADOLAN dataset, and I’m having difficulties in figuring out how to proceed to do the same with my file.
How exactly should this tutorial be adapted for a Rainbow file?

Thanks in advance,
Danilo

[kmuehlbauer]

Hi Danilo,

The example you are referring works on cartesian data. Sou you would have to interpolate your polar data to cartesian before.

There is no combined example (cartesian interpolation and output to GIS) so far. It’s on my todo-list, but not as immediate priority. If you are up to it, I’ll gladly assist.

[d-siden]

UPDATE: First try

In my first attempt I exported the values to a .csv file with their projected coordinates.
This file is loaded into QGIS, then made into a points vector layer, which is converted into raster layers with different resolutions.

Drawbacks:

1. The raster layer has a equally dimensioned area of value validity (like a square), which does not represent the actual shape of coverage of a radar bin.
2. A raster layer has fixed resolution, which again does not represent the variable resolution of a radar sweep data (decreasing with range).
3. Data loss: The value of the pixel comes from a single point inside it, other values are discarded.

To ‘fix’ drawback no. 2, the resolution is lowered until the pixels cover the entire area of the desired region.

This first image show the points vector layer loaded in QGIS.

This one shows the previous layer rasterized to range resolution, the inner points become an uniform layer, but outer points still miss data between them.

This one show the layer rasterized to 10 times the range resolution, the filled area reaches outer as the resolution becomes poorer.

I’ll be trying new methods and keeping the community updated.

[d-siden]

UPDATE: Second try

The second method I’m trying is to save an black and white PPI plot without any extra borders or image elements to PNG format, alongside a World File (.pngw), which QGIS can load to georreference the image.

The positive effect is that the plotted image has a good representation of the area that a bin represents.

Unfortunately the PNG seems to be somewhat misaligned when seen under the point layer in QGIS.

The points come from wradlib’s spherical_to_centroids, as described in the First Try Update above.

Red and blue color scale shows the points.

Zoomed in image shows raster (PNG, black and white) under points. Notice that both layers show the same variations on their respective colors.

But looking closer at the bins, the point (centroid) does not lay on the center of the area. Radar is located to the lower right side. Points are always closer to the radar than a pixelated area in any direction.

[d-siden]

This misalignment was fixed by subtracting half of the range from the range ndarray when plotting.
The next step is saving the values in png, because python png saving functions only support saving 8 bit integer values, whose range is between 0 and 255. So, rounding the number and adding the minimum value can temporarily solve the issue (if the data is reflectivity, whose range is between -31.5 and ~100 dBZ). Just remember to subtract again when reading.

[d-siden]

Update: Still haven’t found a way to achieve this.

I’m currently using Wradlib 2.1.0 and I still haven’t managed to make it work.
I’m trying this function:

wradlib.io.gdal.write_raster_dataset("geotiff.tif", dados, driver="GTiff")

Where “dados” comes from a correct implementation of:

wradlib.georef.create_xarray_dataarray(...)

Followed by

dados = wradlib.georef.georeference(dados, crs=epsg_to_osr(4326))

An attempt of writing the dataset gives

TypeError: in method ‘Driver_CreateCopy’, argument 3 of type ‘GDALDatasetShadow *’

I still can’t find a solution for this, Is there an actual way to solve it?

[kmuehlbauer]

Hi @d-siden,

Danilo, would it be possible for you to create a MCVE (How to create a Minimal, Reproducible Example - Help Center - Stack Overflow) showing the problem? It would help tremendously in getting to the bottom of the issue here, if people can start out from working with the code.

Thanks,
Kai

[d-siden]

Of course! Here is the code:

import wradlib # 2.1.0
from os import environ
from numpy import arange # numpy 2.0.1

environ["WRADLIB_DATA"] = "/home/danilo/Documentos/Python/radar"
sample_file = environ["WRADLIB_DATA"]+"/outros-dados/2025030100000400dBZ.vol"
output_file = environ["WRADLIB_DATA"]+"/geotiff.tif"

content = wradlib.io.read_rainbow(sample_file)
azi = content["volume"]["scan"]["slice"][0]["slicedata"]["rayinfo"][0]["data"]
azi_bits = float(content["volume"]["scan"]["slice"][0]["slicedata"]["rayinfo"][0]["@depth"])
n_rays = float(content["volume"]["scan"]["slice"][0]["slicedata"]["rayinfo"][0]["@rays"])
angstep = float(content["volume"]["scan"]["slice"][0]["anglestep"])
azi = (azi * n_rays / 2 ** azi_bits) * angstep

maxrange = float(content["volume"]["scan"]["slice"][0]["stoprange"])
rangstep = float(content["volume"]["scan"]["slice"][0]["rangestep"])
# convert range from kilometers to meters
radar_range = arange(rangstep * 1000, maxrange * 1000 + rangstep * 1000, rangstep * 1000)

data = content["volume"]["scan"]["slice"][0]["slicedata"]["rawdata"]["data"]
data_bits = float(content["volume"]["scan"]["slice"][0]["slicedata"]["rawdata"]["@depth"])
mindata = float(content["volume"]["scan"]["slice"][0]["slicedata"]["rawdata"]["@min"])
maxdata = float(content["volume"]["scan"]["slice"][0]["slicedata"]["rawdata"]["@max"])
data = mindata + data * (maxdata - mindata) / 2 ** data_bits

elevation = float(content["volume"]["scan"]["slice"][0]["posangle"])

lat_radar = float(content["volume"]["sensorinfo"]["lat"])
lon_radar = float(content["volume"]["sensorinfo"]["lon"])
alt_radar = float(content["volume"]["sensorinfo"]["alt"])

data = wradlib.georef.create_xarray_dataarray(data, r=radar_range, phi=azi, theta=elevation,
                                              site=(lon_radar, lat_radar, alt_radar))
data = wradlib.georef.georeference(data, crs=wradlib.georef.epsg_to_osr(4326))

wradlib.io.gdal.write_raster_dataset(output_file, data, driver="GTiff")

And the sample file.

Thanks in advance,
Danilo

[kmuehlbauer]

@d-siden Try this as starting point

import numpy as np
import wradlib as wrl
import xarray as xr

sample_file = "/home/kai/data/daten/daten/2025030100000400dBZ.vol"
output_file = "geotiff.tif"

# use xarray plus xradar engine
ds = xr.open_dataset(sample_file, engine="rainbow", group="sweep_0")

# georeference (might use some more fitting projection!)
crs = wrl.georef.epsg_to_osr(4326)
ds = ds.set_coords("sweep_mode").wrl.georef.georeference(crs=crs)

# create grid (better use a fitting projection above) this is just to showcase
xgrid = np.linspace(ds.x.min().values, ds.x.max().values, 1000)
ygrid = np.linspace(ds.y.min().values, ds.y.max().values, 1000)
cart = xr.Dataset(coords={"x": (["x"], xgrid), "y": (["y"], ygrid)})

# use nearest neighbour gridder
gridded = ds.DBZH.wrl.comp.togrid(cart, radius=200000.0, center=(0,0), interpol=wrl.ipol.Nearest)

# output to raster file, needs rasterio/rioxarray packages
gridded.rio.to_raster(output_file)

[d-siden]

That works!

I’ll be studying this method, I thought that the rainbow tutorial steps were the only way to get those files read.

Thank you so much!

[kmuehlbauer]

Hi @d-siden,

great it works for you. And yes, we could improve the documentation over at wradlib. All polar radar data related code which existed in wradlib has been ported over to xradar — xradar 0.9.0 documentation. There is a section on rainbow, too.

Would you mind opening an issue over at GitHub - wradlib/wradlib: weather radar data processing - python package with some details how the particular part of the documentation could be improved?

Thanks,
Kai