Как извлечь оценку плотности ядра, сгенерированную функцией smooth_map? - PullRequest
Новая
       15

Как извлечь оценку плотности ядра, сгенерированную функцией smooth_map?

0 голосов
/ 05 июля 2019

Мне нравится вычислять оценки плотности ядра, используя функцию smooth_map пакета tmaptools.Я следую процедуре Криса Брансдона и Лекса Комбера из их учебника "Введение в R для пространственного анализа и картографирования".Из главы 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')

1 Ответ

0 голосов
/ 08 июля 2019

Кажется, я нашел то, что искал: 

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