企业官网建设流程全解析

做大气碳浓度这块的人，手里大概率都躺着一堆.nc文件。它们可能来自一颗专门测温室气体的遥感卫星 OCO-2，文件名长得像oco2_L3CO2_day_20150101_v10r.nc这种。你双击打不开，拖进 GIS 软件要么一片空白，要么只在角落里出现一堆奇怪的散点。想把它们变成大家熟知的.tif，好做切片、好嵌到 Web 地图、好丢给同事在别的软件里看一眼——这就是"OCO-2 的 NC 格式转 tif"这件事的全部动机。

但它远没有听起来那么简单。OCO-2 的数据家族里，有的产品本身是规则格网，有的产品本质上是沿轨道的成千上万个观测点；有的变量是浓度柱，有的是廓线，有的是质量标记。闭着眼睛写一句gdal_translate in.nc out.tif，得到的十有八九是错的，而且错得很隐蔽——图能出来，位置偏了半个像元，或者值域整个错位，你不细看根本发现不了。这篇内容我会把自己在这一路上摸清的结构、算过的参数、翻过的车，从头到尾捋一遍，写成可以直接抄的流程，适合刚接触这类数据的学生、做遥感应用的工程师，以及需要把科学数据接进业务系统的朋友。

1. 动手之前，先认清手里的 NC 到底是什么

1.1 OCO-2 的 NC 文件并不都是同一种东西

很多人对"遥感数据"有一个默认假设：既然是卫星拍的，那它一定是一张图。这个假设在 OCO-2 上会让你吃大亏。它的产品线大体分三层，转换方式完全不同。

**一级产品（L1B）**记录的是仪器探测到的光谱，每一个"footprint"对应一条光谱曲线，附带观测几何、太阳天顶角、地面高程等信息。这类数据转 tif 意义不大，因为它的核心内容根本不是一个空间分布的物理量。

二级产品（L2）是从光谱反演出来的结果，最有名的是 Lite 版本的柱浓度产品，比如oco2_L2_Lite_FP_*.nc。它的核心变量是xco2，单位是 ppm，但关键在于——它的维度是sounding_id，也就是"第多少次观测"，不是lat × lon。换句话说，它是一个点表，不是栅格。每个 sounding 有自己的经纬度、时间、质量标记。

**三级产品（L3）**才是真正的格网化结果，官方已经把散点按一定分辨率分箱平均过了，比如oco2_L3CO2_day_*.nc里，xco2的维度是(lat, lon)，经纬度也是规整的等间距向量。

我见过不止一个人把 L2 Lite 硬往gdal_translate里塞，然后抱怨"出来的图全是噪声"。原因很简单：你把一列 10 万个按时间排序的点，按顺序铺成了一张 316×316 的方阵，那不是地图，那是一张随机噪声图。

提示：转换前第一件事，永远是确认维度名。看到sounding_id、footprint、obs_id这种维度，就当点数据对待；看到成对的lat/lon二维结构，才走栅格路线。

1.2 NetCDF 装的是"多维数组 + 元数据"，GeoTIFF 装的是"一张带坐标的图"

理解这两个格式的差别，能帮你避免后面 80% 的错误。

NetCDF 本质上是一个自描述的容器：里面可以塞任意多个变量，每个变量有自己的维度、类型、单位、缩放因子、缺测值。地理位置以"坐标变量"（coordinate variable）的形式存在，比如一个叫lat的一维数组，和一个叫lon的一维数组。软件需要自己把这些拼起来，才知道第 i 行第 j 列对应地球上的哪个点。

GeoTIFF 则强制要求一件事：你必须提供一个仿射变换（geotransform），用 6 个数字把像素坐标和地理坐标线性地绑死。这 6 个数分别是：

问题的核心就在 GT(0)、GT(3)、GT(5) 这三个数上。NC 文件里给你的lat/lon数组，绝大多数情况下是像素中心坐标，你需要往回推半个像元才能得到边界；而且lat数组多半是从南到北升序排列，GeoTIFF 的第一行却是最北边，这个方向差必须处理。

