如何确定一个点是否在缓冲区内
How to determine if a point is within a buffer
我正在处理一个如下所示的 SpatialPointsDataFrame：
new("SpatialPointsDataFrame", data = structure(list(date = c("01dec2013",
"01jul2003", "01nov2008", "01dec2017", "01dec2017", "01dec2003"
), company = c("Shwe Taung", "PetroChina Exploration and Development",
"Repsol SA", "Repsol SA", "Ipsen Pharmaceutical", "Ceva Laval"
), parent_company = c("Shwe Taung", "China National Petroleum (CNPC)",
"Repsol SA", "Repsol SA", "Ipsen Pharmaceutical", "Ceva Sante Animale"
), Website = c("www.shwetaunggroup.com", "www.cnpc.com.cn", "www.repsol.com",
"www.repsol.com", "www.ipsen.com", "www.ceva.com"), revenues_usd_ml = c(NA,
394554.53, 53215.45, 53215.45, 1760.671, 967.152), Headcount = c(NA,
1396144L, 24634L, 24634L, NA, 3500L), r_d_exp = c(NA, NA, 77.67,
77.67, NA, NA), est_year = c(NA, 1988L, 1927L, 1927L, 1929L,
1989L), o_country = c("Myanmar", "China", "Spain", "Spain", "France",
"France"), o_state = c("Rangoon (Yangon)", "Beijing Municipality",
"Comunidad de Madrid", "Comunidad de Madrid", "Ile-de-France",
"Sud-Ouest (FR)"), o_admin = c("Not Specified", "Not Specified",
"Madrid", "Madrid", "Ile-de-France", "Not Specified"), o_city = c("Rangoon (Yangon)",
"Beijing", "Madrid", "Madrid", "Paris", "Not Specified"), country = c("Algeria",
"Algeria", "Algeria", "Algeria", "Algeria", "Algeria"), state = c("Adrar",
"Adrar", "Adrar", "Adrar", "Adrar", "Adrar"), region = c("Not Specified",
"Not Specified", "Not Specified", "Not Specified", "Not Specified",
"Not Specified"), city = c("Adrar", "Adrar", "Reggane", "Reggane",
"Sidi Abdallah", "Sidi Abdallah"), free_zone = c("", "", "",
"", "", ""), relocation = c("", "", "", "", "", ""), sector = c("Building materials",
"Coal, oil & gas", "Coal, oil & gas", "Coal, oil & gas", "Pharmaceuticals",
"Healthcare"), sub_sector = c("Cement & concrete products", "Oil & gas extraction",
"Oil & gas extraction", "Oil & gas extraction", "Pharmaceutical preparations",
"Other (Healthcare)"), cluster = c("Construction", "Energy",
"Energy", "Energy", "Life sciences", "Life sciences"), activity = c("Manufacturing",
"Extraction", "Extraction", "Extraction", "Manufacturing", "Manufacturing"
), fdi_jobs = c(351L, 145L, 235L, 227L, 150L, 45L), est_fdi_jobs = c("Yes",
"Yes", "Yes", "Yes", "No", "No"), capital = c(139.9, 350, 565,
299.7, 29.55, 2.5), est_capital = c("Yes", "No", "No", "Yes",
"No", "No"), fdi_type = c("New", "New", "New", "Expansion", "New",
"New"), fdi_status = c("Announced", "Announced", "Announced",
"Opened", "Announced", "Opened"), g_lon = c(-0.287488, -0.287488,
0.1751507, 0.1751507, 5.0152099, 5.0152099), year = c(2013L,
2003L, 2008L, 2017L, 2017L, 2003L), code_d = c(12L, 12L, 12L,
12L, 12L, 12L), income_d = c("MIDLW", "MIDLW", "MIDLW", "MIDLW",
"MIDLW", "MIDLW"), continent_d = c("Africa", "Africa", "Africa",
"Africa", "Africa", "Africa"), lang_d = c("Arabic", "Arabic",
"Arabic", "Arabic", "Arabic", "Arabic"), landlocked = c(0L, 0L,
0L, 0L, 0L, 0L), iso_d = c("DZA", "DZA", "DZA", "DZA", "DZA",
"DZA"), isic = c(26L, 11L, 11L, 11L, 24L, 85L), isic4 = c(2695L,
1110L, 1110L, 1110L, 2411L, 8519L), sector_eora = c("Petroleum, Chemical and Non-Metallic Mineral Products",
"Mining and Quarrying", "Mining and Quarrying", "Mining and Quarrying",
"Petroleum, Chemical and Non-Metallic Mineral Products", "Mining and Quarrying"
), oc_lat = c(26.4888155, 26.4888155, 26.7207267, 26.7207267,
32.5966492, 32.5966492), oc_lng = c(-1.3582442, -1.3582442, 0.1751507,
0.1751507, -7.8308492, -7.8308492), oc_formatted = c("Adrar, Algeria",
"Adrar, Algeria", "Reggane, Reggane District, Algeria", "Reggane, Reggane District, Algeria",
"Sidi Abdallah, cercle de Rehamna, Morocco", "Sidi Abdallah, cercle de Rehamna, Morocco"
), oc_lat_st = c(26.4888155, 26.4888155, 26.4888155, 26.4888155,
26.4888155, 26.4888155), oc_lng_st = c(-1.3582442, -1.3582442,
-1.3582442, -1.3582442, -1.3582442, -1.3582442), oc_formatted_st = c("Adrar, Algeria",
"Adrar, Algeria", "Adrar, Algeria", "Adrar, Algeria", "Adrar, Algeria",
"Adrar, Algeria")), row.names = c(NA, 6L), class = "data.frame"),
coords.nrs = numeric(0), coords = structure(c(-1.3582442,
-1.3582442, 0.1751507, 0.1751507, -7.8308492, -7.8308492,
26.4888155, 26.4888155, 26.7207267, 26.7207267, 32.5966492,
32.5966492), .Dim = c(6L, 2L), .Dimnames = list(NULL, c("coords.x1",
"coords.x2"))), bbox = structure(c(-7.8308492, 26.4888155,
0.1751507, 32.5966492), .Dim = c(2L, 2L), .Dimnames = list(
c("coords.x1", "coords.x2"), c("min", "max"))), proj4string = new("CRS",
projargs = "+proj=longlat +datum=WGS84 +no_defs"))
我的任务是确定用 oc_lat_st 和 oc_lng_st 列标识的质心是否落在用 oc_lat 和 [=15= 列标识的质心周围 50 公里缓冲区内].
我目前生成的代码是：
library(rgdal)
library(ggplot2)
library(dplyr)
library(tidyr)
library(maptools)
library(rgeos)
FDI <- read.csv("FDI_harmonized.csv")
coords <- FDI[,40:41]
plot(coords)
coords <- drop_na(coords)
FDI <- FDI[!is.na(FDI$oc_lat), ]
coordinates(FDI) <- cbind(coords$oc_lng,coords$oc_lat)
#change coordinate reference system
geo_proj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(FDI) <- CRS(geo_proj)
#build buffer
distInMeters <- 1000
FDI50km <- gBuffer(FDI, width=100*distInMeters, byid=TRUE )
# Add data, and write to shapefile
FDI50km <- SpatialPolygonsDataFrame(FDI50km, data=FDI50km@data )
writeOGR( FDI50km, "FDI50km", "FDI50km", driver="ESRI Shapefile" )
plot(FDI50km)
但输出只是一个没有任何 CRS 的大缓冲区。
你有什么建议吗？理想的输出是一个额外的列，其中包含一个 T/F 或 1/0 值，用于表示第一组坐标是否落在缓冲区内。
通过像这样读取平面数据来使用 sf 包：
library(sf)
geo_proj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
df_new <- st_as_sf(df, coords = c("oc_lng", "oc_lat"), crs = geo_proj)
或者，将您的数据从 sp 转换为 sf：
library(sf)
df_new <- sf::st_as_sf(yourspatialdataframe)
然后将缓冲区定义为单独的功能 class:
# Define polygons
df_buff <- st_buffer(df_new, dist = 1000)
并按照 ?sf::st_intersects() 中的示例与您的点和缓冲区相交：
# Example data
pts = st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5, 2.5)))
pol = st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0))))
# View data
plot(x = NA, xlim = c(0,3), ylim = c(0,3), axes = FALSE, ann = FALSE)
plot(pts, col = "blue", add= TRUE)
plot(pol, add= TRUE)
# Intersect
lst = st_intersects(pts, pol)
lst
#> Sparse geometry binary predicate list of length 3, where the predicate
#> was `intersects'
#> 1: 1
#> 2: 1
#> 3: (empty)
mat = st_intersects(pts, pol, sparse = FALSE)
mat
#> [,1]
#> [1,] TRUE
#> [2,] TRUE
#> [3,] FALSE
# which points fall inside a polygon?
apply(mat, 1, any)
#> [1] TRUE TRUE FALSE
lengths(lst) > 0
#> [1] TRUE TRUE FALSE
# which points fall inside the first polygon?
st_intersects(pol, pts)[[1]]
#> [1] 1 2
我正在处理一个如下所示的 SpatialPointsDataFrame：
new("SpatialPointsDataFrame", data = structure(list(date = c("01dec2013",
"01jul2003", "01nov2008", "01dec2017", "01dec2017", "01dec2003"
), company = c("Shwe Taung", "PetroChina Exploration and Development",
"Repsol SA", "Repsol SA", "Ipsen Pharmaceutical", "Ceva Laval"
), parent_company = c("Shwe Taung", "China National Petroleum (CNPC)",
"Repsol SA", "Repsol SA", "Ipsen Pharmaceutical", "Ceva Sante Animale"
), Website = c("www.shwetaunggroup.com", "www.cnpc.com.cn", "www.repsol.com",
"www.repsol.com", "www.ipsen.com", "www.ceva.com"), revenues_usd_ml = c(NA,
394554.53, 53215.45, 53215.45, 1760.671, 967.152), Headcount = c(NA,
1396144L, 24634L, 24634L, NA, 3500L), r_d_exp = c(NA, NA, 77.67,
77.67, NA, NA), est_year = c(NA, 1988L, 1927L, 1927L, 1929L,
1989L), o_country = c("Myanmar", "China", "Spain", "Spain", "France",
"France"), o_state = c("Rangoon (Yangon)", "Beijing Municipality",
"Comunidad de Madrid", "Comunidad de Madrid", "Ile-de-France",
"Sud-Ouest (FR)"), o_admin = c("Not Specified", "Not Specified",
"Madrid", "Madrid", "Ile-de-France", "Not Specified"), o_city = c("Rangoon (Yangon)",
"Beijing", "Madrid", "Madrid", "Paris", "Not Specified"), country = c("Algeria",
"Algeria", "Algeria", "Algeria", "Algeria", "Algeria"), state = c("Adrar",
"Adrar", "Adrar", "Adrar", "Adrar", "Adrar"), region = c("Not Specified",
"Not Specified", "Not Specified", "Not Specified", "Not Specified",
"Not Specified"), city = c("Adrar", "Adrar", "Reggane", "Reggane",
"Sidi Abdallah", "Sidi Abdallah"), free_zone = c("", "", "",
"", "", ""), relocation = c("", "", "", "", "", ""), sector = c("Building materials",
"Coal, oil & gas", "Coal, oil & gas", "Coal, oil & gas", "Pharmaceuticals",
"Healthcare"), sub_sector = c("Cement & concrete products", "Oil & gas extraction",
"Oil & gas extraction", "Oil & gas extraction", "Pharmaceutical preparations",
"Other (Healthcare)"), cluster = c("Construction", "Energy",
"Energy", "Energy", "Life sciences", "Life sciences"), activity = c("Manufacturing",
"Extraction", "Extraction", "Extraction", "Manufacturing", "Manufacturing"
), fdi_jobs = c(351L, 145L, 235L, 227L, 150L, 45L), est_fdi_jobs = c("Yes",
"Yes", "Yes", "Yes", "No", "No"), capital = c(139.9, 350, 565,
299.7, 29.55, 2.5), est_capital = c("Yes", "No", "No", "Yes",
"No", "No"), fdi_type = c("New", "New", "New", "Expansion", "New",
"New"), fdi_status = c("Announced", "Announced", "Announced",
"Opened", "Announced", "Opened"), g_lon = c(-0.287488, -0.287488,
0.1751507, 0.1751507, 5.0152099, 5.0152099), year = c(2013L,
2003L, 2008L, 2017L, 2017L, 2003L), code_d = c(12L, 12L, 12L,
12L, 12L, 12L), income_d = c("MIDLW", "MIDLW", "MIDLW", "MIDLW",
"MIDLW", "MIDLW"), continent_d = c("Africa", "Africa", "Africa",
"Africa", "Africa", "Africa"), lang_d = c("Arabic", "Arabic",
"Arabic", "Arabic", "Arabic", "Arabic"), landlocked = c(0L, 0L,
0L, 0L, 0L, 0L), iso_d = c("DZA", "DZA", "DZA", "DZA", "DZA",
"DZA"), isic = c(26L, 11L, 11L, 11L, 24L, 85L), isic4 = c(2695L,
1110L, 1110L, 1110L, 2411L, 8519L), sector_eora = c("Petroleum, Chemical and Non-Metallic Mineral Products",
"Mining and Quarrying", "Mining and Quarrying", "Mining and Quarrying",
"Petroleum, Chemical and Non-Metallic Mineral Products", "Mining and Quarrying"
), oc_lat = c(26.4888155, 26.4888155, 26.7207267, 26.7207267,
32.5966492, 32.5966492), oc_lng = c(-1.3582442, -1.3582442, 0.1751507,
0.1751507, -7.8308492, -7.8308492), oc_formatted = c("Adrar, Algeria",
"Adrar, Algeria", "Reggane, Reggane District, Algeria", "Reggane, Reggane District, Algeria",
"Sidi Abdallah, cercle de Rehamna, Morocco", "Sidi Abdallah, cercle de Rehamna, Morocco"
), oc_lat_st = c(26.4888155, 26.4888155, 26.4888155, 26.4888155,
26.4888155, 26.4888155), oc_lng_st = c(-1.3582442, -1.3582442,
-1.3582442, -1.3582442, -1.3582442, -1.3582442), oc_formatted_st = c("Adrar, Algeria",
"Adrar, Algeria", "Adrar, Algeria", "Adrar, Algeria", "Adrar, Algeria",
"Adrar, Algeria")), row.names = c(NA, 6L), class = "data.frame"),
coords.nrs = numeric(0), coords = structure(c(-1.3582442,
-1.3582442, 0.1751507, 0.1751507, -7.8308492, -7.8308492,
26.4888155, 26.4888155, 26.7207267, 26.7207267, 32.5966492,
32.5966492), .Dim = c(6L, 2L), .Dimnames = list(NULL, c("coords.x1",
"coords.x2"))), bbox = structure(c(-7.8308492, 26.4888155,
0.1751507, 32.5966492), .Dim = c(2L, 2L), .Dimnames = list(
c("coords.x1", "coords.x2"), c("min", "max"))), proj4string = new("CRS",
projargs = "+proj=longlat +datum=WGS84 +no_defs"))
我的任务是确定用 oc_lat_st 和 oc_lng_st 列标识的质心是否落在用 oc_lat 和 [=15= 列标识的质心周围 50 公里缓冲区内].
我目前生成的代码是：
library(rgdal)
library(ggplot2)
library(dplyr)
library(tidyr)
library(maptools)
library(rgeos)
FDI <- read.csv("FDI_harmonized.csv")
coords <- FDI[,40:41]
plot(coords)
coords <- drop_na(coords)
FDI <- FDI[!is.na(FDI$oc_lat), ]
coordinates(FDI) <- cbind(coords$oc_lng,coords$oc_lat)
#change coordinate reference system
geo_proj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(FDI) <- CRS(geo_proj)
#build buffer
distInMeters <- 1000
FDI50km <- gBuffer(FDI, width=100*distInMeters, byid=TRUE )
# Add data, and write to shapefile
FDI50km <- SpatialPolygonsDataFrame(FDI50km, data=FDI50km@data )
writeOGR( FDI50km, "FDI50km", "FDI50km", driver="ESRI Shapefile" )
plot(FDI50km)
但输出只是一个没有任何 CRS 的大缓冲区。
你有什么建议吗？理想的输出是一个额外的列，其中包含一个 T/F 或 1/0 值，用于表示第一组坐标是否落在缓冲区内。
通过像这样读取平面数据来使用 sf 包：
library(sf)
geo_proj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
df_new <- st_as_sf(df, coords = c("oc_lng", "oc_lat"), crs = geo_proj)
或者，将您的数据从 sp 转换为 sf：
library(sf)
df_new <- sf::st_as_sf(yourspatialdataframe)
然后将缓冲区定义为单独的功能 class:
# Define polygons
df_buff <- st_buffer(df_new, dist = 1000)
并按照 ?sf::st_intersects() 中的示例与您的点和缓冲区相交：
# Example data
pts = st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5, 2.5)))
pol = st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0))))
# View data
plot(x = NA, xlim = c(0,3), ylim = c(0,3), axes = FALSE, ann = FALSE)
plot(pts, col = "blue", add= TRUE)
plot(pol, add= TRUE)
# Intersect
lst = st_intersects(pts, pol)
lst
#> Sparse geometry binary predicate list of length 3, where the predicate
#> was `intersects'
#> 1: 1
#> 2: 1
#> 3: (empty)
mat = st_intersects(pts, pol, sparse = FALSE)
mat
#> [,1]
#> [1,] TRUE
#> [2,] TRUE
#> [3,] FALSE
# which points fall inside a polygon?
apply(mat, 1, any)
#> [1] TRUE TRUE FALSE
lengths(lst) > 0
#> [1] TRUE TRUE FALSE
# which points fall inside the first polygon?
st_intersects(pol, pts)[[1]]
#> [1] 1 2