【发布时间】：2020-06-19 11:44:23
【问题描述】：
我有数据可以映射每天的风暴纬度，如下所示（我删除了日期 ID，因此它只有年和月）
lat_lonDat <- structure(list(Latitude = c(9.1, 9.15755, 9.2, 9.21496, 9.06,
8.89169, 8.84286, 9.04619, 9.3, 9.31138, 9.3, 9.45992, 9.65,
9.76127, 9.85, 9.94991, 10.05, 10.1651, 10.25, 10.2725, 10.25,
10.1887, 10.15, 10.21, 10.3, 10.3163, 10.4, 10.6962, 11, 11.0975,
11.15, 11.2587, 11.45, 11.755, 12.15, 12.615, 13.05, 13.3338,
13.55, 13.775, 14, 14.2463, 14.5, 14.765, 15, 15.1537, 15.3,
15.525, 15.75, 15.8963, 16.05, 16.2963, 16.55, 16.7062, 16.9571,
17.4192, 18.04, 18.735, 19.2, 19.7925, 20.4, 20.9275, 21.4),
Longitude = c(88.1, 87.885, 87.7, 87.5575, 87.32, 87.1017,
86.9429, 86.8625, 86.7, 86.3436, 85.95, 85.685, 85.45, 85.16,
84.9, 84.7537, 84.6, 84.315, 84, 83.7237, 83.5, 83.36, 83.25,
83.1113, 82.95, 82.7462, 82.55, 82.4213, 82.3, 82.1137, 81.95,
81.8737, 81.85, 81.8462, 81.85, 81.8475, 81.8, 81.6687, 81.5,
81.3462, 81.2, 81.0488, 80.95, 80.9599, 81, 80.9788, 80.95,
80.9724, 80.95, 80.8087, 80.6, 80.3788, 80.15, 79.9423, 79.7857,
79.6735, 79.64, 79.7275, 79.7, 79.72, 79.8, 79.9348, 80.1),
yearRef = c(1990, 1990, 1990, 1990, 1990, 1990, 1990,
1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990,
1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990,
1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990,
1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990,
1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990, 1990,
1990, 1990, 1990, 1990, 1990, 1990), monthRef = c(5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5)), class = "data.frame", row.names = c(NA, -63L))
将其绘制在 shapefile 上
library(raster)
library(sp)
library(ggplot2)
temp.shp <- getData('GADM', country = 'IND', level = 2)
ggplot() +
geom_polygon(data = temp.shp, aes(x = long, y = lat, group = group), fill = NA, colour = "black") +
geom_point(data = lat_lonDat, aes(x = Longitude, y = Latitude, colour = 'red')) +
coord_quickmap() + theme_classic()
从风暴数据中，我想提取风暴在陆地上（即与多边形相交）时的那些行（纬度经度）。我从这里尝试了这种方法：
检查点是否在多边形 Shapefile 中
coordinates(lat_lonDat) <- ~ Latitude + Longitude
proj4string(lat_lonDat) <- proj4string(temp.shp)
over(lat_lonDat, temp.shp)
但是，我仍然收到 NA。我做错了什么？
【问题讨论】：
标签： r intersection sp rgdal