识别属于纬度和经度坐标的邮政编码 [英] identify zip codes that fall within latitude and longitudinal coordinates

问题描述

我在R中有几个数据帧.第一个数据帧包含按市场计算的一组经度和纬度坐标的凸包(由R中的chull提供).看起来像这样:

I have several data frames in R. The first data frame contains the computed convex hull of a set of lat and long coordinates by market (courtesy of chull in R). It looks like this:

MyGeo<- "Part of Chicago & Wisconsin"
Longitude <- c(-90.31914,  -90.61911,  -89.37842,  -88.0988,  -87.44875)
Latitude <- c(38.45781, 38.80097, 43.07961, 43.0624,41.49182)

dat <- data.frame(Longitude, Latitude, MyGeo)

第二个按其纬度和纵向坐标具有邮政编码(由美国人口普查网站提供).看起来像这样:

The second has zip codes by their latitude and longitudinal coordinates (courtesy of the US census website). It looks like this:

CensuseZip <- c("SomeZipCode1","SomeZipCode2","SomeZipCode3","SomeZipCode4","SomeZipCode5","SomeZipCode6","SomeZipCode7") 
Longitude2 <- c(-131.470425,-133.457924,-131.693453,-87.64957,-87.99734,-87.895,-88.0228)
Latitude2 <- c(55.138352,56.239062,56.370538,41.87485,42.0086,42.04957,41.81055)

cen <- data.frame(Longitude2, Latitude2,   CensuseZip)

现在，我相信第一个数据表为我提供了一个多边形或边界，我应该可以使用该边界或边界来识别该边界内的邮政编码.理想情况下，我想创建一个看起来像这样的第三个数据表:

Now I believe the first data table provides me with a polygon, or a border, that I should be able to use to identify zip codes that fall within that border. Ideally, I would want to create a third data table that looks something like this:

 Longitude2 Latitude2    CensusZip                        MyGeo
-131.470425 55.138352 SomeZipCode1  
-133.457924 56.239062 SomeZipCode2  
-131.693453 56.370538 SomeZipCode3
-87.64957    41.87485 SomeZipCode4  Part of Chicago & Wisconsin 
-87.99734     42.0086 SomeZipCode5  Part of Chicago & Wisconsin 
-87.895      42.04957 SomeZipCode6  Part of Chicago & Wisconsin 
-88.0228     41.81055 SomeZipCode7  Part of Chicago & Wisconsin 

从本质上讲，我正在寻找所有位于蓝色长点和纬度点之间的邮政编码(请参见下面的可点击图像).虽然在下面可以看到，但实际上我正在寻找上述表格.

In essence, I am looking to identify all the zip codes that fall between the blue (see clickable image below) long and lat points. While it is visualized below, I am actually looking for the table described above.

但是...我在执行此操作时遇到了麻烦...我尝试使用以下软件包和脚本:

However... I am having trouble doing this... I have tried using the below packages and script:

library(rgeos)
library(sp)
library(rgdal)

coordinates(dat) <- ~ Longitude + Latitude
coordinates(cen) <- ~ Longitude2 + Latitude2

over(cen, dat)

但是我收到所有NA.

推荐答案

我使用library(sf)来解决这种类型的多边形点问题(sf是sp的后继者).

I use library(sf) to solve this type of point-in-polygon problem (sf is the successor to sp).

函数sf::st_intersection()为您提供两个sf对象的交集.根据您的情况，您可以构造单独的POLYGON和POINT sf对象.

The function sf::st_intersection() gives you the intersection of two sf objects. In your case you can construct separate POLYGON and POINT sf objects.

library(sf)

Longitude <- c(-90.31914,  -90.61911,  -89.37842,  -88.0988,  -87.44875)
Latitude <- c(38.45781, 38.80097, 43.07961, 43.0624,41.49182)

## closing the polygon
Longitude[length(Longitude) + 1] <- Longitude[1]
Latitude[length(Latitude) + 1] <- Latitude[1]

## construct sf POLYGON
sf_poly <- sf::st_sf( geometry = sf::st_sfc( sf::st_polygon( x = list(matrix(c(Longitude, Latitude), ncol = 2)))) )

## construct sf POINT
sf_points <- sf::st_as_sf( cen, coords = c("Longitude2", "Latitude2"))

sf::st_intersection(sf_points, sf_poly)

# Simple feature collection with 4 features and 1 field
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: -88.0228 ymin: 41.81055 xmax: -87.64957 ymax: 42.04957
# epsg (SRID):    NA
# proj4string:    NA
# CensuseZip                   geometry
# 4 SomeZipCode4 POINT (-87.64957 41.87485)
# 5 SomeZipCode5  POINT (-87.99734 42.0086)
# 6 SomeZipCode6   POINT (-87.895 42.04957)
# 7 SomeZipCode7  POINT (-88.0228 41.81055)
# Warning message:
#   attribute variables are assumed to be spatially constant throughout all geometries 

结果是多边形内的所有点

The result is all the points which are inside the polygon

您也可以使用sf::st_join(sf_poly, sf_points)给出相同的结果

You can also use sf::st_join(sf_poly, sf_points) to give the same result

然后，函数sf::st_intersects(sf_points, sf_poly)将返回一个列表，说明给定的POINT是否在多边形内

And, the function sf::st_intersects(sf_points, sf_poly) will return a list saying whether the given POINT is inside the polygon

sf::st_intersects(sf_points, sf_poly)

# Sparse geometry binary predicate list of length 7, where the predicate was `intersects'
#  1: (empty)
# 2: (empty)
# 3: (empty)
# 4: 1
# 5: 1
# 6: 1
# 7: 1

您可以将其用作原始sf_points对象的索引/标识符，以在其上添加新列

Which you can use as an index / identifier of the original sf_points object to add a new column on

is_in <- sf::st_intersects(sf_points, sf_poly)

sf_points$inside_polygon <- as.logical(is_in)

sf_points
# Simple feature collection with 7 features and 2 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: -133.4579 ymin: 41.81055 xmax: -87.64957 ymax: 56.37054
# epsg (SRID):    NA
# proj4string:    NA
# CensuseZip                   geometry inside_polygon
# 1 SomeZipCode1 POINT (-131.4704 55.13835)             NA
# 2 SomeZipCode2 POINT (-133.4579 56.23906)             NA
# 3 SomeZipCode3 POINT (-131.6935 56.37054)             NA
# 4 SomeZipCode4 POINT (-87.64957 41.87485)           TRUE
# 5 SomeZipCode5  POINT (-87.99734 42.0086)           TRUE
# 6 SomeZipCode6   POINT (-87.895 42.04957)           TRUE
# 7 SomeZipCode7  POINT (-88.0228 41.81055)           TRUE

这篇关于识别属于纬度和经度坐标的邮政编码的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！