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

From: Lin Pei-Ling <barthealin_at_hotmail.com>
Date: Tue, 12 Apr 2011 16:28:52 -0500

 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)})                                                



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Received on Tue 12 Apr 2011 - 22:07:18 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 12 Apr 2011 - 22:40:28 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive