Calculate TCPW on a grid of netCDF data using MetPy.
import metpy.calc as mpcalc
import xarray as xr
import numpy as np
from metpy.calc import cape_cin, surface_based_cape_cin, dewpoint_from_specific_humidity, parcel_profile,relative_humidity_from_specific_humidity,most_unstable_cape_cin, precipitable_water
from metpy.units import units
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature## opening NetCDF file using xarray
ds = xr.open_dataset("NETCDF_FILE.nc", decode_times=True)dsLoading...
#### making a function to slice the xarray dataset according to our need.
def slicer (data,lat1, lat2, lon1, lon2, time1,time2) :
sliced_data = data.sel(lat =slice(lat1, lat2), lon = slice(lon1, lon2),time = slice(time1, time2))
return sliced_data#slicing the data for CONUS only
new_data = slicer(ds,23.5,50.5,-125.5,-66.5, ds.time[0], ds.time[0])new_dataLoading...
##extracting temperature, pressure and specific humidity from the dataset
p =new_data.lev*units('hPa')
T = new_data.T
sh = new_data.QV*units('dimensionless')
TCPW = np.zeros((1, 28, 60),dtype=float)
###calculating TCPW in multidimension
for h in range (len(new_data.time)):
for i in range (len(new_data.lat)):
for j in range(len(new_data.lon)):
TC = (T[h,:,i,j]-(273.15))*units('degC')
Td = dewpoint_from_specific_humidity(np.array(p)*units('hPa'),TC, sh[h,:,i,j])
TPW_S = precipitable_water(np.array(p)*units('hPa'), np.array(Td)*units('degC'))
TCPW[h,i,j] = TPW_S.magnitude
TCPWarray([[[22.33904328, 22.13514006, 21.65314758, ..., 49.16251238,
48.16421028, 47.52806178],
[21.69776415, 21.15994415, 20.12837013, ..., 45.6115163 ,
45.29600502, 44.49049573],
[20.60370483, 20.2902779 , 19.38095387, ..., 40.41971515,
39.48164667, 40.0572096 ],
...,
[15.19375194, 16.69203318, 16.68316602, ..., 14.95516252,
16.0881911 , 16.27202346],
[15.93827239, 15.43163369, 16.34822017, ..., 13.12592686,
13.85348823, 14.24460194],
[14.43932076, 14.92929048, 11.26527078, ..., 9.54927061,
11.22804845, 13.29490269]]], shape=(1, 28, 60))dataproj = ccrs. PlateCarree ()
# # Plot projection
# # The look you want for the view.
plotproj = ccrs. PlateCarree ()
fig=plt.figure(1, figsize=(15.,12.))
ax=plt.subplot(111,projection=plotproj)
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
ax.add_feature(cfeature.STATES, linewidth=0.5)
plt.title("Total Column Precipitable Water Calculated Using Metpy (mm)",size = 30)
plt.contourf (new_data.lon, new_data.lat, TCPW[0,:,:],levels =np.arange(0,70,5),cmap = "PuBuGn", transform=dataproj,extend = "max")
plt.colorbar (orientation = "horizontal", pad=0.01).ax.tick_params(labelsize=20)
plt. show ()