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

Dylan Beaudette debeaudette at ucdavis.edu
Thu Apr 23 22:33:49 CEST 2009


On Thursday 23 April 2009, Greg King wrote:
> 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


Have a look at GRASS, it runs on UNIX.

Cheers,
Dylan

>
> 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]]
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at stat.math.ethz.ch
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341



More information about the R-sig-Geo mailing list