将经纬度投影到平面坐标系

Question

免责声明：我目前使用的是 {sf} 包的 ppp 分支版本，因为它提供了在 {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 是使用从互联网上抓取的经纬度坐标创建的。我的目标是创建一个 ppp 对象，其中 vz 中的位置作为点，area_one 中的多边形作为 window。问题是我找不到正确的坐标参考系统 (CRS)，因为我的点位于多边形内。我收到一条错误消息，告诉我这些点位于 window 之外。以下是使下面的代码可重现的两个文件：

area_one: download here

vz: download here

# 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)

# 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)

Answer 1

正如@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))

将经纬度投影到平面坐标系

Project longitude and latitude to planar coordinates system

r

coordinate-systems

sf

spatstat