Skip to content
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
315 changes: 175 additions & 140 deletions L1C_Aux/collocate_hirs_avhrr.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,13 @@ def create_filepaths(args):

return hirs_dir,fiduceo_dir,gac_dir,l1c_dir

def closest(num,arr):
curr = arr[0]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a docstring to describe what the function does, inputs and outputs.

for val in arr:
if abs (num - val) < abs (num - curr):
curr = val
return curr

if __name__ in '__main__':

#Pass in the date to process
Expand All @@ -58,143 +65,171 @@ def create_filepaths(args):
string_date_hirs = in_date[0:4]+'/'+in_date[4:6]+'/'+in_date[6:8]+'/'

for f in os.listdir(hirs_dir+string_date_hirs):

fiduceo_filename = rf.construct_fiduceo(fiduceo_dir,string_date_hirs,f)
time_hirs = rf.read_time(hirs_dir+string_date_hirs+f)

plt.figure()
plt.xlabel('Lon')
plt.ylabel('Lat')

fileobj = -1

for date in [prev,in_date,next]:
string_date = date[0:4]+'/'+date[4:6]+'/'+date[6:8]+'/'



for g in os.listdir(gac_dir+string_date):
time_gac = rf.read_time(gac_dir+string_date+g)
#Identify files where the HIRS/AVHRR times overlap
if np.min(time_gac) > np.max(time_hirs):
continue
elif np.max(time_gac) < np.min(time_hirs):
continue
else:
print 'overlap'
#Identify overlap times for GAC and HIRS
gac_min,gac_max,hirs_min,hirs_max,data_check \
= td.time_mask(time_hirs,time_gac)
if data_check == True:
#Read in for GAC and HIRS: lat,lon and l2 flags
gac_lat,gac_lon,gac_flags \
= rf.refine_geo(gac_dir+string_date+g,gac_min,gac_max)
hirs_lat,hirs_lon,hirs_flags \
= rf.refine_geo(hirs_dir+string_date_hirs+f,\
hirs_min,hirs_max)

#Read in GAC and HIRS BT
hirs_bt,hirs_noise,hirs_obs,hirs_ind,hirs_str,hirs_com\
= rf.read_hirs_obs(fiduceo_dir+string_date_hirs \
+fiduceo_filename,hirs_min,hirs_max)
l1c_filename = rf.construct_l1c(l1c_dir,string_date,g,gac_dir)
print l1c_filename
gac_bt,avhrr_obs,avhrr_noise \
= rf.read_avhrr_obs(l1c_dir+string_date+l1c_filename,\
gac_min,gac_max,args.hirs_sensor)

if os.path.isfile(avhrr_sim_dir+string_date+g):
gac_dy = rf.read_dy(avhrr_sim_dir+string_date+g,gac_min,\
gac_max)
else:
gac_dy=np.zeros([1])

#Read in the PClear arrays
gac_prob = rf.read_pclear(gac_dir+string_date+g,gac_min,gac_max)

hirs_prob = rf.read_pclear(hirs_dir+string_date_hirs+f,hirs_min,\
hirs_max)

#### Extra probabilities from different runs - move from
#### core code??
# hirs_prob_lim = rf.read_pclear(comp_dir+string_date_hirs+f,hirs_min,\
# hirs_max)
# hirs_prob_five = rf.read_pclear(all_dir+string_date_hirs+f,hirs_min,\
# hirs_max)

#Read in the dY arrays
dy = rf.read_ffm(hirs_dir+string_date_hirs+f,hirs_min,hirs_max)

#Co-locate the GAC and HIRS arrays
lat_centre,lon_centre,flag_centre,gac_as_hirs_min,\
gac_as_hirs_mean,gac_bt_mean,gac_bt_cloud_std,gac_bt_all_std,\
cloud_frac,gac_as_hirs_dy,gac_as_hirs_n,gac_as_hirs_noise, \
gac_as_hirs_landmask \
= col.collocate_gac_hirs(time_gac,time_hirs,gac_min,gac_max,\
hirs_min,hirs_max,gac_lat,gac_lon,\
gac_flags,hirs_flags,gac_prob,\
avhrr_obs,avhrr_noise,args.hirs_sensor,gac_dy)


#Calculate the cloud height
# try:
# cloud_height = np.zeros([lat_centre.shape[0],lat_centre.shape[1]])
# mask = (dy[7,:,:] <= 0.)
# cloud_height[mask] = abs(dy[7,:,:][mask])*(1./6.)
# try:
# cloud_height[dy[7,:,:].mask] = -np.nan
# except:
# pass
# except:
# pass

#Write out data
outfile = wf.get_output_filename(args,fiduceo_filename)

if fileobj == -1:
fileobj,var_rtm,var_lat,var_lon,var_flag,\
var_hlat,var_hlon,var_hflag,var_a_obs,var_a_noise, \
var_cloud_frac,var_cloud_height,var_cloud_std,\
var_obs_std,var_n,var_likelihood,var_rtm_land, \
var_likelihood_land,var_landmask \
= wf.create_output(hirs_obs,avhrr_obs,\
time_hirs,outfile)


