[R-sig-Geo] External Drift Kriging

[Luke Spadavecchia]

Hello!

I am currently studying my PhD at Edinburgh University, and I have 
recently been moving to R as my main data analysis package. One of the 
main reasons for this is the availability of GeoR, which will allow me 
to do everything I need to do in one environment. However, I am 
currently having a few problems getting up and running with the 
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.

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