<jupyter_start><jupyter_text>** 数据最全的站点 **<jupyter_code>data_dir = './hourly_data'
# 中国环境监测总站, 数据较全
jczz = os.path.join(data_dir, '110000012.txt')

data = pd.read_csv(
    jczz, 
    sep=',', 
    names=[
        'time','stationcode','longitude','latitude',
        'PM10','PM25',
        'Cl-','NO3-','SO42-','Ca2+','Na+','K+','NH4+','Mg2+',
        'OC','EC',
        'Hg','Br','As','Si','Se','Te','V','Ca','Ti','Ba','Sc',
        'Pd','Co','Mo','K','Fe','Pb','TI','Cu','Cr','Cs','Ag',
        'Zn','Sb','Sn','Mn','Cd','Ni','Ga'],
    encoding='utf-8')
# -999.0 表示无效记录，但至少有记录表示，有些时间段是完全没有记录的
data.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 8760 entries, 0 to 8759
Data columns (total 45 columns):
 #   Column       Non-Null Count  Dtype  
---  ------       --------------  -----  
 0   time         8760 non-null   object 
 1   stationcode  8760 non-null   int64  
 2   longitude    8760 non-null   float64
 3   latitude     8760 non-null   float64
 4   PM10         8760 non-null   float64
 5   PM25         8760 non-null   float64
 6   Cl-          8760 non-null   float64
 7   NO3-         8760 non-null   float64
 8   SO42-        8760 non-null   float64
 9   Ca2+         8760 non-null   float64
 10  Na+          8760 non-null   float64
 11  K+           8760 non-null   float64
 12  NH4+         8760 non-null   float64
 13  Mg2+         8760 non-null   float64
 14  OC           8760 non-null   float64
 15  EC           8760 non-null   float64
 16  Hg           8760 non-null   float64
 17  Br           8760 non-null   float64
 18  As           8760 non-null   float64
 19  Si    [...]<jupyter_text>+ 先进行组分分类，分别关联PM10，PM25以便后期分析+ 第一步预处理，替换无效值+ 处理时间+ 单位转化<jupyter_code>obs_metal = data[['time', 'PM10', 'PM25', 'Hg','Br','As','Si','Se','Te','V','Ca','Ti','Ba','Sc',
    'Pd','Co','Mo','K','Fe','Pb','TI','Cu','Cr','Cs','Ag', 'Zn','Sb','Sn','Mn','Cd','Ni','Ga']]
obs_metal.replace(-999.0, np.nan, inplace=True)
obs_metal['time'] = obs_metal['time'].apply(lambda x: datetime.strptime(x, "%Y-%m-%d:%H"))
obs_metal[obs_metal.columns[3:]] = obs_metal[obs_metal.columns[3:]].apply(lambda x: x/1000, axis=1)
obs_metal.sample(3)
obs_metal.info()
# 完整的数据只有16组观测
obs_metal.dropna().shape
plt.figure(figsize=(20, 20))
sns.heatmap(obs_metal.corr(), annot=True)
plt.title('重金属与污染物相关性热力图', fontsize=16)
plt.savefig('../DataVisual/images/重金属与污染物相关性热力图.png')
plt.show()<jupyter_output><empty_output><jupyter_text>+ 看下数据的统计分布<jupyter_code>obs_metal[['PM25', 'Ca', 'Si', 'Fe']].describe(percentiles=[.08, .09, .1, .7, .8, .9, .95, .99])
obs_metal['Ca'].plot()
obs_metal['Si'].plot()
obs_metal['Fe'].plot()
obs_metal[[
    'PM25','Hg','Br','As','Se','Te','V','Ti','Ba','Sc','Pd','Co','Mo','K','Pb',
    'TI','Cu','Cr','Cs','Ag', 'Zn','Sb','Sn','Mn','Cd','Ni','Ga']].describe(percentiles=[.05, .7, .8, .9, .95, .99])
# 计算百分位数

obs_metal['Ca'].quantile(q=0.95), obs_metal['Te'].quantile(q=.05)
# 其他元素 99%分位数 大于 1up/m3 的情况分析
obs_metal['Sn'].plot()
# 其他元素 99%分位数 大于 1up/m3 的情况分析
obs_metal[(obs_metal['Sn']>0)&(obs_metal['Sn']<1)]['Sn'].plot()
fig, ax = plt.subplots(figsize=(6, 6))

ax.pie(
    x=obs_metal.describe()[1:2][obs_metal.columns[3:]], 
    labels=obs_metal.describe().columns[2:], 
    labeldistance=1.01)
