External Drift Kriging
Hi Luke
krige.conv command. I am trying to krige annual average temperatures from 112 met stations located in the Eastern Cascade mountains of Oregon. The data are highly non-stationary, due to the presence of a linear lapse response between the temperature and the altitude (the altitudinal range here is considerable). I have quantified this response at a loss of 3.5 degrees centegrade per 1km rise. I have removed this trend from the data, and produced a zero mean, gaussian distributed stationary process, and from this derived the spatial covarince structure.
OK. This is one option indeed, to modelling covariance after getting "residuals" from some kind of model to take auto the effect of known sources of variation. And indeed the usual variogram will only make sense for stationary data. In the max. lik. estimation you can do both together, getting estimates for a (linear) model for the mean and for the covariance parameters. An example call would be something along the lines
ml <- likfit(foo, trend=~altitude)
provided foo is a geodata object with covariate information on altitude Later for kriging you would need a call of the type:
krige.conv(foo, loc=gr, krige=krige.control(obj=ml, trend.d=~altitude,
trend.l = ~altitude.gr), borders=NULL) where gr is a matrix with coordinates of prediction points and altitude.gr is the value of the covariate at those prediction points,
The problem is that I only want to predict values from a small subset of the area which the stations occupy; although some stations are outwith the required borders, they are more or less within the range of the covariance function. I have an elevation model for the prediction locations, but due to computational limitations, cannot extend this to cover the full area: thus trend.l is only specified within the borders.
So you can only krige within the borders, unless you extrapolate your covariate information somehow. Kriging with external trend (w/ covariates) requires values of the coavriates at all prediction locations. Hope this helps, and if the problempersists feel free to contact me directly best P.J.
All of the stations have elevations attached, and so I am able to
specify trend.d. To summarise I have the following:
grid [a 1km grid which covers the entire domain of
interest, containing all the met stations]
lc [a two dimensional matrix specifying the
borders for prediction]
Normals [a geodata object with data on the annual station
temperature average]
pts [a geodata object with data on the station
elevations]
cressie_cu [a variogram model (cressie weights, cubic model),
derived from the residuals of the linear model average_temp~elevation]
DEM [a geodata object with elevation data at 1km
spacing for all prediction locations within lc]
I am attempting the following:
pars <- krige.control(obj.m=cressie_cu,
trend.d = trend.spatial(~data, pts),
trend.l = trend.spatial(~data,DEM))
kc <- krige.conv(Normals, loc = grid, borders = lc, krige = pars)
On trying to execute the commands, I get an output which has lots of
lines and other artifacts running through it, which seems to bear no
resemblance to the expected output at all. Could you please suggest what
I am getting wrong?
Many thanks
Luke Spadavecchia
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Paulo Justiniano Ribeiro Jr LEG (Laborat?rio de Estat?stica e Geoinforma??o) Departamento de Estat?stica Universidade Federal do Paran? Caixa Postal 19.081 CEP 81.531-990 Curitiba, PR - Brasil Tel: (+55) 41 3361 3573 Fax: (+55) 41 3361 3141 e-mail: paulojus at est.ufpr.br http://www.est.ufpr.br/~paulojus