Skip to content

functions on kernel density estimation in point pattern analysis

8 messages · zhijie zhang, Tim Keitt, Rob Robinson +2 more

Original
Roger Bivand
#
On Sun, 24 Jun 2007, zhijie zhang wrote:

            

      
      
      
    

        
You might consider possibilities in spatstat, and it would be easier to 
answer if your needs were clearer - it is possible that the splancs 
functions could work better with different argument values, but without 
knowing what you used, it is difficult to tell.

Roger

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Tim Keitt
    

    

      
      #
    

  


  

      
I rather like locfit.

THK

        
On 6/24/07, zhijie zhang <epistat at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Dan Putler
    

    

      
      #
    

  


  

      
Hi All,

To add some detail to Roger's earlier, post density.ppp in the  
spatstat seems to be a very good answer to the original post since it  
is specifically designed to estimate a kernel density for a point  
process pattern. This function use a bivariate Gaussian smoother that  
lends itself to user configuration.

Dan

        
On 24-Jun-07, at 5:28 PM, Tim Keitt wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Rob Robinson
    

    

      
      #
    

  


  

      
Have you had a look at package ks? It seems to do rather nice bivariate
kernels, and (I think) generalises to other dimensions.
Cheers
Rob

*** Want to know about Britain's birds? Try  www.bto.org/birdfacts ***

Dr Rob Robinson, Senior Population Biologist
British Trust for Ornithology, The Nunnery, Thetford, Norfolk, IP24 2PU
Ph: +44 (0)1842 750050         E: rob.robinson at bto.org
Fx: +44 (0)1842 750030         W: http://www.bto.org
eSafe scanned this email for viruses, vandals and malicious content (!)

==== "How can anyone be enlightened, when truth is so poorly lit" =====

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      zhijie zhang
    

    

      
      #
    

  


  

      
Dear Roger Bivand,Tim Keitt, and Dan Putler,
  Thanks for your answers. I have tried density.ppp(spatstat) and
kernel2d(splancs), but the results are not very satisfied. I think there
should be a higher density in the blue part of the map in the attachment.
  My dataset has been put in the attachment and programs have been pasted in
the following, so that u can use and check it.
   I want to do the work like "non-parametric estimation of a spatially
varying intensity" in Diggle's book(2003.P.116-121).
  BTW, i'm not familiar with locfit,would u please also check it using
locfit? Thanks very much.
###############################################################################
                           Kernel density estimation--spatstat
################################################################################
library(sp)
library(foreign)
library(maptools)
library(mgcv)
library(spatstat)

guichi<-readShapePoly("d:/deleting/kernel/kernel/guichi2.shp")
W <-as(as(guichi, "SpatialPolygons"), "owin")  #boundary polygons

cases<-coordinates(readShapePoints("d:/deleting/kernel/kernel/cases.shp"))
#points
colnames(cases)<-c("x","y")
cases[1:2,]

#plot(W);points(cases)

pointcase <- ppp(cases[,1], cases[,2], window=W)  #generate the ppp object

kdensity<-density.ppp(pointcase, 0.05)
plot(kdensity)
*Q:*there are almost the same density in the whole area,but in fact it may
have a higher density in the blue part of the attached map? I think the
problem may the inappropriate value of sigma, how to determine its value?

################################################################################
                                 Kernel density estimation--splancs
################################################################################
library(sp)
library(splancs)
library(foreign)
library(maptools)

case<-readShapePoints("d:/deleting/kernel/kernel/cases.shp")
guichi<-readShapePoly("d:/deleting/kernel/kernel/guichi2.shp")

#Conversion
case_pts <- coordinates(case)
case <- as.points(case_pts)
splancs_poly <-
getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(guichi)[[1]])[[1]])


#to unpack the coordinates of the points and the single ring boundary
polymap(splancs_poly,xlab="x(?)",ylab="y(?)")
pointmap(case_pts, add=TRUE)

m<-mse2d(case,splancs_poly,nsmse=1000,range=5)  #plots the estimated mean
square error as a function of h0
plot(m$h[290:1000],m$mse[290:1000],type="l")

n<-which(m$mse==min(m$mse))
h0<-m$h[n]


