[R-sig-Geo] Spatial Interpolation of Regularly Gridded Data

[Roger Bivand]

On Mon, 6 Apr 2009, Greg King wrote:

> Hi,
>
> I'm fairly new to R, but finding the experience a good one.  However I am a
> little overwhelmed by the number of packages and not sure I am necessarily
> using the most appropriate ones.

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