简介这份代码包围绕MCD19A2遥感数据处理展开面向具备一定IDL与Python基础、从事地球科学、环境监测及植被分析的科研人员与学习者重点解决批量预处理、影像融合、镶嵌、重采样、裁剪和栅格点值提取等流程的代码实现问题。压缩包共12个文件约11KB以3个py脚本和3个pro程序为主体分别承担重投影、裁剪、AOD栅格提取点值以及MCTK调用、均值统计、批量镶嵌等任务另含xml、iml等IDE项目配置便于在对应环境中查看与复用。代码覆盖MCD19A2从影像组织到专题值提取的关键环节可帮助读者理解MCTK批量处理的脚本结构、参数设置与流程衔接并据此改造自己的自动化处理方案对刚接触该产品、需要快速搭建批量流程的读者尤其友好。已有1248人学习适合作为遥感数据批量处理的学习参考与代码模板。1. 为什么一份气溶胶数据的处理代码值得单独拆开看第一次拿到 MCD19A2 的人八成会卡在打开文件这一步HDF4 容器、二维正弦投影、带比例因子的整型子数据集、藏在一串二进制位里的质量标记四件事叠在一起能直接跑通的完整代码反而不多。这份资源就是围绕 MCD19A2 数据处理过程整理的一组代码覆盖读取、质量筛选、重投影、裁剪镶嵌到时间聚合的整条链路。它不是那种装完就能点的软件/插件而是一批可读写、可改参数、可挂到批处理任务上的脚本适合要批量出 AOD 栅格的研究生、做大气环境分析的工程师以及手头有一堆 HDF 文件却不知道怎么变成 GeoTIFF 的人。2. 环境与数据准备先把 HDF4 容器和正弦网格摸清楚2.1 依赖装不上的那半天pyhdf 与 GDAL 的编译依赖这套代码要跑起来绕不开两个重量级依赖读 HDF4 需要pyhdf重投影和写出需要gdal。这两个包最恶心的地方在于它们都不是纯 Python底层各自带着 C 库。用pip install pyhdf直接装在没有预编译轮子的机器上大概率会卡在编译阶段报错提示找不到hdf.h。血泪经验是先装系统的 HDF4 开发库再装 Python 包。# 用 conda 装是最省事的路径它会自动把底层 C 库一起解决 conda create -n maiac python3.10 conda activate maiac conda install -c conda-forge pyhdf gdal rasterio pyproj xarray numpy -y # 装完先验一下 GDAL 有没有把 HDF4 驱动编进去 python -c from osgeo import gdal; print(gdal.GetDriverByName(HDF4) is not None)上面这段命令做了三件事建环境、装依赖、验证 GDAL 是否带 HDF4 驱动。最后一行打印True才算过关打印False说明这个 GDAL 是精简编译的读不了 MCD19A2得换一个 conda-forge 的包。这个检查我一般放在所有脚本之前因为驱动缺失的表现是文件明明存在gdal.Open却返回 None排查起来特别费时间。第二条命令里pyhdf和gdal装在同一个 conda 环境里是为了避免两边链接到不同版本的 HDF4 库导致段错误。如果团队里有人用 pip 混装建议让所有人都用同一份环境文件别各装各的。2.2 一个文件里到底装了什么SDS 清单与属性字段MCD19A2 是网格化的二级产品每个 HDF4 文件里装着十几个子数据集SDS常用的其实就三组气溶胶光学厚度Optical_Depth_047、Optical_Depth_055、质量标记AOD_QA、观测几何Solar_Zenith、View_Zenith、Relative_Azimuth。这些 SDS 全是整型真实物理量要靠scale_factor和add_offset反算无效像元用_FillValue填充。在动手写正式代码前先跑一段探测脚本把文件结构打出来比翻文档快得多from pyhdf.SD import SD, SDC path MCD19A2.A2020150.h26v05.006.2020158031655.hdf hdf SD(path, SDC.READ) # datasets() 返回 {SDS 名称: (维度, 类型, 维度索引)} for name, (dims, dtype, idx) in hdf.datasets().items(): sds hdf.select(name) attrs sds.attributes() print(f{name:24s} dims{dims} dtype{dtype}) # 只打印和数值换算相关的关键属性属性全打会刷屏 for key in (scale_factor, add_offset, _FillValue, valid_range): if key in attrs: print(f {key} {attrs[key]}) sds.endaccess() hdf.end()第一个datasets()给出的是全部 SDS 的元信息维度一般是(1200, 1200)也就是单块瓦片 1200 行 1200 列。注意dtype字段气溶胶厚度通常是int16QA 可能是 8 位也可能是 16 位不同批次的数据集版本会不一样所以代码里不要写死位宽先打印出来再决定掩码怎么写。第二段attributes()是重点。scale_factor决定了原始整数乘多少才是物理量气溶胶一般都是 0.001也就是原始值 1000 对应 AOD 1.0_FillValue常见是 -28672所有等于这个值的像元都必须先剔除再做运算否则参与平均会把结果拉得一塌糊涂。valid_range则是产品自己声明的有效区间AOD 允许轻微负值代表低于探测限的洁净大气所以别一上来就把负数全砍掉。还有一点每次select之后必须endaccess()文件句柄会累积批量处理几千个文件时不释放很快就会撞上系统句柄上限。这个细节在单文件测试里完全看不出来一上批量就翻车。2.3 目录约定与配置入口让批处理脚本认识你的数据盘原始数据一般按年/日序或者瓦片号分层存放脚本如果每换一批数据就改一次路径维护成本会越来越高。比较常见的做法是抽一个 JSON 配置文件出来把所有可变的东西集中到一个地方{ input_dir: /data/maiac/2020, output_dir: /data/out/aod_550, boundary_shp: /data/shp/region.shp, dataset_name: Optical_Depth_055, qa_datasets: [AOD_QA], qa_accept_bits: [0, 1], scale: 0.001, fill_value: -28672, valid_min: -100, valid_max: 5000, out_nodata: -9999.0, resampling: bilinear, target_epsg: 4326, max_workers: 8, resume_log: /data/out/processed.jsonl }字段分成三类看input_dir/output_dir/boundary_shp是路径dataset_name、scale、fill_value、valid_*是数据解码参数换了产品版本或者换了波段就要动这里resampling、target_epsg、max_workers是处理策略。把这些全部外置的好处是同一份代码既能出 550nm 也能出 470nm改一个字段就行不用去翻源码。配置里特意留了resume_log它是断点续跑的依据后面第 4 章会用到。目录上建议输出按年/月再分一层单目录塞几万个 GeoTIFFWindows 上资源管理器会直接卡死Linux 上ls也要等半天。3. 核心读取链路从 HDF4 子数据集到带地理信息的数组3.1 读一个 SDS 的完整动作填充值、比例因子与无效值读取这一步看起来简单实际上顺序错了结果就全废。正确的顺序是先取原始整型数组 → 剔除填充值 → 转浮点 → 乘比例因子 → 用valid_range之外的再筛一遍。反过来先乘比例因子再比较填充值会因为浮点精度问题漏掉一部分无效像元。import numpy as np from pyhdf.SD import SD, SDC def read_sds(path, name, scale0.001, fill-28672, vmin-100, vmax5000): 读取单个 SDS返回 float32 数组无效像元为 NaN。 hdf SD(path, SDC.READ) try: sds hdf.select(name) raw sds.get() # 原始整型形状 (1200, 1200) attrs sds.attributes() # 属性里有的以属性为准没有就用传入的默认值 scale attrs.get(scale_factor, scale) fill attrs.get(_FillValue, fill) sds.endaccess() finally: hdf.end() arr raw.astype(np.float32) # 第一步填充值置 NaN必须在乘比例之前做 arr[raw fill] np.nan # 第二步整型原始值按 valid_range 过滤超出即无效 arr[(raw vmin) | (raw vmax)] np.nan # 第三步换算成物理量 arr arr * float(scale) return arrread_sds里有几个点值得说。sds.get()拿到的是完整的二维数组1200×1200 的int16大约 2.8 MB一次读进内存没问题但如果是镶嵌后的更大范围就得考虑分块读pyhdf支持按行列切片读取可以后面再优化。attrs.get(...)这种写法是给不同批次数据留后路。同一个产品不同版本scale_factor偶尔会变写死之后换数据就静默出错——数值看着正常量级差 10 倍这种错误最难发现因为图还是能画出来的。用属性优先、参数兜底的方式两边的坑都能避开。第三步的量纲换算放在最后是为了让前两步的比较都在整型域完成语义最清晰。返回float32而不是float64是因为后面要做多年合成内存占用差一倍float32对 AOD 这种精度完全够用。还有个小细节try/finally保证异常时也关闭文件。批处理跑几千个文件中间只要有一个损坏文件抛异常没关句柄跑到后面就会莫名其妙地失败。3.2 AOD_QA 位运算把一团数字拆成三档可信度AOD 的 QA 不是好/坏两档而是一串二进制位不同位段分别表示云掩膜邻域、地表类型、反演质量等级、耀斑影响等。直接把 QA 当数值比大小比如qa 2是典型误用因为数值大小和位段含义完全不对应。稳妥的做法是用位掩码逐段提取import numpy as np def decode_aod_qa(qa, quality_shift5, quality_mask0b111, accept_levels(2, 3)): 从 AOD_QA 中提取反演质量等级返回布尔掩膜。 qa qa.astype(np.uint16) # 右移到目标位段起始位置再与掩码取出该段的原始值 level (qa quality_shift) quality_mask # 只保留质量等级在可接受集合内的像元 good np.isin(level, accept_levels) # 顺带把未反演level0的像元单独标出来便于统计 return good, levelquality_shift和quality_mask是可调参数不同批次产品对 QA 位段的定义可能不一样改之前一定要对着产品文档的位定义表核一遍别照抄别人博客里的数字。这也是我把 shift 和 mask 做成函数参数而不是写死常量的原因。np.isin那一步等价于只接受质量等级为 2 和 3 的像元如果你想更宽松把accept_levels改成(1, 2, 3)即可。返回的level数组别扔掉后面做质量统计、画 QA 分布图都用得上比重新解一遍省事。实际用的时候通常会再叠一层几何筛选太阳天顶角过大低太阳高度角时反演误差明显上升加一条solar_zenith 70的条件能砍掉不少边缘数据。几何 SDS 的比例因子一般是 0.01读法跟 AOD 一样。最终掩膜是质量合格 几何合格 非 NaN三者取与。这一步做完一整块瓦片大概只剩 30%60% 的有效像元看着心疼但这是让结果能用的前提。3.3 正弦投影到 WGS84自己算仿射参数还是交给 GDALMCD19A2 用的是正弦投影Sinusoidal坐标单位是米原点在赤道与本初子午线交点。每个瓦片覆盖固定的地理范围只要知道瓦片号和网格尺寸就能手算出左上角坐标和像元大小import numpy as np R 6371007.181 # 正弦投影球体半径米 TILE 1111950.519667 # 单块瓦片边长米 GRID 1200 # 每块瓦片行列数 PIXEL TILE / GRID # 像元边长约 926.625 米 X_ORIGIN -20015109.354 # 全球网格最左边界 Y_ORIGIN 10007554.677 # 全球网格最上边界 def tile_transform(h, v): 由瓦片行列号推出 GDAL 风格的仿射六参数。 x_min X_ORIGIN h * TILE y_max Y_ORIGIN - v * TILE return (PIXEL, 0.0, x_min, 0.0, -PIXEL, y_max)tile_transform返回的六个数字就是 GDAL 的 GeoTransform(x方向像元宽, 旋转, 左上角x, 旋转, y方向像元高(负), 左上角y)。h和v从文件名里解析出来第 4 章会给正则。注意 y 方向的像元高是负数因为栅格行号向下增长而地理坐标向上增长符号写错的结果就是图像上下颠倒但看上去还挺像那么回事属于隐蔽性很强的错误。手算仿射参数的好处是不依赖外部库速度快适合先落一个中间 GeoTIFF 再统一重投影的两段式流程。缺点是只对标准网格有效如果数据被裁过或者重采样过公式就对不上了。另一种做法是直接用 GDAL 一步到位from osgeo import gdal, osr import numpy as np def warp_to_wgs84(src_path, dst_path, arr, transform, epsg4326): 把内存数组写成带正弦投影的临时栅格再重投影到 WGS84。 driver gdal.GetDriverByName(GTiff) tmp /vsimem/tmp_sinu.tif # 内存文件不落磁盘 ds driver.Create(tmp, arr.shape[1], arr.shape[0], 1, gdal.GDT_Float32) ds.SetGeoTransform(transform) srs osr.SpatialReference() srs.ImportFromProj4( fprojsinu lon_00 x_00 y_00 a{R} b{R} unitsm no_defs) ds.SetProjection(srs.ExportToWkt()) band ds.GetRasterBand(1) band.SetNoDataValue(np.nan) band.WriteArray(arr) band.FlushCache() ds None gdal.Warp(dst_path, tmp, dstSRSfEPSG:{epsg}, resampleAlgbilinear, srcNodatanp.nan, dstNodatanp.nan, multithreadTrue)/vsimem/是 GDAL 的内存文件系统临时栅格不落磁盘批量处理时能省掉大量 IO。SetNoDataValue(np.nan)配合srcNodata/dstNodata保证无效像元在重采样时不会污染邻域——这一点很关键双线性插值会把 NaN 周围的有效值一起拉低。multithreadTrue在大范围镶嵌时提升明显但单块 1200×1200 的瓦片上收益有限可以按需关掉。两段式先手算仿射写中间文件再gdal.Warp和一步式内存数组直接 warp各有场景前者适合需要保留中间结果的调试后者适合流水线里只想拿到最终产物。我在批量任务里用的是后者因为磁盘上多出来的中间文件清理起来很烦。4. 批处理流水线把单文件脚本变成能跑一整年的任务4.1 从文件名解析瓦片号与日期MCD19A2 的文件名是高度结构化的年份、日序、瓦片行列号、版本、处理时间全在里面。与其一个个手动填参数不如用正则一次性抽出来import re from datetime import datetime, timedelta PATTERN re.compile( r^MCD19A2\.A(?Pyear\d{4})(?Pdoy\d{3}) r\.h(?Ph\d{2})v(?Pv\d{2}), re.IGNORECASE) def parse_name(path): 从文件名解析日期与瓦片号解析失败返回 None。 name path.rsplit(/, 1)[-1] m PATTERN.match(name) if not m: return None year, doy int(m.group(year)), int(m.group(doy)) date datetime(year, 1, 1) timedelta(daysdoy - 1) return { date: date, year: year, month: date.month, h: int(m.group(h)), v: int(m.group(v)), tile: fh{m.group(h)}v{m.group(v)}, }日序转日期那行要注意datetime(year, 1, 1) timedelta(daysdoy - 1)减 1 是因为 1 月 1 日对应日序 001。忘了减 1 会导致全年日期整体偏移一天做时间序列分析时和地面站点对不上查起来得从头翻。正则只匹配到瓦片号就停不去匹配版本和处理时间因为这两个字段在后续流程里用不上匹配得越长越容易因为格式差异失配。返回None而不是抛异常是为了让批处理能跳过不符合命名的杂文件比如产品附带的元数据文件。解析出来的tile直接用于 3.3 的tile_transformdate用于输出路径分年和时间聚合的分组键。4.2 多进程 断点续跑中断之后不用从头再来单文件处理大概要几百毫秒到一两秒一年的全国范围数据动辄几千个文件串行跑要几个小时。用进程池并发是常规操作但真正让批处理能长期用的是断点续跑——跑到一半机器重启、磁盘满、被踢下线重新跑的时候不用从头开始。import json, os from concurrent.futures import ProcessPoolExecutor, as_completed def load_done(log_path): 读取已完成记录返回已处理文件路径集合。 done set() if os.path.exists(log_path): with open(log_path, r, encodingutf-8) as f: for line in f: line line.strip() if not line: continue try: done.add(json.loads(line)[src]) except (json.JSONDecodeError, KeyError): continue # 中断时最后一行可能是半截直接跳过 return done def run_batch(files, worker, log_path, max_workers8): done load_done(log_path) todo [f for f in files if f not in done] print(ftotal{len(files)} done{len(done)} todo{len(todo)}) with ProcessPoolExecutor(max_workersmax_workers) as pool, \ open(log_path, a, encodingutf-8) as log: futures {pool.submit(worker, f): f for f in todo} for fut in as_completed(futures): src futures[fut] try: out fut.result() log.write(json.dumps({src: src, out: out}) \n) log.flush() # 每条立刻落盘进程被杀也不丢 except Exception as e: log.write(json.dumps({src: src, error: str(e)}) \n) log.flush()load_done里对 JSON 解析失败的行直接跳过是因为进程被强杀时最后一行大概率只写了半截。如果这里抛异常整批任务就再也启动不了了非常要命。log.flush()那两行是整个断点续跑的关键。Python 的文件缓冲默认几 KB不主动 flush 的话程序被kill -9时缓冲区里的记录全丢重跑时又得把已经做完的文件再做一遍。加上 flush 之后最多损失正在处理的那几个文件。max_workers不要设得太大。每个进程要加载 GDAL 和 numpy内存占用不低8 个进程大概吃掉 34 GB普通工作站开 16 个很容易触发内存交换速度反而下降。一般取 CPU 物理核数的一半左右比较稳。日志用的是 JSON Lines 格式一行一条记录追加写不用读全文件就能续跑比 CSV 更适合这种场景。4.3 镶嵌、按边界裁剪与时间维度聚合单景处理完之后常见需求有三个把相邻瓦片拼成一整块、按研究区边界裁出来、按时间维度做合成。这三件事各自有坑。镶嵌用 GDAL 的Warp一次性传多个输入最省事from osgeo import gdal def mosaic_and_clip(inputs, out_path, shp_path, epsg4326): 多景镶嵌后用边界裁切一步出结果。 options gdal.WarpOptions( formatGTiff, dstSRSfEPSG:{epsg}, cutlineDSNameshp_path, # 直接按矢量边界裁 cropToCutlineTrue, # 输出范围紧贴边界不留空白 resampleAlgbilinear, srcNodatafloat(nan), dstNodata-9999.0, multithreadTrue, creationOptions[COMPRESSLZW, TILEDYES, BIGTIFFIF_SAFER], ) gdal.Warp(out_path, inputs, optionsoptions)cutlineDSName直接吃矢量文件省掉先镶嵌再裁剪两步。cropToCutlineTrue让输出范围紧贴边界外接框不保留整景范围文件能小很多。BIGTIFFIF_SAFER是给大范围数据留的后路超过 4 GB 的 GeoTIFF 不加这个参数会直接写失败报错信息还很不直观。COMPRESSLZW对浮点栅格压缩率一般但聊胜于无。时间聚合建议用xarray它对带时间维的数据立方体处理起来比手写循环清晰得多import xarray as xr import glob def monthly_composite(pattern, out_nc): 按天文件聚合成月均值输出 NetCDF。 files sorted(glob.glob(pattern)) ds xr.open_mfdataset( files, combineby_coords, preprocesslambda d: d.where(d -9000), # 剔除 nodata parallelTrue, ) monthly ds.resample(time1MS).mean(skipnaTrue) monthly.to_netcdf(out_nc, encoding{__xarray_dataarray_variable__: {zlib: True, complevel: 4}})preprocess里那行where(d -9000)是专门对付 nodata 的。如果输入文件的 nodata 设成了 -9999 而xarray没识别出来月均值会被拉低一大截图上看就是大片异常的暗区。这一步加上去之后就干净了。resample(time1MS)表示按月起算skipnaTrue保证有效像元少的日子不会把整月拖成 NaN。编码里的complevel4是在压缩率和速度之间取的折中设到 9 会慢好几倍。5. 避坑记录五个让我重跑整夜的排查现场5.1 读出来的数组全是 -28672现象脚本没报错输出的 GeoTIFF 打开后一片均匀值统计一下全是 -28672 或者 0。原因sds.get()返回的原始数组里无效像元就是-28672代码里只乘了比例因子没有先剔除填充值或者剔除用了np.isnan(raw)——原始数组是整型isnan判断永远返回 False。解决先确认原始数组的 dtype整型的话用raw fill_value做布尔索引不要用isnan。填充值从_FillValue属性里取取不到再用默认值。这个坑我在第一次跑批的时候踩过一整晚的机时全浪费在生成垃圾文件上从那以后每个新数据集我都先跑一遍小样本统计。5.2 重投影后图像整体偏移了半个像元现象把结果和研究区矢量叠在一起图像整体往左上角偏了差不多一个像元边缘尤其明显。原因犯的是 GeoTransform 的经典错误——把像元中心坐标当成了像元左上角坐标或者把正弦投影的球体半径写成了 WGS84 的椭球参数。MCD19A2 用的是球体a和b都应该是 6371007.181混用椭球参数会让整幅图产生系统性位移。解决检查ImportFromProj4里的a和b是否都写成了球半径检查 GeoTransform 的第三个和第六个参数是不是左上角外边界。验证方法很简单取一个已知经纬度的地面站点看它在栅格上的位置和图上的位置是否一致偏半个像元在图上肉眼就能看出来。5.3 多进程跑到一半内存爆掉现象批处理跑到几百个文件后系统开始疯狂交换速度骤降最后被 OOM Killer 干掉。原因一是max_workers开太大每个进程都加载了完整 GDAL 和 numpy基础内存就有几百 MB二是sds.endaccess()或者ds None漏了GDAL 数据集对象没释放内存只增不减三是用open_mfdataset一次性打开了上千个文件每个都保持句柄。解决把max_workers降到物理核数的一半用psutil监控每个进程的 RSS所有 HDF 和 GDAL 对象都放在try/finally里或者显式置 Noneopen_mfdataset加chunks{time: 30}做分块读取让xarray按需加载而不是全读进内存。5.4 镶嵌结果出现明显的接缝现象相邻两景拼起来之后交界处有一条颜色突变的线AOD 值在接缝两侧差 0.1 以上。原因接缝处的像元在两景里都是边缘像元重采样时双线性插值把 nodata 混了进来导致边缘值被拉低。另一种可能是两景的观测时间不同气溶胶本身就在变化属于真实差异而不是处理错误。解决先确认是不是 nodata 污染把srcNodata正确设上或者在镶嵌前对每景做一次边缘腐蚀把最外圈几个像元置为无效。如果确认是时间差异就改成分时间段镶嵌或者接受这个差异并在论文里说明。我一般会做一个接缝剖面图沿垂直于接缝的方向采样一眼就能看出是插值问题还是真实变化。5.5 时间聚合后数值整体偏低现象月均值和单日值对比均值明显偏低而且有效像元比例越高偏低越明显。原因xarray在做mean时如果没有正确识别 nodata会把-9999当成真实值参与平均另一种情况是老版本open_mfdataset的combineby_coords在时间坐标有微小差异时产生了重复时间点导致权重错乱。解决在preprocess里显式把 nodata 转成 NaN并且用ds.where(ds ! nodata)而不是ds.where(ds -9000)这种模糊判断聚合前打印一下ds.time的长度和文件数是否一致不一致说明合并出了问题。还有一个习惯聚合完成后随机抽一个像元把它所有日的原始值和均值手动算一遍对一下这招比看统计量管用得多。6. 二次开发把脚本改成一个能长期维护的处理库脚本跑通只是起点真正要长期用得把它拆成配置 处理类 命令行入口三层。我一般的做法是定义一个AodProcessor类把解码参数和处理策略都在__init__里接住process_one只负责单个文件run负责调度。这样换产品、换波段、换研究区都只动配置不动逻辑。class AodProcessor: def __init__(self, cfg: dict): self.cfg cfg self.cache {} # 缓存边界矢量避免每个进程重复读 def process_one(self, src): meta parse_name(src) if meta is None: return None aod read_sds(src, self.cfg[dataset_name], scaleself.cfg[scale], fillself.cfg[fill_value]) qa_ok, _ decode_aod_qa(self._read_qa(src)) valid np.isfinite(aod) qa_ok arr np.where(valid, aod, np.nan) transform tile_transform(meta[h], meta[v]) dst self._out_path(meta) warp_to_wgs84(src, dst, arr, transform, self.cfg[target_epsg]) return dst改造完之后验证正确性我固定走三步。第一步拿一个月的数据跑新旧两版代码逐像元做差最大绝对差应该为 0第二步抽三景结果和手算的单像元值对照确保量纲没变第三步把输出拖进 QGIS 和研究区矢量叠一下看边界有没有偏移。这三步做完基本能排除 90% 的改造事故。最容易翻车的地方是顺手优化。看到read_sds里有三个循环就想合并成一个看到float32想改成float64提精度结果要么改变了无效值的判断逻辑要么内存翻倍。我的血泪经验是重构只做结构搬运不改数值路径任何涉及计算顺序的改动都单独提一次用上面三步验一遍。还有一个实用技巧给处理类加一个dry_run开关只解析文件名、算输出路径、检查目标文件是否已存在不真正读数据。批量任务前面先跑一遍 dry run能在几分钟内发现路径拼错、目录没权限、瓦片号越界这类问题比跑到一半再报错省事太多。从那以后我每次改这套代码都强制先跑一遍 dry run 加一个月的小样本对照确认没问题再挂到全量任务上。希望帮到你。本文还有配套的精品资源点击获取
