Destaggering U and V wind components

Introduction and aims

  • Atmospheric models are often run on a “staggered” grid for computational stability. That is, some model variables are offset by half a grid length relative to other variables.

  • This includes the ACCESS model and AUS2200. In the case of AUS2200, U and V wind components have been saved on their staggered grids, and therefore need to be de-staggered prior to analysis.

  • This notebook demonstrates two simple methods for destaggering U and V wind components from AUS2200 data, and visualises the output with a quiver plot.

[1]:

#Assumes access to xp65 conda environment

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
from dask.distributed import Client

[2]:

#Start a distributed dask client. You can click "Launch dashboard in JupyterLab" in the cell output to see progress

client = Client(threads_per_worker=1)
client

[2]:

Client

Client-c7d89de2-4f66-11f1-9282-0000018dfe80

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: /proxy/8787/status

Cluster Info

 
[4]:

# Lets use intake to grab a couple of AUS2200 files to play with
import intake


datastore = intake.cat.access_nri["AUS2200"]
datastore = datastore.search(realm='atmos')
datastore = datastore.search(experiment_id='ecoastlow-evolvsst')
datastore = datastore.search(frequency='1hr')
datastore = datastore.search(variable_id=["uas","vas"])

uas = datastore.search(variable_id="uas").to_dask(xarray_open_kwargs={'chunks' : {}})
vas = datastore.search(variable_id="vas").to_dask(xarray_open_kwargs={'chunks' : {}})

[5]:

uas = uas.sel(time='2016-06-05 18:00',method="nearest").squeeze()
vas = vas.sel(time='2016-06-05 18:00',method="nearest").squeeze()

[6]:

print("Experiment description: \n")
print(uas.attrs["exp_description"])

Experiment description:

A limited area model study of the entire Australian continent at 2.2 km resolution, using the UM atmospheric model. ERA5+ERA5Land reanalysis data was used to provide initial and boundary conditions, sea-surface temperatures are evolving throughout the simulation. The study covers the period from 2016-06-03 to 2016-06-07, including the East Coast Low event which gave rise to widespread flooding in many areas stretching from southeast Queensland, eastern New South Wales, eastern Victoria, and large areas of northern Tasmania.

You can see more information on this runs, including how to cite it,here

Chunking

We have loaded in in 10-metre u and v wind data, with 4 chunks in the lat and lon directions. We will return to this later when we consider different destaggering methods

[7]:

uas.uas

[7]:

<xarray.DataArray 'uas' (lat: 2120, lon: 2600)> Size: 22MB
dask.array<getitem, shape=(2120, 2600), dtype=float32, chunksize=(1060, 1300), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 17kB -48.79 -48.77 -48.75 ... -6.871 -6.852 -6.832
  * lon      (lon) float64 21kB 114.3 114.3 114.3 114.3 ... 165.7 165.7 165.7
    time     datetime64[ns] 8B 2016-06-05T18:30:00.000000256
    height   float64 8B ...
