经验首页 前端设计 程序设计 Java相关 移动开发 数据库/运维 软件/图像 大数据/云计算 其他经验
当前位置:技术经验 » 程序设计 » Python » 查看文章
Python基于Excel数据加以反距离加权空间插值并掩膜图层
来源:cnblogs  作者:疯狂学习GIS  时间:2024/4/10 15:15:08  对本文有异议

本文介绍基于PythonArcPy模块,实现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得到的,如下图所示。

image

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

而记录有北京市部分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插值与掩膜。

具体代码如下。

  1. # -*- coding: utf-8 -*-
  2. # @author: ChuTianjia
  3. import copy
  4. import arcpy
  5. import pandas as pd
  6. from arcpy.sa import *
  7. arcpy.env.workspace=arcpy.GetParameterAsText(0)
  8. csv_path=arcpy.GetParameterAsText(1)
  9. shape_file_path=arcpy.GetParameterAsText(2)
  10. idw_result_dir=arcpy.GetParameterAsText(8)
  11. boundary_file_path=arcpy.GetParameterAsText(3)
  12. mask_result_dir=arcpy.GetParameterAsText(9)
  13. spatial_resolution=arcpy.GetParameterAsText(4)
  14. power=arcpy.GetParameterAsText(5)
  15. look_point=arcpy.GetParameterAsText(6)
  16. max_distance=arcpy.GetParameterAsText(7)
  17. csv_data=pd.read_csv(csv_path,header=0,encoding="gbk")
  18. column_name_list=list(csv_data)
  19. hour_column=csv_data["hour"]
  20. pm_25_list=[[0]*len(csv_data) for i in range(csv_data.shape[1]-3)]
  21. for i in range(3,csv_data.shape[1]):
  22. for index,data in csv_data.iterrows():
  23. pm_25_list[i-3][index]=data[i]
  24. 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"]
  25. field_list_use=copy.deepcopy(field_list)
  26. field_list_use.insert(0,"Name")
  27. # Update the columns in the attribute table
  28. for i in range(len(field_list)):
  29. arcpy.AddField_management(shape_file_path,field_list[i],"SHORT")
  30. delete_num=0
  31. delete_name=[]
  32. with arcpy.da.UpdateCursor(shape_file_path,field_list_use) as cursor:
  33. for row in cursor:
  34. for column_name in column_name_list:
  35. if column_name==row[0]:
  36. for i in range(len(csv_data[column_name])):
  37. row[i+1]=csv_data[column_name][i]
  38. cursor.updateRow(row)
  39. # Find stations that without any data
  40. if row[0] not in column_name_list:
  41. cursor.deleteRow()
  42. delete_num+=1
  43. delete_name.append(row[0])
  44. arcpy.AddWarning("Delete {0} site(s) that do not contain any data, and the site(s) name is(are):".format(delete_num))
  45. for i in delete_name:
  46. arcpy.AddMessage(i)
  47. arcpy.AddMessage("\n")
  48. # Perform IDW interpolation
  49. arcpy.env.extent=boundary_file_path
  50. for i in range(len(field_list)):
  51. idw_result=Idw(shape_file_path,field_list[i],spatial_resolution,power,RadiusVariable(look_point,max_distance))
  52. idw_result_path=idw_result_dir+"\\"+"BJ_"+field_list[i]+".tif"
  53. idw_result.save(idw_result_path)
  54. arcpy.AddMessage("{0} has completed IDW interpolation.".format(field_list[i]))
  55. # Perform mask
  56. tif_file_list=arcpy.ListRasters("BJ_hour_*","TIF")
  57. for raster in tif_file_list:
  58. mask_result=ExtractByMask(raster,boundary_file_path)
  59. mask_result_path=mask_result_dir+"\\"+raster.strip(".tif")+"_Mask.tif"
  60. mask_result.save(mask_result_path)
  61. arcpy.AddMessage("{0} has been masked.".format(raster.strip(".tif")))

3 运行结果

执行上述代码,如果是在ArcMap中直接通过工具箱运行,则可以看到代码运行过程中出现的提示。

例如,下图所示提示可以知道有哪几个站点是没有数据、从而被剔除的。

下图则可以显示出目前代码的运行情况。

同时,在我们设定的结果文件夹中可以看到,23小时的插值图与掩膜图都将自动生成并保存在指定文件夹。

再来看看具体的图片长什么样子。

首先查看IDW插值结果图;我们以当日10时的插值结果图为例进行查看。可以看到其已对北京市边界矢量数据所包含的矩形范围完成了插值。

接下来,查看IDW插值结果图经过掩膜后的图像;我们同样以当日10时的插值、掩膜结果图为例进行查看。可以看到,经过掩膜操作后的图像已经完全符合北京市边界矢量数据的范围。

至此,大功告成。

原文链接:https://www.cnblogs.com/fkxxgis/p/18125409

 友情链接:直通硅谷  点职佳  北美留学生论坛

本站QQ群:前端 618073944 | Java 606181507 | Python 626812652 | C/C++ 612253063 | 微信 634508462 | 苹果 692586424 | C#/.net 182808419 | PHP 305140648 | 运维 608723728

W3xue 的所有内容仅供测试,对任何法律问题及风险不承担任何责任。通过使用本站内容随之而来的风险与本站无关。
关于我们  |  意见建议  |  捐助我们  |  报错有奖  |  广告合作、友情链接(目前9元/月)请联系QQ:27243702 沸活量
皖ICP备17017327号-2 皖公网安备34020702000426号