映射点和多边形
我是一个R初学者第一次处理与R和空间数据。所以我希望我能说清楚。
我有一个意大利地区的档案，分为几个人口普查区。另一方面，我有一个csv文件，有一个案例列表和每个案例的地址。我想要地图点和形状文件，并得到有多少点是在每个人口普查部门。
以下是我从现在开始所做的事情：
#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，但不能将它们添加到相同的绘图中。同时，正如我所说，我想计算每个多边形(即“区域”的人口普查划分)中有多少个点。
我必须承认我对地理不太感兴趣。我怀疑不同的坐标和/或投影参考系统可能是问题所在，因此我进行了检查。
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"形状文件坐标是否可能是以度为单位，而大小写坐标是十进制的呢？如果是的话，如何将其转化？
谢谢你，萨罗
Stack Overflow用户
发布于 2017-06-27 13:14:49
您有两个不同的投影-您的“区域”形状文件投影在UTM区32和您告诉R使用Web Mercator的“案例”。然而，如果你只是从谷歌下载案例数据为lat/long，那么你就不想告诉R他们的投影是Web Mercator，因为它不是--它是未投影的WGS 84，所以你会想要EPSG 4326。what是谷歌用于显示地图的工具，但如果下载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()的help部分复制的。
region$pointCount <- sapply(over(region, geometry(cases.points),
returnList = TRUE), length)https://stackoverflow.com/questions/44780107
复制相似问题
腾讯云开发者
