Блог

Главная - Вопросы по программированию - Как определить, находится ли точка в буфере

Как определить, находится ли точка в буфере

#r #buffer #coordinates #spatial #geosphere

Вопрос:

Я обрабатываю пространственно-временной кадр, который выглядит следующим образом:

 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 amp; gas", "Coal, oil amp; gas", "Coal, oil amp; gas", "Pharmaceuticals", 
"Healthcare"), sub_sector = c("Cement amp; concrete products", "Oil amp; gas extraction", 
"Oil amp; gas extraction", "Oil amp; 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 в пределах 50-километрового буфера вокруг центроида, идентифицированного со столбцами oc_lat и oc_lng .

Код, который я создал до сих пор, таков:

 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 для определения того, попадает ли первый набор координат в буфер.

Комментарии:

1. Таким образом, gBuffer() буферизует объект, чтобы сделать его больше. Вам нужно будет использовать gIntersects() для возврата логического вектора, геометрия которого пересекается, или gIntersection() для возврата фактических пересекающихся областей. В качестве альтернативы, если вы используете более новый sf пакет (который заменяет sp его ), вы используете st_intersects() и st_intersection для той же работы. Пояснения к sf функциям можно найти здесь geocompr.robinlovelace.net/geometric-operations.html , особенно Раздел 5.2

2. Хорошо, так help(st_intersects) говорит мне, что x и y должно быть классно sf или sfc . У меня есть sf объект с координатами буфера и пространственный кадр Pointsdataframe с координатами точек. Должен ли я использовать первый как x, а второй как y и наоборот?

3. st_* Функции взяты из sf пакета и используются с sf объектами, в то время g* как функции взяты из rgeos пакета и работают для sp объектов. Таким образом, вы должны либо использовать gIntersects() свои текущие sp объекты, либо преобразовать свой рабочий процесс в sf (рекомендуется) и использовать st_intersects() .

Ответ №1:

Используйте 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)

Затем определите свой буфер как отдельный класс объектов:

 # 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
Метки: ,