How can you convert the output of a kernel density output function, for example bkde2D, into units per square km?

I.e. there must be a relationship between the bandwidth, gridsize and real density, or maybe the fhat value? Does the projection also impact the "real" density?

For example, in the code below, it would be nice to know that the range in the legend (~0-50) was 0 to 50 thefts per square km.


inurl <- "https://data.cityofchicago.org/api/views/22s8-eq8h/rows.csv?accessType=DOWNLOAD"
infile <- "mvthefts.csv"

# Load and clean up variable names, and convert dates
download.file(url = inurl, destfile = infile)
dat <- data.table::fread(infile)
setnames(dat, tolower(colnames(dat)))
setnames(dat, gsub(" ", "_", colnames(dat)))
dat <- dat[!is.na(longitude)]
dat[ , date := as.IDate(date, "%m/%d/%Y")]

# Create kernel density output
kde <- bkde2D(dat[ , list(longitude, latitude)],
          bandwidth=c(.0045, .0068), gridsize = c(1000,1000))

# Create Raster from Kernel Density output
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))

# Set low density cells as NA so we can make them transparent with the colorNumeric function
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA

# Create pal function for coloring the raster
palRaster <- colorNumeric("YlOrRd", domain = KernelDensityRaster@data@values, na.color = "transparent")

# Draw on a map
leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
             colors = palRaster, 
             opacity = .8) %>%
  addLegend(pal = palRaster, 
        values = KernelDensityRaster@data@values, 
        title = "Kernel density of points")

Output of MVE

  • Please read up on Kernel Density Estimates. What you are describing is a simple point density [n/area] and not KDE. – Jeffrey Evans Feb 14 at 17:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.