In [1]:
# Load Biospytial modules and etc.
%matplotlib inline
import sys
sys.path.append('/apps')
import django
django.setup()
import pandas as pd
import numpy as np
In [2]:
plotdata_path = "/RawDataCSV/idiv_share/plotsClimateData_11092017.csv"
empirical_data_path = "/apps/external_plugins/spystats/HEC_runs/results/logbiomas_logsppn_res.csv"
In [4]:
from external_plugins.spystats import tools
%run ../HEC_runs/fit_fia_logbiomass_logspp_GLS.py "/RawDataCSV/idiv_share/plotsClimateData_11092017.csv" "/apps/external_plugins/spystats/HEC_runs/results/logbiomas_logsppn_res.csv" -85 -80 30 35
In [5]:
from external_plugins.spystats import tools
hx = np.linspace(0,800000,100)
In [6]:
new_data = results['dataframe']
gvg = results['variogram']
model1 = results['modelGLS']
The object new_data has been reprojected to Alberts and a linear model have been fitted with residuals stored as residuals
In [7]:
new_data.residuals[:10]
Out[7]:
In [8]:
gvg.plot(refresh=False,legend=False,percentage_trunked=20)
plt.title("Semivariogram of residuals $log(Biomass) ~ log(SppR)$")
## HERE we can cast a model (Whittle) and fit it inside the global variogram
whittle_model = tools.WhittleVariogram(sill=0.345,range_a=100000,nugget=0.33,alpha=1.0)
#tt = gvg.fitVariogramModel(whittle_model)
#plt.plot(hx,gvg.model.f(hx),'--',lw=4,c='black')
#print(whittle_model)
In [9]:
## This section is an example for calculating GLS. Using a small section because of computing intensity
In [10]:
whole_var = tools.Variogram(new_data,'logBiomass',model=whittle_model)
In [12]:
MMdist = whole_var.distance_coordinates.flatten()
In [13]:
distances = pd.DataFrame(MMdist)
In [38]:
xx = distances.loc[ (distances[0] < 1000) & (distances[0] > 0)]
In [51]:
xx.sort([0])
Out[51]:
In [9]:
minx = -85
maxx = -80
miny = 30
maxy = 35
section = tools._subselectDataFrameByCoordinates(new_data,'LON','LAT',minx,maxx,miny,maxy)
secvg = tools.Variogram(section,'logBiomass',model=whittle_model)
MMdist = secvg.distance_coordinates.flatten()
CovMat = secvg.model.corr_f(MMdist).reshape(len(section),len(section))
plt.imshow(CovMat)
Out[9]:
In [ ]:
k = resum.as_text
In [52]:
matm = tools.MaternVariogram(sill=0.34,range_a=100000,nugget=0.33,kappa=0.5)
expmm = tools.ExponentialVariogram(sill=0.34,range_a=100000,nugget=0.33)
gausms = tools.GaussianVariogram(sill=0.34,range_a=100000,nugget=0.33)
sphmm = tools.SphericalVariogram(sill=0.34,range_a=100000,nugget=0.33)
wm = tools.WhittleVariogram(sill=0.34,range_a=100000,nugget=0.33,alpha=1)
map(lambda l : l.fit(gvg), [matm,expmm,gausms,sphmm,wm])
print(matm)
print(expmm)
print(gausms)
print(sphmm)
print(wm)
In [53]:
gvg.plot(refresh=False,legend=False,percentage_trunked=20)
plt.title("Semivariogram of residuals $log(Biomass) ~ log(SppR)$")
## HERE we can cast a model (Whittle) and fit it inside the global variogram
#whittle_model = tools.WhittleVariogram(sill=0.345,range_a=100000,nugget=0.33,alpha=1.0)
#tt = gvg.fitVariogramModel(whittle_model)
#plt.plot(hx,matm.f(hx),'--',lw=4)
plt.plot(hx,expmm.f(hx),'--',lw=4)
plt.plot(hx,gausms.f(hx),'--',lw=4)
plt.plot(hx,sphmm.f(hx),'--',lw=4)
plt.plot(hx,wm.f(hx),'--',lw=4)
#print(whittle_model)
Out[53]:
I think the Matern is the right way to go but I need to recheck the model to see if it´s correctly written.
In [54]:
newhx = np.linspace(0,2000,100)
plt.plot(newhx,matm.corr_f(newhx),'--',lw=4)
#plt.plot(hx,expmm.corr_f(hx),'--',lw=4)
#plt.plot(hx,gausms.corr_f(hx),'--',lw=4)
#plt.plot(hx,sphmm.corr_f(hx),'--',lw=4)
plt.plot(newhx,wm.corr_f(newhx),'--',lw=4)
plt.title("Correlation for Matern and Wittle models")
Out[54]:
In [63]:
wm.corr_f(np.linspace(0,100000,100))
Out[63]:
In [57]:
xx
Out[57]:
In [ ]: