Grid structure of radar composite

General
1

Hi all,

as I am trying to compare my ICON simulations with radar observation, I want to keep the underlying grid structure in mind. Therefore, I am trying to plot the grid structure of each.

I am reading the data with rasterio, but I am not sure how I could plot the underlying radar grid structure. Could someone help me how I could do this?

The code for plotting my radar data is the following:

# Erstellen der Figuren und Achsen
fig = plt.figure(figsize=(20, 6))

ax = fig.add_subplot(1, 1, 1, projection=ccrs.Mercator())
ax.set_extent([10.3, 38.3, 55.83, 67.83], crs=ccrs.PlateCarree())
ax.coastlines(linewidth=.5)
ax.add_feature(cf.LAND, color='white', alpha=.15)
ax.add_feature(cf.BORDERS, linestyle='-', linewidth=.6)
gl = ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
gl.left_labels = False

# Load data
dataset = rasterio.open(flist_rasterio[3])

# read the data
data = dataset.read()

# mask all NoData-values with True, convert the data with the meta-data values and replace all as true-marked nodata-values with a nan
nodata_mask = data == nodata
data = gain * data + offset
data[nodata_mask] = np.nan

# Reshaped the data
data_reshaped = np.reshape(data, (dataset.height, dataset.width))

# Plotting GeoTIFF-Data
extent = (dataset.bounds.left, dataset.bounds.right, dataset.bounds.bottom, dataset.bounds.top)
boundaries = [0, 0.01, 0.1, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 2, 2.5, 3, 5, 10, 20]
colors = ('#FF000000', '#c4c4e4','#8989c9','#4e4eaf','#141494','#272776','#626256','#9d9d35','#d8d815','#f6eb00','#dbb100','#c17600','#a63b00','#8b0000') 
cmap = mcolors.ListedColormap(colors)
norm = mcolors.BoundaryNorm(boundaries, cmap.N)
im = ax.imshow(data_reshaped, extent=extent, cmap=cmap, norm=norm, transform=get_crs_of_data(dataset))
    
# Add colorbar to figure
divider = make_axes_locatable(ax)
cax = divider.new_horizontal(size="2%", pad=0.05, axes_class=plt.Axes)
fig.add_axes(cax)
cbar = fig.colorbar(im, cax=cax, orientation='vertical', ticks=boundaries, label='Precipitation rate [kg/(m$^2$h)]')

#Adjust the colour bar limits to remove the masked values
cbar.ax.set_ylim(0.01, cbar.ax.get_ylim()[1])

fig.tight_layout()

# Zeige das gesamte Bild an
plt.show()

while my output looks like this:

rasterio_data
rasterio_data704×590 67.5 KB

2

Hey

What do you mean by the “underlying radar grid structure” (sorry for not understanding the terminology)?

If you mean the lat-lon grid, then you can use the bounds from your dataset (dataset.bounds.left and so on) and form meshgrid with numpy.meshgrid.

From the example:

# You could also use np.arange()
lons = np.linspace(dataset.bounds.left, dataset.bounds.right, *insert your unit*)
lats = np.linspace(dataset.bounds.bottom, dataset.bounds.top, *insert your unit*)
lons_grid, lats_grid = np.meshgrid(lons, lats)
# Maybe you want to flip latitudes from top to bottom with:
# lats_grid = np.flipud(lats_grid)

Again, sorry if I misunderstood your question.

3

Hey @jorahu,

yes, I mean the lat/lon-grid. I was just not sure if the data is plottet on a lat/lon grid. Is that with radar data always the case?

Best,
Julian

4

Good. Sorry for being so insecure with my answers. I am not a high-end expert in the field. But sure, having the radar data pixels referenced in the correct geographical location is essential for this kind of data.

Using cartopy allows for different ways to visualise data, and using imshow with extent keyword is definitely one of them for simple visualisation.

If you want to use different data sources and utilise the lat-lon grid, then the meshgrid is solid option to get the grid. But then, instead of imshow, use pcolormesh (or contourf, whatever your needs are).

Instead of the line:
im = ax.imshow(data_reshaped, extent=extent, cmap=cmap, norm=norm, transform=get_crs_of_data(dataset))

use this:
im = ax.pcolormesh(lons_grid, lats_grid, data_reshaped, cmap=cmap, norm=norm, transform=get_crs_of_data(dataset))
The first 2 arguments are the coordinate grids (could also be the 1D lons/lats vectors).

This should give exactly the same output as in your original post, but now you have (and use) the coordinate grids.

5

Hi @jorahu,

thanks again. So I used now the following code:

# define lon, lat for plotting the grid structure
lons = np.linspace(dataset.bounds.left, dataset.bounds.right, dataset.width)
lats = np.linspace(dataset.bounds.top, dataset.bounds.bottom, dataset.height)
lons_grid, lats_grid = np.meshgrid(lons, lats)

# plot the data
im = ax.pcolormesh(lons_grid, lats_grid, data_reshaped, cmap=cmap, norm=norm, ec='black', linewidth=.005, transform=get_crs_of_data(dataset))

where I replaced your *insert your unit* with the width and height of my dataset. I can now see the grid structure and that’s perfect. But, I have another question here: Are the (lon, lat)-tupels the center point of the rectangles?

Furthermore, I am wondering, how I could draw the grid manually as I just need only a smaller domain. Do you knwo how I could translate my lons, lats from my source src (dataset.crs) to the typical lon,lat WGS84 crs?

Best,
Julian

6

Hey

Good to hear that the suggestions worked.

About the reference point of the coordinates, I think it’s the ‘shading’ keyword in pcolormesh documentation that controls the origin. Check the shading explanation for more details.
If you want the reference to be in the centre, then you should explicitly use shading='nearest' in the pcolormesh arguments. Also, the shapes of input matrices need to be the same.

To zoom into some region, I think the safest and easiest option would be to just change the visible extent of the image (ax.set_extent()). This ensures the correct georeference values straight from the files.

Or go dirty and mask the values using (for example) numpy.where(lons > *some value*, np.nan, lons), or directly using numpy array lons[lons>*some value*] = np.nan. You would need to do this for every “side” of the extent and also apply the same filtering to the data matrix (using numpy.isnan()). Maybe someone else has better suggestions here instead.

For reprojecting coordinates, wradlib has a nice function just for that: wradlib.georef.projection.reproject — wradlib
Here is tutorial how to do it also for the xarray structures. But as you use rasterio import, then you probably need to use the function directly.