Spatial Interpolation of Regularly Gridded Data

This reply is a little late, but after much searching, I cannot find an R
package/function that fits my needs.

I am trying to perform either bilinear or  bicubic interpolation using
great circle distances (ie using a spherical model of the earth).  Here is
what I have tried:

gstat: Can only perform kriging or IDW
akima: Cannot interpolate regularly gridded data and does not support
Coordinate Reference Systems
fields: Thin plate spline function does not seem to support regularly
gridded data or Coordinate Refernec Systems
rsaga: Only supported on windows (I need unix based)

I'm now at a bit of a loss where to look.  It seems the regular R functions
are very good at bilinear and bicubic interpolation in cartesian geometry,
but I can't find a way of using great circle distances for interpolation.
Are there any more avenues I should investigate?

Thanks,

Greg

On Tue, Apr 7, 2009 at 11:06 PM, Scott W Mitchell <

smitch at connect.carleton.ca> wrote:

Greg, you just need to come back here!  We just hired Murray Richardson,
and he's a whiz at R-SAGA.

;)

Cheers,
Scott



On 7-Apr-09, at 16:32 , Greg King wrote:

 Thanks for your reply Roger.  I little while back I got kriging working
in

the gstat package with my data, but found that for my purposes it was
too computationally expensive (read slow).  Therefore I am looking for
something
a little quicker.  I now appreciate I cannot use interp due to it
treating coordinates as planar.  I found a good tutorial for kriging in
gstat at http://casoilresource.lawr.ucdavis.edu/drupal/node/442.  From a
quick play
around I can't immediately see how functions like Tps in the field
package meet my needs.  Could anyone provide directions to primers on
the web?. Googling for R is quite tricky (even with rseek).

Essentially I want a fast simple method for interpolating a single point
value from its surounding points.  It shouldn't be too difficult but I
must
admit I am struggling a little!

Greg

On Tue, Apr 7, 2009 at 8:26 AM, Roger Bivand <Roger.Bivand at nhh.no>
wrote:

 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

       [[alternative HTML version deleted]]

	[[alternative HTML version deleted]]

This reply is a little late, but after much searching, I cannot find an R
package/function that fits my needs.

I am trying to perform either bilinear or  bicubic interpolation using great
circle distances (ie using a spherical model of the earth).  Here is what I
have tried:

gstat: Can only perform kriging or IDW
akima: Cannot interpolate regularly gridded data and does not support
Coordinate Reference Systems
fields: Thin plate spline function does not seem to support regularly
gridded data or Coordinate Refernec Systems
rsaga: Only supported on windows (I need unix based)

I'm now at a bit of a loss where to look.  It seems the regular R functions
are very good at bilinear and bicubic interpolation in cartesian geometry,
but I can't find a way of using great circle distances for interpolation.
Are there any more avenues I should investigate?

Thanks,

Greg

On Tue, Apr 7, 2009 at 11:06 PM, Scott W Mitchell <
smitch at connect.carleton.ca> wrote:

Greg, you just need to come back here!  We just hired Murray Richardson,
and he's a whiz at R-SAGA.

;)

Cheers,
Scott



On 7-Apr-09, at 16:32 , Greg King wrote:

 Thanks for your reply Roger.  I little while back I got kriging working in

the gstat package with my data, but found that for my purposes it was too
computationally expensive (read slow).  Therefore I am looking for
something
a little quicker.  I now appreciate I cannot use interp due to it treating
coordinates as planar.  I found a good tutorial for kriging in gstat at
http://casoilresource.lawr.ucdavis.edu/drupal/node/442.  From a quick
play
around I can't immediately see how functions like Tps in the field package
meet my needs.  Could anyone provide directions to primers on the web?.
Googling for R is quite tricky (even with rseek).

Essentially I want a fast simple method for interpolating a single point
value from its surounding points.  It shouldn't be too difficult but I
must
admit I am struggling a little!

Greg

On Tue, Apr 7, 2009 at 8:26 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:

 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

       [[alternative HTML version deleted]]

	[[alternative HTML version deleted]]