<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Ideas2&nbsp;&nbsp;Libraries3&nbsp;&nbsp;Load data3.1&nbsp;&nbsp;PM2.53.2&nbsp;&nbsp;Observed data from Noibai (ISD data)3.3&nbsp;&nbsp;MERRA-2 SLV group3.4&nbsp;&nbsp;MERRA-2 FLX group4&nbsp;&nbsp;Concatinate data5&nbsp;&nbsp;Correlation and filter out5.1&nbsp;&nbsp;let get wind speed data in different formats5.2&nbsp;&nbsp;Chart6&nbsp;&nbsp;Nudgets are here!  Ideas - select parameters influence concentration of PM2.5 mostly in transport- data should be indenpendent to PM2.5, and should be a physical entities rather a prepresentation of such data- I explicitly select data from MERRA-2 groups (SLV: single level, and some from FLX: flux diagnosis)  Libraries<jupyter_code>import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['figure.figsize'] = (14,6)
import datetime
# plt.rcParamsDefault<jupyter_output><empty_output><jupyter_text>- let visit some graph we hare in the previous exercise - with observed ground data      - they are all basic parameter, let use all of them. For wind speeds, we will have more options with MERRA-2 - with SLV groups in MERRA-2 reanalysis      -  selected columnsT2MWET, T2MDEW, T2M, QV10M, PS, TQV, TQL, H1000, DISPH, wind (2-50, and 1400m) - FLX group select columns:      - QV10M 10-meter_specific_humidity      - QLML surface_specific_humidity      - TQL total_precipitable_liquid_water      - FRCAN areal_fraction_of_anvil_showers      - HLML surface_layer_height      - H1000 height_at_1000_mb      - DISPH zero_plane_displacement_height       - U2M 2-meter_eastward_wind      - V2M 2-meter_northward_wind      - U50M eastward_wind_at_50_meters      - V50M northward_wind_at_50_meters      - U850 eastward_wind_at_850_hPa      - V850 northward_wind_at_850_hPa      - SPEED surface_wind_speed, m s-1      - SPEEDMAX surface_wind_speed, m s-1<jupyter_code># let build a function to look up the name and unit from netcdf4 file
import netCDF4 as nc
def collect_meta_name():
    meta_name = dict()
    files = ['data/nc4/MERRA2_400.tavg1_2d_slv_Nx.20180722.nc4',
             'data/nc4/MERRA2_400.tavg1_2d_flx_Nx.20180102.nc4',
             'data/nc4/MERRA2_400.tavg1_2d_aer_Nx.20180101.nc4']
    for file in files:
        ds = nc.Dataset(file)
        for k, v in ds.variables.items():
            meta_name[k] = f'{v.long_name}, {v.units}'
    return meta_name
name_ = collect_meta_name()
def look_kw(kw=None):
    name_ = collect_meta_name()
    
    for k,v in name_.items():
        if kw in v:
            print(k,v)
    return None
look_kw(kw='wind')
look_kw(kw='plane')
look_kw(kw='height')<jupyter_output>H250 height_at_250_hPa, m
DISPH zero_plane_displacement_height, m
H1000 height_at_1000_mb, m
H850 height_at_850_hPa, m
H500 height_at_500_hPa, m
TCZPBL transcom_planetary_boundary_layer_height, m
HLML surface_layer_height, m
PBLH planetary_boundary_layer_height, m<jupyter_text>Load data  PM2.5<jupyter_code>pm25 = pd.read_csv('data/cleaned_pm25_Hanoi_PM2.5_2018_YTD.csv',
                   parse_dates=['Date (LT)'], 
                   index_col=['Date (LT)'])
pm25.head()<jupyter_output><empty_output><jupyter_text>Observed data from Noibai (ISD data)<jupyter_code>nb = pd.read_csv('data/cleaned_noibai_noaa_isd_2018.csv',
                parse_dates=['DATE'],
                index_col=['DATE'])
nb.head()
 # since it use UTC time, convert to the local time (GMT+7)
