In [ ]:
%matplotlib inline
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import pymepps
The first step is to load the dataset. This will be performed with pymepps.open_model_dataset. The NetCDF4 backend is also supporting opendap paths. So we could specify nc as data type.
In [ ]:
metno_path = 'http://thredds.met.no/thredds/dodsC/meps25files/' \
'meps_det_pp_2_5km_latest.nc'
metno_ds = pymepps.open_model_dataset(metno_path, 'nc')
The resulting dataset is a SpatialDataset. The dataset has several methods to load a xr.DataArray from the path. It also possible to print the content of the dataset. The content contains the dataset type, the number of file handlers within the dataset and all available data variables.
In [ ]:
print(metno_ds)
The next step is to select/extract a variable from the Dataset. We will select the air temperature in 2 metre height and print the content of the resulting data
In [ ]:
metno_t2m = metno_ds.select('air_temperature_2m')
print(metno_t2m)
metno_t2m.isel(validtime=0).plot()
plt.show()
We could see that the resulting data is a normal xarray.DataArray and all of the DataArray methods could be used. The coordinates of the DataArray are normalized. The DataArray is expanded with an accessor. Also the coordinates are normalized. We could access the accessor with metno_t2m.pp. The main methods of the accessor are allowing a grid handling. So our next step is to explore the grid of the DataArray.
In [ ]:
print(metno_t2m.pp.grid)
We could see that the grid is a grid with a defined projection. In our next step we will slice out an area around Hamburg. We will see that a new DataArray with a new grid is created.
In [ ]:
hh_bounds = [9, 54, 11, 53]
t2m_hh = metno_t2m.pp.sellonlatbox(hh_bounds)
print(t2m_hh.pp.grid)
print(t2m_hh)
We sliced a longitude and latitude box around the given grid. So we sliced the data in a longitude and latitude projection. Our original grid was in another projection with unstructured lat lon coordinates. So it is not possible to create a structured grid based on this slice. So the grid becomes an unstructured grid. In the next step we will show the remapping capabilities of the pymepps grid structure.
If we slice the data we have seen that the structured grid could not maintained. So in the next step we will create a structured LonLatGrid from scratch. After the grid building we will remap the raw DataArray basen on the new grid.
The first step is to calculate the model resolution in degree.
In [ ]:
res = 2500 # model resolution in metre
earth_radius = 6371000 # Earth radius in metre
res_deg = np.round(res*360/(earth_radius*2*np.pi), 4)
# rounded model resolution equivalent in degree if it where on the equator
print(res_deg)
Our next step is to build the grid. The grid implementation is inspired by the climate data operators. So to build the grid we will use the same format.
In [ ]:
grid_dict = dict(
gridtype='lonlat',
xsize=int((hh_bounds[2]-hh_bounds[0])/res_deg),
ysize=int((hh_bounds[1]-hh_bounds[3])/res_deg),
xfirst=hh_bounds[0],
xinc=res_deg,
yfirst=hh_bounds[3],
yinc=res_deg,
)
Now we use our grid dict together with the GridBuilder to build our grid.
In [ ]:
builder = pymepps.GridBuilder(grid_dict)
hh_grid = builder.build_grid()
print(hh_grid)
Now we created the grid. The next step is a remapping of the raw DataArray to the new Grid. We will use th enearest neighbour approach to remap the data.
In [ ]:
t2m_hh_remapped = metno_t2m.pp.remapnn(hh_grid)
print(t2m_hh_remapped)
To plot the data in a map, we have to slice the data. We will select the first validtime as plotting parameter.
In [ ]:
t2m_hh_remapped.isel(validtime=0).plot()
plt.show()
In the map around Hamburg we could see the north and baltic sea in the top edges. But with the nearest enighbour approach we retain some of the sharp edges at the map. Our last step is a second remap plot, this time with a bilinear approach.
In [ ]:
# sphinx_gallery_thumbnail_number = 3
metno_t2m.pp.remapbil(hh_grid).isel(validtime=0).plot()
plt.show()