This is a python / R implementation for spatial analysis of radar rainfall fields. All courtesy for the R code implementation goes to Marc Schleiss

Notes before running:

o) make sure to have installed R properly

o) install the python library rpy2

linux/mac users should be fine by just pip install rpy2
for windows, consider using anaconda and install rpy2

o) install the R libraries sp, gstat and intamap inside the R environment (best as sudo/adminstrator):

install.packages("sp")
install.packages("gstat")
install.packages("intamap")

To plot data in R:

png("myvariogram_map_24h.png",height=600,width=600)
print(plot(myisomap))
dev.off()

Import python / R interface packages



In [ ]:

    
from rpy2.robjects.packages import importr
from rpy2.robjects import r
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

Inside the R environment import these geospatial packages



In [ ]:

    
sp = importr('sp')
gstat = importr('gstat')
intamap = importr('intamap')

Set colors within the R environment



In [ ]:

    
r('jet.colors  <- c("#00007F","blue","#007FFF","cyan","#7FFF7F","yellow","#FF7F00","red","#7F0000")')
r('col.palette <- colorRampPalette(jet.colors)')

Read pandas coordinates of the radar grid



In [ ]:

    
coords = pd.read_csv('./radar_xy.csv', header=None)
coords.columns = ['x', 'y']
coords.head()

Read the 24h dataset and (re)arange the pandas DataFrame



In [ ]:

    
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()

Activate the pandas to R conversion interface



In [ ]:

    
from rpy2.robjects import pandas2ri
pandas2ri.activate()

Select only gridpoints > 0mm rain (wet mask) and assign it in the R environment



In [ ]:

    
mask = rainfall.R>0
rainfall = rainfall[mask]

r_df = pandas2ri.py2ri(rainfall)
r.assign('mydata', r_df)

Transform dataset in R to geospatial dataset



In [ ]:

    
r('''
mydata <- data.frame(mydata)
coordinates(mydata) <- ~x+y
''')

Plot map with python



In [ ]:

    
cur_cmap = plt.cm.jet



In [ ]:

    
plt.scatter(rainfall['x'], rainfall['y'], marker='.', c=rainfall['R'], cmap=cur_cmap)
plt.colorbar()

Bypass plot R map



In [ ]:

    
r('''
RAD24 <- read.table("./radar_sent/radar_snap_24h_2011_08_05-00_00.csv",sep=",",colClasses="numeric")
RAD24 <- as.numeric(as.vector(RAD24))
RAD24 <- RAD24[6:length(RAD24)]
png("map_24h.png",height=900,width=900)
ncuts <- 20
cuts <- seq(min(RAD24,na.rm=TRUE),max(RAD24,na.rm=TRUE),length=ncuts)
print(spplot(mydata["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()
''')

Generate a isotropic variogram (2km separated lags, max 100km)



In [ ]:

    
p_myiso = r('myiso <- variogram(R~1,mydata,width=2,cutoff=100)')



In [ ]:

    
p_myiso.head()



In [ ]:

    
plt.plot(p_myiso['dist'], p_myiso['gamma'], '-o')

Generate and save the 2D variogram map



In [ ]:

    
p_myiso_map = r('myisomap <- variogram(R~1,mydata,width=2,cutoff=50,map=TRUE)')

r('''
png("myvariogram_map_24h.png",height=600,width=600)
print(plot(myisomap))
dev.off()
''')

Investigate the (an)isotropy of the dataset

Only possible with up to 1499 values. Therefore we sort the rainfall values descendingly and assign the sorted dataset to the R environment.



In [ ]:

    
rain_sorted = rainfall.sort_values('R', ascending=False)
rain_sorted = rain_sorted.iloc[0:1499]
rs_df = pandas2ri.py2ri(rain_sorted)
r.assign('data_sorted', rs_df)
r('''
data_sorted <- data.frame(data_sorted)
coordinates(data_sorted) <- ~x+y
''')

Returns the direction of the minumum variablity clockwise from North and the anisotropy ratio



In [ ]:

    
r('''
hat.anis <- estimateAnisotropy(data_sorted,"R")
anis <- c(90-hat.anis$direction,1/hat.anis$ratio)
''')

Compute directional variograms with anisotropy direction



In [ ]:

    
dir_var = r('directional_variograms <- variogram(R~1,mydata,width=2,cutoff=100,alpha=c(99.9,189.9),tol.hor=5)')

dir_1 = dir_var['dir.hor']==99.9
plt.figure()
plt.subplot(121)
plt.plot(dir_var.dist[dir_1], dir_var.gamma[dir_1])
plt.subplot(122)
plt.plot(dir_var.dist[~dir_1], dir_var.gamma[~dir_1])

Fit initial spherical variogram to isotropic variogram



In [ ]:

    
r('initial_vario_sph <- vgm(psill=500,model="Sph",range=40,nugget=0)')
sph_fitted = r('fitted_vario_sph <- fit.variogram(myiso,initial_vario_sph)')
print(sph_fitted)

Save image of fitted variogram



In [ ]:

    
r('''
png("fitted_isotropic_variogram_sph_24h.png",height=600,width=900)
print(plot(myiso,fitted_vario_sph))
dev.off()
''')

fitted range



In [ ]:

    
r('range  <- fitted_vario_sph$range[2]')

fitted nugget



In [ ]:

    
r('nugget <- fitted_vario_sph$psill[1]')

fitted sill



In [ ]:

    
r('sill   <- sum(fitted_vario_sph$psill)')

sum of squared errors



In [ ]:

    
r('SSErr_sph <- attributes(fitted_vario_sph)$SSErr')

Fit exponential model



In [ ]:

    
r('initial_vario_exp <- vgm(psill=500,model="Exp",range=40/3,nugget=0)')
exp_fitted = r('fitted_vario_exp <- fit.variogram(myiso,initial_vario_exp)')

sum of squared errors



In [ ]:

    
r('SSErr_exp  <- attributes(fitted_vario_exp)$SSErr')

Save image



In [ ]:

    
r('''
png("fitted_isotropic_variogram_exp_24h.png",height=600,width=900)
print(plot(myiso,fitted_vario_exp))
dev.off()
''')