在ggmap上的点之间添加网格线

Question

*编辑：图像在第一篇文章中不可见

我是 ggmap 新手。我试图将网格中的特定点绘制到地图上，但我只有网格角点的 GPS 坐标。有没有办法在这些点之间构建 10 X 10 网格，提取这些网格中正方形的中心，并绘制这些点的子集？

这是我的地图现在的样子，具有三个不同网格的四个角。我不想包含代码，因为它包含位置数据。

我不确定如何解决这个问题，或者是否可能。

Answer 1

由于您尚未提供示例数据集，因此此表示假定您的网格坐标位于 WGS84/EPSG:4326 中。如果您的数据有预计的 CRS，您可以跳过

st_transform()

步骤。无论哪种方式，您都需要在

org_crs

定义数据的原始 CRS。

使用投影数据对于准确性很重要。在此示例中，WGS84/EPSG:4326 点与首次投影然后转换回 WGS84/EPSG:4326 点之间的最大差异约为 43mm，因此可能不是什么大问题，但遵循最佳实践是很好的。

根据您的问题和评论，我还假设您希望在“统一的行和列导致不规则大小的单元格”的 11 x 11 网格内有一个 10 x 10 的质心网格。由于网格代表不规则四边形，因此每个质心之间的距离是“相对”恒定的，但仅相对于它所属的横断面而言。我说相对是因为投影回 WGS84/EPSG:4326（我假设您的实地工作需要它？）会导致一些变化。

对于本示例中使用的数据，当数据为 WGS84/EPSG:4326 时，任何给定样线上的点之间的最大距离变化为 7.74e-06m，因此可能不是问题。如果您选择将数据保留在预计的 CRS 中，例如UTM 区域 31N/EPSG:32631，最大变化为 3.55e-21m。

工作流程：

确保包含网格坐标的数据集是一个 sf 对象，每个网格都有一个 id 列（本例中为 id）
定义所需数量的样线质心，例如
```
n <- 10
```
确定坐标的 UTM 区域，如果您的坐标未投影，则需要获得更高的精度。如果您的数据属于多个 UTM 区域，请寻找一种适当的方法来管理这一问题。或者，最简单的解决方案是不投影数据并接受更大程度的不准确性
- 投影网格点（如果认为有必要）
- 定义每个网格的左侧和右侧
- 使用
```
mutate()
```
  和
```
complete()
```
  ，计算左侧点和右侧点之间的坐标序列
- ```
pivot_wider()
```
  并计算从上到下的坐标
```
seq(n+1)
```
- ```
filter()
```
  和
```
select()
```
  所以只保留样线质心
- ```
pivot_longer()
```
  并创建 transect_idx 和 transect_idy 列来定义经度和纬度样线。然后，您可以使用 id、transect_idx 和 transect_idy 列的任意组合来对点进行子集化
- 转换为sf对象并投影到原始CRS
- 创建精度为 n 位的 lon/lat 列

library(dplyr)
library(tidyr)
library(sf)
library(ggplot2)

# Example df with WGS86/EPSG:4326 coordinates
p <- data.frame(id = rep(1:3, each = 4),
                lon = c(1.86, 1.91, 1.9, 1.85, 
                        1.83, 1.79, 1.80, 1.835, 
                        1.72, 1.67, 1.65, 1.68),
                lat = c(1.84, 1.85, 1.9, 1.88,
                        1.87, 1.855, 1.825, 1.835, 
                        1.63, 1.67, 1.64, 1.60))

# Define original CRS (change if your coordinates are not WGS84)
org_crs <- 4326

# Convert p to an sf object (skip if your data are alreadly an sf object with CRS)
p <- st_as_sf(p, coords = c("lon", "lat"), crs = org_crs, remove = FALSE)

# Set number of desired transects
n <- 10

# Function to return UTM Zone EPSG code
utm_epsg <- function(x, y) {
  
  utm_zone <- floor((x + 180) / 6) + 1
  if (y >= 0) {
    epsg <- 32600 + utm_zone
  } else {
    epsg <- 32700 + utm_zone
  }
  
}

