In [1]:
import rpy2
print(rpy2.__version__)
In [2]:
from rpy2.robjects.packages import importr
In [3]:
sp = importr('sp')
gstat = importr('gstat')
intamap = importr('intamap')
In [4]:
from rpy2 import robjects
In [5]:
robjects.r('jet.colors <- c("#00007F","blue","#007FFF","cyan","#7FFF7F","yellow","#FF7F00","red","#7F0000")')
robjects.r('col.palette <- colorRampPalette(jet.colors)')
Out[5]:
In [6]:
xy = robjects.r('xy <- read.table("radar_xy.csv",sep=",")')
In [7]:
nrow = robjects.r('nrow <- dim(xy)[1]')
In [8]:
RAD24 = robjects.r('RAD24 <- read.table("./radar_sent/radar_snap_24h_2011_08_05-00_00.csv",sep=",",colClasses="numeric")')
In [9]:
RAD24 = robjects.r('RAD24 <- as.numeric(as.vector(RAD24))')
In [10]:
RAD24 = robjects.r('RAD24 <- RAD24[6:length(RAD24)]')
In [11]:
RAD24
Out[11]:
In [12]:
data = robjects.r('data <- data.frame(xy,RAD24)')
In [13]:
data
Out[13]:
In [14]:
data.colnames = robjects.r('names(data) <- c("x","y","R")')
In [15]:
data
Out[15]:
In [16]:
robjects.r('coordinates(data) <- ~x+y')
Out[16]:
In [17]:
robjects.r('''
png("map.png",height=900,width=900)
ncuts <- 20
cuts <- seq(min(RAD24,na.rm=TRUE),max(RAD24,na.rm=TRUE),length=ncuts)
print(spplot(data["R"],xlab="East [m]",ylab="North [m]",key.space="right",cuts=cuts,region=TRUE,col.regions=col.palette(ncuts),main="Rainfall [mm]",scales=list(draw=TRUE)))
dev.off()
''')
Out[17]:
In [18]:
wet = robjects.r('wet <- which(RAD24>0)')
In [19]:
robjects.r('isotropic_variogram <- variogram(R~1,data[wet,],width=2,cutoff=100) # isotropic variogram with a class width of 2 km and a max. distance of 100 km (wet pixels only)')
Out[19]:
In [20]:
robjects.r('''
dist <- isotropic_variogram$dist # range bins (inter-distances)
gam <- isotropic_variogram$gamma # semivariance estimates for each range bin
np <- isotropic_variogram$np # number of pairs for each range bin
''')
Out[20]:
In [21]:
robjects.r('''
variogram_map <- variogram(R~1,data[wet,],width=2,cutoff=50,map=TRUE)
png("variogram_map.png",height=600,width=600)
print(plot(variogram_map))
dev.off()
''')
Out[21]:
In [ ]:
In [ ]:
In [ ]:
In [37]:
import pandas as pd
import numpy as np
In [40]:
coords = pd.read_csv('./radar_xy.csv', header=None)
coords.columns = ['x', 'y']
coords.head()
Out[40]:
In [43]:
rainfall = pd.read_csv('./radar_sent/radar_snap_24h_2011_08_05-00_00.csv', header=None)
rainfall = pd.DataFrame(rainfall.iloc[0,5::])
rainfall.index = np.arange(0,len(rainfall),1)
rainfall.columns = ['R']
rainfall['x'] = coords['x']
rainfall['y'] = coords['y']
rainfall.head()
Out[43]:
In [44]:
from rpy2.robjects import pandas2ri
In [58]:
mask = rainfall.R>0
In [61]:
rainfall = rainfall[mask]
In [49]:
pandas2ri.activate()
In [62]:
r_df = pandas2ri.py2ri(rainfall)
In [63]:
robjects.r.assign('mydata', r_df)
Out[63]:
In [64]:
robjects.r('head(mydata)')
Out[64]:
In [70]:
robjects.r('''
mydata <- data.frame(mydata)
coordinates(mydata) <- ~x+y
''')
Out[70]:
In [71]:
robjects.r('myiso <- variogram(R~1,mydata,width=2,cutoff=100)')
Out[71]:
In [72]:
robjects.r('''
png("myisotropic_variogram.png",height=600,width=900)
print(plot(myiso))
dev.off()
''')
Out[72]:
In [75]:
asd = robjects.r['myiso']
In [79]:
type(asd)
Out[79]:
In [80]:
print(asd)
In [82]:
asdf = robjects.r('myisomap <- variogram(R~1,mydata,width=2,cutoff=50,map=TRUE)')
In [83]:
import matplotlib.pyplot as plt
%matplotlib inline
In [90]:
robjects.r('''
png("myvariogram_map.png",height=600,width=600)
print(plot(myisomap))
dev.off()
''')
Out[90]:
In [94]:
rainfall.head()
Out[94]:
In [131]:
asdff = pd.DataFrame(rainfall.sort_values(ascending=False, by='R'))
asdff = asdff[0:1499]
robjects.r.assign('sorted_rainfall', asdff)
Out[131]:
In [114]:
robjects.r('head(sorted_rainfall)')
Out[114]:
In [122]:
robjects.r('''
S <- sort(RAD24,index.return=TRUE,decreasing=TRUE)
''')
Out[122]:
In [142]:
robjects.r('I')
Out[142]:
In [124]:
robjects.r('head(RAD24)')
Out[124]:
In [132]:
robjects.r.assign('myrad24', np.array(asdff.R))
Out[132]:
In [138]:
robjects.r('''
S <- sort(RAD24,index.return=TRUE,decreasing=TRUE)
I <- S$ix[1:1499]
hat.anis <- estimateAnisotropy(data[I,],"R") # anisotropy estimates (direction and ratio)
anis <- c(90-hat.anis$direction,1/hat.anis$ratio)
''')
Out[138]:
In [139]:
robjects.r('''
directional_variograms <- variogram(R~1,data[wet,],width=2,cutoff=100,alpha=c(99.9,189.9),tol.hor=5)
png("directional_variograms.png",height=600,width=600)
plot(directional_variograms)
dev.off()
''')
Out[139]:
In [140]:
robjects.r('''
initial_vario_sph <- vgm(psill=500,model="Sph",range=40,nugget=0)
fitted_vario_sph <- fit.variogram(isotropic_variogram,initial_vario_sph)
range <- fitted_vario_sph$range[2] # fitted range
nugget <- fitted_vario_sph$psill[1] # fitted nugget
sill <- sum(fitted_vario_sph$psill) # fitted sill (total sill)
SSErr_sph <- attributes(fitted_vario_sph)$SSErr # sum of squared errors (goodness of fit)
png("fitted_isotropic_variogram_sph.png",height=600,width=900)
print(plot(isotropic_variogram,fitted_vario_sph))
dev.off()
''')
Out[140]:
In [141]:
robjects.r('''
initial_vario_exp <- vgm(psill=500,model="Exp",range=40/3,nugget=0) # pseudo range = range at which you reach 95\% of the sill.
fitted_vario_exp <- fit.variogram(isotropic_variogram,initial_vario_exp)
SSErr_exp <- attributes(fitted_vario_exp)$SSErr
# error is 2.4 times larger than for spherical model
# exponential is clearly not a good choice here.
png("fitted_isotropic_variogram_exp.png",height=600,width=900)
print(plot(isotropic_variogram,fitted_vario_exp))
dev.off()
''')
Out[141]:
In [ ]: