Неправильные результаты при пересечении данных точек и полигонов с sf в R

Я получаю странные результаты, когда пересекаю некоторые точечные данные с многоугольниками с помощью пакета sf в R, и хотя я думаю, что нашел решение, я настолько озадачен этой проблемой, что был бы очень признателен, если бы кто-нибудь мог помогите мне понять, что происходит и, следовательно, упускаю ли я какие-либо потенциальные проблемы.

У меня есть один набор данных только с координатами; Мне нужно знать, в каком участке 2000 находится каждая точка, поэтому я хотел бы пересечь его с помощью шейп-файла участка 2000, чтобы назначить правильные идентификаторы участка. Мои реальные данные конфиденциальны, поэтому в качестве примера я сейчас использую карту метро Вашингтона, доступную здесь. Мои настоящие данные — это просто координаты, а не шейп-файл, поэтому изначально у них нет проекции. Я удалю проекцию со карты метро, ​​чтобы имитировать реальные данные (и продублирую данные, чтобы показать свое решение):

library(sf)

# open metro shapefile using file paths not shown here
metro <- read_sf(dsn = metro.fp, layer = metro.filename)

dim(metro)
[1] 98 15

# remove metro projection to mirror initial state of real data
st_crs(metro) <- ""

# duplicate this data for later example of my solution
metro2 <- metro

Я использую шейп-файл с полигонами участков переписи населения США 2000 года (в геометрии TIGER/Line 2010 года), загруженный с NHGIS.org. (Если есть способ проще поделиться этим примером, дайте мне знать; я мог бы использовать dput(), но даже в сочетании с head() результат огромен.) Здесь используется проекция ESRI, которая, я думаю, имеет отношение к проблеме. под рукой. Я не буду показывать полный вывод, когда смотрю на его проекцию, потому что он очень длинный, но когда вы используете st_crs(), он показывает, что на этой карте используется USA_Contiguous_Albers_Equal_Area_Conic с ID["ESRI",102003] внизу.

Когда я конвертирую свой набор точечных данных в ту же проекцию, что и многоугольники, используя st_crs() вместо st_transform(), потому что в моих точечных данных изначально нет проекции, а затем пересекаю их, я получаю странную проблему, когда все мои точки предположительно находятся в Канзасе (20 штатов), хотя все это станции метро в округе Колумбия, Мэриленде и Вирджинии:

# open tract map using file paths not shown here
tract00 <- read_sf(dsn = fp, layer = file.name)

# give metro same projection as tract map
st_crs(metro) <- st_crs(tract00)

# intersect the two
metro.tract <- st_intersection(x = metro,
                               y = tract00)

# now all states are kansas (state fips of 20):
table(metro.tract$STATEFP00)

20 
98

Если вместо этого я дам своему исходному набору данных метро случайную проекцию, а затем использую st_transform, чтобы поместить его в ту же проекцию, что и шейп-файл участка, то я правильно получу, что станции существуют в Вашингтоне, Мэриленде и Вирджинии (DC — 11, MD — 24 , а ВА – 51):

# now try again, first giving metro2 a random projection so I can use st_transform
st_crs(metro2) <- "WGS84" 

# give metro2 same projection as tract map
metro2 <- st_transform(metro2, crs = st_crs(tract00))

# intersect the two
metro.tract2 <- st_intersection(x = metro2,
                                y = tract00)

# now the states are correct:
table(metro.tract2$STATEFP00)

11 24 51 
40 26 32 

Кто-нибудь знает, что происходит? Хотя мое решение здесь дает результаты, которые выглядят правильно, поскольку я не понимаю, в чем именно проблема, я беспокоюсь, что могу упустить другие потенциальные проблемы с моими реальными данными.

Я думаю, что это как-то связано с проекцией ESRI шейп-файла тракта, поскольку в sf документации st_transform говорится: «С форматом шейп-файла ESRI требуется особая осторожность, поскольку WKT1 не сохраняет однозначно порядок осей. » Но это, по сути, мои предположения, и я не уверен, почему мое решение сработает в данном случае.

По сути, если кто-нибудь может сказать мне, что происходит или упускаю ли я какие-либо другие потенциальные проблемы, я буду признателен.

Ваш второй пример (или st_as_sf(..., crs = "WGS84")) — это именно то, чему вам следует следовать. WGS84 не случайна; если вы работаете с географическими координатами, скорее всего, это CRS, который вам следует установить при создании объекта sf из CSV или какого-либо другого источника, который не содержит CRS (если метаданные набора данных не указывают что-то другое). Вы получите значимые результаты с помощью st_crs(metro) <- st_crs(tract00) только в том случае, если координаты metro действительно находятся в той же системе, что и tract00 — координаты в CSV также могут использовать ESRI:102003 или любую другую CRS.

margusl 24.04.2024 18:37

При работе с пространственными данными не должно быть «случайного» выбора. Вы должны знать, с чего начать, чтобы получить правильный результат. Вы не можете угадать, и это особенно верно в отношении данных по Северной Америке. Это потому, что, например, координаты NAD83/EPSG:4269 и WGS84/EPSG:4326 выглядят очень похоже. Другая распространенная ошибка — применение координат PCS к координатам GCS, что может привести к тому, что данные будут отображаться на расстоянии миллионов метров от их правильного местоположения. Подробности смотрите в моих ответах здесь и здесь.

L Tyrone 24.04.2024 20:46

Скорее всего проблема в том, что вы поменяли местами широту и долготу.

