Skip to content
Snippets Groups Projects
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
manage_GSNdata.py 36.91 KiB
'''
Data manager for PermaSense data. Initially created for the paper "A decade of detailed observations (2008--2018) in steep 
bedrock permafrost at Matterhorn Hoernligrat (Zermatt, CH)"

Copyright (c) 2020, Swiss Federal Institute of Technology (ETH Zurich)
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

* Neither the name of the copyright holder nor the names of its
  contributors may be used to endorse or promote products derived from
  this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

@date: October 31, 2020
@author: Samuel Weber, Jan Beutel and Matthias Meyer

usage: manage_GSNdata.py [-h] [--gsn2raw GSN2RAW] [--load LOAD]
                              [--filter FILTER] [--clean CLEAN]
                              [--aggregationInterval AGGREGATIONINTERVAL]
                              [--export2csv EXPORT2CSV]
                              [--sanityPlot SANITYPLOT] [--gsn2img GSN2IMG]
                              [--nImg NIMG] [-p PATH] [-y1 YEARBEGIN]
                              [-m1 MONTHBEGIN] [-d1 DAYBEGIN] [-y2 YEAREND]
                              [-m2 MONTHEND] [-d2 DAYEND]

optional arguments:
  -h, --help            show this help message and exit
  --gsn2raw GSN2RAW     Get raw data from GSN and store locally (step1,
                        default: True)
  --load LOAD           Load locally stored data (step2, default: True)
  --filter FILTER       Filter data (according reference values, step3,
                        default: True)
  --clean CLEAN         Clean data manually (step4, default: True)
  --aggregationInterval AGGREGATIONINTERVAL
                        Aggregate data over timewindows of X minutes (step5,
                        default: 60)
  --export2csv EXPORT2CSV
                        Export data in yearly csv files (step6, default: True)
  --sanityPlot SANITYPLOT
                        Plot for sanity check (step7, default: True)
  --gsn2img GSN2IMG     Downlod images from GSN and convert NEF to JPEG
                        (default: True)
  --nImg NIMG           Number of images to download and convert (default: 10)
  -p PATH, --path PATH  Relative path for the data output (default: ./data)
  -y1 YEARBEGIN, --yearBegin YEARBEGIN
                        Year begin (default: 2008)
  -m1 MONTHBEGIN, --monthBegin MONTHBEGIN
                        Month begin (default: 1)
  -d1 DAYBEGIN, --dayBegin DAYBEGIN
                        Day begin (default: 1)
  -y2 YEAREND, --yearEnd YEAREND
                        Year end (default: today)
  -m2 MONTHEND, --monthEnd MONTHEND
                        Month end (default: today)
  -d2 DAYEND, --dayEnd DAYEND
                        Day end (default: today)
'''
def main():
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    ## SETUP PYTHON ENVIRONMENT
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    
    # Load functions
    import sys
    import socket
    import argparse
    import os as os
    import numpy as np
    import pandas as pd
    import datetime as dt
    import matplotlib as mpl

    # For proper plotting on server
    if socket.gethostname() == 'tik43x':
        mpl.use('Agg')
    import matplotlib.pyplot as plt

    # faster plotting
    import matplotlib.style as mplstyle
    mplstyle.use('fast')
    # remove matplotlib warning 
    from pandas.plotting import register_matplotlib_converters
    register_matplotlib_converters()
    mpl.rc('figure', max_open_warning = 0)
    
    #adjust display width for variables
    pd.set_option('display.max_rows', 500)
    pd.set_option('display.max_columns', 500)
    pd.set_option('display.width', 1000)

    
    # Load functions in relative path
    sys.path.append('~/permasense/')
    from permasense.GSNdata import get_GSNdata, save_data, load_data, filter_data, clean_data, get_GSNimg
    
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    ## ASSIGN ARGUMENTS TO VARIABLES
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    args = parseArguments()
    
    # Define path to data
    PATH_DATA = args.path
    
    # Define time range
    TBEG = dt.datetime(args.yearBegin,args.monthBegin,args.dayBegin)
    TEND = dt.datetime(args.yearEnd,args.monthEnd,args.dayEnd)
    
    # Get raw data from GSN (step1) and produce derived data product (step2 - step7)
    GSN2RAW = args.gsn2raw
    
    LOAD = args.load
    FILTER = args.filter
    CLEAN = args.clean
    AGGREGATION_INTERVAL = args.aggregationInterval
    EXPORT2CSV = args.export2csv
    SANITY_PLOT = args.sanityPlot
    if FILTER | CLEAN | EXPORT2CSV | SANITY_PLOT:
        LOAD = True
    
    # Get high-resolution images (.NEF) and convert to .JPG
    GSN2IMG = args.gsn2img
    N_IMG = args.nImg
    KEEP_NEF = False
    # max number of images to download
    DOWNLOAD_LIMIT = 30000
    #DOWNLOAD_MAX_SIZE =  20000
    DOWNLOAD_MAX_SIZE = 1000000
    RESIZE = False
    #resize_w, resize_h = 356, 536
    RESIZE_WIDTH = 712
    RESIZE_HEIGHT = 1072  
    
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    ## DICTIONARIES ASSIGNING 'VIRTUAL_SENSOR' TO 'POSITION'
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    DEPO_VS = {
#            'AD01': {'AD01': {'displacement', 'temperature_rock'}},      #OK
#            'AD02': {'AD02': {'displacement'}},                          #OK
#            'AD03': {'AD03': {'displacement'}},                          #OK
#            'AD04': {'AD04': {'displacement', 'temperature_rock'}},      #OK
#            'AD05': {'AD05': {'displacement'}},                          #OK
#            'AD06': {'AD06': {'displacement', 'temperature_rock'}},      #OK
#            'AD07': {'AD07': {'temperature_rock'}},                      #OK
#            'AD08': {'AD08': {'temperature_rock'}},                      #OK
#            'AD09': {'AD09': {'temperature_rock'}},                      #OK
            
           'DH00': {'RD01': {'gps_differential__batch__daily','gps_differential__rtklib__daily'}},                    #                    OK
           'DH05': {'DI55': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #                    OK
           'DH06': {'RA01': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #000deg,             OK
           'DH07': {'DI57': {'gps_differential__batch__daily','gps_differential__rtklib__daily'}},                    #                    OK
           'DH09': {'RA02': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #000deg,             OK
           'DH12': {'BH03': {'gps_differential__batch__daily','gps_differential__rtklib__daily'}},                    #                    OK
           'DH15': {'DI02': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #225/131/281/246deg, OK
           'DH17': {'DI07': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #105/277/321/334deg, OK
           'DH21': {'LS01': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #110deg,             OK
           'DH23': {'LS04': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #225deg,             OK
           'DH25': {'BH07': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #240deg,             OK
           'DH27': {'BH09': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #130deg,             OK
           'DH29': {'ST02': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #004deg,             OK
           'DH31': {'ST05': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #124deg,             OK 
           'DH33': {'GU02': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #335deg,             OK
           'DH35': {'GU03': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #010deg,             OK
           'DH39': {'RG01': {'gps_differential__batch__daily','gps_differential__rtklib__daily'}},                    #                    OK
           'DH41': {'GG52': {'gps_differential__batch__daily','gps_differential__rtklib__daily'}},                    #                    OK
           'DH43': {'GG01': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #325deg,             OK
           'DH44': {'GG02': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #205/041deg          OK
           'DH55': {'BH10': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #280/270/294/312deg  OK
           'DH56': {'RL01': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #309deg              OK
           'DH57': {'GU04': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #014deg,             OK
           'DH62': {'BH12': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #140deg,             OK
           'DH63': {'BH13': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #090/000deg,         OK
           'DH64': {'LS05': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #318/330deg,         OK
           'DH66': {'GG66': {'gps_differential__batch__daily','gps_differential__rtklib__daily'}}, #,                                      OK
           'DH67': {'GG67': {'gps_differential__batch__daily','gps_differential__rtklib__daily'}}, #,                                      OK
           'DH70': {'RA03': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #000deg,             OK
           'DH81': {'RAND': {'gps_differential__batch__daily','gps_differential__rtklib__daily'}},                    #,                   OK
           'DH82': {'WYS1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #170deg,             OK
           'DH83': {'LS11': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #262deg,             OK
           'DH84': {'LS12': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #006deg,             OK
           'DH86': {'DI03': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #000deg,             OK
           'DH87': {'DI04': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #000deg,             OK
           'DH88': {'LS06': {'gps_differential__rtklib__daily','gps_inclinometer'}},                                  #000deg,             OK
           'DH89': {'SA01': {'gps_differential__rtklib__daily','gps_inclinometer'}},                                  #000deg
           'DH91': {'SATT': {'gps_differential__rtklib__daily'}},                                                     #
           'DH92': {'SA92': {'temperature_rock'}},                                                                    #

           
           'DH13': {'DH13': {'vaisalawxt520prec', 'vaisalawxt520windpth'}},                                           #
           'DH42': {'DH42': {'vaisalawxt520prec', 'vaisalawxt520windpth'}},                                           #
           'DH68': {'DH68': {'vaisalawxt520prec', 'vaisalawxt520windpth'}},                                           #
           'DH69': {'DH69': {'vaisalawxt520prec', 'vaisalawxt520windpth'}},                                           #
           'DH73': {'DH73': {'radiometer__conv'}},                                                                    #

           'PE01': {'DIS1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE02': {'DIS2': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE03': {'RIT1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE04': {'GRU1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE05': {'JAE1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE06': {'SCH1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE07': {'MUA1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE08': {'LAR1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE09': {'LAR2': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK
           'PE10': {'COR1': {'gps_differential__batch__daily','gps_differential__rtklib__daily','gps_inclinometer'}}, #,                   OK

        # #    'JJ01': {'JJ01': {'resistivity_rock', 'temperature_rock'}},      #OK
        # #    'JJ02': {'JJ02': {'resistivity_rock', 'temperature_rock'}},      #OK
        # #    'JJ03': {'JJ03': {'resistivity_rock', 'temperature_rock'}},      #OK
        # #    'JJ04': {'JJ04': {'resistivity_rock', 'temperature_rock'}},      #OK
        # #    'JJ05': {'JJ05': {'resistivity_rock', 'temperature_rock'}},      #OK
        # #    'JJ06': {'JJ06': {'resistivity_rock', 'temperature_rock'}},      #OK
        # #    'JJ07': {'JJ07': {'resistivity_rock', 'temperature_rock'}},      #OK
        # #    'JJ08': {'JJ08': {'temperature_fracture'}},                      #OK
        # #    'JJ09': {'JJ09': {'resistivity_rock', 'temperature_rock'}},      #OK
        # #    'JJ10': {'JJ10': {'temperature_fracture'}},                      #OK
        # #    'JJ13': {'JJ13': {'temperature_rock'}},                          #
        # #    'JJ14': {'JJ14': {'temperature_rock'}},                          #
        # #    'JJ18': {'JJ18': {'temperature_rock'}},                          #
        # #    'JJ22': {'JJ22': {'temperature_rock'}},                          #
        # #    'JJ25': {'JJ25': {'temperature_rock'}},                          #


            'MH01': {'MH01': {'displacement', 'temperature_fracture', 'temperature_rock'}},                                    #ok
            'MH02': {'MH02': {'displacement', 'temperature_fracture', 'temperature_rock'}},                                    #ok
            'MH03': {'MH03': {'displacement', 'temperature_fracture', 'temperature_rock'}},                                    #ok
            'MH04': {'MH04': {'displacement', 'temperature_fracture', 'temperature_rock'}},                                    #ok
            'MH05': {'MH05': {'resistivity_fracture', 'temperature_fracture'}},                                                #ok but resistivity
            'MH06': {'MH06': {'displacement', 'temperature_rock'}},                                                            #ok
            'MH07': {'MH07': {'resistivity_fracture', 'temperature_fracture'}},                                                #ok but resistivity
            'MH08': {'MH08': {'displacement', 'temperature_rock'}},                                                            #ok
            'MH09': {'MH09': {'displacement', 'temperature_fracture'}},                                                        #ok
            'MH10': {'MH10': {'resistivity_rock', 'temperature_rock'}},                                                        #ok but resistivity
            'MH11': {'MH11': {'resistivity_rock', 'temperature_rock'}},                                                        #ok but resistivity_rock
            'MH12': {'MH12': {'resistivity_rock', 'temperature_rock'}},                                                        #ok but resistivity
            'MH18': {'MH18': {'displacement'}},                                                                                #ok
            'MH20': {'MH20': {'displacement'}},                                                                                #ok
            'MH21': {'MH21': {'displacement'}},                                                                                #ok
            'MH22': {'MH22': {'displacement'}},                                                                                #ok
            'MH27': {'MH27': {'temperature_rock'}},                                                                            #ok
            'MH30': {'MH30': {'temperature_rock'}},                                                                            #ok
            'MH46': {'MH46': {'temperature_rock'}},                                                                            #ok
            'MH47': {'MH47': {'temperature_rock'}},                                                                            #ok
   
            'MH15': {'MH15': {'radiometer__conv'}},                                                                            #ok
            'MH25': {'MH25': {'vaisalawxt520prec', 'vaisalawxt520windpth'}},                                                   #ok
            'MH51': {'MH51': {'vaisalawxt520prec', 'vaisalawxt520windpth'}},                                                   #ok

            'MH33': {'MH33': {'gps_differential__batch__daily', 'gps_differential__rtklib__daily', 'gps_inclinometer'}},         #ok
            'MH34': {'MH34': {'gps_differential__batch__daily', 'gps_differential__rtklib__daily', 'gps_inclinometer'}},         #ok
            'MH35': {'MH35': {'gps_differential__batch__daily', 'gps_differential__rtklib__daily', 'gps_inclinometer'}},         #ok
            'MH40': {'MH40': {'gps_differential__batch__daily', 'gps_differential__rtklib__daily'}},                             #ok
            'MH42': {'MH42': {'gps_differential__batch__daily', 'gps_differential__rtklib__daily'}},                             #ok
            'MH43': {'MH43': {'gps_differential__rtklib__daily'}},
            }

    DEPO_GPS = ['MH33', 'MH34', 'MH35', 'MH40', 'MH42', 'MH43']

    DEPO_IMG = {
            # 'MH19': {'binary__mapped'}
            }
        
    
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    ## Create PATH_DATA if it does not exist yet
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    
    os.makedirs(PATH_DATA, exist_ok=True)
    os.makedirs(PATH_DATA + '/gnss_data_raw', exist_ok=True)
    os.makedirs(PATH_DATA + '/gnss_derived_data_products', exist_ok=True)
    os.makedirs(PATH_DATA + '/timelapse_images', exist_ok=True)
    os.makedirs(PATH_DATA + '/timeseries_data_raw', exist_ok=True)
    os.makedirs(PATH_DATA + '/timeseries_derived_data_products', exist_ok=True)
    os.makedirs(PATH_DATA + '/timeseries_sanity_plots', exist_ok=True)
    os.makedirs(PATH_DATA + '/timeseries_sanity_plots/' + 'channelwise', exist_ok=True)
        
    
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    ## GSN timeseries data
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    
    # Loop through keys (depo) in dictionary DEPO_VS
    for depo, depo_line in DEPO_VS.items():
#         print(depo)
#         print(depo_line)
        for label in depo_line:
            print('Working on data of label {:s} and depo {:s}.'.format(label, depo))

        if depo[0:2] == 'JJ':
            deployment = 'jungfraujoch'
        elif depo[0:2] == 'AD':
            deployment = 'adm'
        elif depo[0:2] == 'MH':
            deployment = 'matterhorn'
        elif depo[0:2] == 'DH':
            deployment = 'dirruhorn'
        elif depo[0:2] == 'PE':
            deployment = 'permos'
        position = int(depo[2:4])
        
    # Step 1: Query data from GSN server and save it locally
    # Loop through values for a given key of dictionary DEPO_VS
        for vsensor in sorted(list(depo_line[label])):
#             print(label, depo, deployment, vsensor, position)
            if GSN2RAW:
                print('Data from label {:s}, position {:s} and virtual sensor {:s} is imported.'.format(label, depo, vsensor))
                print('Time interval: ',TBEG, ' to ', TEND)
                df = get_GSNdata(deployment, position, vsensor, TBEG, TEND, download_max_size = DOWNLOAD_MAX_SIZE)
                if len(df) > 1:
                    print('Saving raw data')
                    if vsensor == 'gps_differential__batch__daily':
                        save_data(df, PATH_DATA + '/gnss_data_raw/', label, vsensor, deployment, file_type='csv')
                        save_data(df, PATH_DATA + '/gnss_data_raw/', label, vsensor, deployment, yearly=True, iso=True, file_type='csv')
                    elif vsensor == 'gps_differential__rtklib__daily':
                        save_data(df, PATH_DATA + '/gnss_data_raw/', label, vsensor, deployment, file_type='csv')
                        save_data(df, PATH_DATA + '/gnss_data_raw/', label, vsensor, deployment, yearly=True, iso=True, file_type='csv')
                    else:
                        save_data(df, PATH_DATA + '/timeseries_data_raw/', label, vsensor, deployment, file_type='csv')
                        save_data(df, PATH_DATA + '/timeseries_data_raw/', label, vsensor, deployment, yearly=True, iso=True, file_type='csv')
                print()
                
    # Loop through values for a given key of dictionary DEPO_VS
        for vsensor in list(vs for vs in sorted(list(depo_line[label]))):
            df={}

            # Step 2: Load data
            if LOAD:
                print('Loading data', label, depo, deployment, vsensor)
                if vsensor == 'gps_differential__batch__daily':
                    df = load_data(PATH_DATA + '/gnss_data_raw/', label, vsensor, file_type='csv')
                    df_raw = df.copy()
                elif vsensor == 'gps_differential__rtklib__daily':
                    df = load_data(PATH_DATA + '/gnss_data_raw/', label, vsensor, file_type='csv')
                    df_raw = df.copy()
                else:
                    df = load_data(PATH_DATA + '/timeseries_data_raw/', label, vsensor, file_type='csv')
                    df_raw = df.copy()
            
    # Step 3 (optional): Filter data (according to reference values)
            if FILTER:
                print('Filtering data', label, depo, deployment, vsensor)
                df = filter_data(df, depo)
            df_filt = df.copy()
            
    # Step 4 (optional): Clean data manually
            if CLEAN:
                print('Cleaning data', label, depo, deployment, vsensor)
                try:
                    df = clean_data(df, depo, vsensor)
                except Exception as i:
                    print('An issue with data cleaning raised - please check the data', i)
                    pass
            df_filt_clean = df.copy()
    
    # Step 5: Correct GPS data for moving reference HOGR
            if vsensor == 'gps_differential__rtklib__daily':
#                 print(df)
#                 print(label)
                print(df.index[0])
                print(df.index[-1])
                t1 = df.index[0]
                t2 = df.index[-1]
                trng_reindex = pd.date_range(t1, t2, freq='24h')

                label_list = ['MH33', 'MH34', 'MH35', 'MH40', 'MH43']
                if label in label_list:
                    print("Correcting GPS for moving reference HOGR")

                    year_list = list(range(2008, dt.datetime.today().year + 1))
                    depo_HOGR = 'MH42'
                    vsensor = 'gps_differential__rtklib__daily'
                    # #vsensor = 'gps_differential__batch__daily'

                    df_HOGR = load_data(PATH_DATA + '/gnss_derived_data_products/', depo_HOGR, vsensor, year=year_list, file_type='csv')
                    # print(df_HOGR)
                    
                    df_HOGR = df_HOGR.reindex(trng_reindex)


                    df['e'] = df['e'] + df_HOGR['e [m]'] - df_HOGR['e [m]'][0]
                    df['n'] = df['n'] + df_HOGR['n [m]'] - df_HOGR['n [m]'][0]
                    df['h'] = df['h'] + df_HOGR['h [m]'] - df_HOGR['h [m]'][0]
                df_gps_corr = df.copy
                        
    # Step 6 (optional): Aggregate data
            if AGGREGATION_INTERVAL > 0:
                print('Aggregating data', label, depo, deployment, vsensor)
                try:
                    if vsensor == 'vaisalawxt520prec':
                        f = {'position':'mean', 'device_id':'mean', 'rain_accumulation':'sum', 'rain_duration':'sum', 'rain_intensity':'mean', 'rain_peak_intensity':'max', 'hail_accumulation':'sum', 'hail_duration':'sum', 'hail_intensity':'mean', 'hail_peak_intensity':'mean'}
                        df = df.groupby(pd.Grouper(freq=str(AGGREGATION_INTERVAL) + 'Min')).agg(f)
                    elif vsensor == 'vaisalawxt520windpth':
                        f = {'position':'mean', 'device_id':'mean', 'wind_direction_minimum':'min', 'wind_direction_average':'mean', 'wind_direction_maximum':'max', 'wind_speed_minimum':'min', 'wind_speed_average':'mean', 'wind_speed_maximum':'max', 'temp_air':'mean', 'temp_internal':'mean', 'relative_humidity':'mean', 'air_pressure':'mean'}
                        df = df.groupby(pd.Grouper(freq=str(AGGREGATION_INTERVAL) + 'Min')).agg(f)
                    elif vsensor == 'gps_differential__batch__daily':
                        print('Copying GPS data, no aggregation.')
#                         df = df
                    elif vsensor == 'gps_differential__rtklib__daily':
                        print('Copying GPS data, no aggregation.')
#                         df = df
                    else:
                        df = df.groupby(pd.Grouper(freq=str(AGGREGATION_INTERVAL) + 'Min')).aggregate(np.mean)
                        #df = df.groupby(pd.TimeGrouper(freq=str(AGGREGATION_INTERVAL) + 'Min')).aggregate(np.median)
                        
                except:
                    print('There was an issue with time aggregation')
                    pass
            df_agg = df.copy()

    # Step 7: Export to csv
            if EXPORT2CSV:
                print('Exporting data', label, depo, deployment, vsensor)
                df_csv = df.copy()
                for col in ['device_id',  'ref1',  'ref2',  'ref3',  'ref4',  'ref5',  'ref6']:
                    try:
                        df_csv.drop([col], axis=1, inplace=True)
                    except:
                        pass
    
                try:
                    if vsensor == 'gps_differential__batch__daily':
                        save_data(df_csv, PATH_DATA + '/gnss_derived_data_products/', label, vsensor, deployment, yearly=True, iso=True, file_type='csv', float_format=4)
                    elif vsensor == 'gps_differential__rtklib__daily':
                        save_data(df_csv, PATH_DATA + '/gnss_derived_data_products/', label, vsensor, deployment, yearly=True, iso=True, file_type='csv', float_format=4)
                    else:
                        save_data(df_csv, PATH_DATA + '/timeseries_derived_data_products/', label, vsensor, deployment, yearly=True, iso=True, file_type='csv', float_format=4)
                except:
                    pass
                
    # Step 8: Plot for sanity check
            if SANITY_PLOT:
                print('Plotting data', label, depo, deployment, vsensor)
                for i, col in enumerate(df_agg):
                    if col not in ['device_id', 'device_type', 'position', 'label', 'ref1', 'ref2', 'ref3', 'ref4', 'ref5', 'ref6', 'sd_e', 'sd_n', 'sd_h', 'version', 'reference_label', 'processing_time', 'ratio_of_fixed_ambiguities']:
                        #print(col)
                        #plt.figure(figsize=(24/2.54, 18/2.54), facecolor='w', edgecolor='k')
                        plt.figure(figsize=(40/2.54, 18/2.54), facecolor='w', edgecolor='k')
                        plt.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.85) # left=0.125, right=0.9, bottom=0.1, top=0.9, wspace=0.2, hspace=0.2
                        plt.axhline(y=0, linestyle='dotted', linewidth=1, color='grey')

                        plt.plot(df_raw[col][~df_raw[col].isnull()] , color='C0', label='raw')
                        plt.plot(df_filt[col], color='C1', label='filtered')
                        plt.plot(df_filt_clean[col], color='C2', label='cleaned')
                        plt.plot(df_agg[col], color='C3', label='aggregated')
                        diff = df_agg[col].max() - df_agg[col].min()
                        try:
                            plt.ylim([df_agg[col].min() - diff/10, df_agg[col].max() + diff/10])
                        except:
                            pass
                        plt.title(' Label: {:s} \n Deployment: {:s} \n Position: {:s} \n VSENSOR: {:s}\n VARIABLE: {:s}'.format(label,deployment,depo[2:4],vsensor,col), loc='left')
                        #plt.legend(['raw','+ filtered', '+ cleaned', '+ aggregated'], loc="upper left")
                        plt.legend(loc="upper left")
                        #print(mpl.rcParams['agg.path.chunksize'])
                        mpl.rcParams['agg.path.chunksize'] = 10000
                        # start from fixed date
                        if deployment == 'dirruhorn' or deployment == 'permos':
                            plt.xlim([TBEG, TEND])
#                             plt.xlim([dt.datetime(2011,1,1), TEND])
                        else :
                            plt.xlim([TBEG, TEND])

                        plt.savefig('{:s}/timeseries_sanity_plots/channelwise/{:s}_{:s}_{:s}.png'.format(PATH_DATA , label, vsensor, col), dpi=300)
                        plt.close('all')
                        mpl.rcParams['agg.path.chunksize'] = 0
                
                #plt.figure(figsize=(24/2.54, 18/2.54), facecolor='w', edgecolor='k')
                fig = plt.figure(figsize=(40/2.54, 18/2.54), facecolor='w', edgecolor='k')
                plt.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.85) # left=0.125, right=0.9, bottom=0.1, top=0.9, wspace=0.2, hspace=0.2

                for i, col in enumerate(df_agg):
                    if col not in ['device_id', 'device_type', 'position', 'label', 'ref1', 'ref2', 'ref3', 'ref4', 'ref5', 'ref6', 'sd_e', 'sd_n', 'sd_h', 'version', 'reference_label', 'processing_time', 'ratio_of_fixed_ambiguities']:
                        #print(df_agg[col])
                        if ~df_agg[col].isnull().all():
                            # print(col)
                            if col == 'e':
                                col_label = 'East'
                            elif col == 'n':
                                col_label = 'North'
                            elif col == 'h':
                                col_label = 'Altitude'
                            else:
                                col_label = col
                            plt.axhline(y=0, linestyle='dotted', linewidth=1, color='grey')
                            plt.plot(df_agg[col], label=col_label)
                            #print(p1)
                            #print(col_label)	
                plt.title(' Label: {:s} \n Deployment: {:s} \n Position: {:s} \n VSENSOR: {:s}\n All channels'.format(label,deployment,depo[2:4],vsensor), loc='left')
                # any legend location is too slow
                # #plt.legend()
                plt.legend(loc="upper left")
                
                # start from fixed date
                if deployment == 'dirruhorn' or deployment == 'permos':
                    plt.xlim([TBEG, TEND])
#                     plt.xlim([dt.datetime(2011,1,1), TEND])
                else :
                    plt.xlim([TBEG, TEND])
                #plt.tight_layout()
                #if fig.get_axes(): # Do stuff when the figure isn't empty.
                plt.savefig('{:s}/timeseries_sanity_plots/{:s}_{:s}_all.png'.format(PATH_DATA , label, vsensor), dpi=300)
                plt.close('all')
                print()



    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## 
    ## GSN images
    ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##
    print(GSN2IMG)
    if GSN2IMG:
        # Loop through keys (depo) in dictionary DEPO_IMG
        for depo in DEPO_IMG:
            if depo[0:2] == 'MH':
                deployment = 'matterhorn'
            position = int(depo[2:4])
            
            # Loop through values for a given key of dictionary DEPO_IMG
            for vsensor in list(DEPO_IMG[depo]):
                print(deployment, position, vsensor, TBEG, TEND)
                df = get_GSNimg(PATH_DATA + '/timelapse_images/', deployment, position, vsensor, TBEG, TEND, n_img = N_IMG, keep_nef = KEEP_NEF, download_limit = DOWNLOAD_LIMIT, download_max_size = DOWNLOAD_MAX_SIZE, resize = RESIZE, resize_w = RESIZE_WIDTH, resize_h = RESIZE_HEIGHT)



def parseArguments():
    import argparse
    import datetime as dt
    
    def true_or_false(arg):
        ua = str(arg).upper()
        if 'TRUE'.startswith(ua):
           return True
        elif 'FALSE'.startswith(ua):
           return False
        else:
           pass  #error condition maybe?
    
    # Create argument parser
    parser = argparse.ArgumentParser()

    # Positional (mandatory) arguments
    #parser.add_argument("positionalArgument1", help="Positional argument 1", type=float)

    # Optional arguments
    parser.add_argument("--gsn2raw", help="Get raw data from GSN and store locally (step1, default: %(default)s)", type=true_or_false, default=True)
    parser.add_argument("--load", help="Load locally stored data (step2, default: %(default)s)", type=true_or_false, default=True)
    parser.add_argument("--filter", help="Filter data (according reference values, step3, default: %(default)s)", type=true_or_false, default=True)
    parser.add_argument("--clean", help="Clean data manually (step4, default: %(default)s)", type=true_or_false, default=True)
    parser.add_argument("--aggregationInterval", help="Aggregate data over timewindows of X minutes (step5, default: %(default)s)", type=int, default=60)
    parser.add_argument("--export2csv", help="Export data in yearly csv files (step6, default: %(default)s)", type=true_or_false, default=True)
    parser.add_argument("--sanityPlot", help="Plot for sanity check (step7, default: %(default)s)", type=true_or_false, default=True)
    parser.add_argument("--gsn2img", help="Downlod images from GSN and convert NEF to JPEG (default: %(default)s)", type=true_or_false, default=True)
    parser.add_argument("--nImg", help="Number of images to download and convert (default: %(default)s)", type=int, default=10)
    parser.add_argument("-p", "--path", help="Relative path for the data output (default: %(default)s)", type=str, default='./data')
    parser.add_argument("-y1", "--yearBegin", help="Year begin (default: %(default)s)", type=int, default=2008)
    parser.add_argument("-m1", "--monthBegin", help="Month begin (default: %(default)s)", type=int, default=1)
    parser.add_argument("-d1", "--dayBegin", help="Day begin (default: %(default)s)", type=int, default=1)
    parser.add_argument("-y2", "--yearEnd", help="Year end (default: today)", type=int, default=dt.datetime.now().year)
    parser.add_argument("-m2", "--monthEnd", help="Month end (default: today)", type=int, default=dt.datetime.now().month)
    parser.add_argument("-d2", "--dayEnd", help="Day end (default: today)", type=int, default=dt.datetime.now().day)

    # Print version
    #parser.add_argument("--version", action="version", version='%(prog)s - Version 1.0')

    # Parse argument
    return parser.parse_args()

if __name__ == "__main__":
    main()