Process a matlab .mat file generated by proc_lisst.m

proc_lisst.py

# -*- coding: utf-8 -*-
"""
:DESCRIPTION:code to process LISST VC data

:REQUIRES:

:TODO:

:AUTHOR: John Franco Saraceno
:ORGANIZATION: U.S. Geological Survey, United States Department of Interior
:CONTACT: saraceno@usgs.gov
:VERSION: 0.1
Mon Aug 08 18:04:29 2016
"""
# =============================================================================
# IMPORT STATEMENTS
# =============================================================================
import datetime as datetime
import numpy as np
import scipy.io as sio
import pandas as pd

# =============================================================================
# METHODS
# =============================================================================


def is_empty(any_structure):
    if any_structure:
        #        print('Structure is not empty.')
        return False
    else:
        #        print('Structure is empty.')
        return True


def is_leap_year(year):
    """ if year is a leap year return True
        else return False """
    if year % 100 == 0:
        return year % 400 == 0
    return year % 4 == 0


def doy(Y, M, D):
    """ given year, month, day return day of year
        Astronomical Algorithms, Jean Meeus, 2d ed, 1998, chap 7 """
    if is_leap_year(Y):
        K = 1
    else:
        K = 2
    N = int((275 * M) / 9.0) - K * int((M + 9) / 12.0) + D - 30
    return N


def ymd(Y, N):
    """ given year = Y and day of year = N, return year, month, day
        Astronomical Algorithms, Jean Meeus, 2d ed, 1998, chap 7 """
    if is_leap_year(Y):
        K = 1
    else:
        K = 2
    M = int((9 * (K + N)) / 275.0 + 0.98)
    if N < 32:
        M = 1
    D = N - int((275 * M) / 9.0) + K * int((M + 9) / 12.0) + 30
    return Y, M, D


def datenumfromdata(data3940, year=2012):
    # PORTED TO PYTHON FROM SEQUOIA MATLAB FUNC OF SAME NAME
    #   compute days, hours, minutes, seconds from list date/time cols
    years = year*np.ones((len(data3940), 1))  # create vec of the initial year
    days = np.floor(data3940[:, 1] / 100)
    hours = data3940[:, 1] - (100 * days)
    minutes = np.floor(data3940[:, 0] / 100)
    seconds = np.floor((data3940[:, 0] - (100*minutes))/100*60)

    NewYear = np.where(np.diff(days) < 0)
    # find negative differences in doy (indicate deployment over new year)
    if not is_empty(NewYear[0]):  # If we have a new year deployment...
        years[NewYear+1:] = year+1
        """...the year after new year is one higher than
        the inital year (we assume that the deployment
        doesn't span more than one new year)!"""
    realdatenum = []
    for i in np.arange(len(data3940)):
        ymd_tuple = ymd(int(years[i]), int(days[i]))
        realdatenum.append(datetime.datetime(ymd_tuple[0],
                                             ymd_tuple[1], ymd_tuple[2],
                                             int(hours[i]),
                                             int(minutes[i]),
                                             int(seconds[i])))
#    return
    return realdatenum
# =============================================================================
# MAIN METHOD AND TESTING AREA
# =============================================================================


def proc_data_structure(oct_struct):
    #    filenames = []
    #    vd = []
    #    ts = []
    diameters_dataframes = []
    vd_dataframes = []

    dims = (oct_struct[0, 0].dias[0])
    diam_cols = [u'Total Volume Concentration in µl/l',
                 u'Mean Size in µm',
                 u'Standard Deviation in µm',
                 u'transmission ',
                 u'D10 in µm',
                 u'D16 in µm',
                 u'D50 in µm',
                 u'D60 in µm',
                 u'D84 in µm',
                 u'D90 in µm',
                 u'D60/D10',
                 u'Surface area in cm2/l',
                 u'Silt Density',
                 u'Silt Volume']

    for i in np.arange(len(oct_struct)):
        if i is 3:
            year = 2012
        else:
            year = 2013
        print "processing file: {}".format(oct_struct[i, 0].filename)
    #    ts.append(datenumfromdata(oct_struct[i, 0].data[:,37:39], year=year))
    #    ts.append(datenumfromdata(oct_struct[i, 0].data[:,37:39], year=year))
        index = datenumfromdata(oct_struct[i, 0].data[:, 37:39], year=year)
    #    filenames.append(str(oct_struct[i, 0].filename))
    #    vd.append(oct_struct[i, 0].vd_corrected)
        diameters_dataframes.append(pd.DataFrame(oct_struct[i, 0].diamters,
                                    columns=diam_cols, index=index))
        vd_dataframes.append(pd.DataFrame(oct_struct[i, 0].vd_corrected,
                             columns=dims,
                             index=index))

    return diameters_dataframes, vd_dataframes


def main(filename, struct_name='output'):

    mat_contents = sio.loadmat(filename,
                               struct_as_record=False, squeeze_me=False)
    oct_struct = mat_contents[struct_name]
    diameters_dataframes, vd_dataframes = proc_data_structure(oct_struct)
    return diameters_dataframes, vd_dataframes

# i = 0
# year = 2013
# data3940 = oct_struct[i, 0].data[:, 37:39]
fname = 'cache_lisst_random_v7.mat'
diameters_dataframes, vd_dataframes = main(fname)