nb.index = nb.index + datetime.timedelta(hours=7)<jupyter_output><empty_output><jupyter_text>MERRA-2 SLV group<jupyter_code># T2MWET, T2MDEW, T2M, QV10M, PS, TQV, TQL, H1000, DISPH, winds
cols = ['T2MWET', 'T2MDEW', 'T2M', 'QV10M', 'PS', 'TQV', 'TQL', 'H1000', 'DISPH', 'U2M', 'V2M', 'U50M', 'V50M', 'U850', 'V850']
# read and filter out by the column name
slv = pd.read_csv('data/merra2_slv_hanoi_2018.csv',
                 parse_dates=['time'],
                 index_col=['time'])
print(slv.columns)
slv = slv[cols]
slv.head()<jupyter_output>Index(['U2M', 'V250', 'TROPT', 'TROPPB', 'T2M', 'TQL', 'T500', 'U850', 'PS',
       'V850', 'H250', 'Q250', 'T2MDEW', 'V50M', 'Q500', 'DISPH', 'H1000',
       'TS', 'T10M', 'TROPPT', 'SLP', 'U250', 'Q850', 'ZLCL', 'TQV', 'V2M',
       'T250', 'TROPQ', 'V10M', 'H850', 'T850', 'U50M', 'U10M', 'TROPPV',
       'H500', 'V500', 'T2MWET', 'U500', 'QV10M'],
      dtype='object')<jupyter_text>MERRA-2 FLX group<jupyter_code># - QV10M 10-meter_specific_humidity
# - QLML surface_specific_humidity
# - TQL total_precipitable_liquid_water
# - FRCAN areal_fraction_of_anvil_showers
# - HLML surface_layer_height
# - H1000 height_at_1000_mb
# - DISPH zero_plane_displacement_height
# - HLML surface_layer_height
# - U2M 2-meter_eastward_wind
# - V2M 2-meter_northward_wind
# - U50M eastward_wind_at_50_meters
# - V50M northward_wind_at_50_meters
# - U850 eastward_wind_at_850_hPa
# - V850 northward_wind_at_850_hPa
# - SPEED surface_wind_speed, m s-1
# - SPEEDMAX surface_wind_speed, m s-1
cols = ['QLML','FRCAN', 'HLML',  'SPEED', 'SPEEDMAX', 'RHOA']
# similar to FLX group, select and filter out
flx = pd.read_csv('data/merra2_flx_hanoi_2018.csv',
                 parse_dates=['time'],
                 index_col=['time'])
flx.columns
flx = flx[cols]
flx.head()
flx.index = flx.index.tz_localize(None)<jupyter_output><empty_output><jupyter_text>Concatinate data<jupyter_code># now will combine each dataframe using merge 
# if we have the same column name for the index, say 'time', pd.concat() can combine a list of dataframe
flx.head()
# PM2.5 and SLV group
df = pd.merge(pm25, slv, right_index=True, left_index=True)
# then FLX group to the two previous dfs
df = pd.merge(df, flx, right_index=True, left_index=True, how='left')
# finally data from Noibai Intl Airpot
df = pd.merge(df, nb, right_index=True, left_index=True, how='left')
df.describe()
df.columns<jupyter_output><empty_output><jupyter_text>Correlation and filter out<jupyter_code># we will do several steps to filter out dependent columns (well correlated with existing columns)
fig, ax = plt.subplots(figsize=(7,6))
sns.heatmap(df.corr(), cmap='seismic')
fig.tight_layout()
fig.savefig('img/2020Aug-corr-heatmap.png', dpi=120, optimize=True)
fig, ax = plt.subplots(figsize=(6,10))
df.corr()['pm25'].sort_values().to_frame().drop('pm25').plot.barh(ax=ax)<jupyter_output><empty_output><jupyter_text>let get wind speed data in different formats<jupyter_code># let detour and get wind in speed and direction (degree from 2.4 exercise)
dfw = pd.read_csv('data/merra2_hanoi_2018_wind_converted.csv', 
                  parse_dates=['Unnamed: 0'],
                 index_col=['Unnamed: 0'])
