首页 / 资讯中心 / 文章详情

电路理论与Omniscape生态廊道连通性批量分析实践

电路理论与Omniscape生态廊道连通性批量分析实践 ★ FEATURED ARTICLE
简介压缩包面向生态学、景观连通性与保护规划方向的研究者、规划师及学生围绕电路理论模拟物种迁移和栖息地连通性解决生态廊道识别、瓶颈区定位与障碍阻隔评估等问题。包内共264个文件、约95.49MB含Python脚本、地理数据库表与索引、ADF栅格、MXD工程、PDF说明及SHP矢量数据分别用于批量处理、空间数据组织、制图展示和结果核查。已有2612人学习下载具备较高参考价值。覆盖Linkage Mapper的廊道识别与连通性制图、Pinchpoint Mapper的关键“掐点”保护分析、Barriers Mapper的障碍量化并涉及物种分布、地形高程到栅格/矢量输入的完整流程帮助读者理解电路理论如何转化为连通性矩阵和地图为保护区布局、生态走廊规划及障碍修复提供可复用的工具与案例参照。1. 从一张阻力面到一张电流密度图这套包解决的是“多条路都算”的问题做生态廊道的人常被一个问题噎住最小成本路径给你一条线评审会上有人问“为什么不是旁边那条”你很难答。基于电路理论的 Circuitscape 换了个问法——把景观当成电阻面把栖息地斑块当成电极通上电看电流在哪儿密。电流密度高的地方就是个体移动概率高的地方而且是所有路径同时解出来的不是一条最优线。这套Circuitscape软件包集.zip收的正是这么一套东西计算引擎、Python 接口、图形界面外加一个做全景连通性的 Omniscape。它适合两类人手上已经有阻力面、想把廊道结论从“一条线”升级成“一张概率图”的从业者以及要批量跑几十组情景、必须靠脚本而不是点鼠标的人。2. 解包后先认清四件套引擎、接口、界面和 Omniscape 各管什么压缩包解开之后不要急着双击。这套东西不是一个大程序而是四个定位完全不同的部件拼在一起的。搞混了就会出现“我在界面上调好的参数脚本里跑出来完全不一样”这种事——不是软件有 bug是你用错了入口。下面按“谁在什么场景下该用哪个”拆开讲再落到三份输入数据的合格标准最后跑通一个最小算例。2.1 四个部件的分工与选择部件形态什么场景用依赖与注意计算引擎命令行入口可用python -m circuitscape 配置.ini调起服务器批量跑、接调度系统、复现实验需要 Python 环境和包内依赖路径别带空格和中文Python 接口circuitscape模块在自己的脚本里循环、读结果、做统计与引擎同源参数完全一致图形界面桌面程序第一次摸参数、检查输入对不对、单次出图生成的 ini 是极好的模板直接抄Omniscape独立程序区域尺度、不分源汇的全景连通运行时独立配置格式和 Circuitscape 不是一套我一般的工作流是先在界面里用一份小范围数据把参数试对把界面生成的 ini 存下来当模板然后把模板里的路径改成变量扔进 Python 脚本批量跑只有需要全景图的时候才切到 Omniscape。这样既能享受界面的即时反馈又不会把几十组情景的重复劳动压在鼠标上。有一点要提前说清楚包里各组件的具体目录名、启动方式不同打包版本会有出入。别硬套网上看来的路径先看解压目录里有没有 README 或示例配置照着示例改比自己猜快得多。2.2 三份输入数据的最低合格标准Circuitscape 的输入只有三样一张阻力面、一份节点文件、一个配置文件。听起来简单翻车基本都翻在这三样上。阻力面。必须是数值型栅格而且值是“阻力”不是“适宜性”。值越大表示越难穿过这是反直觉的地方——很多做物种分布的人手上拿的是适宜性图层值越大越好直接喂进去结果会整个反过来廊道全落在本该避开的区域。如果只有适宜性数据常见做法是做一次单调递减变换再转成阻力比如取倒数或做线性反转变换方式要在报告里写清楚因为它直接决定结论。坐标系也值得单独提一句。优先用投影坐标系像元单位是米用地理坐标系单位是度时八邻域的对角线距离换算会失真横竖方向一度和斜向一度的实际地面距离不一样电阻累积就偏了。分辨率则决定了结论的粒度30 米和 300 米跑出来的廊道宽度不可比跨研究对比时这是硬伤。NoData 是第三个雷。多来源图层拼接之后NoData 值经常是 -9999、255、nan 混着用软件只认配置里声明的那个。没统一的话要么整片区域被当成“极高阻力”把电流全挡死要么被当成 0 阻力变成一条假的高速通道。节点文件。纯文本两列x y或三列ID x y分隔符空格或制表符都行。三个硬性要求坐标必须和栅格同一坐标系节点不能落在 NoData 像元上任意两个节点不能落进同一个像元——最后这条最容易被忽略落进同一像元后两点间电阻是 0输出矩阵里一堆零你还以为是软件坏了。另外use_pixel_centers这个开关要留意。它决定节点坐标被解释成像元中心还是像元角点。节点是从栅格导出来的就设 True是从别的矢量数据来的、又做过坐标系转换就得实际验证一下落点对不对。配置文件。ini 格式分节。核心是四节模式、栅格、连接方案、输出。下面直接给一份可改的。2.3 用 pairwise 模式跑通第一个闭环先把模式定在 pairwise——所有节点两两配对每个配对单独解一次电流。这是最容易验证的模式结果也最好解释。[Circuitscape mode] data_type raster scenario pairwise [Habitat raster or graph] habitat_file /data/proj/inputs/resistance_30m.tif habitat_map_is_resistances True [Options for pairwise and one-to-all and all-to-one modes] point_file /data/proj/inputs/nodes.txt use_included_pairs False [Connection scheme for raster habitat data] connect_using_avg_resistances True connect_four_neighbors_only False [Options for raster habitat data] use_pixel_centers True [Output options] output_file /data/proj/out/scenario_A.out write_cumulative_current True write_current True write_volt_maps False set_null_currents_to_nodata True compress_grids False [Options for parallel processing] parallelize True max_parallel 8几个键值得解释。habitat_map_is_resistances True声明栅格值是阻力如果你的图层是电导值越大越好走这里要改成 False否则方向整个翻过来。connect_four_neighbors_only False表示用八邻域电流可以走对角线更接近真实的扩散改成 True 就只能走横竖廊道会显得方正保守粗分辨率数据上我倾向用四邻域抑制那种贴着对角线钻出来的假捷径。parallelize配max_parallel开并行但要注意瓶颈通常不是 CPU 而是内存核数开太大反而容易被系统杀掉。跑起来python -m circuitscape /data/proj/cfg/scenario_A.ini跑完之后以output_file的主名为前缀会落下一组东西scenario_A_resistances.out是节点两两之间的电阻距离矩阵这是定量分析的主表同目录下还有一个电流图目录里面是每对节点的电流密度栅格再加一张累计电流图把所有配对的电流叠起来。如果只想要电阻矩阵、不要那一堆栅格把write_current关掉磁盘能省下一大截。第一次跑建议只放三到五个节点跑完先看电阻矩阵是不是对称、对角是不是 0。这两个检查花不了一分钟能挡掉后面八成的排查时间。3. 把一次点击变成一条流水线Python 侧批量驱动与结果解析界面能跑通不代表能交付。真实项目里往往是十几个阻力情景乘上几组节点配置再乘上敏感性分析的几个变体几十次计算起步。这个量级靠点鼠标既慢又不可复现——三个月后你根本说不清当时那个图是用哪一版阻力面跑出来的。这一章把 ini 生成、调用、结果解析串成一条能直接抄的流水线。3.1 为什么批处理必须绕开图形界面界面的问题不只是慢。它把参数状态保存在窗口里而不是保存成一个可以进版本库的文本文件。你调了五次参数最后只有最后一次的结果中间那四次的过程信息丢了。而脚本化的核心价值就是配置文件即实验记录。还有两个现实原因。一是内存。大栅格跑 pairwise单次计算的内存峰值可能好几个 G用 subprocess 起独立进程进程退出内存立刻还给系统在同一个 Python 进程里反复调 API内存碎片会一点点攒起来跑到第十几个情景就 OOM。二是失败隔离某个情景因为节点落进 NoData 挂了subprocess 抛异常你捕获掉就行不会把整个批处理进程带崩。所以我的习惯是用 Python 生成配置、用 subprocess 调用、用 numpy 读结果。三段解耦每段都能单独测。3.2 用模板生成 ini 再逐场景调用先写一个模板函数。核心思路是把上一章那份 ini 变成带占位符的字符串路径一律resolve()成绝对路径——相对路径是批量脚本里最隐蔽的坑工作目录一变读到的就是另一份文件。from pathlib import Path import subprocess import sys INI_TMPL [Circuitscape mode] data_type raster scenario pairwise [Habitat raster or graph] habitat_file {habitat} habitat_map_is_resistances True [Options for pairwise and one-to-all and all-to-one modes] point_file {points} use_included_pairs False [Connection scheme for raster habitat data] connect_using_avg_resistances True connect_four_neighbors_only {four_nb} [Options for raster habitat data] use_pixel_centers True [Output options] output_file {out} write_cumulative_current True write_current {write_cur} write_volt_maps False set_null_currents_to_nodata True compress_grids False [Options for parallel processing] parallelize True max_parallel {nproc} def build_ini(ini_path, habitat, points, out, four_nbFalse, write_curTrue, nproc4): ini_path Path(ini_path) ini_path.parent.mkdir(parentsTrue, exist_okTrue) ini_path.write_text(INI_TMPL.format( habitatPath(habitat).resolve(), pointsPath(points).resolve(), outPath(out).resolve(), four_nbstr(four_nb), write_curstr(write_cur), nprocint(nproc), ), encodingutf-8) return ini_pathfour_nb和write_cur用字符串True/False而不是 Python 布尔值因为 ini 解析器读的是字面量。nproc建议先给 4跑通再往上加加到机器开始频繁换页就说明过头了。resolve()那一行不是洁癖是保命。然后是调用器。这里用 subprocess 而不是直接 import理由上面说了内存回收干净、失败可隔离。def run_cs(ini_path, timeout7200): cmd [sys.executable, -m, circuitscape, str(Path(ini_path).resolve())] r subprocess.run(cmd, capture_outputTrue, textTrue, timeouttimeout) if r.returncode ! 0: # 只留尾部报错信息通常就在最后几百行 raise RuntimeError(ffailed: {ini_path}\n{r.stderr[-3000:]}) return r.stdouttimeout一定要给。我遇到过一次栅格范围填错软件在解一个几乎无穷大的连通区域进程安安静静挂了六个小时日志一个字没有。给了超时之后这种事最多浪费一次设定值。批量的主循环就很简单了把情景表写成字典列表一个个过scenarios [ {name: A_4nb, habitat: inputs/res_A.tif, four_nb: True}, {name: A_8nb, habitat: inputs/res_A.tif, four_nb: False}, {name: B_8nb, habitat: inputs/res_B.tif, four_nb: False}, ] for sc in scenarios: ini build_ini(fcfg/{sc[name]}.ini, sc[habitat], inputs/nodes.txt, fout/{sc[name]}.out, four_nbsc[four_nb], nproc4) run_cs(ini)注意每个情景的输出主名必须不同否则后一次会静默覆盖前一次——这是我最不想承认踩过的一个坑跑了半天发现只有最后一个情景的结果。3.3 解析电阻矩阵与电流图结果文件有两种形态。矩阵形态是方阵首行和首列是节点 ID三列形态每行是节点A 节点B 电阻解析起来更省事做统计分析我基本都用它。import numpy as np def load_pairs(path, skip_headerFalse): 读 *_resistances_3columns.out返回 {(a,b): r} 字典 arr np.loadtxt(path, ndmin2, skiprows1 if skip_header else 0, comments#) out {} for a, b, r in arr[:, :3]: key (int(a), int(b)) if a b else (int(b), int(a)) out[key] float(r) return out def to_matrix(pairs, ids): 把字典摊回方阵方便做对称性检查 idx {nid: i for i, nid in enumerate(ids)} m np.zeros((len(ids), len(ids))) for (a, b), r in pairs.items(): m[idx[a], idx[b]] m[idx[b], idx[a]] r return mskip_header这个参数留着是有用的不同版本写出来的三列文件有的带一行列名有的不带读之前用head -2看一眼比猜快。键统一排序成(小, 大)是为了后面做汇总时不重复计数。电流图那边是标准栅格用常规的栅格读取库打开就行。真正要注意的是累计电流图的量纲——它是所有配对电流的简单叠加节点数一多数值就跟着涨所以跨情景比较之前必须先归一化除以配对数或者除以自身的分位数。直接拿两张累计电流图的绝对值比高低是无效比较。def normalize_current(arr, n_pairs): 累计电流除以配对数得到可比的平均电流密度 return np.where(np.isfinite(arr), arr / max(n_pairs, 1), np.nan)配对数就是n*(n-1)/2n 是节点数。这一步做完不同节点配置、不同情景的电流图才站在同一条起跑线上。3.4 参数改动对照表参数不是随便调的每个都有自己的物理含义。下面这张表是我实际会动的几个其余保持默认。参数常见取值改了会怎样什么情况下动它connect_four_neighbors_onlyFalse八邻域True 时电流只走横竖廊道更方正、更保守分辨率粗、想压掉对角假捷径connect_using_avg_resistancesTrueFalse 时取相邻两格里更好走的那个值整体连通性偏高一般不动属于少数派用法use_pixel_centersTrue改变节点坐标的解释方式节点来自矢量、做过投影转换时要验证write_currentTrue关掉只留电阻矩阵磁盘省很多大批量、只要定量指标时set_null_currents_to_nodataTrue背景 0 变 NoData出图不糊成一片出图必开max_parallel视机器加到内存打满为止不是越多越快大栅格先看内存再看核数connect_using_avg_resistances这一行值得多说两句。它决定相邻两格之间的连接电阻取平均值还是取较优值。取平均值更接近“整个像元都在均匀阻挡”取较优值相当于假设个体总能找到格子里最好走的那一小块连通性会系统性偏高。两者都能用但必须在方法里写清楚不然不同研究之间没法比。4. Omniscape 的滑动窗口不分源汇的全景电流怎么算pairwise 模式有个隐含前提你得先知道哪些斑块是源、哪些是汇。区域尺度的规划项目里这个前提经常不成立——你手上只有一张全区的阻力面不知道哪块是“出发点”。Omniscape 就是为这个场景准备的它在一个滑动窗口里把窗口内所有低阻力像元都当成潜在源逐块算电流再把所有窗口的结果叠起来。不需要预先指定节点。4.1 pairwise 在区域尺度上的两个硬伤第一个是组合数。n 个节点要算 n(n-1)/2 次五十个节点就是一千两百多次计算每个节点对都要解一遍全网电流时间随节点数平方增长。省域尺度的斑块动辄上百个这条路走不通。第二个更本质源汇假设本身就是结论的一部分。你指定了哪几块是源结果就绕着这几块转。换个专家来指定结论可能就变了。而廊道规划要回答的恰恰是“哪一带整体上更容易连通”不该被少数几个节点的选择绑架。Omniscape 的解法是把“源”这件事弱化成一种连续权重每个像元都有一个源强度可以是二值的也可以是栖息地质量归一化后的连续值。窗口滑动时窗口内强度高的像元自然承担更多电流注入。这样既保留了电路理论的物理意义又不需要拍板指定节点。代价是计算量。滑动窗口意味着每个像元被重复计算很多次参数没调好能把机器跑死。4.2 半径、块大小、缓冲区怎么定这三个参数是 Omniscape 的全部难点其余都是常规配置。radius窗口半径单位像元。它等于你设定的最大扩散距离。填小了廊道只在一个街区范围内成立长距离扩散被截断填大了相当于假设个体一步能跨过半个省算出来的电流会把所有区域糊成一片没有分辨力。定这个值的依据应该是目标物种的扩散能力而不是“看起来好看”。我的经验是先按文献里该物种的最大扩散距离折成像元数再用一半和两倍各跑一次做敏感性看结论稳不稳。block_size分块边长。控制单次计算的内存占用。窗口滑动本身是逐像元推进的但实现上会切块并行块太大内存炸块太小调度开销上来。buffer块间重叠宽度。必须大于等于 radius否则相邻两块的边缘会各自按自己的边界条件解一遍拼起来就出现一条条规律的高电流或低电流缝——这个伪影非常像“廊道”第一次见很容易当成重要发现实际上是算法边界。配置大概长这样字段名以包内示例为准我这里给的是结构[Inputs] resistances_file resistance_30m.tif source_file source_strength.tif [Outputs] output_dir out/omni_A [Omniscape] radius 40 block_size 12 buffer 60 mode rank source_threshold 0.0 calc_parallel true write_raw false correct_artifacts truemode一般用rank按窗口内电流排序做归一化对极端值不敏感跨区域可比性好。correct_artifacts保持开启它专门处理窗口边缘的伪影。source_threshold用来把低于阈值的源强度当 0避免大片低质量像元持续注入电流把背景抬高。4.3 输出后处理从原始电流到可比强度Omniscape 输出的原始累计电流图不能直接当廊道图用。它是每个窗口结果的叠加数值大小和窗口数、半径都强相关换个参数就不可比。两件事必须做。第一是归一化。按窗口覆盖次数或者直接按分位数折算让结果落在可解释的区间里。import numpy as np def omni_normalize(raw, valid_mask): 按有效像元的分位数归一压掉极端值的影响 v raw[valid_mask np.isfinite(raw)] p99 np.percentile(v, 99) return np.clip(raw / p99, 0, 1) if p99 0 else raw用 99 分位而不是最大值是因为裸露地和窄谷地形会偶尔冒出极高的单点用最大值归一会把整张图压扁。第二是分级。归一化之后按分位数切档比如前 10% 划为核心廊道、10%–30% 为次级通道、其余为背景。分档线要在方法里写明因为它是结论的一部分。我还会做一步交叉检查把核心廊道区叠回阻力面看它有没有大量穿过高阻力核心区。穿过去说明参数太松廊道被“推”上了不该走的地方。5. 避坑与排查五类翻车现场和对应的检查动作前面几章讲的是怎么跑对这一章讲怎么跑错之后快速定位。下面五条都是按“现象→原因→解决”写的遇到问题直接对号入座。5.1 数据层面的三个坑现象跑完的电流图里廊道全落在本该避开的区域。原因阻力面方向反了。手头那份图层是适宜性值越大越好却按阻力值越大越难走喂了进去于是“好地方”变成了“最难穿的地方”电流全跑到边缘去了。 解决先确认图层语义再改habitat_map_is_resistances。如果源数据本来就是适宜性做一次单调递减变换取倒数或线性反转并且把这个变换写进方法。判断方法很简单找出已知的强阻力区比如建设用地、裸岩看它的值是不是全图最高的一档。现象电阻矩阵里大量 0 值或者干脆全是 NaN。原因节点坐标和栅格不在同一坐标系或者use_pixel_centers设反了导致多个节点落进同一个像元距离 0或者落到栅格范围外NaN。 解决把节点坐标和栅格的范围叠在一起看一眼五个节点画四个点就能发现。落进同一像元的节点要么合并成一个要么在像元尺度上重新布点。这一步我会在跑之前固定做一次比跑完之后回溯省事得多。现象电流全部绕着一个巨大的空白区走像有一堵墙。原因NoData 没统一。拼接图层里这个用 -9999、那个用 255软件只认声明的那个另一种就被当成真实数值读进去了可能是极高阻力也可能接近 0。 解决跑之前统一重投影并统一 NoData 值用一次栅格运算把所有非法值刷成同一个。同时把set_null_currents_to_nodata打开出图时才不会把背景和真电流混在一起。5.2 运行与输出层面的两个坑现象进程安安静静挂了几个小时日志一动不动最后被系统杀掉。原因pairwise 模式组合数爆炸——节点数一多计算量按平方涨或者栅格范围大且阻力值分布极端单次求解的迭代次数暴涨。CPU 并不是瓶颈内存才是。 解决三件事一起做。用use_included_pairs只算关心的节点对把组合数压下来把max_parallel降到 4 或 2先确认不是内存问题给 subprocess 加timeout让失败快速暴露而不是静默挂死。区域尺度就直接换 Omniscape不要硬扛。现象同一份配置今天跑和昨天跑结果不一样。原因ini 里用了相对路径两次运行的工作目录不同读到的其实是两份不同的输入文件。或者节点文件被改动过而没记录下来。 解决所有路径在生成 ini 时resolve()成绝对路径配置文件本身进版本库和输出目录一一对应。这条听起来像洁癖但它是唯一能让三个月后的你复现当时结论的办法。我现在连中间试验的 ini 都留着占不了多少空间。6. 结果自检拿一条只有唯一解法的栅格当标尺参数调对了没有不能靠图好看判断。我习惯在跑真实数据之前先拿一个人为构造、答案已知的小栅格做一次标定相当于给软件做体检。构造方法很直接造一张 100×100 的栅格背景阻力全填 100再从 (10,50) 到 (90,50) 横向拉一条一像元宽的低阻力通道值填 1。两个节点分别放在通道两端。因为四周全是高阻力电流没有别的路可走通道就是唯一路径。这时候总电阻应该约等于通道长度乘单元阻力——从第 10 列到第 90 列是 80 个连接理论值约 80。跑出来如果落在 78–82 这个量级说明阻力语义、连接方案、坐标解释三件事都是对的。如果差得离谱往下查的方向就明确了数值翻倍通常是八邻域在通道两侧通过高阻力像元形成了并联路径或者connect_using_avg_resistances的设置和你以为的不一样数值差一个数量级基本是阻力单位搞错了是每像元还是每米数值接近 0 则是两个节点被塞进了同一个像元。这个标定还有个附带好处它能帮你判断四邻域和八邻域在你的数据分辨率下差多少。同一条通道四邻域会给出接近 80 的干净结果八邻域略低一些差值就是对角并联路径的贡献。这个差值有多大直接决定你在正式数据上该选哪个方案——差值小就随大流用八邻域差值大就说明你的像元太粗对角捷径会实质性改变结论。三个定量自检指标跑完真实数据我都会过一遍。第一电阻矩阵必须对称对角为 0且满足三角不等式——电阻距离是一个度量任意三点之间R_ij ≤ R_ik R_kj必须成立不成立就是数值出问题了。第二电流守恒把所有源注入的电流加起来和汇流出的总量应当抵消累计电流图的全局量级应该随配对数线性增长如果不是多半有配对没跑完。第三把电流高值区和最小成本路径叠在一起前者应当覆盖或包住后者如果两者完全不重叠先别怀疑电路理论回头查坐标系。这套体检流程我走了很多次之后已经成了肌肉记忆新数据进来先跑标定栅格再看对称性最后才看真实结果。有一次赶进度跳过标定直接拿一张投影定义错误的阻力面跑了三天出来的廊道漂亮得不像话全部落在研究区外边——从那以后我每次换数据集都强制走一遍这条标定加三项检查的流程宁可多花二十分钟。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助?
咨询建站