利用一个shp数据中多个面要素同时裁剪栅格的ArcGIS和Python方法实现过程

一、实验目的

  地理学研究过程中，用一个矢量面去提取栅格的方法较为常见，但空间分析中需要将一个栅格数据上的多个空间位置分别进行提取时，逐个提取的方法不够快捷，例如，分别提取中国所有省份的DEM数据时，逐个提取很浪费时间。因此，本实验的目的是利用ArcGIS或Python方法批量提取一个矢量图层下的多个面要素的栅格数据。

二、实验数据和方法

1、实验数据：

（1）2001年中国东部地区物候返青期SOS的遥感影像（MCD12Q1 500M）

（2）中国东北部省份shp图层

（3）中国东北部省份的excel表格
2、代码
（2）提取不同省份的shp图层

（3）剪裁不同省份的SOS图像

3、软件

（1）ArcGIS10.8

（2）Pycharm

三、实现步骤

1、ArcGIS实现步骤

第一步：加载数据

第二步：打开工具箱（ArcToolbox），找到“数据管理工具/栅格/栅格处理/分割栅格（Data Management Tools.tbx / Raster /Raster Processing / Split Raster）”

第三步：设置环境参数

处理完成结果：

软件实现的注意事项：

（1）必须保证矢量数据和栅格数据的投影一致，如果不一致，需要使用“porject”工具手动调整。

（2）设置参数时必须设置环境中的平行处理为“0”

2、代码实现步骤

其内部的逻辑是：先按每个面要素（每个省份）的唯一代码把一个shp数据分割成多个shp，再逐个进行按掩膜提取。

第一步：分割shp

import arcpy
from arcpy import env
from arcpy.sa import *
import pandas as pd
import os
input_path = r"C:\Users\Administrator\Desktop\TEST\shp\NMG_84.shp"    #SHP layer storage path
province_xls = r"C:\Users\Administrator\Desktop\TEST\shp\province.xls"
output_path = r"C:\Users\Administrator\Desktop\TEST\split_shp"

if os.path.exists(output_path)==False:
   os.mkdir(output_path)

env.workspace = r"C:\Users\Administrator\Desktop\TEST\shp"
#Read xls file
df = pd.read_excel(province_xls)
prov_code = df.iloc[:,0]   #To ensure the uniqueness of the province code, duplicate codes cannot be extracted
#Extract shp layers from different provinces
for each_province in prov_code:
   SQL = """ "DZM" = '{}' """.format(each_province)  # Set SQL extraction code
   out_name = output_path + '\\' + 'prov_' + str(each_province) + ".shp"
   print(out_name)
   arcpy.Select_analysis(input_path,out_name,SQL)   # Extract by attribute

第二步：按掩膜提取

# 2.Extract by Mask - Cut
import glob
input_path = r"C:\Users\Administrator\Desktop\TEST\split_shp"    # Obtain the path where all provinces' shps are located
raster_path = r"D:\text\data-text3_new\phenology\MCD12Q2\2_convert\2001_SOS.tif"    #Path to obtain raster images
output_path = r"C:\Users\Administrator\Desktop\TEST\split_tif"
#creat folder
if os.path.exists(output_path)==False:
  os.mkdir(output_path)
#Defining Workspaces
arcpy.env.workspace = input_path
#Obtain all shp format images
shplist = arcpy.ListFiles("*.shp")
shps = glob.glob(os.path.join(input_path,"*.shp"))
for shp in shps:
   name_input = os.path.basename(shp).split(".")[0] + ".tif"
   name_output = os.path.join(output_path, name_input)
   print(name_output)
   arcpy.gp.ExtractByMask_sa(raster_path,shp,name_output)

结果显示：

代码的实现过程中需要特别注意几个事项：

（1）province.xls文件是利用ArcGIS中的转换工具获得

（2）转化完成后需要调整表内的内容，或者根据表内的内容修改代码中的“DZM”

（3）批量注释和批量解除注释的快捷键：“Ctrl+/”