1.3 为什么一句 gdal_translate 经常不灵

GDAL 的 netCDF 驱动确实能读这类文件，但它依赖一套约定：变量得带coordinates属性指向坐标变量，得有grid_mapping描述投影，维度顺序得是(time, lat, lon)这种标准排列。文件里稍有偏差，它要么只返回一个莫名其妙的 band，要么把整个文件识别成一堆 subdataset。

更麻烦的是，即便它能读出来，缺测值、缩放因子、单位这些信息在转换过程中常常被丢掉，你最终拿到的 tif 数值可能完全对不上原始物理量。

我的建议很直接：别指望一行命令搞定，老老实实用脚本显式地构造 geotransform。多写二十行代码，换来的是可复现、可核查、可批量。

2. 转换前的侦察：三分钟摸清文件结构

2.1 命令行体检，先看骨架

拿到一个新文件，我第一步永远是gdalinfo：

gdalinfo oco2_L3CO2_day_20150101_v10r.nc

输出里会列出 subdataset 清单，每一个都长这样：

SUBDATASET_1_NAME=NETCDF:"oco2_L3CO2_day_20150101_v10r.nc":xco2 SUBDATASET_1_DESC=[180][360] xco2 (32-bit floating-point)

方括号里的[180][360]就是维度形状。看到二维的[180][360]，基本可以判定这是个规则格网；看到[1][123456]这种一维长条，就是点数据。

如果gdalinfo报错或者输出过于混乱，就换ncdump：

ncdump -h oco2_L3CO2_day_20150101_v10r.nc | head -n 60

-h只打印头信息（维度和变量声明），不会刷屏。这一眼能看清维度、变量名、类型、属性，效率极高。

2.2 Python 里做交互式探查更灵活

命令行只能看结构，要看实际数值还得进 Python。我一般用xarray，一句话就能把整个文件的"数据模型"摊开：

import xarray as xr ds = xr.open_dataset("oco2_L3CO2_day_20150101_v10r.nc") print(ds)

输出会明确告诉你Dimensions: (lat: 180, lon: 360)，以及Coordinates: lat, lon，还有每个变量的dtype、单位。

接着验证坐标向量是不是等间距：

import numpy as np lat = ds["lat"].values lon = ds["lon"].values print(lat[:3], lat[-3:]) print(lon[:3], lon[-3:]) print(np.diff(lat)[:5]) # 确认步长一致 print(np.diff(lon)[:5])

这一步别省。我就遇到过一次，某个中间版本的文件里lat用float32存储，边缘值有极微小的抖动，np.diff出来的结果差了 1e-6 量级。数值上无所谓，但如果你的代码里写死了assert res == 0.5，脚本直接崩。

2.3 必须盯住的四类变量

一个 L3 文件里变量动辄几十上百个，转换时其实只需要盯四类：

第一类是坐标，lat和lon，决定 geotransform。

第二类是物理量本身，比如xco2（干空气柱平均二氧化碳混合比，ppm）。这类变量是你要真正转出去的东西。

第三类是质量标记，xco2_quality_flag最常用，0表示通过质量筛选，非零表示有各种问题（云污染、地表反照率不足、气溶胶过厚等）。做正式出图时不过滤质量，出来的图会有一大片明显离谱的高值。

第四类是缺测值标识。OCO-2 家族用的填充值是-999999.0，不是-9999，不是1e30，也不是NaN。这个值如果在转换时没清掉，会直接把你的色带拉伸拉到天上去，整张图变成一片单色。

提示：-999999这个数字很大，用data > -1e30或者data != -999999.0都能筛，但后者更稳，因为有些变量实际上会把填充值存成-999999.的浮点近似。建议统一用where(data > -900000, np.nan)这种宽松阈值，避免边界情况。

3. 格网型产品的转换全流程