ax.set_title('在线重金属 均值 分布', fontsize=30)
ax.legend(loc='upper left', fontsize=15)
# plt.savefig('../DataVisual/images/在线重金属 均值 分布.png')
plt.show()
obs_metal.describe()[1:2][obs_metal.columns[3:]]<jupyter_output><empty_output><jupyter_text>+ 当所有元素浓度加合大于 PM2.5浓度时，则该小时所有项目监测结果判定为异常数据。---+ 当单一项目监测结果突然出现单点异常高值或低值（离群，脱离总体趋势），则该数据异常。---+ 正常情况下，Ca 和 Si 监测值范围为 0.1-10μg/m3，Al 不大于 5μg/m3，Fe 不大于 2μg/m3，---+ 其它元素一般不大于 1μg/m3。重污染天气或沙尘天气可能会大于各项目限定值，优良天气可能会小于各项目限定值，具体需要结合颗粒物变化趋势判断是否异常。<jupyter_code>for _, x in enumerate(['Ca', 'Si', 'Fe']):
    print("{} not na : {}".format(x, obs_metal[x].notna().sum()))
obs_metal[obs_metal[['Fe']]>=0].shape, obs_metal[obs_metal[['Fe']]<2].shape
obs_metal[obs_metal['Ca']<=10].shape, obs_metal[obs_metal['Ca']>10].shape
obs_metal['total'] = obs_metal[obs_metal.columns[3:]].apply(lambda x: x.sum(), axis=1)
obs_metal.shape
obs_metal.dropna().shape
# total < PM25
obs_metal[obs_metal['PM25']<obs_metal['total']].shape
obs_metal[obs_metal['PM25']<obs_metal['total']].sample(3)
obs_metal[obs_metal.columns].isna().sum().plot(kind='bar')
obs_metal[['PM25','total']].corr()
# miss obs
obs_metal1 = obs_metal[obs_metal.isnull().values==True].drop_duplicates()
obs_metal1.shape
def plot_line(df):
    fig, ax = plt.subplots(figsize=(20, 6))
    ax.plot(df['time'], df['total'], color='g', label='total(μg/m3)')
    ax.plot(df['time'], df['PM25'], color='b', label='PM25(μg/m3)')
    
    ax.set_xlabel('时间', fontsize=15)
    ax.set_ylabel('浓度', fontsize=15)
    ax.set_title('PM25 vs. Total 逐小时变化', fontsize=18)

    ax.xaxis.set_tick_params(rotation=50, labelsize=12)
    xstart, xend = ax.get_xlim()
    ax.xaxis.set_ticks(np.arange(xstart, xend, 10))
    ax.yaxis.set_tick_params(labelsize=14)
    ax.legend(fontsize=15)
    plt.show()
plot_line(obs_metal)
sns.pairplot(obs_metal[['PM10', 'PM25', 'Ca', 'Si', 'Fe']])
sns.heatmap(obs_metal[['PM10', 'PM25', 'Ca', 'Si', 'Fe']].corr(), annot=True)
# patay 方式建立线性模型
import past
from patsy import dmatrices

y, X = dmatrices('y ~ x', obs_metal[['Ca', 'Si']].rename(
    columns={'Ca': 'x', 'Si': 'y'}))
# Patsy对象可以直接传递到算法（比如numpy.linalg.lstsq）中，它执行普通最小二乘回归
coef, resid, _, _ = np.linalg.lstsq(X, y, rcond=-1)
# 模型的元数据保留在design_info属性中，因此你可以重新附加列名到拟合系数
coef = pd.Series(coef.squeeze(), index=X.design_info.column_names)
coef
# 可视化结果
fig, ax = plt.subplots(figsize=(8, 8))

ax.scatter(obs_metal['Ca'], obs_metal['Si'], color='g', marker='.')
ax.plot(
    np.arange(0, obs_metal['Ca'].max()), 
    coef.x*np.arange(0, obs_metal['Ca'].max())+coef.Intercept,
    color='r', label='拟合曲线\nCoef:{}\nInterp:{}'.format(round(coef.x,4), round(coef.Intercept,4)))
ax.set_title('Si vs. Ca 回归分析', fontsize=15)
ax.set_xlabel('Ca', fontsize=12)
ax.set_ylabel('Si', fontsize=12)
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2238285144566301, 'fraction_numerical': 0.058325024925224324, 'mean_word_length': 4.461887477313975, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '663039', 'license': 'MIT', 'n_tokens_mistral': 3101, 'n_tokens_neox': 2869, 'n_words': 381, 'path': 'analysis/analysis_metal.ipynb', 'repo_name': 'JoeBuzh/WeatherMining', 'size': '6018'}