污染观测数据处理说明

质量控制+转换成ppb

• Concentrations were reported by the MEE in micrograms per cubic meter (µg m 3) under standard conditions (273 K, 1013 hPa) until 31 August 2018.
• This reference state was changed on 1 September 2018 to 298 K and 1013 hPa for gases and to local ambient state for PM2.5

from scipy.stats import zscore

def qualitycontrol(x):
    temp = zscore(x,nan_policy='omit')
    x1 = temp[0:-2]
    x2 = temp[1:-1]
    x3 = temp[2:]
    RM3 = 3*x2/(x1+x2+x3)
    mask1 = abs(temp>4)
    mask2 = RM3>2
    mask1[1:-1] = (mask2|mask1[1:-1])
    return np.where(mask1,np.nan,x)

• 单位转换

  def convert_unit(x,T,P,M):
      # Convert from ug/m3 to ppb
      return x*22.41*T/273.*1013/P/M
  
  T, P = 273, 1013
  
  for year in [2013,2014,2015,2016,2017]:
      inpdir  = f"处理后数据/NoQC/{year}"
      outdir1 = f"处理后数据/QC+PPB/hour/{year}"
      outdir2 = f"处理后数据/QC+PPB/day/{year}"
      if( not os.path.exists(outdir1)):
          os.makedirs(outdir1)
      if( not os.path.exists(outdir2)):
          os.makedirs(outdir2)
  
      for ProvinceCode,CityCode,StationCode in zip(allsite.ProvinceCode, allsite.CityCode, allsite.StationCode):
          fname = f"{inpdir}/{ProvinceCode}_{CityCode}_{StationCode}_{year}.csv"
          if(os.path.exists(fname)):
              data = pd.read_csv(fname)
              data.columns = cols_order
              data['o3']  = convert_unit(data['o3'],T,P,48)
              data['co']  = convert_unit(data['co'],T,P,28)
              data['so2'] = convert_unit(data['so2'],T,P,64)
              data['no2'] = convert_unit(data['no2'],T,P,46)
              data[cols_order].to_csv(f"{outdir1}/{ProvinceCode}_{CityCode}_{StationCode}_{year}.csv",index=None)
  
              data.index = pd.to_datetime(data.datetime,format="%Y-%m-%d %H:%M:%S")
              dfday = data[['pm25','pm10','o3','co','so2','no2']].resample('D').apply(strongmean)
              dfday['MDA8'] = data["o3"].resample('D').apply(MDA8)
              dfday['MDA8'] = np.where(np.isnan(dfday['o3']),np.nan,dfday['MDA8'])
              
              dfday['MaxO3'] = data["o3"].resample('D').apply(strongmax)
              dfday[cols_order_day].to_csv(f"{outdir2}/{ProvinceCode}_{CityCode}_{StationCode}_{year}.csv")
          else:
              print("%s doesn't exist" %(fname))
             
  year = 2018
  T1, P1 = 273, 1013
  T2, P2 = 298, 1013
  
  inpdir  = f"处理后数据/NoQC/{year}"
  outdir1 = f"处理后数据/QC+PPB/hour/{year}"
  outdir2 = f"处理后数据/QC+PPB/day/{year}"
  if( not os.path.exists(outdir1)):
      os.makedirs(outdir1)
  if( not os.path.exists(outdir2)):
      os.makedirs(outdir2)
  
  for ProvinceCode,CityCode,StationCode in zip(allsite.ProvinceCode, allsite.CityCode, allsite.StationCode):
      fname = f"{inpdir}/{ProvinceCode}_{CityCode}_{StationCode}_{year}.csv"
      if(os.path.exists(fname)):
          data = pd.read_csv(fname)
          data.columns = cols_order
          data.index = pd.to_datetime(data.datetime,format="%Y-%m-%d %H:%M:%S")
          
          data.loc['2018-01-01':'2018-08-31','o3']  = convert_unit(data.loc['2018-01-01':'2018-08-31','o3'] ,T1,P1,48)
          data.loc['2018-01-01':'2018-08-31','co']  = convert_unit(data.loc['2018-01-01':'2018-08-31','co'] ,T1,P1,28)
          data.loc['2018-01-01':'2018-08-31','so2'] = convert_unit(data.loc['2018-01-01':'2018-08-31','so2'],T1,P1,64)
          data.loc['2018-01-01':'2018-08-31','no2'] = convert_unit(data.loc['2018-01-01':'2018-08-31','no2'],T1,P1,46)
          
          data.loc['2018-09-01':'2018-12-31','o3']  = convert_unit(data.loc['2018-09-01':'2018-12-31','o3'] ,T2,P2,48)
          data.loc['2018-09-01':'2018-12-31','co']  = convert_unit(data.loc['2018-09-01':'2018-12-31','co'] ,T2,P2,28)
          data.loc['2018-09-01':'2018-12-31','so2'] = convert_unit(data.loc['2018-09-01':'2018-12-31','so2'],T2,P2,64)
          data.loc['2018-09-01':'2018-12-31','no2'] = convert_unit(data.loc['2018-09-01':'2018-12-31','no2'],T2,P2,46)
          
          data[cols_order].to_csv(f"{outdir1}/{ProvinceCode}_{CityCode}_{StationCode}_{year}.csv",index=None)
  
          dfday = data[['pm25','pm10','o3','co','so2','no2']].resample('D').apply(strongmean)
          dfday['MDA8'] = data["o3"].resample('D').apply(MDA8)
          dfday['MDA8'] = np.where(np.isnan(dfday['o3']),np.nan,dfday['MDA8'])
          
          dfday['MaxO3'] = data["o3"].resample('D').apply(strongmax)
          dfday[cols_order_day].to_csv(f"{outdir2}/{ProvinceCode}_{CityCode}_{StationCode}_{year}.csv")
      else:
          print("%s doesn't exist" %(fname))

-