本文介紹基於Python中ArcPy模塊，實現Excel數據讀取並生成矢量圖層，同時進行IDW插值與批量掩膜的方法。 1 任務需求 首先，我們來明確一下本文所需實現的需求。 現有一個記錄有北京市部分PM2.5濃度監測站點在2019年05月18日00時至23時（其中不含19時）等23個逐小時PM2.5 ...
本文介紹基於Python中ArcPy模塊，實現Excel數據讀取並生成矢量圖層，同時進行IDW插值與批量掩膜的方法。
1 任務需求
首先，我們來明確一下本文所需實現的需求。
現有一個記錄有北京市部分PM2.5濃度監測站點在2019年05月18日00時至23時（其中不含19時）等23個逐小時PM2.5濃度數據的Excel表格文件，我們需要將其中的數據依次讀入一個包含北京市各PM2.5濃度監測站點的矢量點要素圖層中；隨後，基於這些站點導入的23個逐小時PM2.5濃度數據，逐小時對北京市PM2.5濃度加以反距離加權（IDW）方法的插值，即共繪製23幅插值圖；最後，基於已有的北京市邊界矢量數據，分別對這23幅插值圖加以掩膜。
在這裡，包含北京市各PM2.5濃度監測站點的矢量點要素圖層是基於Python基於Excel生成矢量圖層及屬性表信息：ArcPy得到的，如下圖所示。
其中，該矢量圖層還包括屬性表，屬性表內容包括每一個站點的編號、地理位置與中文名稱，如下圖所示。
而記錄有北京市部分PM2.5濃度監測站點在2019年05月18日00時至23時（其中不含19時）等23個逐小時PM2.5濃度數據的Excel表格文件則如下所示，其中包括各站點在23個整點時所監測到的PM2.5濃度。
2 代碼實現
瞭解了需求後，我們就基於Python中的ArcPy模塊，進行詳細代碼的撰寫與介紹。
這裡需要說明的是：在編寫代碼的時候，為了方便執行，所以希望代碼後期可以在ArcMap中直接通過工具箱運行，即用到Python程式腳本新建工具箱與自定義工具的方法；因此，代碼中對於一些需要初始定義的變數，都用到了arcpy.GetParameterAsText()函數。大家如果只是希望在IDLE中運行代碼，那麼直接對這些變數進行具體賦值即可。關於Python程式腳本新建工具箱與自定義工具，大家可以查看ArcMap將Python寫的代碼轉為工具箱與自定義工具詳細瞭解。
上面提到需要初始定義的變數一共有十個，其中arcpy.env.workspace參數表示當前工作空間，csv_path參數表示存儲有北京市2019年05月18日00時至23時（其中不含19時）逐小時PM2.5濃度數據的.csv文件，shape_file_path參數表示站點信息矢量數據文件，boundary_file_path參數表示投影後北京市邊界矢量數據文件，spatial_resolution參數表示IDW插值結果柵格圖的像元大小，power參數表示IDW插值時所用距離的冪指數，look_point參數表示IDW插值時所用最鄰近輸入採樣點數量的整數值，max_distance參數表示IDW插值時對最鄰近輸入採樣點的限制距離，單位依據地圖坐標系確定；idw_result_dir參數表示IDW插值結果圖層保存路徑，mask_result_dir參數表示IDW插值結果圖層經掩膜後保存路徑。
代碼的整體思路為：首先利用pd.read_csv函數讀取記錄有北京市部分PM2.5濃度監測站點在2019年05月18日00時至23時（其中不含19時）等23個逐小時PM2.5濃度數據的Excel表格文件數據，隨後在北京市各PM2.5濃度監測站點的矢量點要素圖層的屬性表中新建23個列，每一個列表示該監測站點在某一時刻的濃度數據（共有23個時刻，因此共有23個列）；其次，由於矢量要素圖層中的部分站點在Excel文件中並沒有數據，因此需要將這些站點從矢量要素圖層中刪除；最後，分別利用Idw函數與ExtractByMask函數進行IDW插值與掩膜。
具體代碼如下。
# -*- coding: utf-8 -*-
# @author: ChuTianjia
import copy
import arcpy
import pandas as pd
from arcpy.sa import *
arcpy.env.workspace=arcpy.GetParameterAsText(0)
csv_path=arcpy.GetParameterAsText(1)
shape_file_path=arcpy.GetParameterAsText(2)
idw_result_dir=arcpy.GetParameterAsText(8)
boundary_file_path=arcpy.GetParameterAsText(3)
mask_result_dir=arcpy.GetParameterAsText(9)
spatial_resolution=arcpy.GetParameterAsText(4)
power=arcpy.GetParameterAsText(5)
look_point=arcpy.GetParameterAsText(6)
max_distance=arcpy.GetParameterAsText(7)
csv_data=pd.read_csv(csv_path,header=0,encoding="gbk")
column_name_list=list(csv_data)
hour_column=csv_data["hour"]
pm_25_list=[[0]*len(csv_data) for i in range(csv_data.shape[1]-3)]
for i in range(3,csv_data.shape[1]):
for index,data in csv_data.iterrows():
pm_25_list[i-3][index]=data[i]
field_list=["hour_00","hour_01","hour_02","hour_03","hour_04","hour_05",\
"hour_06","hour_07","hour_08","hour_09","hour_10",\
"hour_11","hour_12","hour_13","hour_14","hour_15",\
"hour_16","hour_17","hour_18","hour_20",\
"hour_21","hour_22","hour_23"]
field_list_use=copy.deepcopy(field_list)
field_list_use.insert(0,"Name")
# Update the columns in the attribute table
for i in range(len(field_list)):
arcpy.AddField_management(shape_file_path,field_list[i],"SHORT")
delete_num=0
delete_name=[]
with arcpy.da.UpdateCursor(shape_file_path,field_list_use) as cursor:
for row in cursor:
for column_name in column_name_list:
if column_name==row[0]:
for i in range(len(csv_data[column_name])):
row[i+1]=csv_data[column_name][i]
cursor.updateRow(row)
# Find stations that without any data
if row[0] not in column_name_list:
cursor.deleteRow()
delete_num+=1
delete_name.append(row[0])
arcpy.AddWarning("Delete {0} site(s) that do not contain any data, and the site(s) name is(are):".format(delete_num))
for i in delete_name:
arcpy.AddMessage(i)
arcpy.AddMessage("\n")
# Perform IDW interpolation
arcpy.env.extent=boundary_file_path
for i in range(len(field_list)):
idw_result=Idw(shape_file_path,field_list[i],spatial_resolution,power,RadiusVariable(look_point,max_distance))
idw_result_path=idw_result_dir+"\\"+"BJ_"+field_list[i]+".tif"
idw_result.save(idw_result_path)
arcpy.AddMessage("{0} has completed IDW interpolation.".format(field_list[i]))
# Perform mask
tif_file_list=arcpy.ListRasters("BJ_hour_*","TIF")
for raster in tif_file_list:
mask_result=ExtractByMask(raster,boundary_file_path)
mask_result_path=mask_result_dir+"\\"+raster.strip(".tif")+"_Mask.tif"
mask_result.save(mask_result_path)
arcpy.AddMessage("{0} has been masked.".format(raster.strip(".tif")))
3 運行結果
執行上述代碼，如果是在ArcMap中直接通過工具箱運行，則可以看到代碼運行過程中出現的提示。
例如，下圖所示提示可以知道有哪幾個站點是沒有數據、從而被剔除的。
下圖則可以顯示出目前代碼的運行情況。
同時，在我們設定的結果文件夾中可以看到，23小時的插值圖與掩膜圖都將自動生成並保存在指定文件夾。
再來看看具體的圖片長什麼樣子。
首先查看IDW插值結果圖；我們以當日10時的插值結果圖為例進行查看。可以看到其已對北京市邊界矢量數據所包含的矩形範圍完成了插值。
接下來，查看IDW插值結果圖經過掩膜後的圖像；我們同樣以當日10時的插值、掩膜結果圖為例進行查看。可以看到，經過掩膜操作後的圖像已經完全符合北京市邊界矢量數據的範圍。
至此，大功告成。