3.1 反推 geotransform 的六个参数

假设lat从 -89.75 到 89.75，步长 0.5；lon从 -179.75 到 179.75，步长 0.5。这些都是像素中心坐标。

那么：

• 经度左边界 =lon[0] - Δlon/2= -179.75 - 0.25 =-180.0
• 纬度上边界 =lat[-1] + Δlat/2= 89.75 + 0.25 =90.0
• GT(1) =0.5（像元宽）
• GT(5) =-0.5（像元高，负值表示从上往下）

这里有两个很容易搞混的点，我单独拎出来说。

第一个是为什么要用lat[-1]而不是lat[0]。因为lat是升序的（南到北），lat[-1]才是最大的那个纬度，也就是北边界。如果你写成lat[0] + Δ/2，得到的会是 -89.5，整个图会被平移到南半球去，误差刚好是 179 度。这种错误在缩略图上不容易看出来，一旦做叠加分析就全乱了。

第二个是为什么 GT(5) 必须为负。GeoTIFF 的像素行号从上往下增大，而地理纬度从下往上增大，两者方向相反。把 GT(5) 设成负数，等于告诉读取软件"每往下一行，纬度减 0.5"。

3.2 行序对齐：数组得翻个身

有了正确的 geotransform，还差一步：数组的行顺序。

xarray读进来的xco2，第 0 行对应lat[0]（最南边），最后一行对应lat[-1]（最北边）。而 GeoTIFF 要求第 0 行是最北边。所以必须做一次上下翻转：

arr = np.flipud(xco2)

这一步漏了，图能不能出来？能。位置对不对？完全反了——南极的数据跑到北极，北极的跑到南极。而且因为 XCO2 的全球分布南北差异不算特别剧烈（大致 395 到 415 ppm 之间），有些人在调色板拉伸下根本看不出来，直到去核对某个具体坐标才傻眼。

顺带说一句，翻转和 geotransform 的符号是两件事，不能互相替代。有人想着"我把 GT(5) 设成正数不就等于翻了"，那是把像元方向定义改错了，输出的 tif 在别的软件里会被当成上下颠倒的"负分辨率"文件，大多数阅读器直接拒绝。

3.3 一份可以照着敲的转换脚本

下面是我现在常用的一版，用rasterio显式写出，可控性最好：

import numpy as np import rasterio from rasterio.transform import from_origin import xarray as xr SRC = "oco2_L3CO2_day_20150101_v10r.nc" DST = "xco2_20150101.tif" VAR = "xco2" ds = xr.open_dataset(SRC) lat = ds["lat"].values.astype("float64") lon = ds["lon"].values.astype("float64") # 分辨率取中位数，避免浮点抖动导致的边界误差 res_x = float(np.median(np.diff(lon))) res_y = float(np.median(np.diff(lat))) data = ds[VAR].values.astype("float32") # 缺测值转 NaN data = np.where(data > -900000, data, np.nan) # 上下翻转，让第 0 行对齐北边界 data = np.flipud(data) left = lon[0] - res_x / 2.0 top = lat[-1] + res_y / 2.0 transform = from_origin(left, top, res_x, res_y) profile = { "driver": "GTiff", "height": data.shape[0], "width": data.shape[1], "count": 1, "dtype": "float32", "crs": "EPSG:4326", "transform": transform, "nodata": np.nan, "compress": "deflate", "predictor": 3, "tiled": True, "blockxsize": 256, "blockysize": 256, } with rasterio.open(DST, "w", **profile) as dst: dst.write(data, 1) dst.update_tags( UNITS="ppm", LONG_NAME="Monthly/DC XCO2 from gridded product", SOURCE_VAR=VAR, FILL_VALUE="-999999.0", )

几个参数选择背后的理由值得说一下。

dtype选float32而不是float64，是因为 XCO2 的物理意义只精确到小数点后两三位，float32的有效数字完全够用，文件体积直接减半。

