<jupyter_start><jupyter_text>---How to make BSH data ATMODAT-compliant version 1 <NAME>, DKRZ; last edited Oct 22, 2021Tested on DKRZ jupyterhub with using Python 3 bleeding edge kernel (module anaconda3/bleeding_edge) TODO:* Wind direction as numeric with bounds* Add geometry  Load required libraries<jupyter_code>import os, glob
import pandas as pd
import numpy as np
import xarray as xr
from datetime import datetime
import matplotlib.pyplot as plt 
%matplotlib inline
%config InlineBackend.print_figure_kwargs = {'bbox_inches':None}
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"              #-- to show all outputs in a single cell<jupyter_output><empty_output><jupyter_text>Input Data  Get annual ascii files **Data description:**https://www.bsh.de/DE/DATEN/Klima-und-Meer/Wetterlagen-und-Stuerme/wetterlagen-und-stuerme_node.html  Define input data  Hard-code file information<jupyter_code>folder = "data/BSH"        #-- folder where input data are stored
fnroot = "lwtns06"         #-- file name root
ifiles = sorted(glob.glob(folder+'/'+fnroot+'_????.dat'))<jupyter_output><empty_output><jupyter_text>Hard-code missing header info<jupyter_code>lwtvars  = [ 'lamb', 'flowdir', 'winddir', 'stormsev' ]
datevars = ['year', 'month', 'day']
headers  = datevars + lwtvars
headers
colspecs = [ (0,4), (5,7), (8,10), (11,13), (14,16), (17,19), (20,23)]<jupyter_output><empty_output><jupyter_text>Hard-code variable metadata: long_names and units<jupyter_code># long name
lwtvars_ln = ['lamb weather type', 'flow direction', 'wind direction', 'storm severity']
# units
lwtvars_un = ['1', '1', '1', '1']
# https://doi.org/10.1088/0026-1394/52/1/40: in SI, dimensionless quantities are deemed to have no unit or the coherent derived unit 'one'. 
# fill value
lwtvars_fl = [-9999]*len(lwtvars)<jupyter_output><empty_output><jupyter_text>Define with what numeric values the variables will be encoded<jupyter_code>lwtvars_cd = [[1, 1], [-1, 2], [10, 10], [0, 100]]   #-- parameters start and step in range(start, stop, step)<jupyter_output><empty_output><jupyter_text>Read in all input files into a single pandas dataframe<jupyter_code>idf         = pd.concat([pd.read_fwf(ifile,index_col=None,header=None, colspecs=colspecs) 
                  for ifile in ifiles],sort=False)
idf.columns = headers
idf['time'] = pd.to_datetime(idf[datevars])    #-- alternatively idf.index=
df          = idf.drop(datevars, axis = 1)
df.head(4)
#print(df.lamb.dtype)   #-- dtype('O'), i.e. pandas objects, which corresponds to Numpy string_, or unicode_ types.<jupyter_output><empty_output><jupyter_text>Recode the data using flag_values and flag_meanings  Plot frequency distribution of classes in each variable<jupyter_code>figsiz = 5
nrows, ncols = 1, len(lwtvars)
fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(ncols*figsiz,nrows*figsiz))
fig.suptitle('Frequency distribution as histograms (period)', fontsize=16)
for varindex in range(len(lwtvars)):
    stdout = df[lwtvars[varindex]].value_counts().plot(kind='bar', ax=axes[varindex], title=lwtvars_ln[varindex])<jupyter_output><empty_output><jupyter_text>Create dictionaries with flag_values and flag_meanings for variable and recode data<jupyter_code>vardict = dict((lwtvars[varindex], {}) for varindex in range(len(lwtvars)))

for varindex in range(len(lwtvars)):
    freqbyclass = df[lwtvars[varindex]].value_counts()
    #-- NB: arange would create int64 object which cannot be written to json
    codes       = range(lwtvars_cd[varindex][0], lwtvars_cd[varindex][0] + len(freqbyclass) * lwtvars_cd[varindex][1], 
                        lwtvars_cd[varindex][1])
    flagdict    = dict(zip(freqbyclass.index, codes))  
    vardict[lwtvars[varindex]].update(flagdict)
    
    print(vardict[lwtvars[varindex]].values())
    df = df.replace({ lwtvars[varindex]: flagdict })  #-- the replace argument is a nested directory

    print(f'\
The {fnroot} BSH {lwtvars[varindex]} variable contains in total {len(freqbyclass)} classes.\n\
These are: {freqbyclass.index.values}')
    print(f'\
The following encoding will be used {flagdict}')

    
df.head(4)<jupyter_output><empty_output><jupyter_text>Save nested dictionary with flag encodings as json file<jupyter_code>with open(fnroot+'_flagcodes.json', 'w') as fp:
    json.dump(vardict, fp)<jupyter_output><empty_output><jupyter_text>Edit metadata  Convert pandas dataframe to xarray dataset<jupyter_code>ds = df.set_index(['time']).to_xarray()<jupyter_output><empty_output><jupyter_text>Encoding for data variables<jupyter_code>varenc     = dict(zlib=False, complevel=1, dtype='i4', _FillValue= -9999)
