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