[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