Конвертировать XY в lat long в R
Я пытаюсь нанести на карту листовки карту портов, которые попадают в определенный регион. Эти 30 странных областей были созданы с использованием mapedit
пакет как sf
Многоугольники. Код для одного такого региона WestCoastIndiax
(Западное побережье Индии) выглядит следующим образом:
structure(list(X_leaflet_id = 629L, feature_type = "polygon",
geometry = structure(list(structure(list(structure(c(-298.1964,
-298.4601, -298.3722, -297.6855, -293.9062, -293.4009, -292.8296,
-292.1265, -291.7529, -289.3359, -288.1934, -289.248, -289.0723,
-287.7539, -284.8535, -286.04, -283.0957, -282.3157, -282.4255,
-282.9529, -283.4253, -284.0295, -284.9634, -285.8643, -286.7651,
-286.7432, -286.875, -286.9519, -288.2922, -288.6328, -290.05,
-291.8958, -292.8955, -293.4009, -294.6863, -298.1964, 25.6861,
24.8565, 24.682, 24.8765, 25.1453, 24.5671, 23.8256, 23.5237,
22.6343, 19.6012, 20.0559, 12.5975, 7.1881, -1.1425, -0.791,
8.4941, 8.0484, 7.8416, 8.2441, 9.2973, 10.4338, 11.62, 14.094,
16.1197, 19.3215, 20.0456, 21.8411, 22.6748, 23.433, 24.4572,
24.7768, 25.0657, 25.3043, 25.8592, 26.116, 25.6861), .Dim = c(36L,
2L))), class = c("XY", "POLYGON", "sfg"))), class = c("sfc_POLYGON",
"sfc"), precision = 0, bbox = structure(c(xmin = -298.4601,
ymin = -1.1425, xmax = -282.3157, ymax = 26.116), class = "bbox"), crs = structure(list(
epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L)), row.names = 1L, sf_column = "geometry", agr = structure(c(X_leaflet_id = NA_integer_,
feature_type = NA_integer_), class = "factor", .Label = c("constant",
"aggregate", "identity")), class = c("sf", "data.frame"))
Область правильно отображается в leaflet
,
leaflet() %>% addTiles() %>% addPolygons(data=WestCoastIndiax)
Затем я пытаюсь увидеть, попадает ли конкретное место в эти регионы, используя st_within()
функция от sf
пакет. Пример города Мумбаи, который попадает в этот регион:
st_within(st_point(c(72.49,18.54)),WestCoastIndiax)
который возвращает:
Sparse geometry binary predicate list of length 1, where the predicate was `within'
1: (empty)
Когда я пытаюсь построить полигон и координаты Мумбаи, я получаю следующее:
leaflet() %>% addTiles() %>% addPolygons(data=WestCoastIndiax) %>% addMarkers(lng = 72.49, lat = 18.54)
Я посмотрел на это и понял, что мне нужно правильно crs
и затем преобразовать его снова в 4326. Однако я не знаю, какой должна быть система координат для этих 30 нечетных геометрий. Я попробовал это:
st_crs(WestCoastIndiax) <- "+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
WestCoastIndiax <- st_transform(WestCoastIndiax, 4326)
И хотя сейчас bbox
единицы поменялись, но они не верны:
Geometry set for 1 feature
geometry type: POLYGON
dimension: XY
bbox: xmin: -0.002681113 ymin: -1.026325e-05 xmax: -0.002536085 ymax: 0.000234604
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
POLYGON ((-0.002678744 0.0002307422, -0.0026811...
Буду признателен за помощь здесь, пожалуйста.
1 ответ
Ваши данные в CSR 4326 хорошо! То, что произошло, очень вероятно, что при создании многоугольника с mapedit объекты были нарисованы не в том месте (то есть вне -180/180 долготы). Вы можете проверить это, вычтя 360° из вашего точечного местоположения, которое затем находится в пределах вашего многоугольника.
mapview::mapview(tst) + sf::st_point(c(72.49-360,18.54))
Следовательно, вам нужно либо обновить свои координаты X, добавив 360°, либо заново нарисовать многоугольник в правильном положении. Вы можете обновить координаты вашего многоугольника следующим образом:
st_geometry(WestCoastIndiax) = st_geometry(WestCoastIndiax) + c(360, 0)
st_crs(WestCoastIndiax) = 4326
Обратите внимание, что мы теряем информацию crs при добавлении смещений, поэтому нам нужно переопределить crs.
У меня есть openend запрос проблемы / функции в mapedit, чтобы убедиться, что рисование разрешено только в пределах разумных длинных границ.