# Return EPSG code/s. If more than one, check online for appropriate strategy or
# search for alternative CRS EPSG
unique(mapply(utm_epsg, st_coordinates(p)[,1], st_coordinates(p)[,2]))
# [1] 32631

# Define projection CRS
prj_crs <- 32631

# Project p, return centroids between each x/y transect, project back to original CRS
p1 <- p %>%
  st_transform(prj_crs) %>%
  mutate(lon = st_coordinates(.)[, 1],
         lat = st_coordinates(.)[, 2]) %>%
  group_by(id) %>%
  mutate(pid1 = if_else(lon == min(lon) | lat == min(lat), "l", "r")) %>%
  arrange(id, pid1, desc(lat)) %>%
  mutate(pid2 = rep(c(1, n * 2 + 1), 2)) %>%
  ungroup() %>%
  complete(id, pid1, pid2 = full_seq(pid2, 1)) %>%
  mutate(lon = seq(first(lon), last(lon), length.out = n()),
         lat = seq(first(lat), last(lat), length.out = n()),
         .by = c(id, pid1)) %>%
  pivot_wider(id_cols = c(id, pid1),
              names_from = pid2,
              values_from = c(lon, lat),
              names_sep = "_") %>%
  mutate(pid2 = c(1, n * 2 + 1), .by = id) %>%
  complete(id, pid2 = full_seq(pid2, 1)) %>%
  mutate(across(matches("^l.*_\\d+"), ~ seq(first(.), last(.), length.out = n())),
         .by = id) %>%
  filter(pid2 %% 2 == 0) %>%
  select(-matches("_\\d*[13579]$")) %>%
  pivot_longer(cols = starts_with(c("lon_", "lat_")),
               names_to = c(".value", "transect_idx"),
               names_pattern = "(lon|lat)_(\\d+)") %>%
  select(-starts_with("pid")) %>%
  mutate(transect_idx = as.integer(transect_idx) / 2,
         transect_idy = 1:n(),
         .by = c(id, transect_idx)) %>%
  arrange(id, transect_idx) %>%
  st_as_sf(coords = c("lon", "lat"), crs = prj_crs) %>%
  st_transform(org_crs) %>%
  mutate(lon = st_coordinates(.)[, 1],
         lat = st_coordinates(.)[, 2]) %>%
  mutate(across(c(lat, lon), ~ formatC(., digits = 8, format = "f")))

p1
# Simple feature collection with 300 features and 5 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 1.65255 ymin: 1.6035 xmax: 1.907 ymax: 1.896525
# Geodetic CRS:  WGS 84
# # A tibble: 300 × 6
#       id transect_idx transect_idy         geometry lon        lat       
#  * <int>        <dbl>        <int>      <POINT [°]> <chr>      <chr>     
#  1     1            1            1 (1.853 1.878975) 1.85299998 1.87897503
#  2     1            1            2 (1.858 1.880925) 1.85799991 1.88092506
#  3     1            1            3 (1.863 1.882875) 1.86299985 1.88287508
#  4     1            1            4 (1.868 1.884825) 1.86799982 1.88482510
#  5     1            1            5 (1.873 1.886775) 1.87299980 1.88677511
#  6     1            1            6 (1.878 1.888725) 1.87799980 1.88872511
#  7     1            1            7 (1.883 1.890675) 1.88299982 1.89067510
#  8     1            1            8 (1.888 1.892625) 1.88799986 1.89262508
#  9     1            1            9 (1.893 1.894575) 1.89299991 1.89457506
# 10     1            1           10 (1.898 1.896525) 1.89799999 1.89652503
# # ℹ 290 more rows
# # ℹ Use `print(n = ...)` to see more rows

在ggmap上的点之间添加网格线

问题描述投票：0回答：1

1个回答

最新问题

在ggmap上的点之间添加网格线

问题描述 投票：0回答：1

1个回答

最新问题

问题描述投票：0回答：1