Большую часть пространственной работы теперь можно выполнить довольно легко, используя пакет sf
.
Пример кода для аналогичной проблемы приведен ниже. Комментарии объясняют большую часть того, что он делает.
Трудной частью может быть понимание проекций карты (crs). Некоторые используют единицы (метры, футы и т. Д. c), а другие используют широту / долготу. Какой из них вы выберете, зависит от того, с какой частью земного шара вы работаете и чего пытаетесь достичь sh. Большинство веб-картографов используют crs 4326, но это не включает в себя легко используемое измерение расстояния.
На карте ниже показаны точки за пределами ~ 3 миль от Херефорда в красном цвете, а точки внутри - темно-бордовые. Синяя точка используется в качестве центра для Херефорда и буферной зоны.
library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
library(mapview)
set.seed(4)
#hereford approx location, ggmap requires api key
hereford <- data.frame(place = 'hereford', lat = -2.7160, lon = 52.0564) %>%
st_as_sf(coords = c('lat', 'lon')) %>% st_set_crs(4326)
#simulation of data points near-ish hereford
random_points <- data.frame(point_num = 1:20,
lat = runif(20, min = -2.8, max = -2.6),
lon = runif(20, min = 52, max = 52.1)) %>%
st_as_sf(coords = c('lat', 'lon')) %>% st_set_crs(4326) %>%st_transform(27700)
#make a buffer of ~3miles (4800m) around hereford
h_buffer <- hereford %>% st_transform(27700) %>% #change crs to one measured in meters
st_buffer(4800)
#only points inside ~3mi buffer
points_within <- random_points[st_within( random_points, h_buffer, sparse = F), ]
head(points_within)
#> Simple feature collection with 6 features and 1 field
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: 346243.2 ymin: 239070.3 xmax: 355169.8 ymax: 243011.4
#> CRS: EPSG:27700
#> point_num geometry
#> 1 1 POINT (353293.1 241673.9)
#> 3 3 POINT (349265.8 239397)
#> 4 4 POINT (349039.5 239217.7)
#> 6 6 POINT (348846.1 243011.4)
#> 7 7 POINT (355169.8 239070.3)
#> 10 10 POINT (346243.2 239690.3)
#shown in mapview
mapview(hereford, color = 'blue') +
mapview(random_points, color = 'red', legend = F, col.regions = 'red') +
mapview(h_buffer, legend = F) +
mapview(points_within, color = 'black', legend = F, col.regions = 'black')
Создано в 2020-04-12 пакетом представ. (v0.3.0)