1 2 3 4 | import xarray as xr import numpy as np import matplotlib.pyplot as plt import xesmf as xe |
1 2 | ### dr4292 import iris |
"targetgrid" is the N48 grid
"era_clima_hour" is the climatology calculated from ERA-I files on 241x480 grid, used here to downscale to N48 grid
1 | targetgrid = xr.open_dataset('/g/data/w40/pf4000/UM/um_experiments/n48/era_n48_tuv.nc') |
1 | era_clima_hour = xr.open_dataset('/g/data/w40/pf4000/UM/um_experiments/clima_erai/ta_ERAI_climate.nc') |
1 | hhrT = era_clima_hour.ta.sel(lev=60, method='nearest') #same as above for highres data |
1 2 3 4 5 6 7 8 9 | ### dr4292 - calculate area weights for hi-res grid hc=iris.load_cube('/g/data/w40/pf4000/UM/um_experiments/clima_erai/ta_ERAI_climate.nc',"ta") hr_lat=hc.dim_coords[2] hr_lat.guess_bounds() hr_lon=hc.dim_coords[3] hr_lon.guess_bounds() hr_weights=iris.analysis.cartography.area_weights(hc[0][0],normalize=True) hr_weights=xr.DataArray(hr_weights).rename(dim_0="lat",dim_1="lon") hr_weights |
/g/data/hh5/public/apps/cms_conda/envs/analysis3-23.01/lib/python3.9/site-packages/iris/fileformats/cf.py:584: UserWarning: Missing CF-netCDF formula term variable 'hyam', referenced by netCDF variable 'lev' warnings.warn( /g/data/hh5/public/apps/cms_conda/envs/analysis3-23.01/lib/python3.9/site-packages/iris/fileformats/cf.py:584: UserWarning: Missing CF-netCDF formula term variable 'hybm', referenced by netCDF variable 'lev' warnings.warn( /g/data/hh5/public/apps/cms_conda/envs/analysis3-23.01/lib/python3.9/site-packages/iris/fileformats/cf.py:584: UserWarning: Missing CF-netCDF formula term variable 'aps', referenced by netCDF variable 'lev' warnings.warn( /g/data/hh5/public/apps/cms_conda/envs/analysis3-23.01/lib/python3.9/site-packages/iris/fileformats/_nc_load_rules/helpers.py:660: UserWarning: Ignoring netCDF variable 'lev' invalid units 'level' warnings.warn(msg) /g/data/hh5/public/apps/cms_conda/envs/analysis3-23.01/lib/python3.9/site-packages/iris/analysis/cartography.py:412: UserWarning: Using DEFAULT_SPHERICAL_EARTH_RADIUS. warnings.warn("Using DEFAULT_SPHERICAL_EARTH_RADIUS.")
<xarray.DataArray (lat: 241, lon: 480)> array([[2.23107665e-08, 2.23107665e-08, 2.23107665e-08, ..., 2.23107665e-08, 2.23107665e-08, 2.23107665e-08], [1.78480398e-07, 1.78480398e-07, 1.78480398e-07, ..., 1.78480398e-07, 1.78480398e-07, 1.78480398e-07], [3.56930214e-07, 3.56930214e-07, 3.56930214e-07, ..., 3.56930214e-07, 3.56930214e-07, 3.56930214e-07], ..., [3.56930214e-07, 3.56930214e-07, 3.56930214e-07, ..., 3.56930214e-07, 3.56930214e-07, 3.56930214e-07], [1.78480398e-07, 1.78480398e-07, 1.78480398e-07, ..., 1.78480398e-07, 1.78480398e-07, 1.78480398e-07], [2.23107665e-08, 2.23107665e-08, 2.23107665e-08, ..., 2.23107665e-08, 2.23107665e-08, 2.23107665e-08]]) Dimensions without coordinates: lat, lon
1 2 3 | ### dr4292 - add weights here # hhrT_timeseries = hhrT.mean(dim=["lat","lon"]) #calculate mean to plot for comparison hhrT_timeseries = hhrT.weighted(hr_weights).mean(dim=["lat","lon"]) #calculate mean to plot for comparison |
The data pre regridding
1 2 | hhrT_timeseries.plot(color='tab:blue',label='highres ERA-I',alpha=0.5) plt.legend() |
<matplotlib.legend.Legend at 0x150eebd44f40>
From here the high resolution grid as the input and the low resolution grid as the output are prepared to calculate the regridder
1 | targetgrid
|
<xarray.Dataset> Dimensions: (longitude: 96, latitude: 73, hybrid: 1, t: 1, longitude_1: 96, latitude_1: 72) Coordinates: * longitude (longitude) float32 0.0 3.75 7.5 11.25 ... 348.8 352.5 356.2 * latitude (latitude) float32 -90.0 -87.5 -85.0 -82.5 ... 85.0 87.5 90.0 * hybrid (hybrid) float32 60.0 * t (t) datetime64[ns] 1988-11-01 * longitude_1 (longitude_1) float32 1.875 5.625 9.375 ... 350.6 354.4 358.1 * latitude_1 (latitude_1) float32 -88.75 -86.25 -83.75 ... 83.75 86.25 88.75 Data variables: T (t, hybrid, latitude, longitude) float32 ... U (t, hybrid, latitude, longitude_1) float32 ... V (t, hybrid, latitude_1, longitude) float32 ... Attributes: history: Fri May 22 15:33:09 BST 2020 - XCONV V1.94 03-May-2018
1 | lrgrid = targetgrid.T[0,0,:,:] |
1 2 3 4 5 6 7 8 9 | ### dr4292 - calculate weights for low-res grid tc=iris.load_cube('/g/data/k10/pf4000/UM/um_experiments/climatology/climatology_hourofyear.nc',"T") lat=tc.dim_coords[2] lat.guess_bounds() lon=tc.dim_coords[3] lon.guess_bounds() lr_weights=iris.analysis.cartography.area_weights(tc[0][0],normalize=True) lr_weights=xr.DataArray(lr_weights).rename(dim_0="latitude",dim_1="longitude") lr_weights |
/g/data/hh5/public/apps/cms_conda/envs/analysis3-23.01/lib/python3.9/site-packages/iris/analysis/cartography.py:412: UserWarning: Using DEFAULT_SPHERICAL_EARTH_RADIUS. warnings.warn("Using DEFAULT_SPHERICAL_EARTH_RADIUS.")
<xarray.DataArray (latitude: 73, longitude: 96)> array([[1.23943683e-06, 1.23943683e-06, 1.23943683e-06, ..., 1.23943683e-06, 1.23943683e-06, 1.23943683e-06], [9.91198529e-06, 9.91198529e-06, 9.91198529e-06, ..., 9.91198529e-06, 9.91198529e-06, 9.91198529e-06], [1.98051801e-05, 1.98051801e-05, 1.98051801e-05, ..., 1.98051801e-05, 1.98051801e-05, 1.98051801e-05], ..., [1.98051801e-05, 1.98051801e-05, 1.98051801e-05, ..., 1.98051801e-05, 1.98051801e-05, 1.98051801e-05], [9.91198529e-06, 9.91198529e-06, 9.91198529e-06, ..., 9.91198529e-06, 9.91198529e-06, 9.91198529e-06], [1.23943683e-06, 1.23943683e-06, 1.23943683e-06, ..., 1.23943683e-06, 1.23943683e-06, 1.23943683e-06]]) Dimensions without coordinates: latitude, longitude
1 | hrgrid = hhrT[0] #reduce to grid only |
1 | hrgrid = hrgrid.drop_vars(['time','lev']) #drop extra info |
1 | hrgrid
|
<xarray.DataArray 'ta' (lat: 241, lon: 480)> array([[247.12662, 247.12662, 247.12662, ..., 247.12662, 247.12662, 247.12662], [246.78316, 246.77962, 246.77637, ..., 246.79301, 246.78947, 246.78632], [246.63676, 246.63649, 246.63593, ..., 246.63731, 246.63751, 246.63722], ..., [248.33058, 248.3356 , 248.34175, ..., 248.31961, 248.32268, 248.32622], [249.1894 , 249.19536, 249.20233, ..., 249.17192, 249.1776 , 249.18317], [250.28491, 250.28491, 250.28491, ..., 250.28491, 250.28491, 250.28491]], dtype=float32) Coordinates: * lon (lon) float64 -180.0 -179.2 -178.5 -177.8 ... 177.8 178.5 179.2 * lat (lat) float64 90.0 89.25 88.5 87.75 ... -87.75 -88.5 -89.25 -90.0 Attributes: standard_name: air_temperature long_name: Temperature units: K code: 130 table: 128 MD5: d41d8cd98f00b204e9800998ecf8427e
1 | regridder_con = xe.Regridder(hrgrid, lrgrid, 'conservative', periodic=True) |
/g/data/hh5/public/apps/cms_conda/envs/analysis3-23.01/lib/python3.9/site-packages/xesmf/backend.py:56: UserWarning: Latitude is outside of [-90, 90] warnings.warn('Latitude is outside of [-90, 90]') /g/data/hh5/public/apps/cms_conda/envs/analysis3-23.01/lib/python3.9/site-packages/xesmf/backend.py:56: UserWarning: Latitude is outside of [-90, 90] warnings.warn('Latitude is outside of [-90, 90]')
1 | regridder_bil = xe.Regridder(hrgrid, lrgrid, 'bilinear', periodic=True) |
1 | regridder_con
|
xESMF Regridder Regridding algorithm: conservative Weight filename: conservative_241x480_73x96.nc Reuse pre-computed weights? False Input grid shape: (241, 480) Output grid shape: (73, 96) Periodic in longitude? False
1 | regridder_bil
|
xESMF Regridder Regridding algorithm: bilinear Weight filename: bilinear_241x480_73x96_peri.nc Reuse pre-computed weights? False Input grid shape: (241, 480) Output grid shape: (73, 96) Periodic in longitude? True
Apply both regridders and calculate mean to plot
1 | hr2lr_con = regridder_con(hhrT) |
1 2 3 | ### dr4292 - add weights here # hr2lr_con_mean = hr2lr_con.mean(dim=["latitude","longitude"]) hr2lr_con_mean = hr2lr_con.weighted(lr_weights).mean(dim=["latitude","longitude"]) |
1 | hr2lr_bil = regridder_bil(hhrT) |
1 2 3 | ### dr4292 - add weights here # hr2lr_bil_mean = hr2lr_bil.mean(dim=["latitude","longitude"]) hr2lr_bil_mean = hr2lr_bil.weighted(lr_weights).mean(dim=["latitude","longitude"]) |
1 2 3 4 | hhrT_timeseries.plot(color='tab:blue',label='orig highres ERA-I', alpha=0.5) hr2lr_con_mean.plot(color='tab:blue',label='regridded ERA-I') plt.title('Regridding method conservative') plt.legend() |
<matplotlib.legend.Legend at 0x150ee20c8a00>
1 2 3 4 | hhrT_timeseries.plot(color='tab:blue',label='orig highres ERA-I', alpha=0.5) hr2lr_bil_mean.plot(color='tab:blue',label='regridded ERA-I') plt.title('Regridding method bilinear') plt.legend() |
<matplotlib.legend.Legend at 0x150ee21a8820>
1 |