"""This program was used to calibrate the CAM-chem PM2.5 results with TWEAP station data
usage: python YYYY """
import numpy as np
from pandas import *
from PseudoNetCDF.camxfiles.Memmaps import uamiv
import netCDF4
import matplotlib.pyplot as plt
from scipy.io import FortranFile
import sys,os,datetime
# input the YYYY
iyr=int(sys.argv[1])

#model pm results, bundled by grd04_yd.py, or the following "except" block
fnameO = 'PMf13_12_124_137_83.bin'
year=[i for i in range(2007,2020)]
try:
  with FortranFile(fnameO, 'r') as f:
    pmm = f.read_record(dtype=np.float64)  
except:
  pmm=np.zeros(shape=(max(year)-min(year)+1,12,124,137, 83))-1
  for y in year:
    yr=str(y)
    for m in range(12):
      mm='{:02d}'.format(m+1)
      try:
        grd04=uamiv('/nas1/CAM-chem/'+yr+'/'+yr+mm+'/output/'+yr[2:4]+mm+'IC.S.grd04L','r')
        nt,nlay,nrow,ncol=grd04.variables['PM25'].shape
        pmm[y-min(year),m,:nt,:,:]=grd04.variables['PM25'][:,0,:,:]
      except:
        continue
  with FortranFile(fnameO, 'w') as f:
    f.write_record(pmm)
# observation PM2.5 (21 year,13 month, 24 hours, 608 stations
# the bin file was bundled by dfpm.py
fnameO = 'PMf21_13_32_24_608.bin'
with FortranFile(fnameO, 'r') as f:
  pmf = f.read_record(dtype=np.float64)

# the township id and fractions(CMAQ landuse formatted, generated by mk_townNew.py)
nc = netCDF4.Dataset('20160101.ncT','r')
V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
nt,nlay,nrow,ncol=nc.variables[V[3][0]].shape

# the township name, ID and stations
df_twnaq=read_csv('town_aqstEnew.csv')
df_ll=read_csv('sta_ll.csv')
intv=1
delH=6
AQD={'Northern':[1,11,12,31,32,33,35],'Central':[17,22,36,37,38,39,40],'Southern':[2,21,41,42,43]}
NCS=[i for i in AQD]
pmm=pmm.reshape(max(year)-min(year)+1,12, 124, 137, 83)
pmf=pmf.reshape(2020-2000+1,13,32,24,608)
# observation and simulation values for each township
obs=np.ones(shape=(2020-2000+1,12,124,372))
sim=np.ones(shape=(2020-2007+1,12,124,372))
# dict between town code and AQ station number
si={}
for ii in range(len(df_twnaq)):
  i=df_twnaq.loc[ii,'new_code']
  if i==0:continue
  a=df_twnaq.loc[ii,'aq_st'].split(';')
  if len(a)==0:continue
  s=set(a)
  if len(s)==0 or s==set(['0']):continue
  si.update({i:[int(j) for j in s if len(j)>0]})

