用网格覆盖空间多边形，并使用R检查在其中定位特定网格元素的坐标

# Load packages 
library(maptools)

# Download shapefile for NYC
# OLD URL (no longer working)
# shpurl <- "http://www.nyc.gov/html/dcp/download/bytes/nybb_13a.zip"
shpurl <- "https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nybb_13a.zip"

tmp    <- tempfile(fileext=".zip")
download.file(shpurl, destfile=tmp)
files <- unzip(tmp, exdir=getwd())

# Load & plot shapefile
shp <- readShapePoly(files[grep(".shp$", files)])
plot(shp)

# Define coordinates 
points_of_interest <- data.frame(y=c(919500, 959500, 1019500, 1049500, 1029500, 989500), 
                 x =c(130600, 150600, 180600, 198000, 248000, 218000),
                 id  =c("A"), stringsAsFactors=F)

# Plot coordinates
points(points_of_interest$y, points_of_interest$x, pch=19, col="red")

#1 楼

这是一个使用SpatialGrid对象的示例：

### read shapefile
library("rgdal")
shp <- readOGR("nybb_13a", "nybb")

proj4string(shp)  # units us-ft
# [1] "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 
# +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
# +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"

### define coordinates and convert to SpatialPointsDataFrame
poi <- data.frame(x=c(919500, 959500, 1019500, 1049500, 1029500, 989500),
                  y=c(130600, 150600, 180600, 198000, 248000, 218000),
                  id="A", stringsAsFactors=F)
coordinates(poi) <- ~ x + y
proj4string(poi) <- proj4string(shp)

### define SpatialGrid object
bb <- bbox(shp)
cs <- c(3.28084, 3.28084)*6000  # cell size 6km x 6km (for illustration)
                                # 1 ft = 3.28084 m
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd
# cellcentre.offset 923018 129964
# cellsize           19685  19685
# cells.dim              8      8

sp_grd <- SpatialGridDataFrame(grd,
                               data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(shp)))
summary(sp_grd)
# Object of class SpatialGridDataFrame
# Coordinates:
#      min     max
# x 913175 1070655
# y 120122  277602
# Is projected: TRUE
# ...

现在您可以使用已实现的over-方法获取单元格ID：

要绘制带有单元格ID的shapefile和网格，请执行以下操作：彩色键：

over(poi, sp_grd)
#   id
# 1 57
# 2 51
# 3 38
# 4 39
# 5 14
# 6 28

#2 楼

问题中提供的纽约数据集不再可供下载。我使用sf包中的nc数据集来演示使用sf包的解决方案：

#3 楼

如果您没有看过R栅格数据包，则它具有用于与矢量GIS对象之间进行转换的工具，因此您应该能够a）创建具有200x200m像元的栅格（网格），以及b）使用以下方法将其转换为一组多边形某种逻辑ID。从那里，我将看一下sp程序包，以帮助将点和多边形网格相交。此http://cran.r-project.org/web/packages/sp/vignettes/over.pdf页面可能是一个不错的开始。在sp包文档中徘徊，您也许可以从SpatialGrid类开始，而完全跳过栅格部分。

#4 楼

“ GIS世界”很复杂，并且有许多必须符合您的数据标准。所有“ GIS工具”均按GIS标准进行互操作。今天（2014年）的所有“严重GIS数据”都存储在数据库中。

在SQL环境中，与其他FOSS工具一起，在GIS上下文中“使用R”的最佳方法已嵌入。最好的工具是PostgreSQL 9.X（请参阅PL / R）和PostGIS。


您回答：


导入/导出形状文件：请使用shp2pgsql和pgsql2shp。
要“将形状文件分割为200米见方/子多边形”：请参见ST_SnapToGrid()和ST_AsRaster()等。我们需要更好地了解您表达为“食谱”的需求。 >您说需要“定位地理坐标” ..也许ST_Centroid()是正方形（？）...您可以表达“更多数学式”，所以我理解。任何栅格转换，仅是规则采样点的矩阵。


一种原始方法是在通常的外部编译器中使用不带PL / R的R：仅转换多边形并将其导出为形状或为WKT（请参阅ST_AsText），然后使用awk或其他过滤器将数据转换为R格式。