dfw.head()
dfw.columns
cols = ['v_2m', 'd_2m', 'v_50m', 'd_50m', 'v_850', 'd_850']
dfw = dfw[cols]
# and combine all data in one place
df = pd.merge(df, dfw, right_index=True, left_index=True, how='left')
df.columns
fig, ax = plt.subplots(figsize=(6,10))
df.corr()['pm25'].sort_values().to_frame().drop('pm25').plot.barh(ax=ax)
df.index.rename('DATE', inplace=True)
df.rename(columns={'pm25': 'PM2.5'}, inplace=True)
# save it, and we can use it later one (if needed)
df.to_csv('data/comb_PM25_Hanoi_2018.csv')<jupyter_output><empty_output><jupyter_text>Chart<jupyter_code># and recreate the figure above with the standard name
plt.style.use('seaborn-whitegrid')
fig, ax = plt.subplots(figsize=(8,10))
df.corr()['PM2.5'].sort_values().to_frame().dropna().drop('PM2.5').plot.barh(ax=ax)
ax.legend(['$PM_{2.5}$'], frameon=True)
labels = [item.get_text() for item in ax.get_yticklabels()]

# looking for a standard name for each abbreviation
new_label = dict()
for label in labels:
    if label in list(name_.keys()):
        new_label[label] = name_[label]
    else:
        new_label[label] = label
ax.set_yticklabels(new_label.values())
plt.title('Correlation of $PM_{2.5}$ with selected data\nfor Hanoi 2018', loc='right')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>- too many, and some duplicated data, let define a list of columns to be dropped<jupyter_code>df.drop(columns=['U2M', 'V2M', 'U50M', 'V50M', 'U850', 'V850'], inplace=True)
# save data as always
df.to_csv('data/comb_PM25_wind_Hanoi_2018.csv')
df.columns
# and recreate the figure above with the standard name
fig, ax = plt.subplots(figsize=(8,8))
df.corr()['PM2.5'].sort_values().to_frame().dropna().drop('PM2.5').plot.barh(ax=ax)
ax.legend(['$PM_{2.5}$'], frameon=True)
labels = [item.get_text() for item in ax.get_yticklabels()]

# looking for a standard name for each abbreviation
new_label = dict()
for label in labels:
    if label in list(name_.keys()):
        new_label[label] = name_[label]
    else:
        new_label[label] = label
ax.set_yticklabels(new_label.values())
ax.set_title('Correlation of $PM_{2.5}$ with selected data\nfor Hanoi 2018', 
          y=1.05, fontsize=15, weight='bold')
fig.subplots_adjust(bottom=0.1)
fig.text(0.9,.04,s='Datasource:\n airnow.gov ($PM_{2.5}$)\nncei.noaa.gov (ISD)\nnasa.gov (MERRA-2)', 
         ha='right', color='gray')
