[R-sig-Geo] Plot variograms ensembles in a single plot
subash
subbu_843 at yahoo.co.in
Thu Dec 20 08:31:01 CET 2012
I am beginner in R.I would like to see all the varigrams ensembles of 8
months rainfall data in a single plot. I am using autofitVariogram and
autoKrige functions.
The following code adapted using the Paul Hiemstra trick at
https://stat.ethz.ch/pipermail/r-sig-geo/2009-December/007109.html
<https://stat.ethz.ch/pipermail/r-sig-geo/2009-December/007109.html> does
the job, but it plots them in * (panels)/subplots*. I want all the
variograms in a single plot.
I tried to see the parameters in *xyplot * but could'nt find the solution.
Thanks in advance
# Import the required libraries
library(rgdal)
library(maptools)
library(gstat)
library(sp)
library(automap)
library(XLConnect)
remove(list=ls())
workdir = "/home/subash/MyData/KSNDMC/Homogenous Regions/2011/JJAS/"
fileName = "grd1_dly_Rainfall_JJAS2011.xls"
path_fileName = paste(workdir,fileName,sep='') #Concatenation
rawdata = readWorksheetFromFile(path_fileName,"Sheet3", region = "D1:DC123",
header = FALSE)
N = nrow(rawdata)
rain = rawdata[4:N,]
lat = as.numeric(t(rawdata[2,]))
lon = as.numeric(t(rawdata[3,]))
lonlat = cbind(lon,lat)
sp = SpatialPoints(lonlat,proj4string = CRS("+proj=longlat"))
sp_utm = spTransform(sp, CRS("+proj=utm +zone=43N +datum=WGS84"))
av_list = list()
nRows = nrow(rain)
for (i in 1:nRows)
{
irain = rain[i,]
indx = (irain == "NaN")
irain = irain[!indx]
irain = as.numeric(irain)
isallzeros = max(irain) # To take care of the cases of dry day(irain =0)
at all stations
irain = as.data.frame(irain)
M = nrow(irain)
if (M > 5) # To avoid cases of NaN across many stations
{
if (!(isallzeros == 0))
{
print(i)
foo_utm = sp_utm[!indx]# Removing the locations with NaN values
data = data.frame(foo_utm,irain)
names(data) = c("Easting","Northing","rain")
coordinates(data) = c("Easting","Northing")
variogram = autofitVariogram(rain~1,data)
# plot(variogram)
av_list[[i]] = variogram
}
}
}
sv_list = do.call("rbind", lapply(1:nRows, function(x) {
sv = av_list[[x]]$exp_var
sv$identification = as.character(x)
return(sv)
}))
# Make a list of fitted models
model_list = lapply(av_list, function(x) x$var_model)
xyplot(gamma ~ dist|identification,sv_list,
panel = function(...) {
ret = variogramLine(model_list[[panel.number()]], maxdist =
max(sv_list$dist))
llines(ret$dist, ret$gamma, col = "red")
})
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Plot-variograms-ensembles-in-a-single-plot-tp7582030.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
More information about the R-sig-Geo
mailing list