[R-sig-Geo] functions on kernel density estimation in pointpattern analysis

Dan Putler putler at sauder.ubc.ca
Mon Jun 25 18:35:49 CEST 2007 

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:

> 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: 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:
>
> > I rather like locfit.
> >
> > THK
> >
> > On 6/24/07, zhijie zhang <epistat at gmail.com> wrote:
> >> Dear Friends,
> >>  Except kernel2d(splancs) function, are there any other functions
> >> on kernel
> >> density estimation in point pattern analysis? I use the kernel2d
> >> (splancs)
> >> function to anayze my dataset, and the result seems not to be very
> >> good.
> >>   Any suggestions or help in kernel density estimation of
> >> univariate or
> >> multivariate point process are greatly appreciated.
> >>   BTW, i mainly want ot do the kernel density estimation in both
> >> univariate and multivariate point process.
> >> --
> >> With Kind Regards,
> >>
> >> oooO:::::::::
> >> (..):::::::::
> >> :\.(:::Oooo::
> >> ::\_)::(..)::
> >> :::::::)./:::
> >> ::::::(_/::::
> >> :::::::::::::
> >>  
> [********************************************************************
> >> ***]
> >> Zhi Jie,Zhang ,PHD
> >> Tel:86-21-54237149
> >> Dept. of Epidemiology,School of Public Health,Fudan University
> >> Address:No. 138 Yi Xue Yuan Road,Shanghai,China
> >> Postcode:200032
> >> Email:epistat at gmail.com
> >> Website: www.statABC.com
> >>  
> [********************************************************************
> >> ***]
> >> oooO:::::::::
> >> (..):::::::::
> >> :\.(:::Oooo::
> >> ::\_)::(..)::
> >> :::::::)./:::
> >> ::::::(_/::::
> >> :::::::::::::
> >>
> >>         [[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
> >>
> >
> >
> > --
> > Timothy H. Keitt, University of Texas at Austin
> > Contact info and schedule at http://www.keittlab.org/tkeitt/
> > Reprints at http://www.keittlab.org/tkeitt/papers/
> > ODF attachment? See http://www.openoffice.org/
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
>
> -- 
> With Kind Regards,
>
> oooO:::::::::
> (..):::::::::
> :\.(:::Oooo::
> ::\_)::(..)::
> :::::::)./:::
> ::::::(_/::::
> :::::::::::::
> [********************************************************************* 
> **]
> Zhi Jie,Zhang ,PHD
> Tel:86-21-54237149
> Dept. of Epidemiology,School of Public Health,Fudan University
> Address:No. 138 Yi Xue Yuan Road,Shanghai,China
> Postcode:200032
> Email:epistat at gmail.com
> Website: www.statABC.com
> [********************************************************************* 
> **]
> oooO:::::::::
> (..):::::::::
> :\.(:::Oooo::
> ::\_)::(..)::
> :::::::)./:::
> ::::::(_/::::
> :::::::::::::
> <data.rar>
> <map.jpg>

More information about the R-sig-Geo mailing list