Skip to content

How to set the dimension of a matrix correctly?

2 messages · Lin Pei-Ling, Ben Bolker

Lin Pei-Ling
# 
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

------------------
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 = F,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		pptpred <- kc$predict		}else{	        pptpred <- 0*(1:pdim[1])	}}

newpptpred <- matrix(pptpred, nrow=62, ncol=165) # need to fit the lat/lon frame
plat <- seq(37.5,42,by=0.07273)plon <- seq(-105.5,-93.5,by=0.07273)
lat2 <- which((plat> 37)&(plat< 42))lon2 <- which((plon> -106)&(plon< -93))
newpptpred <- tpptpred[lat2,pon2]
nLevel <- 60
quartz()filled.contour(plon,plat,t(newpredppt),col=rainbow(nLevel),plot.axes={axis(1);axis(2)})


 		 	   		  
-------------- next part --------------
A non-text attachment was scrubbed...
Name: predtest.pdf
Type: application/pdf
Size: 969114 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20110412/a6d1b775/attachment-0001.pdf>
Ben Bolker
# 
Lin Pei-Ling <barthealin <at> hotmail.com> writes:
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)})