In [1]:
%run basics
%matplotlib
In [149]:
def perdelta(start, end, delta):
"""
Yields an iterator of datetime objects from start to end with time step delta.
"""
curr = start
while curr <= end:
yield curr
curr += delta
In [54]:
# fire events
fire_events = {}
fire_events["2001-08-06"] = {"intensity":3550,"uncertainty":640,"units":"kW/m"}
fire_events["2002-08-18"] = {"intensity":3000,"uncertainty":None,"units":"kW/m"}
fire_events["2003-08-29"] = {"intensity":3200,"uncertainty":None,"units":"kW/m"}
fire_events["2004-08-18"] = {"intensity":3610,"uncertainty":745,"units":"kW/m"}
fire_events["2005-06-08"] = {"intensity":1399,"uncertainty":837,"units":"kW/m"}
fire_events["2006-05-26"] = {"intensity":1047,"uncertainty":232,"units":"kW/m"}
fire_events["2009-08-06"] = {"intensity":2000,"uncertainty":None,"units":"kW/m"}
fire_events["2011-05-18"] = {"intensity":964,"uncertainty":345,"units":"kW/m"}
fire_events["2012-08-19"] = {"intensity":2500,"uncertainty":None,"units":"kW/m"}
fire_events["2013-06-07"] = {"intensity":1000,"uncertainty":None,"units":"kW/m"}
fire_events["2014-09-12"] = {"intensity":2200,"uncertainty":None,"units":"kW/m"}
fire_events["2015-04-29"] = {"intensity":500,"uncertainty":None,"units":"kW/m"}
fire_dates = [dateutil.parser.parse(fd) for fd in fire_events.keys()]
In [49]:
# DINGO version of interpolated and smoothed EVI data
dingo_name = "../../Sites/HowardSprings/Data/MODIS/HowardSprings_MODIS.nc"
dingo_file = netCDF4.Dataset(dingo_name,"r")
dingo_time = dingo_file.variables["time"][:]
dingo_time_units = getattr(dingo_file.variables["time"],"units")
dingo_dt = netCDF4.num2date(dingo_time,dingo_time_units)
dingo_evi_interp = dingo_file.variables["EVI_MODIS_interp"][:]
dingo_evi_smooth = dingo_file.variables["EVI_MODIS_smooth"][:]
In [7]:
dap_url = "http://www.auscover.org.au/thredds/dodsC/"
dap_folder = "auscover/lpdaac-aggregates/c5/v2-nc4/aust/MOD13Q1.005/"
dap_name = "MOD13Q1.aggregated.aust.005.enhanced_vegetation_index.ncml"
dap_file = dap_url+dap_folder+dap_name
print dap_file
In [8]:
nc_file = netCDF4.Dataset(dap_file,"r")
In [9]:
lat_resolution = getattr(nc_file,"geospatial_lat_resolution")
lon_resolution = getattr(nc_file,"geospatial_lon_resolution")
print lat_resolution,lon_resolution
In [155]:
site_name = "HowardSprings"
site_latitude = -12.495200
site_longitude = 131.150050
ts = 30
cutout_width = 3
cutout_height = 3
delta_width = cutout_width*lat_resolution/2
delta_height = cutout_height*lon_resolution/2
In [136]:
lat_bound_lower = site_latitude - delta_width
lat_bound_upper = site_latitude + delta_width
lon_bound_lower = site_longitude - delta_height
lon_bound_upper = site_longitude + delta_height
In [137]:
lat = nc_file.variables["latitude"][:]
lon = nc_file.variables["longitude"][:]
In [150]:
modis_time = nc_file.variables["time"][:]
modis_time_units = getattr(nc_file.variables["time"],"units")
modis_dt = netCDF4.num2date(modis_time,modis_time_units)
print modis_dt[0],modis_dt[-1]
In [139]:
lat_index = numpy.where((lat>=lat_bound_lower)&(lat<=lat_bound_upper))[0]
lon_index = numpy.where((lon>=lon_bound_lower)&(lon<=lon_bound_upper))[0]
# get the EVI
evi = nc_file.variables["evi"][:,lat_index,lon_index]
# get the quality flags
quality = nc_file.variables["quality"][:,lat_index,lon_index]
# get the "typical mask"
#typical_mask = nc_file.variables["typical_mask"][:,lat_index,lon_index]
In [188]:
ok_mask = numpy.ones_like(evi)
ok_list = [2048,2049,2052,2053,2112,2113,2116,2117,2560,2561,2564,2565,2624,2625,2628,2629]
for item in ok_list:
index = numpy.ma.where(quality==item)[0]
ok_mask[index] = 0
evi_masked = numpy.ma.masked_where(ok_mask!=0,evi)
evi_masked_median = numpy.ma.median(evi_masked.reshape(evi_masked.shape[0],-1),axis=1)
In [152]:
fig=plt.figure()
for i in range(evi.shape[1]):
for j in range(evi.shape[2]):
plt.plot(modis_dt,evi[:,i,j],'b.')
plt.plot(dingo_dt,dingo_evi_interp[:,0,0],'r-')
plt.plot(dingo_dt,dingo_evi_smooth[:,0,0],'g-')
for item in fire_dates: plt.axvline(item)
plt.show()
In [190]:
start = modis_dt[0]
end = modis_dt[-1]
modis_dt_interp = [result for result in perdelta(start,end,datetime.timedelta(minutes=ts))]
modis_time_interp = netCDF4.date2num(modis_dt_interp,modis_time_units)
modis_time_masked = numpy.ma.masked_where(numpy.ma.getmaskarray(evi_masked_median)==True,modis_time)
modis_time_comp = numpy.ma.compressed(modis_time_masked)
evi_masked_median_comp = numpy.ma.compressed(evi_masked_median)
x_org = modis_time_comp
y_org = evi_masked_median_comp
In [ ]:
# linear interpolation followed by Savitsky-Golay filter
f = scipy.interpolate.interp1d(x_org,y_org,bounds_error=False)
x_int = modis_time_interp
evi_interp = f(x_int)
evi_interp_smooth= scipy.signal.savgol_filter(evi_interp, 10001, 4, 0, 1)
In [212]:
# smoothed spline interpolation
tck = scipy.interpolate.splrep(x_org, y_org, s=0.03)
evi_interp_smooth = scipy.interpolate.splev(x_int, tck, der=0)
In [192]:
s = scipy.interpolate.InterpolatedUnivariateSpline(x_org,y_org,k=3)
evi_interp_smooth = s(x_int)
In [213]:
fig=plt.figure()
for i in range(evi_masked.shape[1]):
for j in range(evi_masked.shape[2]):
plt.plot(modis_dt,evi_masked[:,i,j],'b.')
plt.plot(dt,evi_masked_median,'ro')
plt.plot(modis_dt_interp,evi_interp,'y-')
plt.plot(modis_dt_interp,evi_interp_smooth,'r-')
#plt.plot(dingo_dt,dingo_evi_interp[:,0,0],'r-')
plt.plot(dingo_dt,dingo_evi_smooth[:,0,0],'g-')
for item in fire_dates: plt.axvline(item)
plt.show()
In [85]:
flag_masks = getattr(nc_file.variables["quality"],"flag_masks")
flag_values = getattr(nc_file.variables["quality"],"flag_values")
flag_meanings = getattr(nc_file.variables["quality"],"flag_meanings")
comment = getattr(nc_file.variables["quality"],"comment")
In [82]:
flag_meanings_list = flag_meanings.split(" ")
print len(flag_masks),len(flag_values),len(flag_meanings_list)
In [84]:
for i in range(len(flag_values)):
print flag_masks[i],",",flag_values[i]," - ",flag_meanings_list[i]
In [86]:
print comment
In [108]:
bin(2048)[2:]
Out[108]:
In [112]:
int('100000000100',2)
Out[112]:
In [ ]: