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

Question

我喜欢使用tmaptools包的smooth_map函数来计算内核密度估计。我遵循Chris Brunsdon和Lex Comber在他们的教科书“用于空间分析和绘图的R入门”中的过程。我从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')

Answer 1

我想我找到了我正在寻找的东西：

# 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]