Total Column Precipitable Water (TCPW)

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)

ds

#### 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_data

##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

TCPW

array([[[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 ()