<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 ICARTT Evaluation
Author
Step1:  Evaluation
Step one pairing data
Homogonize dimensions
Use test metrics
Pairing Data - Better Interpolation
In PNC_02Figures , we masked data instead of matching data.
Ideally, we would have a multi-dimensional interpolation.
For this we will use scipy KDTree
First, load the data again
Step2:  Read Model Data
Step3:  Read Observations
Step4:  Scipy has a spatial library
Step5:  Perform 2-d nearest neighbor interpolation
Step6:  Simple Regression
Step7:  1 Include More Variables
Step8:  2 Nearest Neighbor using pressure/lon/lat as peers
Step9:  Plot Paired Results
Step10:  Lon/Lat First
Step11:  Plot the results
Step12:  Remove Assymetric Variability
Benchmark has 1 mean-value all months
DC3 has 1 per minute
Step13:  Plot Paired Time-Means
Step14:  Interpolation not KDTree
KDTree is just finding a cell
Step15:  Questions
Are KDTree and NearestNDInterpolator the same?
How different would an analysis be if we had
Step16:  Non-Parametric Tests
Spearman Correlation (spearmanr)
Mann Whitney U (mannwhitneyu)
Wilcoxon signed-rank test (wilcoxon)
<ASSISTANT_TASK:>
Python Code:
# Prepare my slides
%pylab inline
%cd working
Explanation: ICARTT Evaluation
Author: Barron H. Henderson
This presentation will teach basics of spatial interpolation and statistical evaluation using Python with PseudoNetCDF, numpy, and scipy.
End of explanation
from PseudoNetCDF import PNC
Explanation: Evaluation
Step one pairing data
Homogonize dimensions
Use test metrics
Pairing Data - Better Interpolation
In PNC_02Figures , we masked data instead of matching data.
Ideally, we would have a multi-dimensional interpolation.
For this we will use scipy KDTree
First, load the data again
End of explanation
margs = PNC('--format=bpch,nogroup=("IJ-AVG-$","PEDGE-$"),vertgrid="GEOS-5-NATIVE"', \
            "bpch/ctm.bpch.v10-01-public-Run0.2013050100")
mfile = margs.ifiles[0]
mlon = mfile.variables['longitude']
mlat = mfile.variables['latitude']
MLON, MLAT = np.meshgrid(mlon, mlat)
mO3 = mfile.variables['O3']
Explanation: Read Model Data
End of explanation
oargs = PNC("-f", "ffi1001",
            "--rename=v,O3_ESRL,O3",
            "--variables=Fractional_Day,O3_ESRL,LATITUDE,LONGITUDE,PRESSURE",
            "--rename=v,LATITUDE,latitude",
            "--expr=longitude=np.where(LONGITUDE>180,LONGITUDE-360,LONGITUDE);",
            "--expr=time=Fractional_Day*24*3600;time.units=\"seconds since 2011-12-31\"",
            "icartt/dc3-mrg60-dc8_merge_20120518_R7_thru20120622.ict")
ofile = oargs.ifiles[0]
olons = ofile.variables['longitude']
olats = ofile.variables['latitude']
oO3 = ofile.variables['O3']
Explanation: Read Observations
End of explanation
import scipy.spatial
?scipy.spatial
?scipy.spatial.KDTree
Explanation: Scipy has a spatial library
End of explanation
from scipy.spatial import KDTree
tree = KDTree(np.ma.array([MLAT.ravel(), MLON.ravel()]).T)
dists, idxs = tree.query(np.ma.array([olats, olons]).T)
latidxs, lonidxs = np.unravel_index(idxs, MLAT.shape)
plt.close()
mpO3 = mO3[:, :, latidxs, lonidxs]
mpO3.shape
moO3 = mpO3[:, 22]
moO3.shape
Explanation: Perform 2-d nearest neighbor interpolation
End of explanation
from scipy.stats.mstats import linregress
?linregress
lr = linregress(oO3[:], moO3[:])
lr
x = np.array(oO3[:].min(), oO3[:].max())
plt.scatter(oO3[:], moO3)
plt.plot(x, lr.slope * x + lr.intercept, ls = '--')
plt.text(x.mean(), plt.ylim()[1], 'r=%.2f; p=%.2f' % (lr.rvalue, lr.pvalue))
Explanation: Simple Regression
End of explanation
from scipy.spatial import KDTree
from PseudoNetCDF import PNC
# Read Model Data
margs = PNC('-c', "layer73,valid,0.5,0.5", \
            '--format=bpch,nogroup=("IJ-AVG-$","PEDGE-$"),'+
            'vertgrid="GEOS-5-NATIVE"', \
            "bpch/ctm.bpch.v10-01-public-Run0.2013050100")
mfile = margs.ifiles[0]
mlon = mfile.variables['longitude']
mlat = mfile.variables['latitude']
mpres = mfile.variables['PSURF']
mO3 = mfile.variables['O3']
# Read Observations
oargs = PNC("-f", "ffi1001", "--rename=v,O3_ESRL,O3", \
            "--variables=Fractional_Day,O3_ESRL,LATITUDE,LONGITUDE,PRESSURE", \
            "--expr=latitude=LATITUDE;"+\
            "longitude=np.where(LONGITUDE>180,LONGITUDE-360,LONGITUDE);"+
            "time=Fractional_Day*24*3600;time.units=\"seconds since 2011-12-31\"", 
            "icartt/dc3-mrg60-dc8_merge_20120518_R7_thru20120622.ict")
