<jupyter_start><jupyter_text>Calculate historical Vapor Pressure Deficit (VPD) grids.Do this for gridded ERA-Interim data using the estimated RH% field. VPD is a simple forula. \begin{align}VPD & = e_{sat} - e_{actual}\end{align}where: \begin{align}e_{actual}=\frac{e_{sat}*RH)}{100}\end{align} where the saturation water vopor pressure ($e_{sat}$) is expressed with the Teten's formula:Teton's formula, eq. 7.5 https://www.ecmwf.int/sites/default/files/elibrary/2015/9211-part-iv-physical-processes.pdf\begin{align}e_{sat}(T) & = a_{1}exp(a_{3}\frac{T - T_{0}}{T-a_{4}})\\\\a_{1} & = 611.21Pa \\a_{3} & = 17.502 \\a_{4} & = 32.19 \\T_{0} & = 273.16\end{align}<jupyter_code>import numpy as np
from netCDF4 import Dataset
import os
import glob

def calculate_es(T):
    """Calculate saturation vapor pressure using Tetons formula"""
    a1 = 611.21 # Pa
    a3 = 17.502  
    a4 = 32.19
    T0 = 273.16 # K
    e_sat = a1 * np.exp( a3 * (T-T0)/(T-a4) )
    
    return e_sat
def get_nc(var, year, dataDir): 
    loadFile = os.path.join(dataDir, var + "_" + str(year) + ".nc")
    nc = Dataset(loadFile)
    vals = nc.variables[var][:]
    t = nc.variables["time"][:]
    lon = nc.variables["longitude"][:]
    lat = nc.variables["latitude"][:]
    nc.close()
    return vals, t, lon, lat
def write_VPD_nc(VALS, t, x, y, year, dataDir):
    
    outputFile = os.path.join(dataDir, "VPD_" + str(year) + ".nc")

    ncFile = Dataset(outputFile, 'w', format='NETCDF4')
    ncFile.description = 'VPD (saturation pressure relative to water, Tetons eq.)'
    ncFile.location = 'Global'
    ncFile.createDimension('time',  len(t) )
    ncFile.createDimension('latitude', len(y) )
    ncFile.createDimension('longitude', len(x) )

    VAR_ = ncFile.createVariable("VPD", 'f4',('time', 'latitude','longitude'))
    VAR_.long_name = "Vapor Pressure Deficit"
    VAR_.units = "Pa"

    # Create time variable
    time_ = ncFile.createVariable('time', 'i4', ('time',))
    time_.units = "hours since 1900-01-01 00:00:0.0"
    time_.calendar = "gregorian"
    
    # create lat variable
    latitude_ = ncFile.createVariable('latitude', 'f4', ('latitude',))
    latitude_.units = 'degrees_north'

    # create longitude variable
    longitude_ = ncFile.createVariable('longitude', 'f4', ('longitude',))
    longitude_.units = 'degrees_east'

    # Write the actual data to these dimensions
    VAR_[:]       = VALS
    time_[:]      = t
    latitude_[:]  = y
    longitude_[:] = x

    ncFile.close()<jupyter_output><empty_output><jupyter_text>Use above functions to create VPD file for era-interim.<jupyter_code>dataDir = os.path.join("..","..", "metSpreadData", "ERA-Interim", 'merged_time')
year = "1983-2017"
    
# Get temperature to calculate es
t2m,t,x,y = get_nc("t2m", year, dataDir)
    
# Get RH% to calculate e_actual
RH,t,x,y = get_nc("RH", year, dataDir)

# get saturation and actual vapor pressure 
e_sat = calculate_es(t2m)
e_actual = e_sat * (RH / 100.)

# Calculate VPD using what is known as subtraction
VPD = e_sat - e_actual

write_VPD_nc(VPD, t, x, y, year, dataDir)
import cdo as cdo
cdo = cdo.Cdo()

dataDir = os.path.join("..","..", "metSpreadData", "ERA-Interim", 'merged_time')
common_grid_txt = os.path.join("..","..", "metSpreadData", "ERA-Interim", "COMMON_GRID.txt")
f_in = os.path.join(dataDir, "VPD_1983-2017.nc")
f_out = os.path.join("..","..", "metSpreadData", "ERA-Interim", 'merged_t_COMMON_GRID', "VPD_1983-2017.nc")
cdo.remapbil(common_grid_txt, input=f_in, output=f_out, options="-b F64")<jupyter_output><empty_output><jupyter_text>Use cdo remapbil,COMMON_GRID.txt to create the needed file for comparing to CMIP5. Move that file to the metSpread directory. - Done most recently for 1983-2017 VPD data via bash cdo on Nov 24 2018  Do the same for CMIP5 output.  NOTE:This script creates VPD files for merged_common_grid variables only, so only run after Python/format_cmip5_data.py has been most recently executed.  TODO: - Update for new files and the fact that we have huss now and do longer need to estimate using RH%- Be more careful in the VPD file creation such that the month spaned that the data represent time span is better preserved. The output and calculations are correct but there is something sloppy happenning with the time stamp when I quickly rewrite the netcdf file using write_VPD_nc which takes t.units and time.data from a tas file. - Use available files, loop through them, vs. checking for a T and RH file for each attempted model.<jupyter_code>def get_CMIP_nc(var, loadFile): 
    
    nc = Dataset(loadFile)
    vals = np.array(nc.variables[var])
    t = np.array(nc.variables["time"])
    t_units = nc.variables["time"].units
    lon = np.array(nc.variables["lon"])
    lat = np.array(nc.variables["lat"])
    nc.close()
    
    return vals, t, t_units, lon, lat
