[R-sig-Geo] Interpolation / smoothing of points (without raster grids)

[mike.elliott at openreach.co.uk]

My original post on this topic was far more detailed but I got very
little response so my second attempt was much more focussed. Here it is
to see if you are able to offer more help:

Sent: 08 May 2009 12:46
To: r-sig-geo at stat.math.ethz.ch
Subject: [R-sig-Geo] Calculating Risk at point locations
(withoutrasterising)

I would appreciate advice on how to do a spatial statistics task:
detecting points at a high risk of Cases (and cluster analysis). I have
semi-randomly located vector point data with 2 attributes (Cases and
Controls) at each point, i.e. marked point data. I want to apply spatial
smoothing as I hypothesize that the underlying process is such that
clusters of cases are of interest (although I also know that the process
operates on a small enough scale that isolated points with high numbers
of cases are also of interest). Having read the book by Bivand et al
"Applied Spatial Data Analysis with R" I was following the approach of
an inhomogeneous Poisson process analysis (p.175) to calculate a risk
image of log(cases/controls), which incidentally gives a few challenges
with divide by zeros. This risk image (I used data class ppp in R
library SpatStat as I had marked point data) is a raster interpolated
image, even though I originated with vector point data. However, the aim
of my analysis is to score each of the (vector) points with its risk
value, i.e. give each of my original vector points a new attribute
called Risk. For my application I am not interested in what happens in
between my point locations (as the point locations have a meaning). So
my problem is how to get the Risk value at the original point locations.
My current approach involves converting the points to raster, so I lose
the precise point location and I have a large data set so in practice
cannot set a very fine resolution. So my key question is: please, can
someone offer me advice on how to get the Risk value at the original
point locations? 

I wonder if I can avoid raster and do all processing in a vector map
(perhaps by interpolation or spline approximation). This would mean I
can avoid resampling a raster map back to the original vector point
locations as occasionally  the locations are so close (less than
resolution) this will create inaccuracy. 

(Background: I have not done much spatial statistics before so am trying
to keep it simple; I have used GRASS; not much experience with R).

When I originally import the map into R it has class
SpatialPointsDataFrame. Below are a couple of code examples:
> Case <- smooth.ppp(Cases,75, eps=50)   # smooth my Case data with a
kernel width of 75m; 50m pixel dimension in the resulting smoothed ppp.
> plot(eval.im( log((Case+9e-6)/(Control+1e-5)) + 0.3521755)) # to get
the Risk map; the 9e6- and 1e-5 are how I avoid divide by zero type
problems. 

-----Original Message-----
From: Paul Hiemstra [mailto:p.hiemstra at geo.uu.nl] 
Sent: 12 May 2009 10:37
To: Elliott,MR,Mike,BMK6 R
Cc: r-sig-geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Interpolation / smoothing of points (without
raster grids)

mike.elliott at openreach.co.uk wrote:
> Hello - can anyone suggest functions/packages that allow interpolation

> / smoothing of point data (without using raster grids)? i.e. which 
> give the results at the original point locations.
>
> Many thanks, Mike Elliott.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>   
If you could tell us more in detail what you want to do, we could offer
you more guidance. To add to the promotion of packages that we wrote
ourselves (edzer :)), you can also use automap, a package for automatic
interpolation based on gstat and sp (two other packages). For
cross-validation (if this is what you need) the following script will do
the trick:

library(automap)
# Load the data
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y

# Perform cross-validation
kr.cv = autoKrige.cv(log(zinc)~1, meuse, model = c("Exp")) kr_dist.cv =
autoKrige.cv(log(zinc)~sqrt(dist), meuse,
            model = c("Exp"))
kr_dist_ffreq.cv = autoKrige.cv(log(zinc)~sqrt(dist)+ffreq,
           meuse, model = c("Exp"))
# Compare the results
compare.cv(kr.cv, kr_dist.cv, kr_dist_ffreq.cv) compare.cv(kr.cv,
kr_dist.cv, kr_dist_ffreq.cv,
           bubbleplots = TRUE)
compare.cv(kr.cv, kr_dist.cv, kr_dist_ffreq.cv,
           bubbleplots = TRUE, col.names = c("OK","UK1","UK2"))
compare.cv(kr.cv, kr_dist.cv, kr_dist_ffreq.cv,
           bubbleplots = TRUE, col.names = c("OK","UK1","UK2"),
           plot.diff = TRUE)

automap is available from CRAN.

cheers and good luck,
Paul

--
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul