[R] How to set the dimension of a matrix correctly?

[Ben Bolker]

Lin Pei-Ling <barthealin <at> hotmail.com> writes:

>  Hi all, I use kriging to interpolate the precipitation from
> stations, but the map of this results show lots of stripes. (please
> see the attachment)I think there's something wrong with the setting
> of the dimension of this matrix, however, I have no idea how to know
> or test to see if this setting is correct or not.I've tried to
> switch the latitude and longitude, but still got the same results
> (stripes). Hope anyone can take a look at it and give me some
> suggestion.  Thanks for help,Peiling


  My main suggestion is that you try to boil your example down
to something smaller, simpler, and reproducible.  The attention
span of R-helpers is only about 10 minutes, maximum (maybe a little
longer for problems they are particularly interested in), and
it took me almost that long just to reformat your R code so
it was readable.

   Often the process of trying to simplify the problem leads
you to discover your own problem.

  good luck,
   Ben Bolker


  
> 
> ------------------
library(geoR) #functions for geostatistical data analysis
coords<- as.matrix(read.table('/Users/R/Code/stncoords.dat'))
ppt<- as.matrix(read.table('/Users/R/Code/ppt_15day.dat'))
xx <- dim(ppt) # (77,528)
plat <- seq(37.5,42,by=0.07273)
plon <- seq(-105.5,-93.5,by=0.07273)
pgrid <- expand.grid(x=plon,y=plat)
pdim <- dim(pgrid) # (10230,2)
#plot(pgrid, cex=0.5)
lat <- coords[,1] 
lon <- coords[,2] 
ppt1 <- ppt[,1:xx[2]] # 1:528
pptpred <- matrix(0,ncol=xx[2],nrow=1)
 
############# Only test one period data##################
ptemp <- ppt1[,3]		
ll <- which(ptemp>0)
ppt2 <- matrix(0,nrow=length(ll),ncol=3) # (lon,lat,ptemp)		
ppt2[,1] <- lat[ll] # y-axis
ppt2[,2] <- lon[ll] # x-axis
ppt2[,3] <- ptemp[ll] # ppt
pptd <- as.geodata(ppt2)
bin1 <- variog(pptd) 
## 	plot(bin1) # fig1
bin2 <- variog(pptd,estimator.type="modulus")
##	plot(bin2) # fig2
ini1 <- max(bin1$v) 		
ols <- variofit(bin1, fix.nugget = FALSE,
   weights="cressie",ini.cov.pars=c(ini1,4))
kc <- krige.conv(pptd,loc=pgrid,
  krige=krige.control(type.krige="OK",trend.d="2nd",
   trend.l="2nd",cov.pars=ols[2]$cov.pars))
pvalxx <- which(kc$predict < 0)		
kc$predict[pvalxx] <- 0		
## ?? something got mangled here ??
## pptpred <- kc$predict		}else{	  
##      pptpred <- 0*(1:pdim[1])	}}
 
## need to fit the lat/lon frame
newpptpred <- matrix(pptpred, nrow=62, ncol=165) 
plat <- seq(37.5,42,by=0.07273)  ## repeat from above?
plon <- seq(-105.5,-93.5,by=0.07273)  ## repeat from above?
lat2 <- which((plat> 37)&(plat< 42))
lon2 <- which((plon> -106)&(plon< -93))
## fixed ? typos ?
newpptpred <- pptpred[lat2,lon2]
nLevel <- 60
quartz()
filled.contour(plon,plat,t(newpredppt),
   col=rainbow(nLevel),plot.axes={axis(1);axis(2)})