如何提取 smooth_map 函数生成的核密度估计值?

How can I extract the kernel density estimate generated by the smooth_map function?

我喜欢使用 tmaptools 包的 smooth_map 函数来计算内核密度估计值。我遵循 Chris Brunsdon 和 Lex Comber 教科书中的程序 "An introduction to R for spatial analysis and mapping"。来自 chapter 6 我改编了下面的代码示例。

我想为构成 smooth_map 的每个坐标对提取核密度估计,但到目前为止无法实现。

# Load GISTools (for the data) and tmap (for the mapping)

require(GISTools)
require(tmap)
require(tmaptools)

# Get the data
data(newhaven)
# look at it
# select 'view' mode
tmap_mode('view')

# Create the map of blocks and incidents
tm_shape(blocks) + tm_borders() + tm_shape(breach) +
  tm_dots(col='navyblue')

# Function to choose bandwidth according Bowman and Azzalini / Scott's rule
# for use with <smooth_map> in <tmaptools>

choose_bw <- function(spdf) {
  X <- coordinates(spdf)
  sigma <- c(sd(X[,1]),sd(X[,2]))  * (2 / (3 * nrow(X))) ^ (1/6)
  return(sigma/1000)
}

# Calculate kernel density

breach_dens <- smooth_map(breach,cover=blocks, bandwidth = choose_bw(breach))

# Plot resulting polygon map

tm_shape(breach_dens$polygons) +  tm_fill(col='level',alpha=0.8)+
  tm_shape(blocks) + tm_borders() + tm_shape(breach) +
  tm_dots(col='navyblue')

我想我找到了我要找的东西:

# Transform to sp spatialpoints data frame with common projection

  breach <- spTransform(breach, CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")) 

# Extract polygon from smooth_map object and transform to sp spatial polygon dataframe with common projection

  breach_poly <- as(breach_dens$polygons, 'Spatial')
  breach_poly <- spTransform(breach_poly, CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")) 

# Generate export dataframe

  # extract point coordinates

    export <- data.frame(coordinates(breach))

  # extract kernel density estimate for each coordinate

    export$kde <-  over(breach,breach_poly, returnList = FALSE)

  # Save vertical and horizontal bandwiths

    export$bw_x <- breach_dens$bandwidth[1:1]
    export$bw_y <- breach_dens$bandwidth[2:2]