【空气质量数据分析专题八】污染物浓度日变化分析_python

前言

对空气质量日级别五年数据进行日变化分析，可以看出污染物浓度逐日变化的特征。

分析流程

对数据进行专题二的预处理后，计算出各污染物全时段的各日平均浓度，最后进行可视化分析。

日变化分析方式有多种，这里通过日历图进行分析。

核心代码

这部分使用Python处理数据，然后使用R进行绘图，其中a_January代表1月，后面的以此类推，在月份英文名前加字母主要是确保月份有序，从而峰峦图才能按月份排列，更为直观。

（1）处理数据（Python代码）

def daily_trend_analysis(self, df_station, year_list):
    """
    日均浓度趋势分析
    :param df_station: 站点日均浓度数据
    :param year_list: 年份列表
    :return:
    """
    result2 = pd.pivot_table(df_station, index=['month', 'day'], aggfunc=np.mean, values=['PM10', 'PM2.5', 'SO2', 'NO2', 'O3', 'CO', 'AQI'])
    for i in result2.index[:]:
        result2.loc[i, 'date'] = datetime(year=2020, month=i[0], day=i[1])
    result = result2.copy()
    result.set_index('date', inplace=True)
    pic_loc0 = Path(self.cf_info['output']['picture']).joinpath(df_station['city'].values[0])
    pic_loc = pic_loc0.joinpath('污染物日变化特征')
    if not os.path.exists(pic_loc):
        os.mkdir(pic_loc)
    result.to_excel(pic_loc / (df_station['station'].values[0] + str(year_list[0]) + '-' + str(year_list[-1]) + '年各污染物浓度日变化.xls'), encoding='gbk')

（2）可视化分析（R代码，省去了路径），这里划分了更加细致的浓度及AQI等级主要是为了使可视化区分度更明显，效果更好

library('xlsx')
library('openair')
Sys.setlocale(category  = "LC_ALL", locale = "C")

micepath <- "数据所在文件夹"
micefiles <- list.files(micepath, full.names = TRUE)

breaks <- c(0,49,99,149,199,249,299)
labels <- c("优","良","轻度污染","中度污染","重度污染","严重污染")
cols=c("green","yellow","orange","red","purple","maroon")

for(i in 1:length(micefiles)){
  file_full <- strsplit(micefiles[i], '/')
  f_names <- strsplit(file_full[[1]][2], '各')
  name <- f_names[[1]][1]
  
  datas <- read.xlsx(micefiles[i], sheetIndex=1)
  all_col_name = names(datas)
  all_col_name_list = strsplit(all_col_name, ' ')
  len <- 2:length(all_col_name_list)
  for(j in len){
    x1 = paste(name, unlist(all_col_name_list[j]), sep = "")
    data <- datas[ , c(unlist(all_col_name_list[1]),unlist(all_col_name_list[j]))]
    if (unlist(all_col_name_list[j])== "O3"){
      jpeg(file=paste(x1,"jpeg",sep="."),units = "px",width = 3*600,height = 3*600,res = 3*72)
      calendarPlot(data, pollutant = "O3",breaks=c(0,20,40,60,80,100,120,140,160,180,200), labels=c("一级Ⅰ","一级Ⅱ","一级Ⅲ","一级Ⅳ","一级Ⅴ","二级Ⅰ","二级Ⅱ","二级Ⅲ","超标Ⅰ","超标Ⅱ"),year=2020,annotate="value",
                   lim=160,layout=c(3,4), cols = "jet",col.lim=c("black","cyan"),key.footer='0μg/m^3',key.header='200μg/m^3')
      dev.off()
    }

    if (unlist(all_col_name_list[j])== "CO"){
      jpeg(file=paste(x1,"jpeg",sep="."),units = "px",width = 3*600,height = 3*600,res = 3*72)
      calendarPlot(data, pollutant = "CO",breaks=c(0,0.3,0.6,0.9,1.2,1.5,4,5), labels=c("达标Ⅰ","达标Ⅱ","达标Ⅲ","达标Ⅳ","达标Ⅴ","达标Ⅵ","超标"),year=2020,annotate="value",
                   lim=4,layout=c(3,4),col.lim=c("black","cyan"),key.footer='0mg/m^3',key.header='5mg/m^3')
      dev.off()
    }

    if (unlist(all_col_name_list[j])== "NO2"){
      jpeg(file=paste(x1,"jpeg",sep="."),units = "px",width = 3*600,height = 3*600,res = 3*72)
      calendarPlot(data, pollutant = "NO2",breaks=c(0,10,20,30,40,50,60,70,80,90,100), labels=c("达标Ⅰ","达标Ⅱ","达标Ⅲ","达标Ⅳ","达标Ⅴ","达标Ⅵ","达标Ⅶ","达标Ⅷ","超标Ⅰ","超标Ⅱ"),year=2020,annotate="value",
                   lim=80,layout=c(3,4),col.lim=c("black","cyan"),key.footer='0μg/m^3',key.header='100μg/m^3')
      dev.off()
    }

    if (unlist(all_col_name_list[j])== "SO2"){
      jpeg(file=paste(x1,"jpeg",sep="."),units = "px",width = 3*600,height = 3*600,res = 3*72)
      calendarPlot(data, pollutant = "SO2",breaks=c(0,25,50,75,100,125,150,175), labels=c("一级Ⅰ","一级Ⅱ","二级Ⅰ","二级Ⅱ","二级Ⅲ","二级Ⅳ","超标"),year=2020,annotate="value",
                   lim=150,layout=c(3,4),col.lim=c("black","cyan"),key.footer='0μg/m^3',key.header='175μg/m^3')
      dev.off()
    }

    if (unlist(all_col_name_list[j])== "PM10"){
      jpeg(file=paste(x1,"jpeg",sep="."),units = "px",width = 3*600,height = 3*600,res = 3*72)
      calendarPlot(data, pollutant = "PM10",breaks=c(0,25,50,75,100,125,150,175,200), labels=c("一级Ⅰ","一级Ⅱ","二级Ⅰ","二级Ⅱ","二级Ⅲ","二级Ⅳ","超标Ⅰ","超标Ⅱ"),year=2020,annotate="value",
                   lim=150,layout=c(3,4),col.lim=c("black","cyan"),key.footer='0μg/m^3',key.header='200μg/m^3')
      dev.off()
    }

    if (unlist(all_col_name_list[j])== "PM2.5"){
      jpeg(file=paste(x1,"jpeg",sep="."),units = "px",width = 3*600,height = 3*600,res = 3*72)
      calendarPlot(data, pollutant = "PM2.5",breaks=c(0,15,25,35,45,55,65,75,85,100), labels=c("一级Ⅰ","一级Ⅱ","一级Ⅲ","二级Ⅰ","二级Ⅱ","二级Ⅲ","二级Ⅳ","超标Ⅰ","超标Ⅱ"),year=2020,annotate="value",
                   lim=75,layout=c(3,4),col.lim=c("black","cyan"),key.footer='0μg/m^3',key.header='100μg/m^3')
      dev.off()
    }

    if (unlist(all_col_name_list[j])== "AQI"){
      jpeg(file=paste(x1,"jpeg",sep="."),units = "px",width = 3*600,height = 3*600,res = 3*72)
      calendarPlot(data, pollutant = "AQI",breaks=c(0,50,100,150,200,250,300), labels=c("优","良","轻度污染","中度污染","重度污染","严重污染"),year=2020,annotate="value",cols=cols,
                   lim=100,layout=c(3,4),col.lim=c("black","cyan"),key.footer='0μg/m^3',key.header='300μg/m^3')
      dev.off()
    }

  }
}