MATLAB与遥感影像处理：geotiffwrite 写入带地理参考的影像时反复报 proj2gtif 第 17 行错误，SpatialRef 和 R 到底怎么构造才能让 R.RasterSize 与数据矩阵维度对得上？

一、问题的真实位置：报错行号不等于故障点

在 MATLAB 中调用 geotiffwrite 写入带地理参考的影像时，用户常看到类似如下的错误栈：

Error using proj2gtif (line 17)
Error using geotiffwrite>writeGeoTIFF (line ...)
Error using geotiffwrite (line ...)

很多人的第一反应是"投影字符串写错了"或"PROJ 库版本不兼容"。但把 proj2gtif 第 17 行单独拎出来看，它做的是把 MATLAB 的投影结构体翻译成 GeoTIFF 的 GeoKey 目录。这一步失败，往往不是投影本身非法，而是上游传入的元数据自相矛盾，导致翻译器拿不到可用的栅格几何描述。

本文评述：行号只是异常抛出的位置，不是因果链的起点。把 proj2gtif:17 当作"投影库 bug"去排查，方向就偏了。真正需要检查的是 geotiffwrite 在进入投影翻译之前，对 R 与 SpatialRef 做的一致性检查是否被绕过。

从工程视角看，这个报错可以拆成三个可独立验证的命题：

• 尺寸命题：R.RasterSize 是否严格等于 size(A) 的前两维？
• 方向命题：行方向对应北→南还是南→北？列方向对应西→东还是东→西？
• 语义命题：像素代表点（PixelIsPoint）还是面（PixelIsArea）？这决定半像素偏移。

三者任一不满足，投影翻译阶段就可能因为"没有合法的栅格尺寸/方向描述"而中断。下文将围绕这条主线逐层展开。

二、背景：MATLAB 地理参考对象模型与 GeoTIFF 标签体系

2.1 三个核心对象

MATLAB 的 Mapping Toolbox 用三个对象描述一幅地理参考栅格：

需要强调的是，georefcells 与 georefpostings 生成的是两类不同语义的参考对象：前者对应 PixelIsArea（像元是面），后者对应 PixelIsPoint（像元是采样点）。这一区别在写入 GeoTIFF 时会直接映射到 GeoKey 中的 GTRasterTypeGeoKey。

2.2 GeoTIFF 的标签映射

GeoTIFF 规范（OGC 19-008r4，即 GeoTIFF 1.1）把地理参考信息编码进一组 GeoKey，核心包括：

• ModelPixelScaleTag (33550)：像元在 X/Y/Z 方向的地面尺寸；
• ModelTiepointTag (33922)：栅格坐标 (I,J,K) 与世界坐标 (X,Y,Z) 的锚点；
• ModelTransformationTag (34264)：完整的 4×4 仿射矩阵（与上两者互斥）；
• GeoKeyDirectoryTag (34735)：投影、基准面、单位、栅格类型等键值。

MATLAB 的 geotiffwrite 内部会依据 R 生成 ModelPixelScale 与 ModelTiepoint，依据 SpatialRef 生成 GeoKeyDirectory。若 R 的尺寸与数据矩阵不符，ModelPixelScale 的推导会得到零或负值，进而在 proj2gtif 组装 GeoKey 时触发异常。

笔者认为：把 GeoTIFF 写入理解为"元数据序列化"是常见误解。它本质上是一次几何一致性校验 + 两次编码：先校验 R 与矩阵的几何关系，再分别编码像素尺度与投影键值。proj2gtif 处在第二次编码，但它依赖第一次校验的产物。

三、三元一致性：SpatialRef、R、数据矩阵的隐式契约

3.1 尺寸契约

最基础也最常被破坏的契约是：

R.RasterSize == size(A, [1 2])

注意 size(A,[1 2]) 而非 size(A)。多波段影像常见 [H W B] 三维数组，若直接拿 size(A) 去构造 R，会把波段数误当第三空间维，导致 RasterSize 变成 [H W B]，与 GeoTIFF 二维栅格模型冲突。

更隐蔽的是单波段被写成 [H W 1]，此时 size(A,[1 2]) 仍为 [H W]，但若用户用 georasterref('RasterSize', size(A)) 构造，就会得到三维尺寸而报错。

3.2 范围契约

给定 RasterSize = [H W] 与 XWorldLimits = [xmin xmax]、YWorldLimits = [ymin ymax]，分辨率由下式唯一确定：

dx = (xmax - xmin) / W
dy = (ymax - ymin) / H

