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


```python
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
```


```python
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
```


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

### Use above functions to create VPD file for era-interim. 


```python
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)
```


```python
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")
```




    '../../metSpreadData/ERA-Interim/merged_t_COMMON_GRID/VPD_1983-2017.nc'



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. 


```python
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
```


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

### Loop through the unique models based on the data that are actually available at the common_grid_level 


```python
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)
```


```python
unique_model_names
```




    array(['ACCESS1-0', 'ACCESS1-3', 'CCSM4', 'CMCC-CESM', 'CMCC-CM',
           'CMCC-CMS', 'CNRM-CM5', 'CSIRO-Mk3-6-0', 'CanESM2', 'FGOALS-g2',
           'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H',
           'GISS-E2-H-CC', 'GISS-E2-R', 'GISS-E2-R-CC', 'HadGEM2-AO',
           'HadGEM2-CC', 'HadGEM2-ES', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR',
           'IPSL-CM5B-LR', 'MIROC-ESM', 'MIROC-ESM-CHEM', 'MIROC5',
           'MPI-ESM-LR', 'MPI-ESM-MR', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M',
           'NorESM1-ME', 'bcc-csm1-1-m', 'inmcm4'], dtype='<U14')




```python
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(" ")
```

    ../../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
     
    start----------------------------
    CNRM-CM5
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    CSIRO-Mk3-6-0
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    CanESM2
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    FGOALS-g2
    *********************************
    Missing a file: 
    T_file present: True
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    GFDL-CM3
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GFDL-ESM2G
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GFDL-ESM2M
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GISS-E2-H
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GISS-E2-H-CC
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GISS-E2-R
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GISS-E2-R-CC
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    HadGEM2-AO
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    HadGEM2-CC
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    HadGEM2-ES
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    IPSL-CM5A-LR
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    IPSL-CM5A-MR
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    IPSL-CM5B-LR
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MIROC-ESM
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MIROC-ESM-CHEM
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MIROC5
    *********************************
    Missing a file: 
    T_file present: True
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    MPI-ESM-LR
    *********************************
    Missing a file: 
    T_file present: True
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    MPI-ESM-MR
    *********************************
    Missing a file: 
    T_file present: True
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    MRI-CGCM3
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MRI-ESM1
    *********************************
    Missing a file: 
    T_file present: False
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    NorESM1-M
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    NorESM1-ME
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    bcc-csm1-1-m
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    inmcm4
    Writing VPD esetimate
    ------------------------------end
     
    Working on: rcp85
    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: True
    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
     
    start----------------------------
    CNRM-CM5
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    CSIRO-Mk3-6-0
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    CanESM2
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    FGOALS-g2
    *********************************
    Missing a file: 
    T_file present: True
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    GFDL-CM3
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GFDL-ESM2G
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GFDL-ESM2M
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GISS-E2-H
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GISS-E2-H-CC
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GISS-E2-R
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    GISS-E2-R-CC
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    HadGEM2-AO
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    HadGEM2-CC
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    HadGEM2-ES
    *********************************
    Missing a file: 
    T_file present: False
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    IPSL-CM5A-LR
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    IPSL-CM5A-MR
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    IPSL-CM5B-LR
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MIROC-ESM
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MIROC-ESM-CHEM
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MIROC5
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MPI-ESM-LR
    *********************************
    Missing a file: 
    T_file present: True
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    MPI-ESM-MR
    *********************************
    Missing a file: 
    T_file present: True
    RH_file present: False
    *********************************
    ------------------------------end
     
    start----------------------------
    MRI-CGCM3
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    MRI-ESM1
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    NorESM1-M
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    NorESM1-ME
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    bcc-csm1-1-m
    Writing VPD esetimate
    ------------------------------end
     
    start----------------------------
    inmcm4
    Writing VPD esetimate
    ------------------------------end
     


c3d1c77c54c7e391abf6d5e8d2b9fe1625e5efb0
ipynb
Jupyter Notebook
Python/calculate_gridded_data_VPD.ipynb
stevenjoelbrey/metSpread
38d667f0e2f58563345fed14132bb5a6eb7022af
MIT
Python/calculate_gridded_data_VPD.ipynb
stevenjoelbrey/metSpread
38d667f0e2f58563345fed14132bb5a6eb7022af
MIT
Python/calculate_gridded_data_VPD.ipynb
stevenjoelbrey/metSpread
38d667f0e2f58563345fed14132bb5a6eb7022af
MIT
Qwen/Qwen-72B
1. YES
2. YES