compress="deflate"加predictor=3（浮点预测器）是比较稳的组合。全球 0.5 度分辨率的单日文件，不加压缩大概 180×360×4 字节，也就 260 KB；加了压缩通常能到 80 KB 左右。如果是 0.25 度或者月平均产品，压缩带来的收益会更明显。

tiled=True是为了让 QGIS、GeoServer 这类软件局部读取时不用加载整幅图。不过对于这种小尺寸的全球图，分块其实意义有限，如果文件小于 1 MB，我更倾向不分块，减少文件头开销。

nodata=np.nan是 GDAL 支持的特性，写入后 tif 头的 GDAL_NODATA 标签就是nan。有些老旧软件对 nan 的 nodata 处理不好，会把它当有效值参与统计，如果你的下游是这类软件，改成-9999.0更保险。

3.4 批量处理与命名规范

真实场景里你很少只转一个文件。约束条件是：一年的日产品就是 365 个文件，转完还要能按时间排序、能看出变量名和分辨率。

我的命名习惯是<变量>_<产品类型>_<时间>_<分辨率>_<版本>.tif，比如xco2_L3day_20150101_0p5deg_v10.tif。好处是丢进任何文件管理器，按名称排序天然就是时间序列；用脚本做glob匹配也方便。

批量脚本的骨架：

from pathlib import Path files = sorted(Path("./raw").glob("oco2_L3CO2_day_*_v10r.nc")) for f in files: try: date_str = f.stem.split("_")[3] # 从文件名取日期段 out = Path("./tif") / f"xco2_L3day_{date_str}_0p5deg_v10.tif" convert_one(f, out) except Exception as e: print(f"[FAIL] {f.name}: {e}")

注意那个try/except。批处理最怕的就是中途一个坏文件把整个流程打断，等你回来发现只跑了一半。捕获异常、打印文件名、继续下一个，是最省心的做法。

4. 点产品：别硬塞进栅格

4.1 点数据的本质决定了它不能被"直接转换"

L2 Lite 产品的结构长这样：xco2的维度是sounding_id，长度可能是 80000 到 150000 不等（取决于当天的观测覆盖）。每个 sounding 对应一条窄条带内的一个观测足迹，足迹在地面上的直径大概 1.3 公里左右，但相邻足迹之间的间距、以及不同轨道之间的间距，都是高度不均匀的。

这种数据结构有两层含义。

第一层：它天然适合做矢量点，每个点带着自己的坐标和属性，GIS 里可以直接按属性符号化、可以做空间查询。

第二层：如果你想得到一张"全球 XCO2 分布图"，必须自己做空间聚合——把散落的点按规则格网分箱，然后对每个箱内的点求统计量。这个过程在遥感里叫 binning，本质上是降采样。

把这两件事混在一起，就会出现开头说的"噪声图"。我见过的最典型错误是有人写data.reshape(360, 720)，然后抱怨图不对——他把时间序列切成了空间网格。

4.2 路线一：转成矢量点，保真度最高

如果你下游要做的是空间分析、按区域统计、或者给每个点保留全部属性，那矢量点是最合适的选择。

import geopandas as gpd from shapely.geometry import Point import pandas as pd import xarray as xr ds = xr.open_dataset("oco2_L2_Lite_FP_20150101.nc") df = pd.DataFrame({ "xco2": ds["xco2"].values, "qf": ds["xco2_quality_flag"].values, "lat": ds["latitude"].values, "lon": ds["longitude"].values, "time": ds["time"].values, }) # 只保留通过质量筛选的记录 df = df[(df["qf"] == 0) & (df["xco2"] > 0)] gdf = gpd.GeoDataFrame( df, geometry=[Point(xy) for xy in zip(df["lon"], df["lat"])], crs="EPSG:4326", ) gdf.to_file("xco2_points_20150101.gpkg", driver="GPKG", layer="xco2")

