??本文介紹基于Python中ArcPy模塊，實(shí)現(xiàn)Excel數(shù)據(jù)讀取并導(dǎo)入圖層，同時(shí)進(jìn)行IDW插值與批量掩膜的方法。

1 任務(wù)需求

??首先，我們來明確一下本文所需實(shí)現(xiàn)的需求。

??現(xiàn)有一個(gè)記錄有北京市部分PM2.5濃度監(jiān)測(cè)站點(diǎn)在2019年05月18日00時(shí)至23時(shí)（其中不含19時(shí)）等23個(gè)逐小時(shí)PM2.5濃度數(shù)據(jù)的Excel表格文件，我們需要將其中的數(shù)據(jù)依次讀入一個(gè)包含北京市各PM2.5濃度監(jiān)測(cè)站點(diǎn)的矢量點(diǎn)要素圖層中；隨后，基于這些站點(diǎn)導(dǎo)入的23個(gè)逐小時(shí)PM2.5濃度數(shù)據(jù)，逐小時(shí)對(duì)北京市PM2.5濃度加以反距離加權(quán)（IDW）方法的插值，即共繪制23幅插值圖；最后，基于已有的北京市邊界矢量數(shù)據(jù)，分別對(duì)這23幅插值圖加以掩膜。

??在這里，包含北京市各PM2.5濃度監(jiān)測(cè)站點(diǎn)的矢量點(diǎn)要素圖層是基于這篇博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/122849279）得到的，如下圖所示。

??其中，該矢量圖層還包括屬性表，屬性表內(nèi)容包括每一個(gè)站點(diǎn)的編號(hào)、地理位置與中文名稱，如下圖所示。

??而記錄有北京市部分PM2.5濃度監(jiān)測(cè)站點(diǎn)在2019年05月18日00時(shí)至23時(shí)（其中不含19時(shí)）等23個(gè)逐小時(shí)PM2.5濃度數(shù)據(jù)的Excel表格文件則如下所示，其中包括各站點(diǎn)在23個(gè)整點(diǎn)時(shí)所監(jiān)測(cè)到的PM2.5濃度。

2 代碼實(shí)現(xiàn)

??了解了需求后，我們就基于Python中的ArcPy模塊，進(jìn)行詳細(xì)代碼的撰寫與介紹。

??這里需要說明的是：在編寫代碼的時(shí)候，為了方便執(zhí)行，所以希望代碼后期可以在ArcMap中直接通過工具箱運(yùn)行，即用到Python程序腳本新建工具箱與自定義工具的方法；因此，代碼中對(duì)于一些需要初始定義的變量，都用到了arcpy.GetParameterAsText()函數(shù)。大家如果只是希望在IDLE中運(yùn)行代碼，那么直接對(duì)這些變量進(jìn)行具體賦值即可。關(guān)于Python程序腳本新建工具箱與自定義工具，大家可以查看這篇博客（https://blog.csdn.net/zhebushibiaoshifu/article/details/121518404）詳細(xì)了解。

??上面提到需要初始定義的變量一共有十個(gè)，其中arcpy.env.workspace參數(shù)表示當(dāng)前工作空間，csv_path參數(shù)表示存儲(chǔ)有北京市2019年05月18日00時(shí)至23時(shí)（其中不含19時(shí)）逐小時(shí)PM2.5濃度數(shù)據(jù)的.csv文件，shape_file_path參數(shù)表示站點(diǎn)信息矢量數(shù)據(jù)文件，boundary_file_path參數(shù)表示投影后北京市邊界矢量數(shù)據(jù)文件，spatial_resolution參數(shù)表示IDW插值結(jié)果柵格圖的像元大小，power參數(shù)表示IDW插值時(shí)所用距離的冪指數(shù)，look_point參數(shù)表示IDW插值時(shí)所用最鄰近輸入采樣點(diǎn)數(shù)量的整數(shù)值，max_distance參數(shù)表示IDW插值時(shí)對(duì)最鄰近輸入采樣點(diǎn)的限制距離，單位依據(jù)地圖坐標(biāo)系確定；idw_result_dir參數(shù)表示IDW插值結(jié)果圖層保存路徑，mask_result_dir參數(shù)表示IDW插值結(jié)果圖層經(jīng)掩膜后保存路徑。

??代碼的整體思路為：首先利用pd.read_csv函數(shù)讀取記錄有北京市部分PM2.5濃度監(jiān)測(cè)站點(diǎn)在2019年05月18日00時(shí)至23時(shí)（其中不含19時(shí)）等23個(gè)逐小時(shí)PM2.5濃度數(shù)據(jù)的Excel表格文件數(shù)據(jù)，隨后在北京市各PM2.5濃度監(jiān)測(cè)站點(diǎn)的矢量點(diǎn)要素圖層的屬性表中新建23個(gè)列，每一個(gè)列表示該監(jiān)測(cè)站點(diǎn)在某一時(shí)刻的濃度數(shù)據(jù)（共有23個(gè)時(shí)刻，因此共有23個(gè)列）；其次，由于矢量要素圖層中的部分站點(diǎn)在Excel文件中并沒有數(shù)據(jù)，因此需要將這些站點(diǎn)從矢量要素圖層中刪除；最后，分別利用Idw函數(shù)與ExtractByMask函數(shù)進(jìn)行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 運(yùn)行結(jié)果

??執(zhí)行上述代碼，如果是在ArcMap中直接通過工具箱運(yùn)行，則可以看到代碼運(yùn)行過程中出現(xiàn)的提示。

??例如，下圖所示提示可以知道有哪幾個(gè)站點(diǎn)是沒有數(shù)據(jù)、從而被剔除的。

??下圖則可以顯示出目前代碼的運(yùn)行情況。

??同時(shí)，在我們?cè)O(shè)定的結(jié)果文件夾中可以看到，23小時(shí)的插值圖與掩膜圖都將自動(dòng)生成并保存在指定文件夾。

??再來看看具體的圖片長(zhǎng)什么樣子。

??首先查看IDW插值結(jié)果圖；我們以當(dāng)日10時(shí)的插值結(jié)果圖為例進(jìn)行查看?？梢钥吹狡湟褜?duì)北京市邊界矢量數(shù)據(jù)所包含的矩形范圍完成了插值。

??接下來，查看IDW插值結(jié)果圖經(jīng)過掩膜后的圖像；我們同樣以當(dāng)日10時(shí)的插值、掩膜結(jié)果圖為例進(jìn)行查看?？梢钥吹剑?jīng)過掩膜操作后的圖像已經(jīng)完全符合北京市邊界矢量數(shù)據(jù)的范圍。