fileobj = wf.write_data(hirs_obs,dy,lat_centre,lon_centre,\
flag_centre,hirs_lat,hirs_lon,\
hirs_flags,gac_bt_mean,gac_as_hirs_noise,\
fileobj,hirs_min,hirs_max,cloud_frac,\
gac_bt_cloud_std,\
gac_bt_all_std,gac_as_hirs_n,hirs_prob,\
gac_as_hirs_landmask,var_rtm,var_lat,var_lon,\
var_flag,var_hlat,var_hlon,\
var_hflag,var_a_obs,var_a_noise,\
var_cloud_frac,var_cloud_height,\
var_cloud_std,var_obs_std,var_n,var_likelihood,\
var_rtm_land,var_likelihood_land,var_landmask)


#Make some comparison plots

mask_gac,mask_hirs \
= pr.plot_geoloc(lat_centre,lon_centre,flag_centre,hirs_lat,\
hirs_lon,hirs_flags)
# pr.plot_cloud_char(cloud_frac,gac_bt_std,gac_bt_mean,ch8)


# pr.plot_sim_comp(gac_bt_mean,gac_as_hirs_dy,hirs_bt,ch8)
# cm.explore_bias(hirs_bt,gac_bt_mean,hirs_noise,gac_bt_all_std,\
# gac_as_hirs_n)

# cm.calc_stats(hirs_bt,hirs_noise,gac_bt_mean,gac_as_hirs_mean,\
# hirs_prob,gac_bt_std,cloud_frac,gac_bt_mean_ext,\
# gac_as_hirs_mean_ext,gac_bt_std_ext,cloud_frac_ext)

# pr.plot_cloudmasks(gac_as_hirs_min,hirs_prob,mask_gac,mask_hirs,\
# ch7,ch8,ch10,ch11,ch12,hirs_prob_lim,\
# hirs_prob_five,gac_as_hirs_mean,cloud_frac)


#Close output file
fileobj.close()
try:
fiduceo_filename = rf.construct_fiduceo(fiduceo_dir,string_date_hirs,f)

time_hirs = rf.read_time(hirs_dir+string_date_hirs+f)

plt.figure()
plt.xlabel('Lon')
plt.ylabel('Lat')

fileobj = -1

for date in [prev,in_date,next]:
string_date = date[0:4]+'/'+date[4:6]+'/'+date[6:8]+'/'



for g in os.listdir(gac_dir+string_date):
time_gac = rf.read_time(gac_dir+string_date+g)
#Identify files where the HIRS/AVHRR times overlap
if np.min(time_gac) > np.max(time_hirs):
continue
elif np.max(time_gac) < np.min(time_hirs):
continue
else:
print 'overlap'
#Identify overlap times for GAC and HIRS
gac_min,gac_max,hirs_min,hirs_max,data_check \
= td.time_mask(time_hirs,time_gac)
if data_check == True:
#Read in for GAC and HIRS: lat,lon and l2 flags

gac_lat,gac_lon,gac_flags \
= rf.refine_geo(gac_dir+string_date+g,gac_min,gac_max)
hirs_lat,hirs_lon,hirs_flags \
= rf.refine_geo(hirs_dir+string_date_hirs+f,\
hirs_min,hirs_max)

#Read in GAC and HIRS BT
hirs_bt,hirs_noise,hirs_obs,hirs_ind,hirs_str,hirs_com\
= rf.read_hirs_obs(fiduceo_dir+string_date_hirs \
+fiduceo_filename,hirs_min,hirs_max)

l1c_filename = rf.construct_l1c(l1c_dir,string_date,g,gac_dir)

# Code now attempts to find file with matching hour start time
try:
os.path.isfile(l1c_dir+string_date+l1c_filename)
#l1c_filename = rf.construct_l1c(l1c_dir,string_date,g,gac_dir)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This commented line should be removed as it is copied above.

gac_bt,avhrr_obs,avhrr_noise \
= rf.read_avhrr_obs(l1c_dir+string_date+l1c_filename,\
gac_min,gac_max,args.hirs_sensor)
if os.path.isfile(avhrr_sim_dir+string_date+g):
gac_dy = rf.read_dy(avhrr_sim_dir+string_date+g,gac_min,\
gac_max)
else:
gac_dy=np.zeros([1])
gac_prob = rf.read_pclear(gac_dir+string_date+g,gac_min,gac_max)

# If hour-matched start time is missed, it now finds the closest in time and replaces
# l1c_filename with the a new version which contains the time of the nearest file
except:

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't like this change because if problem is that the file is on the next day, rather than the current day due to a timing issue, your code will now read the previous file (closest in time) which doesn't necessarily overlap. It should run in the sense that no overlap should be found but it is an unnecessary processing overhead. You could also perform a check for a maximum time delta on an accepted filename which would improve things.

print "Mismatched file found in l2p and l1c for: ",g
old_g_time = int(g[8:14])

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would be inclined to move these ~five lines of code to a separate function for clarity.

