r - 来自多边形输入的点网格
r - Grid of points from polygon input
我正在编写一个脚本，该脚本采用在 google earth 中创建的输入 KML 文件，并在多边形内部绘制坐标点网格。
到目前为止，我有多边形输入和多边形边界框的点网格，但我只想在多边形内部有点。
我尝试使用 over() 函数执行此操作，但它不起作用。建议？
您可以下载我的测试 KML 文件 HERE。
library(rgdal)
library(sp)
library(maptools)
# ogrInfo() to find layer name... not as labelled in Google Earth?!
my.poly = readOGR(ds = "PolyNYC.kml", layer = "PolyNYC")
proj4string(my.poly) <- "+proj=longlat +datum=WGS84 +no_defs"
# Creating grid of points
grdpts <- makegrid(my.poly)
# Converting from df to spdf
coords = cbind(grdpts$x1, grdpts$x2)
sp = SpatialPoints(coords)
spdf = SpatialPointsDataFrame(coords, grdpts, proj4string = CRS(proj4string(my.poly)))
# Using over() to select only those points in the polygon
inPoly = over(spdf, my.poly)
# This is not working
# Plotting the polygon with the points overlaid.
plot(my.poly)
points(spdf, pch = 3, col = "red")
#kmlPoints(obj = spdf, kmlfile = "BBoxFromPoly.kml", kmlname = "Testing123")
我将展示一个使用 library(sf) 的解决方案，它是 library(sp)
的继承者
正在加载数据
library(sf)
## read the kml
my.poly <- sf::st_read("~/Downloads/PolyNYC.kml")
## create a grid of points
grdpts <- sf::st_make_grid(my.poly, what = "centers")
## convert it to an `sf` object, as opposed to an `sfc`
my.points <- sf::st_sf(grdpts)
查看数据
要查看地图上的对象，我使用我的 googleway 包将其绘制在 Google 地图上（因此您需要一个 API 键），但您可以使用 leaflet 或者任何你想要的地图
library(googleway)
set_key("your_api_key_here")
google_map() %>%
add_polygons(my.poly) %>%
add_markers(my.points)
多边形中的点
您可以使用函数sf::st_join()连接几何
pointsInside <- sf::st_join(x = my.points, y = my.poly, left = FALSE)
# Simple feature collection with 59 features and 2 fields
# geometry type: POINT
# dimension: XY
# bbox: xmin: -74.1754 ymin: 40.63513 xmax: -73.75675 ymax: 40.8514
# epsg (SRID): 4326
# proj4string: +proj=longlat +datum=WGS84 +no_defs
# First 10 features:
# Name Description geometry
# 1 TestLayerNYC POINT (-74.08237 40.63513)
# 2 TestLayerNYC POINT (-74.03585 40.63513)
# 3 TestLayerNYC POINT (-74.1754 40.65916)
# 4 TestLayerNYC POINT (-74.12889 40.65916)
# 5 TestLayerNYC POINT (-74.08237 40.65916)
# 6 TestLayerNYC POINT (-74.03585 40.65916)
# 7 TestLayerNYC POINT (-73.80326 40.65916)
# 8 TestLayerNYC POINT (-74.1754 40.68319)
# 9 TestLayerNYC POINT (-74.12889 40.68319)
# 10 TestLayerNYC POINT (-74.08237 40.68319)
这里，pointsInside是多边形内的所有点
查看结果
google_map() %>%
add_polygons(my.poly) %>%
add_markers(pointsInside)
我正在编写一个脚本，该脚本采用在 google earth 中创建的输入 KML 文件，并在多边形内部绘制坐标点网格。
到目前为止，我有多边形输入和多边形边界框的点网格，但我只想在多边形内部有点。
我尝试使用 over() 函数执行此操作，但它不起作用。建议？
您可以下载我的测试 KML 文件 HERE。
library(rgdal)
library(sp)
library(maptools)
# ogrInfo() to find layer name... not as labelled in Google Earth?!
my.poly = readOGR(ds = "PolyNYC.kml", layer = "PolyNYC")
proj4string(my.poly) <- "+proj=longlat +datum=WGS84 +no_defs"
# Creating grid of points
grdpts <- makegrid(my.poly)
# Converting from df to spdf
coords = cbind(grdpts$x1, grdpts$x2)
sp = SpatialPoints(coords)
spdf = SpatialPointsDataFrame(coords, grdpts, proj4string = CRS(proj4string(my.poly)))
# Using over() to select only those points in the polygon
inPoly = over(spdf, my.poly)
# This is not working
# Plotting the polygon with the points overlaid.
plot(my.poly)
points(spdf, pch = 3, col = "red")
#kmlPoints(obj = spdf, kmlfile = "BBoxFromPoly.kml", kmlname = "Testing123")
我将展示一个使用 library(sf) 的解决方案，它是 library(sp)
正在加载数据
library(sf)
## read the kml
my.poly <- sf::st_read("~/Downloads/PolyNYC.kml")
## create a grid of points
grdpts <- sf::st_make_grid(my.poly, what = "centers")
## convert it to an `sf` object, as opposed to an `sfc`
my.points <- sf::st_sf(grdpts)
查看数据
要查看地图上的对象，我使用我的 googleway 包将其绘制在 Google 地图上（因此您需要一个 API 键），但您可以使用 leaflet 或者任何你想要的地图
library(googleway)
set_key("your_api_key_here")
google_map() %>%
add_polygons(my.poly) %>%
add_markers(my.points)
多边形中的点
您可以使用函数sf::st_join()连接几何
pointsInside <- sf::st_join(x = my.points, y = my.poly, left = FALSE)
# Simple feature collection with 59 features and 2 fields
# geometry type: POINT
# dimension: XY
# bbox: xmin: -74.1754 ymin: 40.63513 xmax: -73.75675 ymax: 40.8514
# epsg (SRID): 4326
# proj4string: +proj=longlat +datum=WGS84 +no_defs
# First 10 features:
# Name Description geometry
# 1 TestLayerNYC POINT (-74.08237 40.63513)
# 2 TestLayerNYC POINT (-74.03585 40.63513)
# 3 TestLayerNYC POINT (-74.1754 40.65916)
# 4 TestLayerNYC POINT (-74.12889 40.65916)
# 5 TestLayerNYC POINT (-74.08237 40.65916)
# 6 TestLayerNYC POINT (-74.03585 40.65916)
# 7 TestLayerNYC POINT (-73.80326 40.65916)
# 8 TestLayerNYC POINT (-74.1754 40.68319)
# 9 TestLayerNYC POINT (-74.12889 40.68319)
# 10 TestLayerNYC POINT (-74.08237 40.68319)
这里，pointsInside是多边形内的所有点
查看结果
google_map() %>%
add_polygons(my.poly) %>%
add_markers(pointsInside)