Проект долготы и широты для плоской системы координат - PullRequest
Новая
       26

Проект долготы и широты для плоской системы координат

1 голос
/ 07 января 2020

Отказ от ответственности: в настоящее время я использую ветвь ppp версии пакета {sf}, поскольку в нем доступны новые функции для преобразования объектов между {sf} и {spatstat} (см. https://github.com/r-spatial/sf/issues/1233). Чтобы он работал правильно, мне пришлось вручную удалить пакет {sf} со своего жесткого диска, а затем переустановить его с Github. Я также использую версию разработки {spatstat} без особой причины. 

# install.packages("remotes")
install_github("r-spatial/sf@ppp")
install_github("spatstat/spatstat")

У меня есть два геопространственных объекта: area_one, который представляет собой объединение многоугольников нескольких округов в Техасе и vz, который является точечным расположением нескольких магазинов в Техасе, и оба они являются объектами семейства sf vz было создано с использованием координат долготы и широты, вырезанных из целого rnet. Моя цель - создать объект ppp с местоположениями в vz в качестве точек и многоугольником в area_one в качестве окна. Проблема в том, что я не могу найти правильную систему координат (CRS) для моих точек на l ie внутри многоугольника. Я получаю сообщение об ошибке, сообщающее, что точки l ie за окном. Вот два файла, чтобы сделать код ниже воспроизводимым:

area_one: скачать здесь

vz: скачать здесь 

# Load packages

library(sf) # Development version in the ppp branch
library(spatstat) # Development version in the master branch
library(tmap)
library(here)

# Read the geospatial data (CRS = 4326)

area_one <- st_read(dsn = here("area_one/area_one.shp"), layer = "area_one")
vz <- st_read(dsn = here("vz/vz.shp"), layer = "vz")

# Plot a quick map

tm_shape(area_one) +
  tm_borders() +
  tm_shape(vz) +
  tm_dots(col = "red", size = 0.1)

enter image description here

# Create a planar point pattern

vz_lonlat <- st_coordinates(vz)
area_one_flat <- st_transform(area_one, crs = 6345)

p <- ppp(x = vz_lonlat[, "X"], y = vz_lonlat[, "Y"], window = as.owin(area_one_flat)) # Note that the reason why this line of code does not throw an error (especially despite the window argument being an sf object) is because of the version of the {sf} package I am using

Warning message:
49 points were rejected as lying outside the specified window 

plot(p)

enter image description here

1 Ответ

1 голос
/ 07 января 2020

Как указывает @spacedman, сначала вы должны преобразовать vz в ту же систему координат, что и область наблюдения. Я думаю, вы могли бы сделать что-то вроде (не проверено): 

vz_flat <- st_coordinates(st_transform(vz, crs = 6345))
area_one_flat <- st_transform(area_one, crs = 6345)
p <- ppp(x = vz_flat[, "X"], y = vz_flat[, "Y"], window = as.owin(area_one_flat))
...