映射点和多边形

Question

我是 R 初学者，第一次处理 R 和空间数据。所以我希望我说清楚。

我有一个意大利地区的 shapefile，分为几个人口普查区。另一方面，我有一个 csv 文件，其中包含每个案例的案例列表和地址。我想映射点和 shapefile 并计算每个人口普查部门有多少点。

这是我从现在开始所做的：

#get cases file
cases <- read.csv("cases.csv", sep =';', header = TRUE)
names(cases)
[1] "name"    "address"

#geocode addresses from Google Maps
library(GISTools)
library(rgeos)
library(ggmap)
geolocalize <- geocode(as.character(cases$address))

# bind latitude and longitude to the previous cases data frame
cases <- data.frame(cases, geolocalize)
names(cases)
[1] "name"    "address" "lon"     "lat"

#make cases a SpatialPointDataFrame
#since addresses were retrieved using GoogleMaps, I set proj4string as follows
cases.points <- SpatialPointsDataFrame(cases[,3:4], cases, proj4string = CRS("+init=EPSG:3857"))

#get the shapefile
region <- readOGR("R02_11_WGS84.shp")

现在，我可以分别绘制 cases.points 和 shapefile，但不能将它们添加到同一图中。除此之外，正如我所说，我想计算每个多边形中有多少个点（即 "region" 的人口普查区）。

我必须承认我对地理不是很感兴趣。我怀疑不同的坐标 and/or 投影参考系统可能是问题所在，因此我检查过。

head(coordinates(region))
[,1]    [,2]
0 364509.0 5065900
1 363916.3 5056629
2 372585.0 5068078
3 360692.3 5048321
4 356029.7 5062399
5 360012.1 5065663


coordinates(cases)
lon      lat
[1,] 7.323667 45.73664

proj4string(region)
[1] "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

proj4string(cases.points)
[1] "+init=EPSG:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"

是否有可能 shapefile 坐标是度数，而案例坐标是十进制？如果可以，如何转换？

谢谢，萨罗

Answer 1

您有两个不同的投影 - 您的 "region" shapefile 投影在 UTM zone 32 中并且您告诉 R 使用 Web Mercator for "cases"。但是，如果您只是从 Google 下载案例数据 lat/long，那么您不想告诉 R 他们的投影是 Web Mercator，因为它不是 - 它是未投影的 WGS 84，所以您我想要 EPSG 4326。Web Mercator 是 Google 用来显示地图的工具，但如果您下载 lat/long，那只是未投影的坐标。要正确读取 lat/long 数据，请使用：

library(sp)
cases.points <- SpatialPointsDataFrame(cases[,3:4], cases, 
                                   proj4string = CRS("+init=EPSG:4326"))

然后对于投影，您需要 spTransform - 尝试：

cases.points.utm32 <- spTransform(cases.points, CRS(proj4string(region)))

使用不合适或不匹配的投影会产生很多问题，而且它们并不总是立即引人注意。

编辑：要在多边形内选择点，您需要 over() 函数，它也来自 sp 包（这是一个单独的问题）。请仔细阅读 sp 包的基本功能以及 sp 类的工作原理 - 下面基本上是从 over() 的帮助部分复制的。

region$pointCount <- sapply(over(region, geometry(cases.points), 
                             returnList = TRUE), length)

映射点和多边形

Mapping points and polygons

r

spatial

geospatial