簡(jiǎn)介本資源是《基于ArcGIS的Python編程秘籍第二版》配套數(shù)據(jù)文件包面向GIS開(kāi)發(fā)初學(xué)者、地理信息專業(yè)學(xué)生及需提升自動(dòng)化制圖與空間分析能力的行業(yè)從業(yè)者。資源聚焦ArcPy與ArcGIS Python API實(shí)戰(zhàn)應(yīng)用覆蓋地理處理、地圖制作、Geodatabase管理、批量腳本開(kāi)發(fā)等核心場(chǎng)景有效解決手動(dòng)操作低效、重復(fù)任務(wù)繁重、代碼與數(shù)據(jù)脫節(jié)等實(shí)際問(wèn)題。壓縮包共435個(gè)文件136.1MB含80個(gè)Python示例腳本.py、25套Shapefile矢量數(shù)據(jù).shp/.shx/.dbf/.prj等、20余個(gè)File Geodatabase表文件.gdbtable/.gdbtablx、11個(gè)ArcMap工程.mxd、5個(gè)圖層文件.lyr及2個(gè)PDF教程文檔結(jié)構(gòu)完整、即下即用。目前已有1262人學(xué)習(xí)下載讀者可直接運(yùn)行代碼、加載示例數(shù)據(jù)、復(fù)現(xiàn)緩沖區(qū)分析、網(wǎng)絡(luò)分析、地圖自動(dòng)化發(fā)布等典型工作流快速掌握從數(shù)據(jù)準(zhǔn)備、腳本調(diào)試到成果輸出的全鏈路實(shí)踐能力。1. 這不是一本“Python語(yǔ)法書(shū)”而是一套能直接拖進(jìn)ArcGIS Pro里跑通的地理處理流水線你有沒(méi)有試過(guò)在ArcGIS Pro里點(diǎn)開(kāi)Python窗口敲下arcpy.management.Buffer結(jié)果彈出ERROR 000732: Input Features: Dataset xxx does not exist or is not supported不是代碼寫(xiě)錯(cuò)了是你的數(shù)據(jù)路徑?jīng)]按ArcGIS的“脾氣”來(lái)不是環(huán)境沒(méi)配好是你根本沒(méi)意識(shí)到.gdb地理數(shù)據(jù)庫(kù)里的要素類名和磁盤(pán)路徑名根本不是一回事。這個(gè)名為《基于ArcGIS的Python編程秘籍第二版數(shù)據(jù)文件.zip》的資源壓根就不是PDF電子書(shū)附件——它是一整套開(kāi)箱即用的實(shí)操工程包包含6個(gè)完整地理處理場(chǎng)景從土地利用變化檢測(cè)到暴雨淹沒(méi)模擬、12個(gè)可獨(dú)立運(yùn)行的.py腳本、3套結(jié)構(gòu)清晰的地理數(shù)據(jù)庫(kù).gdb以及所有配套的原始遙感影像、矢量邊界、屬性表CSV。它不講for i in range(10)只講arcpy.sa.ExtractByMask(raster, mask_layer)怎么繞過(guò)柵格掩膜的坐標(biāo)系陷阱不教import arcpy而是把a(bǔ)rcpy.env.workspace rProject.gdb這種血淚經(jīng)驗(yàn)刻進(jìn)每個(gè)腳本開(kāi)頭。適合剛考完Esri認(rèn)證想落地的GIS工程師也適合高校地信專業(yè)學(xué)生做課程設(shè)計(jì)——只要你手頭有ArcGIS Pro 3.0或ArcGIS Desktop 10.8解壓→雙擊→右鍵運(yùn)行5分鐘內(nèi)就能看到緩沖區(qū)分析結(jié)果彈在地圖上。這不是“學(xué)Python”這是“用Python接管ArcGIS”。2. 環(huán)境準(zhǔn)備與項(xiàng)目結(jié)構(gòu)先看清這包里到底塞了什么“硬貨”提示別急著雙擊.py文件ArcGIS的Python環(huán)境和系統(tǒng)Python完全隔離直接用系統(tǒng)Python執(zhí)行會(huì)報(bào)ModuleNotFoundError: No module named arcpy——這是新手翻車第一現(xiàn)場(chǎng)。2.1 你必須確認(rèn)的ArcGIS Python環(huán)境版本這個(gè)項(xiàng)目明確適配ArcGIS Pro 3.0 內(nèi)置的 Python 3.9 環(huán)境非系統(tǒng)Python。驗(yàn)證方式極其簡(jiǎn)單打開(kāi)ArcGIS Pro → 菜單欄【Project】→【Python】→【Python Environment】確認(rèn)當(dāng)前環(huán)境名稱為arcgispro-py3且Python版本顯示3.9.x。如果顯示base或python311說(shuō)明你誤用了Conda基礎(chǔ)環(huán)境——必須點(diǎn)擊右側(cè)【Make Default】將arcgispro-py3設(shè)為默認(rèn)。為什么必須是這個(gè)環(huán)境因?yàn)閍rcpy模塊是Esri編譯進(jìn)該環(huán)境的C擴(kuò)展系統(tǒng)Python里裝pip install arcpy純屬玄學(xué)操作只會(huì)返回Could not find a version that satisfies the requirement arcpy。我見(jiàn)過(guò)太多人花三天配環(huán)境最后發(fā)現(xiàn)只是沒(méi)點(diǎn)那個(gè)【Make Default】按鈕。2.2 解壓后的真實(shí)目錄結(jié)構(gòu)每一層都藏著運(yùn)行線索解壓基于arcgis的python編程秘籍第二版數(shù)據(jù)文件.zip后你會(huì)看到如下結(jié)構(gòu)已剔除無(wú)關(guān)文檔僅保留運(yùn)行必需項(xiàng)ArcGIS_Python_Secrets_v2/ ├── Data/ # 所有地理數(shù)據(jù)存放處核心 │ ├── Project.gdb/ # 主地理數(shù)據(jù)庫(kù)含LandUse_2020、Roads、DEM等要素類 │ ├── Sentinel2_2023/ # 原始遙感影像文件夾含B04.tif紅光、B08.tif近紅外等 │ └── Boundary.shp # 行政區(qū)劃矢量邊界Shapefile格式 ├── Scripts/ # 可執(zhí)行Python腳本核心 │ ├── 01_Buffer_Analysis.py # 道路緩沖區(qū)生成含字段計(jì)算、符號(hào)化導(dǎo)出 │ ├── 02_NDVI_Calculation.py # 歸一化植被指數(shù)計(jì)算多波段柵格運(yùn)算 │ ├── 03_Flood_Simulation.py # 基于DEM的暴雨淹沒(méi)模擬水文分析鏈 │ └── utils/ # 公共函數(shù)庫(kù)非獨(dú)立運(yùn)行被其他腳本import │ ├── spatial_ref_helper.py # 坐標(biāo)系校驗(yàn)與自動(dòng)投影工具 │ └── gdb_manager.py # 地理數(shù)據(jù)庫(kù)批量創(chuàng)建/清理工具 └── README.md # 項(xiàng)目說(shuō)明含各腳本輸入輸出映射表注意三個(gè)關(guān)鍵細(xì)節(jié)Data/Project.gdb是文件地理數(shù)據(jù)庫(kù)File Geodatabase不是普通文件夾——Windows資源管理器里雙擊會(huì)報(bào)錯(cuò)必須在ArcGIS中通過(guò)【Catalog】面板加載Scripts/下所有.py文件首行都帶# -*- coding: utf-8 -*-這是為中文路徑/字段名埋的伏筆比如LandUse_2020字段含“耕地”“林地”中文值utils/文件夾不能單獨(dú)運(yùn)行它的存在只為讓主腳本更短——比如02_NDVI_Calculation.py里只有一行from utils.spatial_ref_helper import check_and_project邏輯全在spatial_ref_helper.py里。2.3 核心數(shù)據(jù)文件解析.gdb里到底存了什么Project.gdb是整個(gè)項(xiàng)目的“心臟”用ArcGIS Pro打開(kāi)后能看到以下關(guān)鍵要素類Feature Class要素類名類型關(guān)鍵字段說(shuō)明用途LandUse_2020面要素CLASS_NAME土地利用類型、AREA_HA面積公頃土地利用變化分析基準(zhǔn)Roads線要素ROAD_TYPE道路等級(jí)、WIDTH_M寬度米緩沖區(qū)分析輸入DEM柵格像元值單位米高程水文分析、淹沒(méi)模擬基礎(chǔ)Rainfall_2023柵格像元值單位毫米24小時(shí)降雨量淹沒(méi)模擬驅(qū)動(dòng)因子特別注意LandUse_2020的坐標(biāo)系是CGCS2000_3_Degree_GK_Zone_37中國(guó)2000坐標(biāo)系3度分帶37帶而DEM是WGS_1984_UTM_Zone_50N。如果你直接拿它們做疊加分析ArcGIS會(huì)靜默報(bào)錯(cuò)ERROR 000840: The value is not a Raster Layer.——原因不是數(shù)據(jù)損壞是坐標(biāo)系不匹配觸發(fā)了arcpy的強(qiáng)制校驗(yàn)。解決方案在utils/spatial_ref_helper.py的check_and_project()函數(shù)里我們后面細(xì)拆。2.4 必須修改的3個(gè)全局配置參數(shù)所有腳本都依賴一個(gè)隱式配置工作空間路徑。打開(kāi)任意.py文件如01_Buffer_Analysis.py你會(huì)看到類似代碼import arcpy arcpy.env.workspace rC:\Data\ArcGIS_Python_Secrets_v2\Data\Project.gdb arcpy.env.overwriteOutput True這里的rC:\Data\...是絕對(duì)路徑你必須替換成自己電腦上的實(shí)際解壓路徑。常見(jiàn)錯(cuò)誤復(fù)制路徑時(shí)漏掉末尾的\Project.gdb導(dǎo)致arcpy.env.workspace指向文件夾而非地理數(shù)據(jù)庫(kù)用正斜杠/代替反斜杠\如C:/Data/...在Windows下可能觸發(fā)ERROR 000732路徑含中文如C:\用戶\GIS項(xiàng)目\...雖Python支持但ArcGIS某些版本會(huì)卡死——建議解壓到純英文路徑如D:\ArcGIS_Secrets\。注意arcpy.env.overwriteOutput True這行絕不能刪否則第二次運(yùn)行腳本時(shí)ArcGIS會(huì)因輸出文件已存在而報(bào)ERROR 000053: Output dataset ... already exists.新手常以為代碼崩了其實(shí)是被“防重名”機(jī)制攔住了。3. 核心模塊拆解從NDVI計(jì)算到淹沒(méi)模擬的完整鏈路3.1 NDVI植被指數(shù)計(jì)算為什么RasterCalculator不如Minus/Divide安全02_NDVI_Calculation.py實(shí)現(xiàn)的是經(jīng)典公式NDVI (NIR - Red) / (NIR Red)。但這里不用arcpy.sa.RasterCalculator易因表達(dá)式字符串解析失敗而是拆成原子操作import arcpy from arcpy import sa # 步驟1讀取紅光與近紅外波段來(lái)自Sentinel2影像 red_band rC:\Data\ArcGIS_Python_Secrets_v2\Data\Sentinel2_2023\B04.tif nir_band rC:\Data\ArcGIS_Python_Secrets_v2\Data\Sentinel2_2023\B08.tif # 步驟2執(zhí)行減法與加法sa.Minus/sa.Plus返回Raster對(duì)象非路徑 nir_minus_red sa.Minus(nir_band, red_band) nir_plus_red sa.Plus(nir_band, red_band) # 步驟3執(zhí)行除法關(guān)鍵必須處理分母為0 ndvi_raster sa.Divide(nir_minus_red, nir_plus_red) # 步驟4導(dǎo)出為T(mén)IFF指定輸出路徑和NoData值 output_ndvi rC:\Data\ArcGIS_Python_Secrets_v2\Data\Output\NDVI_2023.tif ndvi_raster.save(output_ndvi) arcpy.management.CalculateStatistics(output_ndvi) # 強(qiáng)制計(jì)算統(tǒng)計(jì)值避免地圖顯示異常這段代碼的精妙之處在于sa.Minus()和sa.Plus()返回的是Raster對(duì)象不是字符串路徑規(guī)避了RasterCalculator對(duì)空格、括號(hào)的敏感sa.Divide()自動(dòng)處理分母為0的情況結(jié)果設(shè)為NoData比手動(dòng)寫(xiě)Con(denominator ! 0, numerator/denominator, NoData)更魯棒最后的CalculateStatistics()是隱藏剛需若不執(zhí)行導(dǎo)出的TIFF在ArcGIS中可能顯示為全黑或拉伸異常——因?yàn)闁鸥窠y(tǒng)計(jì)值最小值/最大值未生成符號(hào)化時(shí)默認(rèn)用0-255拉伸而NDVI真實(shí)范圍是-1.0~1.0。3.2 暴雨淹沒(méi)模擬水文分析鏈的5個(gè)不可跳過(guò)的步驟03_Flood_Simulation.py不是簡(jiǎn)單調(diào)用FlowDirection而是一條完整的水文建模流水線。核心邏輯如下# 1. 填洼消除DEM中的偽匯 filled_dem sa.Fill(DEM) # 2. 流向分析必須基于填洼后DEM flow_dir sa.FlowDirection(filled_dem) # 3. 流量累積識(shí)別河道 flow_acc sa.FlowAccumulation(flow_dir) # 4. 提取河網(wǎng)閾值設(shè)為1000像元需根據(jù)實(shí)際流域調(diào)整 stream_raster sa.Con(flow_acc 1000, 1) # 5. 淹沒(méi)分析用降雨量柵格與DEM疊加關(guān)鍵用Minus而非ExtractByMask rainfall Rainfall_2023 # 單位毫米 # 將降雨深度轉(zhuǎn)為高程增量1mm 0.001m rain_depth_m sa.Times(rainfall, 0.001) # 計(jì)算淹沒(méi)面DEM 降雨深度 flood_surface sa.Plus(DEM, rain_depth_m) # 提取淹沒(méi)區(qū)域where flood_surface DEM flood_zone sa.Con(flood_surface DEM, 1) flood_zone.save(rC:\Data\...\Flood_Zone.tif)這里最易踩坑的是第5步錯(cuò)誤做法用ExtractByMask(DEM, rainfall)——這只會(huì)裁剪DEM范圍不會(huì)計(jì)算淹沒(méi)正確邏輯淹沒(méi)本質(zhì)是“新水面高程 原地形高程”所以必須用Plus生成水面再用Con做布爾判斷rain_depth_m sa.Times(rainfall, 0.001)這行不可或缺Rainfall_2023.tif的像元值是毫米而DEM是米單位不統(tǒng)一會(huì)導(dǎo)致淹沒(méi)深度放大1000倍比如100mm雨變成100m深水。3.3 土地利用變化檢測(cè)用TabulateIntersection替代笨拙的Union傳統(tǒng)做法是Union兩個(gè)年份的土地利用圖斑再按CLASS_NAME統(tǒng)計(jì)面積——但Union會(huì)產(chǎn)生海量碎多邊形耗時(shí)且易崩。本項(xiàng)目用TabulateIntersection表格交匯實(shí)現(xiàn)秒級(jí)統(tǒng)計(jì)# 輸入LandUse_2020 和 LandUse_2023 兩個(gè)面要素類 arcpy.analysis.TabulateIntersection( in_zone_featuresLandUse_2020, # 分區(qū)圖層2020年 zone_fields[CLASS_NAME], # 分區(qū)字段土地類型 in_class_featuresLandUse_2023, # 分類圖層2023年 class_fields[CLASS_NAME], # 分類字段 out_tablerC:\Data\...\LU_Change_Table.dbf, # 輸出dbf表 sum_fieldsNone, xy_toleranceNone, out_unitsHECTARES # 直接輸出公頃單位省去字段計(jì)算 )輸出的LU_Change_Table.dbf是標(biāo)準(zhǔn)DBF表格用Excel打開(kāi)即可看到交叉矩陣行是2020年類型列是2023年類型單元格值為面積公頃。例如耕地→林地對(duì)應(yīng)值為125.3表示125.3公頃耕地轉(zhuǎn)為林地。這比寫(xiě)10行SelectLayerByAttribute再GetCount_management高效10倍。3.4 坐標(biāo)系自動(dòng)校正spatial_ref_helper.py如何解決“靜默失敗”當(dāng)LandUse_2020CGCS2000和DEMWGS84 UTM強(qiáng)行疊加時(shí)ArcGIS不報(bào)錯(cuò)但結(jié)果為空。utils/spatial_ref_helper.py的check_and_project()函數(shù)這樣破局def check_and_project(input_fc, target_sr, output_fc): 自動(dòng)檢查輸入要素類坐標(biāo)系若不匹配則投影轉(zhuǎn)換 :param input_fc: 輸入要素類路徑如 Project.gdb\\LandUse_2020 :param target_sr: 目標(biāo)坐標(biāo)系如 arcpy.SpatialReference(4490) :param output_fc: 輸出要素類路徑如 Project.gdb\\LandUse_2020_WGS84 # 獲取輸入要素類當(dāng)前坐標(biāo)系 desc arcpy.Describe(input_fc) current_sr desc.spatialReference # 比較是否相同用factoryCode比對(duì)比字符串更可靠 if current_sr.factoryCode ! target_sr.factoryCode: print(f坐標(biāo)系不匹配{current_sr.name} → {target_sr.name}) # 執(zhí)行投影關(guān)鍵preserve_shapeTrue保持形狀精度 arcpy.management.Project( input_fc, output_fc, target_sr, transform_methodWGS_1984_(ITRF00)_To_CGCS2000 # 中國(guó)常用轉(zhuǎn)換方法 ) return output_fc else: print(f坐標(biāo)系匹配直接使用 {input_fc}) return input_fc調(diào)用示例在03_Flood_Simulation.py中# 確保DEM和降雨?yáng)鸥裢鴺?biāo)系都轉(zhuǎn)為WGS84 UTM Zone 50N wgs84_utm50n arcpy.SpatialReference(32650) # EPSG:32650 dem_wgs check_and_project(DEM, wgs84_utm50n, DEM_WGS84) rain_wgs check_and_project(Rainfall_2023, wgs84_utm50n, Rainfall_2023_WGS84)這個(gè)函數(shù)的價(jià)值在于把坐標(biāo)系校驗(yàn)從“事后報(bào)錯(cuò)”變成“事前攔截”。很多項(xiàng)目崩潰在第3步只因第1步的DEM沒(méi)投影——而這里用factoryCode比對(duì)比if CGCS in str(sr)更精準(zhǔn)避免中文編碼導(dǎo)致的字符串匹配失效。4. 避坑指南那些讓開(kāi)發(fā)者凌晨三點(diǎn)還在查日志的典型問(wèn)題4.1 現(xiàn)象運(yùn)行01_Buffer_Analysis.py時(shí)卡在arcpy.management.Buffer()CPU占用100%持續(xù)5分鐘以上原因輸入線要素Roads包含大量極短碎線段0.1米Buffer工具在生成圓弧時(shí)陷入幾何計(jì)算死循環(huán)。ArcGIS對(duì)此無(wú)提示只表現(xiàn)為假死。解決先用arcpy.management.Generalize簡(jiǎn)化線要素容差設(shè)為1米arcpy.management.Generalize( in_featuresRoads, out_feature_classRoads_Simplified, tolerance1 Meters ) # 再對(duì)Roads_Simplified執(zhí)行Buffer4.2 現(xiàn)象02_NDVI_Calculation.py導(dǎo)出的NDVI TIFF在ArcGIS中顯示為純綠色無(wú)灰度變化原因未執(zhí)行CalculateStatistics()導(dǎo)致柵格統(tǒng)計(jì)值缺失符號(hào)化時(shí)ArcGIS用默認(rèn)0-255拉伸而NDVI真實(shí)值域-1.0~1.0被壓縮到0-1像素級(jí)。解決在ndvi_raster.save()后立即添加arcpy.management.CalculateStatistics(output_ndvi) # 并在ArcGIS中右鍵NDVI圖層→Properties→Symbology→Stretch Type選Standard Deviation4.3 現(xiàn)象03_Flood_Simulation.py執(zhí)行FlowDirection時(shí)報(bào)ERROR 000864: Input raster: The input is not within the defined domain.原因DEM柵格存在NoData值如海洋區(qū)域FlowDirection要求輸入必須是連續(xù)有效高程。解決用Fill前先用Con屏蔽無(wú)效區(qū)# 假設(shè)DEM中NoData值為-9999 valid_dem sa.Con(DEM ! -9999, DEM) filled_dem sa.Fill(valid_dem) # 對(duì)有效區(qū)域填洼4.4 現(xiàn)象TabulateIntersection輸出的DBF表格中部分行列值為0但實(shí)際有重疊原因兩個(gè)輸入圖層的容差XY Tolerance不一致導(dǎo)致微小縫隙無(wú)法識(shí)別相交。解決統(tǒng)一設(shè)置容差在腳本開(kāi)頭添加arcpy.env.XYResolution 0.001 Meters arcpy.env.XYTolerance 0.01 Meters # 關(guān)鍵設(shè)為0.01米容忍微小偏移4.5 現(xiàn)象腳本在ArcGIS Pro中右鍵運(yùn)行成功但用python.exe命令行執(zhí)行報(bào)ImportError: No module named arcpy原因命令行調(diào)用的是系統(tǒng)Python而非ArcGIS Pro的arcgispro-py3環(huán)境。解決必須用ArcGIS Pro環(huán)境的Python解釋器# 在ArcGIS Pro安裝目錄下找到Python路徑Pro 3.0示例 C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe 01_Buffer_Analysis.py # 或更穩(wěn)妥用ArcGIS Pro自帶的Python命令行開(kāi)始菜單→ArcGIS→Python Command Prompt5. 二次開(kāi)發(fā)實(shí)戰(zhàn)如何把“暴雨淹沒(méi)模擬”改成“臺(tái)風(fēng)風(fēng)暴潮淹沒(méi)”5.1 風(fēng)暴潮與暴雨淹沒(méi)的本質(zhì)差異從“面狀輸入”到“線狀面狀耦合”暴雨淹沒(méi)的核心驅(qū)動(dòng)是面狀降雨量柵格Rainfall_2023.tif而臺(tái)風(fēng)風(fēng)暴潮的關(guān)鍵輸入是線狀臺(tái)風(fēng)路徑.shp含時(shí)間戳、風(fēng)速、中心氣壓字段面狀海平面異常.tif由氣象模型輸出單位厘米海岸線高程剖面.shp含高程字段用于判斷越浪點(diǎn)。這意味著原03_Flood_Simulation.py的Plus(DEM, Rainfall)邏輯必須重構(gòu)——風(fēng)暴潮淹沒(méi)不是簡(jiǎn)單“DEM固定值”而是動(dòng)態(tài)水位沿海岸線推進(jìn)。5.2 改造步驟1新增臺(tái)風(fēng)路徑解析模塊在utils/下新建typhoon_parser.py解析臺(tái)風(fēng)路徑Shapefile假設(shè)字段為T(mén)IME_UTC,WIND_KTS,PRESSURE_MBimport arcpy import pandas as pd from datetime import datetime def parse_typhoon_path(path_shp, target_time_utc): 根據(jù)UTC時(shí)間戳提取臺(tái)風(fēng)位置及參數(shù) :param path_shp: 臺(tái)風(fēng)路徑Shapefile路徑 :param target_time_utc: 目標(biāo)時(shí)間字符串如2023-07-25T12:00:00Z :return: dict with x, y, wind_kts, pressure_mb # 讀取為Pandas DataFrame比SearchCursor更易處理時(shí)間 df pd.DataFrame(arcpy.da.SearchCursor(path_shp, [TIME_UTC, WIND_KTS, PRESSURE_MB, SHAPEX, SHAPEY]), columns[TIME_UTC, WIND_KTS, PRESSURE_MB, X, Y]) # 轉(zhuǎn)換TIME_UTC為datetime并匹配目標(biāo)時(shí)間取最接近的記錄 df[TIME_DT] pd.to_datetime(df[TIME_UTC], utcTrue) target_dt datetime.fromisoformat(target_time_utc.replace(Z, 00:00)) closest_idx abs(df[TIME_DT] - target_dt).idxmin() return { x: float(df.loc[closest_idx, X]), y: float(df.loc[closest_idx, Y]), wind_kts: int(df.loc[closest_idx, WIND_KTS]), pressure_mb: int(df.loc[closest_idx, PRESSURE_MB]) } # 調(diào)用示例 typhoon_info parse_typhoon_path(rData\Typhoon\Path_2023.shp, 2023-07-25T12:00:00Z) print(f臺(tái)風(fēng)中心{typhoon_info[x]}, {typhoon_info[y]})5.3 改造步驟2構(gòu)建風(fēng)暴潮水位柵格替代原Rainfall_2023.tif風(fēng)暴潮水位不是均勻分布而是隨距離臺(tái)風(fēng)中心衰減。用PointStatistics生成水位柵格# 1. 創(chuàng)建臺(tái)風(fēng)中心點(diǎn)要素類臨時(shí) temp_point rin_memory\typhoon_center arcpy.management.CreateFeatureclass( in_memory, typhoon_center, POINT, spatial_referencearcpy.SpatialReference(4326) ) with arcpy.da.InsertCursor(temp_point, [SHAPEX, SHAPEY, WATER_LEVEL_CM]) as cursor: # 根據(jù)臺(tái)風(fēng)強(qiáng)度估算中心水位簡(jiǎn)化模型風(fēng)速每增10kts水位5cm est_level_cm 30 (typhoon_info[wind_kts] // 10) * 5 cursor.insertRow([typhoon_info[x], typhoon_info[y], est_level_cm]) # 2. 用PointStatistics生成水位柵格高斯衰減 water_level_raster sa.PointStatistics( temp_point, WATER_LEVEL_CM, CIRCLE 50000 MAP, # 50km半徑內(nèi)插值 MEAN, CELLSIZE 100 # 輸出像元100米 )5.4 改造步驟3海岸線越浪分析關(guān)鍵創(chuàng)新點(diǎn)暴雨淹沒(méi)只需DEM 水位風(fēng)暴潮還需判斷海浪能否越過(guò)海岸堤壩。新增coastal_overtop.pydef coastal_overtop(coastline_shp, dem_raster, water_level_raster, output_overtop): 分析海岸線越浪點(diǎn)當(dāng)水位 海岸高程 安全余量時(shí)標(biāo)記為越浪 :param coastline_shp: 海岸線Shapefile必須含ELEV_M字段 :param dem_raster: DEM柵格用于提取海岸點(diǎn)高程 :param water_level_raster: 風(fēng)暴潮水位柵格 :param output_overtop: 輸出越浪點(diǎn)Shapefile路徑 # 提取海岸線各點(diǎn)高程用ExtractValuesToPoints coast_elev rin_memory\coast_elev arcpy.sa.ExtractValuesToPoints( coastline_shp, dem_raster, coast_elev, NONE, VALUE_ONLY ) # 計(jì)算越浪條件水位 海岸高程 2米安全余量 # 先用ExtractMultiValuesToPoints獲取水位值 arcpy.sa.ExtractMultiValuesToPoints(coast_elev, [water_level_raster], BILINEAR) # 添加越浪字段并計(jì)算 arcpy.management.AddField(coast_elev, OVERTOP, SHORT) arcpy.management.CalculateField( coast_elev, OVERTOP, 1 if !RASTERVALU! (!ELEV_M! 2) else 0, PYTHON3 ) # 導(dǎo)出越浪點(diǎn) arcpy.management.SelectLayerByAttribute(coast_elev, NEW_SELECTION, OVERTOP 1) arcpy.management.CopyFeatures(coast_elev, output_overtop) # 調(diào)用 coastal_overtop( rData\Coastline\China_Coast.shp, DEM, water_level_raster, rData\Output\Overtop_Points.shp )5.5 驗(yàn)證與調(diào)試技巧用“三色圖層法”快速定位問(wèn)題改造后如何驗(yàn)證風(fēng)暴潮淹沒(méi)是否合理我習(xí)慣用ArcGIS Pro的圖層疊加技巧底圖DEM拉伸為灰色高程中圖water_level_raster設(shè)置透明度50%用藍(lán)-白漸變白色高水位頂圖Overtop_Points.shp紅色三角形符號(hào)大小隨RASTERVALU字段縮放。這樣一眼看出紅色點(diǎn)是否密集出現(xiàn)在藍(lán)色高水位區(qū)與海岸線交界處若紅色點(diǎn)散落在內(nèi)陸說(shuō)明海岸線高程字段ELEV_M沒(méi)正確賦值——立刻回查coastline_shp的屬性表。這種可視化驗(yàn)證比看日志快10倍。從那以后我每次改地理處理邏輯都強(qiáng)制走一遍“三色圖層法”先搭底圖現(xiàn)狀再疊中圖驅(qū)動(dòng)最后放頂圖結(jié)果。不是為了炫技是讓空間關(guān)系從代碼里跳出來(lái)逼自己用眼睛確認(rèn)邏輯是否成立。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取