new_data.crs = {'init':'epsg:4326'}

"""
Explanation: Let´s reproject to Alberts or something with distance
End of explanation
"""


#new_data =  new_data.to_crs("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")

"""
Explanation: Uncomment to reproject
proj string taken from: http://spatialreference.org/
End of explanation
"""


##### OLD #######
len(data.lon)
#X = data[['AET','StandAge','lon','lat']]
#X = data[['SppN','lon','lat']]
X = data[['lon','lat']]
#Y = data['plotBiomass']
Y = data[['SppN']]
## First step in spatial autocorrelation
#Y = pd.DataFrame(np.zeros(len(Y)))
## Let´s take a small sample only for the spatial autocorrelation
#import numpy as np
#sample_size = 2000
#randindx = np.random.randint(0,X.shape[0],sample_size)
#nX = X.loc[randindx]
#nY = Y.loc[randindx]
nX = X
nY = Y



## Small function for systematically selecting the k-th element of the data.
#### Sughgestion use for now a small k i.e. 10
systematic_selection = lambda k : filter(lambda i : not(i % k) ,range(len(data))) 

idx = systematic_selection(50)
print(len(idx))
nX = X.loc[idx]
nY = Y.loc[idx]
new_data = data.loc[idx]

len(new_data)

# Import GPFlow
import GPflow as gf
k = gf.kernels.Matern12(2,ARD=False,active_dims = [0,1]) + gf.kernels.Bias(1)
#k = gf.kernels.Matern12(2, lengthscales=1, active_dims = [0,1] ) + gf.kernels.Bias(1)
l = gf.likelihoods.Poisson()
model = gf.gpmc.GPMC(nX.as_matrix(), nY.as_matrix().reshape(len(nY),1).astype(float), k, l)


#model = gf.gpr.GPR(nX.as_matrix(),nY.as_matrix().reshape(len(nY),1).astype(float),k)
## If priors
#model.kern.matern12.lengthscales.prior = gf.priors.Gaussian(25.0,3.0)
#model.kern.matern32.variance.prior = GPflow.priors.Gamma(1.,1.)
#model.kern.bias.variance.prior = gf.priors.Gamma(1.,1.)
## Optimize

%time model.optimize(maxiter=100) # start near MAP


model.kern

samples = model.sample(50, verbose=True, epsilon=1, Lmax=100)

"""
Explanation: Model Fitting Using a GLM
The general model will have the form:
$$ Biomass(x,y) = \beta_1 AET + \beta_2 Age + Z(x,y) + \epsilon $$
Where:
$\beta_1$ and $\beta_2$ are model parameters, $Z(x,y)$ is the Spatial Autocorrelation process and $\epsilon \sim N(0,\sigma^2)$
End of explanation
"""


model.kern.lengthscales = 25.4846122373
model.kern.variance = 10.9742076021
model.likelihood.variance = 4.33463026664

%time mm = k.compute_K_symm(X.as_matrix())

import numpy as np
Nn = 500
dsc = data
predicted_x = np.linspace(min(dsc.lon),max(dsc.lon),Nn)
predicted_y = np.linspace(min(dsc.lat),max(dsc.lat),Nn)
Xx, Yy = np.meshgrid(predicted_x,predicted_y)
## Fake richness
fake_sp_rich = np.ones(len(Xx.ravel()))
predicted_coordinates = np.vstack([ Xx.ravel(), Yy.ravel()]).transpose()
#predicted_coordinates = np.vstack([section.SppN, section.newLon,section.newLat]).transpose()

len(predicted_coordinates)

#We will calculate everything with the new model and parameters
#model = gf.gpr.GPR(X.as_matrix(),Y.as_matrix().reshape(len(Y),1).astype(float),k)

%time means,variances = model.predict_y(predicted_coordinates)

-

"""
Explanation: Fitted parameters (From HEC)
End of explanation
"""


#Using k-partition = 7
import cartopy
plt.figure(figsize=(17,11))

proj = cartopy.crs.PlateCarree()
ax = plt.subplot(111, projection=proj)


ax = plt.axes(projection=proj)
#algo = new_data.plot(column='SppN',ax=ax,cmap=colormap,edgecolors='')


#ax.set_extent([-93, -70, 30, 50])
ax.set_extent([-125, -60, 20, 50])
#ax.set_extent([-95, -70, 25, 45])

#ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.add_feature(cartopy.feature.LAKES, alpha=0.9)
ax.stock_img()
#ax.add_geometries(new_data.geometry,crs=cartopy.crs.PlateCarree())
#ax.add_feature(cartopy.feature.RIVERS)
mm = ax.pcolormesh(Xx,Yy,means.reshape(Nn,Nn) + (2* np.sqrt(variances).reshape(Nn,Nn)),transform=proj )
#cs = plt.contour(Xx,Yy,np.sqrt(variances).reshape(Nn,Nn),linewidths=2,cmap=plt.cm.Greys_r,linestyles='dotted')
cs = plt.contour(Xx,Yy,means.reshape(Nn,Nn) + (2 * np.sqrt(variances).reshape(Nn,Nn)),linewidths=2,colors='k',linestyles='dotted',levels=range(1,20))
plt.clabel(cs, fontsize=16,inline=True,fmt='%1.1f')
#ax.scatter(new_data.lon,new_data.lat,edgecolors='',color='white',alpha=0.6)
plt.colorbar(mm)
plt.title("Predicted Species Richness + 2stdev")