new_g_list = os.listdir(l1c_dir+string_date)
newlist = np.array([int(x[8:14]) for x in new_g_list])
new_g_time = closest(old_g_time,newlist)
new_g = str(g[:8])+str(new_g_time)+str(g[14:])
new_l1c_filename = str(l1c_filename[:8])+str(new_g_time)+str(l1c_filename[14:])
print "New l1c_filename: ",new_l1c_filename
gac_bt,avhrr_obs,avhrr_noise \
= rf.read_avhrr_obs(l1c_dir+string_date+new_l1c_filename,\
gac_min,gac_max,args.hirs_sensor)
if os.path.isfile(avhrr_sim_dir+string_date+g):
gac_dy = rf.read_dy(avhrr_sim_dir+string_date+g,gac_min,\
gac_max)
else:
gac_dy=np.zeros([1])
gac_prob = rf.read_pclear(gac_dir+string_date+g,gac_min,gac_max)


hirs_prob = rf.read_pclear(hirs_dir+string_date_hirs+f,hirs_min,\
hirs_max)

#### Extra probabilities from different runs - move from
#### core code??
# hirs_prob_lim = rf.read_pclear(comp_dir+string_date_hirs+f,hirs_min,\
# hirs_max)
# hirs_prob_five = rf.read_pclear(all_dir+string_date_hirs+f,hirs_min,\
# hirs_max)

#Read in the dY arrays
dy = rf.read_ffm(hirs_dir+string_date_hirs+f,hirs_min,hirs_max)

#Co-locate the GAC and HIRS arrays
lat_centre,lon_centre,flag_centre,gac_as_hirs_min,\
gac_as_hirs_mean,gac_bt_mean,gac_bt_cloud_std,gac_bt_all_std,\
cloud_frac,gac_as_hirs_dy,gac_as_hirs_n,gac_as_hirs_noise, \
gac_as_hirs_landmask \
= col.collocate_gac_hirs(time_gac,time_hirs,gac_min,gac_max,\
hirs_min,hirs_max,gac_lat,gac_lon,\
gac_flags,hirs_flags,gac_prob,\
avhrr_obs,avhrr_noise,args.hirs_sensor,gac_dy)


#Calculate the cloud height
# try:
# cloud_height = np.zeros([lat_centre.shape[0],lat_centre.shape[1]])
# mask = (dy[7,:,:] <= 0.)
# cloud_height[mask] = abs(dy[7,:,:][mask])*(1./6.)
# try:
# cloud_height[dy[7,:,:].mask] = -np.nan
# except:
# pass
# except:
# pass

#Write out data
outfile = wf.get_output_filename(args,fiduceo_filename)

if fileobj == -1:
fileobj,var_rtm,var_lat,var_lon,var_flag,\
var_hlat,var_hlon,var_hflag,var_a_obs,var_a_noise, \
var_cloud_frac,var_cloud_height,var_cloud_std,\
var_obs_std,var_n,var_likelihood,var_rtm_land, \
var_likelihood_land,var_landmask \
= wf.create_output(hirs_obs,avhrr_obs,\
time_hirs,outfile)


fileobj = wf.write_data(hirs_obs,dy,lat_centre,lon_centre,\
flag_centre,hirs_lat,hirs_lon,\
hirs_flags,gac_bt_mean,gac_as_hirs_noise,\
fileobj,hirs_min,hirs_max,cloud_frac,\
gac_bt_cloud_std,\
gac_bt_all_std,gac_as_hirs_n,hirs_prob,\
gac_as_hirs_landmask,var_rtm,var_lat,var_lon,\
var_flag,var_hlat,var_hlon,\
var_hflag,var_a_obs,var_a_noise,\
var_cloud_frac,var_cloud_height,\
var_cloud_std,var_obs_std,var_n,var_likelihood,\
var_rtm_land,var_likelihood_land,var_landmask)


#Make some comparison plots

mask_gac,mask_hirs \
= pr.plot_geoloc(lat_centre,lon_centre,flag_centre,hirs_lat,\
hirs_lon,hirs_flags)
# pr.plot_cloud_char(cloud_frac,gac_bt_std,gac_bt_mean,ch8)


# pr.plot_sim_comp(gac_bt_mean,gac_as_hirs_dy,hirs_bt,ch8)
# cm.explore_bias(hirs_bt,gac_bt_mean,hirs_noise,gac_bt_all_std,\
# gac_as_hirs_n)

# cm.calc_stats(hirs_bt,hirs_noise,gac_bt_mean,gac_as_hirs_mean,\
# hirs_prob,gac_bt_std,cloud_frac,gac_bt_mean_ext,\
# gac_as_hirs_mean_ext,gac_bt_std_ext,cloud_frac_ext)

# pr.plot_cloudmasks(gac_as_hirs_min,hirs_prob,mask_gac,mask_hirs,\
# ch7,ch8,ch10,ch11,ch12,hirs_prob_lim,\
# hirs_prob_five,gac_as_hirs_mean,cloud_frac)


#Close output file
fileobj.close()
except:
print "Fiduceo file match failed - ",f