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

Ben Bolker bbolker at gmail.com
Wed Apr 13 14:21:14 CEST 2011

```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

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
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] # 1:528
pptpred <- matrix(0,ncol=xx,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\$cov.pars))
pvalxx <- which(kc\$predict < 0)
kc\$predict[pvalxx] <- 0
## ?? something got mangled here ??
## pptpred <- kc\$predict		}else{
##      pptpred <- 0*(1:pdim)	}}

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

```