"""
Explanation: Predictions with +2std.Dev
End of explanation
"""


#Using k-partition = 7
import cartopy
plt.figure(figsize=(17,11))

proj = cartopy.crs.PlateCarree()
ax = plt.subplot(111, projection=proj)


ax = plt.axes(projection=proj)
#algo = new_data.plot(column='SppN',ax=ax,cmap=colormap,edgecolors='')


#ax.set_extent([-93, -70, 30, 50])
ax.set_extent([-125, -60, 20, 50])
#ax.set_extent([-95, -70, 25, 45])

#ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.add_feature(cartopy.feature.LAKES, alpha=0.9)
ax.stock_img()
#ax.add_geometries(new_data.geometry,crs=cartopy.crs.PlateCarree())
#ax.add_feature(cartopy.feature.RIVERS)
mm = ax.pcolormesh(Xx,Yy,means.reshape(Nn,Nn),transform=proj )
#cs = plt.contour(Xx,Yy,np.sqrt(variances).reshape(Nn,Nn),linewidths=2,cmap=plt.cm.Greys_r,linestyles='dotted')
cs = plt.contour(Xx,Yy,means.reshape(Nn,Nn),linewidths=2,colors='k',linestyles='dotted',levels=range(1,20))
plt.clabel(cs, fontsize=16,inline=True,fmt='%1.1f')
#ax.scatter(new_data.lon,new_data.lat,edgecolors='',color='white',alpha=0.6)
plt.colorbar(mm)
plt.title("Predicted Species Richness")


#Using k-partition = 7
import cartopy
plt.figure(figsize=(17,11))

proj = cartopy.crs.PlateCarree()
ax = plt.subplot(111, projection=proj)


ax = plt.axes(projection=proj)
#algo = new_data.plot(column='SppN',ax=ax,cmap=colormap,edgecolors='')


#ax.set_extent([-93, -70, 30, 50])
ax.set_extent([-125, -60, 20, 50])
#ax.set_extent([-95, -70, 25, 45])

#ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.add_feature(cartopy.feature.LAKES, alpha=0.9)
ax.stock_img()
#ax.add_geometries(new_data.geometry,crs=cartopy.crs.PlateCarree())
#ax.add_feature(cartopy.feature.RIVERS)
mm = ax.pcolormesh(Xx,Yy,means.reshape(Nn,Nn) - (2* np.sqrt(variances).reshape(Nn,Nn)),transform=proj )
#cs = plt.contour(Xx,Yy,np.sqrt(variances).reshape(Nn,Nn),linewidths=2,cmap=plt.cm.Greys_r,linestyles='dotted')
cs = plt.contour(Xx,Yy,means.reshape(Nn,Nn) - (2 * np.sqrt(variances).reshape(Nn,Nn)),linewidths=2,colors='k',linestyles='dotted',levels=[4.0,5.0,6.0,7.0,8.0])
plt.clabel(cs, fontsize=16,inline=True,fmt='%1.1f')
#ax.scatter(new_data.lon,new_data.lat,edgecolors='',color='white',alpha=0.6)
plt.colorbar(mm)
plt.title("Predicted Species Richness - 2stdev")

"""
Explanation: Predicted means
End of explanation
"""


model.get_parameter_dict()

"""
Explanation: Model Analysis
End of explanation
"""


X_ = data[['LON','LAT']]
%time Y_hat = model.predict_y(X_)

pred_y = pd.DataFrame(Y_hat[0])
var_y = pd.DataFrame(Y_hat[1])

new_data['pred_y'] = pred_y
new_data['var_y'] = var_y

new_data= new_data.assign(error=lambda y : (y.SppN - y.pred_y)**2 )

new_data.error.hist(bins=50)

print(new_data.error.mean())
print(new_data.error.std())

"""
Explanation: Let's calculate the residuals
End of explanation
"""


import raster_api.tools as rt
from raster_api.models import MeanTemperature,ETOPO1,Precipitation,SolarRadiation
from sketches.models import Country


## Select US
us_border = Country.objects.filter(name__contains='United States')[1]

from django.db import close_old_connections

close_old_connections()

#Get Raster API
us_meantemp = rt.RasterData(Precipitation,us_border.geom)
us_meantemp.getRaster()

us_meantemp.display_field()

%time coords = us_meantemp.getCoordinates()

"""
Explanation: Experiment
In this section we will bring a Raster Data from the US, using Biospytial Raster API. 
1. First select a polygon, then get A Raster FRom there, say Mean Temperature.
End of explanation
"""


