Export a dataset in GIS-compatible format

Contents

Export a dataset in GIS-compatible format#

In this notebook, we demonstrate how to export a gridded dataset in GeoTIFF and ESRI ASCII format. This will be exemplified using RADOLAN data from the German Weather Service.

You have two options for output:

rioxarray.to_raster
builtin GDAL functionality

import matplotlib.pyplot as plt
import os
import wradlib as wrl
import wradlib_data
import xarray as xr
import numpy as np
import warnings
from pyproj.crs import CRS
from IPython.display import display

warnings.filterwarnings("ignore")

Step 2a (output with rioxarray)#

drop encoding

ds.SF.encoding = {}

ds = ds.rio.write_crs(crs)
ds.SF.rio.to_raster(wdir / "geotiff_rio.tif", driver="GTiff")

ds.SF.rio.to_raster(
    wdir / "aaigrid_rio.asc",
    driver="AAIGrid",
    profile_kwargs=dict(options=["DECIMAL_PRECISION=2"]),
)

Step 2b: (output with GDAL)#

Get the projected coordinates of the RADOLAN grid#

# Get projected RADOLAN coordinates for corner definition
xy_raw = wrl.georef.get_radolan_grid(900, 900)
xy_raw.shape

(900, 900, 2)

Check Origin and Row/Column Order#

We know, that wradlib.io.read_radolan_composite() returns a 2D-array (rows, cols) with the origin in the lower left corner. Same applies to wradlib.georef.get_radolan_grid(). For the next step, we need to flip the data and the coords up-down. The coordinate corner points also need to be adjusted from lower left corner to upper right corner.

data, xy = wrl.georef.set_raster_origin(ds.SF.values, xy_raw, "upper")
print(data.shape)

(900, 900)

Export as GeoTIFF#

For RADOLAN grids, this projection will probably not be recognized by ESRI ArcGIS.

# create 3 bands
data = np.stack((data, data + 100, data + 1000), axis=0)
print(data.shape)
gds = wrl.georef.create_raster_dataset(data, xy, crs=proj_osr)
wrl.io.write_raster_dataset(wdir / "geotiff.tif", gds, driver="GTiff")

(3, 900, 900)

Export as ESRI ASCII file (aka Arc/Info ASCII Grid)#

# Export to Arc/Info ASCII Grid format (aka ESRI grid)
#     It should be possible to import this to most conventional
# GIS software.
# only use first band
proj_esri = proj_osr.Clone()
proj_esri.MorphToESRI()
ds = wrl.georef.create_raster_dataset(data[0], xy, crs=proj_esri)
wrl.io.write_raster_dataset(
    wdir / "aaigrid.asc", ds, driver="AAIGrid", options=["DECIMAL_PRECISION=2"]
)

Step 3b: Read with GDAL#

fig = plt.figure(figsize=(15, 6))
ax1 = fig.add_subplot(121)
ds1 = wrl.io.open_raster(wdir / "geotiff.tif")
data1, xy1, proj1 = wrl.georef.extract_raster_dataset(ds1, nodata=-9999.0)
ax1.pcolormesh(xy1[..., 0], xy1[..., 1], data1[0])

ax2 = fig.add_subplot(122)
ds2 = wrl.io.open_raster(wdir / "geotiff_rio.tif")
data2, xy2, proj2 = wrl.georef.extract_raster_dataset(ds2, nodata=-9999.0)
ax2.pcolormesh(xy2[..., 0], xy2[..., 1], data2)

<matplotlib.collections.QuadMesh at 0x7991f4d3e0d0>

../../../_images/68a0d935f28f9d91184ebea8561d3c210c9cd0d59aca752072ff7b9ccd57d01d.png

fig = plt.figure(figsize=(15, 6))
ax1 = fig.add_subplot(121)
ds1 = wrl.io.open_raster(wdir / "aaigrid.asc")
data1, xy1, proj1 = wrl.georef.extract_raster_dataset(ds1, nodata=-9999.0)
ax1.pcolormesh(xy1[..., 0], xy1[..., 1], data1)

ax2 = fig.add_subplot(122)
ds2 = wrl.io.open_raster(wdir / "aaigrid_rio.asc")
data2, xy2, proj2 = wrl.georef.extract_raster_dataset(ds2, nodata=-9999.0)
ax2.pcolormesh(xy2[..., 0], xy2[..., 1], data2)

<matplotlib.collections.QuadMesh at 0x7992025439d0>

../../../_images/68a0d935f28f9d91184ebea8561d3c210c9cd0d59aca752072ff7b9ccd57d01d.png