若用户同时指定 CellExtentInWorldX 与 XWorldLimits，两者必须满足上式，否则 MATLAB 会在构造阶段报"inconsistent"错误；但若用户绕过构造函数、直接改字段，则错误会推迟到写入阶段，表现为 proj2gtif 异常。

3.3 方向契约

MATLAB 图像约定是"第 1 行在顶部"，而地理坐标中"顶部"通常是北。因此对北向上的影像，YWorldLimits 应满足 ymax > ymin 且行索引递增对应纬度递减。若数据本身是南向上（例如某些 HDF/NetCDF 产品），必须先 flipud 再构造 R，否则写出的 GeoTIFF 会出现南北翻转。

本文评述：方向问题不会总是抛异常——它可能"成功写入但内容翻转"，这比报错更危险。proj2gtif 报错反而是幸运情形，因为它阻止了错误数据的落盘。工程上应把方向校验做成写入前的强制断言，而不是依赖异常兜底。

四、构造规则：从地理范围反推 R 与 SpatialRef 的标准流程

4.1 推荐构造顺序

笔者的实践是"先定尺寸、再定范围、后定投影"，避免字段互相覆盖：

% 1) 明确数据矩阵的空间维
[H, W, ~] = size(A);

% 2) 用地理范围 + 尺寸构造 R（PixelIsArea 语义）
latlim = [ymin ymax];      % 注意 ymin < ymax
lonlim = [xmin xmax];
R = georefcells(latlim, lonlim, [H W]);

% 3) 断言尺寸契约
assert(isequal(R.RasterSize, [H W]), 'RasterSize mismatch');

% 4) 构造 SpatialRef（与 R 同语义）
SpatialRef = refmatToGeoRasterReference(R, [H W]);
% 或直接复用 R（georefcells 返回的已是地理参考对象）

% 5) 写入
geotiffwrite('out.tif', A, R, 'GeoKeyDirectoryTag', ...);

关键点在于第 2 步：georefcells(latlim, lonlim, [H W]) 会一次性把尺寸、范围、分辨率三者绑定，避免手工赋值时的不一致。若使用 georasterref 逐字段赋值，务必在最后调用 validate 或手动断言。

4.2 从投影坐标构造

若数据处于投影坐标系（如 UTM），应使用 maprefcells：

xlim = [xmin xmax];   % 投影坐标，单位米
ylim = [ymin ymax];
R = maprefcells(xlim, ylim, [H W], 'ColumnsStartFrom', 'north');

% 投影信息通过 info 结构体传入
info = geotiffinfo('template.tif');   % 复用模板投影
geotiffwrite('out.tif', A, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);

这里 'ColumnsStartFrom','north' 明确声明列索引从北开始，与 MATLAB 图像约定一致。若数据源是 'south'，必须同步翻转矩阵。

4.3 复用已有影像的参考

最常见的稳健做法是从一幅已知正确的 GeoTIFF 读取参考，再套用到新矩阵：

[~, Rref] = geotiffread('template.tif');
info = geotiffinfo('template.tif');

