前阵子做跨省生态区的Sentinel-2影像覆盖统计被一堆T50TLH、T50TMH、T51RTL之类的图幅号折腾得够呛。明明只想查清楚研究区到底涉及多少景数据、分别落在哪几条轨道条带上结果手动对编号对了整整一下午中间还漏了两幅后来在整理数据清单时才发现。所以我就花了一点时间用Python把这套流程彻底自动化了输入一个研究区范围自动计算相交的哨兵2号图幅关联轨道条带信息输出一张统计表和一份矢量边界文件。这篇文章就把这个过程中的概念、代码和坑一次性讲清楚适合经常和Sentinel-2数据打交道、被图幅与轨道编号折磨过的遥感从业者和GIS开发朋友。1. 先把两个概念理清楚轨道条带与图幅格网刚开始接触Sentinel-2数据的人最容易被两套编号系统弄晕一套是轨道号一套是图幅号。这两者经常被混着说但实际上是两码事。这不怪大家因为欧空局ESA的官方数据目录里不管你在哪个平台查数据展示出来的都是图幅编号但数据在卫星上真正拍摄时又是沿轨道连续成像并没有按图幅切分。搞清楚这一层关系后面做统计才不会懵。1.1 Sentinel-2卫星与143条相对轨道Sentinel-2是哥白尼计划里的光学遥感卫星目前业务运行的是2A和2B两颗。单颗卫星的重访周期是10天两颗卫星组网后同一地点最快5天就能再看一次。两颗星运行在同一类太阳同步轨道上轨道高度大约786公里幅宽290公里降交点地方时控制在10:30左右。这个时间点选得很讲究光照角度合适既不太早导致阴影太长也不太晚导致云热对流发展起来对地表反射率产品很友好。关于“143条相对轨道”这个数字可能有人好奇是怎么来的。Sentinel-2的回归周期是10天卫星每天大约绕地球飞行14.3圈10天就是143圈。每一圈在地面上投影成一个固定的轨道条带编号从001到143。同一个相对轨道号在每天的同一时刻会经过地球上的同一个位置。也就是说轨道条带是卫星飞行的“路”图幅是路两边划分好的“格子”卫星沿着既定的路飞过按格子把影像切好卖给用户。从空间关系上看Sentinel-2的成像主要发生在降轨段也就是卫星从北往南飞的时候这时候是白天地面有光照。升轨段大多是夜间一般不作为光学观测使用。所以在做图幅统计时如果手头的轨道信息里带方向字段通常你会发现大部分图幅落在DESCENDING方向这符合预期不用大惊小怪。1.2 MGRS图幅数据产品怎么切片有了轨道路径还需要一套“切格子”的标准这就是MGRS军用网格参考系统。经纬度网格在低纬度和高纬度地区的面积差异很大直接拿它切图会出现有的图块大、有的图块小的问题。MGRS的方案更聪明先把全球按6度经度分成60个UTM投影带再从南纬80度到北纬84度按8度纬度分成纬度带每个UTM带内再按100公里边长切方格子。Sentinel-2的产品图幅ID就是基于这套MGRS规则命名的。一个完整图幅ID由五部分组成比如T50TLH拆开看就很清晰位置示例含义第1-2位50UTM带号全球共60个带从西经180度起算每6度一个带第3位T纬度带从C到X跳过I和O每8度一个带T对应北纬40-48度第4位L列号在UTM带内按100公里宽度从西往东编号第5位H行号从南往北按100公里高度编号说句题外话MGRS的字母表故意去掉了I和O目的是避免和数字1、0混淆所以你会看到图幅ID里永远没有这两个字母。这个细节在做字符串校验时很管用如果从某个系统导出的图幅编号里出现了I或O基本可以断定数据有问题。1.3 轨道条带和图幅的关系轨道条带和图幅之间并不是严格的“一对应一”关系。一条轨道在地面上扫过大几百公里长的范围会跨越多个MGRS图幅反过来由于地球自转和UTM带边界的影响有些图幅的边缘可能被不同轨道的光学影像重叠覆盖。但在任务设计层面每个图幅都对应一个主用的相对轨道号这个对应关系是固定的并记录在欧空局发布的图幅参数文件里。做图幅统计的时候不能只靠图幅ID的字母硬推轨道号最稳妥的办法是查轨道映射表。但现实情况是很多小伙伴手里的图幅网格shp或者GeoJSON里直接带着orbit字段那就能省掉查表的步骤。如果碰巧你的数据里没有这个字段再去下载TILPARTile Parameter文件它本质上是一份XML里面记录了全球每个图幅对应的轨道号解析出来就能关联。2. 动手前需要准备什么写这段之前我先声明一下我不推荐上来就写代码。先把数据和依赖准备好顺手理一下流程后面写起来会顺很多。这个项目看起来只是“做个统计”但数据源、坐标系统、输出格式这些细节一旦没想清楚后面返工是必然的。2.1 数据源清单要做图幅统计至少需要三个东西研究区边界、图幅网格、轨道映射信息。研究区自己准备比如一块林地、一个县、一个流域导成GeoJSON或者Shapefile都行。图幅网格可以从欧空局官网、部分开源遥感项目仓库或者美国USGS的数据分发页面找到全球范围的Sentinel-2图幅格网通常是一个包含所有图幅面要素的矢量文件字段里最好自带TILE_ID和ORBIT字段。如果你的图幅网格里只有TILE_ID没有ORBIT字段那么需要额外准备一份轨道映射表。官方渠道是下载TILPAR文件下载后不需要人工看XML直接用Python解析出每个TILE_ID对应的ORBIT_NUMBER即可。有些第三方工具也直接提供了类似的功能比如一些Sentinel-2数据下载客户端在搜索时就会返回每条影像对应的轨道号但那种方式依赖外部服务不如本地准备一份映射表来得稳定。2.2 环境与依赖这套统计逻辑的技术栈很常规基于GeoPandas就能搞定。GeoPandas封装了空间数据结构和空间计算接口底层的空间索引用的是Shapely和PyGEOS。还需要Pandas做表格聚合如果要导出Excel再补一个openpyxl。装一次就行pip install geopandas pandas openpyxl如果你的环境里没有GeoPandas安装过程可能会自动拉GDAL、Fiona等二进制依赖Windows用户建议直接用官方预编译的wheel或者用conda创建环境会省掉很多编译相关的麻烦。2.3 整体处理流程设计这个工具的处理流程不算复杂一共五步读取研究区矢量统一坐标系读取图幅网格统一坐标系空间关联找出与研究区相交的所有图幅关联轨道信息按轨道分组聚合输出统计表和矢量文件每一步都可以独立调试尤其是空间关联那一步结果直接决定后面统计对不对。我建议在写代码的时候把这五步拆成四个函数load_roi、find_tiles、aggregate_by_orbit、export_outputs。这样测试时单独跑第一步就能发现数据读取问题不用每次从头跑一遍。3. Python实现从研究区到统计表与矢量下面开始写代码。我这里用的是比较典型的GeoPandas写法代码段可以直接复制改成自己的文件路径运行。需要提前说明我假设你的图幅网格文件里已经存在TILE_ID和ORBIT字段。如果字段名不一样把代码里的字段替换成你手头数据的实际字段名即可。3.1 第一步读取研究区并统一坐标系所有空间运算前先统一坐标系。我习惯统一到WGS84经纬度也就是EPSG:4326。原因很简单图幅网格的原始坐标系通常是WGS84研究区的坐标系可能多种多样。统一到4326做空间关系判断图幅边界不会因为投影变形产生大的偏差。import geopandas as gpd # 读取研究区支持GeoJSON/Shapefile roi gpd.read_file(research_area.geojson) roi roi.to_crs(EPSG:4326) print(研究区范围, roi.total_bounds) print(研究区坐标系, roi.crs)total_bounds会返回一个四元组顺序是左、下、右、上也就是经度和纬度的最小、最大值。这一步很有用可以快速判断数据有没有读错。比如你本来研究的是黑龙江结果total_bounds显示东经120度北纬30度那八成是研究区文件本身有问题或者坐标单位不对早点发现能少走很多弯路。3.2 第二步空间关联找出相交图幅图幅网格数据通常很大全球范围可能有几千上万个图幅多边形。直接逐个图幅做相交判断速度会非常慢。GeoPandas的sjoin内部会自动构建R树空间索引先初筛再精算效率高得多。这也是我建议用GeoPandas而不是自己写循环的根本原因。# 读取全球图幅网格 tiles gpd.read_file(S2_tiles_world.geojson) tiles tiles.to_crs(EPSG:4326) # 空间连接找出与研究区相交的图幅 joined gpd.sjoin(tiles, roi, howinner, predicateintersects) # 一个图幅可能与研究区多个子要素相交这里去重 joined joined.drop_duplicates(subset[TILE_ID]).copy() print(相交图幅数量, len(joined)) print(joined[[TILE_ID, ORBIT]].head())这段代码里有几个细节值得说。sjoin的how参数用inner只保留两边能匹配上的记录predicate参数指定判断方式intersects是碰到就算。如果你的研究区特别大比如占了半个省那么边缘上只要擦到一点点皮的图幅也会被算进来这对统计来说一般没问题因为下载影像时你也得覆盖这些边缘图幅。如果只想统计那些中心点落在研究区里的图幅后续可以用图幅中心点和研究区再做一次within判断二选一即可。3.3 第三步关联轨道号与统计如果图幅网格里没有ORBIT字段这一步就需要额外做一次映射表关联。这里我特意演示一下用TILPAR映射表的情况因为真正用起来你会发现网上能找到的很多图幅网格数据并不带轨道信息。# 从TILPAR解析出来的映射表读取成DataFrame orbit_map pd.read_csv(tile_orbit_map.csv) # 取需要的字段去掉重复 orbit_map orbit_map[[TILE_ID, ORBIT, PASS_DIR]].drop_duplicates() # 与图幅统计结果关联 joined joined.merge(orbit_map, onTILE_ID, howleft, suffixes(, _MAP)) # 优先使用原有的ORBIT字段如果没有则使用映射表里的 if ORBIT_MAP in joined.columns: joined[ORBIT] joined[ORBIT].fillna(joined[ORBIT_MAP])关联完成后就可以按轨道号分组统计了。统计维度我建议至少包含轨道号、图幅数量、图幅ID列表。如果映射表里带了升降轨方向也一起分组进去后面做数据筛选时能直接看出哪些轨道是降轨成像。summary joined.groupby([ORBIT, PASS_DIR]).agg( tile_count(TILE_ID, count), tile_ids(TILE_ID, lambda x: ; .join(sorted(x))), geom(geometry, unary_union) ).reset_index() print(summary)这里用到的unary_union是Shapely的操作会把同一轨道下的所有图幅几何合并成一个MultiPolygon。这个合并后的几何就是该轨道在研究区内的实际覆盖范围后续输出矢量时会非常有用。3.4 第四步输出统计表与Excel统计表的输出可以做成两种一份汇总表一份明细表。汇总表按轨道分组方便快速看整体分布明细表保留每一个相交图幅的记录方便下一步逐景下载或者筛选。# 明细表每个图幅一条记录 detail joined[[TILE_ID, ORBIT, PASS_DIR, geometry]] detail detail.sort_values([ORBIT, TILE_ID]).reset_index(dropTrue) # 汇总表按轨道统计 summary_display summary[[ORBIT, PASS_DIR, tile_count, tile_ids]] summary_display.columns [轨道号, 方向, 图幅数量, 图幅列表] # 导出Excel两个sheet with pd.ExcelWriter(sentinel2_tile_stat.xlsx, engineopenpyxl) as writer: summary_display.to_excel(writer, sheet_name轨道汇总, indexFalse) detail[[TILE_ID, ORBIT, PASS_DIR]].to_excel(writer, sheet_name图幅明细, indexFalse)导出来的Excel表格大概长这样轨道号方向图幅数量图幅列表65DESCENDING8T50TLH; T50TMH; T50TNH; ...97DESCENDING6T51RTL; T51RML; T51RNL; ...这里我用了中文列名方便同事直接看。如果你要保留英文再搞后续程序处理建议还是保留英文列名中文列名在部分代码里会带来编码小麻烦这个看个人习惯。3.5 第五步输出矢量成果矢量输出是整个流程里最容易被忽略的部分。很多人做完图幅统计觉得表格就够了但实际做项目汇报时领导或者甲方根本不想看那几行编号他们要看图。把图幅和高亮条带输出成矢量文件叠加在底图上一眼就能看出研究区覆盖了哪些轨道范围。# 输出相交图幅面 detail.to_file(sentinel2_tile_result.gpkg, layertiles, driverGPKG) # 输出按轨道合并后的条带范围 orbit_geom gpd.GeoDataFrame(summary, geometrygeom, crsEPSG:4326) orbit_geom orbit_geom[[ORBIT, PASS_DIR, tile_count, geom]] orbit_geom.columns [轨道号, 方向, 图幅数量, geometry] orbit_geom.to_file(sentinel2_tile_result.gpkg, layerorbit_footprints, driverGPKG) # 输出研究区本身 roi.to_file(sentinel2_tile_result.gpkg, layerroi, driverGPKG)GPKGGeoPackage这种格式最大的好处是可以把多个图层放进同一个文件里不用像Shapefile那样每个图层单独一个文件夹。你最终会得到一个sentinel2_tile_result.gpkg里面有tiles相交的图幅边界、orbit_footprints按轨道合并的条带范围、roi研究区边界三个图层QGIS打开直接叠加看非常直观。4. 实际跑数据时最容易翻车的几个坑这套流程看起来简单但实际运行时会踩到各种奇奇怪怪的坑。有些坑我是在项目汇报前一天晚上才发现的那时候真想摔键盘。下面挨个说说基本能覆盖大部分人的问题。4.1 坐标系不统一导致统计结果为空最常见的问题是空间关联结果为空。你确认研究区明明在境内但sjoin一条记录都返回不了。这种时候先打印研究区和图幅网格的CRS如果一个是EPSG:4326另一个是EPSG:32650那大概率坐标系不一致。另一个坑是坐标单位导致的范围异常研究区坐标可能是米图幅坐标可能是度两者范围差几个数量级空间索引直接匹配不上。解决办法就一句所有数据在运算前统一to_crs(EPSG:4326)。还有一种隐藏的坐标系问题发生在图幅网格本身。部分非官方来源的图幅网格虽然标注EPSG:4326但几何顶点坐标其实是经纬度乘以了比例因子比如除以1000这种情况total_bounds一眼就能看出来数值明显不在正常经纬度范围内。遇到这种数据建议直接换数据源比你去纠偏划算得多。4.2 研究区跨UTM带时需要注意什么一个研究区如果很大比如跨两三个省它必然跨越多个UTM投影带也就对应着一长串不同UTM带号的图幅。MGRS图幅ID前两位就是UTM带号所以这种情况下统计结果里会出现T50、T51、T52等多个带号的图幅这是正常的。但要注意在跨UTM带的区域轨道条带和图幅边界关系会变得比较复杂。同一个轨道条带在某些纬度可能正好跨在UTM带边界上导致一条轨道对应两列图幅而不是通常的一列图幅。这种情况在空间关联后按轨道聚合时会发现某条轨道的图幅列表里出现了跨带编号不要觉得是程序错了这是MGRS系统本身的特点。4.3 图幅网格和映射表字段衔接不上轨道映射表里的TILE_ID格式可能和你图幅网格里的不一致。比如图幅网格里是“T50TLH”映射表里却是“50TLH”少了开头的T或者映射表里的TILE_ID是数字编码。这些字段不一致的问题在merge之后会出现大量缺失值。建议在关联前先统一格式比如统一去掉开头的T或者统一转为大写再做merge。轨道号本身也可能有前导零问题。比如017和17在Excel里显示不一样但程序里可能是同一个值。如果你在Excel里手工处理映射表Excel可能会自动把017变成17这时候再导回CSV做关联轨道号就对不上了。解决方法是在读取CSV时指定dtype把ORBIT列当作字符串读取。4.4 大数据量时的内存和性能优化如果做全国范围的图幅统计图幅网格有几千个要素研究区也有大量子要素sjoin虽然有效率但如果不注意也能跑半天。这里有两个实用优化技巧。第一提前对研究区做简化。很多研究区边界极其精细几万个顶点但做空间关系判断时根本用不到这么精细。可以先用shapely的simplify方法把容差设到0.01度能大幅减少计算量同时不会影响图幅关联结果。第二如果图幅网格特别大可以先按研究区的total_bounds做一个矩形框过滤只保留bounds与研究区bounds相交的图幅再进入sjoin阶段。这一下就能把几千个图幅过滤到几十个速度提升立竿见影。4.5 没有TILPAR文件时的备用方案如果实在找不到TILPAR文件也不是完全没法做。有一个粗略的近似办法利用图幅ID命名规律。同一轨道条带上的图幅通常在同一个UTM带内且列号字母相同行号字母连续变化。先按这个规则把所有图幅分成若干“列链”再结合卫星轨道方向去推断可能的轨道顺序。但这个方法有较大误差尤其在UTM带边缘和极高纬度区域图幅的列号会跳变所以只能作为应急方案不建议在正式项目里用。更推荐的做法是换一个思路很多在线数据目录平台其实会返回影像的orbit信息。你只需要检索某一时间范围内研究区经过的所有影像从返回结果里直接提取每个tile对应的轨道号然后和本地图幅网格做一次关联。这样虽然多了一步在线检索但拿到的轨道号是完全真实的比TILPAR文件更贴近实际数据分发情况。5. 最后再说点实操中的心得写了这么多最后分享几个我个人的使用习惯。首先是关于图幅统计结果的后续应用不要只停留在输出一张表。我在实际工作里经常把这个工具和后续的影像下载脚本串起来先做图幅统计得到图幅ID列表和轨道号然后直接用这个列表去批量检索或者下载对应的Sentinel-2影像。这样整个从“研究区确定”到“数据下载”的流程就闭环了中间不需要人工干预。其次是关于成果沉淀。我建议把统计表和矢量输出保存成带版本号的文件比如sentinel2_tile_stat_v20250101.xlsx因为研究区边界可能会随时调整边界一变统计结果就变了。保留历史版本可以追溯你最终用了哪一批图幅做分析写报告的时候能解释清楚为什么选了这些图幅而不是那些。最后想说轨道条带与图幅的关系在Sentinel-2数据体系里算是最基础的一环看起来不起眼但对下游所有工作都有影响。把这个逻辑理清楚后面无论是批量下载、云覆盖统计还是轨道内时序拼接都能少走很多弯路。你现在如果正好有一个研究区要处理不妨拿这套思路跑一遍相信会省下不少时间。