Spatial Interpolation of Regularly Gridded Data
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