#smooth variation
smooth<-kernel2d(case, splancs_poly, h0=h0, nx=100, ny=100)
polymap(splancs_poly)  #sets the axes correctly and draws the polygon
image(smooth,add=T)   #the smoothed image is superimposed
polymap(splancs_poly,add=T) #redrawn the polymap in order not to be obsured
by smooth image
*Q*:The result is still not satisfied,there must be something wrong with my
programs.

        
On 6/25/07, Dan Putler <putler at sauder.ubc.ca> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Dan Putler
    

    

      
      #
    

  


  

      
Hi Zhi Jie,

Below are some changes to your code which should make you much  
happier. It does involve the use of the RColorBrewer package to  
create a color palette that makes the plot a bit easier to see. As  
you can guess, the size of the standard deviation given to  
density.ppp was the problem. Your data is in the Xian 1980/3-degree  
Gauss-Kruger CM 117E projection (which is EPSG code 2384). Based on  
looking at things, the units of this projection are meters. Using the  
defaults of density.ppp, the standard deviation of the bandwidth of  
the smoother is 0.05 meters, far smaller than you wanted. The  
implicit assumption in density.ppp appears to be that the window of  
the study area is a unit square. The width of your study area is  
about 60Km, so to get a comparable bandwidth for your study area  
relative to the unit square, I upped the standard deviation to 3000  
meters.

Here is my altered version of your code, you will need to change  
things back to the correct paths to the shapefile sets.
________
library(spatstat)
library(maptools)
library(RColorBrewer)

myPal=brewer.pal(12,"Paired") # An easily seen color palette

guichi<-readShapePoly("~/Research/data/guichi2.shp")
W <-as(as(guichi, "SpatialPolygons"), "owin")  #boundary polygons

cases<-coordinates(readShapePoints("~/Research/data/cases.shp"))   
#points
colnames(cases)<-c("x","y")
cases[1:2,]

#plot(W);points(cases)

pointcase <- ppp(cases[,1], cases[,2], window=W)  #generate the ppp  
object

kdensity<-density.ppp(pointcase, 3000)
plot(kdensity, col=myPal)

rm(list=c("guichi","W","cases","pointcase","kdensity","myPal"))
_______

Dan

        
On 25-Jun-07, at 7:49 AM, zhijie zhang wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Mon, 25 Jun 2007, Dan Putler wrote:

        
Thanks, Dan, I agree that the bandwidths are the problem here. Often the 
image breaks settings are also very misleading for 2D density plots, so 
one needs to take care, as you show.

Roger

            

      
      
      
    

        
Below are some changes to your code which should make you much 
happier. It does involve the use of the RColorBrewer package to 
create a color palette that makes the plot a bit easier to see. As 
you can guess, the size of the standard deviation given to 
density.ppp was the problem. Your data is in the Xian 1980/3-degree 
Gauss-Kruger CM 117E projection (which is EPSG code 2384). Based on 
looking at things, the units of this projection are meters. Using the 
defaults of density.ppp, the standard deviation of the bandwidth of 
the smoother is 0.05 meters, far smaller than you wanted. The 
implicit assumption in density.ppp appears to be that the window of 
the study area is a unit square. The width of your study area is 
about 60Km, so to get a comparable bandwidth for your study area 
relative to the unit square, I upped the standard deviation to 3000 
meters.

Here is my altered version of your code, you will need to change 
things back to the correct paths to the shapefile sets.
________
library(spatstat)
library(maptools)
library(RColorBrewer)

myPal=brewer.pal(12,"Paired") # An easily seen color palette

guichi<-readShapePoly("~/Research/data/guichi2.shp")
W <-as(as(guichi, "SpatialPolygons"), "owin")  #boundary polygons

cases<-coordinates(readShapePoints("~/Research/data/cases.shp")) 
#points
colnames(cases)<-c("x","y")
cases[1:2,]

#plot(W);points(cases)

pointcase <- ppp(cases[,1], cases[,2], window=W)  #generate the ppp 
object

kdensity<-density.ppp(pointcase, 3000)
plot(kdensity, col=myPal)

rm(list=c("guichi","W","cases","pointcase","kdensity","myPal"))
_______

Dan

        
On 25-Jun-07, at 7:49 AM, zhijie zhang wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo