检查点是否与多边形重叠

Question

我有如下数据，每天映射风暴的经纬度（我删除了日期 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。我做错了什么？

Answer 1

long/lat 坐标应交换：

  coordinates(lat_lonDat) <- ~ Longitude + Latitude

结果：

R> tail(over(lat_lonDat, temp.shp), 15)
   GID_0 NAME_0    GID_1         NAME_1 NL_NAME_1       GID_2     NAME_2
49   IND  India  IND.2_1 Andhra Pradesh      <NA>   IND.2.5_1    Krishna
50   IND  India  IND.2_1 Andhra Pradesh      <NA>   IND.2.4_1     Guntur
51   IND  India  IND.2_1 Andhra Pradesh      <NA>   IND.2.4_1     Guntur
52   IND  India  IND.2_1 Andhra Pradesh      <NA>   IND.2.4_1     Guntur
53   IND  India  IND.2_1 Andhra Pradesh      <NA>   IND.2.4_1     Guntur
54   IND  India IND.32_1      Telangana      <NA>  IND.32.7_1   Nalgonda
55   IND  India IND.32_1      Telangana      <NA>  IND.32.7_1   Nalgonda
56   IND  India IND.32_1      Telangana      <NA>  IND.32.7_1   Nalgonda
57   IND  India IND.32_1      Telangana      <NA> IND.32.10_1   Warangal
58   IND  India IND.32_1      Telangana      <NA>  IND.32.1_1   Adilabad
59   IND  India IND.32_1      Telangana      <NA>  IND.32.1_1   Adilabad
60   IND  India IND.20_1    Maharashtra      <NA>  IND.20.8_1 Chandrapur
61   IND  India IND.20_1    Maharashtra      <NA>  IND.20.8_1 Chandrapur
62   IND  India IND.20_1    Maharashtra      <NA>  IND.20.5_1   Bhandara
63   IND  India IND.20_1    Maharashtra      <NA> IND.20.11_1    Gondiya
   VARNAME_2 NL_NAME_2   TYPE_2 ENGTYPE_2 CC_2   HASC_2
49    Kistna      <NA> District  District <NA> IN.AD.KR
50      <NA>      <NA> District  District <NA> IN.AD.GU
51      <NA>      <NA> District  District <NA> IN.AD.GU
52      <NA>      <NA> District  District <NA> IN.AD.GU
53      <NA>      <NA> District  District <NA> IN.AD.GU
54      <NA>      <NA> District  District <NA> IN.TG.NA
55      <NA>      <NA> District  District <NA> IN.TG.NA
56      <NA>      <NA> District  District <NA> IN.TG.NA
57      <NA>      <NA> District  District <NA> IN.TG.WA
58      <NA>      <NA> District  District <NA> IN.TG.AD
59      <NA>      <NA> District  District <NA> IN.TG.AD
60    Chanda      <NA> District  District <NA> IN.MH.CH
61    Chanda      <NA> District  District <NA> IN.MH.CH
62      <NA>      <NA> District  District <NA> IN.MH.BH
63      <NA>      <NA> District  District <NA> IN.MH.GO

检查点是否与多边形重叠

check if points overlap with polygon

intersection

r

rgdal

sp