diff --git a/src/mintpy/cli/load_data.py b/src/mintpy/cli/load_data.py index 7585fd9d7..d38481828 100755 --- a/src/mintpy/cli/load_data.py +++ b/src/mintpy/cli/load_data.py @@ -5,14 +5,21 @@ # Author: Antonio Valentino, Zhang Yunjun, Aug 2022 # ############################################################ +# os 用来处理文件路径，例如取当前脚本文件名、判断模板文件名。 import os +# sys 用来读取命令行参数和退出程序；sys.exit(0) 表示正常退出，sys.exit(1) 表示出错退出。 import sys +# auto_path 保存不同处理软件（ISCE/GAMMA 等）的默认输入路径规则。 from mintpy.defaults import auto_path +# get_template_content() 会读取 MintPy 内置的模板说明文本，显示在命令行帮助中。 from mintpy.defaults.template import get_template_content +# create_argument_parser() 是 MintPy 封装的命令行解析器创建函数。 from mintpy.utils.arg_utils import create_argument_parser ################################################################# +# DEFAULT_TEMPLATE 是 `load_data.py -H` 时显示的示例模板。 +# """...""" 是多行字符串；format(...) 把不同处理器的自动路径说明拼进去。 DEFAULT_TEMPLATE = """template: ########## 1. Load Data (--load to exit after this step) {}\n @@ -23,14 +30,17 @@ auto_path.AUTO_PATH_ISCE_TOPS, ) +# TEMPLATE 是 load_data 对应的完整模板配置说明。 TEMPLATE = get_template_content('load_data') +# NOTE 是补充说明，提醒用户哪些数据集是必需的、文件名里需要包含什么日期信息。 NOTE = """NOTE: For interferogram, unwrapPhase is required, the other dataset are optional, including coherence, connectComponent, wrapPhase, etc. The unwrapPhase metadata file requires DATE12 attribute in YYMMDD-YYMMDD format. All path of data file must contain the reference and secondary date, either in file name or folder name. """ +# EXAMPLE 是命令行示例，会显示在 `load_data.py -h` 的帮助信息里。 EXAMPLE = """example: # MUST run in the mintpy working directory! @@ -54,6 +64,7 @@ def create_parser(subparsers=None): """Create command line parser.""" + # 创建命令行解析器：把终端输入的 `load_data.py -t xxx.cfg` 转换成 Python 变量。 synopsis = 'Load stacks of interferograms to HDF5 files' epilog = TEMPLATE + '\n' + NOTE + '\n' + EXAMPLE name = __name__.split('.')[-1] @@ -61,18 +72,23 @@ def create_parser(subparsers=None): name, synopsis=synopsis, description=synopsis, epilog=epilog, subparsers=subparsers) # extra help + # action='store_true' 表示用户输入 -H 时，inps.print_example_template 变成 True。 parser.add_argument('-H', dest='print_example_template', action='store_true', help='Print/Show the example template file for loading.') # input files + # nargs='+' 表示 -t 后面可以跟一个或多个模板文件。 parser.add_argument('-t', '--template', dest='template_file', type=str, nargs='+', help='template file(s) with path info.') + # --geom 只加载几何文件，不加载干涉图栈；常用于先检查几何输入是否正确。 parser.add_argument('--geom','--geometry', dest='only_load_geometry', action='store_true', help='Load the geometry file(s) ONLY.') # options from template file name & content + # PROJECT_NAME 会写入输出 HDF5 元数据，便于后续产品识别。 parser.add_argument('--project', type=str, dest='PROJECT_NAME', help='project name of dataset for INSARMAPS Web Viewer') + # --enforce 会把 updateMode 设为 False，意思是强制重新加载，不用已有文件跳过。 parser.add_argument('--enforce', '-f', dest='updateMode', action='store_false', help='Disable the update mode, or skip checking dataset already loaded.') @@ -83,10 +99,12 @@ def cmd_line_parse(iargs=None): """Command line parser.""" # parse parser = create_parser() + # parse_args() 解析参数；iargs 为 None 时读取真实命令行，传列表时方便其它代码或测试调用。 inps = parser.parse_args(args=iargs) # check: -H option if inps.print_example_template: + # 只打印示例模板并退出，不执行真正的数据加载。 print(DEFAULT_TEMPLATE) sys.exit(0) @@ -94,6 +112,7 @@ def cmd_line_parse(iargs=None): # -t option is required AND # smallbaselineApp.cfg file is required if not inps.template_file: + # 如果没有 -t/--template，就打印用法并退出。 parser.print_usage() script_name = os.path.basename(__file__) print(f'{script_name}: error: -t/--template option is required.') @@ -101,10 +120,12 @@ def cmd_line_parse(iargs=None): sys.exit(1) elif all(not x.endswith('smallbaselineApp.cfg') for x in inps.template_file): + # load_data 至少需要默认模板 smallbaselineApp.cfg，因为很多默认路径和选项都来自它。 script_name = os.path.basename(__file__) print(f'{script_name}: error: at least smallbaselineApp.cfg file is required for -t/--template option.') sys.exit(1) + # 返回解析后的参数对象，后面的 main() 会把它交给核心 load_data() 函数。 return inps @@ -112,12 +133,15 @@ def cmd_line_parse(iargs=None): ################################################################# def main(iargs=None): # parse + # main() 是 CLI 入口函数：先解析命令行参数。 inps = cmd_line_parse(iargs) # import + # 延迟导入核心函数，只有真正要运行时才加载 src/mintpy/load_data.py。 from mintpy.load_data import load_data # run + # 调用核心模块的 load_data(inps)，把外部数据读入 MintPy 标准 HDF5 文件。 load_data(inps) diff --git a/src/mintpy/cli/smallbaselineApp.py b/src/mintpy/cli/smallbaselineApp.py index c0de95c0e..0c91eaadd 100755 --- a/src/mintpy/cli/smallbaselineApp.py +++ b/src/mintpy/cli/smallbaselineApp.py @@ -6,15 +6,24 @@ ############################################################ +# datetime 是 Python 标准库，专门用来处理日期和时间；这里用它打印程序开始运行的时间。 import datetime +# os 是 Python 标准库，专门用来和操作系统交互，比如判断文件是否存在、拼接路径、获取当前目录。 import os +# sys 是 Python 标准库，专门用来访问 Python 解释器相关信息；这里用 sys.argv 读取命令行参数，用 sys.exit 退出程序。 import sys +# mintpy 是本项目的主包；这里用它读取软件版本、默认配置文件所在位置等信息。 import mintpy +# STEP_LIST 是 MintPy 默认模板中定义好的处理步骤顺序，例如 load_data、modify_network、velocity 等。 from mintpy.defaults.template import STEP_LIST +# create_argument_parser 是 MintPy 自己封装的 argparse 工具，用来创建风格统一的命令行参数解析器。 from mintpy.utils.arg_utils import create_argument_parser ########################################################################## +# 命令行帮助文本：说明可通过 --start/--end/--dostep 选择运行的处理步骤。 +# """...""" 是 Python 的多行字符串，适合保存较长的帮助说明。 +# .format(...) 会把 STEP_LIST 分成几段填入上面的 {}，这样 help 页面里能显示所有步骤名称。 STEP_HELP = """Command line options for steps processing with names chosen from the following list: {} @@ -33,6 +42,7 @@ doi:10.1016/j.cageo.2019.104331. """ +# EXAMPLE 保存命令行使用示例，会显示在 -h/--help 帮助信息的末尾。 EXAMPLE = """example: smallbaselineApp.py # run with default template 'smallbaselineApp.cfg' smallbaselineApp.py # run with default and custom templates @@ -49,18 +59,29 @@ def create_parser(subparsers=None): + # 创建 smallbaselineApp.py 的命令行参数解析器。 + # “解析器”可以理解为：把用户在终端输入的文字参数，转换成 Python 里可访问的变量。 synopsis = 'Routine Time Series Analysis for Small Baseline InSAR Stack' + # epilog 是帮助信息最后显示的内容，这里把参考文献和示例命令拼在一起。 epilog = REFERENCE + '\n' + EXAMPLE + # __name__ 是 Python 自动提供的模块名；split('.')[-1] 取最后一段，作为命令名称。 name = __name__.split('.')[-1] + # 调用 MintPy 的 create_argument_parser() 方法，得到一个 parser 对象。 + # 后面多次调用 parser.add_argument()，就是不断告诉它“这个命令支持哪些参数”。 parser = create_argument_parser( name, synopsis=synopsis, description=synopsis, epilog=epilog, subparsers=subparsers) + # 位置参数：用户自定义模板文件；可选参数：工作目录。 + # nargs='?' 表示这个位置参数可以不填；如果用户不填，就使用当前目录里的 smallbaselineApp.cfg。 parser.add_argument('customTemplateFile', nargs='?', help='custom template with option settings.\n' + "ignored if the default smallbaselineApp.cfg is input.") + # --dir/--work-dir 是同一个参数的两个名字；dest='workDir' 表示解析后保存到 inps.workDir。 parser.add_argument('--dir', '--work-dir', dest='workDir', default='./', help='work directory, (default: %(default)s).') + # 辅助选项：生成模板、打印模板、查看版本或仅绘图。 + # action='store_true' 表示用户写了这个选项时，对应变量就是 True；没写时就是 False。 parser.add_argument('-g', dest='generate_template', action='store_true', help='generate default template (if it does not exist) and exit.') parser.add_argument('-H', dest='print_template', action='store_true', @@ -70,47 +91,68 @@ def create_parser(subparsers=None): parser.add_argument('--plot', dest='plot', action='store_true', help='plot results [only] without running smallbaselineApp.') + # 步骤控制选项：指定起止步骤，或只运行某一个步骤。 + # add_argument_group() 会在 help 信息中单独分组显示这几个步骤相关参数。 step = parser.add_argument_group('steps processing (start/end/dostep)', STEP_HELP) + # default=STEP_LIST[0] 表示默认从第一个步骤开始运行。 step.add_argument('--start', dest='startStep', metavar='STEP', default=STEP_LIST[0], help='start processing at the named step (default: %(default)s).') + # default=STEP_LIST[-1] 表示默认运行到最后一个步骤；[-1] 是 Python 中“最后一个元素”的写法。 step.add_argument('--end','--stop', dest='endStep', metavar='STEP', default=STEP_LIST[-1], help='end processing at the named step (default: %(default)s)') + # --dostep 用来只运行一个步骤，后面会把 startStep 和 endStep 都改成这个步骤。 step.add_argument('--dostep', dest='doStep', metavar='STEP', help='run processing at the named step only') + # 返回 parser，供 cmd_line_parse() 使用。 return parser def cmd_line_parse(iargs=None): """Command line parser.""" + # 解析命令行输入，返回后续工作流需要的参数对象。 # parse parser = create_parser() + # parser.parse_args() 会真正解析命令行参数。 + # 如果 iargs 是 None，表示从真实终端命令读取；如果传入列表，常用于测试或其它 Python 代码调用。 inps = parser.parse_args(args=iargs) + # 保存原始参数列表，用于判断用户是否只指定了特殊选项。 # save argv (to check the manually specified arguments) # use iargs for python call # use sys.argv[1:] for command line call + # sys.argv 是完整命令行参数列表，第 0 个通常是脚本名，所以用 [1:] 去掉脚本名。 inps.argv = iargs or sys.argv[1:] # check + # os.path.dirname(mintpy.__file__) 得到 mintpy 包所在目录。 + # os.path.join(...) 用正确的路径分隔符拼接路径，避免手写 '/' 在不同系统上出错。 template_file = os.path.join(os.path.dirname(mintpy.__file__), 'defaults/smallbaselineApp.cfg') + # -H 只打印默认模板内容并退出，不进入处理流程。 # check: -H option (print default template) if inps.print_template: + # with open(...) as f 是安全打开文件的写法；代码块结束后文件会自动关闭。 with open(template_file) as f: lines = f.read() try: + # 如果安装了 rich，则用语法高亮显示 cfg 文件。 # syntax highlight via rich from rich.console import Console from rich.syntax import Syntax + # Console 是 rich 的输出对象；Syntax 负责按 cfg 语法给文本上色。 console = Console() console.print(Syntax(lines, "cfg", background_color='default')) except ImportError: + # 如果没有安装 rich，就退回普通 print，保证功能不依赖额外包。 print(lines) + # sys.exit(0) 表示正常退出程序，不再继续执行后面的处理流程。 sys.exit(0) + # -v 只打印 MintPy 版本信息并退出。 # check: -v option (print software version) if inps.version: + # mintpy.version.version_description 是 MintPy 的版本说明字符串。 print(mintpy.version.version_description) sys.exit(0) @@ -118,26 +160,36 @@ def cmd_line_parse(iargs=None): if (not inps.customTemplateFile and not os.path.isfile(os.path.basename(template_file)) and not inps.generate_template): + # 没有用户模板、当前目录也没有默认模板，并且不是生成模板模式时，无法继续运行。 + # os.path.basename(template_file) 只取文件名 smallbaselineApp.cfg，不带前面的目录。 + # os.path.isfile(...) 用来判断这个文件在当前目录下是否真实存在。 parser.print_usage() print(EXAMPLE) msg = "no template file found! It requires:" msg += "\n 1) input a custom template file, OR" msg += "\n 2) there is a default template 'smallbaselineApp.cfg' in current directory." + # raise SystemExit 会终止程序，并把错误信息显示给用户。 raise SystemExit(f'ERROR: {msg}') # check: custom input template file if inps.customTemplateFile: + # 自定义模板必须存在，并转换为绝对路径，避免后续切换工作目录后路径失效。 # check the existence if not os.path.isfile(inps.customTemplateFile): + # FileNotFoundError 是 Python 内置异常，用来表示文件不存在。 raise FileNotFoundError(inps.customTemplateFile) + # os.path.abspath(...) 把相对路径转换成绝对路径，例如 a.txt 变成 /home/.../a.txt。 inps.customTemplateFile = os.path.abspath(inps.customTemplateFile) + # 如果用户传入的其实是默认模板，则按无自定义模板处理。 # ignore if it is the default template file (smallbaselineApp.cfg) + # os.path.basename(...) 只比较文件名，不比较目录。 if os.path.basename(inps.customTemplateFile) == os.path.basename(template_file): inps.customTemplateFile = None + # 单独使用 --plot 时只绘图，不运行任何处理步骤。 # check: --plot option if inps.argv == ['--plot']: plot_only = True @@ -145,35 +197,46 @@ def cmd_line_parse(iargs=None): else: plot_only = False + # 根据 start/end/dostep 选项计算本次需要执行的步骤列表。 # check: --start/end/dostep options + # read_inps2run_steps() 是本文件下面定义的函数，负责把用户输入转换成真正要运行的步骤列表。 inps.runSteps = read_inps2run_steps(inps, step_list=STEP_LIST, plot_only=plot_only) + # 返回 inps；后面的 main() 会把它交给真正的处理工作流。 return inps def read_inps2run_steps(inps, step_list, plot_only=False): """read/get run_steps from input arguments.""" + # 将命令行指定的步骤范围转换成有序步骤列表。 # check: if start/end/do step input is valid for key in ['startStep', 'endStep', 'doStep']: + # vars(inps) 会把 argparse 解析出来的对象转换成字典，便于用字符串 key 取值。 value = vars(inps)[key] if value and value not in step_list: + # 用户输入的步骤名必须来自默认模板定义的 STEP_LIST。 msg = f'Input step not found: {value}' msg += f'\nAvailable steps: {step_list}' raise ValueError(msg) + # --dostep 优先级最高：忽略 --start/--end，只执行指定步骤。 # check: ignore --start/end input if --dostep is specified if inps.doStep: inps.startStep = inps.doStep inps.endStep = inps.doStep # get list of steps to run + # list.index(x) 返回 x 在列表中的位置编号，用它确定起始步骤和结束步骤的位置。 idx0 = step_list.index(inps.startStep) idx1 = step_list.index(inps.endStep) if idx0 > idx1: + # 起始步骤不能位于结束步骤之后。 msg = f'start step "{inps.startStep}" is CAN NOT be after the end step "{inps.endStep}"' raise ValueError(msg) + # Python 切片 step_list[idx0:idx1+1] 表示从 idx0 到 idx1，包括 idx1。 run_steps = step_list[idx0:idx1+1] + # 生成模板或仅绘图时，不执行处理步骤。 # empty the step list if: # a) -g OR # b) iargs == ['--plot'] @@ -182,30 +245,40 @@ def read_inps2run_steps(inps, step_list, plot_only=False): # print mssage - processing steps if len(run_steps) > 0: + # 多步骤运行时打印 logo，单步骤运行时打印紧凑版本信息。 # for single step - compact version info if len(run_steps) == 1: print(mintpy.version.version_description) else: print(mintpy.version.logo) + # datetime.datetime.now() 返回当前日期时间；os.getcwd() 返回当前工作目录。 print(f'--RUN-at-{datetime.datetime.now()}--') print(f'Current directory: {os.getcwd()}') + # os.path.basename(__file__) 只取当前脚本文件名，用于日志输出。 print(f'Run routine processing with {os.path.basename(__file__)} on steps: {run_steps}') print(f'Remaining steps: {step_list[idx0+1:]}') print('-'*50) + # 返回最终步骤列表，主工作流会逐个执行这些步骤。 return run_steps ########################################################################## def main(iargs=None): + # CLI 入口：解析参数后调用主工作流实现。 # parse inps = cmd_line_parse(iargs) + # 延迟导入主工作流，减少只看帮助/版本时的启动开销。 # import + # from ... import ... 表示从另一个 Python 文件/模块里拿到一个函数。 + # 这里导入的是 src/mintpy/smallbaselineApp.py 里的 run_smallbaselineApp()。 from mintpy.smallbaselineApp import run_smallbaselineApp + # 将命令行参数对象交给 src/mintpy/smallbaselineApp.py 中的主流程执行。 # run + # 这里是 CLI 文件和主工作流文件的连接点：CLI 只负责解析参数，真正处理数据在 run_smallbaselineApp() 中。 run_smallbaselineApp(inps) diff --git a/src/mintpy/dem_error.py b/src/mintpy/dem_error.py index d415f53c3..c01c2d6c5 100644 --- a/src/mintpy/dem_error.py +++ b/src/mintpy/dem_error.py @@ -5,17 +5,25 @@ ############################################################ +# os 用来处理路径、判断文件是否存在、比较文件修改时间。 import os +# time 用来统计 DEM 误差校正耗时。 import time +# h5py 用来直接读取 HDF5 数据集大小，判断输出文件是否完整写完。 import h5py +# numpy 用于数组计算；时序数据、设计矩阵、掩膜都用 numpy 表示。 import numpy as np +# scipy.linalg 提供最小二乘求解器，用来估计 DEM 残差参数。 from scipy import linalg +# cluster 负责分块和并行；geometry/timeseries 是 MintPy 读取几何文件和时序文件的对象。 from mintpy.objects import cluster, geometry, timeseries +# ptime 处理日期；time_func 构建设计矩阵；readfile/writefile 负责 HDF5 读写；ut 是通用工具。 from mintpy.utils import ptime, readfile, time_func, utils as ut, writefile # key configuration parameter name +# 这些 topographicResidual 配置会写入输出文件属性，用于 update 模式判断是否需要重跑。 key_prefix = 'mintpy.topographicResidual.' config_keys = [ 'polyOrder', @@ -30,6 +38,7 @@ ############################################################################ def run_or_skip(inps): + # update 模式判断：输出存在、完整、比输入新、关键配置没变，则跳过。 print('-'*50) print('update mode: ON') flag = 'skip' @@ -42,6 +51,7 @@ def run_or_skip(inps): # check if time-series file is partly written using file size # since time-series file is not compressed with h5py.File(inps.ts_cor_file, 'r') as f: + # timeseries 是 float32，理论字节数约等于元素个数 * 4。 fsize_ref = f['timeseries'].size * 4 fsize = os.path.getsize(inps.ts_cor_file) if fsize <= fsize_ref: @@ -82,12 +92,14 @@ def run_or_skip(inps): ############################################################################ def read_template2inps(template_file, inps): """Read input template file into inps.excludeDate""" + # 从 smallbaselineApp.cfg 中读取 DEM 残差校正相关配置，写入 inps 对象。 iDict = vars(inps) print('read options from template file: '+os.path.basename(template_file)) template = readfile.read_template(template_file, skip_chars=['[', ']']) template = ut.check_template_auto_value(template) # Read template option + # key_prefix+i 对应模板项，例如 mintpy.topographicResidual.polyOrder。 keyList = [i for i in list(iDict.keys()) if key_prefix+i in template.keys()] for key in keyList: value = template[key_prefix+key] @@ -100,6 +112,7 @@ def read_template2inps(template_file, inps): iDict[key] = ptime.yyyymmdd(value.split(',')) # computing configurations + # 同时读取并行/内存相关配置，例如 mintpy.compute.maxMemory。 dask_key_prefix = 'mintpy.compute.' keyList = [i for i in list(iDict.keys()) if dask_key_prefix+i in template.keys()] for key in keyList: @@ -123,12 +136,14 @@ def read_exclude_date(ex_date_list, date_list_all, print_msg=True): Returns: date_flag : 1D array of bool in size of (num_date,) """ # Read exclude date input + # ptime.read_date_list() 可以读取日期列表，也可以读取包含日期的文本文件。 ex_date_list = ptime.read_date_list(ex_date_list) if ex_date_list and print_msg: print(('exclude the following dates for DEM error estimation:' ' ({})\n{}').format(len(ex_date_list), ex_date_list)) # convert to mark array + # date_flag=True 表示该日期参与 DEM 误差估计；False 表示排除。 date_flag = np.array([i not in ex_date_list for i in date_list_all], dtype=np.bool_) return date_flag, ex_date_list @@ -141,6 +156,7 @@ def get_design_matrix4defo(inps): """ # key msg + # 这里构建的是“真实地表形变”的时间函数模型，用来和 DEM 误差项一起联合拟合。 msg = '-'*80 msg += '\ncorrect topographic phase residual (DEM error) (Fattahi & Amelung, 2013, IEEE-TGRS)' msg += '\nordinal least squares (OLS) inversion with L2-norm minimization on: phase' @@ -155,17 +171,20 @@ def get_design_matrix4defo(inps): print(msg) # prepare temporal deformation model + # model 字典描述时间函数：多项式、阶跃项、周期项等。 model = dict() model['polynomial'] = inps.polyOrder model['stepDate'] = inps.stepDate model['periodic'] = inps.periodic # prepare SAR info + # 时间函数设计矩阵需要知道每期 SAR 影像日期和成像时间。 ts_obj = timeseries(inps.ts_file) date_list = ts_obj.get_date_list() seconds = ts_obj.get_metadata().get('CENTER_LINE_UTC', 0) # compose design matrix + # G_defo 每一列是一个时间函数参数，例如常数项、速度项、阶跃项等。 G_defo = time_func.get_design_matrix4time_func(date_list, model, seconds=seconds) return G_defo @@ -180,16 +199,19 @@ def read_geometry(ts_file, geom_file=None, box=None): range_dist - 0/2D array, slant range distance in meter pbase - 0/3D array, perp baseline in meter """ + # DEM 误差相位与入射角、斜距、垂直基线有关；这个函数读取这些几何量。 ts_obj = timeseries(ts_file) ts_obj.open(print_msg=False) # 0/2/3D geometry if geom_file: + # 有 geometry 文件时，优先读取二维/三维几何数据，精度更高。 geom_obj = geometry(geom_file) geom_obj.open() # 0/2D incidence angle / slant range distance if 'incidenceAngle' not in geom_obj.datasetNames: + # 如果几何文件里没有 incidenceAngle，就从时序元数据估算平均值。 inc_angle = ut.incidence_angle(ts_obj.metadata, dimension=0) range_dist = ut.range_distance(ts_obj.metadata, dimension=0) else: @@ -200,6 +222,7 @@ def read_geometry(ts_file, geom_file=None, box=None): # 0/3D perp baseline if 'bperp' in geom_obj.datasetNames: + # bperp 是每期影像的垂直基线；三维 bperp 表示每期、每个像元都可能不同。 print(f'read 3D bperp from {geom_obj.name} file: {os.path.basename(geom_obj.file)} ...') dset_list = [f'bperp-{d}' for d in ts_obj.dateList] pbase = geom_obj.read(datasetName=dset_list, box=box, print_msg=False).reshape((ts_obj.numDate, -1)) @@ -210,6 +233,7 @@ def read_geometry(ts_file, geom_file=None, box=None): # 0D geometry else: + # 没有 geometry 文件时，只能使用元数据中的平均几何值，空间变化会被忽略。 print(f'read mean incidenceAngle, slantRangeDistance, bperp value from {ts_obj.name} file') inc_angle = ut.incidence_angle(ts_obj.metadata, dimension=0) range_dist = ut.range_distance(ts_obj.metadata, dimension=0) @@ -235,17 +259,21 @@ def estimate_dem_error(ts0, G0, tbase, date_flag=None, phase_velocity=False, residual timeseries = tsOrig - delta_z_phase - defModel Example: delta_z, ts_cor, ts_res = estimate_dem_error(ts, G, tbase, date_flag) """ + # 这个函数是 DEM 误差估计的核心：解线性方程 G0 * X = ts0。 + # X 的第一个参数 delta_z 就是每个像元估计出的 DEM 高程残差。 if len(ts0.shape) == 1: ts0 = ts0.reshape(-1, 1) if date_flag is None: date_flag = np.ones(ts0.shape[0], np.bool_) # skip noisy acquisitions + # date_flag 会排除用户指定不参与估计的日期。 G = G0[date_flag, :] ts = ts0[date_flag, :] if phase_velocity: # adjust from phase to phase velocity + # phase_velocity=True 时，不拟合位移本身，而拟合相邻日期之间的速度变化。 tbase_diff = np.diff(tbase[date_flag], axis=0).reshape(-1,1) ts = np.diff(ts, axis=0) / np.repeat(tbase_diff, ts.shape[1], axis=1) G = np.diff(G, axis=0) / np.repeat(tbase_diff, G.shape[1], axis=1) @@ -260,9 +288,11 @@ def estimate_dem_error(ts0, G0, tbase, date_flag=None, phase_velocity=False, # X = [delta_z, constC, vel, acc, deltaAcc, ..., step1, step2, ...] # equivalent to X = np.dot(np.dot(np.linalg.inv(np.dot(G.T, G)), G.T), ts) # X = np.dot(np.linalg.pinv(G), ts) + # linalg.lstsq() 进行普通最小二乘估计，返回最能解释观测时序的参数。 X = linalg.lstsq(G, ts, cond=cond)[0] # prepare outputs + # ts_cor 是扣除 DEM 误差后的时序；ts_res 是扣除 DEM 误差和形变模型后的残差。 delta_z = X[0, :] ts_cor = ts0 - np.dot(G0[:, 0].reshape(-1, 1), delta_z.reshape(1, -1)) ts_res = ts0 - np.dot(G0, X) @@ -297,6 +327,7 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None, box - tuple of 4 int in (x0, y0, x1, y1) for the area of interest """ + # patch 表示图像的一块区域；大图按块处理以降低内存占用。 ts_obj = timeseries(ts_file) ts_obj.open(print_msg=False) @@ -320,12 +351,15 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None, num_date = ts_obj.numDate # 1.1 read time-series + # 读取当前 patch 的位移时序，并展开成 num_date x num_pixel 的二维矩阵。 ts_data = readfile.read(ts_file, box=box)[0].reshape(num_date, -1) # 1.2 read geometry + # 读取当前 patch 的入射角、斜距、垂直基线，用于构建 DEM 误差项。 sin_inc_angle, range_dist, pbase = read_geometry(ts_file, geom_file, box=box) # 1.3 mask of pixels to invert + # 下面几步构建有效像元掩膜，跳过全零、NaN、低相干或几何异常的像元。 print('skip pixels with ZERO in ALL acquisitions') mask = np.nanmean(ts_data, axis=0) != 0. @@ -367,6 +401,7 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None, # 2.2 estimate if range_dist.size == 1: + # 几何量是单个平均值时，所有像元共享同一个设计矩阵，可以一次性反演所有像元。 print('estimating DEM error ...') # compose design matrix G_geom = pbase / (range_dist * sin_inc_angle) @@ -387,6 +422,7 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None, ts_res[:, mask] = ts_res_i else: + # 几何量逐像元变化时，每个像元的 G_geom 不同，需要逐像元估计。 print('estimating DEM error pixel-wisely ...') prog_bar = ptime.progressBar(maxValue=num_pixel2inv) for i in range(num_pixel2inv): @@ -430,6 +466,7 @@ def correct_dem_error_patch(G_defo, ts_file, geom_file=None, box=None, def correct_dem_error(inps): """Correct DEM error of input timeseries file""" + # correct_dem_error() 是本模块主入口：准备输出文件、切块、调用 patch 函数、写入结果。 start_time = time.time() # limit the number of threads to 1 @@ -453,6 +490,7 @@ def correct_dem_error(inps): inps.polyOrder, np.sum(date_flag))) # 1.2 design matrix part 1 - time func for surface deformation + # G_defo 只包含地表形变时间函数；DEM 几何项会在每个 patch/像元中再拼进去。 G_defo = get_design_matrix4defo(inps) @@ -465,16 +503,19 @@ def correct_dem_error(inps): meta[key_prefix+key] = str(vars(inps)[key]) # 2.2 instantiate est. DEM error + # dem_err_file 保存估计出的 DEM 高程残差，单位米。 meta['FILE_TYPE'] = 'dem' meta['UNIT'] = 'm' ds_name_dict = {'dem' : [np.float32, (length, width), None]} writefile.layout_hdf5(inps.dem_err_file, ds_name_dict, metadata=meta) # 2.3 instantiate corrected time-series + # ts_cor_file 保存扣除 DEM 残差相位后的校正时序。 meta['FILE_TYPE'] = 'timeseries' writefile.layout_hdf5(inps.ts_cor_file, metadata=meta, ref_file=inps.ts_file) # 2.4 instantiate residual phase time-series + # timeseriesResidual.h5 保存拟合残差，后续 residual_RMS 会用它评估噪声。 ts_res_file = os.path.join(os.path.dirname(inps.ts_cor_file), 'timeseriesResidual.h5') writefile.layout_hdf5(ts_res_file, metadata=meta, ref_file=inps.ts_file) @@ -482,6 +523,7 @@ def correct_dem_error(inps): ## 3. run the estimation and write to disk # 3.1 split ts_file into blocks to save memory + # 根据日期数、图像大小和内存上限估算需要切成多少块。 # 1st dimension size: ts (obs / cor / res / step) + dem_err/inc_angle/rg_dist (+pbase) num_epoch = num_date * 3 + num_step + 3 if inps.geom_file: @@ -499,6 +541,7 @@ def correct_dem_error(inps): ) # 3.2 prepare the input arguments for *_patch() + # data_kwargs 是传给 correct_dem_error_patch() 的固定参数；循环里只更新 box。 data_kwargs = { 'G_defo' : G_defo, 'ts_file' : inps.ts_file, @@ -522,10 +565,12 @@ def correct_dem_error(inps): # invert if not inps.cluster: # non-parallel + # 单进程模式：直接处理当前 patch。 delta_z, ts_cor, ts_res = correct_dem_error_patch(**data_kwargs)[:-1] else: # parallel + # 并行模式：用 Dask 对当前 patch 内计算进一步并行。 print('\n\n------- start parallel processing using Dask -------') # initiate the output data @@ -552,6 +597,7 @@ def correct_dem_error(inps): # write the block to disk # with 3D block in [z0, z1, y0, y1, x0, x1] # and 2D block in [y0, y1, x0, x1] + # 分块写入 HDF5，保证每个 patch 只写自己负责的空间范围。 # DEM error - 2D block = [box[1], box[3], box[0], box[2]] diff --git a/src/mintpy/geocode.py b/src/mintpy/geocode.py index ed6d5cf76..aa196a1be 100644 --- a/src/mintpy/geocode.py +++ b/src/mintpy/geocode.py @@ -5,29 +5,41 @@ ############################################################ +# os 用来处理路径、创建输出目录、判断文件大小等。 import os +# time 用来统计地理编码处理耗时。 import time +# numpy 用来处理数组，例如创建输出数组、判断最大文件索引。 import numpy as np +# resample 是 MintPy 的重采样对象，负责根据查找表把数据从一个坐标系重采样到另一个坐标系。 from mintpy.objects.resample import resample +# attribute(attr) 负责更新文件元数据；readfile/writefile 负责读写数据；ut 是通用工具函数。 from mintpy.utils import attribute as attr, readfile, utils as ut, writefile ############################################################################################ def auto_output_filename(in_file, inps): + # 自动生成输出文件名。 + # 如果只处理一个输入文件，并且用户明确指定了输出文件名，就直接使用用户指定的名字。 if len(inps.file) == 1 and inps.outfile: return inps.outfile + # os.path.basename() 取文件名；os.path.splitext() 把文件名拆成“主名”和“扩展名”。 fbase, fext = os.path.splitext(os.path.basename(in_file)) + # radar2geo=True 表示雷达坐标转地理坐标，输出名前缀用 geo_；否则用 rdr_。 prefix = 'geo_' if inps.radar2geo else 'rdr_' + # 如果用户只处理某个数据集 dset，就用数据集名做输出名后缀；否则用原文件主名。 suffix = inps.dset if inps.dset else fbase out_file = f'{prefix}{suffix}{fext}' if inps.out_dir: + # 如果指定输出目录但目录不存在，就先创建目录。 if not os.path.isdir(inps.out_dir): os.makedirs(inps.out_dir) print(f'create directory: {inps.out_dir}') + # os.path.join() 把目录和文件名拼成完整路径。 out_file = os.path.join(inps.out_dir, out_file) return out_file @@ -35,12 +47,17 @@ def auto_output_filename(in_file, inps): def run_geocode(inps): """geocode all input files""" + # geocode 的核心任务：把输入文件中的二维/三维栅格数据重采样到另一个坐标系。 + # 常见方向是 radar2geo：雷达坐标 -> 经纬度地理坐标。 start_time = time.time() # feed the largest file for resample object initiation + # 选择最大的输入文件初始化 resample 对象，通常可以代表最完整的数据尺寸和结构。 + # os.path.getsize(i) 返回文件大小；np.argmax() 返回最大值的索引。 ind_max = np.argmax([os.path.getsize(i) for i in inps.file]) # prepare geometry for geocoding + # kwargs 是传给 resample 对象的配置字典，包括插值方法、填充值、并行数和内存限制。 kwargs = dict(interp_method=inps.interpMethod, fill_value=inps.fillValue, nprocs=inps.nprocs, @@ -48,42 +65,52 @@ def run_geocode(inps): software=inps.software, print_msg=True) if inps.latFile and inps.lonFile: + # 有些数据用经纬度文件而不是 lookupFile 表示坐标转换关系。 kwargs['lat_file'] = inps.latFile kwargs['lon_file'] = inps.lonFile + # 创建 resample 对象。lut_file 是 lookup table，描述源坐标和目标坐标的对应关系。 res_obj = resample(lut_file=inps.lookupFile, src_file=inps.file[ind_max], SNWE=inps.SNWE, lalo_step=inps.laloStep, **kwargs) + # open() 读取查找表和源文件的基本信息；prepare() 预计算重采样需要的块和索引。 res_obj.open() res_obj.prepare() # resample input files one by one for infile in inps.file: print('-' * 50+f'\nresampling file: {infile}') + # 读取输入文件属性；datasetName=inps.dset 表示可只读取某个指定数据集的属性。 atr = readfile.read_attribute(infile, datasetName=inps.dset) outfile = auto_output_filename(infile, inps) # update_mode if inps.updateMode: print('update mode: ON') + # 如果输出文件已存在且比输入文件/查找表新，就跳过，避免重复计算。 if ut.run_or_skip(outfile, in_file=[infile, inps.lookupFile]) == 'skip': continue ## prepare output # update metadata if inps.radar2geo: + # 雷达坐标转地理坐标时，需要把元数据中的坐标信息更新成经纬度网格。 atr = attr.update_attribute4radar2geo(atr, res_obj=res_obj) else: + # 地理坐标转雷达坐标时，需要把元数据更新成雷达坐标尺寸和坐标定义。 atr = attr.update_attribute4geo2radar(atr, res_obj=res_obj) # instantiate output file + # 判断输出是否是 HDF5 文件；.h5 和 .he5 需要先创建文件结构，再分块写数据。 hdf5_file = os.path.splitext(outfile)[1] in ['.h5', '.he5'] if hdf5_file: # grab metadata from input file: dataset compression and UNIT + # 保留原 HDF5 文件的数据压缩方式和每个数据集的 UNIT 属性。 compression = readfile.get_hdf5_compression(infile) ds_unit_dict = readfile.get_hdf5_dataset_attrs(infile, key='UNIT') # initiate output file + # layout_hdf5() 先创建空的 HDF5 文件和数据集布局，后面再按块写入数据。 writefile.layout_hdf5( outfile, metadata=atr, @@ -92,18 +119,22 @@ def run_geocode(inps): compression=compression, ) else: + # 非 HDF5 输出先用字典在内存里收集所有数据集，最后统一写出。 dsDict = dict() ## run + # get_dataset_list() 获取要处理的数据集列表；如果指定 inps.dset，就只返回相关数据集。 dsNames = readfile.get_dataset_list(infile, datasetName=inps.dset) maxDigit = max(len(i) for i in dsNames) for dsName in dsNames: if not hdf5_file: + # 非 HDF5 文件需要先创建目标大小的空数组。 dsDict[dsName] = np.zeros((res_obj.length, res_obj.width), dtype=atr['DATA_TYPE']) # loop for block-by-block IO for i in range(res_obj.num_box): + # 为了节省内存，大文件不会一次性读入，而是按 res_obj 预先划分的小块逐块处理。 src_box = res_obj.src_box_list[i] dest_box = res_obj.dest_box_list[i] @@ -118,42 +149,51 @@ def run_geocode(inps): print_msg=False)[0] # resample + # run_resample() 对当前数据块执行坐标转换和插值。 data = res_obj.run_resample(src_data=data, box_ind=i) # write / save block data if data.ndim == 3: + # 三维数据通常是 time/date x row x col，因此 block 要包含第 0 维范围。 block = [0, data.shape[0], dest_box[1], dest_box[3], dest_box[0], dest_box[2]] else: + # 二维数据只需要行列范围。 block = [dest_box[1], dest_box[3], dest_box[0], dest_box[2]] if hdf5_file: print(f'write data in block {block} to file: {outfile}') + # HDF5 文件直接把当前块写入目标文件对应位置。 writefile.write_hdf5_block(outfile, data=data, datasetName=dsName, block=block, print_msg=False) else: + # 非 HDF5 文件先把当前块放进内存数组对应位置。 dsDict[dsName][block[0]:block[1], block[2]:block[3]] = data # for binary file: ensure same data type if not hdf5_file: + # 保持输出数组数据类型和实际重采样结果一致。 dsDict[dsName] = np.array(dsDict[dsName], dtype=data.dtype) # write binary file if not hdf5_file: + # 非 HDF5 格式最后一次性写出，并附带更新后的元数据。 atr['BANDS'] = len(dsDict.keys()) writefile.write(dsDict, out_file=outfile, metadata=atr, ref_file=infile) # create ISCE XML and GDAL VRT file if using ISCE lookup table file if inps.latFile and inps.lonFile: + # ISCE 生态常用 XML/VRT 描述二进制文件，这里一并生成。 writefile.write_isce_xml(atr, fname=outfile) # used time + # 统计所有输入文件地理编码总耗时。 m, s = divmod(time.time()-start_time, 60) print(f'time used: {m:02.0f} mins {s:02.1f} secs.\n') diff --git a/src/mintpy/ifgram_inversion.py b/src/mintpy/ifgram_inversion.py index 16c8e1dc2..ee79210a3 100644 --- a/src/mintpy/ifgram_inversion.py +++ b/src/mintpy/ifgram_inversion.py @@ -8,18 +8,27 @@ # from mintpy import ifgram_inversion as ifginv +# os 用来处理文件路径、判断文件是否存在、比较文件修改时间。 import os +# time 用来统计反演运行耗时。 import time +# h5py 用来直接读写 HDF5 文件中的数据集和属性。 import h5py +# numpy 用来做数组计算；反演过程中相位栈、时间序列、掩膜都以 numpy 数组表示。 import numpy as np +# scipy.linalg 提供线性代数求解器；这里比 numpy.linalg 更适合大型最小二乘问题。 from scipy import linalg # more effieint than numpy.linalg +# cluster 负责分块/并行；ifgramStack 是读取 ifgramStack.h5 的对象。 from mintpy.objects import cluster, ifgramStack +# decorrelation 模块把相干性转换成权重，用于加权最小二乘反演。 from mintpy.simulation import decorrelation as decor +# ptime 处理日期；readfile/writefile 读写文件；ut 是通用工具函数。 from mintpy.utils import ptime, readfile, utils as ut, writefile # key configuration parameter name +# 这些配置项会写入输出文件属性；下次 update 模式用它们判断是否需要重跑。 key_prefix = 'mintpy.networkInversion.' config_keys = [ 'obsDatasetName', @@ -36,6 +45,7 @@ def run_or_skip(inps): + # update 模式判断：如果输出文件已存在、比输入新、关键配置没变，就跳过反演。 print('-'*50) print('update mode: ON') flag = 'skip' @@ -48,6 +58,7 @@ def run_or_skip(inps): # check if time-series file is partly written using file size # since time-series file is not compressed with h5py.File(inps.outfile[0], 'r') as f: + # timeseries 数据集是 float32，理论数据字节数约为元素个数 * 4。 fsize_ref = f['timeseries'].size * 4 fsize = os.path.getsize(inps.outfile[0]) if fsize <= fsize_ref: @@ -58,6 +69,7 @@ def run_or_skip(inps): print(f'1) output files already exist: {inps.outfile}.') # check modification time with h5py.File(inps.ifgramStackFile, 'r') as f: + # MODIFICATION_TIME 是输入数据集的修改时间；没有这个属性时退回文件修改时间。 ti = float(f[inps.obsDatasetName].attrs.get('MODIFICATION_TIME', os.path.getmtime(inps.ifgramStackFile))) to = min(os.path.getmtime(i) for i in inps.outfile) if ti > to: @@ -70,6 +82,7 @@ def run_or_skip(inps): if flag == 'skip': atr_ifg = readfile.read_attribute(inps.ifgramStackFile) atr_ts = readfile.read_attribute(inps.tsFile) + # 当前实际参与反演的干涉图数量，受 dropIfgram 网络筛选影响。 inps.numIfgram = len(ifgramStack(inps.ifgramStackFile).get_date12_list(dropIfgram=True)) meta_keys = [i for i in ['REF_Y', 'REF_X'] if i in atr_ts.keys()] @@ -141,8 +154,10 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity used during the inversion """ + # y 是观测值矩阵：行是干涉图，列是像元。reshape 保证它是二维数组。 y = y.reshape(A.shape[0], -1) if weight_sqrt is not None: + # weight_sqrt 是权重的平方根，用于加权最小二乘：sqrt(W) * G * x = sqrt(W) * y。 weight_sqrt = weight_sqrt.reshape(A.shape[0], -1) num_date = A.shape[1] + 1 num_pixel = y.shape[1] @@ -156,6 +171,7 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity num_inv_obs = 0 ##### skip invalid phase/offset value [NaN] + # 对单个像元，如果某些干涉图是 NaN，就同时从观测值和设计矩阵中删掉对应行。 y, [A, B, weight_sqrt] = skip_invalid_obs(y, mat_list=[A, B, weight_sqrt]) # check 1 - network redundancy: skip inversion if < threshold @@ -177,6 +193,7 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity try: if min_norm_velocity: ##### min-norm velocity + # 最小范数速度模式使用 B 矩阵求解速度增量，再累加成位移时序。 if weight_sqrt is not None: X, e2 = linalg.lstsq(np.multiply(B, weight_sqrt), np.multiply(y, weight_sqrt), @@ -192,11 +209,13 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity print_msg=print_msg) # assemble time-series + # X 是速度；乘以相邻日期时间间隔 tbase_diff 得到位移增量，再 cumsum 累加成时序。 ts_diff = X * np.tile(tbase_diff, (1, num_pixel)) ts[1:, :] = np.cumsum(ts_diff, axis=0) else: ##### min-norm displacement + # 最小范数位移模式直接使用 A 矩阵求解每个日期相对参考日期的位移。 if weight_sqrt is not None: X, e2 = linalg.lstsq(np.multiply(A, weight_sqrt), np.multiply(y, weight_sqrt), @@ -215,6 +234,7 @@ def estimate_timeseries(A, B, y, tbase_diff, weight_sqrt=None, min_norm_velocity ts[1: ,:] = X except linalg.LinAlgError: + # 线性方程不可解或矩阵数值异常时，保留默认输出，避免整个批块崩溃。 pass # number of observations used for inversion @@ -238,6 +258,7 @@ def estimate_timeseries_cov(G, y, y_std, rcond=1e-5, min_redundancy=1.0): y_std - 2D np.ndarray in size of (num_pair, 1), stack of obs std. dev. Returns: ts_cov - 2D np.ndarray in size of (num_date-1, num_date-1), time-series obs std. dev. """ + # 这个函数用于不确定性传播：把干涉图观测标准差传播成时间序列协方差。 y = y.reshape(G.shape[0], -1) y_std = y_std.reshape(G.shape[0], -1) @@ -259,7 +280,9 @@ def estimate_timeseries_cov(G, y, y_std, rcond=1e-5, min_redundancy=1.0): ## TS var. --> TS std. dev. #ts_cov = np.sqrt(ts_var) Gplus = linalg.pinv(G) + # stack_cov 是观测协方差矩阵；这里假设不同干涉图之间相互独立，所以是对角矩阵。 stack_cov = np.diag(np.square(y_std.flatten())) + # 线性传播公式：Cov(x) = G+ * Cov(y) * (G+)^T。 ts_cov = np.linalg.multi_dot([Gplus, stack_cov, Gplus.T]) return ts_cov @@ -275,6 +298,7 @@ def skip_invalid_obs(obs, mat_list): """ if np.any(np.isnan(obs)): # get flag matrix + # flag=True 表示该干涉图观测值有效；False 表示当前像元在这幅干涉图中无效。 flag = (~np.isnan(obs[:, 0])).flatten() # update obs @@ -303,12 +327,14 @@ def calc_inv_quality(G, X, y, e2, inv_quality_name='temporalCoherence', weight_s Returns: inv_quality - 1D np.ndarray in size of (num_pixel), temporalCoherence / residual """ + # 反演质量用于判断时序结果可信度：相位用 temporalCoherence，偏移用 residual。 num_pair, num_pixel = y.shape inv_quality = np.zeros(num_pixel, dtype=np.float32) # chunk_size as the number of pixels chunk_size = int(ut.round_to_1(2e5 / num_pair)) if num_pixel > chunk_size: + # 像元很多时分块计算质量指标，避免一次性创建巨大残差矩阵占满内存。 num_chunk = int(np.ceil(num_pixel / chunk_size)) if print_msg: print('calculating {} in chunks of {} pixels: {} chunks in total ...'.format( @@ -321,6 +347,7 @@ def calc_inv_quality(G, X, y, e2, inv_quality_name='temporalCoherence', weight_s if inv_quality_name == 'temporalCoherence': #for phase + # e 是残差相位；exp(1j*e) 把残差转为单位复数，平均长度越接近 1 表示越一致。 e = y[:, c0:c1] - np.dot(G, X[:, c0:c1]) inv_quality[c0:c1] = np.abs(np.sum(np.exp(1j*e), axis=0)) / num_pair @@ -370,6 +397,7 @@ def check_design_matrix(ifgram_file, weight_func='var'): Check Rank of Design matrix for weighted inversion """ + # 设计矩阵必须满秩，才能把干涉图网络稳定反演成时间序列。 date12_list = ifgramStack(ifgram_file).get_date12_list(dropIfgram=True) A = ifgramStack.get_design_matrix4timeseries(date12_list)[0] if weight_func == 'no': @@ -399,6 +427,7 @@ def read_stack_obs(stack_obj, box, ref_phase, obs_ds_name='unwrapPhase', dropIfg Returns: stack_obs - 2D array of unwrapPhase in size of (num_pair, num_pixel) """ # Read unwrapPhase + # dropIfgram=True 表示只读取网络筛选后保留的干涉图。 num_pair = stack_obj.get_size(dropIfgram=dropIfgram)[0] if print_msg: print(f'reading {obs_ds_name} in {box} * {num_pair} ...') @@ -426,6 +455,7 @@ def read_stack_obs(stack_obj, box, ref_phase, obs_ds_name='unwrapPhase', dropIfg # reference unwrapPhase for i in range(num_pair): + # 每幅干涉图减去参考像元相位，使所有相位相对于同一个参考点。 mask = stack_obs[i, :] != 0. stack_obs[i, :][mask] -= ref_phase[i] return stack_obs @@ -436,6 +466,7 @@ def mask_stack_obs(stack_obs, stack_obj, box, mask_ds_name=None, mask_threshold= """Mask input unwrapped phase by setting them to np.nan.""" # Read/Generate Mask + # 掩膜的目的：把低质量像元设为 NaN，让最小二乘反演自动忽略它们。 num_pair = stack_obj.get_size(dropIfgram=dropIfgram)[0] if mask_ds_name and mask_ds_name in stack_obj.datasetNames: if print_msg: @@ -451,11 +482,13 @@ def mask_stack_obs(stack_obs, stack_obj, box, mask_ds_name=None, mask_threshold= msk = np.ones(msk_data.shape, dtype=np.bool_) if mask_ds_name in ['connectComponent']: + # connectComponent 为 0 通常表示该像元解缠不可靠。 msk *= msk_data != 0 if print_msg: print(f'mask out pixels with {mask_ds_name} == 0 by setting them to NaN') elif mask_ds_name in ['coherence', 'offsetSNR']: + # coherence/offsetSNR 小于阈值表示质量低，设为无效。 msk *= msk_data >= mask_threshold if print_msg: print(f'mask out pixels with {mask_ds_name} < {mask_threshold} by setting them to NaN') @@ -514,6 +547,7 @@ def calc_weight_sqrt(stack_obj, box, weight_func='var', dropIfgram=True, chunk_s print('calculating weight from spatial coherence ...') # read coherence + # 相干性越高，干涉图通常越可靠；加权反演会给高相干观测更大权重。 weight = read_coherence(stack_obj, box=box, dropIfgram=dropIfgram) num_pixel = weight.shape[1] @@ -532,6 +566,7 @@ def calc_weight_sqrt(stack_obj, box, weight_func='var', dropIfgram=True, chunk_s ': {n} chunks in total ...').format(c=chunk_size, n=num_chunk)) for i in range(num_chunk): + # 分块把 coherence 转成 weight，减少内存峰值。 c0 = i * chunk_size c1 = min((i + 1) * chunk_size, num_pixel) if i == 0: @@ -564,6 +599,7 @@ def get_design_matrix4std(stack_obj): """ # get ref_date from template file + # 时间序列标准差计算需要明确参考日期，因为协方差矩阵要相对于该日期组织。 mintpy_dir = os.path.dirname(os.path.dirname(stack_obj.file)) cfg_file = os.path.join(mintpy_dir, 'smallbaselineApp.cfg') ref_date = readfile.read_template(cfg_file)['mintpy.reference.date'] @@ -618,6 +654,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam Example: run_ifgram_inversion_patch('ifgramStack.h5', box=(0,200,1316,400)) """ + # patch 表示图像的一块区域。大图会被切成多个 patch 分别反演，降低内存占用。 stack_obj = ifgramStack(ifgram_file) stack_obj.open(print_msg=False) stack_dir, stack_base = os.path.split(ifgram_file) @@ -630,6 +667,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam ## 1. input info # size + # box 是 (x0, y0, x1, y1)。如果 box=None，就处理整幅图。 if box: num_row = box[3] - box[1] num_col = box[2] - box[0] @@ -639,12 +677,14 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam num_pixel = num_row * num_col # get tbase_diff in the unit of year + # tbase_diff 是相邻日期之间的时间间隔，单位为年，用于把速度积分成位移。 date_list = stack_obj.get_date_list(dropIfgram=True) num_date = len(date_list) tbase = np.array(ptime.date_list2tbase(date_list)[0], np.float32) / 365.25 tbase_diff = np.diff(tbase).reshape(-1, 1) # design matrix + # A/B 是由干涉图网络生成的设计矩阵，是“干涉图相位差 -> 日期时序”的数学桥梁。 date12_list = stack_obj.get_date12_list(dropIfgram=True) A, B = stack_obj.get_design_matrix4timeseries(date12_list=date12_list)[0:2] @@ -655,6 +695,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam if obs_ds_name.startswith(('unwrapPhase', 'ion')): # calculate weight if weight_func not in ['no', 'sbas']: + # 相位观测可以根据 coherence 计算权重。 weight_sqrt = calc_weight_sqrt(stack_obj, box, weight_func=weight_func, dropIfgram=True, @@ -675,6 +716,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam elif 'offset' in obs_ds_name.lower(): if calc_cov or weight_func == 'var': # calculate weight for offset + # 偏移观测的权重来自对应的 Std 数据集，标准差越小权重越大。 print('reading {} in {} * {} ...'.format(obs_ds_name+'Std', box, len(date12_list))) weight_sqrt = stack_obj.read(datasetName=obs_ds_name+'Std', box=box, @@ -704,6 +746,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam raise ValueError(f'un-recognized observation dataset name: {obs_ds_name}') # 1.2 read / mask unwrapPhase and offset + # 读取观测栈数据，并减去参考相位，得到可用于反演的观测矩阵。 stack_obs = read_stack_obs(stack_obj, box, ref_phase, obs_ds_name=obs_ds_name, dropIfgram=True) @@ -722,10 +765,12 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam dropIfgram=True) # 1.3 mask of pixels to invert + # mask 是一维布尔数组，长度等于当前 patch 的像元数；True 表示该像元需要反演。 mask = np.ones(num_pixel, np.bool_) # 1.3.1 - Water Mask if water_mask_file: + # 水体区域通常没有稳定 InSAR 相位，因此跳过 waterMask 为 0 的像元。 print(f'skip pixels (on the water) with zero value in file: {os.path.basename(water_mask_file)}') atr_msk = readfile.read_attribute(water_mask_file) len_msk, wid_msk = int(atr_msk['LENGTH']), int(atr_msk['WIDTH']) @@ -739,6 +784,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam del waterMask # 1.3.2 - Mask for NaN value in ALL ifgrams + # 如果一个像元在所有干涉图里都是 NaN，没有任何观测，无法反演。 print(f'skip pixels with {obs_ds_name} = NaN in all interferograms') mask *= ~np.all(np.isnan(stack_obs), axis=0) @@ -757,6 +803,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam stack_quality_file = os.path.join(stack_dir, '../avgSpatialCoh.h5') if os.path.isfile(stack_quality_file): + # 如果有平均相干性/SNR 文件，则把质量指标为 0 的像元也跳过。 atr_stack = readfile.read_attribute(stack_quality_file) len_stack, wid_stack = int(atr_stack['LENGTH']), int(atr_stack['WIDTH']) if (len_stack, wid_stack) == (stack_obj.length, stack_obj.width): @@ -775,6 +822,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam ## 2. inversion # 2.1 initiale the output matrices + # 先创建当前 patch 的输出数组，后面只把有效像元位置填进去。 ts = np.zeros((num_date, num_pixel), np.float32) ts_cov = np.zeros((num_date, num_date, num_pixel), np.float32) if calc_cov else None inv_quality = np.zeros(num_pixel, np.float32) @@ -802,6 +850,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam # 2.2 un-weighted inversion (classic SBAS) if weight_sqrt is None: + # 非加权反演是经典 SBAS。对“所有干涉图都有有效观测”的像元，可以一次性矩阵求解。 msg = f'estimating time-series for pixels with valid {obs_ds_name} in' # a. split mask into mask_all/part_net @@ -829,6 +878,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam # c. pixel-by-pixel for pixels with obs not in all ifgrams if np.sum(mask_part_net) > 0: + # 有些像元只在部分干涉图有效，它们的有效观测行不同，需要逐像元反演。 num_pixel2inv_part = int(np.sum(mask_part_net)) idx_pixel2inv_part = np.where(mask_part_net)[0] print(f'{msg} some ifgrams ({num_pixel2inv_part} pixels; {num_pixel2inv_part/num_pixel2inv*100:.1f}%) ...') @@ -853,6 +903,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam # 2.3 weighted inversion - pixel-by-pixel else: + # 加权反演中每个像元权重可能不同，因此逐像元求解。 print('estimating time-series via WLS pixel-by-pixel ...') prog_bar = ptime.progressBar(maxValue=num_pixel2inv) for i in range(num_pixel2inv): @@ -875,6 +926,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam # 2.4 time-series std. dev. - pixel-by-pixel if calc_cov: + # 如果要求协方差，就把干涉图观测标准差传播到时间序列协方差。 print('propagating std. dev. from network of interferograms to time-series (Yunjun et al., 2021, FRINGE) ...') prog_bar = ptime.progressBar(maxValue=num_pixel2inv) for i in range(num_pixel2inv): @@ -907,12 +959,14 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam # 3.2 convert displacement unit to meter if obs_ds_name.startswith(('unwrapPhase','ion')): + # 相位单位是弧度，按雷达波长转换为 LOS 位移，单位米。 phase2range = -1 * float(stack_obj.metadata['WAVELENGTH']) / (4.*np.pi) ts *= phase2range ts_cov = ts_cov * np.abs(phase2range)**2 if calc_cov else ts_cov print('converting LOS phase unit from radian to meter') elif (obs_ds_name == 'azimuthOffset') & (stack_obj.metadata['PROCESSOR'] != 'cosicorr'): + # 方位向偏移原本单位是像素，需要乘以地面分辨率转换为米。 az_pixel_size = ut.azimuth_ground_resolution(stack_obj.metadata) az_pixel_size /= float(stack_obj.metadata['ALOOKS']) ts *= az_pixel_size @@ -920,6 +974,7 @@ def run_ifgram_inversion_patch(ifgram_file, box=None, ref_phase=None, obs_ds_nam print(f'converting azimuth offset unit from pixel ({az_pixel_size:.2f} m) to meter') elif (obs_ds_name == 'rangeOffset') & (stack_obj.metadata['PROCESSOR'] != 'cosicorr'): + # 距离向偏移原本单位是像素，需要乘以距离向像素大小转换为米。 rg_pixel_size = float(stack_obj.metadata['RANGE_PIXEL_SIZE']) rg_pixel_size /= float(stack_obj.metadata['RLOOKS']) ts *= -1 * rg_pixel_size @@ -937,6 +992,8 @@ def run_ifgram_inversion(inps): run_ifgram_inversion(inps) """ + # run_ifgram_inversion() 是整个 ifgram_inversion.py 的主入口。 + # 它负责准备输出文件、把图像切块、逐块调用 run_ifgram_inversion_patch()、再把结果写回磁盘。 start_time = time.time() ## limit the number of threads in numpy/scipy to 1 @@ -963,16 +1020,19 @@ def run_ifgram_inversion(inps): stack_obj = ifgramStack(inps.ifgramStackFile) stack_obj.open(print_msg=False) + # 只读取 dropIfgram=True 的日期对，也就是 modify_network 后保留下来的干涉图。 date12_list = stack_obj.get_date12_list(dropIfgram=True) date_list = stack_obj.get_date_list(dropIfgram=True) length, width = stack_obj.length, stack_obj.width # 1.1 read values on the reference pixel + # refPhase 是每幅干涉图在参考像元上的相位，用于把所有干涉图统一参考到同一点。 inps.refPhase = stack_obj.get_reference_phase(unwDatasetName=inps.obsDatasetName, skip_reference=inps.skip_ref, dropIfgram=True) # 1.2 design matrix + # A 的行数是干涉图数量，列数是日期数-1；它描述每个干涉图由哪些日期位移相减得到。 A = stack_obj.get_design_matrix4timeseries(date12_list)[0] num_pair, num_date = A.shape[0], A.shape[1]+1 inps.numIfgram = num_pair @@ -1022,6 +1082,7 @@ def run_ifgram_inversion(inps): ## 2. prepare output # 2.1 metadata + # 输出文件继承输入 ifgramStack 的元数据，并额外写入本次反演的关键配置。 meta = dict(stack_obj.metadata) for key in config_keys: meta[key_prefix+key] = str(vars(inps)[key]) @@ -1031,6 +1092,7 @@ def run_ifgram_inversion(inps): meta['REF_DATE'] = date_list[0] # 2.2 instantiate time-series + # timeseries.h5 保存反演得到的位移时序，包含 date、bperp 和 timeseries 三个核心数据集。 dates = np.array(date_list, dtype=np.bytes_) pbase = stack_obj.get_perp_baseline_timeseries(dropIfgram=True) ds_name_dict = { @@ -1041,6 +1103,7 @@ def run_ifgram_inversion(inps): writefile.layout_hdf5(inps.tsFile, ds_name_dict, metadata=meta) if inps.calcCov: + # 可选输出：时间序列协方差文件，用来量化每期位移的不确定性。 fbase = os.path.splitext(inps.tsFile)[0] fbase += 'Decor' if inps.obsDatasetName.startswith('unwrapPhase') else '' tsStdFile = f'{fbase}Cov.h5' @@ -1050,6 +1113,7 @@ def run_ifgram_inversion(inps): writefile.layout_hdf5(tsStdFile, ds_name_dict, meta) # 2.3 instantiate invQualifyFile: temporalCoherence / residualInv + # invQualityFile 保存反演质量：相位时通常是 temporalCoherence，偏移时通常是 residual。 if 'residual' in os.path.basename(inps.invQualityFile).lower(): inv_quality_name = 'residual' meta['UNIT'] = 'pixel' @@ -1062,6 +1126,7 @@ def run_ifgram_inversion(inps): writefile.layout_hdf5(inps.invQualityFile, ds_name_dict, metadata=meta) # 2.4 instantiate number of inverted observations + # numInvFile 保存每个像元实际参与反演的干涉图数量，可用于筛选低约束像元。 meta['FILE_TYPE'] = 'mask' meta['UNIT'] = '1' # ignore NO_DATA_VALUE from ifgram stack file here as 1) it makes sense @@ -1075,9 +1140,11 @@ def run_ifgram_inversion(inps): ## 3. run the inversion / estimation and write to disk # 3.1 split ifgram_file into blocks to save memory + # 大图像直接一次性反演会占用大量内存，所以按 maxMemory 把图像切成多个 box。 box_list, num_box = stack_obj.split2boxes(max_memory=inps.maxMemory) # 3.2 prepare the input arguments for *_patch() + # data_kwargs 是传给 run_ifgram_inversion_patch() 的固定参数；循环中只更新 box。 data_kwargs = { "ifgram_file" : inps.ifgramStackFile, "ref_phase" : inps.refPhase, @@ -1104,10 +1171,12 @@ def run_ifgram_inversion(inps): data_kwargs['box'] = box if not inps.cluster: # non-parallel + # 单进程模式：直接调用 patch 函数。 ts, ts_cov, inv_quality, num_inv_obs = run_ifgram_inversion_patch(**data_kwargs)[:-1] else: # parallel + # 并行模式：用 Dask 把同一个 patch 内部的计算再分给多个 worker。 print('\n\n------- start parallel processing using Dask -------') # initiate the output data @@ -1134,6 +1203,7 @@ def run_ifgram_inversion(inps): # write the block to disk # with 3D block in [z0, z1, y0, y1, x0, x1] # and 2D block in [y0, y1, x0, x1] + # HDF5 分块写入时必须指定写入区域 block，避免覆盖其它 patch 的结果。 # time-series - 3D block = [0, num_date, box[1], box[3], box[0], box[2]] writefile.write_hdf5_block(inps.tsFile, @@ -1167,6 +1237,7 @@ def run_ifgram_inversion(inps): # 3.4 update output data on the reference pixel (for phase) if not inps.skip_ref: + # 参考像元理论上位移为 0、质量为最佳；这里显式修正输出文件中的参考像元值。 # grab ref_y/x ref_y = int(stack_obj.metadata['REF_Y']) ref_x = int(stack_obj.metadata['REF_X']) @@ -1183,6 +1254,7 @@ def run_ifgram_inversion(inps): f['mask'][ref_y, ref_x] = num_pair # roll back to the original number of threads + # 恢复进入函数前的线程设置，避免影响后续其它模块。 cluster.roll_back_num_threads(num_threads_dict) m, s = divmod(time.time() - start_time, 60) diff --git a/src/mintpy/load_data.py b/src/mintpy/load_data.py index 1ec38a044..d4a648b74 100644 --- a/src/mintpy/load_data.py +++ b/src/mintpy/load_data.py @@ -5,14 +5,22 @@ ############################################################ +# glob 用来按通配符搜索文件，例如 '/path/*.unw' 会找到所有 .unw 文件。 import glob +# importlib 用来按字符串动态导入模块；这里根据 processor 名称导入 prep_isce/prep_aria 等模块。 import importlib +# os 用来处理路径、当前目录、判断文件是否存在等操作。 import os +# time 用来统计 load_data 整体运行耗时。 import time +# warnings 用来打印警告信息；警告不会立刻中断程序，和 raise Exception 不同。 import warnings +# subset 模块负责从模板中读取裁剪范围，例如按像素范围或经纬度范围裁剪数据。 from mintpy import subset +# auto_path 模块保存不同 InSAR 处理软件的默认路径规则，用于把模板里的 auto 转换成真实路径。 from mintpy.defaults import auto_path +# 这些对象和常量定义了 MintPy HDF5 文件中常见的数据集名称，以及几何/干涉图栈的数据结构。 from mintpy.objects import ( GEOMETRY_DSET_NAMES, IFGRAM_DSET_NAMES, @@ -20,15 +28,23 @@ ifgramStack, sensor, ) +# geometryDict、ifgramDict、ifgramStackDict 是临时字典对象： +# 它们先收集“外部文件路径 + 元数据”，最后统一写成 MintPy 标准 HDF5 文件。 from mintpy.objects.stackDict import geometryDict, ifgramDict, ifgramStackDict +# ptime 处理日期格式；readfile 读取模板/属性/数据；ut 是 MintPy 通用工具函数集合。 from mintpy.utils import ptime, readfile, utils as ut ################################################################# +# PROCESSOR_LIST 列出 load_data 支持的上游处理软件名称。 +# 用户模板中的 mintpy.load.processor 必须是这些名字之一。 PROCESSOR_LIST = ['isce', 'aria', 'hyp3', 'gmtsar', 'snap', 'gamma', 'roipac', 'cosicorr', 'nisar'] # primary observation dataset names +# OBS_DSET_NAMES 是主要观测数据集名称：相位、距离向偏移、方位向偏移。 OBS_DSET_NAMES = ['unwrapPhase', 'rangeOffset', 'azimuthOffset'] +# 下面几个 *_DSET_NAME2TEMPLATE_KEY 字典是“数据集名称 -> 模板配置项”的对应表。 +# 例如 unwrapPhase 在 HDF5 里叫 unwrapPhase，但它的输入文件路径来自模板项 mintpy.load.unwFile。 IFG_DSET_NAME2TEMPLATE_KEY = { 'unwrapPhase' : 'mintpy.load.unwFile', 'coherence' : 'mintpy.load.corFile', @@ -79,27 +95,34 @@ def read_inps2dict(inps): Returns: iDict - dict, input arguments from command line & template file """ # Read input info into iDict + # vars(inps) 把 argparse.Namespace 对象转换成普通字典，方便用 iDict['key'] 读写。 iDict = vars(inps) + # 先给一些关键字段设置默认值，后面模板或项目名识别会覆盖它们。 iDict['PLATFORM'] = None iDict['processor'] = 'isce' # Read template file template = {} for fname in inps.template_file: + # read_template() 把 cfg 文件读成字典；check_template_auto_value() 会处理 yes/no/auto 等值。 temp = readfile.read_template(fname) temp = ut.check_template_auto_value(temp) + # update() 会把当前模板内容合并到 template；后读的模板可以覆盖先读的同名配置。 template.update(temp) for key, value in template.items(): iDict[key] = value # group - load prefix = 'mintpy.load.' + # 找出所有以 mintpy.load. 开头的模板配置，并去掉前缀得到短 key。 key_list = [i.split(prefix)[1] for i in template.keys() if i.startswith(prefix)] for key in key_list: value = template[prefix+key] if key in ['processor', 'autoPath', 'updateMode', 'compression']: + # 这些配置在后面经常使用，所以同时保存成短 key，例如 iDict['processor']。 iDict[key] = template[prefix+key] elif value: + # 其它 load 配置保留完整 key，避免和其它配置冲突。 iDict[prefix+key] = template[prefix+key] print('processor : {}'.format(iDict['processor'])) @@ -117,6 +140,7 @@ def read_inps2dict(inps): # PROJECT_NAME --> PLATFORM if not iDict['PROJECT_NAME']: + # 如果命令行没有给项目名，就尝试从自定义模板文件名中推断项目名和卫星平台。 cfile = [i for i in list(inps.template_file) if os.path.basename(i) != 'smallbaselineApp.cfg'] iDict['PROJECT_NAME'] = sensor.project_name2sensor_name(cfile)[1] @@ -132,6 +156,7 @@ def read_inps2dict(inps): # update file path with auto if iDict.get('autoPath', False): print('use auto path defined in mintpy.defaults.auto_path for options in auto') + # auto_path.get_auto_path() 会把模板里的 auto 路径替换成当前 processor 的默认文件路径模式。 iDict = auto_path.get_auto_path(processor=iDict['processor'], work_dir=os.getcwd(), template=iDict) @@ -148,8 +173,10 @@ def read_subset_box(iDict): the geo bounding box of the box above. """ # Read subset info from template + # box 是用于普通观测/雷达几何文件的裁剪范围；box4geo 是用于地理坐标查找表的裁剪范围。 iDict['box'] = None iDict['box4geo'] = None + # subset.read_subset_template2box() 从模板读取裁剪设置，返回像素框 pix_box 和经纬度框 geo_box。 pix_box, geo_box = subset.read_subset_template2box(iDict['template_file'][0]) # Grab required info to read input geo_box into pix_box @@ -165,6 +192,7 @@ def read_subset_box(iDict): # 2) it is in the same coordinate type as observation files dem_files = glob.glob(iDict['mintpy.load.demFile']) if len(dem_files) > 0: + # 用 DEM 的属性判断输入数据是地理坐标还是雷达坐标；Y_FIRST 是地理坐标文件常见属性。 atr = readfile.read_attribute(dem_files[0]) else: atr = dict() @@ -192,8 +220,10 @@ def read_subset_box(iDict): return iDict # geo_box --> pix_box + # ut.coordinate() 创建坐标转换对象，可以在经纬度框和像素框之间转换。 coord = ut.coordinate(atr, lookup_file=lookup_file) if geo_box is not None: + # bbox_geo2radar() 把经纬度范围转换为雷达坐标下的像素范围。 pix_box = coord.bbox_geo2radar(geo_box) pix_box = coord.check_box_within_data_coverage(pix_box) print(f'input bounding box of interest in lat/lon: {geo_box}') @@ -204,6 +234,7 @@ def read_subset_box(iDict): if lookup_file is not None: atrLut = readfile.read_attribute(lookup_file[0]) if not geocoded and 'Y_FIRST' in atrLut.keys(): + # 对于 Gamma/ROI_PAC 等情况，查找表可能是地理坐标，需要单独计算它的裁剪框。 geo_box = coord.bbox_radar2geo(pix_box) box4geo_lut = ut.coordinate(atrLut).bbox_geo2radar(geo_box) print(f'box to read for geocoded lookup file in y/x: {box4geo_lut}') @@ -221,10 +252,12 @@ def update_box4files_with_inconsistent_size(fnames): Returns: pix_box - None if all files are in same size (0, 0, min_width, min_length) if not. """ + # 读取每个输入干涉图的 LENGTH/WIDTH 属性，检查行列数是否一致。 atr_list = [readfile.read_attribute(fname) for fname in fnames] length_list = [int(atr['LENGTH']) for atr in atr_list] width_list = [int(atr['WIDTH']) for atr in atr_list] if any(len(set(i)) for i in [length_list, width_list]): + # 如果尺寸不一致，选择所有文件都共有的最小范围，避免读取越界。 min_length = min(length_list) min_width = min(width_list) pix_box = (0, 0, min_width, min_length) @@ -254,6 +287,7 @@ def update_box4files_with_inconsistent_size(fnames): def skip_files_with_inconsistent_size(dsPathDict, pix_box=None, dsName='unwrapPhase'): """Skip files by removing the file path from the input dsPathDict.""" + # dsPathDict 保存每类数据集对应的文件路径列表，例如 {'unwrapPhase': [...], 'coherence': [...]}。 atr_list = [readfile.read_attribute(fname) for fname in dsPathDict[dsName]] length_list = [int(atr['LENGTH']) for atr in atr_list] width_list = [int(atr['WIDTH']) for atr in atr_list] @@ -261,6 +295,7 @@ def skip_files_with_inconsistent_size(dsPathDict, pix_box=None, dsName='unwrapPh # Check size requirements drop_inconsistent_files = False if any(len(set(size_list)) > 1 for size_list in [length_list, width_list]): + # 如果没有指定裁剪框，尺寸不同的文件无法放进同一个 HDF5 栈，只能跳过异常尺寸文件。 if pix_box is None: drop_inconsistent_files = True else: @@ -271,6 +306,7 @@ def skip_files_with_inconsistent_size(dsPathDict, pix_box=None, dsName='unwrapPh # update dsPathDict if drop_inconsistent_files: + # most_common() 找出最常见的行列数；非这个尺寸的干涉图会被移除。 common_length = ut.most_common(length_list) common_width = ut.most_common(width_list) @@ -288,6 +324,7 @@ def skip_files_with_inconsistent_size(dsPathDict, pix_box=None, dsName='unwrapPh if length != common_length or width != common_width: dates = ptime.yyyymmdd(date12.split('-')) # update file list for all datasets + # 对某一个异常 date12，要从 unwrapPhase/coherence/connectComponent 等所有数据集中一起移除。 for dsName in dsNames: fnames = [i for i in dsPathDict[dsName] if all(d[2:8] in i for d in dates)] @@ -312,8 +349,10 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key): Returns: stackObj - ifgramStackDict object or None """ if iDict['only_load_geometry']: + # 用户指定 --geom 时，只加载几何文件，不创建干涉图/电离层/偏移栈对象。 return None + # 根据传入的映射表判断这次要加载的是普通干涉图栈、电离层栈还是偏移栈。 if 'mintpy.load.unwFile' in ds_name2template_key.values(): obs_type = 'interferogram' elif 'mintpy.load.ionUnwFile' in ds_name2template_key.values(): @@ -330,12 +369,14 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key): for dsName in [i for i in IFGRAM_DSET_NAMES if i in ds_name2template_key.keys()]: key = ds_name2template_key[dsName] if key in iDict.keys(): + # glob.glob() 按模板路径搜索真实文件；sorted() 保证顺序稳定。 files = sorted(glob.glob(str(iDict[key]))) if len(files) > 0: dsPathDict[dsName] = files print(f'{dsName:<{max_digit}}: {iDict[key]}') # Check 1: required dataset + # 普通干涉图必须有 unwrapPhase；偏移栈必须至少有 rangeOffset/azimuthOffset 之一。 dsName0s = [x for x in OBS_DSET_NAMES if x in ds_name2template_key.keys()] dsName0 = [i for i in dsName0s if i in dsPathDict.keys()] if len(dsName0) == 0: @@ -360,6 +401,7 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key): dsNumList = list(dsNumDict.values()) if any(i != dsNumList[0] for i in dsNumList): + # 不同数据集数量不一致时继续运行，但后面会按 date12 匹配，跳过缺文件的干涉对。 msg = 'WARNING: NOT all types of dataset have the same number of files.' msg += ' -> skip interferograms with missing files and continue.' print(msg) @@ -409,6 +451,7 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key): else: # 2nd guess: any file in the list + # 如果同一索引位置不是对应日期，就在整个文件列表里搜索包含同一对日期的文件。 dsPath2 = [p for p in dsPathDict[dsName] if (all(d6 in p for d6 in date6s) or (date12MJD and date12MJD in dsPath1))] @@ -419,13 +462,16 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key): print(f'WARNING: {dsName:>18} file missing for pair {date6s}') # initiate ifgramDict object + # ifgramDict 表示一个干涉对的所有数据文件路径集合。 ifgramObj = ifgramDict(datasetDict=ifgramPathDict) # update pairsDict object + # pairsDict 的 key 是日期对，例如 ('20180101', '20180113')。 date8s = ptime.yyyymmdd(date6s) pairsDict[tuple(date8s)] = ifgramObj if len(pairsDict) > 0: + # ifgramStackDict 表示整个干涉图栈，后面可以一次性写入 ifgramStack.h5。 stackObj = ifgramStackDict(pairsDict=pairsDict, dsName0=dsName0) else: stackObj = None @@ -441,6 +487,7 @@ def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): """ # eliminate lookup table dsName for input files in radar-coordinates + # 不同上游处理软件输出的查找表坐标系不同，因此需要去掉不适用的数据集名称。 if iDict['processor'] in ['isce', 'doris']: # for processors with lookup table in radar-coordinates, remove azimuth/rangeCoord dset_name2template_key.pop('azimuthCoord') @@ -465,9 +512,11 @@ def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): for dsName in [i for i in GEOMETRY_DSET_NAMES if i in dset_name2template_key.keys()]: key = dset_name2template_key[dsName] if key in iDict.keys(): + # 对每个几何数据集，按模板路径找到真实输入文件。 files = sorted(glob.glob(str(iDict[key]))) if len(files) > 0: if dsName == 'bperp': + # bperp 是垂直基线，通常每个日期一个文件，因此保存成 date -> file 的字典。 bperpDict = {} for file in files: date = ptime.yyyymmdd(os.path.basename(os.path.dirname(file))) @@ -476,6 +525,7 @@ def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): print(f'{dsName:<{max_digit}}: {iDict[key]}') print(f'number of bperp files: {len(list(bperpDict.keys()))}') else: + # 其它几何量通常只需要一个文件，例如 DEM、高程角、阴影掩膜。 dsPathDict[dsName] = files[0] print(f'{dsName:<{max_digit}}: {files[0]}') @@ -486,6 +536,7 @@ def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): # extra metadata from observations # e.g. EARTH_RADIUS, HEIGHT, etc. + # 几何文件可能缺少部分元数据，因此从观测文件中补充一些通用元数据。 obsMetaGeo = None obsMetaRadar = None for obsName in OBS_DSET_NAMES: @@ -499,6 +550,7 @@ def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): break # dsPathDict --> dsGeoPathDict + dsRadarPathDict + # 按文件属性中是否存在 Y_FIRST，把几何文件分成地理坐标和雷达坐标两类。 dsNameList = list(dsPathDict.keys()) dsGeoPathDict = {} dsRadarPathDict = {} @@ -515,11 +567,13 @@ def read_inps_dict2geometry_dict_object(iDict, dset_name2template_key): geomGeoObj = None geomRadarObj = None if len(dsGeoPathDict) > 0: + # geometryDict 是“待写入 geometryGeo.h5 的几何数据集合”。 geomGeoObj = geometryDict( processor=iDict['processor'], datasetDict=dsGeoPathDict, extraMetadata=obsMetaGeo) if len(dsRadarPathDict) > 0: + # geometryDict 是“待写入 geometryRadar.h5 的几何数据集合”。 geomRadarObj = geometryDict( processor=iDict['processor'], datasetDict=dsRadarPathDict, @@ -550,17 +604,21 @@ def run_or_skip(outFile, inObj, box, updateMode=True, xstep=1, ystep=1, geom_obj # skip if there is no dict object to write if not inObj: + # 没有输入对象，说明没有对应数据需要写入。 flag = 'skip' return flag # run if not in update mode if not updateMode: + # updateMode 关闭时，总是重新写入。 return flag if ut.run_or_skip(outFile, readable=True) == 'skip': + # 只有输出文件已经存在且可读时，才进一步检查里面的数据集是否完整。 kwargs = dict(box=box, xstep=xstep, ystep=ystep) if inObj.name == 'ifgramStack': + # 对干涉图栈，要同时检查尺寸、数据集名称和 date12 日期对是否完整。 in_size = inObj.get_size(geom_obj=geom_obj, **kwargs)[1:] in_dset_list = inObj.get_dataset_list() in_date12_list = inObj.get_date12_list() @@ -579,6 +637,7 @@ def run_or_skip(outFile, inObj, box, updateMode=True, xstep=1, ystep=1, geom_obj flag = 'skip' elif inObj.name == 'geometry': + # 对几何文件，只需要检查尺寸和数据集名称是否完整。 in_size = inObj.get_size(**kwargs) in_dset_list = inObj.get_dataset_list() @@ -598,6 +657,8 @@ def run_or_skip(outFile, inObj, box, updateMode=True, xstep=1, ystep=1, geom_obj def prepare_metadata(iDict): """Prepare metadata via prep_{processor}.py scripts.""" + # metadata 是描述数据的重要属性，例如日期、轨道、波长、坐标范围等。 + # 不同处理器的原始产品格式不同，所以这里先调用对应 prep_xxx.py 标准化元数据。 processor = iDict['processor'] script_name = f'prep_{processor}.py' print('-'*50) @@ -609,16 +670,19 @@ def prepare_metadata(iDict): raise ValueError(msg) # import prep_{processor} + # importlib.import_module() 通过字符串动态导入模块，例如 processor='isce' 时导入 mintpy.cli.prep_isce。 prep_module = importlib.import_module(f'mintpy.cli.prep_{processor}') if processor in ['gamma', 'hyp3', 'roipac', 'snap', 'cosicorr']: # run prep_module + # 这些 processor 通常需要对每类输入文件逐个补充/转换元数据。 for key in [i for i in iDict.keys() if (i.startswith('mintpy.load.') and i.endswith('File') and i != 'mintpy.load.metaFile')]: if len(glob.glob(str(iDict[key]))) > 0: # print command line + # iargs 模拟命令行参数，传给 prep_xxx.py 的 main()。 iargs = [iDict[key]] if processor == 'gamma': if iDict['PLATFORM']: @@ -633,6 +697,7 @@ def prepare_metadata(iDict): prep_module.main(iargs) elif processor == 'nisar': + # NISAR 的 GUNW 产品需要 DEM、频率等信息，且这里要求 DEM 必须真实存在。 dem_file = iDict['mintpy.load.demFile'] gunw_files = iDict['mintpy.load.unwFile'] water_mask = iDict['mintpy.load.waterMaskFile'] @@ -655,6 +720,7 @@ def prepare_metadata(iDict): ) # run prep_*.py + # prep_nisar.py 会读取 GUNW 文件并准备 MintPy 后续加载所需的元数据。 iargs = ['-i', gunw_files, '-d', dem_file, '--frequency', frequency] if str(water_mask).lower() not in ['auto', 'none', 'no', ''] and os.path.exists(water_mask): @@ -674,9 +740,11 @@ def prepare_metadata(iDict): try: prep_module.main(iargs) except: + # 这里捕获所有异常并给警告，是为了兼容“元数据已经存在”的情况继续往下走。 warnings.warn('prep_nisar.py failed. Assuming its result exists and continue...') elif processor == 'isce': + # ISCE 产品需要 meta 文件、baseline 目录、geometry 目录以及观测文件路径。 from mintpy.utils import isce_utils, s1_utils # --meta-file @@ -725,6 +793,7 @@ def prepare_metadata(iDict): warnings.warn('prep_isce.py failed. Assuming its result exists and continue...') # [optional] for topsStack: SAFE_files.txt --> S1A/B_date.txt + # Sentinel-1 TOPS 数据有时需要 SAFE_files.txt 来生成 S1A/B 日期列表，用于后续范围偏差修正等。 if os.path.isfile(meta_file) and isce_utils.get_processor(meta_file) == 'topsStack': safe_list_file = os.path.join(os.path.dirname(os.path.dirname(meta_file)), 'SAFE_files.txt') if os.path.isfile(safe_list_file): @@ -736,6 +805,7 @@ def prepare_metadata(iDict): elif processor == 'aria': ## compose input arguments # use the default template file if exists & input + # ARIA 的 prep 阶段会直接生成部分 MintPy 输入文件，因此后面 load_data() 会跳过重复写入。 default_temp_files = [fname for fname in iDict['template_file'] if fname.endswith('smallbaselineApp.cfg')] if len(default_temp_files) > 0: @@ -745,6 +815,7 @@ def prepare_metadata(iDict): iargs = ['--template', temp_file] # file name/dir/path + # ARG2OPT_DICT 把 prep_aria.py 的命令行参数名映射到 MintPy 模板配置项。 ARG2OPT_DICT = { '--stack-dir' : 'mintpy.load.unwFile', '--unwrap-stack-name' : 'mintpy.load.unwFile', @@ -760,6 +831,7 @@ def prepare_metadata(iDict): for arg_name, opt_name in ARG2OPT_DICT.items(): arg_value = iDict.get(opt_name, 'auto') if arg_value.lower() not in ['auto', 'no', 'none']: + # 对目录、文件名、普通文件路径三类参数分别处理。 if arg_name.endswith('dir'): iargs += [arg_name, os.path.dirname(arg_value)] elif arg_name.endswith('name'): @@ -779,6 +851,7 @@ def prepare_metadata(iDict): elif processor == 'gmtsar': # use the custom template file if exists & input + # GMTSAR 的 prep 脚本依赖自定义模板；只有默认模板时信息不够。 custom_temp_files = [fname for fname in iDict['template_file'] if not fname.endswith('smallbaselineApp.cfg')] if len(custom_temp_files) == 0: @@ -804,6 +877,7 @@ def get_extra_metadata(iDict): """ extraDict = {} # all keys in MACRO_CASE + # 约定：全大写 key 通常是要写入输出 HDF5 的额外元数据。 upper_keys = [i for i in iDict.keys() if i.isupper()] for key in upper_keys: value = iDict[key] @@ -821,14 +895,17 @@ def load_data(inps): ## 0. read input start_time = time.time() + # 第一步：把命令行参数和模板文件合并成一个总字典 iDict。 iDict = read_inps2dict(inps) ## 1. prepare metadata + # 第二步：调用 prep_xxx.py 为不同处理器产品准备/补全元数据。 prepare_metadata(iDict) extraDict = get_extra_metadata(iDict) # skip data writing as it is included in prep_aria/nisar if iDict['processor'] in ['aria', 'nisar']: + # ARIA/NISAR 的准备脚本已经生成 MintPy 可用文件，因此这里直接返回。 return ## 2. search & write data files @@ -840,15 +917,18 @@ def load_data(inps): kwargs = dict(updateMode=iDict['updateMode'], xstep=iDict['xstep'], ystep=iDict['ystep']) # read subset info [need the metadata from above] + # 第三步：读取裁剪范围，决定后面只加载全图还是子区域。 iDict = read_subset_box(iDict) # geometry in geo / radar coordinates + # 第四步：准备几何文件。这里把几何数据集和观测数据集的模板映射合并，是为了能从观测文件补充尺寸/坐标信息。 geom_dset_name2template_key = { **GEOM_DSET_NAME2TEMPLATE_KEY, **IFG_DSET_NAME2TEMPLATE_KEY, **OFF_DSET_NAME2TEMPLATE_KEY, } geom_geo_obj, geom_radar_obj = read_inps_dict2geometry_dict_object(iDict, geom_dset_name2template_key) + # 几何输出统一写到 ./inputs 目录。 geom_geo_file = os.path.abspath('./inputs/geometryGeo.h5') geom_radar_file = os.path.abspath('./inputs/geometryRadar.h5') @@ -856,6 +936,7 @@ def load_data(inps): compression = 'lzf' if iDict['compression'] == 'default' else iDict['compression'] if run_or_skip(geom_geo_file, geom_geo_obj, iDict['box4geo'], **kwargs) == 'run': + # write2hdf5() 是 geometryDict 对象的方法，会把外部几何文件读取并写成标准 HDF5。 geom_geo_obj.write2hdf5( outputFile=geom_geo_file, access_mode='w', @@ -865,6 +946,7 @@ def load_data(inps): compression=compression) if run_or_skip(geom_radar_file, geom_radar_obj, iDict['box'], **kwargs) == 'run': + # radar 几何文件会额外写入 PROJECT_NAME、PLATFORM 等元数据。 geom_radar_obj.write2hdf5( outputFile=geom_radar_file, access_mode='w', @@ -876,6 +958,7 @@ def load_data(inps): # observations: ifgram, ion or offset # loop over obs stacks + # 第五步：依次加载普通干涉图栈、电离层栈、偏移栈。 stack_ds_name2tmpl_key_list = [ IFG_DSET_NAME2TEMPLATE_KEY, ION_DSET_NAME2TEMPLATE_KEY, @@ -887,15 +970,18 @@ def load_data(inps): for ds_name2tmpl_opt, stack_file in zip(stack_ds_name2tmpl_key_list, stack_files): # initiate dict objects + # 根据模板路径搜索文件，并构建 ifgramStackDict 对象。 stack_obj = read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2tmpl_opt) # use geom_obj as size reference while loading ionosphere geom_obj = None if os.path.basename(stack_file).startswith('ion'): + # 电离层栈的尺寸需要参考几何对象，避免和主观测栈坐标/尺寸不一致。 geom_obj = geom_geo_obj if iDict['geocoded'] else geom_radar_obj # write dict objects to HDF5 files if run_or_skip(stack_file, stack_obj, iDict['box'], geom_obj=geom_obj, **kwargs) == 'run': + # write2hdf5() 会真正读取所有外部干涉图文件，写入 ifgramStack.h5/ionStack.h5/offsetStack.h5。 stack_obj.write2hdf5( outputFile=stack_file, access_mode='w', diff --git a/src/mintpy/modify_network.py b/src/mintpy/modify_network.py index 7e603c2b9..dd08db6af 100644 --- a/src/mintpy/modify_network.py +++ b/src/mintpy/modify_network.py @@ -5,25 +5,36 @@ ############################################################ +# h5py 用来直接读写 HDF5 文件；这里会直接修改 ifgramStack.h5 里的 dropIfgram 数据集。 import h5py +# numpy 用来做数组计算和布尔筛选，例如按阈值选出要删除的干涉图。 import numpy as np +# matplotlib 用来交互式画图；manual_select_pairs_to_remove() 里用户可以点击网络图手动删边。 from matplotlib import dates as mdates, pyplot as plt +# ifgramStack 是 MintPy 表示干涉图栈 HDF5 文件的对象。 from mintpy.objects import ifgramStack +# network(pnet) 负责干涉网算法；plot(pp) 负责绘图；ptime 处理日期；ut 是通用工具函数。 from mintpy.utils import network as pnet, plot as pp, ptime, utils as ut ########################### Sub Function ############################# def reset_network(stackFile): """Reset/restore all pairs within the input file by set all DROP_IFGRAM=no""" + # 这个函数用于“恢复网络”：把所有干涉图都标记为保留。 + # stackFile 通常是 ./inputs/ifgramStack.h5。 print("reset dataset 'dropIfgram' to True for all interferograms for file: "+stackFile) + # 创建 ifgramStack 对象，并打开 HDF5 文件读取元数据和 dropIfgram 数组。 obj = ifgramStack(stackFile) obj.open(print_msg=False) if np.all(obj.dropIfgram): print('All dropIfgram are already True, no need to reset.') else: + # h5py.File(..., 'r+') 表示以可读写模式打开 HDF5 文件。 with h5py.File(stackFile, 'r+') as f: + # dropIfgram 是布尔数组；True 表示该干涉图参与后续反演，False 表示丢弃。 f['dropIfgram'][:] = True + # touch() 更新 coherenceSpatialAvg.txt 的修改时间，让后续步骤知道网络状态变了。 ut.touch('coherenceSpatialAvg.txt') return stackFile @@ -36,7 +47,9 @@ def nearest_neighbor(x, y, x_array, y_array): Output: idx : int, index of min distance - nearest neighbour """ + # 计算鼠标点击点 (x, y) 和所有网络节点 (x_array, y_array) 的欧氏距离。 dist = np.sqrt((x_array - x)**2 + (y_array - y)**2) + # np.argmin() 返回最小距离所在位置，也就是离鼠标点击最近的日期节点。 idx = np.argmin(dist) #idx = dist==np.min(dist) return idx @@ -44,6 +57,7 @@ def nearest_neighbor(x, y, x_array, y_array): def manual_select_pairs_to_remove(stackFile): """Manually select interferograms to remove""" + # 这个函数会打开交互式网络图，让用户用鼠标点击两个日期来选择一条要删除的干涉边。 print('\n-------------------------------------------------------------') print('Manually select interferograms to remove') print('1) click two dates/points to select one pair of interferogram') @@ -52,9 +66,12 @@ def manual_select_pairs_to_remove(stackFile): print('-------------------------------------------------------------\n') obj = ifgramStack(stackFile) obj.open() + # date12ListAll 是所有干涉图日期对，例如 20180101_20180113。 date12ListAll = obj.date12List + # pbase 是每个日期对应的垂直基线时间序列，用来绘制干涉网络纵轴。 pbase = obj.get_perp_baseline_timeseries(dropIfgram=False) dateList = obj.dateList + # matplotlib 内部用浮点数表示日期，date2num() 把日期转换成可画图的数字。 datesNum = mdates.date2num(np.array(ptime.date_list2vector(dateList)[0])) date12ListKept = obj.get_date12_list(dropIfgram=True) @@ -69,10 +86,12 @@ def manual_select_pairs_to_remove(stackFile): date12_click = [] def onclick(event): + # onclick 是鼠标点击事件回调函数：用户每点一次图，matplotlib 会调用它。 idx = nearest_neighbor(event.xdata, event.ydata, datesNum, pbase) print('click at '+dateList[idx]) date_click.append(dateList[idx]) if len(date_click) % 2 == 0 and date_click[-2] != date_click[-1]: + # 每点击两个不同日期，就组成一个干涉图日期对。 [mDate, sDate] = sorted(date_click[-2:]) mIdx = dateList.index(mDate) sIdx = dateList.index(sDate) @@ -80,14 +99,17 @@ def onclick(event): if date12 in date12ListAll: print('select date12: '+date12) date12_click.append(date12) + # 在图上用红色粗线标出用户选择删除的干涉边。 ax.plot([datesNum[mIdx], datesNum[sIdx]], [pbase[mIdx], pbase[sIdx]], 'r', lw=4) else: print(date12+' does not exist in input file') plt.draw() fig.canvas.mpl_connect('button_press_event', onclick) + # plt.show() 打开图形窗口；用户关闭窗口后代码继续往下执行。 plt.show() + # yes_or_no() 会询问用户是否确认删除刚才选中的干涉图。 if not ut.yes_or_no('Proceed to drop the ifgrams/date12?'): date12_click = None @@ -96,15 +118,18 @@ def onclick(event): def get_aoi_pix_box(meta, lookup_file, pix_box, geo_box): """Get pix_box for AOI.""" + # AOI 是 area of interest，感兴趣区域。这个函数把用户输入的区域统一转换成像素框。 coord = ut.coordinate(meta, lookup_file=lookup_file) # geo_box -> pix_box if geo_box and lookup_file: + # 如果用户输入的是经纬度范围，并且有查找表，就转换成雷达坐标像素范围。 print(f'input AOI in (lon0, lat1, lon1, lat0): {geo_box}') pix_box = coord.bbox_geo2radar(geo_box) # check pix_box if pix_box: + # check_box_within_data_coverage() 保证像素框不会超出数据覆盖范围。 pix_box = coord.check_box_within_data_coverage(pix_box) print(f'input AOI in (x0,y0,x1,y1): {pix_box}') @@ -113,6 +138,8 @@ def get_aoi_pix_box(meta, lookup_file, pix_box, geo_box): def get_mst_date12(keep_mst, par_list_all, date12_list_all, date12_to_drop, min_par, par_name='average coherence'): """Get the date12_list of the MST network for the given parameter.""" + # MST 是 minimum spanning tree，最小生成树。 + # 在干涉网里保留 MST 可以避免网络被删断，保证时间序列反演仍有基本连通性。 if keep_mst: print(f'Get minimum spanning tree (MST) of interferograms with inverse of {par_name}.') msg = ('Drop ifgrams with ' @@ -120,11 +147,13 @@ def get_mst_date12(keep_mst, par_list_all, date12_list_all, date12_to_drop, min_ '2) not in MST network: '.format(par_name, min_par)) # get the current remaining network (after all the above criteria and before data-driven) + # 先从全部日期对里去掉前面规则已经决定删除的日期对。 date12_to_keep = sorted(list(set(date12_list_all) - set(date12_to_drop))) par_to_keep = [par for par, date12 in zip(par_list_all, date12_list_all) if date12 in date12_to_keep] # get MST from the current remaining network + # threshold_coherence_based_mst() 根据参数反比构建 MST，返回需要保留的日期对。 mst_date12_list = pnet.threshold_coherence_based_mst(date12_to_keep, par_to_keep) mst_date12_list = ptime.yyyymmdd_date12(mst_date12_list) @@ -140,6 +169,7 @@ def get_date12_to_drop(inps): Return [] if no ifgram to drop, thus KEEP ALL ifgrams; None if nothing to change, exit without doing anything. """ + # 这是本文件最核心的函数：综合所有筛选规则，计算哪些 date12 干涉图要删除。 obj = ifgramStack(inps.file) obj.open() date12ListAll = obj.date12List @@ -147,10 +177,12 @@ def get_date12_to_drop(inps): print(f'number of interferograms: {len(date12ListAll)}') # Get date12_to_drop + # date12_to_drop 是待删除干涉图日期对列表，后面每条规则都会往里面追加。 date12_to_drop = [] # reference file if inps.referenceFile: + # 如果提供参考文件，就只保留参考文件中仍然保留的干涉图。 date12_to_keep = pnet.get_date12_list(inps.referenceFile, dropIfgram=True) print('--------------------------------------------------') print(f'use reference pairs info from file: {inps.referenceFile}') @@ -161,6 +193,7 @@ def get_date12_to_drop(inps): # temp baseline threshold if inps.tempBaseMax: + # tbaseIfgram 是时间基线；超过阈值的干涉图会被删除。 tempIndex = np.abs(obj.tbaseIfgram) > inps.tempBaseMax tempList = list(np.array(date12ListAll)[tempIndex]) date12_to_drop += tempList @@ -170,6 +203,7 @@ def get_date12_to_drop(inps): # perp baseline threshold if inps.perpBaseMax: + # pbaseIfgram 是垂直空间基线；超过阈值的干涉图会被删除。 tempIndex = np.abs(obj.pbaseIfgram) > inps.perpBaseMax tempList = list(np.array(date12ListAll)[tempIndex]) date12_to_drop += tempList @@ -179,6 +213,7 @@ def get_date12_to_drop(inps): # connection number threshold if inps.connNumMax: + # connNumMax 控制每个日期最多连接多少个相邻日期，常用于构建较稀疏的顺序网络。 seq_date12_list = pnet.select_pairs_sequential(dateList, inps.connNumMax) seq_date12_list = ptime.yyyymmdd_date12(seq_date12_list) tempList = [i for i in date12ListAll if i not in seq_date12_list] @@ -192,6 +227,7 @@ def get_date12_to_drop(inps): # excludeIfgIndex if inps.excludeIfgIndex: + # excludeIfgIndex 按干涉图在文件里的索引编号删除。 tempList = [date12ListAll[i] for i in inps.excludeIfgIndex] date12_to_drop += tempList print('--------------------------------------------------') @@ -201,6 +237,7 @@ def get_date12_to_drop(inps): # excludeDate12 if inps.excludeDate12: + # excludeDate12 直接按日期对字符串删除，例如 20180101_20180113。 tempList = [i for i in inps.excludeDate12 if i in date12ListAll] date12_to_drop += tempList print('--------------------------------------------------') @@ -210,6 +247,7 @@ def get_date12_to_drop(inps): # excludeDate if inps.excludeDate: + # excludeDate 删除所有包含指定单景日期的干涉图。 tempList = [i for i in date12ListAll if any(j in inps.excludeDate for j in i.split('_'))] date12_to_drop += tempList print('-'*50+'\nDrop ifgrams including the following dates: ({})\n{}'.format( @@ -218,6 +256,7 @@ def get_date12_to_drop(inps): # startDate if inps.startDate: + # startDate 删除早于指定日期的干涉图。 minDate = int(inps.startDate) tempList = [i for i in date12ListAll if any(int(j) < minDate for j in i.split('_'))] date12_to_drop += tempList @@ -227,6 +266,7 @@ def get_date12_to_drop(inps): # endDate if inps.endDate: + # endDate 删除晚于指定日期的干涉图。 maxDate = int(inps.endDate) tempList = [i for i in date12ListAll if any(int(j) > maxDate for j in i.split('_'))] date12_to_drop += tempList @@ -236,6 +276,7 @@ def get_date12_to_drop(inps): # coherence file if inps.coherenceBased: + # coherenceBased 根据平均相干性筛选低质量干涉图。 print('--------------------------------------------------') print('use coherence-based network modification') @@ -243,6 +284,7 @@ def get_date12_to_drop(inps): pix_box = get_aoi_pix_box(obj.metadata, inps.lookupFile, inps.aoiYX, inps.aoiLALO) # calculate spatial average coherence + # ut.spatial_average() 计算每个干涉图在 AOI 或掩膜区域内的空间平均值。 cohList = ut.spatial_average(inps.file, datasetName='coherence', maskFile=inps.maskFile, @@ -250,6 +292,7 @@ def get_date12_to_drop(inps): saveList=True)[0] # get coherence-based network + # 只保留平均相干性大于等于阈值的日期对。 coh_date12_list = list(np.array(date12ListAll)[np.array(cohList) >= inps.minCoherence]) # get MST network @@ -268,6 +311,7 @@ def get_date12_to_drop(inps): # area ratio file if inps.areaRatioBased: + # areaRatioBased 根据“高相干像元面积比例”筛选干涉图。 print('--------------------------------------------------') print('use area-ratio-based network modification') @@ -275,6 +319,7 @@ def get_date12_to_drop(inps): pix_box = get_aoi_pix_box(obj.metadata, inps.lookupFile, inps.aoiYX, inps.aoiLALO) # calculate average coherence in masked out areas as threshold + # reverseMask=True 表示在掩膜外区域计算平均相干性，作为动态阈值参考。 meanMaskCoh = np.nanmean(ut.spatial_average(inps.file, datasetName='coherence', maskFile=inps.maskFile, @@ -283,6 +328,7 @@ def get_date12_to_drop(inps): print(f'Average coherence of {inps.maskFile} reverse is {meanMaskCoh:.2f}') # calculate area-ratio with pixels greater than meanMaskCoh + # threshold=meanMaskCoh 表示统计大于这个阈值的像元比例。 areaRatioList = ut.spatial_average(inps.file, datasetName='coherence', maskFile=inps.maskFile, @@ -310,6 +356,7 @@ def get_date12_to_drop(inps): # Manually drop pairs if inps.manual: + # manual=True 时进入交互式手动选边删除流程。 tempList = manual_select_pairs_to_remove(inps.file) if tempList is None: return None @@ -319,6 +366,7 @@ def get_date12_to_drop(inps): ## summary # drop duplicate date12 and sort in order + # set() 去重，sorted() 排序，得到最终稳定的删除列表。 date12_to_drop = sorted(list(set(date12_to_drop))) date12_to_keep = sorted(list(set(date12ListAll) - set(date12_to_drop))) print('--------------------------------------------------') @@ -335,6 +383,7 @@ def get_date12_to_drop(inps): # checking: # 1) no new date12 to drop against existing file # 2) no date12 left after dropping + # 如果计算结果和当前文件已有 dropIfgram 状态完全一致，就不需要重复写文件。 date12ListKept = obj.get_date12_list(dropIfgram=True) date12ListDropped = sorted(list(set(date12ListAll) - set(date12ListKept))) if date12_to_drop == date12ListDropped: @@ -349,18 +398,23 @@ def get_date12_to_drop(inps): ######################################################################## def modify_network(inps): """Run network modification.""" + # modify_network() 是核心入口：计算要删除的干涉图，并写回 ifgramStack.h5 的 dropIfgram。 if inps.reset: + # reset 模式直接恢复所有干涉图，不执行其它筛选规则。 print('--------------------------------------------------') reset_network(inps.file) return # get the list of date12 to drop/exclude + # get_date12_to_drop() 返回 None 表示没有变化；返回 [] 表示保留全部；返回列表表示要删除这些日期对。 inps.date12_to_drop = get_date12_to_drop(inps) # update date dataset in the ifgram stack h5 file if inps.date12_to_drop is not None: + # update_drop_ifgram() 会把指定日期对的 dropIfgram 标记改为 False。 ifgramStack(inps.file).update_drop_ifgram(inps.date12_to_drop) + # 更新辅助文件时间戳，提示后续网络图/平均相干性可能需要重新生成。 ut.touch('coherenceSpatialAvg.txt') print('Done.') diff --git a/src/mintpy/smallbaselineApp.py b/src/mintpy/smallbaselineApp.py index 9e607093a..b6e64721e 100644 --- a/src/mintpy/smallbaselineApp.py +++ b/src/mintpy/smallbaselineApp.py @@ -5,17 +5,27 @@ ############################################################ +# glob 是 Python 标准库，用来按照通配符查找文件，例如 '*.png' 表示所有 png 图片。 import glob +# os 是 Python 标准库，用来和操作系统交互，例如拼接路径、创建目录、切换目录、读取环境变量。 import os +# shutil 是 Python 标准库，用来做高级文件操作，例如复制文件 copy2()、移动文件 move()。 import shutil +# time 是 Python 标准库，用来记录时间；这里用于统计某个处理步骤花了多久。 import time +# numpy 是常用的科学计算库；这里用 np.sum() 等函数对数组像元进行统计。 import numpy as np +# mintpy 是当前项目的主包；mintpy.__file__ 可以告诉我们 MintPy 安装目录在哪里。 import mintpy +# RAMP_LIST 是允许的相位坡度类型列表；cluster 管理并行计算；sensor 处理卫星传感器名称。 from mintpy.objects import RAMP_LIST, cluster, sensor +# readfile/writefile 负责读写 HDF5/模板等文件；utils as ut 是把 utils 模块起短名 ut，方便后面调用。 from mintpy.utils import readfile, utils as ut, writefile +# 说明：本文件负责把 MintPy 小基线时序分析的各个命令行模块串成完整工作流。 +# “工作流”可以理解为：按固定顺序调用 load_data、modify_network、ifgram_inversion 等工具。 # comment out the workflow import in favor of individual imports for speed and robustness # import mintpy.workflow # dynamic import of modules for smallbaselineApp @@ -25,24 +35,37 @@ def get_the_latest_default_template_file(work_dir): """Get the latest version of default template file. If an obsolete file exists in the working directory, the existing option values are kept. """ + # 获取软件自带的最新默认模板，以及当前工作目录中的默认模板路径。 # get the latest and current versions of the default template files + # os.path.dirname(mintpy.__file__) 得到 mintpy 包所在目录。 + # os.path.join(...) 用系统认可的路径分隔符拼出 defaults/smallbaselineApp.cfg 的完整路径。 lfile = os.path.join(os.path.dirname(mintpy.__file__), 'defaults/smallbaselineApp.cfg') + # work_dir 是用户要运行项目的目录；cfile 是这个目录里的默认配置文件路径。 cfile = os.path.join(work_dir, 'smallbaselineApp.cfg') + # os.path.isfile(cfile) 判断 cfile 这个路径是否存在且确实是一个文件。 if not os.path.isfile(cfile): + # 工作目录中没有模板时，直接复制一份默认模板。 print(f'copy default template file {lfile} to work directory') + # shutil.copy2(src, dst) 复制文件，同时尽量保留原文件的时间戳等元信息。 shutil.copy2(lfile, work_dir) else: + # 工作目录已有模板时，检查是否缺少新版模板中的配置项。 #cfile is obsolete if any key is missing + # readfile.read_template() 是 MintPy 的函数，会把 cfg 模板读成 Python 字典 dict。 ldict = readfile.read_template(lfile) cdict = readfile.read_template(cfile) + # any([...]) 表示只要列表里有一个 True，整体就是 True；这里用来判断旧模板是否缺 key。 if any([key not in cdict.keys() for key in ldict.keys()]): print('obsolete default template detected, update to the latest version.') shutil.copy2(lfile, work_dir) + # 更新模板结构后，保留旧模板中用户已经设置过的参数值。 #keep the existing option value from obsolete template file + # ut.update_template_file() 会把 cdict 中已有的用户配置写回新版模板文件。 ut.update_template_file(cfile, cdict) + # 返回当前工作目录里最终可用的模板文件路径。 return cfile @@ -52,8 +75,13 @@ class TimeSeriesAnalysis: """ def __init__(self, customTemplateFile=None, workDir=None): + # 记录用户模板、工作目录以及启动程序时所在目录，便于流程结束后返回。 + # __init__ 是 Python 类的初始化方法；创建 TimeSeriesAnalysis(...) 对象时会自动运行。 + # self 表示“当前这个对象本身”；self.xxx 是把变量保存到对象里，供其它方法继续使用。 self.customTemplateFile = customTemplateFile + # os.path.abspath() 把相对路径转换为绝对路径，避免后续切换目录后找不到文件。 self.workDir = os.path.abspath(workDir) + # os.getcwd() 表示 get current working directory，即获取当前终端所在目录。 self.cwd = os.path.abspath(os.getcwd()) def open(self): @@ -64,29 +92,39 @@ def open(self): """ #1. Get projectName + # projectName 默认没有；如果用户提供了模板，就用模板文件名当项目名。 self.projectName = None if self.customTemplateFile: + # os.path.basename() 取文件名，os.path.splitext() 把文件名和扩展名分开。 self.projectName = os.path.splitext(os.path.basename(self.customTemplateFile))[0] print('Project name:', self.projectName) #2. Go to work directory + # os.makedirs(..., exist_ok=True) 创建目录；目录已存在时不会报错。 os.makedirs(self.workDir, exist_ok=True) + # os.chdir() 切换 Python 当前工作目录；后面相对路径都会以 self.workDir 为起点。 os.chdir(self.workDir) print('Go to work directory:', self.workDir) #3. Read templates + # 调用本文件上面定义的函数，确保工作目录里有最新版 smallbaselineApp.cfg。 self.templateFile = get_the_latest_default_template_file(self.workDir) + # 调用当前类自己的 _read_template() 方法，读取并合并模板参数。 self._read_template() def _read_template(self): + # 读取并合并用户模板与默认模板，是后续所有处理步骤的参数来源。 + # 方法名前面的下划线 _ 表示这个方法主要给类内部使用，不是给用户直接调用的入口。 # 1) update default template self.customTemplate = None if self.customTemplateFile: # customTemplateFile --> customTemplate print('read custom template file:', self.customTemplateFile) + # 把用户模板文件读成字典，例如 {'mintpy.geocode': True, ...}。 cdict = readfile.read_template(self.customTemplateFile) + # 将常见的简写/布尔写法规范化，减少模板输入格式差异造成的问题。 # correct some loose type errors standardValues = { 'def':'auto', 'default':'auto', @@ -94,43 +132,56 @@ def _read_template(self): 'n':'no', 'off':'no', 'false':'no', } for key, value in cdict.items(): + # cdict.items() 会逐个返回字典中的 key 和 value。 if value in standardValues.keys(): cdict[key] = standardValues[value] for key in ['mintpy.deramp', 'mintpy.troposphericDelay.method']: if key in cdict.keys(): + # lower() 转小写，replace('-', '_') 把连字符替换成下划线，统一方法名格式。 cdict[key] = cdict[key].lower().replace('-', '_') + # 兼容旧模板里直接使用 processor 的写法。 if 'processor' in cdict.keys(): cdict['mintpy.load.processor'] = cdict['processor'] + # 以下偏移元数据只在 load_data.py 中临时使用，合并模板后不再需要保留。 # these metadata are used in load_data.py only, not needed afterwards # (in order to manually add extra offset when the lookup table is shifted) # (seen in ROI_PAC product sometimes) for key in ['SUBSET_XMIN', 'SUBSET_YMIN']: if key in cdict.keys(): + # pop() 会从字典中删除这个 key，并返回被删除的值。 cdict.pop(key) + # dict(cdict) 复制一份字典，保存为对象属性，后面写元数据时还会用到。 self.customTemplate = dict(cdict) + # 用用户模板中的配置覆盖默认模板对应项，生成本次运行使用的模板文件。 # customTemplate --> templateFile print('update default template based on input custom template') self.templateFile = ut.update_template_file(self.templateFile, self.customTemplate) + # 将模板备份到 inputs 和 pic 目录，便于追踪本次处理使用的参数。 # 2) backup custom/default template file in inputs/pic folder flen = len(os.path.basename(self.templateFile)) if self.customTemplateFile: + # max(a, b) 取较大值；这里用于对齐打印信息中的文件名宽度。 flen = max(flen, len(os.path.basename(self.customTemplateFile))) for backup_dirname in ['inputs', 'pic']: + # 循环创建 inputs 和 pic 两个备份目录。 backup_dir = os.path.join(self.workDir, backup_dirname) # create directory os.makedirs(backup_dir, exist_ok=True) # back up to the directory + # 列表推导式 [x for x in ... if x] 会过滤掉 None 或空字符串。 tfiles = [x for x in [self.customTemplateFile, self.templateFile] if x] for tfile in tfiles: + # out_file 是备份后的目标文件路径。 out_file = os.path.join(backup_dir, os.path.basename(tfile)) + # ut.run_or_skip() 是 MintPy 的工具函数，用来判断目标文件是否需要重新生成/复制。 if ut.run_or_skip(out_file, in_file=tfile, readable=False, print_msg=False) == 'run': shutil.copy2(tfile, backup_dir) print('copy {f:<{l}} to {d:<8} directory for backup.'.format( @@ -138,9 +189,12 @@ def _read_template(self): # 3) read default template file print('read default template file:', self.templateFile) + # self.template 保存最终生效的完整模板配置，后面每个处理步骤都会读取它。 self.template = readfile.read_template(self.templateFile) + # check_template_auto_value() 会把模板里的 auto/yes/no 等文本进一步转换成程序需要的值。 self.template = ut.check_template_auto_value(self.template) + # 若用户要求输出 HDF-EOS5/KMZ，则需要先开启地理编码。 # correct some loose setup conflicts if self.template['mintpy.geocode'] is False: for key in ['mintpy.save.hdfEos5', 'mintpy.save.kmz']: @@ -157,91 +211,123 @@ def run_load_data(self, step_name): 3) check loading result 4) add custom metadata (optional, for HDF-EOS5 format only) """ + # load_data 是工作流入口步骤：把外部 InSAR 产品整理成 MintPy HDF5 输入。 # 1) copy aux files (optional) + # self._copy_aux_file() 调用当前类中的辅助方法，尝试复制一些可选的外部辅助文件。 self._copy_aux_file() # 2) loading data # compose list of input arguments # instead of using command line then split # to support path with whitespace + # iargs 是“参数列表”，等价于用户在终端里输入 load_data.py 后面的那些参数。 + # 用列表而不是一整串字符串，可以正确处理路径中包含空格的情况。 iargs = ['--template', self.templateFile] if self.customTemplateFile: + # 如果用户提供了自定义模板，也把它传给 load_data.py。 iargs += [self.customTemplateFile] if self.projectName: + # 如果能从模板名得到项目名，就通过 --project 参数传给 load_data.py。 iargs += ['--project', self.projectName] # run command line print('\nload_data.py', ' '.join(iargs)) + # import mintpy.cli.load_data 表示导入 MintPy 里的 load_data 命令模块。 + # 这里不是开启新的终端进程，而是在 Python 代码里直接调用它的 main() 函数。 import mintpy.cli.load_data + # main(iargs) 会按 iargs 里的参数执行 load_data 的功能。 mintpy.cli.load_data.main(iargs) # come back to working directory + # load_data 执行过程中可能切换目录，所以这里强制切回本次项目工作目录。 os.chdir(self.workDir) # 3) check loading result + # ut.check_loaded_dataset() 会检查 inputs 目录中的关键 HDF5 文件是否已经正确生成。 + # [:4] 表示只取返回结果的前 4 个：干涉图栈、几何文件、查找表、电离层文件。 stack_file, geom_file, _, ion_file = ut.check_loaded_dataset(self.workDir, print_msg=True)[:4] + # 若提供了自定义模板，把其中的元数据补写到加载后的数据文件中。 # 4) add custom metadata (optional) if self.customTemplateFile: # use ut.add_attribute() instead of add_attribute.py because of # better control of special metadata, such as SUBSET_X/YMIN + # msg 是日志文字，f'...' 是 Python 的格式化字符串，可以把变量值嵌入文本。 msg = f'updating metadata based on custom template file {os.path.basename(self.customTemplateFile)}' for fname in [stack_file, ion_file]: #, geom_file]: if fname: print(f'{msg} for file: {os.path.basename(fname)}') + # ut.add_attribute() 会把 self.customTemplate 中的键值对写入 HDF5 文件属性。 ut.add_attribute(fname, self.customTemplate) def _copy_aux_file(self): + # 针对 University of Miami 环境的辅助文件拷贝逻辑，普通运行环境通常不会触发。 # for Univ of Miami + # os.getenv('SCRATCHDIR') 用来读取名为 SCRATCHDIR 的环境变量；没有设置时返回 None。 if os.getenv('SCRATCHDIR') and self.projectName: + # proj_dir 是外部临时/共享目录中当前项目的位置。 proj_dir = os.path.join(os.getenv('SCRATCHDIR'), self.projectName) flist = ['PROCESS/unavco_attributes.txt', 'PROCESS/bl_list.txt', 'SLC/summary*slc.jpg'] + # ut.get_file_list() 会把通配符路径展开成真实文件列表；abspath=True 表示返回绝对路径。 flist = ut.get_file_list([os.path.join(proj_dir, i) for i in flist], abspath=True) for fname in flist: if ut.run_or_skip(os.path.basename(fname), in_file=fname, readable=False) == 'run': + # 把辅助文件复制到当前工作目录。 shutil.copy2(fname, self.workDir) print(f'copy {os.path.basename(fname)} to work directory') def run_network_modification(self, step_name): """Modify network of interferograms before the network inversion.""" + # 在反演前裁剪或调整干涉网络，并准备网络图所需的输入。 # check the existence of ifgramStack.h5 + # stack_file 是干涉图栈文件；geom_file 是几何信息文件。 stack_file, geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[:2] + # coherenceSpatialAvg.txt 是空间平均相干性文本文件，plot_network 会用到。 coh_txt = os.path.join(self.workDir, 'coherenceSpatialAvg.txt') + # net_fig 同时检查工作目录和 pic 目录下是否已有 network.pdf。 net_fig = [os.path.join(self.workDir, i, 'network.pdf') for i in ['', 'pic']] net_fig = [x for x in net_fig if os.path.isfile(x)] # 1) output waterMask.h5 to simplify the detection/use of waterMask water_mask_file = os.path.join(self.workDir, 'waterMask.h5') + # readfile.get_dataset_list() 会列出 HDF5 文件里有哪些数据集，例如 waterMask、latitude。 if 'waterMask' in readfile.get_dataset_list(geom_file): if ut.run_or_skip(out_file=water_mask_file, in_file=geom_file) == 'run': print(f'generate {water_mask_file} from {geom_file} for conveniency') + # readfile.read() 从 HDF5 文件读取指定数据集；返回值通常是数据数组和属性字典。 water_mask, atr = readfile.read(geom_file, datasetName='waterMask') + # 经纬度为 0 的像元通常是几何文件中的无效区域，也同步标记为水体/无效。 # ignore no-data pixels in geometry files ds_name_list = readfile.get_dataset_list(geom_file) for ds_name in ['latitude','longitude']: if ds_name in ds_name_list: print(f'set pixels with 0 in {ds_name} to 0 in waterMask') ds = readfile.read(geom_file, datasetName=ds_name)[0] + # water_mask[ds == 0] = 0 是 numpy 数组筛选写法：把经纬度为 0 的位置设为 0。 water_mask[ds == 0] = 0 atr['FILE_TYPE'] = 'waterMask' + # writefile.write() 把数组和元数据写成新的 HDF5 文件。 writefile.write(water_mask, out_file=water_mask_file, metadata=atr) # 2) modify network iargs = [stack_file, '-t', self.templateFile] print('\nmodify_network.py', ' '.join(iargs)) + # 调用 modify_network.py，根据模板配置删除低质量干涉图或修改网络连接。 import mintpy.cli.modify_network mintpy.cli.modify_network.main(iargs) # 3) plot network iargs = [stack_file, '-t', self.templateFile, '--nodisplay'] + # 相位栈和偏移栈使用不同的数据集作为网络质量指标。 dsNames = readfile.get_dataset_list(stack_file) + # any(...) 判断数据集名称中是否包含 phase/offset，从而决定绘图时用 coherence 还是 offsetSNR。 if any('phase' in i.lower() for i in dsNames): iargs += ['-d', 'coherence', '-v', '0.2', '1.0'] elif any('offset' in i.lower() for i in dsNames): @@ -254,6 +340,7 @@ def run_network_modification(self, step_name): if ut.run_or_skip(out_file=net_fig, in_file=[stack_file, coh_txt, self.templateFile], readable=False) == 'run': + # 调用 plot_network.py 生成干涉网络图。 import mintpy.cli.plot_network mintpy.cli.plot_network.main(iargs) else: @@ -262,8 +349,10 @@ def run_network_modification(self, step_name): def generate_ifgram_aux_file(self): """Generate auxiliary files from ifgramStack file""" + # 生成参考点选择和后续反演常用的干涉图辅助文件。 stack_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[0] dsNames = readfile.get_dataset_list(stack_file) + # 下面三个文件分别保存连接分量掩膜、平均空间相干性、平均空间信噪比。 mask_file = os.path.join(self.workDir, 'maskConnComp.h5') coh_file = os.path.join(self.workDir, 'avgSpatialCoh.h5') snr_file = os.path.join(self.workDir, 'avgSpatialSNR.h5') @@ -272,10 +361,12 @@ def generate_ifgram_aux_file(self): if any('phase' in i.lower() for i in dsNames): iargs = [stack_file, '--nonzero', '-o', mask_file, '--update'] print('\ngenerate_mask.py', ' '.join(iargs)) + # generate_mask.py 根据非零像元生成掩膜文件。 import mintpy.cli.generate_mask mintpy.cli.generate_mask.main(iargs) # 2) generate average spatial coherence + # temporal_average.py 虽然名字是时间平均，也可用于对干涉图栈中的 coherence/offsetSNR 求平均。 if any('phase' in i.lower() for i in dsNames): iargs = [stack_file, '--dataset', 'coherence', '-o', coh_file, '--update'] elif any('offset' in i.lower() for i in dsNames): @@ -291,10 +382,13 @@ def run_reference_point(self, step_name): 2) generate average spatial coherence and its mask 3) add REF_X/Y and/or REF_LAT/LON attribute to stack file """ + # 选择参考像元，并将参考像元信息写入干涉图栈属性。 # 1-2) aux files: maskConnComp and avgSpatialCoh + # 先确保参考点选择所需的辅助文件已经生成。 self.generate_ifgram_aux_file() # 3) add REF_X/Y(/LAT/LON) of the reference point + # lookup_file 是雷达坐标到地理坐标的查找表；没有查找表时可能为 None。 stack_file, _, lookup_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[:3] coh_file = os.path.join(self.workDir, 'avgSpatialCoh.h5') @@ -302,6 +396,7 @@ def run_reference_point(self, step_name): if lookup_file is not None: iargs += ['--lookup', lookup_file] print('\nreference_point.py', ' '.join(iargs)) + # reference_point.py 会根据模板和相干性文件选择参考点，并写入 REF_X/REF_Y 等属性。 import mintpy.cli.reference_point mintpy.cli.reference_point.main(iargs) @@ -312,6 +407,7 @@ def run_quick_overview(self, step_name): 2) numTriNonzeroIntAmbiguity.h5: phase unwrapping errors through the integer ambiguity of phase closure """ + # 快速生成两个诊断产品，用于初步检查形变信号和解缠错误。 # check the existence of ifgramStack.h5 stack_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[0] @@ -319,6 +415,7 @@ def run_quick_overview(self, step_name): pha_vel_file = os.path.join(self.workDir, 'avgPhaseVelocity.h5') iargs = [stack_file, '--dataset', 'unwrapPhase', '-o', pha_vel_file, '--update'] print('\ntemporal_average.py', ' '.join(iargs)) + # 对 unwrapPhase 做时间平均，得到一个快速查看整体形变趋势的文件。 import mintpy.cli.temporal_average mintpy.cli.temporal_average.main(iargs) @@ -326,12 +423,15 @@ def run_quick_overview(self, step_name): water_mask_file = os.path.join(self.workDir, 'waterMask.h5') iargs = [stack_file, '--water-mask', water_mask_file, '--action', 'calculate', '--update'] print('\nunwrap_error_phase_closure.py', ' '.join(iargs)) + # unwrap_error_phase_closure.py 这里用于计算相位闭合诊断量，而不是直接改正。 import mintpy.cli.unwrap_error_phase_closure mintpy.cli.unwrap_error_phase_closure.main(iargs) def run_unwrap_error_correction(self, step_name): """Correct phase-unwrapping errors""" + # 根据模板中指定的方法执行解缠误差改正；未指定方法则跳过。 + # self.template 是模板字典；用方括号读取某个配置项的值。 method = self.template['mintpy.unwrapError.method'] if not method: print('phase-unwrapping error correction is OFF.') @@ -344,17 +444,22 @@ def run_unwrap_error_correction(self, step_name): iargs_bridge = [stack_file, '--template', self.templateFile, '--update'] iargs_closure = iargs_bridge + ['--cc-mask', mask_file] + # 先导入两个可能会用到的 MintPy 命令模块，后面根据 method 选择调用哪个 main()。 import mintpy.cli.unwrap_error_bridging import mintpy.cli.unwrap_error_phase_closure if method == 'bridging': + # 基于连接分量桥接的解缠误差修正。 print('\nunwrap_error_bridging.py', ' '.join(iargs_bridge)) mintpy.cli.unwrap_error_bridging.main(iargs_bridge) elif method == 'phase_closure': + # 基于相位闭合的解缠误差修正。 print('\nunwrap_error_phase_closure.py', ' '.join(iargs_closure)) mintpy.cli.unwrap_error_phase_closure.main(iargs_closure) elif method == 'bridging+phase_closure': + # 组合方法：先桥接修正，再进行相位闭合修正。 + # -i 指定输入数据集名，-o 指定输出数据集名，避免覆盖原始 unwrapPhase。 iargs_bridge += ['-i', 'unwrapPhase', '-o', 'unwrapPhase_bridging'] print('\nunwrap_error_bridging.py', ' '.join(iargs_bridge)) @@ -366,6 +471,7 @@ def run_unwrap_error_correction(self, step_name): mintpy.cli.unwrap_error_phase_closure.main(iargs_closure) else: + # raise ValueError 表示输入值不合法，程序会停止并显示错误信息。 raise ValueError(f'un-recognized method: {method}') @@ -374,34 +480,42 @@ def run_network_inversion(self, step_name): 1) network inversion --> timeseries.h5, temporalCoherence.h5, numInvIfgram.h5 2) temporalCoherence.h5 --> maskTempCoh.h5 """ + # 将干涉网络反演为原始时序，并生成可靠像元掩膜。 # check the existence of ifgramStack.h5 stack_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[0] # 1) invert ifgramStack for time-series iargs = [stack_file, '-t', self.templateFile, '--update'] print('\nifgram_inversion.py', ' '.join(iargs)) + # ifgram_inversion.py 是核心反演步骤，把干涉图网络转换为位移时间序列。 import mintpy.cli.ifgram_inversion mintpy.cli.ifgram_inversion.main(iargs) # 2) get reliable pixel mask: maskTempCoh.h5 + # 调用本类自己的方法，根据 temporalCoherence.h5 生成可靠像元掩膜。 self.generate_temporal_coherence_mask() def generate_temporal_coherence_mask(self): """Generate reliable pixel mask from temporal coherence""" + # 根据时间相干性阈值生成可靠像元掩膜，供后续时序分析使用。 geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1] + # temporalCoherence.h5 是网络反演产生的时间相干性文件。 tcoh_file = os.path.join(self.workDir, 'temporalCoherence.h5') mask_file = os.path.join(self.workDir, 'maskTempCoh.h5') tcoh_min = self.template['mintpy.networkInversion.minTempCoh'] # compose list of arguments + # -m 是最小阈值；-o 是输出文件名。 iargs = [tcoh_file, '-m', tcoh_min, '-o', mask_file] # exclude pixels in shadow if shadowMask dataset is available if (self.template['mintpy.networkInversion.shadowMask'] is True and 'shadowMask' in readfile.get_dataset_list(geom_file)): + # 如果几何文件中有 shadowMask，就把阴影区域排除在可靠像元之外。 iargs += ['--base', geom_file, '--base-dataset', 'shadowMask', '--base-value', '1'] print('\ngenerate_mask.py', ' '.join(iargs)) + # update 模式只在输出过期或关键配置变化时重新生成掩膜。 # update mode: run only if: # 1) output file exists and newer than input file, AND # 2) all config keys are the same @@ -412,6 +526,7 @@ def generate_temporal_coherence_mask(self): flag = 'run' else: print(f'1) output file: {mask_file} already exists and newer than input file: {tcoh_file}') + # readfile.read_attribute() 只读取文件属性，不读取大数组数据，速度更快。 atr = readfile.read_attribute(mask_file) if any(str(self.template[i]) != atr.get(i, 'False') for i in config_keys): flag = 'run' @@ -424,6 +539,7 @@ def generate_temporal_coherence_mask(self): import mintpy.cli.generate_mask mintpy.cli.generate_mask.main(iargs) + # 将关键配置写入掩膜文件属性，便于下次判断是否需要重新运行。 # update config_keys atr = {} for key in config_keys: @@ -431,11 +547,13 @@ def generate_temporal_coherence_mask(self): ut.add_attribute(mask_file, atr) # check number of pixels selected in mask file for following analysis + # readfile.read(mask_file)[0] 读取掩膜数组；!= 0 找出有效像元；np.sum() 统计数量。 num_pixel = np.sum(readfile.read(mask_file)[0] != 0.) print(f'number of reliable pixels: {num_pixel}') min_num_pixel = float(self.template['mintpy.networkInversion.minNumPixel']) if num_pixel < min_num_pixel: + # 可靠像元过少时，后续分析结果不可信，因此提前停止并生成诊断图。 msg = f"Not enough reliable pixels (minimum of {int(min_num_pixel)}). " msg += "Try the following:\n" msg += "1) Check the reference pixel and make sure it's not in areas with unwrapping errors\n" @@ -443,6 +561,7 @@ def generate_temporal_coherence_mask(self): print(f'ERROR: {msg}') # populate the pic folder to facilate the trouble shooting + # 如果可靠像元太少，先自动绘图，方便用户排查参考点、网络连接等问题。 self.plot_result(print_aux=False) # terminate the program @@ -455,9 +574,11 @@ def get_timeseries_filename(template, work_dir='./'): Parameters: template : dict, content of smallbaselineApp.cfg Returns: steps : dict of dicts, input/output filenames for each step """ + # 根据模板中开启的改正项，推导每一步的输入/输出时序文件名。 work_dir = os.path.abspath(work_dir) fname0 = os.path.join(work_dir, 'timeseries.h5') fname1 = os.path.join(work_dir, 'timeseries.h5') + # 读取 timeseries.h5 的属性，用来判断卫星平台、是否已地理编码等信息。 atr = readfile.read_attribute(fname0) phase_correction_steps = [ @@ -472,21 +593,25 @@ def get_timeseries_filename(template, work_dir='./'): # loop for all steps steps = dict() for sname in phase_correction_steps: + # 默认情况下输入文件与输出文件相同，表示该步骤无需生成新文件。 # fname0 == fname1 if no valid correction method is set. fname0 = fname1 if sname == 'correct_LOD': + # Envisat 平台才需要 LOD 改正，输出文件名增加 _LOD 后缀。 if atr['PLATFORM'].lower().startswith('env'): fname1 = f'{os.path.splitext(fname0)[0]}_LOD.h5' elif sname == 'correct_SET': method = template['mintpy.solidEarthTides'] if method: + # 如果模板开启固体地球潮汐改正，输出文件名增加 _SET 后缀。 fname1 = f'{os.path.splitext(fname0)[0]}_SET.h5' elif sname == 'correct_ionosphere': method = template['mintpy.ionosphericDelay.method'] if method: if method == 'split_spectrum': + # 电离层改正输出文件名增加 _ion 后缀。 fname1 = f'{os.path.splitext(fname0)[0]}_ion.h5' else: msg = f'un-recognized ionospheric correction method: {method}' @@ -497,12 +622,15 @@ def get_timeseries_filename(template, work_dir='./'): model = template['mintpy.troposphericDelay.weatherModel'] if method: if method == 'height_correlation': + # 高程相关对流层改正输出文件名增加 _tropHgt 后缀。 fname1 = f'{os.path.splitext(fname0)[0]}_tropHgt.h5' elif method == 'gacos': + # GACOS 对流层改正输出文件名增加 _GACOS 后缀。 fname1 = f'{os.path.splitext(fname0)[0]}_GACOS.h5' elif method == 'pyaps': + # PyAPS 会把天气模型名写入输出文件名，例如 _ERA5.h5。 fname1 = f'{os.path.splitext(fname0)[0]}_{model}.h5' else: @@ -513,6 +641,7 @@ def get_timeseries_filename(template, work_dir='./'): method = template['mintpy.deramp'] if method: if method in RAMP_LIST: + # 去相位坡度后的文件名增加 _ramp 后缀。 fname1 = f'{os.path.splitext(fname0)[0]}_ramp.h5' else: msg = f'un-recognized phase ramp type: {method}' @@ -522,9 +651,11 @@ def get_timeseries_filename(template, work_dir='./'): elif sname == 'correct_topography': method = template['mintpy.topographicResidual'] if method: + # DEM 残差改正后的文件名增加 _demErr 后缀。 fname1 = f'{os.path.splitext(fname0)[0]}_demErr.h5' step = dict() + # 每个步骤记录 input 和 output，后面的 run_xxx 方法会根据它们判断是否需要运行。 step['input'] = fname0 step['output'] = fname1 steps[sname] = step @@ -532,17 +663,20 @@ def get_timeseries_filename(template, work_dir='./'): # step - reference_date fnames = [steps[sname]['output'] for sname in phase_correction_steps] fnames += [steps[sname]['input'] for sname in phase_correction_steps] + # set 去重，sorted 排序；这样 reference_date 可以处理所有可能产生过的时序文件。 fnames = sorted(list(set(fnames))) step = dict() step['input'] = fnames steps['reference_date'] = step + # velocity 和 geocode 都使用参考日期处理后的最终时序文件作为输入。 # step - velocity / geocode step = dict() step['input'] = steps['reference_date']['input'][-1] steps['velocity'] = step steps['geocode'] = step + # 若当前数据仍处于雷达坐标，HDF-EOS5 输出应使用 geo 目录中的地理编码结果。 # step - hdfeos5 if 'Y_FIRST' not in atr.keys(): step = dict() @@ -558,7 +692,9 @@ def run_local_oscillator_drift_correction(self, step_name): Automatically applied for Envisat data. Automatically skipped for all the other data. """ + # 本地振荡器漂移改正主要针对 Envisat 数据，其它卫星通常跳过。 geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1] + # get_timeseries_filename() 返回一个字典，里面告诉当前步骤的输入文件和输出文件。 fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name] in_file = fnames['input'] out_file = fnames['output'] @@ -566,9 +702,11 @@ def run_local_oscillator_drift_correction(self, step_name): iargs = [in_file, geom_file, '-o', out_file] print('\nlocal_oscilator_drift.py', ' '.join(iargs)) if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run': + # 调用 local_oscilator_drift.py 执行 Envisat LOD 改正。 import mintpy.cli.local_oscilator_drift mintpy.cli.local_oscilator_drift.main(iargs) else: + # 输入和输出文件名相同，说明这个步骤不需要生成新文件。 atr = readfile.read_attribute(in_file) sat = atr.get('PLATFORM', None) print(f'No local oscillator drift correction is needed for {sat}.') @@ -576,6 +714,7 @@ def run_local_oscillator_drift_correction(self, step_name): def run_solid_earth_tides_correction(self, step_name): """Correct solid Earth tides (SET).""" + # 根据几何文件和模板设置执行固体地球潮汐改正。 geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1] fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name] in_file = fnames['input'] @@ -585,6 +724,7 @@ def run_solid_earth_tides_correction(self, step_name): iargs = [in_file, '-g', geom_file, '-o', out_file, '--update'] print('\nsolid_earth_tides.py', ' '.join(iargs)) if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run': + # 调用 solid_earth_tides.py 给时序文件做固体地球潮汐改正。 import mintpy.cli.solid_earth_tides mintpy.cli.solid_earth_tides.main(iargs) else: @@ -593,6 +733,7 @@ def run_solid_earth_tides_correction(self, step_name): def run_ionospheric_delay_correction(self, step_name): """Correct ionospheric delays.""" + # 电离层改正依赖加载阶段生成的 ionosphere stack。 iono_stack_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[3] fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name] @@ -601,12 +742,14 @@ def run_ionospheric_delay_correction(self, step_name): if in_file != out_file: method = self.template['mintpy.ionosphericDelay.method'] + # 当前工作流支持 split_spectrum 方法。 # range split spectrum (Fattahi et al., 2017; Liang et al. 2018; 2019) if method == 'split_spectrum': print(f'ionospheric delay correction with {method} approach') iargs = ['-t', self.templateFile, '-f', in_file, '-o', out_file, '--iono-stack-file', iono_stack_file] print('\niono_split_spectrum.py', ' '.join(iargs)) + # 调用 iono_split_spectrum.py 使用分裂谱方法改正电离层延迟。 import mintpy.cli.iono_split_spectrum mintpy.cli.iono_split_spectrum.main(iargs) else: @@ -615,19 +758,23 @@ def run_ionospheric_delay_correction(self, step_name): def run_tropospheric_delay_correction(self, step_name): """Correct tropospheric delays.""" + # 对流层延迟改正支持高程相关、GACOS 和 PyAPS 天气再分析三类方法。 geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1] + # maskTempCoh.h5 是可靠像元掩膜，用来限制对流层拟合只在可信像元上进行。 mask_file = os.path.join(self.workDir, 'maskTempCoh.h5') fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name] in_file = fnames['input'] out_file = fnames['output'] if in_file != out_file: + # 从模板读取对流层改正所需参数。 poly_order = self.template['mintpy.troposphericDelay.polyOrder'] tropo_model = self.template['mintpy.troposphericDelay.weatherModel'].upper() weather_dir = self.template['mintpy.troposphericDelay.weatherDir'] method = self.template['mintpy.troposphericDelay.method'] def get_dataset_size(fname): + # 用文件属性中的 LENGTH/WIDTH 判断两个 HDF5 数据集尺寸是否一致。 atr = readfile.read_attribute(fname) return (atr['LENGTH'], atr['WIDTH']) @@ -635,6 +782,7 @@ def get_dataset_size(fname): if method == 'height_correlation': tropo_look = self.template['mintpy.troposphericDelay.looks'] tropo_min_cor = self.template['mintpy.troposphericDelay.minCorrelation'] + # -g 几何文件，-p 多项式阶数，-m 掩膜文件，-l looks，-t 最小相关系数阈值。 iargs = [in_file, '-g', geom_file, '-p', poly_order, @@ -645,44 +793,55 @@ def get_dataset_size(fname): print('tropospheric delay correction with height-correlation approach') print('\ntropo_phase_elevation.py', ' '.join(iargs)) if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run': + # 调用 tropo_phase_elevation.py 用相位-高程相关关系估计对流层延迟。 import mintpy.cli.tropo_phase_elevation mintpy.cli.tropo_phase_elevation.main(iargs) # Weather re-analysis data with iterative tropospheric decomposition (GACOS) # Yu et al. (2018, JGR) elif method == 'gacos': + # GACOS_dir 是 GACOS 数据所在目录，由模板指定。 GACOS_dir = self.template['mintpy.troposphericDelay.gacosDir'] iargs = ['-f', in_file, '-g', geom_file, '-o', out_file, '--dir', GACOS_dir] print('tropospheric delay correction with gacos approach') print('\ntropo_gacos.py', ' '.join(iargs)) if ut.run_or_skip(out_file=out_file, in_file=in_file) == 'run': + # 调用 tropo_gacos.py 使用 GACOS 产品改正对流层延迟。 import mintpy.cli.tropo_gacos mintpy.cli.tropo_gacos.main(iargs) # Weather Re-analysis Data (Jolivet et al., 2011;2014) elif method == 'pyaps': + # PyAPS 需要天气模型名称、几何文件以及天气数据保存目录。 iargs = ['-f', in_file, '--model', tropo_model, '-g', geom_file, '-w', weather_dir] print('Atmospheric correction using Weather Re-analysis dataset (PyAPS, Jolivet et al., 2011)') print('Weather Re-analysis dataset:', tropo_model) + # 预估/缓存的天气模型延迟文件通常放在 inputs 目录下。 tropo_file = f'./inputs/{tropo_model}.h5' if ut.run_or_skip(out_file=out_file, in_file=[in_file, tropo_file]) == 'run': if os.path.isfile(tropo_file) and get_dataset_size(tropo_file) == get_dataset_size(in_file): + # 若已有尺寸匹配的对流层延迟文件，则直接与时序文件相减。 iargs = [in_file, tropo_file, '-o', out_file, '--force'] print('--------------------------------------------') print(f'Use existed tropospheric delay file: {tropo_file}') print('\ndiff.py', ' '.join(iargs)) + # diff.py 这里用于从输入时序中减去已有的对流层延迟文件。 import mintpy.cli.diff mintpy.cli.diff.main(iargs) else: if tropo_model in ['ERA5']: + # ERA5 使用 PyAPS3 接口。 print('\ntropo_pyaps3.py', ' '.join(iargs)) + # tropo_pyaps3.py 会下载/读取 ERA5 并生成对流层改正后的时序。 import mintpy.cli.tropo_pyaps3 mintpy.cli.tropo_pyaps3.main(iargs) elif tropo_model in ['MERRA', 'NARR']: + # MERRA/NARR 保留旧版 PyAPS 处理入口。 print('\ntropo_pyaps.py', ' '.join(iargs)) + # legacy 表示旧版模块；这里为了兼容 MERRA/NARR 仍然调用它。 import mintpy.legacy.tropo_pyaps mintpy.legacy.tropo_pyaps.main(iargs) @@ -695,6 +854,8 @@ def get_dataset_size(fname): def run_phase_deramping(self, step_name): """Estimate and remove phase ramp from each acquisition.""" + # 按模板指定的 ramp 类型估计并移除每期影像中的相位坡度。 + # mask_file 指定哪些像元参与坡度估计；method 指定坡度模型类型。 mask_file = self.template['mintpy.deramp.maskFile'] method = self.template['mintpy.deramp'] @@ -705,6 +866,7 @@ def run_phase_deramping(self, step_name): print(f'Remove for each acquisition a phase ramp: {method}') iargs = [in_file, '-s', method, '-m', mask_file, '-o', out_file, '--update'] print('\nremove_ramp.py', ' '.join(iargs)) + # 调用 remove_ramp.py 去除每期影像中的平面/二次等相位坡度。 import mintpy.cli.remove_ramp mintpy.cli.remove_ramp.main(iargs) else: @@ -715,6 +877,7 @@ def run_topographic_residual_correction(self, step_name): """step - correct_topography Topographic residual (DEM error) correction (optional). """ + # 估计并改正 DEM 误差造成的残余地形相位。 geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1] fnames = self.get_timeseries_filename(self.template, self.workDir)[step_name] in_file = fnames['input'] @@ -723,8 +886,10 @@ def run_topographic_residual_correction(self, step_name): if in_file != out_file: iargs = [in_file, '-t', self.templateFile, '-o', out_file, '--update'] if self.template['mintpy.topographicResidual.pixelwiseGeometry']: + # 如果模板要求逐像元几何，就额外把几何文件通过 -g 传给 dem_error.py。 iargs += ['-g', geom_file] print('\ndem_error.py', ' '.join(iargs)) + # dem_error.py 估计 DEM 残差并从时序中扣除相关相位。 import mintpy.cli.dem_error mintpy.cli.dem_error.main(iargs) @@ -734,10 +899,12 @@ def run_topographic_residual_correction(self, step_name): def run_residual_phase_rms(self, step_name): """Noise evaluation based on the phase residual.""" + # 若残差时序文件存在，则计算 RMS 作为噪声诊断指标。 res_file = 'timeseriesResidual.h5' if os.path.isfile(res_file): iargs = [res_file, '-t', self.templateFile] print('\ntimeseries_rms.py', ' '.join(iargs)) + # timeseries_rms.py 根据残差时序计算 RMS，用于评估噪声水平。 import mintpy.cli.timeseries_rms mintpy.cli.timeseries_rms.main(iargs) else: @@ -746,12 +913,15 @@ def run_residual_phase_rms(self, step_name): def run_reference_date(self, step_name): """Change reference date for all time-series files (optional).""" + # 将所有中间/最终时序文件统一调整到模板指定的参考日期。 if self.template['mintpy.reference.date']: iargs = ['-t', self.templateFile] in_files = self.get_timeseries_filename(self.template, self.workDir)[step_name]['input'] for in_file in in_files: + # reference_date.py 支持一次传入多个时序文件，这里逐个追加到参数列表。 iargs += [in_file] print('\nreference_date.py', ' '.join(iargs)) + # 调用 reference_date.py 修改时序文件的参考日期。 import mintpy.cli.reference_date mintpy.cli.reference_date.main(iargs) else: @@ -760,11 +930,14 @@ def run_reference_date(self, step_name): def run_timeseries2velocity(self, step_name): """Estimate average velocity from displacement time-series""" + # 从最终位移时序估计平均速度。 ts_file = self.get_timeseries_filename(self.template, self.workDir)[step_name]['input'] + # velocity.h5 是 MintPy 默认的平均速度输出文件。 vel_file = os.path.join(self.workDir, 'velocity.h5') iargs = [ts_file, '-t', self.templateFile, '-o', vel_file, '--update'] print('\ntimeseries2velocity.py', ' '.join(iargs)) + # timeseries2velocity.py 会对位移时序拟合速度模型，生成 velocity.h5。 import mintpy.cli.timeseries2velocity mintpy.cli.timeseries2velocity.main(iargs) @@ -772,12 +945,14 @@ def run_timeseries2velocity(self, step_name): tropo_model = self.template['mintpy.troposphericDelay.weatherModel'].upper() tropo_file = os.path.join(self.workDir, f'inputs/{tropo_model}.h5') if os.path.isfile(tropo_file): + # 如果存在估计的对流层延迟文件，也计算其速度趋势，便于评估大气改正影响。 suffix = os.path.splitext(os.path.basename(tropo_file))[0] tropo_vel_file = f'{os.path.splitext(vel_file)[0]}{suffix}.h5' tropo_vel_file = os.path.join(self.workDir, tropo_vel_file) iargs = [tropo_file, '-t', self.templateFile, '-o', tropo_vel_file, '--update'] # add reference info for a meaningful velocity to assess the impact of tropo delay on velocity + # 为了让对流层延迟速度和主速度文件可比较，需要使用相同参考日期和参考像元。 atr = readfile.read_attribute(vel_file) iargs += ['--ref-date', atr['REF_DATE'], '--ref-yx', atr['REF_Y'], atr['REF_X']] print('\ntimeseries2velocity.py', ' '.join(iargs)) @@ -786,16 +961,20 @@ def run_timeseries2velocity(self, step_name): def run_geocode(self, step_name): """geocode data files in radar coordinates into ./geo folder.""" + # 将雷达坐标产品地理编码到 geo 目录；已是地理坐标的数据会跳过。 if self.template['mintpy.geocode']: ts_file = self.get_timeseries_filename(self.template, self.workDir)[step_name]['input'] atr = readfile.read_attribute(ts_file) + # Y_FIRST 是 MintPy 常用的地理坐标属性；没有它通常表示数据还在雷达坐标。 if 'Y_FIRST' not in atr.keys(): # 1. geocode out_dir = os.path.join(self.workDir, 'geo') os.makedirs(out_dir, exist_ok=True) + # 对关键结果文件统一地理编码，保持输出产品坐标一致。 geom_file, lookup_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1:3] in_files = [geom_file, 'temporalCoherence.h5', 'avgSpatialCoh.h5', ts_file, 'velocity.h5'] + # geocode.py 用 lookup_file 把雷达坐标文件转换到地理坐标，并输出到 geo 目录。 iargs = in_files + ['-l', lookup_file, '-t', self.templateFile, '--outdir', out_dir, '--update'] print('\ngeocode.py', ' '.join(iargs)) import mintpy.cli.geocode @@ -811,10 +990,12 @@ def run_geocode(self, step_name): # exclude pixels in shadow if shadowMask dataset is available if (self.template['mintpy.networkInversion.shadowMask'] is True and 'shadowMask' in readfile.get_dataset_list(geom_file)): + # 在地理坐标下重新生成可靠像元掩膜，同时排除阴影区域。 iargs += ['--base', geom_file, '--base-dataset', 'shadowMask', '--base-value', '1'] print('\ngenerate_mask.py', ' '.join(iargs)) if ut.run_or_skip(out_file=mask_file, in_file=tcoh_file) == 'run': + # 调用 generate_mask.py 生成 geo_maskTempCoh.h5。 import mintpy.cli.generate_mask mintpy.cli.generate_mask.main(iargs) @@ -826,26 +1007,31 @@ def run_geocode(self, step_name): def run_save2google_earth(self, step_name): """Save velocity file in geo coordinates into Google Earth raster image.""" + # 按模板开关将速度结果保存为 Google Earth 可读的 KMZ 文件。 if self.template['mintpy.save.kmz'] is True: print('creating Google Earth KMZ file for geocoded velocity file: ...') # input vel_file = os.path.join(self.workDir, 'velocity.h5') atr = readfile.read_attribute(vel_file) if 'Y_FIRST' not in atr.keys(): + # 雷达坐标速度文件需要改用地理编码后的版本。 vel_file = os.path.join(self.workDir, 'geo/geo_velocity.h5') # output + # os.path.splitext(vel_file)[0] 去掉 .h5 后缀，再加 .kmz 得到输出文件名。 kmz_file = f'{os.path.splitext(vel_file)[0]}.kmz' iargs = [vel_file, '-o', kmz_file] print('\nsave_kmz.py', ' '.join(iargs)) # update mode fbase = os.path.basename(kmz_file) + # KMZ 文件可能已经在当前目录、geo 目录或 pic 目录中，这里都检查一遍。 kmz_files = [i for i in [fbase, f'./geo/{fbase}', f'./pic/{fbase}'] if os.path.isfile(i)] kmz_file = kmz_files[0] if len(kmz_files) > 0 else None if ut.run_or_skip(out_file=kmz_file, in_file=vel_file, readable=False) == 'run': + # save_kmz.py 把地理坐标速度文件转换为 Google Earth 可打开的 KMZ 栅格图。 import mintpy.cli.save_kmz mintpy.cli.save_kmz.main(iargs) @@ -855,11 +1041,13 @@ def run_save2google_earth(self, step_name): def run_save2hdfeos5(self, step_name): """Save displacement time-series and its aux data in geo coordinate into HDF-EOS5 format""" + # 按模板开关将时序、相干性、掩膜和几何信息打包为 HDF-EOS5 产品。 if self.template['mintpy.save.hdfEos5'] is True: # input ts_file = self.get_timeseries_filename(self.template, self.workDir)[step_name]['input'] # Add attributes from custom template to timeseries file if self.customTemplate is not None: + # 输出前补充用户模板中的属性，方便产品元数据完整记录。 ut.add_attribute(ts_file, self.customTemplate) tcoh_file = os.path.join(self.workDir, 'temporalCoherence.h5') @@ -867,6 +1055,7 @@ def run_save2hdfeos5(self, step_name): mask_file = os.path.join(self.workDir, 'maskTempCoh.h5') geom_file = ut.check_loaded_dataset(self.workDir, print_msg=False)[1] if 'geo' in ts_file: + # 如果输入时序已经位于 geo 目录，同步使用 geo 目录中的辅助文件。 tcoh_file = os.path.join(self.workDir, 'geo/geo_temporalCoherence.h5') scoh_file = os.path.join(self.workDir, 'geo/geo_avgSpatialCoh.h5') mask_file = os.path.join(self.workDir, 'geo/geo_maskTempCoh.h5') @@ -874,6 +1063,7 @@ def run_save2hdfeos5(self, step_name): # cmd print('--------------------------------------------') + # save_hdfeos5.py 需要同时知道时序、时间相干性、空间相干性、掩膜、几何和模板文件。 iargs = [ts_file, '--tc', tcoh_file, '--asc', scoh_file, @@ -884,7 +1074,9 @@ def run_save2hdfeos5(self, step_name): # output (check existing file) atr = readfile.read_attribute(ts_file) + # sensor.get_unavco_mission_name() 根据文件属性推断 UNAVCO/HDF-EOS5 产品所需的卫星任务名。 SAT = sensor.get_unavco_mission_name(atr) + # ut.get_file_list('SAT_*.he5') 查找是否已经有对应卫星名前缀的 HDF-EOS5 输出文件。 hdfeos5_files = ut.get_file_list(f'{SAT}_*.he5') hdfeos5_file = hdfeos5_files[0] if len(hdfeos5_files) > 0 else None @@ -892,6 +1084,7 @@ def run_save2hdfeos5(self, step_name): in_file=[ts_file, tcoh_file, scoh_file, mask_file, geom_file]) == 'run': + # save_hdfeos5.py 生成标准 HDF-EOS5 格式产品。 import mintpy.cli.save_hdfeos5 mintpy.cli.save_hdfeos5.main(iargs) @@ -901,9 +1094,12 @@ def run_save2hdfeos5(self, step_name): def run(self, steps): """run the chosen steps.""" + # 按解析出的步骤列表顺序执行，每个步骤映射到对应的成员函数。 for sname in steps: print(f'\n\n******************** step - {sname} ********************') + # 这一长串 if/elif 是“调度表”：根据步骤名 sname，调用对应的类方法。 + # 例如 sname == 'load_data' 时，就调用 self.run_load_data(sname)。 if sname == 'load_data': self.run_load_data(sname) @@ -962,13 +1158,16 @@ def run(self, steps): def plot_result(self, print_aux=True): """Plot data files and save to figures in pic folder""" + # 将主要中间结果和最终结果批量绘图，并统一移动到 pic 目录。 print('\n******************** plot & save to pic ********************') + # 从模板读取绘图和计算资源相关参数。 tropo_model = self.template['mintpy.troposphericDelay.weatherModel'].upper() max_plot_memory = abs(float(self.template['mintpy.plot.maxMemory'])) max_memory = abs(float(self.template['mintpy.compute.maxMemory'])) fig_dpi = int(self.template['mintpy.plot.dpi']) + # 再次检查核心输入文件位置，绘图时会用到它们。 stack_file, geom_file, lookup_file, ion_file = ut.check_loaded_dataset( self.workDir, print_msg=False)[:4] @@ -976,6 +1175,7 @@ def plot_result(self, print_aux=True): geo_dir = os.path.join(self.workDir, 'geo') pic_dir = os.path.join(self.workDir, 'pic') + # 打印 view 命令时使用相对路径，输出更短也更容易阅读。 # use relative path for shorter and cleaner printout view command stack_file = os.path.relpath(stack_file) if stack_file else stack_file ion_file = os.path.relpath(ion_file) if ion_file else ion_file @@ -988,7 +1188,10 @@ def plot_result(self, print_aux=True): # for each element list: # the 1st item is the data file # the 2nd item is the dataset if applicable + # opt4ts 是绘制时序文件时常用的一组 view.py 参数。 opt4ts = ['--noaxis', '-u', 'cm', '--wrap', '--wrap-range', '-5', '5'] + # iargs_list0 中每个小列表都代表一次 view.py 调用的参数。 + # 例如 ['velocity.h5', '--dem', geom_file, '--mask', mask_file] 表示绘制 velocity.h5。 iargs_list0 = [ # key files ['velocity.h5', '--dem', geom_file, '--mask', mask_file], @@ -1032,6 +1235,7 @@ def plot_result(self, print_aux=True): ] if ion_file: + # 如果存在电离层栈，也把相关数据集加入绘图列表。 iargs_list0 += [ [ion_file, 'unwrapPhase-', '--zero-mask', '--wrap', '-c', 'cmy'], [ion_file, 'unwrapPhase-', '--zero-mask'], @@ -1043,9 +1247,12 @@ def plot_result(self, print_aux=True): for iargs in iargs_list0: fname, args = iargs[0], iargs[1:] if not fname: + # 如果文件名是 None 或空字符串，说明这个输入不存在，直接跳过。 continue if '*' in fname: + # 展开通配符路径，逐个生成 view.py 参数列表。 + # glob.glob() 会把 'timeseries_*.h5' 这种通配符展开成真实文件名列表。 fnames = sorted(glob.glob(fname)) if len(fnames) > 0: for fname in fnames: @@ -1055,9 +1262,11 @@ def plot_result(self, print_aux=True): iargs_list.append(iargs) # remove element list - file does not exists + # 只保留真实存在的文件，避免 view.py 因找不到文件而报错。 iargs_list = [iargs for iargs in iargs_list if os.path.isfile(iargs[0])] # remote element list - file is ifgramStack and the dataset of interest does not exist + # 对 ifgramStack.h5 这种一个文件包含多个数据集的情况，检查要绘制的数据集是否真的存在。 stack_dset_list = readfile.get_dataset_list(stack_file) stack_dset_list += [f'{x}-' for x in stack_dset_list] iargs_list = [iargs for iargs in iargs_list @@ -1066,8 +1275,10 @@ def plot_result(self, print_aux=True): and iargs[1] in stack_dset_list))] # add the following common options to all element lists + # opt_common 是所有 view.py 调用共用的绘图选项，例如 dpi、无界面显示、内存限制等。 opt_common = ['--dpi', str(fig_dpi), '--noverbose', '--nodisplay', '--update', '--memory', str(max_plot_memory)] + # 把通用参数加到每一次 view.py 调用前面。 iargs_list = [opt_common + iargs for iargs in iargs_list] # run view @@ -1075,7 +1286,9 @@ def plot_result(self, print_aux=True): run_parallel = False cluster_type = self.template['mintpy.compute.cluster'] if cluster_type: + # 若模板指定并行计算，则根据 CPU 数和内存限制确定实际并行度。 num_workers = self.template['mintpy.compute.numWorker'] + # DaskCluster.format_num_worker() 会把模板里的 worker 数量配置转换成可用的数字。 num_workers = cluster.DaskCluster.format_num_worker(cluster_type, num_workers) # limit number of parallel processes based on available CPU @@ -1093,15 +1306,19 @@ def plot_result(self, print_aux=True): import mintpy.cli.view if run_parallel and num_cores > 1: print(f"parallel processing using {num_cores} cores ...") + # Parallel/delayed 来自 joblib；这里把多次 view.py 绘图任务并行运行。 Parallel(n_jobs=num_cores)(delayed(mintpy.cli.view.main)(iargs) for iargs in iargs_list) else: for iargs in iargs_list: + # 串行模式：一个一个调用 view.py 绘图。 mintpy.cli.view.main(iargs) # copy text files to pic print('copy *.txt files into ./pic directory.') + # glob.glob('*.txt') 查找当前目录下所有 txt 文件。 tfiles = glob.glob('*.txt') for tfile in tfiles: + # 把文本诊断文件复制到 pic 目录，方便集中查看。 shutil.copy2(tfile, pic_dir) # move picture files to pic @@ -1111,9 +1328,12 @@ def plot_result(self, print_aux=True): pfiles += glob.glob('*.kmz') pfiles += glob.glob(os.path.join(geo_dir, '*.kmz')) for pfile in pfiles: + # shutil.move() 移动文件到 pic 目录；os.path.basename() 保留原文件名。 shutil.move(pfile, os.path.join(pic_dir, os.path.basename(pfile))) # time info + # time.time() 返回当前时间戳；与 start_time 相减得到绘图耗时秒数。 + # divmod(seconds, 60) 把秒数拆成分钟和剩余秒数。 m, s = divmod(time.time()-start_time, 60) print(f'time used: {m:02.0f} mins {s:02.1f} secs.') @@ -1134,8 +1354,10 @@ def plot_result(self, print_aux=True): def close(self, normal_end=True): + # 流程结束后回到启动程序时所在的目录。 # go back to original directory print('Go back to directory:', self.cwd) + # os.chdir(self.cwd) 切回最开始的目录，避免程序结束后还停留在工作目录中。 os.chdir(self.cwd) # message if normal_end: @@ -1147,23 +1369,30 @@ def close(self, normal_end=True): def run_smallbaselineApp(inps): """Run the small baseline time series analsysis workflow.""" + # smallbaselineApp 的程序化入口：打开工作流、执行步骤、按需绘图并收尾。 start_time = time.time() # open and run + # 创建 TimeSeriesAnalysis 对象；这会自动调用 __init__ 保存模板路径和工作目录。 app = TimeSeriesAnalysis(inps.customTemplateFile, inps.workDir) + # app.open() 负责进入工作目录、准备模板、读取配置。 app.open() + # app.run(...) 根据 CLI 文件解析出来的 runSteps 逐步执行处理流程。 app.run(steps=inps.runSteps) # plot if: # a) --plot in command line OR # b) template['mintpy.plot'] = yes AND runSteps > 1 + # 如果用户明确要求绘图，或模板开启绘图且本次运行了多个步骤，就调用 plot_result()。 if inps.plot or (app.template['mintpy.plot'] and len(inps.runSteps) > 1): app.plot_result() # close + # close() 会切回原目录，并打印正常结束信息。 app.close() # used time + # 统计整个 smallbaselineApp 工作流从开始到结束的总耗时。 m, s = divmod(time.time()-start_time, 60) print(f'Time used: {m:02.0f} mins {s:02.1f} secs\n') diff --git a/src/mintpy/timeseries2velocity.py b/src/mintpy/timeseries2velocity.py index b2148b1fc..7de102d0a 100644 --- a/src/mintpy/timeseries2velocity.py +++ b/src/mintpy/timeseries2velocity.py @@ -7,17 +7,25 @@ # from mintpy import timeseries2velocity as ts2vel +# os 用来处理路径、判断文件是否存在、比较输入输出文件修改时间。 import os +# time 用来统计速度拟合运行耗时。 import time +# numpy 用于数组计算；时序数据和拟合参数都以 numpy 数组表示。 import numpy as np +# scipy.linalg 用于协方差传播等线性代数计算。 from scipy import linalg +# 这些对象用于读取不同类型的时序文件：普通 timeseries、GIAnT 时序、HDF-EOS5 产品。 from mintpy.objects import HDFEOS, cluster, giantTimeseries, timeseries +# ptime 处理日期；time_func 构建/估计时间函数模型；readfile/writefile 读写 HDF5。 from mintpy.utils import ptime, readfile, time_func, writefile +# 输出拟合参数统一使用 float32，兼顾精度和文件大小。 DATA_TYPE = np.float32 # key configuration parameter name +# 这些 timeFunc 配置会写入 velocity.h5 属性，用于 update 模式判断是否需要重跑。 key_prefix = 'mintpy.timeFunc.' config_keys = [ # date @@ -39,6 +47,7 @@ ############################################################################ def run_or_skip(inps): + # update 模式判断：如果输出 velocity 文件存在、比输入新、配置没变，就跳过。 print('update mode: ON') flag = 'skip' @@ -58,6 +67,7 @@ def run_or_skip(inps): # check configuration if flag == 'skip': + # 读取输出文件属性，比对本次拟合模型配置是否与上次一致。 atr = readfile.read_attribute(inps.outfile) if any(str(vars(inps)[key]) != atr.get(key_prefix+key, 'None') for key in config_keys): flag = 'run' @@ -80,6 +90,7 @@ def read_date_info(inps): dropDate - 1D np.ndarray in bool in size of all available dates """ # initiate and open time-series file object + # 根据 FILE_TYPE 选择正确的时序读取对象。 ftype = readfile.read_attribute(inps.timeseries_file)['FILE_TYPE'] if ftype == 'timeseries': ts_obj = timeseries(inps.timeseries_file) @@ -92,6 +103,7 @@ def read_date_info(inps): ts_obj.open() # exclude dates - user inputs + # 根据 startDate/endDate/excludeDate 计算不参与拟合的日期列表。 ex_date_list = ptime.get_exclude_date_list( date_list=ts_obj.dateList, start_date=inps.startDate, @@ -100,6 +112,7 @@ def read_date_info(inps): # exclude dates - no obs data [for offset time-series only for now] if os.path.basename(inps.timeseries_file).startswith('timeseriesRg'): + # 对偏移时序，如果某些日期全为 0，说明没有有效观测，应排除。 data, atr = readfile.read(inps.timeseries_file) flag = np.nansum(data, axis=(1,2)) == 0 flag[ts_obj.dateList.index(atr['REF_DATE'])] = 0 @@ -109,9 +122,11 @@ def read_date_info(inps): ex_date_list = sorted(list(set(ex_date_list))) # dates used for estimation - inps.date_list + # date_list 是真正参与时间函数拟合的日期。 inps.date_list = [i for i in ts_obj.dateList if i not in ex_date_list] # flag array for ts data reading + # dropDate 是布尔数组：True 表示该日期保留，False 表示该日期丢弃。 inps.dropDate = np.array([i not in ex_date_list for i in ts_obj.dateList], dtype=np.bool_) # print out msg @@ -128,9 +143,11 @@ def read_date_info(inps): def run_timeseries2time_func(inps): + # 主入口：读取位移时序，拟合时间函数参数，并写入 velocity.h5。 start_time = time.time() # basic file info + # LENGTH/WIDTH 是图像行列数，用来创建输出数据集。 atr = readfile.read_attribute(inps.timeseries_file) length, width = int(atr['LENGTH']), int(atr['WIDTH']) @@ -147,13 +164,16 @@ def run_timeseries2time_func(inps): print(f' Set "--ref-date {inps.date_list[0]}" and continue.') # get deformation model from inputs + # model 描述要拟合哪些时间函数：速度、多项式、周期、阶跃、指数/对数等。 model = time_func.inps2model(inps, date_list=inps.date_list) + # num_param 是模型参数数量，例如线性速度模型通常有 intercept 和 velocity 两个参数。 num_param = time_func.get_num_param(model) ## output preparation # time_func_param: attributes + # atrV 是输出 velocity.h5 的元数据，继承输入时序并改写文件类型/单位/日期范围。 date0, date1 = inps.date_list[0], inps.date_list[-1] atrV = dict(atr) atrV['FILE_TYPE'] = 'velocity' @@ -168,11 +188,13 @@ def run_timeseries2time_func(inps): atrV['REF_DATE'] = inps.ref_date # time_func_param: config parameter + # 把本次拟合配置写进 HDF5 属性，方便复现和 update 模式检查。 print(f'add/update the following configuration metadata:\n{config_keys}') for key in config_keys: atrV[key_prefix+key] = str(vars(inps)[key]) # time_func_param: instantiate output file + # model2hdf5_dataset() 根据 model 自动生成 velocity、annualAmplitude、stepYYYYMMDD 等输出数据集定义。 ds_name_dict, ds_unit_dict = model2hdf5_dataset(model, ds_shape=(length, width))[1:] # add dataset: residue if inps.uncertaintyQuantification == 'residue': @@ -186,6 +208,7 @@ def run_timeseries2time_func(inps): # timeseries_res: attributes + instantiate output file if inps.save_res: + # 可选输出残差时序：原时序减去拟合时间函数后的剩余部分。 atrR = dict(atr) # remove REF_DATE attribute for key in ['REF_DATE']: @@ -203,6 +226,7 @@ def run_timeseries2time_func(inps): ## estimation # calc number of box based on memory limit + # 按内存上限估算需要把图像切成多少个空间块。 memoryAll = (num_date + num_param * 2 + 2) * length * width * 4 if inps.uncertaintyQuantification == 'bootstrap': memoryAll += inps.bootstrapCount * num_param * length * width * 4 @@ -216,6 +240,7 @@ def run_timeseries2time_func(inps): # loop for block-by-block IO for i, box in enumerate(box_list): + # 当前 box 是空间块范围 (x0, y0, x1, y1)。 box_wid = box[2] - box[0] box_len = box[3] - box[1] num_pixel = box_len * box_wid @@ -225,21 +250,25 @@ def run_timeseries2time_func(inps): print(f'box length: {box_len}') # initiate output + # m 保存模型参数，m_std 保存参数标准差；形状是 num_param x num_pixel。 m = np.zeros((num_param, num_pixel), dtype=DATA_TYPE) m_std = np.zeros((num_param, num_pixel), dtype=DATA_TYPE) # read input + # 读取当前空间块的时序数据。 print(f'reading data from file {inps.timeseries_file} ...') ts_data = readfile.read(inps.timeseries_file, box=box)[0] # referencing in time and space # for file w/o reference info. e.g. ERA5.h5 if inps.ref_date: + # 如果指定参考日期，就把每个像元所有日期减去该参考日期的值。 print(f'referecing to date: {inps.ref_date}') ref_ind = inps.date_list.index(inps.ref_date) ts_data -= np.tile(ts_data[ref_ind, :, :], (ts_data.shape[0], 1, 1)) if inps.ref_yx: + # 如果指定参考点，就把每个日期所有像元减去该参考点的值。 print(f'referencing to point (y, x): ({inps.ref_yx[0]}, {inps.ref_yx[1]})') ref_box = (inps.ref_yx[1], inps.ref_yx[0], inps.ref_yx[1]+1, inps.ref_yx[0]+1) ref_val = readfile.read(inps.timeseries_file, box=ref_box)[0] @@ -248,10 +277,12 @@ def run_timeseries2time_func(inps): ts_data = ts_data[inps.dropDate, :, :].reshape(num_date, -1) if atrV['UNIT'] == 'mm': + # 输出单位需要是米/年；如果输入是毫米，先转成米。 ts_data *= 1./1000. ts_cov = None if inps.uncertaintyQuantification == 'covariance': + # covariance 模式从时间序列协方差文件传播得到模型参数不确定性。 print(f'reading time-series covariance matrix from file {inps.timeSeriesCovFile} ...') ts_cov = readfile.read(inps.timeSeriesCovFile, box=box)[0] if len(ts_cov.shape) == 4: @@ -272,6 +303,7 @@ def run_timeseries2time_func(inps): #ts_cov[ts_cov 2e3: inps.disp_dem_contour = False @@ -90,6 +106,7 @@ def update_inps_with_file_metadata(inps, metadata): inps.multilook = True # Colormap + # 根据文件类型和数据集名称自动选择合适的色带。 inps.colormap = pp.auto_colormap_name( metadata, inps.colormap, @@ -101,6 +118,7 @@ def update_inps_with_file_metadata(inps, metadata): # Convert ref_lalo if existed, to ref_yx, and use ref_yx for the following # ref_yx is referenced to input data coverage, not subseted area for display if inps.ref_lalo: + # 如果用户给的是经纬度参考点，转换成文件中的 y/x 像素坐标。 vprint(f'input reference point in lat/lon: {inps.ref_lalo}') if not inps.geo_box and not coord.lookup_file: print('WARNING: --ref-lalo is NOT supported when 1) file is radar-coded AND 2) no lookup table file found') @@ -120,9 +138,11 @@ def update_inps_with_file_metadata(inps, metadata): inps.ref_lalo = None # Points of interest + # 读取用户指定的点文件或点坐标，用于在图上标记。 inps = pp.read_pts2inps(inps, coord) # Unit and Wrap + # 根据数据单位和 --wrap 设置，决定显示单位、缩放倍数以及是否相位缠绕显示。 inps.disp_unit, inps.wrap = pp.check_disp_unit_and_wrap( metadata, disp_unit=inps.disp_unit, @@ -132,6 +152,7 @@ def update_inps_with_file_metadata(inps, metadata): ) # Map Projection via cartopy + # 根据文件是否地理编码、是否要求 coast/lalo-label，决定是否启用 cartopy 投影。 inps = check_map_projection(inps, metadata, print_msg=inps.print_msg) # Min / Max - Display @@ -144,6 +165,7 @@ def update_inps_with_file_metadata(inps, metadata): inps.vlim = inps.wrap_range # Transparency - Alpha + # alpha 是透明度；叠加 DEM 阴影时默认让主数据稍透明。 if not inps.transparency: # Auto adjust transparency value when showing shaded relief DEM if inps.dem_file and inps.disp_dem_shade: @@ -166,6 +188,7 @@ def update_inps_with_file_metadata(inps, metadata): # Figure output file name if not inps.outfile: + # 没指定输出文件名时，根据图标题自动生成图片名。 # ignore whitespaces in the filename fbase = inps.fig_title.replace(' ', '') if not inps.disp_whitespace: @@ -195,6 +218,7 @@ def check_map_projection(inps, metadata, print_msg=True): """ inps.map_proj_obj = None + # Y_UNIT 是坐标单位，常见为 degrees 或 meter。 inps.coord_unit = metadata.get('Y_UNIT', 'degrees').lower() if (inps.geo_box and inps.coord_unit.startswith(('deg', 'meter')) @@ -205,11 +229,13 @@ def check_map_projection(inps, metadata, print_msg=True): # https://scitools.org.uk/cartopy/docs/latest/crs/projections.html msg = 'initiate cartopy map projection: ' if inps.coord_unit.startswith('deg'): + # PlateCarree 是经纬度坐标最常见的地图投影。 inps.map_proj_obj = ccrs.PlateCarree() print(msg + 'PlateCarree') elif inps.coord_unit.startswith('meter'): if 'UTM_ZONE' in metadata.keys(): + # UTM 是以米为单位的投影坐标；需要知道 UTM 分区。 utm_zone = metadata['UTM_ZONE'] inps.map_proj_obj = ccrs.UTM(utm_zone) print(msg + f'UTM zone {utm_zone}') @@ -228,6 +254,7 @@ def check_map_projection(inps, metadata, print_msg=True): ################################################################################################## def update_data_with_plot_inps(data, metadata, inps): + # 在真正绘图前，根据用户参数对数据做参考点改正、单位转换、缠绕显示和数学操作。 # 1. spatial referencing with respect to the seed point if inps.ref_yx: inps.ref_box = [inps.ref_yx[1], inps.ref_yx[0], @@ -244,6 +271,7 @@ def update_data_with_plot_inps(data, metadata, inps): if len(data.shape) == 2: # read ref_val + # 二维图像直接读取参考像元的一个值，然后整幅图减去它。 if 0 <= ref_y < data.shape[-2] and 0 <= ref_x < data.shape[-1]: ref_val = data[ref_y, ref_x] else: @@ -266,6 +294,7 @@ def update_data_with_plot_inps(data, metadata, inps): elif len(data.shape) == 3: # read ref_val + # 三维数据通常是 date x row x col，每个日期都要减去对应日期的参考像元值。 if 0 <= ref_y < data.shape[-2] and 0 <= ref_x < data.shape[-1]: ref_val = np.squeeze(data[:, ref_y, ref_x]) elif inps.key == 'timeseries': @@ -291,6 +320,7 @@ def update_data_with_plot_inps(data, metadata, inps): inps.ref_yx = None # 2. scale data based on the display unit and re-wrap + # 例如把 m 转成 cm，或者把相位按 [-pi, pi] 重新缠绕显示。 (data, inps.disp_unit, inps.disp_scale, @@ -307,6 +337,7 @@ def update_data_with_plot_inps(data, metadata, inps): # math operation if inps.math_operation: + # 用户可对显示数据做简单数学变换，例如开方、取反、弧度转角度。 vprint(f'Apply math operation: {inps.math_operation}') if inps.math_operation == 'square' : data = np.square(data) elif inps.math_operation == 'sqrt' : data = np.sqrt(data) @@ -317,6 +348,7 @@ def update_data_with_plot_inps(data, metadata, inps): else: raise ValueError(f'un-recognized math operation: {inps.math_operation}') # 4. update display min/max + # 自动计算数据范围和颜色显示范围。 inps.dlim = [np.nanmin(data), np.nanmax(data)] if not inps.vlim: # and data.ndim < 3: (inps.cmap_lut, @@ -356,6 +388,7 @@ def prep_slice(cmd, auto_fig=False): plt.show() """ # parse + # prep_slice() 方便其它 Python 代码用一条 view.py 命令字符串准备绘图数据。 from mintpy.cli.view import cmd_line_parse iargs = cmd.split()[1:] @@ -387,12 +420,14 @@ def prep_slice(cmd, auto_fig=False): iargs.append(temp_iarg) # run parse + # 复用 CLI 的参数解析逻辑，保证 Python 调用和命令行调用行为一致。 inps = cmd_line_parse(iargs) global vprint vprint = print if inps.print_msg else lambda *args, **kwargs: None # read input args + # 读取文件基础信息、数据集列表，并根据元数据更新绘图参数。 inps, atr = read_input_file_info(inps) inps = update_inps_with_file_metadata(inps, atr) @@ -407,6 +442,7 @@ def prep_slice(cmd, auto_fig=False): ) # read data + # 只读取当前要绘制的数据集和裁剪区域。 data, atr = readfile.read( inps.file, datasetName=inps.dset[0], @@ -438,6 +474,7 @@ def prep_slice(cmd, auto_fig=False): # masking if inps.zero_mask: + # zero_mask 表示把值为 0 的像元遮住，不显示。 data = np.ma.masked_where(data == 0., data) if inps.msk is not None: @@ -477,6 +514,7 @@ def plot_slice(ax, data, metadata, inps): cbar : matplotlib.colorbar.Colorbar object Example: See prep_slice() above for example usage. """ + # plot_slice() 绘制一张二维图，是 view.py 可视化的核心绘图函数。 global vprint vprint = print if inps.print_msg else lambda *args, **kwargs: None @@ -486,6 +524,7 @@ def extent2meshgrid(extent: tuple, ds_shape: list): shape - list of int for [length, width] of the data Returns: xx/yy - 1D np.ndarray of the data coordinates """ + # scatter 绘图需要每个像元的 x/y 坐标，因此把 extent 和数据形状转换成网格坐标。 height, width = ds_shape x = np.linspace(extent[0], extent[1], width) y = np.linspace(extent[2], extent[3], height)[::-1] # reverse the Y-axis @@ -494,6 +533,7 @@ def extent2meshgrid(extent: tuple, ds_shape: list): # colormap: str -> object if isinstance(inps.colormap, str): + # 把色带名称字符串转换成 matplotlib colormap 对象。 inps.colormap = pp.ColormapExt( inps.colormap, cmap_lut=inps.cmap_lut, @@ -502,6 +542,7 @@ def extent2meshgrid(extent: tuple, ds_shape: list): # read DEM if inps.dem_file: + # 如果指定 DEM，就读取 DEM 作为背景阴影或等高线。 dem, dem_metadata, dem_pix_box = pp.read_dem( inps.dem_file, pix_box=inps.pix_box, @@ -514,6 +555,7 @@ def extent2meshgrid(extent: tuple, ds_shape: list): lalo_digit = ut.get_lalo_digit4display(metadata, coord_unit=inps.coord_unit) # common options for data visualization + # kwargs 是传给 imshow/scatter 的通用绘图参数。 kwargs = dict(cmap=inps.colormap, vmin=inps.vlim[0], vmax=inps.vlim[1], alpha=inps.transparency, zorder=1) #----------------------- Plot in Geo-coordinate --------------------------------------------# @@ -523,6 +565,7 @@ def extent2meshgrid(extent: tuple, ds_shape: list): vprint('plot in geo-coordinate') # extent info for matplotlib.imshow and other functions + # extent 告诉 matplotlib 图像四条边对应的地理坐标范围。 inps.extent = (inps.geo_box[0], inps.geo_box[2], inps.geo_box[3], inps.geo_box[1]) # (W, E, S, N) SNWE = (inps.geo_box[3], inps.geo_box[1], inps.geo_box[0], inps.geo_box[2]) @@ -543,6 +586,7 @@ def extent2meshgrid(extent: tuple, ds_shape: list): ) # Reference (InSAR) data to a GNSS site + # 可选：把 InSAR 图像参考到某个 GNSS 站点附近的像元值。 coord = ut.coordinate(metadata) if inps.disp_gnss and inps.gnss_component and inps.ref_gnss_site: gnss_obj = gnss.get_gnss_class(inps.gnss_source)(site=inps.ref_gnss_site) @@ -558,14 +602,17 @@ def extent2meshgrid(extent: tuple, ds_shape: list): # Plot data if inps.disp_dem_blend: + # DEM blend 会把数据颜色和 DEM 阴影融合显示。 im = pp.plot_blend_image(ax, data, dem, inps, print_msg=inps.print_msg) elif inps.style == 'image': + # image 模式用 imshow 绘制规则栅格，速度最快。 vprint(f'plotting data as {inps.style} via matplotlib.pyplot.imshow ...') im = ax.imshow(data, extent=inps.extent, origin='upper', interpolation=inps.interpolation, animated=inps.animation, **kwargs) elif inps.style == 'scatter': + # scatter 模式逐像元散点绘制，适合不规则显示但速度较慢。 vprint(f'plotting data as {inps.style} via matplotlib.pyplot.scatter (can take some time) ...') xx, yy = extent2meshgrid(inps.extent, data.shape) im = ax.scatter(xx, yy, c=data.flatten(), marker='o', s=inps.scatter_marker_size, **kwargs) @@ -637,6 +684,7 @@ def extent2meshgrid(extent: tuple, ds_shape: list): ax = pp.plot_gnss(ax, SNWE, inps, metadata, print_msg=inps.print_msg) # Status bar + # format_coord() 自定义鼠标悬停/状态栏显示内容：经纬度、像元值、DEM 高程、x/y。 if inps.dem_file: coord_dem = ut.coordinate(dem_metadata) dem_len, dem_wid = dem.shape @@ -668,6 +716,7 @@ def format_coord(x, y): #------------------------ Plot in Y/X-coordinate ------------------------------------------------# else: + # 非地理坐标绘图：直接使用图像行列号 x/y 作为坐标。 inps.fig_coord = 'yx' vprint('plotting in Y/X coordinate ...') @@ -683,6 +732,7 @@ def format_coord(x, y): ) # extent = (left, right, bottom, top) in data coordinates + # y 轴使用图像行号，origin 在左上角，因此 top/bottom 顺序与地理图不同。 inps.extent = (inps.pix_box[0]-0.5, inps.pix_box[2]-0.5, inps.pix_box[3]-0.5, inps.pix_box[1]-0.5) @@ -742,6 +792,7 @@ def format_coord(x, y): # Status bar # read lats/lons if exist + # 雷达坐标图也可以尝试从 geometryRadar.h5 读取 latitude/longitude，在状态栏显示经纬度。 geom_file = os.path.join(os.path.dirname(metadata['FILE_PATH']), 'inputs/geometryRadar.h5') if os.path.isfile(geom_file): geom_ds_list = readfile.get_dataset_list(geom_file) @@ -783,6 +834,7 @@ def format_coord(x, y): # 3.1 Colorbar cbar = None if inps.disp_cbar: + # 在主图旁边创建一个 colorbar 轴，并绘制颜色条。 divider = make_axes_locatable(ax) cax = divider.append_axes(inps.cbar_loc, inps.cbar_size, pad=inps.cbar_size, axes_class=plt.Axes) inps, cbar = pp.plot_colorbar(inps, im, cax) @@ -835,6 +887,7 @@ def format_coord(x, y): ################################################################################################## def read_input_file_info(inps): + # 读取输入文件基础信息：文件类型、尺寸、数据集列表，并确定要显示哪些数据集。 # File Basic Info atr = readfile.read_attribute(inps.file) msg = 'input file is ' @@ -856,6 +909,7 @@ def read_input_file_info(inps): vprint(f'file size in y/x: {(inps.length, inps.width)}') # File dataset List + # get_slice_list() 返回可以被 view.py 绘制的二维切片名称列表。 inps.sliceList = readfile.get_slice_list(inps.file, no_complex=True) # Read input list of dataset to display @@ -866,6 +920,7 @@ def read_input_file_info(inps): def search_dataset_input(all_list, in_list=[], in_num_list=[], search_dset=True): """Get dataset(es) from input dataset / dataset_num""" + # 这个函数根据用户输入的数据集名称/编号，在文件所有数据集列表中找到匹配项。 # make a copy to avoid weird variable behavior in_num_list = [x for x in in_num_list] @@ -878,6 +933,7 @@ def search_dataset_input(all_list, in_list=[], in_num_list=[], search_dset=True) if search_dset: for ds in in_list: # style of regular expression + # 用户可以输入 velocity 或 *velocity* 这类模式，这里转换成正则表达式搜索。 if '*' not in ds: ds = f'*{ds}*' ds = ds.replace('*','.*') @@ -901,6 +957,7 @@ def search_dataset_input(all_list, in_list=[], in_num_list=[], search_dset=True) def read_dataset_input(inps): """Check input / exclude / reference dataset input with file dataset list""" + # 根据用户输入、文件类型和默认规则，确定最终要绘制的数据集列表 inps.dset。 # read inps.dset + inps.dsetNumList --> inps.dsetNumList if len(inps.dset) > 0 or len(inps.dsetNumList) > 0: # message @@ -909,10 +966,12 @@ def read_dataset_input(inps): # special rule for special file types if inps.key == 'velocity': + # velocity 文件的数据集名通常很明确，关闭模糊搜索避免误匹配。 inps.search_dset = False vprint(f'turning glob search OFF for {inps.key} file') elif inps.key == 'timeseries' and len(inps.dset) == 1 and '_' in inps.dset[0]: + # 支持 view.py timeseries.h5 20180101_20180201 这种写法：显示后者相对前者。 date1, date2 = inps.dset[0].split('_') inps.dset = [date2] inps.ref_date = date1 @@ -926,6 +985,7 @@ def read_dataset_input(inps): else: # default dataset to display for certain type of files + # 用户没有指定数据集时，根据文件类型选择默认绘制内容。 if inps.key == 'ifgramStack': inps.dset = ['unwrapPhase'] elif inps.key == 'HDFEOS': @@ -956,6 +1016,7 @@ def read_dataset_input(inps): search_dset=inps.search_dset) # read inps.plot_drop_ifgram + # ifgramStack 中被 modify_network 丢弃的干涉图，默认不显示。 drop_num_list = [] ftype = readfile.read_attribute(inps.file)['FILE_TYPE'] @@ -1019,11 +1080,13 @@ def update_figure_setting(inps): 2) for multi: figure/row/column number 3) for multi: output file name """ + # 根据数据集数量自动设置单图/多子图布局、字体大小、图像大小和输出文件名。 length = float(inps.pix_box[3]-inps.pix_box[1]) width = float(inps.pix_box[2]-inps.pix_box[0]) # One Plot if inps.dsetNum == 1: + # 单图时优先保证图像比例和 colorbar 空间合适。 if not inps.font_size: inps.font_size = 16 if not inps.fig_size: @@ -1039,6 +1102,7 @@ def update_figure_setting(inps): # Multiple Plots else: + # 多数据集时自动计算行列数，必要时分成多张 figure。 if not inps.fig_size: inps.fig_size = pp.default_figsize_multi vprint(f'figure size : [{inps.fig_size[0]:.2f}, {inps.fig_size[1]:.2f}]') @@ -1117,6 +1181,7 @@ def read_data4figure(i_start, i_end, inps, metadata): """Read multiple datasets for one figure into 3D matrix based on i_start/end""" # initiate output matrix + # data 形状为 子图数量 x 行 x 列；每个子图对应一个二维数据集。 data = np.zeros(( i_end - i_start, int((inps.pix_box[3] - inps.pix_box[1]) / inps.multilook_num), @@ -1137,6 +1202,7 @@ def read_data4figure(i_start, i_end, inps, metadata): ref_kwargs = dict(box=(ref_x, ref_y, ref_x+1, ref_y+1), print_msg=False) # fast reading for single dataset type + # 如果所有子图属于同一种数据集家族，可以一次性读成 3D 矩阵，速度更快。 if (len(inps.dsetFamilyList) == 1 and inps.key in ['timeseries', 'ifgramStack', 'geometry', 'HDFEOS', 'giantTimeseries']): @@ -1158,6 +1224,7 @@ def read_data4figure(i_start, i_end, inps, metadata): # slow reading with one 2D matrix at a time else: + # 数据集类型混杂时，逐个读取二维矩阵更稳妥。 vprint('reading data as a list of 2D matrices ...') kwargs['print_msg'] = False prog_bar = ptime.progressBar(maxValue=i_end-i_start, print_msg=inps.print_msg) @@ -1178,6 +1245,7 @@ def read_data4figure(i_start, i_end, inps, metadata): # ref_date for timeseries if inps.ref_date: + # 对时序数据，可以显示某个日期相对于参考日期的差值。 vprint('consider input reference date: '+inps.ref_date) ref_data = readfile.read(inps.file, datasetName=inps.ref_date, **kwargs)[0] data -= ref_data @@ -1197,10 +1265,12 @@ def read_data4figure(i_start, i_end, inps, metadata): or all(d.endswith('Phase') for d in inps.dsetFamilyList)): same_unit4all_subplots = True else: + # 如果多个子图单位不同，就不能统一做单位转换和统一 colorbar 范围。 same_unit4all_subplots = False # adjust data due to spatial referencing and unit related scaling if same_unit4all_subplots: + # 单位一致时，可以对整个 3D data 一次性做参考点/单位/缠绕处理。 data, inps = update_data_with_plot_inps(data, metadata, inps) else: if any(x in inps.argv for x in ['-u', '--unit']): @@ -1211,6 +1281,7 @@ def read_data4figure(i_start, i_end, inps, metadata): # mask if inps.zero_mask: + # 多子图同样支持把 0 值遮住。 vprint('masking pixels with zero value') data = np.ma.masked_where(data == 0., data) if inps.msk is not None: @@ -1219,6 +1290,7 @@ def read_data4figure(i_start, i_end, inps, metadata): data = np.ma.masked_where(msk == 0., data) # update display min/max + # 多子图如果共享单位且用户没有手动设置 -v，就自动计算统一颜色范围。 inps.dlim = [np.nanmin(data), np.nanmax(data)] if (same_unit4all_subplots and all(arg not in inps.argv for arg in ['-v', '--vlim', '--wrap']) @@ -1239,6 +1311,7 @@ def plot_subplot4figure(i, inps, ax, data, metadata): 1) Plot DEM, data and reference pixel 2) axes setting: tick, ticklabel, title, axis etc. """ + # plot_subplot4figure() 负责多子图中的一个小图：画背景、画数据、画参考点、设置标题和坐标轴。 # Plot DEM if inps.dem_file: pp.plot_dem_background( @@ -1251,6 +1324,7 @@ def plot_subplot4figure(i, inps, ax, data, metadata): print_msg=inps.print_msg) # Plot Data + # 多子图使用 imshow 绘制，每个子图是一个二维矩阵。 vlim = inps.vlim if inps.vlim is not None else [np.nanmin(data), np.nanmax(data)] inps.extent = (inps.pix_box[0]-0.5, inps.pix_box[2]-0.5, inps.pix_box[3]-0.5, inps.pix_box[1]-0.5) @@ -1260,6 +1334,7 @@ def plot_subplot4figure(i, inps, ax, data, metadata): # Plot Seed Point if inps.disp_ref_pixel: + # 在每个子图上标出参考像元位置。 ref_y, ref_x = None, None if inps.ref_yx: ref_y, ref_x = inps.ref_yx[0], inps.ref_yx[1] @@ -1286,6 +1361,7 @@ def format_coord(x, y): # Title if inps.disp_title: # get title + # 时序类数据集通常用日期作为标题，其它数据集用数据集名或索引作为标题。 subplot_title = None if inps.key in TIMESERIES_KEY_NAMES or inps.dset[0].startswith('bperp'): # support "/" in the dataset names, e.g. HDF-EOS5 and py2-mintpy formats @@ -1304,6 +1380,7 @@ def format_coord(x, y): if len(inps.dsetFamilyList) == 1 and '-' in title_str: title_str = title_str.split('-')[1] # for ifgramStack, show index in the date12 list to facilitate the network modification + # 对 ifgramStack，标题中显示索引有助于回到 modify_network 用编号删除干涉图。 if inps.atr['FILE_TYPE'] == 'ifgramStack': title_ind = inps.date12List.index(title_str) @@ -1326,6 +1403,7 @@ def format_coord(x, y): else: kwarg = dict(fontsize=inps.font_size) # mark dropped interferograms in bold crimson + # 被 dropIfgram=False 标记删除的干涉图标题显示为红色粗体。 if inps.dset[i] in inps.dropDatasetList: kwarg['color'] = 'crimson' kwarg['fontweight'] = 'bold' @@ -1359,11 +1437,13 @@ def plot_figure(j, inps, metadata): 3) loop to plot each subplot using plot_subplot4figure() 4) common colorbar and save """ + # plot_figure() 负责生成一整张 figure，里面可能包含很多 subplot。 fig_title = f'Figure {str(j)} - {inps.outfile[j-1]}' vprint('----------------------------------------') vprint(fig_title) # Open a new figure object + # plt.subplots() 创建 figure 和子图网格。 fig, axs = plt.subplots( num=j, figsize=inps.fig_size, nrows=inps.fig_row_num, @@ -1374,6 +1454,7 @@ def plot_figure(j, inps, metadata): axs = axs.flatten() # Read all data for the current figure into 3D np.array + # 每张 figure 只读取自己需要显示的那一段数据集，避免一次读入所有子图。 i_start = (j - 1) * inps.fig_row_num * inps.fig_col_num i_end = min([inps.dsetNum, i_start + inps.fig_row_num * inps.fig_col_num]) data = read_data4figure(i_start, i_end, inps, metadata) @@ -1384,6 +1465,7 @@ def plot_figure(j, inps, metadata): ).colormap # Loop - Subplots + # 逐个子图调用 plot_subplot4figure()。 vprint('plotting ...') prog_bar = ptime.progressBar(maxValue=i_end-i_start, print_msg=inps.print_msg) for i in range(i_start, i_end): @@ -1409,6 +1491,7 @@ def plot_figure(j, inps, metadata): del data # delete empty axes + # 最后一张 figure 可能没有填满所有子图位置，需要删除空白 axes。 for i in range(i_end-i_start, len(axs)): fig.delaxes(axs[i]) @@ -1423,6 +1506,7 @@ def plot_figure(j, inps, metadata): # before fig.add_axes(), which is the case of common colorbar # after fig.colorbar() and fig.set_size_inches(), which is the case of individual/multiple colorbars def adjust_subplots_layout(fig, inps): + # 调整子图间距，减少多子图中的空白区域。 fig.subplots_adjust( left=0.02, right=0.98, bottom=0.02, top=0.98, @@ -1438,12 +1522,14 @@ def adjust_subplots_layout(fig, inps): # Colorbar if inps.disp_cbar: if not inps.vlim: + # 没有统一 vlim 时，每个子图都有自己的颜色范围和 colorbar。 vprint('Note: different color scale for EACH subplot!') vprint('Adjust figsize for the colorbar of each subplot.') fig.set_size_inches(inps.fig_size[0] * 1.1, inps.fig_size[1]) adjust_subplots_layout(fig, inps) else: + # 有统一 vlim 时，整张 figure 共用一个 colorbar。 adjust_subplots_layout(fig, inps) # plot common colorbar cbar_length = 0.4 @@ -1459,6 +1545,7 @@ def adjust_subplots_layout(fig, inps): # Save Figure if inps.save_fig: + # 保存当前 figure 到磁盘；不显示窗口时，保存后清空 figure 释放内存。 vprint(f'save figure to {os.path.abspath(inps.outfile[j-1])} with dpi={inps.fig_dpi}') fig.savefig(inps.outfile[j-1], bbox_inches='tight', transparent=True, dpi=inps.fig_dpi) if not inps.disp_fig: @@ -1473,6 +1560,7 @@ def prepare4multi_subplots(inps, metadata): 3) read dropIfgram info 4) read and prepare DEM for background """ + # 多子图绘制前的准备：确定数据集家族、降采样、参考像元、被删除干涉图、DEM 背景等。 inps.dsetFamilyList = sorted(list({x.split('-')[0] for x in inps.dset})) inps.dsetFamilyList = sorted(list({x.replace('Std','') for x in inps.dsetFamilyList})) if len(inps.dsetFamilyList) == 1 and inps.atr['FILE_TYPE'] == 'ifgramStack': @@ -1495,6 +1583,7 @@ def prepare4multi_subplots(inps, metadata): # multilook mask if inps.msk is not None and inps.multilook_num > 1: + # 如果数据降采样了，掩膜也必须用同样倍率降采样，尺寸才能匹配。 inps.msk = multilook_data( inps.msk, inps.multilook_num, @@ -1503,6 +1592,7 @@ def prepare4multi_subplots(inps, metadata): ) # Reference pixel for timeseries and ifgramStack + # 对 ifgramStack 的 unwrapPhase，多子图显示时要考虑文件已有参考像元。 inps.file_ref_yx = None if inps.key in ['ifgramStack'] and 'REF_Y' in metadata.keys(): ref_y, ref_x = int(metadata['REF_Y']), int(metadata['REF_X']) @@ -1517,6 +1607,7 @@ def prepare4multi_subplots(inps, metadata): inps.ref_marker_size /= 20. # Check dropped interferograms + # 多子图显示 ifgramStack 时，被网络筛选删除的干涉图会在标题中标红。 inps.dropDatasetList = [] if inps.key == 'ifgramStack' and inps.disp_title: obj = ifgramStack(inps.file) @@ -1528,6 +1619,7 @@ def prepare4multi_subplots(inps, metadata): # Read DEM if inps.dem_file: + # 多子图只支持 DEM 与数据文件尺寸一致的情况；否则 DEM 背景会被忽略。 dem_metadata = readfile.read_attribute(inps.dem_file) if all(dem_metadata[i] == metadata[i] for i in ['LENGTH', 'WIDTH']): vprint(f'reading DEM: {os.path.basename(inps.dem_file)} ... ') @@ -1570,32 +1662,39 @@ class viewer(): """ def __init__(self, cmd=None, iargs=None): + # viewer 是面向外部调用的类：可以传命令字符串 cmd，也可以传已经拆好的参数列表 iargs。 if cmd: iargs = cmd.split()[1:] self.cmd = cmd self.iargs = iargs def configure(self, inps): + # configure() 做绘图前准备：读取文件信息、更新参数、判断是否跳过、读取掩膜。 global vprint vprint = print if inps.print_msg else lambda *args, **kwargs: None # matplotlib backend setting if not inps.disp_fig: + # Agg 是无界面绘图后端，适合只保存图片、不弹窗的批处理。 plt.switch_backend('Agg') + # 读取文件元数据和数据集列表，并根据元数据补全绘图配置。 inps, self.atr = read_input_file_info(inps) inps = update_inps_with_file_metadata(inps, self.atr) # --update option self.flag = 'run' if inps.update_mode and run_or_skip(inps) == 'skip': + # 如果 update 模式判断可以跳过，plot() 外层通常就不会继续绘图。 self.flag = 'skip' # copy inps to self object + # setattr(self, key, value) 把 argparse 参数复制成 viewer 对象属性，后面用 self.xxx 访问。 for key, value in inps.__dict__.items(): setattr(self, key, value) # read mask + # 读取用户指定或默认掩膜，后续绘图时遮住无效像元。 self.msk, self.mask_file = pp.read_mask( self.file, mask_file=self.mask_file, @@ -1609,10 +1708,12 @@ def configure(self, inps): def plot(self): + # plot() 是真正绘图的入口：根据 dsetNum 决定画单图还是多子图。 # One Subplot if self.dsetNum == 1: vprint('reading data ...') # read data + # 单图只读取一个数据集和裁剪范围。 data, self.atr = readfile.read( self.file, datasetName=self.dset[0], @@ -1621,6 +1722,7 @@ def plot(self): # reference in time if self.ref_date: + # 时序单图支持减去参考日期，显示两个日期之间的变化。 data -= readfile.read( self.file, datasetName=self.ref_date, @@ -1631,6 +1733,7 @@ def plot(self): if (self.key in ['ifgramStack'] and self.dset[0].split('-')[0].startswith('unwrapPhase') and 'REF_Y' in self.atr.keys()): + # ifgramStack 的 unwrapPhase 需要减去参考像元相位。 ref_y, ref_x = int(self.atr['REF_Y']), int(self.atr['REF_X']) ref_data = readfile.read( self.file, @@ -1641,6 +1744,7 @@ def plot(self): # masking - input options if self.zero_mask: + # 按用户选项遮住 0 值像元。 vprint('masking pixels with zero value') data = np.ma.masked_where(data == 0., data) if self.msk is not None: @@ -1650,6 +1754,7 @@ def plot(self): self.msk = np.ones(data.shape, dtype=np.bool_) # masking - NO_DATA_VALUE + # 如果文件定义了 NO_DATA_VALUE，也把这些像元转成 NaN，不参与显示。 no_data_val = readfile.get_no_data_value(self.file) if self.no_data_value is not None: vprint(f'masking pixels with NO_DATA_VALUE of {self.no_data_value}') @@ -1672,15 +1777,18 @@ def plot(self): self.msk *= ~np.isnan(data) # update data + # 做单位转换、参考点改正、数学操作和颜色范围更新。 data, self = update_data_with_plot_inps(data, self.atr, self) # prepare figure + # 如果启用了 cartopy 地图投影，就把 projection 传给 plt.subplots()。 subplot_kw = dict(projection=self.map_proj_obj) if self.map_proj_obj is not None else {} fig, ax = plt.subplots(figsize=self.fig_size, subplot_kw=subplot_kw) if not self.disp_whitespace: fig.subplots_adjust(left=0,right=1,bottom=0,top=1) # plot + # 调用前面定义的 plot_slice() 绘制单张二维图。 self = plot_slice(ax, data, self.atr, self)[1] # save figure @@ -1695,6 +1803,7 @@ def plot(self): # Multiple Subplots else: + # 多子图模式用于一次查看多个日期、多个干涉图或多个数据集。 # warn single-subplot options opt_names = ['--show-gnss', '--coastline', '--lalo-label', '--lalo-step', '--scalebar', '--pts-yx', '--pts-lalo', '--pts-file'] @@ -1703,9 +1812,11 @@ def plot(self): print(f'WARNING: {opt_name} is NOT supported for multi-subplots, ignore it and continue.') # prepare + # 多子图准备包括自动降采样、DEM 背景、dropIfgram 标题标记等。 self = prepare4multi_subplots(self, metadata=self.atr) # plot + # 可能有多张 figure，因此循环调用 plot_figure()。 self.dlim_all = [0., 0.] for j in range(1, self.fig_num + 1): plot_figure(j, self, metadata=self.atr) @@ -1719,6 +1830,7 @@ def plot(self): # Display Figure if self.disp_fig: + # 如果用户要求显示图形窗口，最后调用 plt.show()。 vprint('showing ...') plt.show() return