MatfiletoNetCDF_CL_ARTofMELT_v1.py

# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# This python script converts the processed .mat file for native Ceilometer observations into a user-friendly NetCDF file. It also converts all time averaged .mat files into their respecitve netCDF files. These final (four) netCDF files are found in the Bolin Centre Database. Data from the ARTofMELT expedition in Spring 2023. This code is written by Sonja Murto (2024).

# +
# import modules
import sys,glob,os
import numpy as np
import math
import pandas as pd
import matplotlib.dates as mdates
import xarray as xr
import time,datetime
import itertools
import matplotlib.pyplot as plt
import string
import scipy.io
import matplotlib.ticker as mticker
from collections import Counter


#metpy functions for thermodynamical variable conversions
import metpy.calc as mpcalc
from metpy.units import units

import warnings
warnings.filterwarnings('ignore')


# +
# define functions
def datestring_fromtuple(date=(2001,2,15,12)):
    '''
    Input: a date tuple (YYYY,M,D,h)
    Output: Date string: 'YYYYMMDD_hh'
    '''
    strdate=str(date[0])+"{:0>2}".format(str(date[1]))+"{:0>2}".format(str(date[2]))+'_'+"{:0>2}".format(str(date[3]))
    return strdate

def datetuple(date):
    '''date is as datetime object, i.e. from pandas using .to_pydatetime()
    '''
    return date.year,date.month,date.day,date.hour

def strtodatetuple(datestr='19900109_12',format_st='%Y%m%d_%H',returntuple=True):
    if returntuple:   
        return datetuple(datetime.datetime.strptime(datestr, format_st))
    else:
        return datetime.datetime.strptime(datestr, format_st)


def returndatetime_fromdoy(doys, year=2023):
    '''
    fucntion to return datetime from DOY dates. 
    Takes the leap year into account if given for a leap year
    Input: list of DOYs
    Output: List of Dates (in datetime) corresponding to the given doys
    Note, returns in microseconds.
    If second rounded: [pd.to_datetime(T).round('1s') for T in date]
    
    '''
    jdate_frac=doys
    date=[]
    for d in jdate_frac:
        year,julian = [year,d]
        date.append(datetime.datetime(year, 1, 1)+datetime.timedelta(days=julian -1))
    return date

def getList(dict): 
    '''
    This function returns a list of keys for a pandas dictionary
    '''
    list = [] 
    for key in dict.keys(): 
        list.append(key) 
    return list

def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100,returncolors=False):
    '''
    This function creates a new colormap from the cmap chosen, taking the limits from minval and maxval, and the number of colors as n.
    The default is a colormap as LinearSegmentedColormap, which can be used as colormap for imshow, contourf...
    If you assign returncolors as True, the there will be an array of n amount of colors from the chosen colormap. This can be used in plots when each event has a color (n=50)
    '''
    new_cmap = mcolors.LinearSegmentedColormap.from_list('new_cmap', plt.get_cmap(cmap)(np.linspace(minval, maxval, n)))
    if returncolors:
        new_cmap=plt.get_cmap(cmap)(np.linspace(minval, maxval, n))
    return new_cmap   



# -

#define directories - change according to your computer settings
cwd = os.getcwd() 
savefigpath = cwd + '/FIGS/Example_figs/' #directory for figures
#directory where to find the data (.mat files)
load_data = r'/Volumes/My Passport for Mac/MISU_sticka/WORK/viktigapapper/ARTofMELT_2023/AoM_data_mod/WX/cl31_output/'


# # Convert ceilometer data from .mat files into .nc files

# Create secondary dimensions
cloud_layer_levs=np.arange(1,4,1).astype('int32') # 3 levels
sc_layer_levs=np.arange(1,6,1).astype('int32') # 5 levels
range_levs=np.arange(1,771,1).astype('int32') # 770 levels

print(cloud_layer_levs,sc_layer_levs)

# ## 30s data

# load 30s native data
CL_SM_org = scipy.io.loadmat(load_data + 'CL31_ceilometer_ARTofMELT_20230507_20230613_30s_v01.mat',struct_as_record=True) #30s

CL_SM_org.keys()

# get the data in the mat file
Names=CL_SM_org['cl31'].dtype.names
ndata = {n: CL_SM_org['cl31'][n][0, 0] for n in Names}
print(Names)

