如何将经纬度点添加至Landsat影像RasterBrick可视化图中？

1. Landsat影像RasterBrick信息 #

landsat_stack_brick
#class      : RasterBrick 
#dimensions : 7711, 7551, 58225761, 7  (nrow, ncol, ncell, nlayers)
#resolution : 30, 30  (x, y)
#extent     : 576885, 803415, 2281485, 2512815  (xmin, xmax, ymin, ymax)
#crs        : +proj=utm +zone=4 +datum=WGS84 +units=m +no_defs 
#source     : r_tmp_2022-11-19_144142_2384_62790.grd 
#names      :    B1,    B2,    B3,    B4,    B5,    B6,    B7 
#min values :     1,     1,     1,     4,  1780,  6138,  7110 
#max values : 60301, 61806, 61527, 62936, 62293, 65454, 65454 

2. 初始可视化尝试代码 #

library(raster)

Long_lat_df  = na.omit(water_long_lat)
coordinates(Long_lat_df ) <- ~ water_long + water_lat
proj4string(Long_lat_df ) <- CRS("+init=epsg:4326")
y <- spTransform(geo_points, crs(landsat_stack_brick))

plot(landsat_stack_brick)
points(y)

3. 问题现象 #

上述代码添加的经纬度点本应分布在岛屿海岸线附近，但实际显示位置异常，未出现在预期区域。

4. 修正后的尝试代码 #

water_chem <- read.csv("data.csv")

###提取所有唯一纬度值
water_lat = unique(water_chem$Latitude)
###提取对应唯一纬度的经度值
water_long= water_chem$Longitude[unique(water_chem$Latitude)]

x <- water_long_lat
x = na.omit(x)
coordinates(x) <- ~ water_long + water_lat
sbux_sf <- st_as_sf(x, coords = c("water_long", "water_lat"),  crs = landsat_stack@crs)

sbux_sd

structure(list(geometry = structure(list(structure(c(-156.63776, 
21.0133265), class = c("XY", "POINT", "sfg")), structure(c(-156.63776, 
21.0134451), class = c("XY", "POINT", "sfg")), structure(c(-156.63776, 
21.0135265), class = c("XY", "POINT", "sfg")), structure(c(-156.63776, [...]

[...], precision = 0, bbox = structure(c(xmin = -156.6407, ymin = 20.77606, 
xmax = -156.63776, ymax = 21.014317), class = "bbox"), crs = structure(list(
    input = "WGS 84", wkt = "GEOGCRS[\"WGS 84\",
    DATUM[\"World Geodetic System 1984\",
        ELLIPSOID[\"WGS 84\",6378137,298.257223563,
            LENGTHUNIT[\"metre\",1]],
        ID[\"EPSG\",6326]],
    PRIMEM[\"Greenwich\",0,
        ANGLEUNIT[\"degree\",0.0174532925199433],
        ID[\"EPSG\",8901]],
    CS[ellipsoidal,2],
        AXIS[\"longitude\",east,
            ORDER[1],
            ANGLEUNIT[\"degree\",0.0174532925199433,
                ID[\"EPSG\",9122]]],
        AXIS[\"latitude\",north,
            ORDER[2],
            ANGLEUNIT[\"degree\",0.0174532925199433,
                ID[\"EPSG\",9122]]],
    USAGE[
        SCOPE[\"unknown\"],
        AREA[\"World.\"],
        BBOX[-90,-180,90,180]]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
124L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(integer(0), class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = character(0)))

问题根源 #

1. 坐标系设置错误：修正代码中错误地将点数据的初始CRS设置为Landsat影像的UTM坐标系，但点数据实际是WGS84地理坐标系（EPSG:4326），直接匹配导致坐标完全错位。
2. 经纬度提取逻辑错误：water_long= water_chem$Longitude[unique(water_chem$Latitude)] 无法保证经纬度一一对应，unique()返回的是纬度唯一值而非索引，会导致坐标配对混乱。

正确实现代码 #

# 加载依赖包
library(raster)
library(sf)
library(dplyr)

# 读取并清理数据：提取唯一且配对的经纬度
water_chem <- read.csv("data.csv")
cleaned_points <- water_chem %>%
  distinct(Latitude, Longitude, .keep_all = TRUE) %>%  # 确保经纬度一一对应
  na.omit()  # 移除缺失值

# 将经纬度转为sf对象，初始CRS设为WGS84（EPSG:4326）
points_sf <- st_as_sf(cleaned_points,
                      coords = c("Longitude", "Latitude"),  # 注意顺序：经度在前，纬度在后
                      crs = 4326)

# 将点转换到Landsat影像的UTM坐标系
points_utm <- st_transform(points_sf, crs = crs(landsat_stack_brick))

# 可视化：先绘制影像，再叠加转换后的点
plot(landsat_stack_brick)
plot(st_geometry(points_utm), add = TRUE, col = "red", pch = 16, cex = 0.8)

关键注意事项 #

• 坐标顺序：sf/sp包要求坐标系转换时，先传入经度，再传入纬度，顺序错误会导致点位置偏移。
• 初始CRS必须匹配点数据的实际坐标系：这里点数据是WGS84地理坐标，必须先设置为EPSG:4326，再转换到影像的UTM坐标系。
• 数据清理时必须保证经纬度配对：使用distinct(Latitude, Longitude, .keep_all = TRUE)确保每一组经纬度都是唯一且对应的，避免索引错误导致坐标混乱。

内容的提问来源于stack exchange，提问作者Pablo Rodriguez