1

Rookie R 用户在这里，我将非常感谢您能给我的任何帮助。

我的项目要求我在我选择的城市周围创建一个矢量边界框，然后过滤大量数据，因此我只有与该区域相关的数据。然而，我使用 R studio 已经有好几年了，公平地说，我对这门语言几乎一无所知。

我最初使用

geocode("Hereford, UK")

bbox <-c(Longitude=-2.72,Latitude=52.1)

myMap <- get_map(location = "Hereford, UK",source="google",maptype="roadmap")

然后，我必须创建一个新的 tibble，它会过滤掉并仅将相关数据提供给该区域。

我不确定如何进行此操作，然后我必须将数据覆盖到我创建的地图上。

由于我只有一个坐标中心点，是否可以在我的位置中心周围创建一个半径为 3 英里的圆，以便我可以过滤这个区域？

感谢大家花时间阅读我的帖子。干杯!

4

1 回答 1

0

大多数空间工作现在可以使用该sf软件包轻松完成。

类似问题的示例代码如下。评论解释了它的大部分功能。

困难的部分可能在于理解地图投影（crs）。有些使用单位（米、英尺等），有些使用纬度/经度。您选择哪一个取决于您正在与地球的哪个区域合作以及您要完成的工作。大多数网络地图使用 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')

由reprex 包于 2020-04-12 创建(v0.3.0)

于 2020-04-12T13:12:32.490 回答