【问题标题】:How to filter a dataframe of geocoordinates using a KML polygon?如何使用 KML 多边形过滤坐标数据框?
【发布时间】:2015-08-17 00:35:15
【问题描述】:

我有一个涵盖整个地球 (the US Geological Service's earthquake feed) 的数据点的 CSV,我想只过滤美国境内的地震。

我从美国人口普查局提取的 KML 文件:

https://www.census.gov/geo/maps-data/data/kml/kml_nation.html

在 R 中,rgdal 库可以加载 KML 文件:

library(rgdal)
kml = readOGR("kmls/cb_2014_us_nation_20m.kml", 'cb_2014_us_nation_20m')

如何使用dplyr/plyr/等。仅针对 KML 指定的边界内的行过滤 data.frame(地理数据的列是 latitudelongitude)?


编辑,后回答:

这是我从@hrbrmstr's answer 使用的用于快速可视化的内容:

library(sp)
library(rgdal)
# download earthquakes
url <- "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_week.csv"
fil <- "all_week.csv"
if (!file.exists(fil)) download.file(url, fil)
quakes <- read.csv("all_week.csv", stringsAsFactors=FALSE)
# create spatial object
sp::coordinates(quakes) <- ~longitude+latitude

# download nation KML
url <- "http://www2.census.gov/geo/tiger/GENZ2014/kml/cb_2014_us_nation_20m.zip"
fil <- "uskml.zip"
if (!file.exists(fil)) download.file(url, fil)
unzip(fil, exdir="uskml")
# read KML file
us <- rgdal::readOGR("./uskml/cb_2014_us_nation_20m.kml", "cb_2014_us_nation_20m")
sp::proj4string(quakes) <- sp::proj4string(us)

length(quakes)
# 1514

usquakes = quakes[us,]
length(usquakes)
# 1260

### visualize
plot(us) 
# plot all quakes
points(quakes$longitude, quakes$latitude)
# plot just US
points(usquakes$longitude, usquakes$latitude)

生成的图像:

谢谢@hrbrmstr!

【问题讨论】:

  • @hrbrmstr:我会满足于任何一个,但我猜你会问,因为与非连续美国,多边形数据结构不同?假设我只关心大陆,尽管我现在使用的文件是人口普查的国家边界:census.gov/geo/maps-data/data/kml/kml_nation.html

标签: r kml geospatial


【解决方案1】:

这将是一个 cpl 方式:

library(sp)
library(maptools)
library(rgeos) # not entirely necessary
library(rgdal) # not entirely necessary

url <- "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_week.csv"
fil <- "all_week.csv"
if (!file.exists(fil)) download.file(url, fil)

quakes <- read.csv("all_week.csv", stringsAsFactors=FALSE)
coordinates(quakes) <- ~longitude+latitude

url <- "http://www2.census.gov/geo/tiger/GENZ2014/shp/cb_2014_us_nation_20m.zip"
fil <- "us.zip"
if (!file.exists(fil)) download.file(url, fil)
unzip(fil, exdir="us")
us <- readShapePoly("us/cb_2014_us_nation_20m.shp")

# alternatively
# us <- rgdal::readOGR("us/cbcb_2014_us_nation_20m.shp", "cb_2014_us_nation_20m")

# TRUE if in
in_us_rgeos <- rgeos::gIntersects(quakes, us, byid=TRUE)

# <NA> if in
in_us_over <- quakes %over% us

gIntersects 需要更长的时间。 rgdalrgeos 可以在某些系统上工作。 YMMV

使用美国 KML 将需要(大部分)rgdal:

# you wanted KML shapefile tho
url <- "http://www2.census.gov/geo/tiger/GENZ2014/kml/cb_2014_us_nation_20m.zip"
fil <- "uskml.zip"
if (!file.exists(fil)) download.file(url, fil)
unzip(fil, exdir="uskml")

us <- rgdal::readOGR("uskml/cb_2014_us_nation_20m.kml", "cb_2014_us_nation_20m")
proj4string(quakes) <- proj4string(us)
rgeos::gIntersects(quakes, us, byid=TRUE)
head(quakes %over% us)

【讨论】:

  • 谢谢!...快速提问...为什么over 看起来比gIntersects 快那么多?
  • rgeos 是 Java 代码的 C++ 端口(这可能意味着它承载了一些 Java 代码包袱)但gIntersects更灵活(它可以处理所有Spatial… 类型)。 over(对于 SpatialPointsSpatialPolygons 中的情况)针对该特殊情况调用 super 高效 C 代码(函数为 R_point_in_polygon_sp)。
猜你喜欢
  • 2016-09-11
  • 2014-08-02
  • 2019-04-26
  • 1970-01-01
  • 2012-10-29
  • 1970-01-01
  • 2012-07-05
  • 2015-04-10
  • 2021-05-29
相关资源
最近更新 更多