# dict between town code and AQ station X and Y index
sx,sy={},{}
for ii in range(len(df_twnaq)):
  i=df_twnaq.loc[ii,'new_code']
  if i==0:continue
  ixx,iyy=[],[]
  for jj in si[i]:
    if jj not in list(df_ll.ID):continue
    lx=list(df_ll.loc[df_ll.ID==jj,'lcp_x'])[0]
    ly=list(df_ll.loc[df_ll.ID==jj,'lcp_y'])[0]
    ix=int((lx-nc.XORIG)//nc.XCELL)
    iy=int((ly-nc.YORIG)//nc.YCELL)
    if ix*(ix-ncol+1)>=0 or iy*(iy-nrow+1)>=0:continue
    ixx.append(ix)
    iyy.append(iy)
  sx.update({i:ixx})
  sy.update({i:iyy})
for y in [iyr]:#year:
  yy=y-min(year) #UTC
  yu=y-2000 #UTC
  for m in range(12):
    bdate=datetime.datetime(y,m+1,1)  
    for t in range(124): 
      tdate=bdate+datetime.timedelta(days=(t*delH+8)/24)
      if tdate>datetime.datetime(2019,6,30):continue
      tdate=tdate.strftime("%Y%m%d%H")
      yr,mn,da,hr=int(tdate[:4]),int(tdate[4:6]),int(tdate[6:8]),int(tdate[8:])    
      yo=yr-2000 #LST
      for ii in range(len(df_twnaq)):
        i=df_twnaq.loc[ii,'new_code']
        if i==0:continue
        if si[i]==['0'] or len(si[i])==0:continue
        try:
          obs[yu,m,t,ii]=np.average(pmf[yo,mn,da,hr,si[i]],weights=(pmf[yo,mn,da,hr,si[i]]>0))
        except:
          obs[yu,m,t,ii]=-1
        try:
          sim[yy,m,t,ii]=np.average(pmm[yy,m,t,sy[i],sx[i]],weights=(pmm[yy,m,t,sy[i],sx[i]]>0))
        except:
          sim[yy,m,t,ii]=-1
  avg=np.average(obs[yu,:,:,:],weights=(obs[yu,:,:,:]>0))
  idx=np.where(obs[yu,:,:,:]<=0)
  obs[yu,idx[0],idx[1],idx[2]]=avg
  avg=np.average(sim[yy,:,:,:],weights=(sim[yy,:,:,:]>0))
  idx=np.where(sim[yy,:,:,:]<=0)
  sim[yy,idx[0],idx[1],idx[2]]=avg

xlabels=[' ']
for y in year:
  if y%2==0:continue
  for d in AQD:
    if d[0]=='N':
      xlabels.append('{:02d}'.format(y%100)+'_'+d[0])
    else:
      xlabels.append(d[0])
fig, ax = plt.subplots()
plt.title(" PM2.5 of CAM-chem Simulations(w/t correction)", loc='center' )
plt.xlabel('Year(2digits)')
plt.ylabel('PM2.5 (microgram/cubic meter)')
data=[]
y_code,s_code,m_code,d_code,h_code,v_code=[],[],[],[],[],[]
for y in [iyr]: #year[:-1]:
#  if y%2==0:continue
  yy=y-min(year) #UTC
  yu=y-2000 #UTC
  for d in AQD:
#['Southern']:
    pm=[]
    for i in AQD[d]:
      dfa=df_twnaq.loc[df_twnaq.code1==i].reset_index(drop=True)
      for ii in dfa.new_code:
        if ii==0:continue
        dct='T'+str(ii)
        if dct not in V[3]:continue
        pp=pmm[yy,:,:,:,:]*nc.variables[dct][0,0,:,:]#).flatten()
        iseq=list(df_twnaq.new_code).index(ii)
        if np.sum(sim[yy,:,:,iseq])<=0:continue #sys.exit('sim<=0 '+dct)
        rat=obs[yu,:,:,iseq]/sim[yy,:,:,iseq]
#        if np.max(rat) >10:sys.exit()
        pr=np.zeros(shape=pp.shape)
        for m in range(12):
          for h in range(124):
#        pr=pp[idx]*rat[idx]
            pr[m,h,:,:]=pp[m,h,:,:]*rat[m,h]
        idx=np.where(pr>0)
        for m in range(12):
          for h in range(124):
            s_idx=[(idx[0][j],idx[1][j],idx[2][j],idx[3][j]) for j in range(len(idx[0])) if (idx[0][j],idx[1][j])==(m,h)]
            if len(s_idx)==0:continue
            y_code.append(y)
            s_code.append(ii)
            m_code.append(m+1)
            d_code.append(h//4+1)
            h_code.append(h%4*6)
            v_code.append(np.mean([pr[s_idx[k]] for k in range(len(s_idx))]))
        pm=pm+list(pr[idx].flatten())
    data.append(pm)
df=DataFrame({'y':y_code,'s':s_code,'m':m_code,'d':d_code,'h':h_code,'v':v_code})
for s in 'ymdh':
  exec('df.'+s+'=[int(i) for i in df.'+s+']')
ymd=[-1 for i in range(len(df))]
for i in range(len(df)):
  y,m,d,h=df.y[i],df.m[i],df.d[i],df.h[i]
  try:
    ymd[i]=int((datetime.datetime(y,m,d,h)+datetime.timedelta(days=8/24.)).strftime('%Y%m%d'))
  except:
    continue
df['ymd']=ymd
df=df.drop(df.loc[df.ymd<0].index).reset_index(drop=True)
dfv=pivot_table(df,index=['ymd','s'],values='v', aggfunc=np.mean).reset_index()
name={i:j for i,j in zip(df_twnaq.new_code,df_twnaq.Name)}
dfv['name']=[name[i] for i in dfv.s]
dfv.set_index('ymd').to_csv('ymd_s_v'+str(iyr)[2:4]+'.csv')
ax.boxplot(data, showfliers=False)
xticks=list(range(0,len(data)+1,intv))
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels, fontsize=8)
plt.savefig('box_AQC.png')


"""This program was used to calibrate the CAM-chem PM2.5 results with TWEAP station data
usage: python YYYY """
import numpy as np
from pandas import *
from PseudoNetCDF.camxfiles.Memmaps import uamiv
import netCDF4
import matplotlib.pyplot as plt
from scipy.io import FortranFile
import sys,os,datetime
# input the YYYY
iyr=int(sys.argv[1])

#model pm results, bundled by grd04_yd.py, or the following "except" block
fnameO = 'PMf13_12_124_137_83.bin'
year=[i for i in range(2007,2020)]
try:
  with FortranFile(fnameO, 'r') as f:
    pmm = f.read_record(dtype=np.float64)  
except:
  pmm=np.zeros(shape=(max(year)-min(year)+1,12,124,137, 83))-1
  for y in year:
    yr=str(y)
    for m in range(12):
      mm='{:02d}'.format(m+1)
      try:
        grd04=uamiv('/nas1/CAM-chem/'+yr+'/'+yr+mm+'/output/'+yr[2:4]+mm+'IC.S.grd04L','r')
        nt,nlay,nrow,ncol=grd04.variables['PM25'].shape
        pmm[y-min(year),m,:nt,:,:]=grd04.variables['PM25'][:,0,:,:]
      except:
        continue
  with FortranFile(fnameO, 'w') as f:
    f.write_record(pmm)
# observation PM2.5 (21 year,13 month, 24 hours, 608 stations
# the bin file was bundled by dfpm.py
fnameO = 'PMf21_13_32_24_608.bin'
with FortranFile(fnameO, 'r') as f:
  pmf = f.read_record(dtype=np.float64)

# the township id and fractions(CMAQ landuse formatted, generated by mk_townNew.py)
nc = netCDF4.Dataset('20160101.ncT','r')
V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]
nt,nlay,nrow,ncol=nc.variables[V[3][0]].shape

# the township name, ID and stations
df_twnaq=read_csv('town_aqstEnew.csv')
df_ll=read_csv('sta_ll.csv')
intv=1
delH=6
AQD={'Northern':[1,11,12,31,32,33,35],'Central':[17,22,36,37,38,39,40],'Southern':[2,21,41,42,43]}
NCS=[i for i in AQD]
pmm=pmm.reshape(max(year)-min(year)+1,12, 124, 137, 83)
pmf=pmf.reshape(2020-2000+1,13,32,24,608)
# observation and simulation values for each township
obs=np.ones(shape=(2020-2000+1,12,124,372))
sim=np.ones(shape=(2020-2007+1,12,124,372))
# dict between town code and AQ station number
si={}
for ii in range(len(df_twnaq)):
  i=df_twnaq.loc[ii,'new_code']
  if i==0:continue
  a=df_twnaq.loc[ii,'aq_st'].split(';')
  if len(a)==0:continue
  s=set(a)
  if len(s)==0 or s==set(['0']):continue
  si.update({i:[int(j) for j in s if len(j)>0]})

# dict between town code and AQ station X and Y index
sx,sy={},{}
for ii in range(len(df_twnaq)):
  i=df_twnaq.loc[ii,'new_code']
  if i==0:continue
  ixx,iyy=[],[]
  for jj in si[i]:
    if jj not in list(df_ll.ID):continue
    lx=list(df_ll.loc[df_ll.ID==jj,'lcp_x'])[0]
    ly=list(df_ll.loc[df_ll.ID==jj,'lcp_y'])[0]
    ix=int((lx-nc.XORIG)//nc.XCELL)
    iy=int((ly-nc.YORIG)//nc.YCELL)
    if ix*(ix-ncol+1)>=0 or iy*(iy-nrow+1)>=0:continue
    ixx.append(ix)
    iyy.append(iy)
  sx.update({i:ixx})
  sy.update({i:iyy})
for y in [iyr]:#year:
  yy=y-min(year) #UTC
  yu=y-2000 #UTC
  for m in range(12):
    bdate=datetime.datetime(y,m+1,1)  
    for t in range(124): 
      tdate=bdate+datetime.timedelta(days=(t*delH+8)/24)
      if tdate>datetime.datetime(2019,6,30):continue
      tdate=tdate.strftime("%Y%m%d%H")
      yr,mn,da,hr=int(tdate[:4]),int(tdate[4:6]),int(tdate[6:8]),int(tdate[8:])    
      yo=yr-2000 #LST
      for ii in range(len(df_twnaq)):
        i=df_twnaq.loc[ii,'new_code']
        if i==0:continue
        if si[i]==['0'] or len(si[i])==0:continue
        try:
          obs[yu,m,t,ii]=np.average(pmf[yo,mn,da,hr,si[i]],weights=(pmf[yo,mn,da,hr,si[i]]>0))
        except:
          obs[yu,m,t,ii]=-1
        try:
          sim[yy,m,t,ii]=np.average(pmm[yy,m,t,sy[i],sx[i]],weights=(pmm[yy,m,t,sy[i],sx[i]]>0))
        except:
          sim[yy,m,t,ii]=-1
  avg=np.average(obs[yu,:,:,:],weights=(obs[yu,:,:,:]>0))
  idx=np.where(obs[yu,:,:,:]<=0)
  obs[yu,idx[0],idx[1],idx[2]]=avg
  avg=np.average(sim[yy,:,:,:],weights=(sim[yy,:,:,:]>0))
  idx=np.where(sim[yy,:,:,:]<=0)
  sim[yy,idx[0],idx[1],idx[2]]=avg

xlabels=[' ']
for y in year:
  if y%2==0:continue
  for d in AQD:
    if d[0]=='N':
      xlabels.append('{:02d}'.format(y%100)+'_'+d[0])
    else:
      xlabels.append(d[0])
fig, ax = plt.subplots()
plt.title(" PM2.5 of CAM-chem Simulations(w/t correction)", loc='center' )
plt.xlabel('Year(2digits)')
plt.ylabel('PM2.5 (microgram/cubic meter)')
data=[]
y_code,s_code,m_code,d_code,h_code,v_code=[],[],[],[],[],[]
for y in [iyr]: #year[:-1]:
#  if y%2==0:continue
  yy=y-min(year) #UTC
  yu=y-2000 #UTC
  for d in AQD:
#['Southern']:
    pm=[]
    for i in AQD[d]:
      dfa=df_twnaq.loc[df_twnaq.code1==i].reset_index(drop=True)
      for ii in dfa.new_code:
        if ii==0:continue
        dct='T'+str(ii)
        if dct not in V[3]:continue
        pp=pmm[yy,:,:,:,:]*nc.variables[dct][0,0,:,:]#).flatten()
        iseq=list(df_twnaq.new_code).index(ii)
        if np.sum(sim[yy,:,:,iseq])<=0:continue #sys.exit('sim<=0 '+dct)
        rat=obs[yu,:,:,iseq]/sim[yy,:,:,iseq]
#        if np.max(rat) >10:sys.exit()
        pr=np.zeros(shape=pp.shape)
        for m in range(12):
          for h in range(124):
#        pr=pp[idx]*rat[idx]
            pr[m,h,:,:]=pp[m,h,:,:]*rat[m,h]
        idx=np.where(pr>0)
        for m in range(12):
          for h in range(124):
            s_idx=[(idx[0][j],idx[1][j],idx[2][j],idx[3][j]) for j in range(len(idx[0])) if (idx[0][j],idx[1][j])==(m,h)]
            if len(s_idx)==0:continue
            y_code.append(y)
            s_code.append(ii)
            m_code.append(m+1)
            d_code.append(h//4+1)
            h_code.append(h%4*6)
            v_code.append(np.mean([pr[s_idx[k]] for k in range(len(s_idx))]))
        pm=pm+list(pr[idx].flatten())
    data.append(pm)
df=DataFrame({'y':y_code,'s':s_code,'m':m_code,'d':d_code,'h':h_code,'v':v_code})
for s in 'ymdh':
  exec('df.'+s+'=[int(i) for i in df.'+s+']')
ymd=[-1 for i in range(len(df))]
for i in range(len(df)):
  y,m,d,h=df.y[i],df.m[i],df.d[i],df.h[i]
  try:
    ymd[i]=int((datetime.datetime(y,m,d,h)+datetime.timedelta(days=8/24.)).strftime('%Y%m%d'))
  except:
    continue
df['ymd']=ymd
df=df.drop(df.loc[df.ymd<0].index).reset_index(drop=True)
dfv=pivot_table(df,index=['ymd','s'],values='v', aggfunc=np.mean).reset_index()
name={i:j for i,j in zip(df_twnaq.new_code,df_twnaq.Name)}
dfv['name']=[name[i] for i in dfv.s]
dfv.set_index('ymd').to_csv('ymd_s_v'+str(iyr)[2:4]+'.csv')
ax.boxplot(data, showfliers=False)
xticks=list(range(0,len(data)+1,intv))
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels, fontsize=8)
plt.savefig('box_AQC.png')


This program was used to calibrate the CAM-chem PM2.5 results with TWEAP station data usage: python YYYY # input the YYYY #model pm results, bundled by grd04_yd.py, or the following "except" block # observation PM2.5 (21 year,13 month, 24 hours, 608 stations # the bin file was bundled by dfpm.py # the township id and fractions(CMAQ landuse formatted, generated by mk_townNew.py) # the township name, ID and stations # observation and simulation values for each township # dict between town code and AQ station number # dict between town code and AQ station X and Y index #year: #UTC #UTC #LST #year[:-1]: #  if y%2==0:continue #UTC #UTC #['Southern']: #).flatten() #sys.exit('sim<=0 '+dct) #        if np.max(rat) >10:sys.exit() #        pr=pp[idx]*rat[idx]