學(xué)gis的人都會遇到柵格影像的批處理，雖說arcgis工具都是可以批處理，但相對于代碼還是繁瑣。本著：怎么簡單省事怎么來，處理代碼如下：

批量掩膜處理：

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

arcpy.CheckOutExtension('Spatial')

## 使用掩膜數(shù)據(jù)對柵格數(shù)據(jù)進(jìn)行批量裁剪
def mask_clip(raster_path, mask, out_path):
? ?n = 0
? ?rasters = os.listdir(raster_path) ?## 遍歷柵格文件，獲得所有柵格數(shù)據(jù)的文件名
? ?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) ?## 執(zhí)行裁剪操作，結(jié)果為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" ?## 目標(biāo)投影信息的文件路徑,可以是shapefile對應(yīng)的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" ?## 柵格數(shù)據(jù)路徑
? ?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" ?## 柵格數(shù)據(jù)路徑
? ?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

這里注意：重投影的信息，是我的研究區(qū)信息，這里要根據(jù)自己的研究區(qū)地理位置來設(shè)置投影帶。

祝大家科研順利！