plt.tight_layout(rect=(0, 0.1, 1,1))
# plt.tight_layout()
plt.savefig('img/2020Aug-PM25-allin.png', dpi=120, optimize=True)
# less inputs
fig, ax = plt.subplots(figsize=(6,10))
df.corr()['PM2.5'].sort_values().to_frame().drop('PM2.5').plot.barh(ax=ax)
# some test to make sure we know they these are highly-dependent, 
df.corr()['T2M']['TMP']
df.corr()['WS']['v_2m']
# some cleaning
df.loc[df.WS>20, 'WS'] = None
# not always as expected, one from observed station, another from a ground observed stations
plt.plot(df.index, df.WS)
plt.plot(df.index, df.v_2m)
# ground station reported with coarse resolution, 0.5m/s increment
plt.scatter(df.WS, df.v_2m)
df.corr()['T2MDEW']['T2MWET']
plt.plot(df.index, df.CIG)
df.corr()['SPEED'].sort_values()
df.corr()['DEW']['T2MDEW']
df.columns
cols_drop = ['T2MWET', 'SPEED', 'SPEEDMAX', 'DEW', 'TMP']
df.columns
df.drop(columns=cols_drop, inplace=True)
df.to_csv('data/comb_PM25_wind_Hanoi_2018_v1.csv')
# let see how the heatmap again
sns.heatmap(df.corr(), cmap='seismic')
# and correlation only with PM2.5
fig, ax = plt.subplots(figsize=(6,10))
df.corr()['PM2.5'].sort_values().to_frame().drop('PM2.5').plot.barh(ax=ax)
colsd = ['QV10M', 'CLDHT', 'QLML']
# Height to the lowest cloud is very weak, let see data before we drop it
plt.plot(df.index, df.CLDHT)
df.drop(columns=colsd, inplace=True)
fig, ax = plt.subplots(figsize=(6,10))
df.corr()['PM2.5'].sort_values().to_frame().drop('PM2.5').plot.barh(ax=ax)
df.to_csv('data/comb_PM25_wind_Hanoi_2018_v2.csv')
more_cols = ['d_850', 'd_50m', 'd_2m', 'VIS', 'WD']
df.drop(columns=more_cols, inplace=True)
fig, ax = plt.subplots(figsize=(6,10))
df.corr()['PM2.5'].sort_values().to_frame().drop('PM2.5').plot.barh(ax=ax)
df['T2MDEW'] = df['T2MDEW'] - 273.15
df['T2M'] = df['T2M'] - 273.15
df.to_csv('data/comb_PM25_wind_Hanoi_2018_v3.csv')
df.describe()
# much smaller dataset now
fig, ax = plt.subplots(figsize=(7,6))
sns.heatmap(df.corr(), cmap='seismic')
# let add name to nameplace holder, 
name_['CIG'] = 'ceiling height dimension, m'
name_['CLDCR'] = 'cloud cover, 1'
name_['v_2m'] = 'wind speed at 2m, m/s'
name_['v_50m'] = 'wind speed at 50m, m/s'
name_['v_850'] = 'wind speed at 850hPa (~1450m)'
name_['WS'] = 'observed ground wind speed, m/s'<jupyter_output><empty_output><jupyter_text>Nudgets are here!<jupyter_code># and recreate the figure above with the standard name
fig, ax = plt.subplots(figsize=(8,6))
df.corr()['PM2.5'].sort_values().to_frame().dropna().drop('PM2.5').plot.barh(ax=ax)
ax.legend(['$PM_{2.5}$'], frameon=True)
labels = [item.get_text() for item in ax.get_yticklabels()]

# looking for a standard name for each abbreviation
new_label = dict()
for label in labels:
    if label in list(name_.keys()):
        new_label[label] = name_[label]
    else:
        new_label[label] = label
ax.set_yticklabels(new_label.values())
ax.set_title('Correlation of $PM_{2.5}$ with selected data\nfor Hanoi 2018', 
          y=1.05, fontsize=15, weight='bold')
fig.subplots_adjust(bottom=0.1)
fig.text(0.7,.03,s='Datasource:\nairnow.gov ($PM_{2.5}$)\nncei.noaa.gov (ISD)\nnasa.gov (MERRA-2)', 
         ha='left', color='gray')
plt.tight_layout(rect=(0, 0.15, 1,1))
# plt.tight_layout()
plt.savefig('img/2020Aug-PM25-selected.png', dpi=120, optimize=True)<jupyter_output><empty_output><jupyter_text>- OK, now we are ready to rock some prediction<jupyter_code>df.shape
new_label<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15426687042548057, 'fraction_numerical': 0.04890211533821282, 'mean_word_length': 5.003862868179623, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '170197', 'license': 'MIT', 'n_tokens_mistral': 5270, 'n_tokens_neox': 4850, 'n_words': 1028, 'path': '3.1 Data Selection.ipynb', 'repo_name': 'binh-bk/air-quality-analysis', 'size': '12433'}