如何提取 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]
我喜欢使用 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]