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

【问题标题】：How can I extract the kernel density estimate generated by the smooth_map function?如何提取由 smooth_map 函数生成的核密度估计？
【发布时间】：2019-07-05 16:32:06
【问题描述】：

我喜欢使用 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')

【问题讨论】：

标签： r dictionary smoothing

【解决方案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]

【讨论】：