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

Lin Pei-Ling barthealin at hotmail.com
Tue Apr 12 23:28:52 CEST 2011


 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>


More information about the R-help mailing list