[R-sig-Geo] Spatial Interpolation of Regularly Gridded Data
Clint Bowman
clint at ecy.wa.gov
Thu Apr 23 22:54:06 CEST 2009
Greg,
I'm quite certain that kriging would not be appropriate with model
output--you won't have any errors to produce the necessary statistics.
Are you looking at smooth terrain (e.g., over the ocean) where topographic
influences are negligible and you may want to interpolate the pressure
field and then diagnose the resulting winds, or is your interest high
enough in the atmosphere where you could use the thermal wind equation?
Recognize that the shortest resolved wavelength will likely be
approximately two grid lengths which will limit the ability of any
interpolation to match a sub-grid scale effect. You can probably get a
better idea of the correct approach by looking in detail at the model
description--is it a spectral model? The interpolation used by the model
(you may need to look at the model source code or communicate with the
model develoopers/maintainers) to determine an appropriate interpolation
method.
Clint
Clint Bowman INTERNET: clint at ecy.wa.gov
Air Dispersion Modeler INTERNET: clint at math.utah.edu
Air Quality Program VOICE: (360) 407-6815
Department of Ecology FAX: (360) 407-7534
USPS: PO Box 47600, Olympia, WA 98504-7600
Parcels: 300 Desmond Drive, Lacey, WA 98503-1274
On Thu, 23 Apr 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
>
> 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
>
More information about the R-sig-Geo
mailing list