1 任务需求

在这里，包含北京市 各PM2.5浓度监测站点 的矢量点要素图层是基于 这篇博客 （https://blog.csdn.net/zhebushibiaoshifu/article/details/122849279）得到的，如下图所示。

其中，该矢量图层还包括属性表，属性表内容包括每一个站点的编号、地理位置与中文名称，如下图所示。

而记录有北京市 部分PM2.5浓度监测站点 在2019年05月18日00时至23时（其中不含19时）等 23个逐小时PM2.5浓度数据 的 Excel 表格文件则如下所示，其中包括各站点在23个整点时所监测到的 PM2.5 浓度。

2 代码实现

这里需要说明的是：在编写代码的时候，为了方便执行，所以希望代码后期可以在 ArcMap 中直接通过工具箱运行，即用到 Python程序脚本新建工具箱与自定义工具 的方法；因此，代码中对于一些需要初始定义的变量，都用到了 arcpy.GetParameterAsText() 函数。大家如果只是希望在 IDLE 中运行代码，那么直接对这些变量进行具体赋值即可。关于 Python程序脚本新建工具箱与自定义工具 ，大家可以查看 这篇博客 （https://blog.csdn.net/zhebushibiaoshifu/article/details/121518404）详细了解。

上面提到需要初始定义的变量一共有十个，其中 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时的插值、掩膜结果图为例进行查看。可以看到，经过掩膜操作后的图像已经完全符合北京市边界矢量数据的范围。
 
欢迎关注公众号/CSDN/知乎/微博：疯狂学习GIS
 
				
注：   全新上架。全新录制，整合3年来学员的常见问题。课程2020年全新上架。将近90课时。内容丰富，可涵盖所有常见工作的问题。同时保留老版本课程，详情请看课程目录及以下课程概述。     课程目录第一章第三节（课程资料（PPt与操作数据等）），是整套课程的课件PPT，资料、课程操作数据等的下载地址     本课程经过全面的再录制，更全面、更系统化，支持随到随学，免费试学。利用ArcGIS10.6中文版教学，试用于ArcGIS10.0、10.1、10.2、10.3、10.4、10.5、10.6、10.7、10.8系列，让零基础或者掌握不全面的人快速系统地了解ArcGIS的应用，让学习者对ArcGIS整体认识、空间数据信息采集、属性表操作、拓扑、空间数据可视化、出图、数据更新、投影变换与格式转换、矢量、栅格数据空间分析有一个全新的认识.  
基于Python空间数据处理课程的Arcpy大作业，分享一些经历和心得。
题目要求：
（1）在“地质调查点基础数据表.xls”中图幅范围内增加100个随机位置的高程点。构建一个shape文件，采用自定义工具的模式，参数有两个：一个是让用户选择excel文件，一个让用户指定新生成的文件名。
（2）采用spline方法插值shape点文件中的高程信息。
（3）将高程图和点图层叠加显示。制作简单的制图模板，将图件按照图幅编号分幅输出为jpg图。
图幅范围计算
要生成随机点，首先应确定随机点的空间范围。
				
本文我们将介绍IDW(反距离加权法（Inverse Distance Weighted）)插值的Python计算方法及插值结果的可视化绘制过程。主要涉及的知识点如下：
	IDW简介
	自定义Python代码计算空间IDW
	分别使用plotnine、Basemap进行IDW插值结果可视化绘制
IDW简介
反距离权重 (IDW)插值假设：彼此距离较近的事物要比彼此距离较远的事物更相似。当为任何未测量的位置预测值时，反距离权重法会采用预测位置周围的测量值与距离预测位置较远的测量值相..
				
前言： 研一，刚开始学python，初衷其实是想以arcgis软件上没有的插值方法进行空间插值，奈何人菜瘾大，只能先从基础的插值方法开始，这次的IDW插值算是一个开始吧，我会慢慢记录自己的成长
				
ArcPy是ArcGIS中的一个Python模块，用于执行GIS数据处理和自动化任务。如果要使用Python写ArcPy分割影像，您可以使用以下代码：
import arcpy
# Set the workspace
arcpy.env.workspace = r"C:\ArcPy\Data"
# Set the input raster
inRaster = "input_raster"
# Set the output feature class
outFeatureClass = "output_feature_class"
# Execute the Feature Class to Raster tool
arcpy.RasterToPolygon_conversion(inRaster, outFeatureClass)
上面的代码首先设置ArcPy的工作空间，然后将输入影像作为变量`inRaster`，设置输出要素类作为变量`outFeatureClass`。最后，使用`RasterToPolygon_conversion`函数执行影像到多边形转换，从而完成分割影像的操作。