在R中处理Sentinel-3 300米分辨率数据提取TSM遇阻求助

• 处理300米分辨率Sentinel-3数据，目标提取每个像元的TSM（总悬浮物浓度）
• 数据以归档形式提供，包含多个独立netCDF文件：geo_coordinates.nc存储经纬度坐标，tsm_nn.nc存储TSM数据
• 使用stars包批量读取文件后，经纬度仅作为对象属性存在，未被识别为空间维度，导致：
  • 无投影信息
  • 绘图时南北/东西方向颠倒，边界框错误，地理区域位置偏移

用户尝试的代码：

library(stars)

nc_files <- c("S3B_OL_2_WFR____20230218T101125_20230218T101425_20230219T222151_0179_076_179_2160_MAR_O_NT_003.SEN3/geo_coordinates.nc", 
              "S3B_OL_2_WFR____20230218T101125_20230218T101425_20230219T222151_0179_076_179_2160_MAR_O_NT_003.SEN3/tsm_nn.nc")

lon_lat_tsm <- read_stars(nc_files)

运行结果摘要：

stars object with 2 dimensions and 5 attributes
attribute(s), summary of first 1e+05 cells:
                   Min.   1st Qu.      Median         Mean     3rd Qu.        Max.  NA's
latitude [°]  39.568894 40.402587  41.0897840  40.99148193  41.6207120   42.046713     0
longitude [°] -7.581946 -3.745773   0.0337230   0.05383067   3.8784562    7.686906     0
TSM_NN        -1.909408 -1.148438  -0.9310176  -0.84433249  -0.5505323    2.438995 56460
dimension(s):
  from   to offset delta refsys point values x/y
x    1 4865      0     1     NA    NA   NULL [x]
y    1 4091   4091    -1     NA    NA   NULL [y]

核心思路是手动将经纬度数据关联为TSM数据的空间维度，并设置正确坐标系。

步骤1：分别读取坐标与TSM数据 #

library(stars)
library(sf)

# 读取坐标文件
geo_data <- read_stars("S3B_OL_2_WFR____20230218T101125_20230218T101425_20230219T222151_0179_076_179_2160_MAR_O_NT_003.SEN3/geo_coordinates.nc")

# 读取TSM文件
tsm_data <- read_stars("S3B_OL_2_WFR____20230218T101125_20230218T101425_20230219T222151_0179_076_179_2160_MAR_O_NT_003.SEN3/tsm_nn.nc")

步骤2：构建带地理坐标的stars对象 #

从坐标文件中提取经纬度矩阵，替换TSM数据的默认索引维度：

# 提取经纬度矩阵
lon_matrix <- geo_data[["longitude"]]
lat_matrix <- geo_data[["latitude"]]

# 用经纬度作为空间维度重构TSM对象
tsm_geo <- st_as_stars(
  list(TSM_NN = tsm_data[["TSM_NN"]]),
  dims = st_dimensions(x = lon_matrix, y = lat_matrix, point = TRUE)
)

步骤3：设置坐标系 #

Sentinel-3的地理坐标默认采用WGS84（EPSG:4326）：

st_crs(tsm_geo) <- 4326

步骤4：验证与提取数据 #

• 查看修正后的对象：

print(tsm_geo)

• 绘图验证地理正确性：

plot(tsm_geo, main = "TSM 浓度分布")

• 提取所有像元的经纬度与TSM值（转为DataFrame）：

tsm_pixel_data <- as.data.frame(tsm_geo, xy = TRUE)
head(tsm_pixel_data)

额外处理：若仍存在南北方向颠倒 #

如果绘图时地理区域仍反向，可反转y轴维度：

# 反转y轴顺序
tsm_geo_fixed <- tsm_geo[, rev(dim(tsm_geo)[[2]]$values)]
plot(tsm_geo_fixed)

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