SODA L3

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