varsenc    = {var: varenc for var in ds.data_vars}<jupyter_output><empty_output><jupyter_text>Encoding for coordinate variables<jupyter_code>dimenc     = dict(zlib=False, complevel=1, dtype='i4', _FillValue=None)
dimsenc    = {var: varenc for var in ds.dims}<jupyter_output><empty_output><jupyter_text>Data variable metadata<jupyter_code>for count, var in enumerate(lwtvars):
    ds[var].attrs = {'flag_meanings': str(list(vardict[var].keys())).replace('[','').replace(']','').replace("'", ''), 
                     'flag_values': list(vardict[var].values()),
                     'long_name:': lwtvars_ln[count],
                     'units:': lwtvars_un[count]}<jupyter_output><empty_output><jupyter_text>Global metadata  Get columns requirements and BSH entry from excel file<jupyter_code>pdglat    = pd.read_excel(open('Metadata_for_atmodat_standard_BSH_PL_AH.xlsx', 'rb'),
              sheet_name='data_files', header=2, nrows=31, usecols=[1, 8])<jupyter_output><empty_output><jupyter_text>Remove all rows that have NaN entry<jupyter_code>glat      = pdglat[pdglat.iloc[:, 1].str.contains("NaN")==0]<jupyter_output><empty_output><jupyter_text>Create dictionary with all global attributes<jupyter_code>dictglat  = dict(zip(glat.iloc[0:, 0].str.strip(),glat['BSH entry']))
#-- add quotation .apply(lambda x: "'" + str(x) + "'"))) or glat.update(glat[['BSH entry']].applymap('"{}"'.format)) # glat.iloc[1:, 1]))<jupyter_output><empty_output><jupyter_text>Write global attributes into xarray dataset<jupyter_code>ds.attrs = dictglat  #-- select first 2 dictionary elements: dict(list(dictglat.items())[:2])<jupyter_output><empty_output><jupyter_text>Metadata of coordinate variables<jupyter_code>ds['time'].attrs={'standard_name':'time', 'long_name':'Time'}<jupyter_output><empty_output><jupyter_text>Write out data as netCDF<jupyter_code>ds.to_netcdf(f"{folder}/{fnroot}.nc", encoding={**varsenc, **dimsenc})<jupyter_output><empty_output><jupyter_text>Write out annual files<jupyter_code>years, datasets = zip(*ds.groupby("time.year"))
# paths = [f"{fnroot}_{y}.nc" for y in years]; xr.save_mfdataset(datasets, paths)  #--> cannot use encoding!! 
for count, year in enumerate(years):
    datasets[count].to_netcdf(f"{folder}/{fnroot}_{year}.nc", encoding={**varsenc, **dimsenc})
!ncdump -h $folder/lwtns06_1952.nc
!cdo infov $folder/lwtns06_1952.nc<jupyter_output>-1 :       Date     Time   Level Gridsize    Miss :     Minimum        Mean     Maximum : Parameter name
     1 : 1952-01-01 00:00:00       0        1       0 :                  2.0000             : lamb          
     2 : 1952-01-01 00:00:00       0        1       0 :                  1.0000             : flowdir       
     3 : 1952-01-01 00:00:00       0        1       0 :                  10.000             : winddir       
     4 : 1952-01-01 00:00:00       0        1       0 :                  0.0000             : stormsev      
     5 : 1952-01-02 00:00:00       0        1       0 :                  4.0000             : lamb          
     6 : 1952-01-02 00:00:00       0        1       0 :                  1.0000             : flowdir       
     7 : 1952-01-02 00:00:00       0        1       0 :                  10.000             : winddir       
     8 : 1952-01-02 00:00:00       0        1       0 :                  0.0000             : stormsev      
     9 : 1[...]
{'fraction_non_alphanumeric': 0.13343291013193925, 'fraction_numerical': 0.038088125466766244, 'mean_word_length': 3.805622009569378, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 82, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '236521', 'license': 'Apache-2.0', 'n_tokens_mistral': 2968, 'n_tokens_neox': 2739, 'n_words': 593, 'path': 'notebooks/Convert_BSH_ascii_to_ATMODAT-netCDF-v1.ipynb', 'repo_name': 'AtMoDat/atmodat_standard_checker', 'size': '8041'}