选 GeoPackage 而不是 Shapefile，理由很实际：Shapefile 字段名限制 10 个字符、单文件不能超过 2 GB、编码容易出问题，而一天的数据量虽然不到 2 GB，但字段名和编码的坑迟早会踩到。GeoPackage 是单文件、支持长字段名、支持 UTF-8，现在主流 GIS 软件都读得动。

如果确实需要 SHP，就加个driver="ESRI Shapefile"，但记得字段名要提前改短。

把time变量存进去的时候要留意：OCO-2 的时间是自 1993 年 1 月 1 日起的秒数（TAI93 时间标准），不是常见的 Unix 时间戳（1970 年）。直接拿去pd.to_datetime会差出 23 年。正确做法是加一个基准偏移：

base = pd.Timestamp("1993-01-01 00:00:00") df["datetime"] = base + pd.to_timedelta(df["time"], unit="s")

这里还有个小细节：TAI 和 UTC 之间有几十秒的跳秒差。对于日尺度的分布图，几十秒完全不影响；但如果你要把 OCO-2 的时间和其他卫星数据做精确对齐，那就得查一下官方的转换说明，不能想当然。

4.3 路线二：分箱聚合，得到规则的栅格

要出"图"而不是"点图层"，就得做分箱。核心思路是给每个点算一个格网索引，然后按索引聚合。

import numpy as np import pandas as pd RES = 0.5 LON_MIN, LON_MAX = -180.0, 180.0 LAT_MIN, LAT_MAX = -90.0, 90.0 n_lon = int(round((LON_MAX - LON_MIN) / RES)) # 720 n_lat = int(round((LAT_MAX - LAT_MIN) / RES)) # 360 df["lon_idx"] = np.floor((df["lon"] - LON_MIN) / RES).astype("int32") df["lat_idx"] = np.floor((df["lat"] - LAT_MIN) / RES).astype("int32") df = df[(df["lon_idx"] >= 0) & (df["lon_idx"] < n_lon)] df = df[(df["lat_idx"] >= 0) & (df["lat_idx"] < n_lat)] grid = np.full((n_lat, n_lon), np.nan, dtype="float32") grouped = df.groupby(["lat_idx", "lon_idx"])["xco2"].mean() for (i, j), v in grouped.items(): grid[i, j] = v

这段代码里有两个设计选择值得解释。

为什么用np.floor而不是四舍五入？因为格网分箱的本质是把连续空间切成一个个左闭右开的区间。用floor能保证一个点只落入一个箱，不会出现边界上的重复或遗漏。用round的话，在箱边界上的点会向两侧飘，虽然最终结果差异不大，但统计上不严谨。

为什么对每个格网求均值，而不是求和或取中位数？均值是最直接的估计，但要注意：一个 0.5 度格网在赤道附近约 55 公里见方，OCO-2 的单日观测在一个格网里可能只有 1 到 2 个有效点，甚至在轨道间隙处是 0 个点。这时候"均值"其实就是那个单点的值，空间代表性很有限。

我一般的做法是同时输出一个count图层，记录每个格网里参与平均的点数。出图时用count >= 3做掩膜，视觉上会更干净，也避免把孤点当成区域平均值。

另外grid初始化用NaN而不是 0，是为了让空格网在 GeoTIFF 里正确显示为 nodata。如果用 0 初始化，那些没有观测的区域会被渲染成 XCO2 = 0，在地图上是一片刺眼的白或黑，完全误导读者。

最后写出去的时候，行序方向还是那个老问题——lat_idx从 0 开始对应最南边，所以输出前依然要np.flipud。

4.4 空间插值：能做，但要想清楚为什么做

有人会问，格网里那么多空值，能不能插值填上？

技术上当然可以，scipy.interpolate.griddata或者pykrige都能做。但我建议你在插值之前先回答一个问题：你要这张图回答什么问题？

如果是"这片区域的平均碳浓度水平是多少"，那插值会把稀疏观测的局部特征平滑掉，反而丢失信息；更合适的做法是明确标注有效格网，让读者知道哪里是实测、哪里是空白。