% 若新矩阵尺寸与模板一致，直接复用
assert(isequal(Rref.RasterSize, size(A,[1 2])));
geotiffwrite('out.tif', A, Rref, ...
    'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);

这条路径把"构造"换成"复制"，从源头消除了尺寸不一致的可能。代价是要求新矩阵与模板同尺寸；若不同，需按比例重采样或重新构造 R。

五、方向与语义：PixelIsPoint / PixelIsArea 与轴向陷阱

5.1 半像素偏移

PixelIsArea 把像元视为覆盖 [i, i+1)×[j, j+1) 的面，锚点在角点；PixelIsPoint 把像元视为采样点，锚点在中心。两者在 ModelTiepoint 上相差半个像元：

% PixelIsArea: 左上角点
X_area = xmin;
Y_area = ymax;

% PixelIsPoint: 左上像元中心
X_point = xmin + dx/2;
Y_point = ymax - dy/2;

若用 georefcells（Area）构造，却在 GeoKey 中写入 PixelIsPoint，或反之，就会产生半像素的系统性偏移。这种偏移在低分辨率影像上可忽略，在 10 m 级影像上则是 5 m 的定位误差。

5.2 轴向与 flipud/fliplr

不同数据源的轴向约定并不统一：

MATLAB 的 geotiffread 默认返回北向上、西向东的矩阵，因此从 GeoTIFF 读入的数据无需翻转。但从 HDF 读入时，若直接套用 georefcells 而不翻转，写出的影像会南北颠倒。

笔者认为：方向校验应做成"数据指纹"式的断言——在写入前随机抽取若干行列，比对已知地标的像素值。仅靠元数据自洽不足以发现翻转，因为翻转后的元数据仍然自洽。

六、错误对照表与最小复现用例

6.1 常见错误对照

6.2 最小复现用例

下面这段代码可稳定复现"尺寸不一致"导致的 proj2gtif 异常（模拟数据）：

% 模拟一幅 512x512 单波段影像
A = uint8(rand(512, 512) * 255);

% 错误构造：用 size(A) 但 A 被扩成三维
A3 = reshape(A, [512 512 1]);
R = georefcells([30 31], [104 105], size(A3));  % 得到 [512 512 1]

% 写入时触发 proj2gtif 异常
try
    geotiffwrite('bad.tif', A3, R);
catch ME
    disp(ME.message);   % 观察错误信息
end

% 正确构造
R2 = georefcells([30 31], [104 105], size(A3, [1 2]));
geotiffwrite('good.tif', A3, R2);

该用例说明：错误并非来自投影字符串，而是 RasterSize 被污染为三维。修复只需在构造时显式取前两维。

七、工程化封装：写入前校验、异常捕获与批量流水线

7.1 写入前校验函数

function validateGeoWrite(A, R)
    % 尺寸契约
    assert(isequal(R.RasterSize, size(A, [1 2])), ...
        'RasterSize %s != data %s', mat2str(R.RasterSize), mat2str(size(A,[1 2])));
    % 范围契约
    assert(R.XWorldLimits(1) < R.XWorldLimits(2), 'X limits reversed');
    assert(R.YWorldLimits(1) < R.YWorldLimits(2), 'Y limits reversed');
    % 分辨率非零
    assert(all(R.CellExtentInWorldX > 0) && all(R.CellExtentInWorldY > 0), ...
        'Non-positive cell extent');
    % 语义一致性
    assert(isa(R, 'map.rasterref.GeographicCellsReference') || ...
           isa(R, 'map.rasterref.MapCellsReference'), 'Unexpected ref type');
end

把上述断言放在 geotiffwrite 之前，可把 90% 以上的 proj2gtif 异常拦截在写入之前，错误信息也更可读。

7.2 批量流水线

处理成百上千景影像时，建议把"读取—校验—写入"封装成统一入口，并对每景记录日志：

files = dir('*.tif');
for k = 1:numel(files)
    try
        [A, R] = geotiffread(files(k).name);
        validateGeoWrite(A, R);
        geotiffwrite(['out_' files(k).name], A, R);
        fprintf('[OK] %s\n', files(k).name);
    catch ME
        fprintf('[FAIL] %s : %s\n', files(k).name, ME.message);
    end
end

这种"逐景隔离"的设计让单景失败不影响整批，且失败原因可追溯。

7.3 拓展资源

以下外部资源对理解 GeoTIFF 标签与 MATLAB 地理参考模型有直接帮助（均在新窗口打开）：

• MathWorks 官方文档：geotiffwrite
• MathWorks 官方文档：georefcells
• OGC GeoTIFF 1.1 标准
• GDAL GeoTIFF 驱动说明
• PROJ 官方文档
• YouTube：MATLAB 地理参考写入教程检索

八、对比与局限：MATLAB 路径 vs GDAL 路径

MATLAB 的 geotiffwrite 与 GDAL 的 gdal_translate / gdalwarp 在语义上存在缝隙：

MATLAB 的优势在于把尺寸与范围绑定成对象，减少手工出错；局限在于对非标准投影或自定义 WKT 的支持不如 GDAL 灵活。当遇到 MATLAB 无法识别的投影时，可先用 GDAL 转成 EPSG 标准码，再回到 MATLAB 处理。

本文评述：两条路径不应被对立看待。工程上更稳妥的做法是"MATLAB 负责数值处理，GDAL 负责格式与投影兜底"，用 system() 调用 GDAL 做最终校验，形成双保险。

九、结论与展望

proj2gtif 第 17 行错误很少是投影库本身的问题，而是 SpatialRef、R、数据矩阵 三元契约被破坏后的下游表现。可操作的结论有三条：

1. 构造 R 时始终用 size(A,[1 2])，并优先用 georefcells/maprefcells 一次性绑定尺寸与范围；
2. 写入前做尺寸、范围、方向、语义四类断言，把异常拦截在编码之前；
3. 对来源不明的数据，先用模板影像的参考对象复用，再按需重采样。

展望来看，随着 GeoTIFF 1.1 与 PROJ 9.x 的普及，MATLAB 与 GDAL 在投影表达上的差距正在缩小；但"对象化参考"与"标签化参考"的语义差异仍会长期存在。工程上应把一致性校验沉淀为团队级工具函数，而非依赖个人经验。