刚拿到一批碳观测卫星的 .nc4 文件很多人第一反应是直接拖进 QGIS 或 ArcGIS——结果要么软件压根不认要么认出来只是一堆散点画不出那种规整的全球浓度分布图。这就是OCO-2 的 NC 格式转 tif这件事的真实起点它不是一个“换个后缀名”的操作而是一次数据结构上的改造——把多维科学数组变成带地理参考的栅格影像。这件事看起来是格式转换实际上卡人的地方全在坐标对齐、缺测值处理、量纲还原和质量筛选上。我前后处理过几万景这类数据踩过的坑够写一本书这里把完整链路和判断逻辑摊开讲适合刚接触遥感数据格式转换的同行也适合已经会写脚本、但总在“影像偏半格”“值全黑”“南北颠倒”上反复返工的人。1. 先看清手里这份 NC 到底是什么产品动手写第一行代码之前有个判断必须先做完你下载的是点数据产品还是格网产品。这两类文件的转换路径完全不同搞错了后面全白干。1.1 点集与栅格两种产品的本质差异这颗卫星的对地观测方式是推扫式的一次过境沿着轨道打出一串脚印footprint每个脚印对应一条独立的观测记录。所谓L2 级别产品本质上就是一张“观测记录清单”每行有 sounding_id、纬度、经度、时间、二氧化碳干空气摩尔分数XCO2、质量标志、不确定度等字段。文件里可能是几万到十几万条记录纬度经度都是逐行变化的没有任何网格规律。而L3 级别产品是由 L2 经过时空聚合生成的规则经纬度格网典型分辨率为 0.5° 或 1°维度一般是(time, lat, lon)。这类文件才是真正意义上的“准栅格”转 tif 属于顺手的事。我见过最常见的误判是拿到 L2 Lite 文件看到文件名里有 nc4、里面也有 latitude/longitude 变量就以为可以直接gdal_translate。结果导出来的影像要么是空的要么是一条极细的斜线——因为 GDAL 把第一个和最后一个坐标值当成影像范围中间所有点都被塞进了错误的像素里。提示判断方法很简单打开文件看一眼纬度的形状。如果lat.shape (N,)且 N 与观测条数一致、相邻值没有固定间隔那就是点数据如果lat是等间距排列的比如从 89.75 到 -89.75 步长 0.5那就是格网数据。1.2 打开文件的第一件事摸清维度和变量名不同版本、不同处理级别的产品变量命名并不统一。有的版本把浓度变量叫xco2有的叫xco2_2020之类带年份后缀的名字有的在 L3 月度产品里会拆成多个变量分别对应不同的质量等级。所以我不建议任何脚本里硬编码变量名先探一遍import xarray as xr ds xr.open_dataset(your_file.nc4, decode_cfFalse, mask_and_scaleFalse) print(ds.dims) for name, var in ds.variables.items(): print(name, var.dims, var.shape, var.dtype, var.attrs.get(units))这里我特意关掉了decode_cf和mask_and_scale原因是自动解码有时候会把_FillValue直接变成 NaN、把整型按 scale_factor 转成浮点表面上看更方便但你后面想核对原始填充值、想复现别人的处理流程时就会对不上。做格式转换这种“底层搬运”的活我倾向于先看清楚原始状态再自己决定怎么处理。1.3 填充值与量纲两个最容易被忽略的字段每个变量几乎都带_FillValue通常是 -999999 或者 9.96921e36 这种天文数字。如果你直接把这个数组写进 tifNoData 会变成影像里的极值拉伸显示的时候整幅图会变成一坨灰因为色带被这个异常值拉爆了。更隐蔽的是打包存储。有些产品为了节省空间把浮点数据存成整型另外用scale_factor和add_offset两个属性记录换算关系真实值 原始值 × scale_factor add_offset。xarray 默认会自动应用pandas 读 csv 版本时却不会。如果你混用了两条路径的结果会发现两批数据的数值差了整整一个数量级还以为是自己算错了。处理原则我总结成一条先还原物理量再处理填充值最后才谈栅格化。顺序反了你可能把填充值也参与了 scale 运算得到一个看起来“合理”但完全错误的数。2. 栅格化之前必须想清楚的坐标问题坐标这一块是整件事里最容易出错、也最难自查的部分。影像看着“差不多对”但偏差半格或者上下颠倒往往要等到和别的底图叠加时才发现。2.1 像元中心与像元边界半格偏移从哪来NC 文件里的lat和lon数组存的是每个格子的中心坐标。而 GeoTIFF 的地理变换geotransform记录的是左上角像元的左上角点的坐标再加上像元宽高。举个例子经度从 -179.75 开始步长 0.5那么第一个像元的中心是 -179.75它的左边界就是 -179.75 - 0.25 -180.0。如果你直接把 -179.75 当成左上角坐标写进去整幅影像就会往东南方向偏半个像元。0.5° 分辨率下这个偏移是 0.25°换算到赤道大约是 27 公里——叠加到行政区划图上肉眼可见地错位。正确做法import numpy as np from rasterio.transform import from_origin lon ds[lon].values lat ds[lat].values dlon abs(lon[1] - lon[0]) dlat abs(lat[1] - lat[0]) lon_left lon[0] - dlon / 2.0 lat_top lat.max() dlat / 2.0 transform from_origin(lon_left, lat_top, dlon, dlat)用 GDAL 命令行时同样要检查gdalinfo输出里的Origin。如果 Origin 的经度正好等于坐标数组的第一个值那基本可以确定没做半格偏移需要用-a_ullr手动把影像范围钉死。2.2 纬度递减一个让所有人翻车的细节很多产品的纬度是从北向南递减排列的90 → -90因为这种顺序符合气象学里“从上到下”的书写习惯。但 GeoTIFF 的存储顺序永远是第一行在最北。两者一撞结果就是南北颠倒。判断和修正的写法flip lat[0] lat[-1] if flip: lat lat[::-1] data data[..., ::-1, :] # 只在纬度维度翻转我见过有人用np.flipud一把梭结果把三维数组的时间维和纬度维一起翻了影像波段顺序全乱。一定按轴索引来翻别图省事。还有一个隐蔽情况少数产品纬度本身是递增的但顶层属性里写了positive: down表示“数值增大方向朝下”。这种元数据与数组实际顺序矛盾的情况以数组实际顺序为准元数据只是给可视化软件看的。2.3 全球网格的边界对齐与覆盖面0.5° 分辨率的全球网格经度方向应该有 720 列-180 到 180纬度方向 360 行-90 到 90。如果你算出来的列数是 719 或者 721说明网格边界没有对齐——多半是产品用了 -179.75 到 179.75 的像元中心共 720 个这没问题但如果它存的是 -180 到 180 共 721 个点那就是把边界点也当成了中心点需要重新推算。这个自检非常值得加进脚本写出 tif 之前打印一下width和height对不上就停下来看。全球 0.5° 是 (360, 720)1° 是 (180, 360)1°×1.25° 是 (180, 288)。这些数字记熟比事后排查省事得多。3. Python 路线从读取到写出影像的完整实现大部分人的诉求是“批量、可控、能塞进自动化流程”Python 是最合适的选择。下面这套写法我用了很久稳定、可迁移。3.1 依赖组合与版本注意点核心依赖就四个xarray负责读numpy负责算rasterio负责写rioxarray作为可选加速层。pip install xarray numpy rasterio netCDF4 rioxarraynetCDF4这个后端必须装否则 xarray 打开 nc4 会退回到 scipy 引擎而 scipy 引擎对 netCDF4 格式HDF5 底层支持不全会直接报错。另外如果文件带有压缩chunkingscipy 引擎也读不了。版本上有个小坑rasterio 2.x 和 3.x 的write接口参数名有差异老脚本迁移时nodata位置参数会报错。新项目直接用 3.x 就行。3.2 L3 格网产品的最小可用脚本import numpy as np import xarray as xr import rasterio from rasterio.transform import from_origin def nc_grid_to_tif(nc_path, tif_path, var_namexco2): ds xr.open_dataset(nc_path, decode_cfFalse, mask_and_scaleFalse) lat ds[lat].values.astype(float64) lon ds[lon].values.astype(float64) data ds[var_name].values.astype(float32) # 1) 还原量纲若存在打包存储 scale ds[var_name].attrs.get(scale_factor, 1.0) offset ds[var_name].attrs.get(add_offset, 0.0) if scale ! 1.0 or offset ! 0.0: data data * scale offset # 2) 填充值转 NaN fill ds[var_name].attrs.get(_FillValue, None) if fill is not None: data np.where(data np.float32(fill), np.nan, data) data np.where(np.abs(data) 1e30, np.nan, data) # 3) 纬度方向修正 if lat[0] lat[-1]: lat lat[::-1] data data[..., ::-1, :] # 4) 地理变换 dlon float(abs(lon[1] - lon[0])) dlat float(abs(lat[1] - lat[0])) transform from_origin(lon[0] - dlon / 2.0, lat[-1] dlat / 2.0, dlon, dlat) # 5) 维度整理把 time 放到第一维其余为 (band, row, col) if data.ndim 2: data data[np.newaxis, :, :] elif data.ndim 3: pass else: raise ValueError(f维度不支持: {data.shape}) profile { driver: GTiff, height: data.shape[1], width: data.shape[2], count: data.shape[0], dtype: float32, crs: EPSG:4326, transform: transform, nodata: np.nan, compress: deflate, zlevel: 6, predictor: 3, tiled: True, blockxsize: 256, blockysize: 256, BIGTIFF: IF_SAFER, } with rasterio.open(tif_path, w, **profile) as dst: dst.write(data) dst.update_tags(1, UNITSppm, SOURCEnc_path.split(/)[-1]) ds.close() nc_grid_to_tif(oco2_l3_day.nc4, oco2_l3_day.tif)这套代码里predictor3是给浮点数据用的配合 DEFLATE 能把文件体积压到原来的三分之一左右如果是整型数据应该用predictor2。BIGTIFFIF_SAFER是保命用的——全球高分辨率多波段文件很容易突破 4GB 的单文件上限不设置会直接写入失败。3.3 质量标志过滤不做这一步等于白忙观测数据里永远混着一部分不可信的记录云层遮挡、地表反射率过低、气溶胶污染、观测几何条件差。产品会提供一个质量标志变量通常0表示“好”非 0 表示各类问题。在栅格化阶段我的处理方式是先生成一个布尔掩膜再把不合格的位置置为 NaN。qf ds[xco2_quality_flag].values mask (qf ! 0) data_masked np.where(mask, np.nan, data)不建议直接删除不合格数据再重排数组因为在格网产品里删除会导致行列对应关系断裂坐标就对不上了。只有点数据L2才适合先删后聚合。还有一种情况是多级质量产品可能提供warn_level0-15 分级0 是最优1-3 是“可用但有告警”。做科研分析时可以根据需求放宽到 ≤1 或 ≤2。这个阈值没有绝对标准取决于你的应用场景——做全球碳通量反演必须严格做教学演示可以放松。我通常会同时输出严格版和宽松版两个文件方便对比。3.4 批量处理的命名与断点续跑一天一个文件的产品一年就是 365 个。批量脚本必须解决两件事命名规范要能让后续按时间排序以及中途失败能续跑。from pathlib import Path src_dir Path(./raw_nc) dst_dir Path(./tif_out) dst_dir.mkdir(exist_okTrue) for nc_file in sorted(src_dir.glob(*.nc4)): out dst_dir / (nc_file.stem _xco2.tif) if out.exists(): continue try: nc_grid_to_tif(str(nc_file), str(out)) print(OK, nc_file.name) except Exception as e: print(FAIL, nc_file.name, repr(e))if out.exists(): continue这一行看起来平平无奇但在处理上千个文件时能救命——凌晨三点脚本因为单个坏文件中断不用重跑前面已经完成的几百个。文件名里的日期我一般单独抽出成子目录2020/01/oco2_20200101.tif一是避免单目录文件过多导致ls卡死二是方便后续用 GDAL 的 VRT 按时间拼接。4. GDAL 命令行不写代码的快速通道有些场景你只是想把几个文件转出来看看写脚本反而麻烦。GDAL 的命令行工具链在这种场景下效率极高但前提是得先摸清它怎么“看”这个文件。4.1 用 gdalinfo 摸清子数据集的真面目netCDF4 文件在 GDAL 眼里是一堆子数据集的集合。直接gdalinfo yourfile.nc4输出里会有一大段SUBDATASET_n_NAMENETCDF:yourfile.nc4:/variable_name的列表。gdalinfo oco2_l3_day.nc4 | grep -i subdataset看完列表后针对性查看某个变量的信息gdalinfo NETCDF:oco2_l3_day.nc4:xco2重点看三处Size is后面的行列数对不对、地理变换里的Origin和Pixel Size是否合理、NoData Value有没有被识别出来。这三项里最容易出问题的是地理变换。GDAL 的 netCDF 驱动会根据坐标变量的units属性比如degrees_north自动推断空间参考但如果产品用的是非标准的单位字符串或者坐标变量存在但维度顺序特殊推断就会失败输出里会出现Coordinate System is:后面空白。遇到这种情况可以用-a_srs和-a_ullr手动指定gdal_translate -of GTiff \ -a_srs EPSG:4326 \ -a_ullr -180 90 180 -90 \ -a_nodata -9999 \ -co COMPRESSDEFLATE -co PREDICTOR3 -co TILEDYES \ -co BIGTIFFIF_SAFER \ NETCDF:oco2_l3_day.nc4:xco2 out.tif-a_ullr的四个参数顺序是左上经度、左上纬度、右下经度、右下纬度很多人会记成左下右上写反了影像会整个翻转。这是命令行路线上最高频的错误。4.2 gdal_translate 与 gdalwarp 的分工gdal_translate只做格式转换和元数据改写不做重采样速度快、不损数据。gdalwarp做重投影和重采样适合需要统一分辨率的场景——比如你要把 0.5° 的产品重采样到 1°或者把不同来源的数据统一到同一个格网上做叠加分析。gdalwarp -t_srs EPSG:4326 \ -tr 1.0 1.0 \ -r average \ -srcnodata -9999 -dstnodata -9999 \ -co COMPRESSDEFLATE -co TILEDYES \ in_05deg.tif out_1deg.tif-r average这个重采样方法在浓度数据上要慎用。XCO2 是浓度量取平均是合理的但如果是通量数据比如净生态系统交换量单位是 µmol/m²/s做面积加权平均才正确直接取算术平均会引入偏差。这一点在做多源数据融合时尤其要注意。4.3 压缩参数与文件体积的现实权衡全球 0.5° 单波段浮点影像未压缩大约 1MB压缩后 300KB 左右随便存。但如果是 0.05° 分辨率3600×7200单波段未压缩就是 100MB 往上几百天的数据累积下来轻松过 TB。我的参数选择逻辑场景压缩方式Predictor分块备注单景快速查看DEFLATE z63关闭简单兼容性最好批量存档DEFLATE z93256×256体积最小写入稍慢频繁切片读取LZW3512×512解压快适合云端对象存储极小文件归档ZSTD z153512×512需要 GDAL 3.x 支持压缩率最高分块TILED这个选项单景小文件不开也没关系但一旦文件超过几十 MB、需要频繁按窗口读取或发布成切片服务不开分块会导致每次读取都要解压整行性能差一个数量级。我吃过这个亏一个 2GB 的多时次文件没开分块做动画渲染时每秒只能读两帧。5. L2 点数据怎么落到格网上如果你的文件是点数据产品前面那套代码完全用不了。这是另一条技术路线核心思路是“先建空网格再把点往里装”。5.1 为什么不能把点直接当栅格有人试过把点数据写成带坐标的点矢量再用gdal_rasterize烧录成栅格。这个思路本身没错但会丢失值——gdal_rasterize默认只看点所在位置的属性如果多个点落在同一个像元只有最后一个生效前面的全被覆盖。正确的做法是做格网统计统计每个格子里有多少个有效观测、它们的均值是多少、标准差是多少。均值出浓度图计数出采样密度图标准差出空间变异图。这三张图配套看才能判断某个区域的浓度异常是真实信号还是采样太少导致的噪声。5.2 用 histogram2d 做双通道统计import numpy as np import pandas as pd import xarray as xr import rasterio from rasterio.transform import from_origin def l2_points_to_grid(nc_path, tif_path, res0.5): ds xr.open_dataset(nc_path, decode_cfFalse, mask_and_scaleFalse) lat ds[latitude].values.astype(float64) lon ds[longitude].values.astype(float64) xco2 ds[xco2].values.astype(float64) qf ds[xco2_quality_flag].values # 有效观测筛选 valid (qf 0) np.isfinite(xco2) (xco2 0) (xco2 1000) lat, lon, xco2 lat[valid], lon[valid], xco2[valid] lat_edges np.arange(-90, 90 res, res) lon_edges np.arange(-180, 180 res, res) sum_grid, _, _ np.histogram2d(lat, lon, bins[lat_edges, lon_edges], weightsxco2) cnt_grid, _, _ np.histogram2d(lat, lon, bins[lat_edges, lon_edges]) with np.errstate(invalidignore, divideignore): mean_grid np.where(cnt_grid 0, sum_grid / cnt_grid, np.nan) # histogram2d 的纬度方向是从南到北GeoTIFF 是从北到南 mean_grid np.flipud(mean_grid) transform from_origin(-180.0, 90.0, res, res) profile { driver: GTiff, height: mean_grid.shape[0], width: mean_grid.shape[1], count: 1, dtype: float32, crs: EPSG:4326, transform: transform, nodata: np.nan, compress: deflate, zlevel: 6, predictor: 3, tiled: True, blockxsize: 256, blockysize: 256, } with rasterio.open(tif_path, w, **profile) as dst: dst.write(mean_grid.astype(float32), 1) l2_points_to_grid(oco2_L2_Lite.nc4, oco2_l2_grid.tif, res0.5)这段代码里有几个细节值得展开说。np.histogram2d的bins参数接受的是边界数组长度要比格子数多 1。np.arange(-90, 90 res, res)生成的是 -90 到 90 共 361 个边界点对应 360 个格子正好。如果写成np.arange(-90, 90, res)就少一格南极点附近的数据会被丢弃。np.flipud不能省。histogram2d第一个维度对应 lat从小到大排列第一行是南纬 90°而 GeoTIFF 第一行必须是北纬 90°。忘了翻转整幅图会南北镜像。np.errstate那段是防止除零告警刷屏。虽然结果会被np.where修正但几万次警告打出来会让日志没法看。5.3 稀疏采样导致的空洞与掩膜策略这类推扫观测有个天然缺陷轨道之间有缝隙加上云层过滤全球 0.5° 网格里的有效采样率通常只有百分之几到十几。直接输出的栅格会像一块破布大片大片是 NoData。处理方式取决于你的目的做全球分布图用 1° 或 2° 的粗网格牺牲空间细节换取覆盖率同时用采样计数做掩膜只显示计数 ≥ 3 的格子。做区域对比限定在特定经纬度范围内用更细的网格接受空洞。做时间聚合把一个月甚至一个季度的点合在一起统计能显著提升覆盖。采样计数的输出值得单独做一份。它是判断某个格子浓度可信度的直接依据——只被 1 个观测覆盖的格子浓度值本身就带着观测误差的量级不能和几十个观测平均出来的值等价对待。6. 转换结果的自检与常见故障排查转出来的文件能不能用不能只凭“能打开”三个字。我养成了一套三步自检的习惯能在几分钟内判断这批数据是否可信。6.1 数值、范围、均值三项硬核对第一步看值域。XCO2 的全球背景值大概在 400 ppm 上下区域差异在几十 ppm 量级。如果你导出的影像最小值是 -999999、最大值是 3.4e38说明填充值没处理干净如果值域是 0 到 1说明 scale_factor 被应用了两次或者根本没用。第二步看影像的行列数。全球 0.5° 必须是 720 列 × 360 行对不上就是网格没对齐。第三步算全图均值。用gdalinfo -stats直接看gdalinfo -stats out.tif | grep -E Minimum|Maximum|Mean|NoData均值落在合理区间、NoData 值显示为nan或者你设定的 -9999基本就过关了。这一步只要三十秒能挡掉八成以上的低级错误。6.2 与原始数据交叉验证更严格的做法是抽几个点回查原始 NC。随机取 10 个有值的像素坐标用它的经纬度去 NC 里找对应的观测比对数值。差值应该为零如果是直接搬运或在一个可解释的范围内如果做过重采样或聚合。这一步在做过任何“优化”之后必须做。比如你为了减小体积把 float64 转成了 float32XCO2 数值在 400 左右float32 的有效位数大约是 7 位精度损失在 0.0001 量级完全可接受。但如果你同时做了压缩和量纲转换误差可能被放大就得仔细核。6.3 报错与故障的排查对照表下面这张表是我这些年攒下来的高频问题清单可以直接当查询手册用。现象可能原因排查与修正gdalinfo显示Coordinate System is:为空坐标变量缺units属性或单位字符串非标准加-a_srs EPSG:4326手动指定影像全黑或全灰填充值未转 NaN色带被极值拉爆检查_FillValue用-a_nodata指定或在脚本里置 NaN影像南北颠倒纬度递减数组未翻转判断lat[0] lat[-1]后 flipud影像整体偏半个像元用像元中心当成了左上角坐标减去半个步长或用-a_ullr精确指定边界导出的波段数不对time 维度被识别为波段或反之检查变量的 dims 顺序必要时手动 transposeCPLE_NotSupportedErrorGDAL 编译时未包含 netCDF/HDF5 驱动gdalinfo --formats确认重装完整版读取时卡死或无响应文件过大且未开分块或网络存储延迟开 TILED或先复制到本地再处理写入失败提示超过 4GB未启用 BigTIFF加-co BIGTIFFIF_SAFERxarray 打开报NetCDF: HDF error文件不完整下载中断或后端引擎不匹配重新下载装netCDF4后端还有一类问题不在表里但很折磨人文件能转、能打开、数值也对就是放进 QGIS 后坐标对不上任何底图。这种情况九成是轴向定义问题——tif 的坐标轴方向和你以为的不一样或者影像被写成了“南上北下”的负 y 分辨率。rasterio 允许生成 dy 为负的 transform部分老版本软件读这种文件会错位。我的习惯是只要源数据纬度递减一律手动翻成递增再写宁可多一次数组复制也不给下游留隐患。我个人在实际操作中的体会是这类工作真正的技术含量不在“转换”两个字上而在转换之前对数据结构的判断和转换之后对结果的质疑。同一份文件用不同的参数组合能转出七八个版本图像看上去都差不多但叠上底图、做个统计差异立刻显形。我现在的习惯是每批数据先花十分钟做三件事打印维度、算一遍值域、随手抽三个像素回查原文件。这十分钟省下来的返工往往是以小时计的。另外一个小技巧是把转换参数分辨率、填充值、压缩方式、质量阈值全部写进日志文件和处理脚本放在一起。半年后回来看某个图为什么是那个样子翻日志比翻记忆靠谱得多。 
长河君 · 2026/10/26 4:38:09 · 阅读 1.2k OCO-2 NC转GeoTIFF：坐标对齐、缺测值与质量筛选实战
刚拿到一批碳观测卫星的 .nc4 文件很多人第一反应是直接拖进 QGIS 或 ArcGIS——结果要么软件压根不认要么认出来只是一堆散点画不出那种规整的全球浓度分布图。这就是OCO-2 的 NC 格式转 tif这件事的真实起点它不是一个“换个后缀名”的操作而是一次数据结构上的改造——把多维科学数组变成带地理参考的栅格影像。这件事看起来是格式转换实际上卡人的地方全在坐标对齐、缺测值处理、量纲还原和质量筛选上。我前后处理过几万景这类数据踩过的坑够写一本书这里把完整链路和判断逻辑摊开讲适合刚接触遥感数据格式转换的同行也适合已经会写脚本、但总在“影像偏半格”“值全黑”“南北颠倒”上反复返工的人。1. 先看清手里这份 NC 到底是什么产品动手写第一行代码之前有个判断必须先做完你下载的是点数据产品还是格网产品。这两类文件的转换路径完全不同搞错了后面全白干。1.1 点集与栅格两种产品的本质差异这颗卫星的对地观测方式是推扫式的一次过境沿着轨道打出一串脚印footprint每个脚印对应一条独立的观测记录。所谓L2 级别产品本质上就是一张“观测记录清单”每行有 sounding_id、纬度、经度、时间、二氧化碳干空气摩尔分数XCO2、质量标志、不确定度等字段。文件里可能是几万到十几万条记录纬度经度都是逐行变化的没有任何网格规律。而L3 级别产品是由 L2 经过时空聚合生成的规则经纬度格网典型分辨率为 0.5° 或 1°维度一般是(time, lat, lon)。这类文件才是真正意义上的“准栅格”转 tif 属于顺手的事。我见过最常见的误判是拿到 L2 Lite 文件看到文件名里有 nc4、里面也有 latitude/longitude 变量就以为可以直接gdal_translate。结果导出来的影像要么是空的要么是一条极细的斜线——因为 GDAL 把第一个和最后一个坐标值当成影像范围中间所有点都被塞进了错误的像素里。提示判断方法很简单打开文件看一眼纬度的形状。如果lat.shape (N,)且 N 与观测条数一致、相邻值没有固定间隔那就是点数据如果lat是等间距排列的比如从 89.75 到 -89.75 步长 0.5那就是格网数据。1.2 打开文件的第一件事摸清维度和变量名不同版本、不同处理级别的产品变量命名并不统一。有的版本把浓度变量叫xco2有的叫xco2_2020之类带年份后缀的名字有的在 L3 月度产品里会拆成多个变量分别对应不同的质量等级。所以我不建议任何脚本里硬编码变量名先探一遍import xarray as xr ds xr.open_dataset(your_file.nc4, decode_cfFalse, mask_and_scaleFalse) print(ds.dims) for name, var in ds.variables.items(): print(name, var.dims, var.shape, var.dtype, var.attrs.get(units))这里我特意关掉了decode_cf和mask_and_scale原因是自动解码有时候会把_FillValue直接变成 NaN、把整型按 scale_factor 转成浮点表面上看更方便但你后面想核对原始填充值、想复现别人的处理流程时就会对不上。做格式转换这种“底层搬运”的活我倾向于先看清楚原始状态再自己决定怎么处理。1.3 填充值与量纲两个最容易被忽略的字段每个变量几乎都带_FillValue通常是 -999999 或者 9.96921e36 这种天文数字。如果你直接把这个数组写进 tifNoData 会变成影像里的极值拉伸显示的时候整幅图会变成一坨灰因为色带被这个异常值拉爆了。更隐蔽的是打包存储。有些产品为了节省空间把浮点数据存成整型另外用scale_factor和add_offset两个属性记录换算关系真实值 原始值 × scale_factor add_offset。xarray 默认会自动应用pandas 读 csv 版本时却不会。如果你混用了两条路径的结果会发现两批数据的数值差了整整一个数量级还以为是自己算错了。处理原则我总结成一条先还原物理量再处理填充值最后才谈栅格化。顺序反了你可能把填充值也参与了 scale 运算得到一个看起来“合理”但完全错误的数。2. 栅格化之前必须想清楚的坐标问题坐标这一块是整件事里最容易出错、也最难自查的部分。影像看着“差不多对”但偏差半格或者上下颠倒往往要等到和别的底图叠加时才发现。2.1 像元中心与像元边界半格偏移从哪来NC 文件里的lat和lon数组存的是每个格子的中心坐标。而 GeoTIFF 的地理变换geotransform记录的是左上角像元的左上角点的坐标再加上像元宽高。举个例子经度从 -179.75 开始步长 0.5那么第一个像元的中心是 -179.75它的左边界就是 -179.75 - 0.25 -180.0。如果你直接把 -179.75 当成左上角坐标写进去整幅影像就会往东南方向偏半个像元。0.5° 分辨率下这个偏移是 0.25°换算到赤道大约是 27 公里——叠加到行政区划图上肉眼可见地错位。正确做法import numpy as np from rasterio.transform import from_origin lon ds[lon].values lat ds[lat].values dlon abs(lon[1] - lon[0]) dlat abs(lat[1] - lat[0]) lon_left lon[0] - dlon / 2.0 lat_top lat.max() dlat / 2.0 transform from_origin(lon_left, lat_top, dlon, dlat)用 GDAL 命令行时同样要检查gdalinfo输出里的Origin。如果 Origin 的经度正好等于坐标数组的第一个值那基本可以确定没做半格偏移需要用-a_ullr手动把影像范围钉死。2.2 纬度递减一个让所有人翻车的细节很多产品的纬度是从北向南递减排列的90 → -90因为这种顺序符合气象学里“从上到下”的书写习惯。但 GeoTIFF 的存储顺序永远是第一行在最北。两者一撞结果就是南北颠倒。判断和修正的写法flip lat[0] lat[-1] if flip: lat lat[::-1] data data[..., ::-1, :] # 只在纬度维度翻转我见过有人用np.flipud一把梭结果把三维数组的时间维和纬度维一起翻了影像波段顺序全乱。一定按轴索引来翻别图省事。还有一个隐蔽情况少数产品纬度本身是递增的但顶层属性里写了positive: down表示“数值增大方向朝下”。这种元数据与数组实际顺序矛盾的情况以数组实际顺序为准元数据只是给可视化软件看的。2.3 全球网格的边界对齐与覆盖面0.5° 分辨率的全球网格经度方向应该有 720 列-180 到 180纬度方向 360 行-90 到 90。如果你算出来的列数是 719 或者 721说明网格边界没有对齐——多半是产品用了 -179.75 到 179.75 的像元中心共 720 个这没问题但如果它存的是 -180 到 180 共 721 个点那就是把边界点也当成了中心点需要重新推算。这个自检非常值得加进脚本写出 tif 之前打印一下width和height对不上就停下来看。全球 0.5° 是 (360, 720)1° 是 (180, 360)1°×1.25° 是 (180, 288)。这些数字记熟比事后排查省事得多。3. Python 路线从读取到写出影像的完整实现大部分人的诉求是“批量、可控、能塞进自动化流程”Python 是最合适的选择。下面这套写法我用了很久稳定、可迁移。3.1 依赖组合与版本注意点核心依赖就四个xarray负责读numpy负责算rasterio负责写rioxarray作为可选加速层。pip install xarray numpy rasterio netCDF4 rioxarraynetCDF4这个后端必须装否则 xarray 打开 nc4 会退回到 scipy 引擎而 scipy 引擎对 netCDF4 格式HDF5 底层支持不全会直接报错。另外如果文件带有压缩chunkingscipy 引擎也读不了。版本上有个小坑rasterio 2.x 和 3.x 的write接口参数名有差异老脚本迁移时nodata位置参数会报错。新项目直接用 3.x 就行。3.2 L3 格网产品的最小可用脚本import numpy as np import xarray as xr import rasterio from rasterio.transform import from_origin def nc_grid_to_tif(nc_path, tif_path, var_namexco2): ds xr.open_dataset(nc_path, decode_cfFalse, mask_and_scaleFalse) lat ds[lat].values.astype(float64) lon ds[lon].values.astype(float64) data ds[var_name].values.astype(float32) # 1) 还原量纲若存在打包存储 scale ds[var_name].attrs.get(scale_factor, 1.0) offset ds[var_name].attrs.get(add_offset, 0.0) if scale ! 1.0 or offset ! 0.0: data data * scale offset # 2) 填充值转 NaN fill ds[var_name].attrs.get(_FillValue, None) if fill is not None: data np.where(data np.float32(fill), np.nan, data) data np.where(np.abs(data) 1e30, np.nan, data) # 3) 纬度方向修正 if lat[0] lat[-1]: lat lat[::-1] data data[..., ::-1, :] # 4) 地理变换 dlon float(abs(lon[1] - lon[0])) dlat float(abs(lat[1] - lat[0])) transform from_origin(lon[0] - dlon / 2.0, lat[-1] dlat / 2.0, dlon, dlat) # 5) 维度整理把 time 放到第一维其余为 (band, row, col) if data.ndim 2: data data[np.newaxis, :, :] elif data.ndim 3: pass else: raise ValueError(f维度不支持: {data.shape}) profile { driver: GTiff, height: data.shape[1], width: data.shape[2], count: data.shape[0], dtype: float32, crs: EPSG:4326, transform: transform, nodata: np.nan, compress: deflate, zlevel: 6, predictor: 3, tiled: True, blockxsize: 256, blockysize: 256, BIGTIFF: IF_SAFER, } with rasterio.open(tif_path, w, **profile) as dst: dst.write(data) dst.update_tags(1, UNITSppm, SOURCEnc_path.split(/)[-1]) ds.close() nc_grid_to_tif(oco2_l3_day.nc4, oco2_l3_day.tif)这套代码里predictor3是给浮点数据用的配合 DEFLATE 能把文件体积压到原来的三分之一左右如果是整型数据应该用predictor2。BIGTIFFIF_SAFER是保命用的——全球高分辨率多波段文件很容易突破 4GB 的单文件上限不设置会直接写入失败。3.3 质量标志过滤不做这一步等于白忙观测数据里永远混着一部分不可信的记录云层遮挡、地表反射率过低、气溶胶污染、观测几何条件差。产品会提供一个质量标志变量通常0表示“好”非 0 表示各类问题。在栅格化阶段我的处理方式是先生成一个布尔掩膜再把不合格的位置置为 NaN。qf ds[xco2_quality_flag].values mask (qf ! 0) data_masked np.where(mask, np.nan, data)不建议直接删除不合格数据再重排数组因为在格网产品里删除会导致行列对应关系断裂坐标就对不上了。只有点数据L2才适合先删后聚合。还有一种情况是多级质量产品可能提供warn_level0-15 分级0 是最优1-3 是“可用但有告警”。做科研分析时可以根据需求放宽到 ≤1 或 ≤2。这个阈值没有绝对标准取决于你的应用场景——做全球碳通量反演必须严格做教学演示可以放松。我通常会同时输出严格版和宽松版两个文件方便对比。3.4 批量处理的命名与断点续跑一天一个文件的产品一年就是 365 个。批量脚本必须解决两件事命名规范要能让后续按时间排序以及中途失败能续跑。from pathlib import Path src_dir Path(./raw_nc) dst_dir Path(./tif_out) dst_dir.mkdir(exist_okTrue) for nc_file in sorted(src_dir.glob(*.nc4)): out dst_dir / (nc_file.stem _xco2.tif) if out.exists(): continue try: nc_grid_to_tif(str(nc_file), str(out)) print(OK, nc_file.name) except Exception as e: print(FAIL, nc_file.name, repr(e))if out.exists(): continue这一行看起来平平无奇但在处理上千个文件时能救命——凌晨三点脚本因为单个坏文件中断不用重跑前面已经完成的几百个。文件名里的日期我一般单独抽出成子目录2020/01/oco2_20200101.tif一是避免单目录文件过多导致ls卡死二是方便后续用 GDAL 的 VRT 按时间拼接。4. GDAL 命令行不写代码的快速通道有些场景你只是想把几个文件转出来看看写脚本反而麻烦。GDAL 的命令行工具链在这种场景下效率极高但前提是得先摸清它怎么“看”这个文件。4.1 用 gdalinfo 摸清子数据集的真面目netCDF4 文件在 GDAL 眼里是一堆子数据集的集合。直接gdalinfo yourfile.nc4输出里会有一大段SUBDATASET_n_NAMENETCDF:yourfile.nc4:/variable_name的列表。gdalinfo oco2_l3_day.nc4 | grep -i subdataset看完列表后针对性查看某个变量的信息gdalinfo NETCDF:oco2_l3_day.nc4:xco2重点看三处Size is后面的行列数对不对、地理变换里的Origin和Pixel Size是否合理、NoData Value有没有被识别出来。这三项里最容易出问题的是地理变换。GDAL 的 netCDF 驱动会根据坐标变量的units属性比如degrees_north自动推断空间参考但如果产品用的是非标准的单位字符串或者坐标变量存在但维度顺序特殊推断就会失败输出里会出现Coordinate System is:后面空白。遇到这种情况可以用-a_srs和-a_ullr手动指定gdal_translate -of GTiff \ -a_srs EPSG:4326 \ -a_ullr -180 90 180 -90 \ -a_nodata -9999 \ -co COMPRESSDEFLATE -co PREDICTOR3 -co TILEDYES \ -co BIGTIFFIF_SAFER \ NETCDF:oco2_l3_day.nc4:xco2 out.tif-a_ullr的四个参数顺序是左上经度、左上纬度、右下经度、右下纬度很多人会记成左下右上写反了影像会整个翻转。这是命令行路线上最高频的错误。4.2 gdal_translate 与 gdalwarp 的分工gdal_translate只做格式转换和元数据改写不做重采样速度快、不损数据。gdalwarp做重投影和重采样适合需要统一分辨率的场景——比如你要把 0.5° 的产品重采样到 1°或者把不同来源的数据统一到同一个格网上做叠加分析。gdalwarp -t_srs EPSG:4326 \ -tr 1.0 1.0 \ -r average \ -srcnodata -9999 -dstnodata -9999 \ -co COMPRESSDEFLATE -co TILEDYES \ in_05deg.tif out_1deg.tif-r average这个重采样方法在浓度数据上要慎用。XCO2 是浓度量取平均是合理的但如果是通量数据比如净生态系统交换量单位是 µmol/m²/s做面积加权平均才正确直接取算术平均会引入偏差。这一点在做多源数据融合时尤其要注意。4.3 压缩参数与文件体积的现实权衡全球 0.5° 单波段浮点影像未压缩大约 1MB压缩后 300KB 左右随便存。但如果是 0.05° 分辨率3600×7200单波段未压缩就是 100MB 往上几百天的数据累积下来轻松过 TB。我的参数选择逻辑场景压缩方式Predictor分块备注单景快速查看DEFLATE z63关闭简单兼容性最好批量存档DEFLATE z93256×256体积最小写入稍慢频繁切片读取LZW3512×512解压快适合云端对象存储极小文件归档ZSTD z153512×512需要 GDAL 3.x 支持压缩率最高分块TILED这个选项单景小文件不开也没关系但一旦文件超过几十 MB、需要频繁按窗口读取或发布成切片服务不开分块会导致每次读取都要解压整行性能差一个数量级。我吃过这个亏一个 2GB 的多时次文件没开分块做动画渲染时每秒只能读两帧。5. L2 点数据怎么落到格网上如果你的文件是点数据产品前面那套代码完全用不了。这是另一条技术路线核心思路是“先建空网格再把点往里装”。5.1 为什么不能把点直接当栅格有人试过把点数据写成带坐标的点矢量再用gdal_rasterize烧录成栅格。这个思路本身没错但会丢失值——gdal_rasterize默认只看点所在位置的属性如果多个点落在同一个像元只有最后一个生效前面的全被覆盖。正确的做法是做格网统计统计每个格子里有多少个有效观测、它们的均值是多少、标准差是多少。均值出浓度图计数出采样密度图标准差出空间变异图。这三张图配套看才能判断某个区域的浓度异常是真实信号还是采样太少导致的噪声。5.2 用 histogram2d 做双通道统计import numpy as np import pandas as pd import xarray as xr import rasterio from rasterio.transform import from_origin def l2_points_to_grid(nc_path, tif_path, res0.5): ds xr.open_dataset(nc_path, decode_cfFalse, mask_and_scaleFalse) lat ds[latitude].values.astype(float64) lon ds[longitude].values.astype(float64) xco2 ds[xco2].values.astype(float64) qf ds[xco2_quality_flag].values # 有效观测筛选 valid (qf 0) np.isfinite(xco2) (xco2 0) (xco2 1000) lat, lon, xco2 lat[valid], lon[valid], xco2[valid] lat_edges np.arange(-90, 90 res, res) lon_edges np.arange(-180, 180 res, res) sum_grid, _, _ np.histogram2d(lat, lon, bins[lat_edges, lon_edges], weightsxco2) cnt_grid, _, _ np.histogram2d(lat, lon, bins[lat_edges, lon_edges]) with np.errstate(invalidignore, divideignore): mean_grid np.where(cnt_grid 0, sum_grid / cnt_grid, np.nan) # histogram2d 的纬度方向是从南到北GeoTIFF 是从北到南 mean_grid np.flipud(mean_grid) transform from_origin(-180.0, 90.0, res, res) profile { driver: GTiff, height: mean_grid.shape[0], width: mean_grid.shape[1], count: 1, dtype: float32, crs: EPSG:4326, transform: transform, nodata: np.nan, compress: deflate, zlevel: 6, predictor: 3, tiled: True, blockxsize: 256, blockysize: 256, } with rasterio.open(tif_path, w, **profile) as dst: dst.write(mean_grid.astype(float32), 1) l2_points_to_grid(oco2_L2_Lite.nc4, oco2_l2_grid.tif, res0.5)这段代码里有几个细节值得展开说。np.histogram2d的bins参数接受的是边界数组长度要比格子数多 1。np.arange(-90, 90 res, res)生成的是 -90 到 90 共 361 个边界点对应 360 个格子正好。如果写成np.arange(-90, 90, res)就少一格南极点附近的数据会被丢弃。np.flipud不能省。histogram2d第一个维度对应 lat从小到大排列第一行是南纬 90°而 GeoTIFF 第一行必须是北纬 90°。忘了翻转整幅图会南北镜像。np.errstate那段是防止除零告警刷屏。虽然结果会被np.where修正但几万次警告打出来会让日志没法看。5.3 稀疏采样导致的空洞与掩膜策略这类推扫观测有个天然缺陷轨道之间有缝隙加上云层过滤全球 0.5° 网格里的有效采样率通常只有百分之几到十几。直接输出的栅格会像一块破布大片大片是 NoData。处理方式取决于你的目的做全球分布图用 1° 或 2° 的粗网格牺牲空间细节换取覆盖率同时用采样计数做掩膜只显示计数 ≥ 3 的格子。做区域对比限定在特定经纬度范围内用更细的网格接受空洞。做时间聚合把一个月甚至一个季度的点合在一起统计能显著提升覆盖。采样计数的输出值得单独做一份。它是判断某个格子浓度可信度的直接依据——只被 1 个观测覆盖的格子浓度值本身就带着观测误差的量级不能和几十个观测平均出来的值等价对待。6. 转换结果的自检与常见故障排查转出来的文件能不能用不能只凭“能打开”三个字。我养成了一套三步自检的习惯能在几分钟内判断这批数据是否可信。6.1 数值、范围、均值三项硬核对第一步看值域。XCO2 的全球背景值大概在 400 ppm 上下区域差异在几十 ppm 量级。如果你导出的影像最小值是 -999999、最大值是 3.4e38说明填充值没处理干净如果值域是 0 到 1说明 scale_factor 被应用了两次或者根本没用。第二步看影像的行列数。全球 0.5° 必须是 720 列 × 360 行对不上就是网格没对齐。第三步算全图均值。用gdalinfo -stats直接看gdalinfo -stats out.tif | grep -E Minimum|Maximum|Mean|NoData均值落在合理区间、NoData 值显示为nan或者你设定的 -9999基本就过关了。这一步只要三十秒能挡掉八成以上的低级错误。6.2 与原始数据交叉验证更严格的做法是抽几个点回查原始 NC。随机取 10 个有值的像素坐标用它的经纬度去 NC 里找对应的观测比对数值。差值应该为零如果是直接搬运或在一个可解释的范围内如果做过重采样或聚合。这一步在做过任何“优化”之后必须做。比如你为了减小体积把 float64 转成了 float32XCO2 数值在 400 左右float32 的有效位数大约是 7 位精度损失在 0.0001 量级完全可接受。但如果你同时做了压缩和量纲转换误差可能被放大就得仔细核。6.3 报错与故障的排查对照表下面这张表是我这些年攒下来的高频问题清单可以直接当查询手册用。现象可能原因排查与修正gdalinfo显示Coordinate System is:为空坐标变量缺units属性或单位字符串非标准加-a_srs EPSG:4326手动指定影像全黑或全灰填充值未转 NaN色带被极值拉爆检查_FillValue用-a_nodata指定或在脚本里置 NaN影像南北颠倒纬度递减数组未翻转判断lat[0] lat[-1]后 flipud影像整体偏半个像元用像元中心当成了左上角坐标减去半个步长或用-a_ullr精确指定边界导出的波段数不对time 维度被识别为波段或反之检查变量的 dims 顺序必要时手动 transposeCPLE_NotSupportedErrorGDAL 编译时未包含 netCDF/HDF5 驱动gdalinfo --formats确认重装完整版读取时卡死或无响应文件过大且未开分块或网络存储延迟开 TILED或先复制到本地再处理写入失败提示超过 4GB未启用 BigTIFF加-co BIGTIFFIF_SAFERxarray 打开报NetCDF: HDF error文件不完整下载中断或后端引擎不匹配重新下载装netCDF4后端还有一类问题不在表里但很折磨人文件能转、能打开、数值也对就是放进 QGIS 后坐标对不上任何底图。这种情况九成是轴向定义问题——tif 的坐标轴方向和你以为的不一样或者影像被写成了“南上北下”的负 y 分辨率。rasterio 允许生成 dy 为负的 transform部分老版本软件读这种文件会错位。我的习惯是只要源数据纬度递减一律手动翻成递增再写宁可多一次数组复制也不给下游留隐患。我个人在实际操作中的体会是这类工作真正的技术含量不在“转换”两个字上而在转换之前对数据结构的判断和转换之后对结果的质疑。同一份文件用不同的参数组合能转出七八个版本图像看上去都差不多但叠上底图、做个统计差异立刻显形。我现在的习惯是每批数据先花十分钟做三件事打印维度、算一遍值域、随手抽三个像素回查原文件。这十分钟省下来的返工往往是以小时计的。另外一个小技巧是把转换参数分辨率、填充值、压缩方式、质量阈值全部写进日志文件和处理脚本放在一起。半年后回来看某个图为什么是那个样子翻日志比翻记忆靠谱得多。