#create time dimension
Time_steps=returndatetime_fromdoy(np.array(list(itertools.chain.from_iterable(ndata['doy'])),dtype=float))
Times_nomicrosec=[pd.to_datetime(T).round('1s') for T in Time_steps]
Time_steps_dt64_org=[np.datetime64(t) for t in Time_steps]
Time_steps_dt64_org=np.array(Time_steps_dt64_org,dtype='datetime64[ns]') #time dimension

# +
# convert data into pandas dataframes
DF_30s=pd.DataFrame(index=range(len(ndata['doy'])),)
DF_30s['doy']=np.array(list(itertools.chain.from_iterable(ndata['doy'])),dtype=float)
DF_30s['year']=np.array([Times_nomicrosec[i].year for i in range(len(Time_steps))],dtype=int)
DF_30s['month']=np.array([Times_nomicrosec[i].month for i in range(len(Time_steps))],dtype=int)
DF_30s['day']=np.array([Times_nomicrosec[i].day for i in range(len(Time_steps))],dtype=int)
DF_30s['hour']=np.array([Times_nomicrosec[i].hour for i in range(len(Time_steps))],dtype=int)
DF_30s['minute']=np.array([Times_nomicrosec[i].minute for i in range(len(Time_steps))],dtype=int)
DF_30s['second']=np.array([Times_nomicrosec[i].second for i in range(len(Time_steps))],dtype=int)
DF_30s['cloudcode']=np.array(list(itertools.chain.from_iterable(ndata['cloudcode'])),dtype=int)
DF_30s['vert_vis']=np.array(list(itertools.chain.from_iterable(ndata['vert_vis'])),dtype=float)
DF_30s['high_sig']=np.array(list(itertools.chain.from_iterable(ndata['high_sig'])),dtype=float)

DF_30s.head()
# -

#check for nans - no nans!
DF_30s.iloc[DF_30s[DF_30s.doy.isna()].index]

#get a list of columns
DF_30s.columns

# +
# save each variable into Data Arrays, with attribute data

#1D variables: dimension time

da_doy=xr.DataArray(data=np.array(DF_30s['doy']).astype('float32'),name="day_of_year",
                dims=["time"],coords=dict(time=Time_steps_dt64_org),
                  attrs=dict(type="float32",dimension="time",units="1",long_name="Day of Year",
                             description="time as decimal day of year"),) #added as 1D

#date and time in separate arrays; microseconds approximated to seconds
da_year=xr.DataArray(data=np.array(DF_30s['year']).astype('int32'),name="year",dims=["time"],
                  coords=dict(time=Time_steps_dt64_org),
                  attrs=dict(type="int32",dimension="time",units="1",long_name="Year"),) #added as 1D

da_month=xr.DataArray(data=np.array(DF_30s['month']).astype('int32'),name="month",dims=["time"],
                  coords=dict(time=Time_steps_dt64_org),
                      attrs=dict(type="int32",dimension="time",units="1", long_name="Month"),) #added as 1D

da_day=xr.DataArray(data=np.array(DF_30s['day']).astype('int32'),name="day",dims=["time"],
                  coords=dict(time=Time_steps_dt64_org),
                    attrs=dict(type="int32",dimension="time",units="1",long_name="Day"),) #added as 1D

da_hour=xr.DataArray(data=np.array(DF_30s['hour']).astype('int32'),name="hour",dims=["time"],
                  coords=dict(time=Time_steps_dt64_org),
                     attrs=dict(type="int32",dimension="time",units="1",long_name="Hour"),) #added as 1D

da_min=xr.DataArray(data=np.array(DF_30s['minute']).astype('int32'),name="minute",dims=["time"],
                  coords=dict(time=Time_steps_dt64_org),
                    attrs=dict(type="int32",dimension="time",units="1",long_name="Minute"),) #added as 1D

da_sec=xr.DataArray(data=np.array(DF_30s['second']).astype('int32'),name="second",dims=["time"],
                  coords=dict(time=Time_steps_dt64_org),
                     attrs=dict(type="int32",dimension="time",units="1",long_name="Second",
                               description="Time averaged to closest second"),) #added as 1D

#add cloudcode flag
da_cloudcode=xr.DataArray(data=np.array(DF_30s['cloudcode']).astype('int32'),name="flag_cloudcode",dims=["time"],