如果是"给某个模型提供空间连续的输入场"，那插值是合理的，但必须做交叉验证，报告插值误差，别默认插出来的值可信。

我自己踩过的坑是：某次用反距离加权对单日的 0.5 度格网做填补，结果在轨道条带之间的空隙处插出了一条明显不合理的高值带，原因是加权半径设得太小，插值几乎退化成了最近邻，把边缘的异常点直接复制过去了。后来把搜索半径放大到 3 个格网以上、同时限制最多用 12 个邻居，才勉强能看。

底线是：插值结果必须和原始观测点一起展示，或者至少在图注里说清楚哪些区域是插的。把插值图当实测图发布，是这类数据最容易出的事故。

5. 踩坑实录与速查

5.1 数值全是对不上的：缺测值和缩放因子

最典型的现象是：图出来了，色带范围从 -1000000 到 420，整张图看起来只有一两个亮斑。

原因基本都是缺测值没处理。OCO-2 的填充值是-999999.0，这个数字比真实值小六个数量级，只要有一个漏网，直方图的拉伸就全毁了。

还有一种更隐蔽的情况：文件里某些变量带scale_factor和add_offset属性，实际值 = 存储值 ×scale_factor+add_offset。用xarray.open_dataset默认会自动应用这两个属性（mask_and_scale=True），但用netCDF4直接读变量时不会，你拿到的就是裸整数。两种读法混用，结果能差几倍。

排查动作很简单，转换后立刻做一次统计：

import rasterio import numpy as np with rasterio.open("out.tif") as src: a = src.read(1) print("min", np.nanmin(a), "max", np.nanmax(a), "nan%", np.isnan(a).mean())

XCO2 的合理范围大致是 380 到 420 ppm（近地面），如果统计出来超出这个范围很多，先回头查缺测值和缩放，别怀疑数据本身。

5.2 位置整体偏移或上下颠倒

上一节说的行序问题是最常见的成因。第二种常见成因是经纬度的"中心 vs 边界"混淆——如果你误把lon[0]当成左边界直接填进 GT(0)，整张图会向右偏半个像元，也就是 0.25 度。这个偏移在全球尺度上不显眼，但做区域裁剪或者和别的栅格做叠加分析时，会出现"明显对不上"的现象。

验证方法：拿一个你熟悉的地标坐标，比如某个高原区域，去 tif 里查对应位置的值，再回头查原始 NC 在同一坐标的值，两者应该一致（同一像元内）。

from rasterio.transform import rowcol with rasterio.open("out.tif") as src: row, col = rowcol(src.transform, lon_test, lat_test) print("tif value:", src.read(1)[row, col]) # 对照原始 nc print("nc value:", ds["xco2"].sel(lon=lon_test, lat=lat_test, method="nearest").values)

5.3 QGIS 打开一片黑或者一片白

这通常不是数据的问题，是渲染的问题。

XCO2 的值域很窄，380 到 415 之间变化不到 10%，默认的线性拉伸如果被个别异常值拉宽到 -1e6 到 420，中间那 35 个单位的差异就被压成了一条线，看起来就是纯色。

处理办法有两个方向。一是数据侧：确保缺测值是 nodata 且被正确识别，QGIS 会自动跳过；二是渲染侧：手动设置最小/最大值，或者用 2%–98% 的分位数拉伸。

QGIS 里也可以直接在图层属性里设置"累积计数剪切"，让软件自动忽略两端的极值。

5.4 常见问题速查表

6. 效率优化与下游衔接

6.1 批处理的并行与内存控制

一年 365 个文件，单个文件的处理时间大概零点几秒，串行跑完不到五分钟，其实没必要上并行。但如果你的目标是多年数据、或者更高分辨率的 L3 产品（0.25 度），文件数量和体积都会上去。

用concurrent.futures.ProcessPoolExecutor做并行是最省事的：

