Плотность населения в пределах полигонов - PullRequest
Новая
       70

Плотность населения в пределах полигонов

1 голос
/ 11 февраля 2020

Итак, у меня есть несколько вопросов, касающихся пакета растров в R. У меня есть растр с оценочной совокупностью в каждой точке сетки. У меня также есть шейп-файл с полигонами областей. Я хочу узнать координаты района с самой высокой плотностью населения в каждом регионе. Предположим, что каждая окрестность представляет собой однородный квадрат 5 на 5 точек сетки.

Следующий игрушечный пример имитирует мою проблему. 

library(raster)
library(maptools)

set.seed(123)

data(wrld_simpl)

wrld_simpl <- st_as_sf(wrld_simpl)
contr_c_am <- wrld_simpl %>%
  filter(SUBREGION ==13) %>%
  filter(FIPS != "MX") %>%
  select(NAME) 


# Create a raster of population (sorry for the bad example spatial distribution)
r <- raster(xmn=-180, xmx=180, ymn=-90, ymx=90, res=0.1)
values(r) <- runif(ncell(r), 0, 100)

# keep only raster around the region of interest
r_small <- crop(r, extent(contr_c_am))
plot(r_small)
plot(st_geometry(contr_c_am), add = T)

raster_contr_c_am <- rasterize(contr_c_am, r)

raster_contr_c_am - это сетка населения, а название региона сохраняется в качестве атрибута.

Каким-то образом мне нужно отфильтровать только точки сетки из одного региона и, возможно, использовать какую-то функцию, например, focal (), чтобы найти общую численность населения поблизости. 

focal(raster_contr_c_am, matrix(1,5,5),sum, pad = T, padValue = 0)

Затем мне нужно найти, какую сетку точка имеет наибольшее значение в каждой области и сохраняет свои координаты.

Надеюсь, мое объяснение не слишком запутанное,

Спасибо за любую помощь!

1 Ответ

2 голосов
/ 20 февраля 2020

Вот пример, который перебирает фигуру, определяющую область, затем использует растровые значения внутри области и функцию focal(), чтобы найти максимум. 

library(raster)
library(maptools)
library(sf)
library(dplyr)

set.seed(123)

data(wrld_simpl)

wrld_simpl <- st_as_sf(wrld_simpl)
contr_c_am <- wrld_simpl %>%
    filter(SUBREGION ==13) %>%
    filter(FIPS != "MX") %>%
    select(NAME) 

# Create a raster of population (sorry for the bad example spatial distribution)
r <- raster(xmn=-180, xmx=180, ymn=-90, ymx=90, res=0.1)
values(r) <- runif(ncell(r), 0, 100)

# keep only raster around the region of interest
r_small <- crop(r, extent(contr_c_am))
raster_contr_c_am <- rasterize(contr_c_am, r_small)

# function to find the max raster value using focal
# in a region 
findMax <- function(region, raster) {
    tt <- trim((mask(raster, region))) # focus on the region
    ff <- focal(tt, w=matrix(1/25,nc=5,nr=5))
    maximumCell <- which.max(ff) # find the maximum cell id
    maximumvalue <- maxValue(ff) # find the maximum value
    maximumx <- xFromCell(ff, maximumCell) # get the coordinates
    maximumy <- yFromCell(ff, maximumCell)
    # return a data frame
    df <- data.frame(maximumx, maximumy, maximumvalue)
    df
}

numberOfShapes <- nrow(contr_c_am)
ll <- lapply(1:numberOfShapes, function(s) findMax(region = contr_c_am[s,], raster = r_small))
merged <- do.call(rbind, ll)
maxpoints <- st_as_sf(merged, coords=c('maximumx', 'maximumy'), crs=crs(contr_c_am))

library(mapview) # optional but nice visualization - select layers to see if things look right
mapview(maxpoints) + mapview(r_small) + mapview(contr_c_am)

Я сделал sf объект, так что он может быть нанесен с другими пространственными объектами. Используя пакет mapview, я получаю это.

enter image description here

...