import xarray as xr
from matplotlib.pyplot import colorbar
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
import numpy as np
import pandas as pd
# Load the wilting point and field capacity datasets
thmin1 = xr.open_dataset('./thmin1.nc')
thmin2 = xr.open_dataset('./thmin2.nc')
thmin3 = xr.open_dataset('./thmin3.nc')
thmax1 = xr.open_dataset('./thmax1.nc')
thmax2 = xr.open_dataset('./thmax2.nc')
thmax3 = xr.open_dataset('./thmax3.nc')
# load the EFAS volumetric soil moisture dataset
ds = xr.open_dataset('./eud_2019013000.nc')
# plotting function
def make_cmap(colors, position=None):
'''
make_cmap takes a list of tuples which contain RGB values.
The RGB values are [0 to 255]
make_cmap returns a cmap with equally spaced colors.
Arrange your tuples so that the first color is the lowest value for the
colorbar and the last is the highest.
position contains values from 0 to 1 to dictate the location of each color.
'''
import matplotlib as mpl
import numpy as np
import sys
bit_rgb = np.linspace(0,1,256)
if position == None:
position = np.linspace(0,1,len(colors))
else:
if len(position) != len(colors):
sys.exit("position length must be the same as colors")
elif position[0] != 0 or position[-1] != 1:
sys.exit("position must start with 0 and end with 1")
for i in range(len(colors)):
colors[i] = (bit_rgb[colors[i][0]],
bit_rgb[colors[i][1]],
bit_rgb[colors[i][2]])
cdict = {'red':[], 'green':[], 'blue':[]}
for pos, color in zip(position, colors):
cdict['red'].append((pos, color[0], color[0]))
cdict['green'].append((pos, color[1], color[1]))
cdict['blue'].append((pos, color[2], color[2]))
cmap = mpl.colors.LinearSegmentedColormap('my_colormap',cdict,256)
return cmap
# There are 3 layers of Volumetric Soil Moisture. We need to get the mean across all timesteps
mean_vsw = ds.vsw.mean('step') # Mean Volumetric Soil Moisture for all time steps
th1 = mean_vsw[0,:,:] # Mean Volumetric Soil Moisture for Layer 1
th2 = mean_vsw[1,:,:] # Mean Volumetric Soil Moisture for Layer 2
# We need to have the values for the first and second soil depths
# Layer 2 includes Layer 1, so we need to remove 2 from 1
# All layers are the same through out so we need to use step 1
sd1 = ds.sod[0,0,:,:].values
sd2 = ds.sod[0,1,:,:].values - ds.sod[0,0,:,:].values
# Compute the mean wilting point
data=(thmin1.wiltingpoint.values * sd1 + thmin2.wiltingpoint.values * sd2)/(sd1 + sd2)
thmin = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='thmin')
ds['thmin'] = thmin
# Compute the field capacity
data=(thmax2.fieldcapacity.values * sd1 + thmax2.fieldcapacity.values * sd2)/(sd1 + sd2)
thmax = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='thmax')
ds['thmax'] = thmax
# Compute the mean soil moisture for layer 1 and 2
data=((th1 * sd1) + (th2 * sd2))/(sd1 +sd2)
smtot = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='smtot')
ds['smtot'] = smtot
# Compute the soil wetness index
data=(((smtot - thmin)/(thmax - thmin)))
SM = xr.DataArray(data,dims=(ds.y.name,ds.x.name),name='SM')
ds['SM'] = SM
# Create a list of RGB colors for the plotting function
colors = [(64,0,3), (133,8,3), (249,0,0), (250,106,0), (249,210,0), (189,248,255), (92,138,255), (0,98,255), (0,0,255), (0,0,88)]
# Create an array or list of positions from 0 to 1.
position = [0, 0.2, 0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
cmap = make_cmap(colors, position=position)
# define the coordinate reference system for the EFAS data
crs = ccrs.LambertAzimuthalEqualArea(central_longitude=10,central_latitude=52,false_easting=4321000,false_northing=3210000)
fig, ax = plt.subplots(1,1,subplot_kw={'projection': crs}, figsize=(15,15) )
ax.gridlines(crs=crs, linestyle="-")
ax.coastlines()
ax.add_feature(cf.BORDERS)
sc = ds["SM"].plot(ax=ax,cmap=cmap,add_colorbar=False)
cbar = plt.colorbar(sc, shrink=.5,)
cbar.set_label("Soil Wetness Index") |