做遥感的人十有八九都被问过这样一类问题项目区在某某省、某某流域哨兵2号Sentinel-2到底走哪几条轨道覆盖哪些图幅以前我也和大家一样打开数据平台筛选界面一块一块把图幅编号复制出来心里还默念“千万别漏”。等到要批量做时序分析、按轨道组织数据时手动核对这种活儿真的能把人逼疯。后来我直接用PythonRS的思路把“轨道条带—图幅统计”做成了完整流程既能生成每个轨道覆盖的图幅列表也能输出表格和矢量文件。这篇文章就把这套东西拆开讲透从Sentinel-2的轨道条带原理、MGRS图幅组织方式到geopandas空间统计的完整代码再到我实测踩过的几个坑一次说清楚。1. 为什么要把“轨道条带—图幅关系”单独拿出来做统计1.1 业务痛点下载数据时最烦的不是算法而是图幅对齐很多人第一次接触哨兵2号数据时会发现一件挺反直觉的事平台是按Tile图幅分发数据的但很多业务却是按轨道来组织的。比如你要做某个区域的NDVI时间序列最理想的做法是锁定一条轨道把这条轨道每次过境的数据都拉下来这样光照角度、成像几何都相对一致时序曲线才平稳。可问题是一个研究区往往横跨好几个Tile而每个Tile又可能被不同轨道拍过到底哪些Tile属于你选的那条轨道没有底图的时候只能一个个去平台里试。另一个更常见的场景是批量下载。假设你要做全省范围的影像拼接先得知道这个省涉及多少Tile、每个Tile被几条轨道覆盖。如果你直接把所有Tile所有轨道的数据全下下来存储量至少翻两三倍但如果只按单一轨道取又可能在某几个Tile上缺数据。所以提前把“轨道—图幅”的映射关系统计出来是很多预处理流程的第一步。1.2 本文的产出物一张表加一个矢量我最终做的流程输出两样东西。第一是表格包括“轨道号→覆盖图幅列表”“图幅号→可用轨道列表”“研究区→涉及图幅与轨道”三类统计结果分别存成CSV和Excel。第二是矢量把图幅边界和轨道信息融合后的图层导成GeoPackage或GeoJSON方便直接在QGIS、ArcGIS里看也能作为后续裁剪、检索、任务分发的空间索引。这套流程不需要联网不依赖商业平台只要本地有两个基础矢量MGRS Tile网格和轨道footprint条带范围剩下用Python的geopandas就能算完。对遥感工程师来说是刚需对刚入门做哨兵2号数据处理的研究生来说也能少走很多弯路。2. Sentinel-2的轨道条带、MGRS图幅到底是怎么对应的2.1 轨道号就是条带的身份证先理一个概念我们说“轨道条带”在哨兵2号的语境里通常指一个相对轨道号Relative Orbit对应的一整条影像覆盖区域。Sentinel-2是太阳同步轨道双星A星B星组网后重访周期5天单星回归周期10天。每个回归周期内卫星在相邻轨道圈之间偏移一定距离从而形成一个固定的轨道网格。这些轨道圈按顺序编号就是相对轨道号范围一般在1到143之间不同平台和资料的起始编号可能略有差异。为什么要强调“相对轨道号”因为绝对轨道号Absolute Orbit会随着卫星发射时间、轨道维持不断累加今天是30000多号明天是30020号不适合做业务标识。而相对轨道号是固定在回归周期内的比如编号37的轨道每次回归都会经过大致相同的地面轨迹位置漂移很小。所以下载数据时平台的筛选条件里有一个“相对轨道号”这个编号就相当于轨道条带的身份证。条带的宽度大约是290公里。这个数值听起来不宽但卫星从北极飞到南极每一条带纵向上能拉出上万公里横向上覆盖接近三个100公里宽度的MGRS网格。所以一条轨道完整跑下来会穿过非常多图幅从高纬到低纬跨好几个UTM投影带。2.2 MGRS网格与产品Tile的关系哨兵2号的产品切分方式用一句话说就是沿轨方向按轨道条带成像横轨方向按MGRS格网分幅。MGRSMilitary Grid Reference System大家可能看着陌生但UTM投影大家都熟。MGRS把全球按UTM分带再在每一个投影带内划出100公里见方的网格每个网格用字母数字组合唯一编码。哨兵2号官方产品就是以这个100公里网格为单元切块的一个网格就是一景Tile。Tile的编号格式需要特别留意因为它存在两种写法。完整写法是T加上5位编码比如T50TLH或T50RQP不少平台和企业内部数据目录里会把开头的T省略直接写成50TLH。代码里处理的时候一定要统一别把带T和不带T混在一起做字符串匹配我见过太多因为这种小问题导致joins不上数据的案例。2.3 条带和Tile的交叠逻辑因为条带宽290公里而单个Tile宽度约100公里所以任意一条轨道过境时横向上至少会切过两到三个Tile。再加上轨道方向与UTM网格不是平行关系而是斜着穿过所以实际交叠情况更复杂一条轨道可能会覆盖同一UTM带内相邻的多个Tile也可能跨到相邻UTM带里。反过来一个Tile也不只被一条轨道覆盖。由于相邻轨道之间有重叠区域很多Tile会同时出现在两条甚至三条轨道的数据里。这就是为什么做时序分析之前必须先搞清楚你只管某一条轨道这个Tile到底有没有数据你希望某个Tile尽量多时相覆盖该选哪条轨道本质上这就是一个空间求交问题交给地理计算来处理比在网页上人工点选靠谱得多。3. 数据准备网格矢量与轨道footprint从哪来3.1 MGRS Tile网格数据要做统计第一步得有全球或区域范围的Tile边界矢量。我常用的获得渠道有三个。第一个渠道是ESA官方发布过Sentinel-2的Tile网格数据通常以KMZ、GeoJSON或Shapefile的形式包含在产品说明页和部分处理软件的资源包里字段里直接带TileID。第二个渠道是开源的sentinelhub相关的工具库或网站它们会把S2 Tile网格做成GeoJSON提供下载精度足够业务使用。第三个渠道是从自己已经下载的L1C或L2A产品里提取footprint打开压缩包里的MTD_MSIL1C.xml找到Product_Footprint节点下的坐标序列解析成多边形再按TileID合并成完整的网格。这个方法不依赖任何外部数据源适合手里已经有大量数据、想自己重建一份区域Tile网格的场景。需要注意不管从哪个渠道拿到的数据拿到之后先打印字段名和几何类型。不同来源的字段命名千差万别有的叫Name有的叫TileID有的直接叫MGRS_TILE代码里计算前统一改列名否则后面merge很容易出问题。3.2 轨道footprint数据轨道footprint就是每一条相对轨道在地面的覆盖范围。这个数据相对难找一点但也不是没有。最正规的渠道是ESA的Sentinel-2任务性能中心或哥白尼数据空间生态系统的元数据接口里面会提供轨道条带的范围信息有的版本是KML有的版本可以在API返回的JSON里找到footprint字段。直接把KML读进QGIS再另存成GeoPackage就完成准备了。如果网络渠道受限还有一个土办法用本地已有数据的Metadata重建。每景L1C产品里不仅写了自己这一景的footprint还会写明它的相对轨道号。把同一轨道号的所有产品footprint按轨道号union起来就得到该轨道在你研究区范围内的实际覆盖范围。这个方法不需要额外下轨道数据缺点是只能得到你“已经下载过”的数据覆盖范围没法覆盖没下载的轨道。但对很多区域级项目来说先知道自己有什么再去补什么已经够用了。3.3 Python环境配置环境方面我建议直接用一个干净的conda环境Python 3.9以上就行。核心依赖是geopandas它会把shapely、fiona、pyproj这些底层库一起带进来不需要单独一个个装。如果还要输出Excel需要加装openpyxl。conda create -n rs_env python3.9 conda activate rs_env pip install geopandas pandas openpyxl matplotlib这里有个小提醒geopandas的版本迭代比较快0.12之后sjoin的op参数改成了predicate网上很多旧教程里写的opintersects在最新版里会报警告。我下面的代码统一用predicate新老版本只要不是太老都能跑。4. 核心代码空间连接完成图幅统计并导出表格和矢量4.1 读取数据并统一坐标系准备阶段假设你手上有两个文件s2_tiles.gpkg存放MGRS Tile网格字段里有TileIDs2_orbits.gpkg存放轨道footprint字段里有Orbit_ID。先把它们读进来看一眼字段和坐标系。import geopandas as gpd import pandas as pd tiles gpd.read_file(s2_tiles.gpkg) orbits gpd.read_file(s2_orbits.gpkg) print(tiles.head()) print(orbits.head()) print(tiles crs:, tiles.crs) print(orbits crs:, orbits.crs) print(tiles columns:, tiles.columns.tolist()) print(orbits columns:, orbits.columns.tolist())实际项目里Tile网格可能是EPSG:4326经纬度也可能是某一种UTM投影轨道footprint大概率是WGS84经纬度。做空间连接前最好统一到同一个坐标系我习惯统一到EPSG:4326。如果只是求“是否相交”4326的变形影响不大效率也高如果后续要精确计算覆盖面积比例那就得临时再转成研究区对应的等面积投影。tiles tiles.to_crs(epsg4326) orbits orbits.to_crs(epsg4326)4.2 按轨道统计图幅一条轨道覆盖了哪些Tile这是最核心的需求。方法很直接用gpd.sjoin把Tile和轨道做空间连接凡是和该轨道footprint相交的Tile都会被匹配上然后按轨道分组汇总。matched gpd.sjoin(tiles, orbits, howinner, predicateintersects) matched matched.drop_duplicates(subset[TileID, Orbit_ID]) summary matched.groupby(Orbit_ID).agg( tile_count(TileID, count), tile_ids(TileID, lambda x: 、.join(sorted(x))) ).reset_index() print(summary.head())输出大致长这样Orbit_IDtile_counttile_ids63849RGP、49RHQ、49RHP……74149REQ、49RER、49RES……374350TLH、50TLJ、50TLK……sjoin这一步是空间计算的主战场它对每个Tile几何去和轨道几何盘算是否相交内部用了空间索引比写双重循环快几个量级。外层加一个drop_duplicates是必要的因为一个Tile可能和轨道footprint有多个接触片段不排重的话统计数量会成倍虚高。如果你想只看某一条轨道比如看37号轨道加一行筛选orbit37 summary[summary[Orbit_ID] 37] print(orbit37[[tile_count, tile_ids]])4.3 按图幅反查轨道一个Tile能被哪些轨道拍过反查逻辑其实就是把4.2的结果透视一下。同一份matched表把索引键倒过来。tile_to_orbits matched.groupby(TileID)[Orbit_ID].apply( lambda x: sorted(x.unique()) ).reset_index() tile_to_orbits[Orbit_List] tile_to_orbits[Orbit_ID].apply( lambda x: 、.join(map(str, x)) )这个表在数据缺失排查时尤其好用。比如你发现某个Tile在某段时间内没有数据先看看这个Tile在不在Orbit_List里如果压根没有轨道覆盖记录那是网格边界问题如果有轨道记录但平台上下不到数据那可能是数据覆盖或云量问题排查方向马上就不一样了。4.4 按研究区范围统计给定一个shp列出该区域的图幅和轨道前面两种统计面向全球或大区域但日常项目里更常见的是“我有一个研究区想知道需要准备哪些Tile和轨道”。实现上先读入研究区边界做一次sjoin把研究区范围内的Tile筛出来再拿这批Tile去和轨道做关联。study gpd.read_file(study_area.shp).to_crs(epsg4326) tiles_in_area gpd.sjoin(tiles, study, howinner, predicateintersects) tiles_in_area tiles_in_area.drop_duplicates(subset[TileID]) area_stats gpd.sjoin(tiles_in_area, orbits, howinner, predicateintersects) area_stats area_stats.drop_duplicates(subset[TileID, Orbit_ID]) print(area_stats[[TileID, Orbit_ID]].sort_values([TileID, Orbit_ID]))这时候你会发现研究区横跨的投影带一目了然TileID里的UTM带号直接写在编号前两位。比如研究区西部落在49带东部落在50带对应数据在后续几何处理时就必须分带投影不能拿一个固定UTM带从头算到尾。4.5 表格与矢量的导出统计结果不落盘等于白干。导出分两步。第一步导出纯表格CSV用utf-8-sig编码这是为了Excel打开不出乱码Excel版直接to_excel。summary.to_csv(orbit_tile_summary.csv, indexFalse, encodingutf-8-sig) tile_to_orbits.to_excel(tile_orbit_lookup.xlsx, indexFalse) with pd.ExcelWriter(sentinel2_area_stats.xlsx) as writer: summary.to_excel(writer, sheet_name轨道覆盖图幅, indexFalse) tile_to_orbits.to_excel(writer, sheet_name图幅可用轨道, indexFalse) area_stats[[TileID, Orbit_ID]].to_excel(writer, sheet_name研究区统计, indexFalse)第二步导出矢量。把轨道信息和Tile边界合并让每个Tile属性里带上能覆盖它的轨道列表。这里建议优先用GeoPackage或GeoJSON不要用Shapefile省去字段名截断和中文乱码的麻烦。tiles_final tiles[[TileID, geometry]].copy() tiles_final tiles_final.merge(tile_to_orbits[[TileID, Orbit_List]], onTileID, howleft) tiles_final.to_file(tiles_with_orbits.gpkg, layertiles, driverGPKG) # 顺便导出一个研究区范围的矢量 area_tiles tiles_final[tiles_final[TileID].isin(tiles_in_area[TileID])] area_tiles.to_file(area_tiles_with_orbits.geojson, driverGeoJSON)如果想把“每条轨道覆盖了哪些Tile”也用矢量表达出来可以把matched表按轨道号dissolve合并几何生成条带级矢量。这一步计算量稍大但结果很直观适合做汇报图。orbit_geom matched[[Orbit_ID, geometry]].dissolve(byOrbit_ID) orbit_geom.to_file(orbit_tiles_dissolve.gpkg, driverGPKG)5. 实测中容易翻车的地方以及效率优化建议5.1 坐标系不一致结果错得离谱我最开始在项目里踩过最大的坑就是坐标系没对齐傻傻地跑sjoin。Tile网格是某个UTM投影坐标轨道footprint是经纬度两套数据画在同一张图上一看都对但一做相交运算要么匹配不上要么匹配出完全错误的结果。geopandas在高版本里面对crs不一致的数据会主动警告甚至报错但如果是老版本它可能直接按数字坐标硬算结果完全不可信。所以代码里统一坐标那一步永远不要省。另外还要提醒一点只是打印crs相同还不够要确认它们都确切是4326的经纬度而不是投影坐标冒充的。最稳妥的办法是输出total_bounds看一眼坐标范围如果x和y都在-180到180、-90到90之间才放心做相交。5.2 intersects和within的取舍不要把擦边的图幅也算进来轨道footprint边缘会斜着擦过不少Tile的角。用predicateintersects时被擦到一小块的Tile也会被统计进来。对“这Tile到底有没有这条轨道的数据”这个问题来说intersects是合理的因为只要擦到一点就有数据。但对“哪些数据是整齐覆盖研究区的”来说这种擦边图幅常常是干扰项下载下来发现只有几行像素有用。更科学的做法是加一个面积占比阈值。比如只保留覆盖面积超过该Tile面积5%的图幅记录减少下游的无效数据量。实现上需要用到overlay求交集面积计算量比sjoin大所以我的做法是先用sjoin粗筛出候选Tile再对候选集做精确面积计算。from shapely.geometry import shape candidates gpd.sjoin(tiles, orbits, howinner, predicateintersects) candidates candidates.drop_duplicates(subset[TileID, Orbit_ID]) overlay_res gpd.overlay(candidates, orbits[[Orbit_ID, geometry]], howintersection) overlay_res[inter_area] overlay_res.geometry.area tile_area candidates[[TileID, geometry]].copy() tile_area[tile_area] tile_area.geometry.area overlay_res overlay_res.merge(tile_area[[TileID, tile_area]], onTileID) overlay_res[cover_ratio] overlay_res[inter_area] / overlay_res[tile_area] valid overlay_res[overlay_res[cover_ratio] 0.05]这里把阈值的判断逻辑写出来了实际用多少需要根据业务定。做无缝拼接时可以放宽到0.01做单轨道时序分析时建议收紧到0.1以上宁可少不能滥。5.3 性能优化少用overlay多用sjoin全球的MGRS Tile数量约五千块上下轨道条带143条理论上笛卡尔积也就是七十万级的判断sjoin很快。但如果你一上来就gpd.overlay(tiles, orbits)让软件把所有Tile裁剪成轨道交集碎片数据量会大好几倍跑起来明显卡。overlay的结果是“重合区几何”sjoin的结果是“属性关联关系”。大部分图幅统计场景只需要属性关联完全没必要去算重合区多边形。实际项目里我一般先用sjoin拿到关系列表再根据需要去重、分组、统计。只有在最后做覆盖面积占比时才对筛选出来的小范围数据做overlay这样速度能差出一个量级。5.4 输出字段与编码问题导出CSV时encodingutf-8-sig这一步很多人会漏。漏掉的后果就是CSV在Windows上用Excel打开直接乱码只能再用记事本转码非常折腾。GeoJSON默认utf-8没有这个问题建议矢量优先导出GeoPackage和GeoJSON。如果用Shapefile导出注意三个坑一是字段名会被截断成10个字符Orbit_List这种名字完全没问题但tile_count_with_orbit_detail这种长字段会被截得乱七八糟二是中文属性值容易乱码除非指定编码否则不建议在SHP里放中文三是轨道号是整数时属性表里排序会按数值来但如果轨道号前面补了0它可能变成字符串排序就乱套了。我习惯于把轨道号直接用整数存储只在拼接字符串展示时才转成字符。6. 一点个人经验以及后续还能怎么玩最后分享几个我实际做项目时的习惯。第一不要凭记忆背“某某区域是哪几条轨道”。轨道号看着稳定但不同平台、不同版本产品里的编号起始位不一样用代码去算永远比脑记可靠。拿到需求第一件事就是跑一遍sjoin用数据说话。第二把这个统计表做成一个缓存文件。轨道—图幅映射关系在短时间内不会剧烈变化做完一次统计后把结果存成CSV和GeoPackage下次做单景数据筛选时直接查表不必每次重新做空间计算。我自己的项目里通常只每季度更新一次轨道footprint全球Tile网格更是可以一直复用。第三把这张表和下载链路串起来。比如通过OData接口批量检索时把统计得出的TileID列表和Orbit_ID作为筛选条件能非常精准地把数据范围框定在研究区相关记录内避免“全下再筛”的低效存储。如果再把云量阈值加进筛选基本就是一个小型自动化数据准备系统。再往深了做这套统计还能延伸两个方向。一个是结合研究区矢量做“按行政边界拆分轨道覆盖面积占比”输出每个行政区涉及的轨道清单和面积比例这个在项目汇报时很有说服力。另一个是做“多时相覆盖度评估”把轨道回归周期叠加上去估算某个Tile在一段时间内能拿到多少有效观测次数这在植被物候、土壤水分监测里很实用。我也是踩过不少坑才把这条流程理顺的。手动核对时代复制粘贴出错率高自动化统计后虽然也需要检查数据但至少把机械劳动交给了代码剩下的精力能花在更值得处理的遥感问题上。希望这套“轨道条带—图幅统计”的流程也能帮你省下半天时间。