如何用R将Sentinel-3 LST NetCDF数据转换为带坐标的GeoTIFF？

问题背景 #

Sentinel-3 SLSTR的LST（地表温度）数据拆分存储在两个NetCDF文件中：LST_in.nc存储温度值，geodetic_in.nc存储经纬度信息。使用R转换为GeoTIFF时结果异常，但SNAP工具处理的栅格结果正常。尝试代码如下：

dir <- "/home/user/S3A_SL_2_LST____20221125T125014_20221125T125314_20221126T214452_0179_092_266_3060_PS1_O_NT_004.SEN3/"

output_raster = "20221126T214452_0179_092_266_3060_PS1_O_NT_004"  

# Loading libraries. 
library(ncdf4)
library(raster)
library(dplyr)
library(ggplot2)
library(terra)

# Creating filepath names
climate_filepath <- paste0(dir, "LST_in.nc")
cart_filepath <- paste0(dir, "geodetic_in.nc")

# Reading them in using nc_open
nc <- nc_open(climate_filepath)
cord <-  nc_open(cart_filepath)

# All three files have a 1200 x 1500 cell matrix. Thus, I collapsed the matrix, and bound them all into a dataframe:
latitude <- ncvar_get(cord, "latitude_in") %>% as.vector()
longitude <- ncvar_get(cord, "longitude_in") %>% as.vector()
lst <- ncvar_get(nc, "LST") %>% as.vector()

LST_DF = data.frame(lon = longitude,
                    lat = latitude,
                    LST = lst) %>% 
  #Convert from Kelvin to Celcius
  dplyr::mutate(LST = LST - 273.15)

# I converted all the variables in a data frame to a matrix
LST_DF_matrix <- data.matrix(LST_DF, rownames.force = NA)
colnames(LST_DF_matrix) <- c('X', 'Y', 'Z')
head(LST_DF_matrix)

# Set up a raster geometry, using terra package
e <- ext(apply(LST_DF_matrix[,1:2], 2, range))
# Set up the raster. I choose this ncol and nrow, however, I don't know if it is correct. 
# The dimension of the data is 1200, 1500, 1800000  (nrow, ncol, ncell) with the 1 km grid resolution. In another attempt, I opened each NetCDF as a raster too.
r <- rast(e, ncol=1000, nrow=1000)

# I apply rasterize
x <- rasterize(LST_DF_matrix[, 1:2], r, LST_DF_matrix[,3]) #, fun=mean)
plot(x) #,  col = rev(rainbow(25)))

# Set CRS    
crs(x) = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# Saving to a GeoTIFF
writeRaster(x = x, filename = paste0(dir, output_raster, "_v3.tif"), overwrite=TRUE)

错误原因 #

1. 栅格维度错误：手动设置ncol=1000, nrow=1000，与原始数据的1200×1500网格不匹配，破坏了数据的空间结构。
2. 错误使用栅格化方法：Sentinel-3 SLSTR LST是规则网格数据，无需用rasterize（针对离散点转栅格），直接基于原始矩阵构建栅格即可。
3. 行列对应关系丢失：将矩阵转为向量合并成数据框，导致经纬度与温度值的网格对应关系错位。

修正后的代码 #

dir <- "/home/user/S3A_SL_2_LST____20221125T125014_20221125T125314_20221126T214452_0179_092_266_3060_PS1_O_NT_004.SEN3/"
output_raster <- "20221126T214452_0179_092_266_3060_PS1_O_NT_004"  

# 加载核心库
library(ncdf4)
library(terra)

# 读取NetCDF文件
nc_lst <- nc_open(paste0(dir, "LST_in.nc"))
nc_geo <- nc_open(paste0(dir, "geodetic_in.nc"))

# 读取原始矩阵（保留行列结构），直接转换为摄氏度
lst_matrix <- ncvar_get(nc_lst, "LST") - 273.15
lat_matrix <- ncvar_get(nc_geo, "latitude_in")
lon_matrix <- ncvar_get(nc_geo, "longitude_in")

# 关闭NetCDF连接，释放资源
nc_close(nc_lst)
nc_close(nc_geo)

# 获取栅格的空间参数
x_range <- range(lon_matrix)
y_range <- range(lat_matrix)
nrows <- nrow(lst_matrix)
ncols <- ncol(lst_matrix)

# 创建匹配原始数据的空栅格
r <- rast(
  nrows = nrows,
  ncols = ncols,
  xmin = x_range[1],
  xmax = x_range[2],
  ymin = y_range[1],
  ymax = y_range[2],
  crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
)

# 转置矩阵后赋值，匹配Terra的行优先存储顺序
values(r) <- as.vector(t(lst_matrix))

# 保存为GeoTIFF
writeRaster(r, filename = paste0(dir, output_raster, "_fixed.tif"), overwrite = TRUE)

# 预览结果
plot(r)

关键修正说明 #

• 保留原始网格结构：直接读取NetCDF中的矩阵数据，避免向量转换导致的坐标与值错位。
• 使用原始维度构建栅格：基于原始数据的1200×1500行列数创建栅格，保证空间分辨率与原始数据一致。
• 直接赋值替代栅格化：将LST矩阵转置后转为向量赋值给栅格，匹配Terra包的栅格存储逻辑。
• 提前设置坐标系：创建栅格时直接指定WGS84坐标系，避免后续修改可能引发的问题。

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