Algorithm Description
This algorithm computes the monthly means of the aerosols and non-opaque (at 333m scale) cloud optical depth retrieved by the SODA 5km products on a 2.5 by 2.5 degrees orthogonal grid. It produces 3 synthesis per month : day orbits only, night orbits only and both. The day/night definition is based on the CALIPSO one.
SODA Level 3 Builder
This is a very basic averaging scheme : an arithmetic mean of observations is computed in each output grid cell. Here is what it is like in pseudo-code :
def get_averaged_dataset(dataset_name, v_files): #---------------------------------------------------------------------# # projects the pixels of the given files in a plane square grid, computes# # the arithmetic mean of the given variable in each cell and return it# #---------------------------------------------------------------------# # set accumulation buffers of size [sz_y, sz_x]. Initialized to zeros outgrid_acc = allocate_zeros_array("float", [sz_y, sz_x]) outgrid_num = allocate_zeros_array("integer", [sz_y, sz_x]) # - accumulate all pixels in the grid cells for f in v_files : # read vectors of(lat,lon)pixel positions v_lat = read_data(f, "Latitude") v_lon = read_data(f, "Longitude") # number of pixels in file "f" npix = v_lat.size() # read vector of retrieved data v_data = read_data(f, dataset_name) fill_value = get_fill_value(f, dataset_name) for ipix in range(npix): # [i, j] are the indices of the grid cell where falls the pixel # (lat, lon) i, j = get_igrid(lat, lon) data = v_data[ipix] if data != fill_value : # accumulate in grid cell x and increment number of values outgrid_acc[i, j] += data outgrid_num[i, j] += 1 # end if # end for # end for # - average the accumulated data outdata = allocate_array("float", [sz_y, sz_x]) for j in range(sz_y): for i in range(sz_x): if outgrid_num > 0 : # compute mean outdata[i, j] = outgrid_acc[i, j] / outgrid_num[i, j] else : # set to fill value outdata[i, j] = fill_value # end if # end for # end for return outdata # end def get_averaged_dataset def build_synthesis(year, month, dn_flag): #-------------------------------------------------------------------------# # build the monthly synthesis of the selected type of orbits # # dn_mode is 'D', 'N' or 'A' and stands for day only, night only or all # #-------------------------------------------------------------------------# v_files = get_infiles(year, month, dn_flag) outfile_name = get_outfilename(year, month, dn_flag) outfile = open(outfile_name, "w")# output HDF file for dataset_name in ["Optical_Depth_532_Cloud", "Optical_Depth_532_Aerosol", "Optical_Depth_1064_Cloud", "Optical_Depth_1064_Aerosol"] : outdata = get_averaged_dataset(dataset_name, v_files) write(outfile, dataset_name, outdata) # end for close(outfile) # end def build_synthesis # --- main --- # for dn_flag in ['D', 'N', 'A'] : build_synthesis(year, month, dn_flag) # end for