2.Python借助Arcpy实现批量掩膜裁剪、重采样、重投影

学gis的人都会遇到栅格影像的批处理，虽说arcgis工具都是可以批处理，但相对于代码还是繁琐。本着：怎么简单省事怎么来，处理代码如下：

批量掩膜处理：

import os
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')

## 使用掩膜数据对栅格数据进行批量裁剪
def mask_clip(raster_path, mask, out_path):
   n = 0
   rasters = os.listdir(raster_path)  ## 遍历栅格文件，获得所有栅格数据的文件名
   for raster in rasters:
       raster_name = raster.split('.')[0]  ## 文件名去掉后缀
       raster_type = raster.split('.')[-1]  ## 文件名后缀，用于类型筛选，只要tif文件
       if raster_type == 'tif':
           in_raster = os.path.join(raster_path, raster)  ## 获取待裁剪的tif文件的路径（raster只是文件名，读取的时候需要加上路径）
           arcpy.CheckExtension("Spatial")
           out_extract_mask = ExtractByMask(in_raster, mask)  ## 执行裁剪操作，结果为out_extract_mask
           out_raster = os.path.join(out_path, raster_name + '.tif')  ## 输出文件路径，raster_name+'_new.tif'为文件名
           out_extract_mask.save(out_raster)  ## 将out_extract_mask保存为out_raster（路径+文件名）
           n = n + 1
   print str(n) + " rasters are processed !!!"


if __name__ == '__main__':
   raster_path = 'D:/G/year_tem/3_5km'  # 原始栅格文件的文件夹路径
   #mask_file = 'G:/CMIP6_pr/GEOTIFF/yanmo.tif'  # 栅格或矢量的掩膜文件路径  矩形掩膜文件
   mask_file = 'G:/CMIP6_pr/GEOTIFF/5km.tif' # 5km栅格掩膜文件
   out_path = 'D:/G/year_tem/4_cut_5'  # 裁剪后的文件保存路径
   mask_clip(raster_path, mask_file, out_path)

批量重投影：

import os
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')

arcpy.env.workspace = 'G:/CMIP6_tas/GEOTIFF/ACCESS-CM2/'


def GetRaster(file_path, out_path):
   n = 0
   files = os.listdir(file_path)
   for file in files:
       type = file.split(".")
       if type[-1] == "tif":
           print file
           raster = os.path.join(file_path, file)
           out_raster = os.path.join(out_path, file)
           Coordinate_System = "PROJCS['UTM_Zone_50N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['false_easting',500000.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',117.0],PARAMETER['scale_factor',0.9996],PARAMETER['latitude_of_origin',0.0],UNIT['Meter',1.0]]"
           cellsize = 110000  # 1km   50000  50km
           #projecttype = "Albers_Province.prj"  ## 目标投影信息的文件路径,可以是shapefile对应的prj文件
           arcpy.ProjectRaster_management(raster, out_raster, Coordinate_System, "BILINEAR", cellsize, "", "", "")
           n = n + 1
   print str(n) + " rasters are processed !!!"


if __name__ == "__main__":
   data_path = "D:/G/year_tem/1_cut"  ## 栅格数据路径
   out_path = "D:/G/year_tem/2_m"  ## 输出文件路径
   GetRaster(data_path, out_path)

批量重采样：

# coding=utf-8
import os
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')

def GetRaster(file_path, out_path):
   n = 0
   files = os.listdir(file_path)
   for file in files:
       type = file.split(".")
       if type[-1] == "tif":
           print file
           raster = os.path.join(file_path, file)
           out_raster = os.path.join(out_path, file)
           cellsize = 5000  # 5km   50000  50km
           arcpy.Resample_management(raster, out_raster, cellsize,"BILINEAR", )
           n = n + 1
   print str(n) + " rasters are processed !!!"


if __name__ == "__main__":
   data_path = "D:/G/year_tem/2_m"  ## 栅格数据路径
   out_path = "D:/G/year_tem/3_5km"  ## 输出文件路径
   GetRaster(data_path, out_path)

代码很简单，更改输入输出路径就ok！

参考的博客：

https://blog.csdn.net/MLH7M/article/details/120977949?spm=1001.2014.3001.5502

这里注意：重投影的信息，是我的研究区信息，这里要根据自己的研究区地理位置来设置投影带。

祝大家科研顺利！