Hi all,

I am trying to create a kernel density estimation map of points representing locations of severe weather events (attached script). The map and the event (point) data are imported from .shp files with LongLat projection. First the map and the point data were converted to a Albers Equal-Area Conic projection and then the kernel density was calculated using: density.ppp(ppp.regular, kernel = "gaussian", sigma=bw.diggle(ppp.regular)). The resulting map looks OK, but I have one question: what are the units of the density plot? Since I am using Albers Equal-Area Conic with units in meter, the resulting density should be events per square meters, but I am not sure if this is right.

Thank you very much for your help.

Cheers,
Bogdan


Bogdan Antonescu

Centre for Atmospheric Science
School of Earth, Atmospheric and Environmental Sciences
The University of Manchester
Simon Building, Oxford Road
Manchester, M13 9PL, United Kingdom
Room 3.16, Simon Building
+44 (0)161 306 3911

--------------------------------
library(sp)
library(maptools)
library(spatstat)
library(rgdal)
library(rgeos)

pointdata <- readShapePoints("~/KDE/points.shp", CRS("+proj=longlat"))
### link to data:
mapdata   <- readShapeLines("~/KDE/GSHHS_l_L1.shp", CRS("+proj=longlat"))
### link to data:

projstr <- "+proj=aea +lon_0=15e +lat_1=30n +lat_2=70n"

xrange <- c(-15,45)
yrange <- c(30,70)
bound <- cbind(xrange,yrange)
boundbox <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(xrange[2],seq(xrange[2],xrange[1],length=200),seq(xrange[1],xrange[2],length=200)),c(yrange[2],seq(yrange[1],yrange[1],length=200),seq(yrange[2],yrange[2],length=200))))),ID="bb")), proj4str=CRS("+proj=longlat")) # works for conical projections

mapdata.clipped <- gIntersection(mapdata, boundbox)
mapdata.proj    <- spTransform(mapdata.clipped, CRS(projstr))
data.clipped    <- gIntersection(pointdata, boundbox)
data.proj       <- spTransform(data.clipped, CRS(projstr))
boundbox.proj   <- spTransform(boundbox, CRS(projstr))

gxrange <- c(-15,45)
gyrange <- c(30,70)
xinterv <- 10
yinterv <- 10
xgrid <- seq(gxrange[1], gxrange[2], xinterv)
ygrid <- seq(gyrange[1], gyrange[2], yinterv)
grid1 <- gridlines(mapdata.clipped, easts=xgrid, norths=ygrid, ndiscr=100)
gridat1 <- gridat(mapdata.clipped, easts=xgrid, norths=ygrid)
grid1.proj    <- spTransform(grid1, CRS(projstr))
gridat1.proj  <- spTransform(gridat1, CRS(projstr))

windows()
par(mar=c(0.5,0.5,0.5,0.5)+0.1) # (bottom, left, top, right)
plot(mapdata.proj, xaxs="i", yaxs="i")
plot(grid1.proj, col="gray60", add=T)
text(coordinates(gridat1.proj), labels=parse(text=as.character(gridat1.proj$labels)), pos=gridat1.proj$pos, offset=-2, col="black")
plot(boundbox.proj, add=T, lwd=2)
w <- as.owin.SpatialPolygons(boundbox.proj)
ppp.regular <- ppp(data.proj@coords[,1], data.proj@coords[,2], window=w)
kde <- density.ppp(ppp.regular, kernel = "gaussian", sigma=bw.diggle(ppp.regular))
plot(kde, add=T)
plot(mapdata.proj, xaxs="i", yaxs="i", add=T)


	[[alternative HTML version deleted]]