ofile = oargs.ifiles[0]
olons = ofile.variables['longitude']
olats = ofile.variables['latitude']
opres = ofile.variables['PRESSURE']
oO3 = ofile.variables['O3']
Explanation: 1 Include More Variables
End of explanation
MLON, MLAT = np.meshgrid(mlon, mlat)
MLAT.shape
MLON = MLON * np.ones_like(mO3[0])
MLAT = MLAT * np.ones_like(mO3[0])
MLAT.shape
tree = KDTree(np.ma.array([mpres.ravel(), MLAT.ravel(), MLON.ravel()]).T)
dists, idxs = tree.query(np.ma.array([opres, olats, olons]).T)
presidxs, latidxs, lonidxs = np.unravel_index(idxs, MLAT.shape)
moO3 = mO3[:, presidxs, latidxs, lonidxs][0]
moO3.shape
Explanation: 2 Nearest Neighbor using pressure/lon/lat as peers
End of explanation
from scipy.stats.mstats import linregress
lr = linregress(oO3[:], moO3[:])
lr
plt.scatter(oO3[:], moO3)
plt.plot([0,500], [0,500], ls = '-', lw = 3, color = 'k')
x = np.array([0, 500])
plt.plot(x, lr.slope * x + lr.intercept, ls = '--')
plt.text(200, 500, 'r=%.2f; p=%.2f' % (lr.rvalue, lr.pvalue))
plt.xlim(0, 500)
plt.ylim(0, 500)
plt.ylabel('Model')
plt.xlabel('Obs');
Explanation: Plot Paired Results
End of explanation
MLON, MLAT = np.meshgrid(mlon, mlat)
tree = KDTree(np.ma.array([MLAT.ravel(), MLON.ravel()]).T)
dists, idxs = tree.query(np.ma.array([olats, olons]).T)
latidxs, lonidxs = np.unravel_index(idxs, MLAT.shape)
meO3 = mO3[:, :, latidxs, lonidxs]
meO3.shape
dp = opres[None, None,:] - mpres[:, :, latidxs, lonidxs]
pidx = np.abs(dp).argmin(1)[0]
moO3 = meO3[:, pidx, np.arange(pidx.size)][0]
moO3.shape
Explanation: Lon/Lat First
End of explanation
from scipy.stats.mstats import linregress
lr = linregress(oO3[:], moO3[:])
plt.scatter(oO3[:], moO3)
x = np.array([0, 500])
plt.plot(x, x, ls = '-', lw = 3, color = 'k')
plt.axis('square')
plt.xlim(*x)
plt.ylim(*x)
plt.ylabel('Model')
plt.xlabel('Obs')
plt.plot(x, lr.slope * x + lr.intercept, ls = '--')
plt.text(200, 500, 'r=%.2f; p=%.2f' % (lr.rvalue, lr.pvalue));
Explanation: Plot the results
End of explanation
idxs = [i for i in zip(pidx, latidxs, lonidxs)]
uniqueset = set(idxs)
idxs = np.array(idxs)
ovals = []
mvals = []
for idx in uniqueset:
    oi = (idxs == idx).all(1)
    ovals.append(oO3[oi].mean())
    mvals.append(mO3[0][idx])
print(len(ovals), len(mvals))
Explanation: Remove Assymetric Variability
Benchmark has 1 mean-value all months
DC3 has 1 per minute
End of explanation
from scipy.stats.mstats import linregress
lr = linregress(ovals[:], mvals[:])
plt.scatter(ovals, mvals)
plt.ylabel('Model'); plt.xlabel('Obs')
plt.axis('square')
x = np.array([0, 350])
plt.xlim(*x); plt.ylim(*x)
plt.plot(x, lr.slope * x + lr.intercept, ls = '--')
plt.title('r=%.2f; p=%.2f' % (lr.rvalue, lr.pvalue));
Explanation: Plot Paired Time-Means
End of explanation
import scipy.interpolate
?scipy.interpolate.NearestNDInterpolator
?scipy.interpolate.LinearNDInterpolator
Explanation: Interpolation not KDTree
KDTree is just finding a cell
End of explanation
from scipy.stats.mstats import ttest_ind
ttr = ttest_ind(ovals, mvals)
ttr
Explanation: Questions
Are KDTree and NearestNDInterpolator the same?
How different would an analysis be if we had:
used LinearNDInterpolator
would time normalization be necessary
Testing Populations
Typical with normal assumption
Non-parametric tests
Parametric Tests
Pearson Correlation (pearsonr)
Students t-test (ttest_ind)
Shapiro-Wilk (shapiro)
End of explanation
from scipy.stats.mstats import mannwhitneyu
mwr = mannwhitneyu(ovals, mvals)
mwr
Explanation: Non-Parametric Tests
Spearman Correlation (spearmanr)
Mann Whitney U (mannwhitneyu)
Wilcoxon signed-rank test (wilcoxon)
End of explanation
<END_TASK>