R - Подгонка сетки к карте города и ввод данных в квадраты сетки

Question

Я пытаюсь разместить сетку над Сан-Хосе следующим образом:

Сетка Сан-Хосе

Вы можете сделать сетку визуально, используя следующий код:

  ca_cities = tigris::places(state = "CA") #using tigris package to get shape file of all CA cities

  sj = ca_cities[ca_cities$NAME == "San Jose",] #specifying to San Jose

  UTM_ZONE = "10" #the UTM zone for San Jose, will be used to convert the proj4string of sj into UTM

  main_sj = sj@polygons[[1]]@Polygons[[5]] #the portion of the shape file I focus on. This is the boundary of san jose

  #converting the main_sj polygon into a spatialpolygondataframe using the sp package
  tst_ps = sp::Polygons(list(main_sj), 1)
  tst_sps = sp::SpatialPolygons(list(tst_ps))
  proj4string(tst_sps) = proj4string(sj)
  df = data.frame(f = 99.9)
  tst_spdf = sp::SpatialPolygonsDataFrame(tst_sps, data = df)

  #transforming the proj4string and declaring the finished map as "map"
  map = sp::spTransform(tst_sps, CRS(paste0("+proj=utm +zone=",UTM_ZONE," ellps=WGS84")))

  #designates the number of horizontal and vertical lines of the grid
  NUM_LINES_VERT = 25
  NUM_LINES_HORZ = 25
  #getting bounding box of map
  bbox = map@bbox
  #Marking the x and y coordinates for each of the grid lines.
  x_spots = seq(bbox[1,1], bbox[1,2], length.out = NUM_LINES_HORZ)
  y_spots = seq(bbox[2,1], bbox[2,2], length.out = NUM_LINES_VERT)

  #creating the coordinates for the lines. top and bottom connect to each other. left and right connect to each other
  top_vert_line_coords = expand.grid(x = x_spots, y = y_spots[1])
  bottom_vert_line_coords = expand.grid(x = x_spots, y = y_spots[length(y_spots)])
  left_horz_line_coords = expand.grid(x = x_spots[1], y = y_spots)
  right_horz_line_coords = expand.grid(x = x_spots[length(x_spots)], y = y_spots)

  #creating vertical lines and adding them all to a list
  vert_line_list = list()
  for(n in 1 : nrow(top_vert_line_coords)){
    vert_line_list[[n]] = sp::Line(rbind(top_vert_line_coords[n,], bottom_vert_line_coords[n,]))
  }

  vert_lines = sp::Lines(vert_line_list, ID = "vert") #creating Lines object of the vertical lines

  #creating horizontal lines and adding them all to a list
  horz_line_list = list()
  for(n in 1 : nrow(top_vert_line_coords)){
    horz_line_list[[n]] = sp::Line(rbind(left_horz_line_coords[n,], right_horz_line_coords[n,]))
  }

  horz_lines = sp::Lines(horz_line_list, ID = "horz") #creating Lines object of the horizontal lines

  all_lines = sp::Lines(c(horz_line_list, vert_line_list), ID = 1) #combining horizontal and vertical lines into a single grid format

  grid_lines = sp::SpatialLines(list(all_lines)) #converting the lines object into a Spatial Lines object
  proj4string(grid_lines) = proj4string(map) #ensuring the projections are the same between the map and the grid lines.

  trimmed_grid = intersect(grid_lines, map) #grid that shapes to the san jose map

  plot(map) #plotting the map of San Jose
  lines(trimmed_grid) #plotting the grid

Однако я изо всех сил пытаюсь превратить каждую «сетку» каждой сетки (некоторые части сетки не являются квадратами, поскольку они соответствуют форме карты Сан-Хосе) в корзину, в которую я мог бы вводить данныев.Другими словами, если бы каждый квадрат «сетки» был пронумерован 1: n, то я мог бы создать такой кадр данных:

  grid_id num_assaults num_thefts
1       1          100         89
2       2           55        456
3       3           12       1321
4       4           48        498
5       5           66          6

и заполнить каждый квадрат «сетки» данными точечной локализацией каждого случая совершения преступления.Надеюсь, что с помощью функции over() из пакета sp.

Я пытался решить эту проблему в течение нескольких недель и не могу ее решить.Я искал простое решение, но я не могу найти его.Любая помощь будет оценена.

Answer 1

Кроме того, вот решение на основе sf и tidyverse:

С помощью sf вы можете создать сетку квадратов с помощью функции st_make_grid().Здесь я сделаю 2-километровую сетку над ограничительной рамкой Сан-Хосе, а затем пересекаю ее с границей Сан-Хосе.Обратите внимание, что я проецируюсь на UTM-зону 10N, поэтому я могу указать размер сетки в метрах.

library(tigris)
library(tidyverse)
library(sf)
options(tigris_class = "sf", tigris_use_cache = TRUE)
set.seed(1234)

sj <- places("CA", cb = TRUE) %>%
  filter(NAME == "San Jose") %>%
  st_transform(26910)

g <- sj %>%
  st_make_grid(cellsize = 2000) %>%
  st_intersection(sj) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(id = row_number())

Далее мы можем сгенерировать некоторые случайные данные о преступлении с помощью st_sample() и нанести их на график, чтобы увидеть, с чем мы работаем.

thefts <- st_sample(sj, size = 500) %>%
  st_sf()

assaults <- st_sample(sj, size = 200) %>%
  st_sf()

plot(g$geometry)
plot(thefts, add = TRUE, col = "red")

Затем данные о преступности можно пространственно присоединить к сетке с помощью st_join().Мы можем подготовить график, чтобы проверить наши результаты.

