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

fnameO = 'PMf13_12_124_137_83.bin'
year=[i for i in range(2007,2020)]
try:
  with FortranFile(fnameO, 'r') as f:
    pmf = f.read_record(dtype=np.float64)  
except:
  pmf=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
        pmf[y-min(year),m,:nt,:,:]=grd04.variables['PM25'][:,0,:,:]
      except:
        continue
  with FortranFile(fnameO, 'w') as f:
    f.write_record(pmf)
fname='/nas1/cmaqruns/2016base/data/land/epic_festc1.4_20180516/gridmask/TWN_CNTY_3X3.nc'
nc = netCDF4.Dataset(fname,'r')
V=[list(filter(lambda x:nc.variables[x].ndim==j, [i for i in nc.variables])) for j in [1,2,3,4]]

df_cnty=read_csv('/nas1/TEDS/teds10_camx/HourlyWeighted/area/cnty2.csv')
df_twnaq=read_csv('/nas1/TEDS/teds10_camx/HourlyWeighted/area/town_aqstE.csv')
df_cnty.no=[int(i) for i in df_cnty.no]
intv=1
xlabels=[' ']
for y in year:
    xlabels.append('{:02d}'.format(y%100))
pmf=pmf.reshape(max(year)-min(year)+1,12, 124, 137, 83)
for i in range(len(df_cnty)):
    fig, ax = plt.subplots()
    n=df_cnty.loc[i,'cnty']
    plt.title(" PM2.5 of CAM-chem Simulations in "+n, loc='center' )
    plt.xlabel('Year(2digits)')
    plt.ylabel('PM2.5 (microgram/cubic meter)')
#   if df_cnty.loc[i,'no']<44:continue
    dct='CNTY_'+'{:02d}'.format(df_cnty.loc[i,'no'])
    if dct not in V[3]:continue
    data=[]
    for y in year:
        pm=[]
        for m in range(12):
            pp=(pmf[y-2007,m,:,:,:]*nc.variables[dct][0,0,:,:]).flatten()
            idx=np.where(pp>0)
            pm=pm+list(pp[idx])
        data.append(pm)
    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_'+n+'.png')