简介一份基于MATLAB的卫星数据读取与可视化示例包面向海洋、气象、遥感等领域需要处理卫星资料的科研人员、工程师和高校学生。资源以真实海表面温度数据为例完整演示了从数据文件读取、变量解析、坐标处理到绘制分布图的流程可帮助初学者快速掌握卫星图像读取与出图的核心思路避开常见的数据格式陷阱和绘图参数误区。压缩包共三个文件一个wmv格式的操作录屏便于直观跟随实际运行环境一个read.m脚本提供可直接修改和复用的出图代码一个doc说明文档补充介绍数据来源、程序逻辑和操作要点。整个压缩包仅1.38MB轻量高效适合碎片化学习。目前已有634人学习下载尤其适合刚接触MATLAB卫星数据处理、希望快速跑通完整流程并理解程序结构的用户。借助该示例包可以大幅缩短从零调试的时间快速搭建自己的卫星数据读取与可视化脚本框架也便于扩展到其他卫星数据的分析场景。1. 从“例子-卫星数据.rar”开始MATLAB 卫星数据读取的真实入口很多人拿到的第一个卫星数据集不是一个目录结构清晰的文件夹而是一个 .rar 压缩包。解压开以后里面有 .tif、.nc、.hdf、.mat既有卫星图也有信噪比序列命名还不太规律。这时候你会发现 MATLAB 的 imread、load 在这里不够用imread 读不了带地理坐标的 GeoTIFFload 只能碰 .matncread 和 h5read 又要先弄清楚文件内部结构。下面按一条完整链路往下走先拆包再盘点格式然后给出每种格式的最小读取代码、批量组装时间序列的思路最后是画图和数值验证。适合做遥感、气象和通信信号分析的工程师目标是把“能打开文件”变成“能拿数据做计算做统计画对图画准数”。2. 先拆包再盘点把“例子-卫星数据.rar”解放成 MATLAB 认识的文件集2.1 .rar 解压MATLAB 没有原生 rar 接口时的处理方式先说一个现实问题MATLAB 自带的 unzip 只支持 zip不支持 rar 格式。所以第一步通常不是写 MATLAB 解析代码而是调用外部解压工具。我一般优先用 7-Zip因为它开源且命令行接口稳定Window 和 Linux 都有对应版本。命令行调用方式如下# Windows 下优先找 7-Zip 安装目录里的命令行程序 C:\Program Files\7-Zip\7z.exe x -oD:\satellite_workspace -y D:\examples\例子-卫星数据.rar # Linux 下 p7zip 提供 7z 命令 7z x -o/home/user/satellite_workspace -y /home/user/example_satdata.rar参数说明x表示解压并保留压缩包内目录结构-o指定输出目录注意-o后面不能有空格-y让解压过程覆盖同名文件时不再询问。路径里有空格或中文时务必用双引号包住否则命令会被错误拆分成多个参数这是最常见的失败原因。在 MATLAB 里自动化这一步推荐先判断外部程序是否存在再拼接命令% 在 MATLAB 中调用 7-Zip 解压 .rar 压缩包 zip7Path C:\Program Files\7-Zip\7z.exe; archive D:\examples\例子-卫星数据.rar; outDir D:\satellite_workspace; if isfile(zip7Path) if ~exist(outDir, dir) mkdir(outDir); end cmd sprintf(%s x -o%s -y %s, zip7Path, outDir, archive); [status, result] system(cmd); if status ~ 0 error(解压失败请查看 system 返回信息%s, result); end fprintf(解压完成: %s\n, outDir); else error(未找到 7z.exe请先安装 7-Zip 并确认安装路径); end这段代码里用isfile检查 7z.exe 是否存在比exist(x,file)更精确不会把同名文件夹误判成可执行文件。system返回的status为 0 表示命令成功非 0 时要把result打出来排查错误信息往往直接指向压缩包损坏或权限不足。如果服务端没有图形界面安装 7-Zip 的字符版本后命令用法完全相同。2.2 用 dir 递归通配符盘点卫星数据格式与文件数量解压完先别急着读数据跑一个递归遍历脚本看看这个数据包里到底有哪些格式各有多少个文件rootDir D:\satellite_workspace; files dir(fullfile(rootDir, **, *)); % 递归列出全部内容 files files(~[files.isdir]); % 只保留文件 exts regexpi({files.name}, \.([a-z0-9])$, tokens, once); exts string(exts); [G, extList] findgroups(exts); counts splitapply(numel, exts, G); tally table(extList(:), counts(:), VariableNames, {Extension, Count}); tally sortrows(tally, Count, descend); disp(tally);这段代码用到 MATLAB 在dir函数中支持的**递归模式可以一次扫出所有子目录里的文件。regexpi提取文件名的最后一段扩展名结果转成 string 数组后用findgroups和splitapply做分组计数。运行后会看到 .tif、.nc、.h5、.mat 各自的数量。这个环节的目的不是欣赏统计数字而是确认该为哪几种格式写读取器。很多数据包里还会出现 .tfw、.hdr、.xml 这类辅助文件它们不是数据本体误当成 GeoTIFF 处理会报错盘点阶段就能提早发现。2.3 按扩展名分工MATLAB 卫星数据读取接口对照数据格式占比明确后下一步是为每种格式选择合适的函数。整理成一张接口对照表开发时照着查即可扩展名典型数据内容首选函数注意事项.tif/.tiffLandsat 等光学影像、地表反射率readgeoraster返回像素矩阵和地理栅格引用对象 R需要 Mapping Toolbox.nc大气、海洋再分析产品ncinfo / ncread内置支持无额外工具箱要求.h5/.he5MODIS 切块数据、部分大气探测产品h5info / h5read读取路径带 group结构比 NetCDF 复杂一层.hdf早期 HDF4 产品hdfread / hdfinfoHDF4 接口维护较保守优先确认能否转换格式.mat信噪比序列、轨道参数、中间处理结果load / matfile注意嵌套结构体和变量名冲突这张表是读取架构的依据。实际项目里我绝大多数情况只处理三类GeoTIFF、NetCDF/HDF5、MAT。.tfw这类辅助文件通常不需要单独读GeoTIFF 的参考信息已经内嵌在文件里。3. 三种高频卫星数据格式的最小读取代码3.1 GeoTIFFreadgeoraster 从像素值换算到经纬度GeoTIFF 是把遥感影像和地理参考信息封装在同一个文件里的格式。MATLAB 读取它的标准做法是[A, R] readgeoraster(LC08_20241105_B10.TIF); % A 是二维灰度矩阵 % R 是地理栅格引用对象保存了坐标转换所需的全部参数 size(A) R.RasterSizeR对象里最常用到的属性包括LatitudeLimits和LongitudeLimits对应整个栅格覆盖的纬度经度范围CellExtentInLatitude和CellExtentInLongitude对应每个像素代表的实际尺寸单位是度ProjectedCRS记录投影坐标系。凡是涉及“把像元位置换算成经纬度”的计算都应该直接依赖 R而不是自己写死角点坐标。同一类产品在不同版本中栅格边界经常微调硬编码会让你的数据产品维护成本直线上升。如果当前环境没有 Mapping Toolboxreadgeoraster会直接报错。这时可以用 MATLAB 内置的Tiff类退回读像素t Tiff(LC08_20241105_B10.TIF, r); A read(t); close(t);这样拿到的仍是没有地理参考的纯矩阵需要自己解析投影模型才能把像素位置对应到经纬度非常容易出错。这里的建议很直接要做卫星数据处理先把 Mapping Toolbox 装上后续省出的排错时间远远大于安装成本。3.2 NetCDF 与 HDF5ncinfo 定位变量ncread 取数NetCDF 和 HDF5 文件的难点在于“不知道读哪个变量”。先用ncinfo看顶层结构f MODIS_Surface_Reflectance_A2004001.nc; info ncinfo(f); for k 1:numel(info.Variables) fprintf(%12s | dims: %s\n, info.Variables(k).Name, ... strjoin(info.Variables(k).Dimensions.Name, , )); end注意ncinfo返回的Variables是结构体数组每个元素包含Name、Dimensions、Attributes等字段。以 MODIS 地表反射率产品为例数据变量通常叫sur_refl_b01、sur_refl_b02经纬度要么漂在变量里要么独立叫Latitude、Longitude。取数时还要顺带读变量属性因为很多产品的缺省值不是纳姆NaN而是-9999并且物理量本身带有缩放data ncread(f, sur_refl_b01); scale ncreadatt(f, sur_refl_b01, scale_factor); offset ncreadatt(f, sur_refl_b01, add_offset); % 先还原物理值再屏蔽填充坏值 data double(data) * scale offset; data(data -9999) NaN;HDF5 的读取逻辑类似区别是变量路径要带 group 层级plist h5info(AIRS_L2_20250101.h5); % 先打印 plist. Groups 结构确认变量所在路径 data h5read(AIRS_L2_20250101.h5, /standard/air_temp); lon h5read(AIRS_L2_20250101.h5, /geolocation/longitude);先通过h5info查看 group 名称再拼完整路径。直接用h5read猜路径失败率极高因为不同产品对变量组织方式的约定差别很大。3.3 .mat 卫星数据load 之后先看层级再取字段.mat 文件通常是卫星信号处理链路里的中间产物比如一段时间的信噪比SNR序列。多数情况下load后直接暴露在工作区的并不是一个只有两层结构的变量而是嵌套的结构体S load(example_snr_series.mat); disp(fieldnames(S));假设输出发现顶层变量名是track而track内部又有Lat、Lon、Snr_dB、TimeUTC几个字段顺序索引的意义就体现出来了。直接load有污染当前工作区的风险把整个文件读进结构体S后不仅能避免变量名覆盖还能用isfield做防御性检查。如果 .mat 文件体积接近 1 GB 以上不建议用load一次性载入优先用matfile对象做切片读取这个后文专门讲。3.4 统一读取函数按扩展名分发到不同解析器为了避免每个脚本里都堆一堆if ... elseif ...通常做法是把上面的分支集中到一个readSatFile函数里function out readSatFile(f) %READSATFILE 按扩展名分发的卫星数据读取器 % 输出 struct包含 raw、lon、lat、units [~, ~, ext] fileparts(f); fprintf(读取 %s\n, f); switch lower(ext) case {.tif, .tiff} [out.raw, R] readgeoraster(f); out.lon linspace(R.LongitudeLimits(1), R.LongitudeLimits(2), R.RasterSize(2)); out.lat linspace(R.LatitudeLimits(1), R.LatitudeLimits(2), R.RasterSize(1)); out.units ; case .nc info ncinfo(f); vname info.Variables(1).Name; % 实际使用时应按产品显式指定 out.raw ncread(f, vname); out.lat ncread(f, lat); out.lon ncread(f, lon); out.units unknown; case {.h5, .he5} out.raw h5read(f, /standard/reflectance); out.lat h5read(f, /geolocation/latitude); out.lon h5read(f, /geolocation/longitude); out.units unknown; case .mat S load(f); flds fieldnames(S); if isfield(S, snr) out.raw S.snr; % 常见字段名按实际项目调整 else out.raw S.(flds{1}); % 保底取第一个字段 end out.units dB; otherwise error(不支持的格式: %s, ext); end end这个函数的核心价值是把“每种格式怎么读”的知识收敛到一处。调用时d readSatFile(D:\satellite_workspace\data_01.nc); imagesc(d.raw);后续的批量循环、画图、统计都复用同一套输出结构。换数据集时只需调整这里面的变量名映射而不需要动上百行的分析代码。实际项目里这个函数的第二个迭代版本会加入_FillValue、scale_factor的自动处理以及把维序统一成 lat × lon 的逻辑。4. 从单景到时间序列批量读取“一段时间的卫星数据”真实业务里很快会面对这样一个场景某个目录下连续存放了几十天甚至几个月的文件每帧是一个样本你需要把它们组装成一个时间序列比如提取某段时间内的信噪比变化或逐像元反射率趋势。这一步的核心不是读取本身而是“文件名解析 维度对齐 数据掩码”。4.1 从文件名中解析观测时间假设文件名是20250101_snr.nc这类带日期戳的命名提取观测时间的做法dataDir D:\satellite_workspace\snr_by_day; files dir(fullfile(dataDir, *.nc)); obsTime datetime(NaT(size(files))); for k 1:numel(files) [~, baseName, ~] fileparts(files(k).name); tok regexp(baseName, (20\d{2})(\d{2})(\d{2}), tokens, once); if ~isempty(tok) obsTime(k) datetime(str2double(tok{1}), str2double(tok{2}), str2double(tok{3})); end end % 剔除解析失败的文件 validTime ~isnat(obsTime); files files(validTime); obsTime obsTime(validTime);正则(20\d{2})(\d{2})(\d{2})一次捕获年、月、日三段数字datetime(Y,M,D)接收数值三元组生成datetime对象。解析失败的先保留为 NaT最后统一过滤避免循环中途判断影响代码可读性。如果文件名里确实没有日期退而求其次用文件系统的修改时间但要注意有些批处理任务会把同一天处理的文件写成“最近修改”逻辑上不够严谨。4.2 把多帧数据堆成三维数组或 timetable拿到的是栅格数据时用三维矩阵堆叠最直接。以 MODIS 地表反射率为例每帧是nLat × nLon几十天的数据就堆成nLat × nLon × nT% 先用第一帧确定尺寸 first readSatFile(fullfile(dataDir, files(1).name)); nLat size(first.raw, 1); nLon size(first.raw, 2); nT numel(files); cube NaN(nLat, nLon, nT); for k 1:nT d readSatFile(fullfile(dataDir, files(k).name)); if isequal(size(d.raw), [nLat, nLon]) cube(:, :, k) d.raw; else warning(第 %d 帧尺寸异常已跳过: %s, k, files(k).name); end end % 逐像素计算多年均值 meanMap mean(cube, 3, omitnan); % 查看每个时刻的有效覆盖率 validCount squeeze(sum(~isnan(cube), [1 2])); plot(obsTime, validCount, o-);cube预先用 NaN 填充好处是后续所有统计天然带有掩膜。尺寸不一致是最容易翻车的点卫星轨道宽度变化、裁切范围不同都会导致某些帧出现偏移。直接跳过用 warning 挑出来比强行塞进去让拼接结果错位要好。如果后期需要更精细的对齐可以用mapresize或imresize把栅格统一到参考网格但做这一步之前先搞清楚到底是谁的坐标系发生变化。对于信噪比这类的点位序列不需要三维矩阵用timetable更顺手T timetable(obsTime, snrData, VariableNames, {SNR}); T sortrows(T, Time); % 把非均匀采样规整到小时 Th retime(T, hourly, mean);retime支持hourly、daily、minutely等多种时间粒度第二参数的聚合方式可以是mean、min、max也可以传函数句柄。4.3 填充值、缩放因子与地理维序的三大坑把读取流程做对本质上是在和三类数据本身的问题做对抗填充值没处理。很多产品缺测区域写的是-9999直接mean会把时间序列拉低到离谱的程度。正确做法是先读_FillValue或missing_value属性再把等值像素替换成 NaN推荐把这步收进统一读取函数而不是散落在各分析脚本里。scale_factor / add_offset 漏乘。MODIS 这类数据为了压缩体积会以整数存储真实物理值读取后必须执行data * scale offset。忘了这步画出来的图明暗层次都在但颜色条上的单位全错。这个错误最具迷惑性因为图像看起来“正常”。lat 和 lon 的维度顺序。有些数据文件变量维度是 lon × lat有些是 lat × lon读取后到底要不要转置取决于你后续怎么索引。建议在统一读取函数的输出中固定为 lat × lon把这个决定做在源头% 假设 info 显示变量维序是 [lon, lat] % 则读取后转置一次后续所有分析代码都按 lat × lon 处理 raw ncread(f, reflectance); % 注意转置前要确保 lon 是第一维避免误用这段转置逻辑同样应该收进读取函数而不是靠每个脚本各自祈祷。5. 卫星图读取之后的地理可视化与交叉验证读对了数据还只是第一步。直接imagesc(A)得到的坐标轴是像素索引图和地图之间没有任何映射关系。要“画出来像一张遥感图”至少要把经纬度坐标编织进去。5.1 worldmap geoshow 画带经纬度的卫星图Mapping Toolbox 里最常用的组合是worldmap设定制图范围和投影方式geoshow绘制网格化数据lat d.lat; lon d.lon; A d.raw; figure; worldmap([min(lat) max(lat)], [min(lon) max(lon)]); geoshow(lat, lon, A, DisplayType, texturemap); % 叠加海岸线作为参照 coast load(coastlines); plotm(coast.lat, coast.lon, k, LineWidth, 0.5); colormap(parula); colorbar;worldmap(latlim, lonlim)创建一个带投影信息的地图坐标轴默认投影适合中低纬度展示。geoshow的DisplayType设为texturemap时二维矩阵会被当作纹理贴到经纬度网格上图像自动处于经纬度坐标体系中。plotm是地图坐标系下的专用绘图函数这里不能用普通plot叠加海岸线否则坐标基准不统一点位会错位。输出图的横纵轴自动显示经纬度颜色条对应物理值单位这才是能拿去做报告的卫星图。5.2 没有 Mapping Toolbox 时的替代画法没有安装 Mapping Toolbox 时也有临时方案只是没有投影变换能力只能按经纬度线性展开figure; imagesc(lon, lat, A); set(gca, YDir, normal); % 把纬度方向反转回北在上 xlabel(Longitude (°)); ylabel(Latitude (°)); axis tight; colorbar;imagesc(lon, lat, A)会自动把横轴映射到经度范围、纵轴映射到纬度范围但 MATLAB 默认的 Y 轴是递增方向朝上因此必须用set(gca,YDir,normal)把南在上翻转为北在上。注意这个方案只在数据本身没有做投影处理时才成立。如果数据产品是 UTM 投影横纵轴代表的是东向和北向坐标必须先把投影坐标转换为经纬度再作图不能直接用像素或米值替代经纬度。5.3 用锚点经纬度做读取结果的数值验证画图只能看出数据大概形态真正证明自己读对了的是数值验证。常见做法是选一个已知地理位置锚点比如沿海站点或城市中心查出经纬度再读取卫星数据对应像素值与外部参考值对比。有 Mapping Toolbox 时可以直接用latlon2pix[row, col] latlon2pix(R, anchorLat, anchorLon); row round(row); col round(col); pixelValue A(row, col);latlon2pix的输入R是readgeoraster返回的地理栅格引用对象纬度在前、经度在后。没有工具箱时手工换算col round((anchorLon - lon(1)) / (lon(2) - lon(1))) 1; row round((anchorLat - lat(1)) / (lat(2) - lat(1))) 1;注意lat数组的方向问题。GeoTIFF 的纬度数组通常从北向南排列lat(1)是北端最大值此时手动计算row要把分子倒过来或者直接用nLat - round(...)变换。这套换算的精度上限是半个像元适合验证数据读取正确性不适合做精确定量定位。比对时相对误差保持在千分位量级基本可以确认整条读取链路没有系统性错误。这一步必须放在批量处理之前而不是之后——数百万像素算到一半才发现坐标系反了返工成本不可接受。6. 进阶给卫星数据读取结果加一个可回查的持久化收尾前五章解决的是“从压缩包到数值”的链路这章解决的是“如何让读出来的东西能反复用”。实际项目里经常遇到的情况是原始数据在服务器上脚本跑完退出工作区第二天想换个参数重新分析又要从头解压、读取、校验。这里有一个更高效的收尾方式把第一次读取并完成基础处理后得到的结果连同元数据一起落盘后续所有分析脚本直接基于这份中间产物运行。6.1 用 matfile 做大文件的增量落盘当数据总量达到几个 GB 时save整个工作区会触发内存复制非常容易把桌面端 MATLAB 顶到内存上限。换成matfile对象按切片写入 v7.3 格式的 .mat 文件matObj matfile(processed_satdata.mat, Writable, true); % 预先分配好维度后续按帧写入 matObj.cube NaN(nLat, nLon, nT, single); for k 1:nT matObj.cube(:, :, k) squeeze(cube(:, :, k)); endmatfile不会把整个文件载入内存写入时只操作对应磁盘区间。首帧写入前必须先通过赋值确定变量的维度之后再逐帧覆盖指定切片这样即使是 30 GB 的三维数组也能在普通工作站上跑完。读取时用matObj.cube(:, :, 5)拉取单个时间帧RAM 占用始终可控。6.2 用 timetable 统一时间戳并做重采样第 4 章里提到过timetable在持久化收尾阶段它是更合适的统一出口。把时间列、经纬度列、物理量列组装进同一个表对象后续筛选、重采样和绘图都基于它完成T timetable(obsTime(:), snr(:), lat(:), lon(:), ... VariableNames, {SNR, Lat, Lon}); T sortrows(T, Time); Th retime(T, daily, mean);retime在带有 NaN 的时间戳上会返回 NaN不会自动插值这对卫星数据这种带云的观测来说反而是优点能天然保留缺失区间。如果业务上必须填补再加fillmissing显式指定插值方法避免隐式填出一个“假趋势”。6.3 一个收尾技巧用函数句柄做读取分派表最后一招是给前面写的readSatFile做一个更可扩展的升级用containers.Map存储扩展名到函数句柄的映射避免以后的格式分支都堆进同一个 switchreaders containers.Map(... {.tif, .nc, .h5, .mat}, ... {readGeoTIFF, readNetCDF, readHDF5, readMAT}, ... UniformValues, true); function out readGeoTIFF(f) [out.raw, R] readgeoraster(f); out.lon linspace(R.LongitudeLimits(1), R.LongitudeLimits(2), R.RasterSize(2)); out.lat linspace(R.LatitudeLimits(1), R.LatitudeLimits(2), R.RasterSize(1)); end调用侧简洁为一个动作[~, ~, ext] fileparts(fileName); d readers(lower(ext))(fileName);以后新增格式只要写一个新函数并在 Map 里注册主流程代码完全不动。收尾阶段最后再执行一次whos(-file, processed_satdata.mat)核对落盘数据的变量名和维度确认无误后再把中间文件接入下游分析。这套做法可以让整个卫星数据读取流程在几周后再回看时仍然能够快速定位到“数据从哪来格式怎么读结果在哪里”。本文还有配套的精品资源点击获取