Spatial Interpolation of Regularly Gridded Data

Tue, Apr 7, 2009 12:26 AM

On Mon, 6 Apr 2009, Greg King wrote:

Have you read the Spatial Task View on your nearest CRAN? There you should 
find some information to help. Firstly, akima treats all coordinates as 
planar, so don't use interp(). Both fields and gstat can interpolate using 
geographical coordinates, gstat with IDW and kriging, fields with thin 
plate splines (I believe using rdist.earth, but have not checked) and 
kriging. You are trying to interpolate, so why not use a geostatistical 
method? You could see whether the new automap package (running gstat 
internally) can handle geographical coordinates - check by comparing the 
output with runs with the proj4string set and unset.

Hope this helps,

Roger

Here is the background to what I am trying to achieve:  I have a CSV file
which contains weather forecast data for latitude and longitude points (see
attached out.csv, the data I understand is on WGS84
http://www.ready.noaa.gov/faq/geodatums.html).  The sample points are at
half degree intervals.  My objective is to work out what the forecast data
is at any specific given latitude/longitude by interpolating data from the
0.5x0.5 degree grid.  I am doing this for a number of different time points
using the following functions:

library(akima)
library(rgdal)

gks <- function(inLat,inLon,dframe,variab) {
   wind.grid <- expand.grid(lat=inLat, lon=inLon)
   coordinates(wind.grid) <- ~ lat+lon
   proj4string(wind.grid) <- CRS("+init=epsg:4326")
   pnt<-interpp(coordinates(dframe)[,1], coordinates(dframe)[,2],
z=as.data.frame(dframe)[,1],
coordinates(wind.grid)[,1],coordinates(wind.grid)[,2],linear=TRUE)
   return(pnt$z)
}

interp_gk <- function(lat, lon) {
   wind<-read.csv(file="/Users/greg/Out.csv",head=TRUE,sep=",",
colClasses=c("POSIXct","numeric","numeric","numeric","numeric"))
   coordinates(wind)=~lat+lon
   proj4string(wind) <- CRS("+init=epsg:4326")

   times<-unique(wind$vt)
   columns<-names(wind)[2:length(names(wind))]

   dOut<-data.frame(dateTime=times[1])

   for (i in 1:length(times)) {
   dOut[i,"dateTime"]<-times[i]
       for (j in 1:length(columns)) {
           sset<-subset(wind, wind$vt==times[i], select=c(columns[j]))
           dOut[i,columns[j]]<-gks(lat,lon,sset,columns[j])
       }
   }
   dOut<-cbind(dOut, mag=sqrt(dOut$ugrd_0^2+dOut$vgrd_0^2))
   return(dOut)
}

However, I have the following concerns:

  - Should I really be using akima?  The documentation states it is not for
  use on regularly spaced grids - what are my alternatives?
  - The interp funcrtion will not work for cubic spline "linear=FALSE"
  interpolation (is it because my data is regularly gridded?).  How can I
  achieve cubic spline interpolation?
  - Is my function really using the Cordinate reference system specified?
  When I comment out the CRS lines, my functions return the same values?

Lots of questions I appreciate, but I am curious!  It seems R can achieve
what I am trying to do... but I may just be missing some vital bits of
information.

Thanks,

Greg

Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no

Spatial Interpolation of Regularly Gridded Data

Thread (3 messages)