Matlab处理遥感影像必看：地理坐标和投影坐标的GeoTIFF读写，别再搞混了！
Matlab遥感影像处理实战：地理坐标与投影坐标的GeoTIFF读写全解析
遥感影像处理中，坐标系的选择与正确读写是许多初学者容易踩坑的环节。今天我们就来深入探讨Matlab环境下如何处理这两种不同坐标系的GeoTIFF文件，从原理到实践，帮你彻底理清思路。
1. 地理坐标与投影坐标的本质区别
在开始代码实操前，我们必须先理解这两种坐标系的根本差异。地理坐标系（Geographic Coordinate System）基于椭球体模型，使用经纬度表示位置，单位为度。而投影坐标系（Projected Coordinate System）则是将三维地球表面投影到二维平面上的坐标系统，常用单位为米。
表：地理坐标系与投影坐标系核心差异对比
| 特征 | 地理坐标系 | 投影坐标系 |
|---|---|---|
| 基础 | 椭球体模型 | 平面投影 |
| 单位 | 度 | 米 |
| 变形 | 存在角度变形 | 存在长度/面积变形 |
| 适用场景 | 全球尺度分析 | 区域尺度制图 |
在Matlab中，这两种坐标系对应的数据结构也不同：
- 地理坐标系：
map.rasterref.GeographicCellsReference - 投影坐标系：
map.rasterref.MapCellsReference
2. 地理坐标系GeoTIFF的读写方法
对于地理坐标系的影像，处理流程相对直接。以下是典型的工作流：
% 读取地理坐标系影像 [A, R] = geotiffread('input_geo.tif'); % 检查坐标系类型 disp(class(R)) % 应显示map.rasterref.GeographicCellsReference % 影像处理（示例：NDVI计算） nir = A(:,:,4); red = A(:,:,1); ndvi = (nir - red) ./ (nir + red); % 输出处理结果 output_filename = 'output_geo.tif'; epsg_code = 4326; % WGS84地理坐标系 geotiffwrite(output_filename, ndvi, R, 'CoordRefSysCode', epsg_code);关键点说明：
geotiffread会自动识别坐标系类型并返回相应的R对象- 地理坐标系输出时必须指定正确的EPSG编码（如WGS84为4326）
- 使用
CoordRefSysCode参数传递EPSG编码
提示：常见地理坐标系EPSG编码
- WGS84: 4326
- CGCS2000: 4490
3. 投影坐标系GeoTIFF的特殊处理
当遇到投影坐标系时，直接套用地理坐标系的写法会导致错误。正确的处理方式如下：
% 读取投影坐标系影像 [A, R] = geotiffread('input_proj.tif'); info = geotiffinfo('input_proj.tif'); % 检查坐标系类型 disp(class(R)) % 应显示map.rasterref.MapCellsReference % 影像处理（示例：波段运算） processed_data = A(:,:,1) * 0.5 + A(:,:,2) * 0.3; % 输出处理结果 output_filename = 'output_proj.tif'; geotiffwrite(output_filename, processed_data, R, ... 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);常见问题解决方案：
- 报错"The input, R, is a map.rasterref.MapCellsReference object..."时：
- 确认是否使用了
geotiffinfo获取原图信息 - 改用
GeoKeyDirectoryTag参数而非EPSG编码 - 确保从原图info中提取正确的GeoKeyDirectoryTag
- 确认是否使用了
4. 高级应用：自定义坐标参考系统
有时我们需要为新建的栅格数据指定全新的坐标参考系统。以下是两种实现方式：
4.1 创建新的地理坐标系参考
% 定义地理范围 latlim = [35 40]; lonlim = [110 120]; rasterSize = [1000 1500]; % 创建地理参考 R = georefcells(latlim, lonlim, rasterSize); % 创建随机数据并保存 data = rand(rasterSize); geotiffwrite('new_geo.tif', data, R, 'CoordRefSysCode', 4326);4.2 创建新的投影坐标系参考
% 定义投影参数 proj = geotiffinfo('sample_proj.tif').GeoTIFFTags.GeoKeyDirectoryTag; % 定义空间范围 xlim = [300000 350000]; ylim = [2500000 2550000]; rasterSize = [800 600]; % 创建投影参考 R = maprefcells(xlim, ylim, rasterSize); % 创建数据并保存 data = rand(rasterSize); geotiffwrite('new_proj.tif', data, R, 'GeoKeyDirectoryTag', proj);5. 实战技巧与性能优化
在实际项目中，我们还需要考虑以下进阶问题：
5.1 批量处理多文件
file_list = dir('*.tif'); for i = 1:length(file_list) [A, R] = geotiffread(file_list(i).name); info = geotiffinfo(file_list(i).name); % 处理逻辑... % 自动判断坐标系类型 if contains(class(R), 'Geographic') geotiffwrite(['proc_' file_list(i).name], result, R, ... 'CoordRefSysCode', 4326); else geotiffwrite(['proc_' file_list(i).name], result, R, ... 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag); end end5.2 内存优化策略
处理大型遥感影像时，可采用分块处理：
block_size = [1000 1000]; [A, R] = geotiffread('large_image.tif', 'PixelRegion', {[1 block_size(1)], [1 block_size(2)]}); % 分块处理逻辑...5.3 坐标系转换
有时需要在两种坐标系间转换：
% 地理坐标转投影坐标 [lat, lon] = meshgrid(35:0.01:36, 110:0.01:111); [x, y] = projfwd(info, lat, lon); % 投影坐标转地理坐标 [lat, lon] = projinv(info, x, y);掌握这些核心技巧后，相信你在处理遥感影像坐标问题时能够更加得心应手。在实际项目中，建议先小范围测试代码，确认坐标系处理正确后再进行批量操作。