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)

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

##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.34989166, 22.14568471, 21.66339844, ..., 49.19934755,
         48.19983052, 47.56431138],
        [21.7082633 , 21.17000119, 20.13788047, ..., 45.64615776,
         45.32941988, 44.5238861 ],
        [20.61434703, 20.30073408, 19.39063505, ..., 40.45091568,
         39.51188793, 40.08718096],
        ...,
        [15.20498414, 16.70405236, 16.6958771 , ..., 14.96512996,
         16.09915845, 16.28375345],
        [15.95014737, 15.44401933, 16.36107178, ..., 13.13411618,
         13.86219912, 14.25409531],
        [14.45032494, 14.94088309, 11.27504841, ...,  9.55777922,
         11.2378108 , 13.30571821]]])

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

../../_images/39f348d2c1cae3b439e128b4b0c5de9093fb22c746d2335721bacdedfb882435.png