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

【问题标题】：Project longitude and latitude to planar coordinates system将经度和纬度投影到平面坐标系
【发布时间】：2020-01-07 07:46:42
【问题描述】：

免责声明：我目前正在使用 {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 中的多边形作为窗口。问题是我找不到正确的坐标参考系统 (CRS)，因为我的点位于多边形内。我收到一条错误消息，告诉我这些点位于窗外。这是使下面的代码可重现的两个文件：

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)

【问题讨论】：

为什么你还没有将vz 转换为你已经将area_one_flat 转换为的坐标系？

标签： r coordinate-systems sf spatstat

【解决方案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))

【讨论】：