本文介绍基于Python中ArcPy模块，实现Excel数据读取并导入图层，同时进行IDW插值与批量掩膜的方法。在这里，包含北京市各PM2.5浓度监测站点的矢量点要素图层是基于Python中ArcPy读取Excel表格数据创建矢量要素图层并生成属性表这篇文章得到的，如下图所示。而记录有北京市部分PM2.5浓度监测站点在2019年05月18日00时至23时等23个逐小时PM2.5浓度数据的Excel表格文件则如下所示，其中包括各站点在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中ArcPy读取Excel表格数据创建矢量要素图层并生成属性表这篇文章得到的，如下图所示。
其中，该矢量图层还包括属性表，属性表内容包括每一个站点的编号、地理位置与中文名称，如下图所示。
而记录有北京市部分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: ChuTianjiaimport copyimport arcpyimport pandas as pdfrom 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 tablefor i in range(len(field_list)):arcpy.AddField_management(shape_file_path,field_list[i],"SHORT")delete_num=0delete_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 dataif row[0] not in column_name_list:cursor.deleteRow()delete_num =1delete_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 interpolationarcpy.env.extent=boundary_file_pathfor 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 masktif_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")))
执行上述代码，如果是在ArcMap中直接通过工具箱运行，则可以看到代码运行过程中出现的提示。
例如，下图所示提示可以知道有哪几个站点是没有数据、从而被剔除的。
下图则可以显示出目前代码的运行情况。
同时，在我们设定的结果文件夹中可以看到，23小时的插值图与掩膜图都将自动生成并保存在指定文件夹。
再来看看具体的图片长什么样子。
首先查看IDW插值结果图；我们以当日10时的插值结果图为例进行查看。可以看到其已对北京市边界矢量数据所包含的矩形范围完成了插值。
接下来，查看IDW插值结果图经过掩膜后的图像；我们同样以当日10时的插值、掩膜结果图为例进行查看。可以看到，经过掩膜操作后的图像已经完全符合北京市边界矢量数据的范围。
