如何将R中插值结果导出为带预定义配色的GeoTIFF用于OpenLayers

项目背景 #

我们在一项生态研究项目中，采用R语言的反距离加权（Inverse Distance Weighting, IDW）方法对雪样污染物分析结果进行插值处理，最终需要将插值结果导入使用OpenLayers库的Web应用。

现有代码实现 #

数据读取与预处理 #

library(sf)
library(stars)
library(tmap)
library(raster)
library(plotly)
library(mapview)
library(tidyverse)
library(ggrepel)
library(dismo) 
library(akima) 
library(gstat) 
library(deldir)

options(scipen=999)

pts = read_table2("D:/www/baikal_so4.dat") %>% st_as_sf(coords = c('x', 'y'), crs = 4326, remove = FALSE)
pts = st_transform(pts, 3857)

coords = st_coordinates(pts)
box = st_bbox(pts)
envelope = box[c(1,3,2,4)]

px_grid = st_as_stars(box, dx = 500, dy = 500)

# 提取网格点坐标
coords_grid = st_coordinates(px_grid)

# 初始化配色映射
stat_colors = colorRampPalette(c('#7EBEDE', '#93C4D9','#A1C9D4','#B2D2CE', '#B9D4C6','#C3D9BB','#D2E3B9','#DDE8B0','#EEF2B1','#F9FBA5','#FDF798','#FDE791','#FBD992','#F9CB8A','#F8BF86','#F5B180','#F5A57C','#F39779','#F18272','#EF7474'))

stat_levels = seq(0, 3, by=0.2)
stat_ncolors = length(stat_levels)-1
           
stat_legend = scale_fill_manual(name = 'mg/l', values = stat_colors(stat_ncolors), guide = guide_legend(label.vjust = -0.3, reverse = TRUE, title.position = "bottom"), labels = stat_levels, na.value = 'white', drop = FALSE)

stat_mapping = aes(fill = cut(values, breaks = stat_levels))

# Delaunay三角化
edges = pts %>% 
  st_union() %>% 
  st_triangulate()

IDW插值与可视化 #

px_grid = px_grid %>% 
  mutate(z_idw2 = gstat::idw(values ~ 1, locations = pts, newdata = px_grid, idp = 2.0) %>% pull(var1.pred) %>% as('vector'),
         z_idw3 = gstat::idw(values ~ 1, locations = pts, newdata = px_grid, idp = 3.0) %>% pull(var1.pred) %>% as('vector'),
         z_idw4 = gstat::idw(values ~ 1, locations = pts, newdata = px_grid, idp = 4.0) %>% pull(var1.pred) %>% as('vector'),
         z_idw5 = gstat::idw(values ~ 1, locations = pts, newdata = px_grid, idp = 5.0) %>% pull(var1.pred) %>% as('vector'))
         
cont_idw2 = st_contour(px_grid['z_idw2'], 
                       breaks = stat_levels, 
                       contour_lines = TRUE)

cont_idw3 = st_contour(px_grid['z_idw3'], 
                       breaks = stat_levels, 
                       contour_lines = TRUE)

cont_idw4 = st_contour(px_grid['z_idw4'], 
                       breaks = stat_levels, 
                       contour_lines = TRUE)

cont_idw5 = st_contour(px_grid['z_idw5'], 
                       breaks = stat_levels, 
                       contour_lines = TRUE)

# 生成可视化图
ggplot() +
  geom_stars(data = cut(px_grid['z_idw5'], breaks = stat_levels)) +
  stat_legend +
  coord_sf(crs = st_crs(pts)) +
  geom_sf(data = cont_idw5, color = 'black', size = 0.2) +
  geom_sf(data = pts, color = 'red', size = 0.5) +
  geom_sf(data = edges, color = 'red', size = 0.1, fill = NA)

问题描述 #

直接使用write_stars导出px_grid['z_idw5']为GeoTIFF时，文件仅保留了插值数值，丢失了预定义的渐变配色方案，导致在Web应用中无法还原R中可视化的颜色效果。

解决方案 #

要导出带有预定义配色的GeoTIFF，需要将连续插值值转换为离散分类值，并为GeoTIFF附加颜色查找表（Color Lookup Table, CLUT）。具体步骤如下：

1. 对插值结果进行分箱分类 #

将连续的z_idw5值按照预定义的stat_levels转换为离散分类：

# 对z_idw5进行分箱，得到分类后的stars对象
px_grid_classified = cut(px_grid['z_idw5'], breaks = stat_levels)

2. 转换为Raster对象并设置颜色表 #

利用raster包处理颜色表，将分类值与预定义颜色绑定：

# 将stars对象转换为Raster对象
r = raster(px_grid_classified)

# 创建颜色表：将每个分类对应到预定义颜色
rat = levels(r)[[1]]
rat$Color = stat_colors(stat_ncolors)
# 为分类添加数值区间标签（可选，方便后续识别）
rat$Level = paste(stat_levels[-length(stat_levels)], stat_levels[-1], sep = "-")
levels(r) = rat

3. 导出带颜色表的GeoTIFF #

使用writeRaster导出，确保颜色表被写入文件：

# 导出为带颜色表的GeoTIFF
writeRaster(
  r, 
  filename = "idw5_classified_colored.tif", 
  format = "GTiff", 
  overwrite = TRUE, 
  datatype = "INT1U", # 分类数<=255，用8位无符号整数即可
  options = c("COMPRESS=LZW") # 可选，压缩文件
)

关键说明 #

• 直接导出连续数值的GeoTIFF时，配色是由可视化工具动态生成的，无法保留R中预定义的配色；而将数值转换为离散分类并附加颜色表后，GeoTIFF会自带颜色映射规则，OpenLayers可以直接读取并应用该配色。
• 选择INT1U数据类型是因为我们的分类数为15（stat_levels从0到3，步长0.2，共15个区间），远小于255，用8位整数足够存储，同时能减小文件体积。

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