def write_VPD_CMIP5_nc(VALS, model, rcp, t, t_units, x, y, dataDir):
    """
    Needs to take rcp files that have been merged with history files, so that these operations 
    can be done only onces and not repeated at several points on the data pipeline. 
        
        DataSource: ../r1i1p1_rcp_COMMON_GRID/
        
    """
    
    outputFile = os.path.join(dataDir, "VPD_Amon_"+model+"_"+rcp+"_r1i1p1_198301-210012.nc")

    ncFile = Dataset(outputFile, 'w', format='NETCDF4')
    ncFile.description = 'VPD (saturation pressure relative to water, Tetons eq.)'
    ncFile.location = 'Global'
    ncFile.createDimension('time',  len(t) )
    ncFile.createDimension('lat', len(y) )
    ncFile.createDimension('lon', len(x) )

    VAR_ = ncFile.createVariable("VPD", 'f4',('time', 'lat','lon'))
    VAR_.long_name = "Vapor Pressure Deficit"
    VAR_.units = "Pa"

    # Create time variable
    time_ = ncFile.createVariable('time', 'i4', ('time',))
    time_.units = t_units
    #time_.calendar = "gregorian"
    
    # create lat variable
    latitude_ = ncFile.createVariable('lat', 'f4', ('lat',))
    latitude_.units = 'degrees_north'
    latitude_.axis = "Y"

    # create longitude variable
    longitude_ = ncFile.createVariable('lon', 'f4', ('lon',))
    longitude_.units = 'degrees_east'
    longitude_.axis = "X"

    # Write the actual data to these dimensions
    VAR_[:]       = VALS
    time_[:]      = t
    latitude_[:]  = y
    longitude_[:] = x

    ncFile.close()<jupyter_output><empty_output><jupyter_text>Loop through the unique models based on the data that are actually available at the common_grid_level<jupyter_code>dataDir = os.path.join("..","..", "metSpreadData", "getCMIP5", 'r1i1p1_rcp_COMMON_GRID')

all_model_names = []
all_files = glob.glob(dataDir+"/tas_Amon*")
for f in all_files :
    all_model_names.append(f.split("/")[5].split("_")[2])
unique_model_names = np.unique(all_model_names)
unique_model_names
dataDir = os.path.join("..","..", "metSpreadData", "getCMIP5", 'r1i1p1_rcp_COMMON_GRID')

print(dataDir)

# Need VPD for both rcps
rcps = ['rcp45', 'rcp85']

# list the files 

for rcp in rcps : 
    
    print("Working on: " + rcp)
    
    for model in unique_model_names :
    
        print("start----------------------------")
        print(model)
                
        # see if Temperature and RH file exist, both are required
        # TODO: replace hurs (RH) with use of estimated huss (e) and that will increase accuracy
        # TODO: and simpliciity though make it slightly different from how things are done for
        # TODO: calculating era-interim VPD. 
        T_file = os.path.join(dataDir, "tas_Amon_"+model+"_"+rcp+"_r1i1p1_198301-210012.nc")
        RH_file = os.path.join(dataDir, "hurs_Amon_"+model+"_"+rcp+"_r1i1p1_198301-210012.nc")

        # If these files exist (they don't for all models and rcps), continue
        if os.path.isfile(T_file) & os.path.isfile(RH_file) : 
            
            # Get the values to calculate e
            tas,t,t_units,x,y = get_CMIP_nc('tas', T_file)
            RH,t,_, x,y  = get_CMIP_nc('hurs', RH_file)

            # get saturation and actual vapor pressure 
            e_sat = calculate_es(tas)
            e_actual = e_sat * (RH / 100.)

            # Calculate VPD using what is known as subtraction
            VPD = e_sat - e_actual

            write_VPD_CMIP5_nc(VPD, model, rcp, t, t_units, x, y,  dataDir)
            
            print("Writing VPD esetimate")
            
        else :
            
            print("*********************************")
            print("Missing a file: ")
            print("T_file present: " + str(os.path.isfile(T_file)))
            print("RH_file present: " + str(os.path.isfile(RH_file)))
            print("*********************************")
            
        print("------------------------------end")
        print(" ")<jupyter_output>../../metSpreadData/getCMIP5/r1i1p1_rcp_COMMON_GRID
Working on: rcp45
start----------------------------
ACCESS1-0
Writing VPD esetimate
------------------------------end
 
start----------------------------
ACCESS1-3
Writing VPD esetimate
------------------------------end
 
start----------------------------
CCSM4
Writing VPD esetimate
------------------------------end
 
start----------------------------
CMCC-CESM
*********************************
Missing a file: 
T_file present: False
RH_file present: False
*********************************
------------------------------end
 
start----------------------------
CMCC-CM
*********************************
Missing a file: 
T_file present: True
RH_file present: False
*********************************
------------------------------end
 
start----------------------------
CMCC-CMS
*********************************
Missing a file: 
T_file present: True
RH_file present: False
*********************************
------------------------------end
 
sta[...]
{'fraction_non_alphanumeric': 0.19158110882956877, 'fraction_numerical': 0.022484599589322384, 'mean_word_length': 2.8624107850911975, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '411506', 'license': 'MIT', 'n_tokens_mistral': 3466, 'n_tokens_neox': 3201, 'n_words': 838, 'path': 'Python/calculate_gridded_data_VPD.ipynb', 'repo_name': 'stevenjoelbrey/metSpread', 'size': '9740'}