Я получаю странные результаты, когда пересекаю некоторые точечные данные с многоугольниками с помощью пакета 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 не сохраняет однозначно порядок осей. » Но это, по сути, мои предположения, и я не уверен, почему мое решение сработает в данном случае.
По сути, если кто-нибудь может сказать мне, что происходит или упускаю ли я какие-либо другие потенциальные проблемы, я буду признателен.
При работе с пространственными данными не должно быть «случайного» выбора. Вы должны знать, с чего начать, чтобы получить правильный результат. Вы не можете угадать, и это особенно верно в отношении данных по Северной Америке. Это потому, что, например, координаты NAD83/EPSG:4269 и WGS84/EPSG:4326 выглядят очень похоже. Другая распространенная ошибка — применение координат PCS к координатам GCS, что может привести к тому, что данные будут отображаться на расстоянии миллионов метров от их правильного местоположения. Подробности смотрите в моих ответах здесь и здесь.
Скорее всего проблема в том, что вы поменяли местами широту и долготу.
Я по очереди рассмотрю ваши вопросы, чтобы прояснить этот вопрос. Сначала загрузите необходимые пакеты и создайте некоторые репрезентативные данные для использования:
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 легла в основу ваших фактических данных, прежде чем приступать к прогнозированию.
Ваш второй пример (или
st_as_sf(..., crs = "WGS84")
) — это именно то, чему вам следует следовать.WGS84
не случайна; если вы работаете с географическими координатами, скорее всего, это CRS, который вам следует установить при создании объектаsf
из CSV или какого-либо другого источника, который не содержит CRS (если метаданные набора данных не указывают что-то другое). Вы получите значимые результаты с помощьюst_crs(metro) <- st_crs(tract00)
только в том случае, если координатыmetro
действительно находятся в той же системе, что иtract00
— координаты в CSV также могут использоватьESRI:102003
или любую другую CRS.