简介:面向需要批量处理NC气象/遥感数据的MATLAB用户,这套脚本可高效将月度或年度NC文件转为GeoTIFF,并支持月度单独导出与年度合成导出两条输出路径。它能灵活应对两类数据组织方式:当单个NC内包含12个月数据时,可逐月导出月平均/月总量tif,也可将12个月求和或平均后生成年度tif;当单个NC内包含三年36个月数据时,既能按月批量导出各月文件,也能每年聚合12个月后输出三张年度tif。压缩包内共2个m文件,整体大小仅2KB,代码注释详尽,使用者只需修改数据输入输出路径即可在MATLAB中直接运行。目前已有548人学习下载;脚本结构清晰、参数改动少,可显著减少逐月手动读取与转换的重复工作,适合气象、水文、遥感等领域的科研人员快速建立NC到tif的批处理流程,并据此灵活生成月尺度或年尺度栅格结果。
1. 被反复问到的NC转TIF,先看懂数据再动手
做气象、遥感或者水文数据处理的人,绕不开NC文件的批量转换。我见过不少同事和同行,拿到一批月度的NC数据,第一反应是用ArcGIS的Make NetCDF Raster Layer工具一张张导出,遇到36个月的文件就手动拉12次,遇到年度合成需求再写一遍栅格计算器公式。整个过程繁琐不说,中间一旦有文件命名不统一或者缺测值处理不一致,后面的分析基本白做。
这个NCtoTIFF工具包做的事情很直接:读取月度或年度NC文件,按月批量导出TIF,或者把12个月的数据做求和、平均后输出年度TIF。素材包里包含NCtoTIFF.m和plNCtoTIFF.m两个脚本,代码注释详细,改输入输出路径就能直接跑。但真正值得关注的是它对时间维的处理方式——这在批量处理里才是最容易踩坑的地方。
2. NC文件的结构解析与读取前的准备
2.1 先搞清楚NC文件里的维度关系
NC格式本质上是一种自描述的数组存储结构,文件内部包含了变量、维度和属性。用MATLAB处理NC文件,第一个动作永远是ncinfo而不是ncread。很多新手拿到文件就急着读数据,结果发现读出来的数组维度顺序和自己想的不一样。
对月度和年度NC文件来说,最常见的变量组织方式是经度(lon)、纬度(lat)、时间(time),以及一个或多个气象变量,比如降水、气温、蒸发等。维度顺序通常是[lon, lat, time]或者[time, lat, lon],这取决于数据生产方的写入习惯。在MATLAB里,ncread按文件中的维度声明顺序返回数组,因此如果不先检查ncinfo里的Dimensions和Variables信息,后续所有循环都会出错。
% 第一步:查看NC文件结构 info = ncinfo('precip_monthly_2001.nc'); disp({info.Variables.Name}); disp(info.Dimensions);这段代码输出的Variables列表告诉你文件里有哪些变量,Dimensions告诉你每个维度的名称和长度。例如time = 12说明这个文件包含12个月的数据,lon = 360、lat = 180说明是0.5度分辨率的网格数据。
2.2 变量名的自动识别
不同来源的NC文件,变量命名差异很大。有的叫precip,有的叫pr,有的叫rainfall。如果代码里写死变量名,换一批数据就得改一次源码。常见的做法是在读取前先扫描所有变量,自动锁定目标变量。
% 自动选择第一个三维变量作为目标变量 varNames = {info.Variables.Name}; varSizes = cellfun(@(v) numel(info.Variables(strcmp(varNames, v)).Size), varNames); targetVar = varNames{find(varSizes == max(varSizes), 1)}; disp(['目标变量: ', targetVar]);这里用varSizes比较每个变量的数组规模,因为气象要素变量通常维度最多、规模最大,坐标变量(lon、lat、time)往往是一维数组。自动选择最大规模的变量,可以避免手动修改变量名参数。如果文件里有多个三维变量,比如降水和温度共存,那就需要根据业务场景手动指定变量名。
2.3 缺测值处理
NC文件里的缺测值通常存放在变量属性_FillValue或missing_value中。读取数据后必须把这个值替换为NaN,否则转出来的TIF里会出现异常极值,ArcGIS或QGIS打开后整幅图可能是花的。
% 读取变量属性并处理缺测值 varData = ncread(ncFile, targetVar); fillVal = ncreadatt(ncFile, targetVar, '_FillValue'); % 有的文件用 missing_value varData(varData == fillVal) = NaN;建议把_FillValue和missing_value两个属性都查一遍,有些数据生产方只写其中一个。查找失败时可以用exist判断属性是否存在,避免报错中断整个批量流程。
3. 月度文件批量转TIF的完整实现
3.1 文件遍历与月份维解耦
月度NC文件的常见组织方式有两种:一种是每个文件包含12个月的逐月数据,一种是每个文件只包含一个月份的数据。工具包支持的是第一种场景,也就是单个NC文件内时间维长度为12。但无论哪种方式,核心逻辑都是先把时间维取出来,再逐月写入TIF。
% 月度批量导出主循环 ncFile = 'precip_mon_2001.nc'; lon = ncread(ncFile, 'lon'); lat = ncread(ncFile, 'lat'); months = 12; for m = 1:months % 按时间索引读取单月数据,lon和lat维度也一并取出 data = squeeze(ncread(ncFile, targetVar, [1 1 m], [Inf Inf 1])); data = data'; % 转置为 [lat, lon],符合GeoTIFF行列习惯 % 构造输出文件名 outFile = sprintf('precip_2001_%02d.tif', m); % 写入TIF,R是地理参考对象 geotiffwrite(outFile, data, R); end这段代码有四个关键参数需要说明。
ncread的第三个参数是起始索引,第四个参数是读取长度。[1 1 m]表示从第一个经度、第一个纬度、第m个月开始读取,[Inf Inf 1]表示经度和纬度维度全读,时间维度只读1个值。用squeeze去掉时间维为1的维度,得到的才是二维矩阵。
data = data'这一步必须注意。NC文件里维度顺序是[lon, lat],而GeoTIFF的行列习惯是[lat, lon],即第一维是行(纬度方向),第二维是列(经度方向)。如果不转置,出来的TIF会发生90度旋转,而且经纬度坐标全错位。
R是地理参考对象,由georasterref构造或者从已有的TIF模板中读取。如果NC文件里没有包含投影信息,直接用georasterref('RasterSize', size(data), 'Latlim', [min(lat) max(lat)], 'Lonlim', [min(lon) max(lon)])构造一个WGS84的参考对象。
%02d保证月份编号补零为两位数,这样文件名排序时不会出现10月排在2月前面的问题。
3.2 坐标数组的重采样问题
NC文件里的经度经常是0到360,而TIF的经度范围是-180到180。如果直接使用georasterref构造的经纬度边界,两者会差180度。常见做法是读取坐标后先判断范围,再做转换。
lon = ncread(ncFile, 'lon'); if max(lon) > 180 lon = lon - 360; % 0-360转-180-180 end转换后要注意数据矩阵的排列顺序。0度经线在中间时,数据矩阵的左右两半需要交换,否则地图上显示的位置仍然不对。处理方式是把数组左右对调,或者用circshift重新按经度排序。这个细节在直接下载现成数据时经常被忽略。
4. 年度合成导出与三年36个月文件的批量处理
4.1 月度求和与平均的业务差异
年度合成有两种常见方式:累加求和与平均。降水、蒸散量这类累积型变量用求和,比如年降水量就是12个月降水之和;气温这类强度型变量用平均,比如年均温是12个月气温的平均值。工具包同时支持两种方式,只需要在参数中指定。
% 年度合成:先循环读取12个月数据,再合成 yearData = zeros(size(data)); % 预分配数组 for m = 1:12 monthData = squeeze(ncread(ncFile, targetVar, [1 1 m], [Inf Inf 1])); monthData(monthData == fillVal) = NaN; yearData = yearData + monthData; end % 如果聚合方式是'average',除以12;'sum'则不用除 if strcmp(aggMethod, 'mean') yearData = yearData / 12; end geotiffwrite(sprintf('precip_2001_annual_%s.tif', aggMethod), yearData', R);这里引入了一个业务问题:平均操作应该除以12还是天数加权平均。对逐月数据来说,每个月的天数不同,严格做法是加权平均:天数作为权重。但气象领域处理逐月气候态时,常用算术平均简化处理。工具包默认算术平均,注释里也写了这个取舍,使用者要根据自己的业务场景判断。
求和操作相对安全,但要特别留意缺测值分布。如果某个月的数据存在大范围缺测,直接累加会污染年度总量。常见做法是统计逐月NaN掩膜,只有12个月都有有效值的像元才参与年度合成。实现上可以在累加时同时累加有效数据计数。
validCount = zeros(size(data)); for m = 1:12 monthData = squeeze(ncread(...)); validMask = ~isnan(monthData); monthData(~validMask) = 0; yearData = yearData + monthData; validCount = validCount + validMask; end yearData(validCount < 12) = NaN; % 有效月份不足12的像元设为缺测4.2 三年36个月文件的按年拆分
部分数据集会把36个月的数据放进一个NC文件里,比如某区域2001年至2003年的逐月序列。此时需要先按时间维索引将数据切分为三年,再套用月度或年度的处理逻辑。切分的关键在于时间维索引的计算。
% 三年36个月批量处理 ncFile = 'precip_2001_2003.nc'; timeLen = 36; for yr = 1:3 % 第几年 % 计算每一年的月份索引范围 startIdx = (yr - 1) * 12 + 1; endIdx = yr * 12; % 年度内逐月导出 for m = startIdx:endIdx data = squeeze(ncread(ncFile, targetVar, [1 1 m], [Inf Inf 1])); yearNum = 2000 + yr; monthNum = m - (yr - 1) * 12; geotiffwrite(sprintf('precip_%d_%02d.tif', yearNum, monthNum), data', R); end end这段代码里,startIdx和endIdx的计算逻辑是固定的:第1年索引1到12,第2年索引13到24,第3年索引25到36。写代码时不要使用mod(timeIndex, 12)来推断月份,因为mod(12, 12)的结果是0,而实际应该是12月。我一般会让月份索引从1到36连续编号,单独用if判断每个循环变量的归属。
4.3 输出文件命名规范
批量处理的项目,文件命名往往是最后最麻烦的一环。我习惯把变量名、年份、月份、聚合方式全部编入文件名,这样后续用ArcGIS批量导入或者Python脚本做时效分析时,文件名本身就携带了完整的元数据信息。
推荐命名格式是变量名_年份_月份_聚合方式.tif,例如precip_2001_05.tif、precip_2001_annual_sum.tif、precip_2001_annual_mean.tif。其中年份用四位数字,月份用两位数字补零,聚合方式在月度文件里省略。这套规则简明,对中文路径有兼容性,也方便Shell脚本或MATLAB的dir函数做后续匹配。
5. 多年度文件路径构建、批处理扩展与结果验证
5.1 多年度文件的自动路径构建
前面的例子都以单文件为输入,实际项目中往往是20年、30年的数据序列。这时需要先构建一个按年份循环的主流程,再在每个年份内层调用月度或年度子函数。这里有一个细节值得注意:NC文件的年份信息往往不在文件内部,而是在文件名里。
% 多年度批量处理框架 baseDir = 'F:/ClimateData/Precip/'; years = 2001:2020; for y = years ncFile = fullfile(baseDir, sprintf('precip_%d.nc', y)); if ~exist(ncFile, 'file') warning('文件不存在: %s', ncFile); continue; end % 调用NCtoTIFF核心函数 NCtoTIFF(ncFile, 'OutputDir', 'F:/ClimateData/Precip/TIF/', ... 'AggMethod', 'sum', 'ExportMonthly', true, 'ExportAnnual', true); end用exist检查文件是否存在是防止循环中断的保底手段。天气再好的数据目录也可能缺文件,一旦ncread找不到文件,整个批处理就会中断在这个年份上,而且不会自动记录失败日志。更稳妥的方案是把处理结果写进日志文件。
另一个容易被忽略的任务是检查年度文件的时间维长度是否为12。某些年份的数据可能缺某个月,导致ncread读取时越界。批量循环里前置一个长度校验,问题文件直接跳过并记录文件名,比运行中途报错再去排查效率高得多。
% 检查时间维长度是否为12,不等于12则记录警告 timeLen = length(ncread(ncFile, 'time')); if timeLen ~= 12 warning('%s 时间维长度为%d,不是12个月,跳过', ncFile, timeLen); continue; end5.2 批量结果的快速验证
转换完成后,不能只靠肉眼打开几张图就确认没问题。尤其是转置、经度范围转换这类操作,往往是整片区域在空间上产生了系统性偏移,单张TIF上几乎看不出来。
我一般用两种方式快速验证。第一种是原始数据与输出TIF的数值一致性比对,在MATLAB里读取NC文件的一个月和对应TIF文件,计算相关系数和最大绝对差。第二种是直接对比一组多个TIF文件的整体统计特征,比如年均合计数据应该约等于12个月的累加值。
% 验证模块:对比NC与TIF的统计值 ncData = squeeze(ncread('precip_2001.nc', 'precip', [1 1 1], [Inf Inf 1])); tifData = geotiffread('precip_2001_01.tif'); % 去掉缺测值后对比 validIdx = ~isnan(ncData) & ~isnan(tifData); corrVal = corr(ncData(validIdx), tifData(validIdx)); maxDiff = max(abs(ncData(validIdx) - tifData(validIdx))); fprintf('相关系数: %.4f, 最大绝对差: %.4f\n', corrVal, maxDiff);如果同一文件直接转出的TIF相关系数都低于0.999,优先检查转置和数据精度。geotiffwrite默认写出单精度浮点型,而NC文件里的数据可能是双精度。转换时的双精度到单精度截断一般不会影响分析,但如果数据范围很大且关注极值差异,可以在写入参数中指定'Precision', 'double'保留精度。
5.3 时间维数据在MATLAB中的另类处理思路
工具包的逻辑是基于[lon, lat, time]逐索引读取,也就是说NC文件整体作为一个三维数据块常驻内存,然后按需要切片写入TIF。如果文件非常大,比如全球0.25度分辨率36个月的数据,整个数组可能超过数GB,ncread一次性读取会比较吃力。
常见替代方案是按需读取,即每次只从磁盘读取一个时间切片,然后立即写入TIF后释放变量。核心就是本文前面代码里的[1 1 m]、[Inf Inf 1]读取方式,这已经是最省内存的做法。需要特别注意的是,循环读取中不要累积临时变量,例如把每个月的二维数据都保存在一个三维数组里,然后再统一写入。这样内存在单月处理上没问题,但后期如果需要把12个月的数组叠加计算,就要考虑用cumsum增量累积,边读边加,而不是先保存12个数组再做二次循环。
另外,如果你手头有多个NC文件需要先拼接再转换,可以考虑在读取时直接按年份顺序加载进一个预留好空间的三维数组里。比如20年数据共240个月,先zeros(lonLen, latLen, 240)分配空间,然后逐文件写入对应年份的索引段,最后再统一循环输出TIF。这种方式适合内存充足、且后续需要跨年统计分析的场景。
关于输出TIF的坐标系,建议统一使用WGS84地理坐标系。如果原始数据是等经纬度网格,直接让R对象携带经纬度范围即可。对于已经包含投影信息的NC文件,可以用geotiffwrite的扩展参数把坐标参考信息一并写入。在ArcGIS里打开后如果发现位置不对,回到第2节检查经度范围和数据转置。这两类问题占了NC转TIF出错场景的八成以上。
本文还有配套的精品资源,点击获取