Calculating Climate Suitable Days using netCDF4 and Python?

Question

This is a follow up to Using netCDF4 Python climate algorithm?

I’m relatively novice in using netCDF files, so I’m wondering if someone can check my logic on my process of delineating what has been defined as “Vector Genesis Days” for epidemiology research. The criteria to be met for each day are:

---

18°C ≤ T ≤ 32°C
 20% ≤ Q ≤ 80%
1.5 mm ≤ R ≤ 20 mm 

Where, 
T is Temperature
Q is Relative Humidity
R is Rainfall

---

I’m attempting to iterate through each of these variables by day to get a count of the number of days that meet the specified parameters per year. I want to output this final annual data (1995-2009) into centroids representative of the current netCDF file grids(if possible, I did it is a “rough” manual manner) or rasters if easier.

Some background on the data: 
NCEP Reanalysis I - downloaded at [ESRL NOAA PSD][2]
Data coverage:
Study Area – Orissa (Odisha), India
2.5° by 2.5° grid … ~69 mi (111 km)
9 x 11 Arrays:
Latitude: 12N – 28N
Longitude: 74E – 94E

---

Variables with units:
Temperature Maximum (tmax) – K (Kelvin)
Temperature Minimum (tmin) – K(Kelvin)
Precipitation Rate (prate) – kg/m2/s
Temperature Daily Average (temp) – K(Kelvin)
Specific Humidity (shum) – kg/kg
Pressure (pres) – Pascals

---

A few conversions were made to create the needed daily mm/day in rainfall from precipitation rate, as well as, conversion of specific humidity to relative humidity (that’s where temperature daily average and pressure needed to be utilized). In addition, a few conversions were used for Celsius to Kelvin, etc.

My current code in Python is:

""" Leah Bowyer Python Code: Calculate Vector Genesis Days >>> Convert to CSV """

import numpy as np
import netCDF4 as ncdf

prate=ncdf.Dataset('N:/NCEP Reanalysis_India/prate.sfc.gauss.1995.nc','r')
tmax=ncdf.Dataset('N:/NCEP Reanalysis_India/tmax.2m.gauss.1995.nc','r')
tmin=ncdf.Dataset('N:/NCEP Reanalysis_India/tmin.2m.gauss.1995.nc','r')
shum=ncdf.Dataset('N:/NCEP Reanalysis_India/shum.2m.gauss.1995.nc','r+')  ## 'r+' when converting humidity file
temp=ncdf.Dataset('N:/NCEP Reanalysis_India/temp.2m.gauss.1995.nc','r')
press=ncdf.Dataset('N:/NCEP Reanalysis_India/pres.2m.gauss.1995.nc','r')

es=6.112 * np.exp((17.67 * (temp.variables['air'][:] - 273.15)) / ((temp.variables['air'][:] - 273.15) + 243.5))
e=shum.variables['shum'][:] *( press.variables['pres'][:] / 100)/ (0.378 * shum.variables['shum'][:] + 0.622)
rhum=(e/es) *100
print shum.variables['shum'][:] ##check values before conversion
shum.variables['shum'][:]=rhum[:]

print shum.variables['shum'][:] ##check values after conversion


##check climate parameters
a=np.logical_and(shum.variables['shum'][:] <=80, shum.variables['shum'][:] >=20)
b=np.logical_and(prate.variables['prate'][:]<= 2.31481481e-03, prate.variables['prate'][:]>= 1.73611111e-05)
c=np.logical_and(tmax.variables['tmax'][:]<=305.15, tmax.variables['tmax'][:]>=291.15)
d=np.logical_and(tmin.variables['tmin'][:]<=305.15, tmin.variables['tmin'][:]>=291.15)

##sum days satisfy daily parameter
Q = sum(a)
R = sum(b)
sum(c)
sum(d)


T=np.logical_and(c,d)
precip_satisfied=np.logical_and(a,b)
criteria_satisfied=np.logical_and(T,precip_satisfied)

sum(criteria_satisfied) ##check logic in output

##output array into .csv file with parameter qualifications (daily value) by column
output = np.column_stack((sum(T).flatten(),Q.flatten(),R.flatten(),sum(criteria_satisfied).flatten()))
np.savetxt('N:\\Python Results\\2011_India.csv',output,delimiter=',')

print sum(criteria_satisfied)

Some things to note:

I’d like to output the individual parameter results to do some controls on the overall “days” count. I’ve also only coded to handle one year at a time, but this data CAN BE downloaded continuously. Also, I’ve exported the resulting array as a .csv file that I joined to an already created centroid file. I know that this is not the most efficient way.

My code does run with somewhat reasonable results, but with some analysis it does not check out correctly.

---