R语言栅格化海岸多边形后采样点无法正常绘制问题求解

我们拟制作Vermeille海岸地图，用于计算两点间路径不穿越陆地的采样点间距，操作流程如下：

已完成操作步骤 #

• 1. 拼接两份海岸shapefile数据
     参考方法：R cran: sf 拼接两类线要素
• 2. 生成外接框构建多边形
     参考方法：sf包：基于复杂线要素闭合生成多边形
• 3. 完成多边形栅格化
     参考方法：R语言sf多边形栅格化方法

示例数据下载地址：海岸测试数据

前期数据读取、拼接、多边形生成代码如下：

library(sf)
library(fasterize)
library(raster)
library(dplyr)
library(tidyverse)

frenchCoast_CoteBanyuls <- st_read("coasts_subnational_France/coasts_subnational.shp")
spainCoast_CoteBanyuls <- st_read("coasts_subnational_Spain/coasts_subnational.shp")

spainurl <- "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:coasts_subnational&outputformat=SHAPE-ZIP&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Emrgid_1%3C%2FPropertyName%3E%3CLiteral%3E3417%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E"
download.file(spainurl, "spain.zip", mode = "wb")
unzip("spain.zip", exdir = "spain", junkpaths = TRUE)

franceurl <- "https://geo.vliz.be/geoserver/wfs?request=getfeature&service=wfs&version=1.0.0&typename=MarineRegions:coasts_subnational&outputformat=SHAPE-ZIP&filter=%3CPropertyIsEqualTo%3E%3CPropertyName%3Emrgid_1%3C%2FPropertyName%3E%3CLiteral%3E19888%3C%2FLiteral%3E%3C%2FPropertyIsEqualTo%3E"
download.file(franceurl, "france.zip", mode = "wb")
unzip("france.zip", exdir = "france", junkpaths = TRUE)

spainCoast_CoteBanyuls <- list.files("spain",
                                     pattern = "shp$",
                                     full.names = TRUE) %>% st_read()
frenchCoast_CoteBanyuls <- list.files("france",
                                      pattern = "shp$",
                                      full.names = TRUE) %>% st_read()


lines_spain <- st_geometry(spainCoast_CoteBanyuls) %>% st_cast("LINESTRING")
spainCoast_l <- st_sf(n = as.character(seq_len(length(lines_spain))), lines_spain)

lines_france <- st_geometry(frenchCoast_CoteBanyuls) %>% st_cast("LINESTRING")
franceCoast_l <- st_sf(n = as.character(seq_len(length(lines_france))), lines_france)

spainmax <- spainCoast_l[which.max(st_length(spainCoast_l)), ]
spainrest <- spainCoast_l[-which.max(st_length(spainCoast_l)), ]

joined <- c(st_geometry(spainmax), st_geometry(franceCoast_l)) %>%
  st_union()

join_end <- st_union(joined, st_geometry(spainrest))

bbox_all <- st_bbox(joined) %>%
  st_as_sfc()

polygon_joined <- bbox_all %>%
  lwgeom::st_split(join_end) %>%
  st_collection_extract("POLYGON")

# 目视检查后移除位置2、3的多余多边形
polygon_end <- polygon_joined[2] # 定义陆地区域多边形而非海域

polyCombin_df <- st_sf(var = 1, polygon_end)
class(polyCombin_df)
st_crs(25831)$units
polyCombin_df_t <- polyCombin_df %>% st_transform(25831)

该步骤得到的多边形结果正常：

随后执行多边形栅格化操作：

r <- raster(polyCombin_df_t, res = 100)
r <- fasterize(polyCombin_df_t, r, fun = "max")
par(mar=c(1,1,1,1))
plot(r)

栅格化结果显示正常：

• 4. 采样点间距计算（待完成）
     后续计划参考R语言距离计算方法（参考教程：R语言三类距离计算方法）计算采样点间距，首先尝试使用points函数在图上添加沿岸采样点坐标，代码如下：

# sites
site_random <- matrix(data = c(3.164887 , 3.123969 , 3.158125 , 3.160378, 42.402158, 
                42.521957, 42.475956, 42.461188), ncol = 2)

site_random <- data.frame(site_random) 
points(site_random$X1, site_random$X2, pch = 19)

执行后无采样点显示，初步怀疑为绘图比例尺或坐标系不匹配问题。

问题原因与修复方案 #

核心原因 #

采样点无法显示完全是坐标系不匹配导致：

1. 栅格化前已将多边形转换为EPSG:25831（UTM 31N投影坐标系），生成的栅格坐标单位为米，绘图窗口的坐标范围为投影后的米制坐标，东距约400000-520000m、北距约4690000-4720000m
2. 直接传入points()的采样点为WGS84经纬度坐标，数值范围仅为经度3°左右、纬度42°左右，和绘图窗口坐标范围差6个数量级，点被绘制在可视区域外，无法显示。

修复代码 #

将采样点转换为与栅格一致的投影坐标系后再绘图即可：

# 整理采样点
site_random <- matrix(data = c(3.164887 , 3.123969 , 3.158125 , 3.160378, 42.402158, 
                               42.521957, 42.475956, 42.461188), ncol = 2)
# 转为sf对象，指定原始坐标为WGS84经纬度(EPSG:4326)
site_random_sf <- st_as_sf(data.frame(lon = site_random[,1], lat = site_random[,2]),
                           coords = c("lon", "lat"), crs = 4326)
# 转换到和栅格一致的EPSG:25831坐标系
site_random_proj <- st_transform(site_random_sf, crs = st_crs(r))
# 提取投影后坐标
site_coords <- st_coordinates(site_random_proj)

# 绘图
par(mar=c(1,1,1,1))
plot(r)
points(site_coords[,1], site_coords[,2], pch = 19, col = "red", cex = 1.2)

后续计算不穿越陆地的最短路径时，直接使用投影后的采样点坐标即可，避免坐标系不一致带来的计算误差。

内容的提问来源于stack exchange，提问作者Florian B.