# Grid structure of radar composite

**URL:** <https://openradar.discourse.group/t/grid-structure-of-radar-composite/464>\
**Category:** General\
**Created:** [April 30, 2024, 12:34pm UTC](https://openradar.discourse.group/t/grid-structure-of-radar-composite/464 "2024-04-30T12:34:56Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![lou](https://avatars.discourse-cdn.com/v4/letter/l/ea666f/32.png) [@lou](https://openradar.discourse.group/u/lou)\
**Post date:** [April 30, 2024, 12:34pm UTC](https://openradar.discourse.group/t/grid-structure-of-radar-composite/464/1 "2024-04-30T12:34:56Z")

</div>

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:

```auto
# 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](https://global.discourse-cdn.com/free1/uploads/openradar/original/1X/a7a2a210dc12fdb0201c96699b273e0533c4cf3c.png)

---

<div class="post-metadata">

**Author:** ![jorahu](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/jorahu/32/19_2.png) [@jorahu](https://openradar.discourse.group/u/jorahu)\
**Post date:** [May 1, 2024, 7:44am UTC](https://openradar.discourse.group/t/grid-structure-of-radar-composite/464/2 "2024-05-01T07:44:50Z")

</div>

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](https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html).

From the example:

```auto
# 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.

---

<div class="post-metadata">

**Author:** ![lou](https://avatars.discourse-cdn.com/v4/letter/l/ea666f/32.png) [@lou](https://openradar.discourse.group/u/lou)\
**Post date:** [May 1, 2024, 8:16am UTC](https://openradar.discourse.group/t/grid-structure-of-radar-composite/464/3 "2024-05-01T08:16:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jorahu](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/jorahu/32/19_2.png) [@jorahu](https://openradar.discourse.group/u/jorahu)\
**Post date:** [May 1, 2024, 8:40am UTC](https://openradar.discourse.group/t/grid-structure-of-radar-composite/464/4 "2024-05-01T08:40:18Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lou](https://avatars.discourse-cdn.com/v4/letter/l/ea666f/32.png) [@lou](https://openradar.discourse.group/u/lou)\
**Post date:** [May 6, 2024, 9:23am UTC](https://openradar.discourse.group/t/grid-structure-of-radar-composite/464/5 "2024-05-06T09:23:53Z")

</div>

Hi @jorahu,

thanks again. So I used now the following code:

```auto
# 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

---

<div class="post-metadata">

**Author:** ![jorahu](https://yyz2.discourse-cdn.com/free1/user_avatar/openradar.discourse.group/jorahu/32/19_2.png) [@jorahu](https://openradar.discourse.group/u/jorahu)\
**Post date:** [May 6, 2024, 1:39pm UTC](https://openradar.discourse.group/t/grid-structure-of-radar-composite/464/6 "2024-05-06T13:39:21Z")

</div>

Hey

Good to hear that the suggestions worked.

About the reference point of the coordinates, I think it’s the ‘shading’ keyword in [pcolormesh](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.pcolormesh.html) documentation that controls the origin. Check the [shading explanation](https://matplotlib.org/stable/gallery/images_contours_and_fields/pcolormesh_grids.html) 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](https://docs.wradlib.org/en/latest/generated/wradlib.georef.projection.reproject.html)  
[Here](https://docs.wradlib.org/en/latest/notebooks/georeferencing/georef.html) 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.