theft_grid <- g %>%
  st_join(thefts) %>%
  group_by(id) %>%
  summarize(num_thefts = n())

plot(theft_grid["num_thefts"])

Затем мы можем сделать то же самое с данными нападений, а затем соединить два набора данных вместе, чтобы получить желаемый результат.Если у вас было много наборов криминальных данных, их можно изменить, чтобы они работали в пределах некоторого изменения purrr::map().

assault_grid <- g %>%
  st_join(assaults) %>%
  group_by(id) %>%
  summarize(num_assaults = n()) 

st_geometry(assault_grid) <- NULL

crime_data <- left_join(theft_grid, assault_grid, by = "id")

crime_data

Simple feature collection with 190 features and 3 fields
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 584412 ymin: 4109499 xmax: 625213.2 ymax: 4147443
epsg (SRID):    26910
proj4string:    +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
# A tibble: 190 x 4
      id num_thefts num_assaults                                                     geometry
   <int>      <int>        <int>                                               <GEOMETRY [m]>
 1     1          2            1 POLYGON ((607150.3 4111499, 608412 4111499, 608412 4109738,…
 2     2          4            1 POLYGON ((608412 4109738, 608412 4111499, 609237.8 4111499,…
 3     3          3            1 POLYGON ((608412 4113454, 608412 4111499, 607150.3 4111499,…
 4     4          2            2 POLYGON ((609237.8 4111499, 608412 4111499, 608412 4113454,…
 5     5          1            1 MULTIPOLYGON (((610412 4112522, 610412 4112804, 610597 4112…
 6     6          1            1 POLYGON ((616205.4 4113499, 616412 4113499, 616412 4113309,…
 7     7          1            1 MULTIPOLYGON (((617467.1 4113499, 618107.9 4113499, 617697.…
 8     8          2            1 POLYGON ((605206.8 4115499, 606412 4115499, 606412 4114617,…
 9     9          5            1 POLYGON ((606412 4114617, 606412 4115499, 608078.2 4115499,…
10    10          1            1 POLYGON ((609242.7 4115499, 610412 4115499, 610412 4113499,…
# ... with 180 more rows

Answer 2

Если вашей целью является только визуальный, а не обязательно весь код и данные агрегирования сетки, вы можете сгенерировать интерактивную карту и сетку в library(mapdeck) (отметив, что вам понадобится токен доступа Mapbox)

Первый шаг для генерации данных заимствован из ответа @kwalkertcu

library(tigris)
library(sf)
options(tigris_class = "sf", tigris_use_cache = TRUE)
set.seed(1234)

sj <- places("CA", cb = TRUE) %>%
  filter(NAME == "San Jose") %>%
  st_transform(26910)

thefts <- st_sample(sj, size = 500) %>%
  st_sf() %>%
  st_transform(crs = 4326)

## some random weight data
thefts$weight <- sample(1:100, size = nrow(thefts), replace = T)

Затем, учитывая объект sf со столбцом weight, вы можете построить его, используяadd_screengrid()

library(mapdeck)

set_token("MAPBOX_TOKEN")

mapdeck(
  style = mapdeck_style("dark")
  , location = c(-121.8, 37.3)
  , zoom = 6
) %>%
  add_screengrid(
    data = thefts
    , cell_size = 15
    , weight = "weight"
  )

Примечания:

• Я использую GitHub-версию mapdeck, где APIизменилось незначительно, но версия CRAN должна дать тот же результат.

Answer 3

Используя объект Spatial * в качестве ваших данных

library(tigris)
ca_cities = tigris::places(state = "CA") #using tigris package to get shape file of all CA cities
sj = ca_cities[ca_cities$NAME == "San Jose",] #specifying to San Jose
sjutm = sp::spTransform(sj, CRS("+proj=utm +zone=10 +datum=WGS84"))

Вы можете создать сетку из многоугольников, подобных этой

library(raster)
r <- raster(sjutm, ncol=25, nrow=25)
rp <- as(r, 'SpatialPolygons')

Показать ее

plot(sjutm, col='red')
lines(rp, col='blue')

Чтобы подсчитать количество наблюдений на ячейку сетки (используя несколько случайных точек здесь), вы не хотите использовать полигоны, а скорее RasterLayer

set.seed(0)
x <- runif(500, xmin(r), xmax(r))
y <- runif(500, ymin(r), ymax(r))
xy1 <- cbind(x, y)
x <- runif(500, xmin(r), xmax(r))
y <- runif(500, ymin(r), ymax(r))
xy2 <- cbind(x, y)

d1 <- rasterize(xy1, r, fun="count", background=0)
d2 <- rasterize(xy2, r, fun="count", background=0)

plot(d1)
plot(sjutm, add=TRUE)

, за которым следует

s <- stack(d1, d2)
names(s) = c("assault", "theft")
s <- mask(s, sjutm)
plot(s, addfun=function()lines(sjutm))

.таблица, за которой вы следите

p <- rasterToPoints(s)
cell <- cellFromXY(s, p[,1:2])
res <- data.frame(grid_id=cell, p[,3:4])
head(res)
#  grid_id assault theft
#1       1       1     1
#2       2       0     1
#3       3       0     3
#4       5       1     1
#5       6       1     0
#6      26       0     0

Вы также можете создать SpatialPolygonsDataFrame из результатов

pp <- as(s, 'SpatialPolygonsDataFrame')
pp
#class       : SpatialPolygonsDataFrame 
#features    : 190 
#extent      : 584411.5, 623584.9, 4109499, 4147443  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=utm +zone=10 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#variables   : 2
#names       : assault, theft 
#min values  :       0,     0 
#max values  :       4,     5