from concurrent.futures import ProcessPoolExecutor from pathlib import Path def job(args): src, dst = args convert_one(src, dst) return src.name if __name__ == "__main__": tasks = [(str(f), str(out_path(f))) for f in files] with ProcessPoolExecutor(max_workers=4) as pool: for name in pool.map(job, tasks): print("done", name)

max_workers别设太高。这类任务瓶颈在磁盘 I/O 和xarray打开文件的开销上，核数开到 4 到 8 就够了，再多反而因为争抢 I/O 变慢。

另外提一句xarray.open_dataset的chunks参数。对于这种几十 MB 的文件，加chunks走 dask 反而增加调度开销，得不偿失。真正需要分块的是多时间步的大文件，那是另一个场景。

6.2 输出格式：普通 GeoTIFF 还是 COG

如果 tif 只是本地看看、丢给同事，普通 GeoTIFF 完全够。

如果目标是发布到服务端、做 Web 切片、或者需要支持 HTTP 范围请求的按需读取，那就该转成COG（Cloud Optimized GeoTIFF）。COG 的特点是把金字塔概览层级和分块按特定顺序排列在文件里，客户端可以先读低分辨率概览快速出图，再按需拉取局部细节。

用 GDAL 转很直接：

gdal_translate -of COG \ -co COMPRESS=DEFLATE \ -co OVERVIEWS=IGNORE_EXISTING \ input.tif output_cog.tif

或者用rasterio的driver="COG"，参数基本一致。注意 COG 要求内部必须是分块的，blockxsize一般设 256 或 512，太大太小都不好——太大导致按需读取拉的数据多，太小导致文件头开销大。

6.3 和下游软件的衔接

转出来的 tif 最终要进什么地方，会反过来影响你的输出选择。

进QGIS / ArcGIS：EPSG:4326、单波段 float32、带 nodata，这三个条件满足就能正常打开。如果要多波段叠在一起，记得波段顺序在元数据里写清楚。

进GeoServer / MapServer：最好转 COG，并且确认 nodata 值在服务端配置里也声明了，否则边界会出现黑边。

进Python 分析流程：其实不必转 tif，直接用xarray读 nc 更省事，元数据更全。转换只在这两种情况下有必要——下游工具只认栅格，或者你需要把数据冻结成一个不可变的中间产物。

进前端可视化：通常不是直接给 tif，而是先切片成 XYZ 或者转成栅格瓦片。这时候 nodata 处理和色带映射要在切片阶段就定好，别留给前端。

还有一点，如果最终用途是"多次复用、需要保留全部科学属性"，我强烈建议在 tif 的元数据里把所有关键信息写全：变量名、单位、缺测值、原始文件名、处理时间、质量筛选条件。这些信息在半年后你自己回头看时会救命，因为到那时你大概率已经忘了当时用的是什么阈值。

dst.update_tags( ORIGIN_FILE="oco2_L3CO2_day_20150101_v10r.nc", VARIABLE="xco2", UNITS="ppm", QUALITY_FILTER="xco2_quality_flag == 0", PROCESS_TIME="2025-01-15T10:30:00", RESOLUTION_DEG="0.5", )

我在实际项目中体会最深的一点是：这类数据转换的代码本身不难，难的是把每一个"隐式约定"都显式化。缺测值是哪个数、坐标是中心还是边界、行序朝哪边、时间基准是哪一年——这些在官方文档里往往散落在不同章节，甚至有的只写在产品说明书的小字里。把它们一次性固化进脚本和元数据，后面所有复用都会轻松很多。

另一个反复验证有效的习惯是：每次转完先画一张缩略图看一眼。不用多精细，matplotlib里imshow一下，或者 QGIS 里拖进去扫一眼，几秒钟的事。我靠这个习惯至少拦下过三次"数据全错但文件正常生成"的事故，其中一次是因为把两个不同版本的文件混在了同一批任务里，变量名改了但脚本没报错，只是静默地输出了全空值。看一眼图，问题立刻暴露。