In [1]:
%run basics
%matplotlib
import pytz
from pysolar import GetAltitude
from scipy.interpolate import InterpolatedUnivariateSpline
In [2]:
def get_Fsd(erai_file,times_3hr,times_tts,fire,erai):
"""
Purpose:
Extracts a time series of incoming solar radiation for a given location (specified
as latitude and longitude) from an ERA Interim re-analysis file.
Usage:
Fsd = get_FSD(erai_file,times_3hr,times_tts,fire,erai)
where
erai_file is the netCDF file retrieved from the ERAI site
times_3hr is a dictionary of times at the ERAI time step (3 hourly)
times_tts is a dictionary of times at the tower time step (30 or 60 minutes)
fire is a dictionary of information about thye fire instance
erai is a dictionary of information about the ERAI data
Author: PRI
Date: March 2015
"""
# local pointers
elat = erai["latitude"]
elon = erai["longitude"]
dt_erai_utc_cor = times_3hr["dt_erai_utc_cor"]
erai_time_3hr = times_3hr["erai_time_3hr"]
dt_erai_utc_tts = times_tts["dt_erai_utc_tts"]
erai_time_tts = times_tts["erai_time_tts"]
flati = fire["lat_index"]
floni = fire["lon_index"]
# get the solar altitude, we will use this later to interpolate the ERA Interim solar
# data from the ERA-I 3 hour time step to the tower time step.
# NOTE: alt_solar is in degrees
alt_solar_3hr = numpy.array([GetAltitude(elat,elon,dt) for dt in dt_erai_utc_cor])
# get the solar altitude at the tower time step
alt_solar_tts = numpy.array([GetAltitude(elat,elon,dt) for dt in dt_erai_utc_tts])
idx = numpy.where(alt_solar_tts<=0)[0]
alt_solar_tts[idx] = float(0)
# Interpolate the 3 hourly accumulated downwelling shortwave to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Fsd_3d = erai_file.variables["ssrd"][:,:,:]
Fsd_accum = Fsd_3d[:,flati,floni]
# Downwelling shortwave in ERA-I is a cummulative value that is reset to 0 at 0300 and 1500 UTC.
# Here we convert the cummulative values to 3 hourly values.
Fsd_erai_3hr = numpy.ediff1d(Fsd_accum,to_begin=0)
idx = numpy.where((times_3hr["hour_utc"]==3)|
(times_3hr["hour_utc"]==15))[0]
Fsd_erai_3hr[idx] = Fsd_accum[idx]
Fsd_erai_3hr = Fsd_erai_3hr/float(10800)
# normalise the ERA-I downwelling shortwave by the solar altitude
# clamp solar altitude to a minimum value to avoid numerical problems
# when alt_solar is close to 0erai_file,times_3hr,times_tts,fire_info
alt_solar_limit = float(5)*numpy.ones(len(alt_solar_3hr))
sa = numpy.where(alt_solar_3hr<=float(5),alt_solar_limit,alt_solar_3hr)
coef_3hr = Fsd_erai_3hr/numpy.sin(numpy.deg2rad(sa))
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, coef_3hr, k=1)
# get the coefficient at the tower time step
coef_tts = s(erai_time_tts)
# get the downwelling solar radiation at the tower time step for this latitude and longitude
Fsd = {}
Fsd["data"] = coef_tts*numpy.sin(numpy.deg2rad(alt_solar_tts))
Fsd["attr"] = qcutils.MakeAttributeDictionary(long_name="Downwelling short wave radiation",units="W/m2")
return Fsd
In [3]:
def get_Ta(erai_file,times_3hr,times_tts,fire):
# local pointers
flati = fire["lat_index"]
floni = fire["lon_index"]
erai_time_3hr = times_3hr["erai_time_3hr"]
erai_time_tts = times_tts["erai_time_tts"]
# Interpolate the 3 hourly air temperature to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Ta_3d = erai_file.variables["t2m"][:,:,:]
Ta_erai_3hr = Ta_3d[:,flati,floni] - 273.15
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Ta_erai_3hr, k=2)
# get the air temperature at the tower time step
Ta = {}
Ta["data"] = s(erai_time_tts)
Ta["attr"] = qcutils.MakeAttributeDictionary(long_name="Air temperature",units="C")
return Ta
In [4]:
def get_RH(erai_file,times_3hr,times_tts,fire):
# local pointers
flati = fire["lat_index"]
floni = fire["lon_index"]
erai_time_3hr = times_3hr["erai_time_3hr"]
erai_time_tts = times_tts["erai_time_tts"]
# Interpolate the 3 hourly air temperature to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Ta_3d = erai_file.variables["t2m"][:,:,:]
Ta_erai_3hr = Ta_3d[:,flati,floni] - 273.15
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Ta_erai_3hr, k=2)
Ta_erai_tts = s(erai_time_tts)
# Interpolate the 3 hourly dew point temperature to the tower time step
Td_3d = erai_file.variables["d2m"][:,:,:]
Td_erai_3hr = Td_3d[:,flati,floni] - 273.15
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, Td_erai_3hr, k=2)
# get the dew point temperature at the towespeedr time step
Td_erai_tts = s(erai_time_tts)
# get the relative humidity
es_erai_tts = mf.es(Ta_erai_tts)
e_erai_tts = mf.es(Td_erai_tts)
RH = {}
RH["data"] = float(100)*e_erai_tts/es_erai_tts
RH["attr"] = qcutils.MakeAttributeDictionary(long_name="Relative humidity",units="percent")
return RH
In [5]:
def get_Ws_and_Wd(erai_file,times_3hr,times_tts,fire):
# local pointers
flati = fire["lat_index"]
floni = fire["lon_index"]
erai_time_3hr = times_3hr["erai_time_3hr"]
erai_time_tts = times_tts["erai_time_tts"]
# Interpolate the 3 hourly U and V components to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
# U first ...
U_3d = erai_file.variables["u10"][:,:,:]
U_erai_3hr = U_3d[:,flati,floni]
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, U_erai_3hr, k=2)
# get the U component at the tower time step
U_erai_tts = s(erai_time_tts)
# ... then V
V_3d = erai_file.variables["v10"][:,:,:]
V_erai_3hr = V_3d[:,flati,floni]
# get the spline interpolation function
s = InterpolatedUnivariateSpline(erai_time_3hr, V_erai_3hr, k=2)
# get the soil moisture at the tower time step
V_erai_tts = s(erai_time_tts)
# now get the wind speed and direction
Ws = {}
Ws["data"] = numpy.sqrt(U_erai_tts*U_erai_tts + V_erai_tts*V_erai_tts)
Ws["attr"] = qcutils.MakeAttributeDictionary(long_name="Wind speed",units="m/s")
Wd = {}
Wd["data"] = float(270) - numpy.arctan2(V_erai_tts,U_erai_tts)*float(180)/numpy.pi
idx = numpy.where(Wd["data"]>360)[0]
if len(idx)>0: Wd["data"][idx] = Wd["data"][idx] - float(360)
Wd["attr"] = qcutils.MakeAttributeDictionary(long_name="Wind direction",units="deg")
return Ws,Wd
In [6]:
def get_Precip(erai_file,times_3hr,times_tts,fire):
# local pointers
flati = fire["lat_index"]
floni = fire["lon_index"]
dt_erai_utc_cor = times_3hr["dt_erai_utc_cor"]
erai_time_3hr = times_3hr["erai_time_3hr"]
dt_erai_utc_tts = times_tts["dt_erai_utc_tts"]
erai_time_tts = times_tts["erai_time_tts"]
hour_utc = times_3hr["hour_utc"]
# Spread the 3 hourly accumulated precipitation to the tower time step
# NOTE: ERA-I variables are dimensioned [time,latitude,longitude]
Precip_3d = erai_file.variables["tp"][:,:,:]
Precip_accum = Precip_3d[:,flati,floni]
Precip_erai_3hr = numpy.ediff1d(Precip_accum,to_begin=0)
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Precip_erai_3hr[idx] = Precip_accum[idx]
Precip_erai_3hr = Precip_erai_3hr*float(1000)
Precip = {}
Precip["data"] = numpy.zeros(len(dt_erai_utc_tts))
idx = qcutils.FindIndicesOfBInA(dt_erai_utc_cor,dt_erai_utc_tts)
Precip["data"][idx] = Precip_erai_3hr
Precip["attr"] = qcutils.MakeAttributeDictionary(long_name="Precipitation",units="mm")
return Precip
In [7]:
def nc_write_series(out_file,label,data):
ncVar = out_file.createVariable(label,"d",("time",))
ncVar[:] = data["data"]
for item in data["attr"]:
setattr(ncVar,item,data["attr"][item])
In [8]:
# get the filename
file_path = "/home/peter/hackfest/"
year = 2014
csv_filename = file_path+"hotspot_"+str(year)+"_short.csv"
erai_name = file_path+"ERAI_"+str(year)+".nc"
# sniff the file to see if we can find the dialect
csv_file = open(csv_filename,'rb')
dialect = csv.Sniffer().sniff(csv_file.readline(), [' ',',','\t'])
# rewind file
csv_file.seek(0)
csv_reader = csv.reader(csv_file,dialect)
# get the header line
header = csv_reader.next()
# close the file
csv_file.close()
In [9]:
# get a list of columns to be read
csv_list = ["load_dt","latitude","longitude"]
col_list = [header.index(item) for item in csv_list]
# read the csv file using numpy's genfromtxt
skip = 1
# define the missing values and the value with which to fill them
missing_values = {}
filling_values = {}
for item in col_list:
missing_values[item] = ["NA","N/A","NAN","#NAME?","#VALUE!","#DIV/0!","#REF!"]
filling_values[item] = c.missing_value
# read the CSV file
data = numpy.genfromtxt(csv_filename,delimiter=dialect.delimiter,skip_header=skip,
names=header,usecols=col_list,missing_values=missing_values,
filling_values=filling_values,dtype=None)
In [10]:
erai_timestep = 180
# dictionary for ERAI times at 3 hour time step
times_3hr = {}
# open the ERAI netCDF file
erai_file = netCDF4.Dataset(erai_name)
# get the latitude and longitude data
latitude = erai_file.variables["latitude"][:]
longitude = erai_file.variables["longitude"][:]
# get the latitude and longitude resolutions
lat_resolution = abs(latitude[-1]-latitude[0])/(len(latitude)-1)
lon_resolution = abs(longitude[-1]-longitude[0])/(len(longitude)-1)
# get the time and convert to Python datetime object
erai_time = erai_file.variables["time"][:]
time_units = getattr(erai_file.variables["time"],"units")
dt_erai = netCDF4.num2date(erai_time,time_units)
times_3hr["hour_utc"] = numpy.array([dt.hour for dt in dt_erai])
# get the datetime in the middle of the accumulation period
erai_offset = datetime.timedelta(minutes=float(erai_timestep)/2)
dt_erai_cor = [x - erai_offset for x in dt_erai]
# get a series of time, corrected for the offset
# NOTE: netCDF4.date2num doesn't handle timezone-aware datetimes
times_3hr["erai_time_3hr"] = netCDF4.date2num(dt_erai_cor,time_units)
# make utc_dt timezone aware so we can generate local times later
times_3hr["dt_erai_utc_cor"] = [x.replace(tzinfo=pytz.utc) for x in dt_erai_cor]
In [11]:
# loop over latitude and longitude of fires
nrecs = len(data["latitude"])
lat_indices = []
lon_indices = []
n = 0
for i in range(100):
# dictionaries for erai and fire information and times
erai = {}; fire = {}; times_tts = {}
# put useful stuff in dictionaries
fire["latitude"] = data["latitude"][i]
fire["longitude"] = data["longitude"][i]
# index of fire in latitude dimension
fire["lat_index"] = int(((latitude[0]-fire["latitude"])/lat_resolution)+0.5)
erai["latitude"] = latitude[fire["lat_index"]]
# index of the fire in longitude dimension
if fire["longitude"]<0: fire["longitude"] = float(360) + fire["longitude"]
fire["lon_index"] = int(((fire["longitude"]-longitude[0])/lon_resolution)+0.5)
erai["longitude"] = longitude[fire["lon_index"]]
# update the latitude and longitude index lists, these are used to reject fire
# locations that resolveto the same ERAI pixels
if i==0:
# add the indices if it is the first time through
lat_indices.append(fire["lat_index"])
lon_indices.append(fire["lon_index"])
# check to see if we have had this combination of indices before
if fire["lat_index"] in lat_indices and fire["lon_index"] in lon_indices:
# if we have, then move on to the next fire location
continue
else:
# if we haven't, then add these indices to the list of processed pixels
lat_indices.append(fire["lat_index"])
lon_indices.append(fire["lon_index"])
# ... and then go and do the processing
n += 1
print "For fire no. ",n
print " Fire coordinates: ",fire["latitude"],fire["longitude"]
print " ERAI coordinates: ",erai["latitude"],erai["longitude"]
# now we get the datetime series at hourly time step
tdts = datetime.timedelta(minutes=60)
# get the start and end datetimes rounded to the nearest time steps
# that lie between the first and last times
start_date = qcutils.rounddttots(times_3hr["dt_erai_utc_cor"][0],ts=60)
if start_date<times_3hr["dt_erai_utc_cor"][0]: start_date = start_date+tdts
end_date = qcutils.rounddttots(times_3hr["dt_erai_utc_cor"][-1],ts=60)
if end_date>times_3hr["dt_erai_utc_cor"][-1]: end_date = end_date-tdts
print " Got data from ",start_date," UTC to ",end_date," UTC"
# UTC datetime series at the tower time step
times_tts["dt_erai_utc_tts"] = [x for x in qcutils.perdelta(start_date,end_date,tdts)]
# UTC netCDF time series at tower time step for interpolation
tmp = [x.replace(tzinfo=None) for x in times_tts["dt_erai_utc_tts"]]
times_tts["erai_time_tts"] = netCDF4.date2num(tmp,time_units)
# get the incoming short wave radiation
Fsd = get_Fsd(erai_file,times_3hr,times_tts,fire,erai)
# get the air temperature
Ta = get_Ta(erai_file,times_3hr,times_tts,fire)
# get the relative humidity
RH = get_RH(erai_file,times_3hr,times_tts,fire)
# get the wind speed
Ws,Wd = get_Ws_and_Wd(erai_file,times_3hr,times_tts,fire)
# get the precipitation
Precip = get_Precip(erai_file,times_3hr,times_tts,fire)
# construct the output file name
out_erai_name = file_path+str(year)+"_"+str(fire["lat_index"])+"_"+str(fire["lon_index"])+".nc"
out_file = netCDF4.Dataset(out_erai_name,"w")
# create the time dimension and the time variable
out_file.createDimension("time",len(times_tts["erai_time_tts"]))
ncVar = out_file.createVariable("time","d",("time",))
ncVar[:] = times_tts["erai_time_tts"]
setattr(ncVar,"long_name","time")
setattr(ncVar,"standard_name","time")
setattr(ncVar,"units",time_units)
setattr(ncVar,"calendar","gregorian")
# write out the incoming short wave radiation
nc_write_series(out_file,"Fsd",Fsd)
# write out the air temperature
nc_write_series(out_file,"Ta",Ta)
# write out the relative humidity
nc_write_series(out_file,"RH",RH)
# write out the wind speed
nc_write_series(out_file,"Ws",Ws)
# write out the wind direction
nc_write_series(out_file,"Wd",Wd)
# write out the precipitation
nc_write_series(out_file,"Precip",Precip)
# close the netCDF file
out_file.close()
In [20]:
times_tts = {}
# UTC datetime series at the tower time step
times_tts["dt_erai_utc_tts"] = [x for x in qcutils.perdelta(start_date,end_date,tdts)]
# UTC netCDF time series at tower time step for interpolation
tmp = [x.replace(tzinfo=None) for x in times_tts["dt_erai_utc_tts"]]
times_tts["erai_time_tts"] = netCDF4.date2num(tmp,time_units)
# local pointers
flati = fire["lat_index"]
floni = fire["lon_index"]
dt_erai_utc_cor = times_3hr["dt_erai_utc_cor"]
erai_time_3hr = times_3hr["erai_time_3hr"]
dt_erai_utc_tts = times_tts["dt_erai_utc_tts"]
erai_time_tts = times_tts["erai_time_tts"]
hour_utc = times_3hr["hour_utc"]
Precip_3d = erai_file.variables["tp"][:,:,:]
Precip_accum = Precip_3d[:,flati,floni]
Precip_erai_3hr = numpy.ediff1d(Precip_accum,to_begin=0)
idx = numpy.where((hour_utc==3)|(hour_utc==15))[0]
Precip_erai_3hr[idx] = Precip_accum[idx]
Precip_erai_3hr = Precip_erai_3hr*float(1000)
Precip = {}
Precip["data"] = numpy.zeros(len(dt_erai_utc_tts))
idx = qcutils.FindIndicesOfBInA(dt_erai_utc_cor,dt_erai_utc_tts)
Precip["data"][idx] = Precip_erai_3hr
In [21]:
fig=plt.figure()
plt.plot(Precip["data"])
plt.show()
In [24]:
idx = qcutils.FindIndicesOfBInA(dt_erai_utc_tts,dt_erai_utc_cor)
print len(idx)
In [25]:
print erai_time_tts[0:24]
print erai_time_cor[0:6]
In [ ]: