This code shows as simple way to:
import xarray as xr
import cf_units as cf
import esmlab
from ctsm_py import utils
import numpy as np
import pandas as pd
from netCDF4 import num2date
# some resources for plotting
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
#import hvplot.xarray
#import holoviews as hv
#hv.extension('bokeh')
%matplotlib inline
## User defined options
model = 'CLM5_GSWP3'
var = 'GPP'
suff = 'lnd/proc/tseries/month_1/'
path = '/glade/p/cgd/tss/people/oleson/CLM_LAND_ONLY_RELEASE/CLM5/clm50_r270_1deg_GSWP3V1_iso_newpopd_hist/'+suff
file = path + 'clm50_r270_1deg_GSWP3V1_iso_newpopd_hist.clm2.h0.'+ var +'.185001-201412.nc'
file
'/glade/p/cgd/tss/people/oleson/CLM_LAND_ONLY_RELEASE/CLM5/clm50_r270_1deg_GSWP3V1_iso_newpopd_hist/lnd/proc/tseries/month_1/clm50_r270_1deg_GSWP3V1_iso_newpopd_hist.clm2.h0.GPP.185001-201412.nc'
# ?? Wrapping the time_set_mid utility corrects time being read in by xr.open_dataset, other options ??
# ?? Would be nice to build a way to read in multiple models here, maybe build on Keith L's Example ??
# data sets have lots of coords, variables, and attributes
#ds = xr.open_dataset(file, decode_times=True)
ds = utils.time_set_mid(xr.open_dataset(file, decode_times=True), 'time')
print(ds)
area = ds.area
landfrac = ds.landfrac
da = ds[var] # strips out the extra coords and variables from the dataset
dataset = da.to_dataset() # convert the previously created unnamed DataArray to Dataset with a variable named "array"
<xarray.Dataset> Dimensions: (hist_interval: 2, lat: 192, levdcmp: 25, levgrnd: 25, levlak: 10, lon: 288, time: 1980) Coordinates: * lat (lat) float32 -90.0 -89.057594 -88.11518 ... 89.057594 90.0 * levdcmp (levdcmp) float32 0.01 0.04 0.09 ... 28.870724 41.998436 * levgrnd (levgrnd) float32 0.01 0.04 0.09 ... 28.870724 41.998436 * levlak (levlak) float32 0.05 0.6 2.1 4.6 ... 18.6 25.6 34.325 44.775 * lon (lon) float32 0.0 1.25 2.5 3.75 ... 355.0 356.25 357.5 358.75 * time (time) object 1850-01-16 11:45:00 ... 2014-12-16 12:00:00 Dimensions without coordinates: hist_interval Data variables: BSW (levgrnd, lat, lon) float32 ... DZLAKE (levlak, lat, lon) float32 ... DZSOI (levgrnd, lat, lon) float32 ... GPP (time, lat, lon) float32 ... HKSAT (levgrnd, lat, lon) float32 ... SUCSAT (levgrnd, lat, lon) float32 ... WATSAT (levgrnd, lat, lon) float32 ... ZLAKE (levlak, lat, lon) float32 ... ZSOI (levgrnd, lat, lon) float32 ... area (lat, lon) float32 ... date_written (time) object ... landfrac (lat, lon) float32 ... landmask (lat, lon) float64 ... mcdate (time) float64 ... mcsec (time) float64 ... mdcur (time) float64 ... mscur (time) float64 ... nbedrock (lat, lon) float64 ... nstep (time) float64 ... pftmask (lat, lon) float64 ... time_bounds (time, hist_interval) object 1849-12-31 23:30:00 ... 2015-01-01 00:00:00 time_written (time) object ... Attributes: _NCProperties: version=1|netcdflibversion=4.4... title: CLM History file information comment: NOTE: None of the variables ar... Conventions: CF-1.0 history: Sat Feb 24 14:09:15 2018: /gla... source: Community Land Model CLM4.0 hostname: cheyenne username: oleson version: unknown revision_id: $Id: histFileMod.F90 42903 201... case_title: UNSET case_id: clm50_r270_1deg_GSWP3V1_iso_ne... Surface_dataset: surfdata_0.9x1.25_78pfts_CMIP6... Initial_conditions_dataset: finidat_interp_dest.nc PFT_physiological_constants_dataset: clm5_params.c171117.nc ltype_vegetated_or_bare_soil: 1 ltype_crop: 2 ltype_UNUSED: 3 ltype_landice_multiple_elevation_classes: 4 ltype_deep_lake: 5 ltype_wetland: 6 ltype_urban_tbd: 7 ltype_urban_hd: 8 ltype_urban_md: 9 ctype_vegetated_or_bare_soil: 1 ctype_crop: 2 ctype_crop_noncompete: 2*100+m, m=cft_lb,cft_ub ctype_landice: 3 ctype_landice_multiple_elevation_classes: 4*100+m, m=1,glcnec ctype_deep_lake: 5 ctype_wetland: 6 ctype_urban_roof: 71 ctype_urban_sunwall: 72 ctype_urban_shadewall: 73 ctype_urban_impervious_road: 74 ctype_urban_pervious_road: 75 cft_c3_crop: 1 cft_c3_irrigated: 2 cft_temperate_corn: 3 cft_irrigated_temperate_corn: 4 cft_spring_wheat: 5 cft_irrigated_spring_wheat: 6 cft_winter_wheat: 7 cft_irrigated_winter_wheat: 8 cft_temperate_soybean: 9 cft_irrigated_temperate_soybean: 10 cft_barley: 11 cft_irrigated_barley: 12 cft_winter_barley: 13 cft_irrigated_winter_barley: 14 cft_rye: 15 cft_irrigated_rye: 16 cft_winter_rye: 17 cft_irrigated_winter_rye: 18 cft_cassava: 19 cft_irrigated_cassava: 20 cft_citrus: 21 cft_irrigated_citrus: 22 cft_cocoa: 23 cft_irrigated_cocoa: 24 cft_coffee: 25 cft_irrigated_coffee: 26 cft_cotton: 27 cft_irrigated_cotton: 28 cft_datepalm: 29 cft_irrigated_datepalm: 30 cft_foddergrass: 31 cft_irrigated_foddergrass: 32 cft_grapes: 33 cft_irrigated_grapes: 34 cft_groundnuts: 35 cft_irrigated_groundnuts: 36 cft_millet: 37 cft_irrigated_millet: 38 cft_oilpalm: 39 cft_irrigated_oilpalm: 40 cft_potatoes: 41 cft_irrigated_potatoes: 42 cft_pulses: 43 cft_irrigated_pulses: 44 cft_rapeseed: 45 cft_irrigated_rapeseed: 46 cft_rice: 47 cft_irrigated_rice: 48 cft_sorghum: 49 cft_irrigated_sorghum: 50 cft_sugarbeet: 51 cft_irrigated_sugarbeet: 52 cft_sugarcane: 53 cft_irrigated_sugarcane: 54 cft_sunflower: 55 cft_irrigated_sunflower: 56 cft_miscanthus: 57 cft_irrigated_miscanthus: 58 cft_switchgrass: 59 cft_irrigated_switchgrass: 60 cft_tropical_corn: 61 cft_irrigated_tropical_corn: 62 cft_tropical_soybean: 63 cft_irrigated_tropical_soybean: 64 time_period_freq: month_1 NCO: "4.5.5"
# ?? Would be nice to use cfunits here ??
landUp = area * landfrac * 1e6 # area in km2, not m2
orig_units = cf.Unit(landUp.attrs['units'])
target_units = cf.Unit('m^2')
landUp.attrs['units'] = target_units# rename units
area_wgt = landUp / landUp.sum() # weighting for each grid cell
spy = 365 * 24 * 3600 # Convert to annual fluxes (gC/m2/y), from gc/m2/s
# apply iterates through all variables in a dataset and applies the function to each variable.
# ?? Would be nice to do this differently for fluxes & stocks ??
mean_dataset = dataset.apply(utils.weighted_annual_mean) * spy
mean_dataset.attrs['units'] = 'gC m-2 y-1'
print(mean_dataset)
<xarray.Dataset> Dimensions: (lat: 192, lon: 288, time: 165) Coordinates: * lat (lat) float32 -90.0 -89.057594 -88.11518 ... 89.057594 90.0 * lon (lon) float32 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.25 357.5 358.75 * time (time) object 1850-12-16 12:00:00 ... 2014-12-16 12:00:00 Data variables: GPP (time, lat, lon) float64 0.0 0.0 0.0 0.0 0.0 ... nan nan nan nan Attributes: units: gC m-2 y-1
# ?? This could also be done regionally (see Paul's comments on google doc) ??
mean_wgt = mean_dataset * landUp
ann_glob = mean_wgt.sum(dim=('lat','lon')) * 1e-15 #convert to Pg/y
ann_glob.attrs['units'] = 'Pg C y^-1'
print(ann_glob)
<xarray.Dataset> Dimensions: (time: 165) Coordinates: * time (time) object 1850-12-16 12:00:00 ... 2014-12-16 12:00:00 Data variables: GPP (time) float64 97.99 97.83 97.77 98.42 ... 124.6 121.3 122.3 122.4 Attributes: units: Pg C y^-1
#quick plot , the `;` on the last line supresses text after executing the cell
plt.plot(ann_glob.coords['time.year'], ann_glob[var], '-' );
# Create a nicer figure
fig = plt.figure(figsize=(10, 6))
# Ask, out of a 1x1 grid, the first axes.
ax = fig.add_subplot(1, 1, 1)
# Plot times as x-variable and temperatures as y-variable
ax.plot(ann_glob.coords['time.year'], ann_glob[var], label=model)
# Add some labels to the plot
ax.set_xlabel('Time')
ax.set_ylabel(var+' (' + ann_glob.attrs['units'] +')')
ax.set_title(var, fontdict={'size':16})
ax.legend(loc='upper left');
# making a data array, because cyclic_dataarray currently set up that way
x1= mean_dataset.isel(time=slice(0,10)).mean('time')[var] # Climatology for first 10 years of data
x = mean_dataset.isel(time=slice(-10,None)).mean('time')[var] # Climatology for last 10 years of data
xdiff = x - x1
xinit = x1 * landUp
xfin = x * landUp
xdelta = xdiff * landUp
glob_init = xinit.sum(dim=('lat','lon')) * 1e-15 #convert to Pg/y
glob_fin = xfin.sum(dim=('lat','lon')) * 1e-15 #convert to Pg/y
glob_change = xdelta.sum(dim=('lat','lon')) * 1e-15 #convert to Pg/y
print(glob_init, glob_fin, glob_change)
# Then add cyclic point, ?? is there a better way to do this ??
x1= utils.cyclic_dataarray(x1)
x = utils.cyclic_dataarray(x)
xdiff = utils.cyclic_dataarray(xdiff)
lat = x.coords['lat']
lon = x.coords['lon']
<xarray.DataArray ()> array(97.88563119) <xarray.DataArray ()> array(121.24235673) <xarray.DataArray ()> array(23.35672554)
/glade/u/home/wwieder/miniconda3/envs/python-tutorial/lib/python3.7/site-packages/xarray/core/nanops.py:140: RuntimeWarning: Mean of empty slice return np.nanmean(a, axis=axis, dtype=dtype)
import matplotlib.colors as colors
import cartopy.feature as cfeature
norm = colors.Normalize(vmin=0., vmax=3.5e3)
norm(x.mean())
# modified from https://unidata.github.io/MetPy/latest/examples/Four_Panel_Map.html
crs = ccrs.Robinson()#central_longitude=305.0)
# Function used to create the map subplots
def plot_background(ax):
ax.set_global()
ax.coastlines()
#ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
return ax
levels = np.linspace(0, 3500, 11)
fig = plt.figure(figsize=(5, 3.5),constrained_layout=True);
fig, axarr = plt.subplots(nrows=2, ncols=1, figsize=(10, 6), constrained_layout=True,
subplot_kw={'projection': crs})
axlist = axarr.flatten()
for ax in axlist:
plot_background(ax)
ax.set_extent([-155, 60, -25, 60]) # I really don't understand how this works??
# ax.set_extent([-150, 60, -25, 60])
# upper plot
cf0 = axlist[0].contourf(lon, lat, x1, levels=levels, norm=norm, cmap='gist_earth_r',
transform=ccrs.PlateCarree());
axlist[0].set_title('Initial '+model+' '+var, fontsize=18)
axlist[0].set_axis_off()
# lower plot
cf1 = axlist[1].contourf(lon, lat, x, levels=levels, norm=norm, cmap='gist_earth_r',
transform=ccrs.PlateCarree());
axlist[1].set_title('Final '+model+' '+var, fontsize=18);
axlist[1].set_axis_off()
# add common colorbar
cbar = fig.colorbar(cf1, ax=axlist.ravel().tolist(),
orientation='vertical', shrink=0.80, pad=0)
#cbar.set_label('gC/m2/y', size=18, orientation='horizontal')
cbar.ax.tick_params(labelsize=16)
ax = cbar.ax
ax.text(0.2,-0.1,mean_dataset.attrs['units'], size=16,rotation=0)
plt.show();
<Figure size 360x252 with 0 Axes>
# Now make a difference map
norm = colors.Normalize(vmin=-1000., vmax=1000.)
norm(xdiff.mean())
levels = np.linspace(-1e3, 1e3, 11)
ax = plt.axes(projection=ccrs.Robinson())
ax.set_global()
ax.coastlines()
ax.set_extent([-155, 60, -25, 60]) # I really don't understand how this works??
cs = ax.contourf(lon, lat, xdiff, levels=levels, norm=norm, cmap='BrBG',
transform=ccrs.PlateCarree());
# add colorbar.
cbar = plt.colorbar(cs, orientation='horizontal')
cbar.set_label('gC/m2/y')
ax.set_title(var + " change")
plt.show();