Attributes:
    standard_name:          eastward_wind
    long_name:              Eastward Near-Surface Wind
    comment:                Eastward component of the near-surface (usually, ...
    units:                  m s-1
    cell_methods:           area: mean time: mean
    cell_measures:          area: areacella
    history:                2024-10-16T01:32:52Z altered by CMOR: Treated sca...
    coverage_content_type:  modelResult

Staggered grid

The AUS2200 u wind data has been saved on a staggered longitude grid, while the v wind data has been saved on a staggered latitude grid. See how the longitude coordinates are different (and offset by 0.5 * grid spacing). If we were to combine u and v now into one dataset using something like xr.Dataset({"u":uas.uas,"v":vas:vas}), then xarray would do it by meshing together the grids with missing values (this is not what we want to do).

[8]:

vas.lon.isel(lon=slice(0,20)).plot(marker="o")
uas.lon.isel(lon=slice(0,20)).plot(marker="o")
plt.gca().grid(ls=":")
../_images/Recipes_AUS2200_destagger_winds_10_0.png

The staggered grid has also resulted in one extra latitude grid point, so the vas variable is on a different size grid

[9]:

uas.uas.shape

[9]:

(2120, 2600)

[10]:

vas.vas.shape

[10]:

(2121, 2600)

Destaggering (centring)

If we want to perform any analysis or make plots, we need u and v to be on the same grid. Our first method will be to centre u in longitude and v in latitude by averaging

[11]:

#Set up the destaggered coordinates.
#Note that our centring approach will result in one less grid point, so we need to slice the lons down accordingly
#The latitude for v is already one grid point longer than u, so we don't need to do this for latitude

destaggered_lon_coord = vas.lon.isel({"lon":slice(0,-1)})
destaggered_lat_coord = uas.lat

[12]:

#Now we can slice and centre/average the U data in longitude.
#We need to manually reassign the longitude coordinates using assign_coords()
uas_destaggered = ((uas.uas.isel({"lon":slice(0,-1)}).assign_coords({"lon":destaggered_lon_coord}) +\
      uas.uas.isel({"lon":slice(1,uas.lon.shape[0])}).assign_coords({"lon":destaggered_lon_coord})) / 2)

[13]:

#Same for V
vas_destaggered = ((vas.vas.isel({"lat":slice(0,-1)}).assign_coords({"lat":destaggered_lat_coord}) +\
      vas.vas.isel({"lat":slice(1,vas.lat.shape[0])}).assign_coords({"lat":destaggered_lat_coord})) / 2)

Combining

Now we can safely combine the u and v wind components into one dataset (useful for plotting later), and calculate wind speed.

Note that out chunking has been preserved in this method

[14]:

destaggered_ds = xr.Dataset({"u":uas_destaggered,"v":vas_destaggered})
destaggered_ds["wind_speed"] = np.sqrt( destaggered_ds["u"]**2 + destaggered_ds["v"]**2 )

[15]:

destaggered_ds["u"]

[15]:

<xarray.DataArray 'u' (lat: 2120, lon: 2600)> Size: 22MB
dask.array<where, shape=(2120, 2600), dtype=float32, chunksize=(1060, 866), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 17kB -48.79 -48.77 -48.75 ... -6.871 -6.852 -6.832
  * lon      (lon) float64 21kB 114.3 114.3 114.3 114.3 ... 165.7 165.7 165.7
    time     datetime64[ns] 8B 2016-06-05T18:30:00.000000256
    height   float64 8B 10.0
Attributes:
    standard_name:          eastward_wind
    long_name:              Eastward Near-Surface Wind
    comment:                Eastward component of the near-surface (usually, ...
    units:                  m s-1
    cell_methods:           area: mean time: mean
    cell_measures:          area: areacella
    history:                2024-10-16T01:32:52Z altered by CMOR: Treated sca...
    coverage_content_type:  modelResult

Destaggering (interp)

It would be a little neater to just use xr.interp to destagger, which is demonstrated below. But note that our chunking is not preserved in this method, because dask needs to access the whole lat/lon dimension to do the spatial interpolation. So we now have only one chunk in longitude for the U variable. This is an important consideration if you need to destagger on really large data (not an issue here for one time step, and surface level data only)

[16]:

uas_destaggered_interp = uas.uas.interp({"lon":destaggered_lon_coord},method="linear")
vas_destaggered_interp = vas.vas.interp({"lat":destaggered_lat_coord},method="linear")
destaggered_ds_interp = xr.Dataset({"u":uas_destaggered_interp,"v":vas_destaggered_interp})
destaggered_ds_interp["wind_speed"] = np.sqrt( destaggered_ds_interp["u"]**2 + destaggered_ds_interp["v"]**2 )
destaggered_ds_interp["u"]

[16]:

<xarray.DataArray 'u' (lat: 2120, lon: 2600)> Size: 22MB
dask.array<where, shape=(2120, 2600), dtype=float32, chunksize=(530, 2599), chunktype=numpy.ndarray>
Coordinates:
  * lat      (lat) float64 17kB -48.79 -48.77 -48.75 ... -6.871 -6.852 -6.832
  * lon      (lon) float64 21kB 114.3 114.3 114.3 114.3 ... 165.7 165.7 165.7
    time     datetime64[ns] 8B 2016-06-05T18:30:00.000000256
    height   float64 8B 10.0
Attributes:
    standard_name:          eastward_wind
    long_name:              Eastward Near-Surface Wind
    comment:                Eastward component of the near-surface (usually, ...
    units:                  m s-1
    cell_methods:           area: mean time: mean
    cell_measures:          area: areacella
    history:                2024-10-16T01:32:52Z altered by CMOR: Treated sca...
    coverage_content_type:  modelResult

Plotting

Visualise the winds using a quiver plot and contour the wind speed. Both methods of destaggering produce the same outcome

[17]:

#Set the spacing of the wind vectors in lat and lon
step = 75

plt.figure(figsize=[16,6])
ax = plt.subplot(1,2,1,projection=ccrs.PlateCarree())
destaggered_ds["wind_speed"].plot(robust=True,cmap="Blues")
destaggered_ds.isel({"lat":slice(0,-1,step),"lon":slice(0,-1,step)}).plot.quiver(
    "lon","lat","u","v")
ax.coastlines()
plt.title("Centring")

ax = plt.subplot(1,2,2,projection=ccrs.PlateCarree())
destaggered_ds_interp["wind_speed"].plot(robust=True,cmap="Blues")
destaggered_ds_interp.isel({"lat":slice(0,-1,step),"lon":slice(0,-1,step)}).plot.quiver(
    "lon","lat","u","v")
ax.coastlines()
plt.title("xr.interp")

[17]:

Text(0.5, 1.0, 'xr.interp')
../_images/Recipes_AUS2200_destagger_winds_24_1.png

Contributors

  • Andrew Brown

    ARC Centre of Excellence for 21st Century Weather, University of Melbourne

  • Charles Turner

    ACCESS-NRI, Australian National University, Canberra

If you have any enquries, suggested improvements or bug reports related to this recipe, please open an issue or start a discussion in this repository.