Ian Turton 24.04.2024 20:54
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
1
3
104
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Я по очереди рассмотрю ваши вопросы, чтобы прояснить этот вопрос. Сначала загрузите необходимые пакеты и создайте некоторые репрезентативные данные для использования:

library(tigris) # For US census tracts
library(sf)
library(dplyr)

# Read tract data for DC. Project to ESRI:102003 as default CRS is NAD83:EPSG4269 
tract00 <- tracts(state = "dc", year = 2000) %>%
  st_transform("ESRI:102003")

tract00$geometry
# Geometry set for 188 features 
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: 1610777 ymin: 306994.4 xmax: 1629412 ymax: 329361
# Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
# First 5 geometries:
# MULTIPOLYGON (((1617395 319019.9, 1617392 31902...
# MULTIPOLYGON (((1617767 319815.2, 1617772 31981...
# MULTIPOLYGON (((1616807 318539.2, 1616794 31854...
# MULTIPOLYGON (((1617619 318298.1, 1617601 31829...
# MULTIPOLYGON (((1618010 318832.6, 1618021 31883...

# DC metro data shp previously unzipped to working directory, downloaded from 
# https://opendata.dc.gov/datasets/DCGIS::metro-stations-regional/explore
metro <- st_read("Metro_Stations_Regional.shp") 

metro$geometry
# Geometry set for 98 features 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -77.49154 ymin: 38.76653 xmax: -76.84455 ymax: 39.11994
# Geodetic CRS:  WGS 84
# First 5 geometries:
# POINT (-76.91147 38.82645)
# POINT (-77.05367 38.81415)
# POINT (-77.06081 38.80659)
# POINT (-77.07088 38.80043)
# POINT (-77.07521 38.79392)

Обратите внимание на разницу в значениях между Metro и Tract00.

Когда я конвертирую свой набор точечных данных в ту же проекцию, что и полигоны... у меня возникает странная проблема: все мои точки предположительно находятся в Канзасе...

Единицами координат для набора данных метро округа Колумбия, с которым вы связались, являются десятичные градусы, тогда как единицами координат для данных участка являются метры. Запустив:

st_crs(metro) <- st_crs(tract00)
# Warning message:
#   st_crs<- : replacing crs does not reproject data; use st_transform for that

вы определяете эти десятичные значения градусов как метры. Именно поэтому вы получите предупреждение, если сделаете это. В результате ваши данные метро будут расположены на расстоянии > 1,6e6 м от того места, где они должны быть по долготе, и на расстоянии > 3,06e5 м по широте, например. Дороти и Тото вернулись в Канзас ;)

Если вместо этого я дам своему исходному набору данных метро случайную проекцию, а затем использую st_transform, чтобы поместить его в ту же проекцию, что и шейп-файл участка, то я правильно пойму, что станции существуют в Вашингтоне...

Оказывается, в этом гипотетическом случае вам «повезло», присвоив WGS84:EPSG4326 данным Metro2. Это потому, что мы знаем, что исходный файл Metro имеет WGS84 в качестве CRS. Итак, когда вы побежали:

st_crs(metro2) <- "WGS84" 

metro2 <- st_transform(metro2, crs = st_crs(tract00))

функция st_transform() сказала: «Я вижу, что метро указано в десятичных градусах, и я знаю, где оно расположено, на основе CRS, поэтому я преобразую эти значения координат в десятичных градусах в метры и найду их, используя CRS тракта 00».

Однако, как я уже сказал, вам «повезло». В случае вашего фактического набора данных вы не указали, какая CRS составляет основу значений координат. Это проблема, поскольку значения географической системы координат (GCS) выглядят одинаково, но могут сильно различаться. Учтите следующее:

# Reload metro dataset
metro <- st_read("Metro_Stations_Regional.shp")

# Convert metro to df with 'unknown' coordinates to replicate your actual data
metro2 <- data.frame(id = 1:nrow(metro),
                     st_coordinates(metro))

head(metro2)
#   id         X        Y
# 1  1 -76.91147 38.82645
# 2  2 -77.05367 38.81415
# 3  3 -77.06081 38.80659
# 4  4 -77.07088 38.80043
# 5  5 -77.07521 38.79392
# 6  6 -76.99537 38.86297

# Coordinates incorrectly assumed to be NAD27/EPSG:4267
p_4267 <- st_as_sf(metro2, coords = c("X", "Y"), crs = 4326) %>%
  st_transform(4267) %>%
  st_set_crs(4326)
# Warning message:
#   st_crs<- : replacing crs does not reproject data; use st_transform for that

p_4267$geometry
# Geometry set for 98 features 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -77.49185 ymin: 38.76649 xmax: -76.84488 ymax: 39.11991
# Geodetic CRS:  WGS 84
# First 5 geometries:
# POINT (-76.91179 38.82642)
# POINT (-77.05399 38.81412)
# POINT (-77.06114 38.80656)
# POINT (-77.0712 38.8004)
# POINT (-77.07553 38.79389)

# Return distance between original and incorrectly projected location
head(st_distance(metro, p_4267 , by_element = TRUE))
# Units: [m]
# [1] 28.61364 28.23506 28.21754 28.19191 28.18171 28.38108

Несмотря на то, что метро имеет тот же CRS, что и p_4267, и значения координат выглядят одинаково, на самом деле они находятся на расстоянии ~ 28 метров друг от друга. Таким образом, вам необходимо знать, какая CRS легла в основу ваших фактических данных, прежде чем приступать к прогнозированию.

Другие вопросы по теме