# Prepare my slides
%pylab inline
%cd working

"""
Explanation: Python Analysis Evaluation
Author: Barron H. Henderson
End of explanation
"""


!pncaqsraw4pnceval.py --help

"""
Explanation: Process AQS for evaluation

Download annual zip file(s)
Unzip
Use sed, grep, or awk to get spatial/temporal subset
Reshape data
Missing data should be masked
Dimensions (time, point)



CHECK POINT


What do you think the dimensions should be for AQS observations?





What meta-data should be present?





Getting AQS Observations

Get AQS observations
Raw outputs from AQS website
Representational State Transfer - good for small amounts of data


REST format was having problems due to its transition.

Raw outputs

Download directly or download inline
End of explanation
"""


from shapely.wkt import loads
geom = loads("POLYGON ((30 10, 40 35, 20 40, 10 20, 30 10))")
x, y = geom.exterior.xy
plt.plot(x, y, ls = '-', marker = 'o')

"""
Explanation: wktpolygon
This option is most relevant for regional extractions.
End of explanation
"""


!pncaqsraw4pnceval.py -O --timeresolution=daily \
    --start-date 2013-05-01 --end-date 2013-07-01 \
    --wktpolygon "POLYGON ((-181.25 0, 178.75 0, 178.75 90, -181.25 90, -181.25 0))"

%ls -l AQS_DATA_20130501-20130701.nc

"""
Explanation: CHECK  POINT
What should the bounding box be as a WKT Polygon?
ANSWERS Hidden
<div style="visibility: hidden">

"POLYGON ((llcrnrlon llcrnrlat, lrcrnrlon lrcrnrlat, urcrnrlon urcrnrlat, ulcrnrlon ulcrnrlat, llcrnrlon llcrnrlat))"

</div>

Download and Process
End of explanation
"""


!pncdump.py --header AQS_DATA_20130501-20130701.nc

"""
Explanation: Review Output
End of explanation
"""


!pncgen -O -f "bpch,vertgrid='GEOS-5-NATIVE',nogroup=('IJ-AVG-$',)" \
        --extract-file AQS_DATA_20130501-20130701.nc --stack=time -v O3 -s layer72,0 \
        bpch/ctm.bpch.v10-01-public-Run0.2013050100 \
        bpch/ctm.bpch.v10-01-public-Run0.2013050100 \
        bpch_aqs_extract.nc

!pncdump.py --header bpch_aqs_extract.nc

!pnceval.py --help

%%bash
pnceval.py --funcs NO,NP,NOP,MO,MP,MB,RMSE,IOA,AC -v O3 \
--pnc " --expr O3=Ozone*1000;O3.units=\'ppb\' -r time,mean AQS_DATA_20130501-20130701.nc" \
--pnc " -r time,mean bpch_aqs_extract.nc"

from PseudoNetCDF import pnceval
help(pnceval)

"""
Explanation: Extract GEOS-Chem at AQS
End of explanation
"""


from PseudoNetCDF import PNC, pnceval

aqs = PNC("--reduce=time,mean", "--expr=O3=Ozone*1000", "AQS_DATA_20130501-20130701.nc")
geos = PNC("--reduce=time,mean", "bpch_aqs_extract.nc")
aqso3 = aqs.ifiles[0].variables['O3']
geoso3 = geos.ifiles[0].variables['O3']
print(aqso3.shape)
print(geoso3.shape)
print(pnceval.RMSE(aqso3, geoso3))

"""
Explanation: Reproduced in Python

Read in files
Get variables
Calculate MB, RMSE
Modify
Repeat except for month specific results
Repeat except for site specific results



Reproduction Provided
End of explanation
"""


