Python遥感开发之批量掩膜和裁剪
- 1.使用arcpy进行批量掩膜
- 1.1 批量掩膜代码
- 1.2 单个掩膜代码
- 2.使用GDAL进行批量掩膜
- 3.使用rasterio进行批量裁剪
前言:主要介绍了使用arcpy、gdal、rasterio对遥感影像进行批量掩膜和裁剪。
1.使用arcpy进行批量掩膜
注意:arcpy是基于python2.7版本的,也就是arcgis里面自带的,使用的时候添加arcgis的python2.7,即可使用arcpy。
环境如图所示:
1.1 批量掩膜代码
- filepath:需要掩膜的tif文件夹
- outfile :掩膜好的tif保存的文件夹
- inMaskData :可以是shp,可以是tif
- unicode:因为python2.7下直接支持中文路径,使用unicode可以支持中文路径
# coding=UTF-8
import arcpy,os,glob
from arcpy import env
if __name__ == '__main__':
filepath = unicode(r"D:\0000Desktop\AAAAAwork\研究方向\data\河南NDVI", "utf-8")
outfile = unicode(r"D:\0000Desktop\AAAAAwork\研究方向\data\鹿邑NDVI", "utf-8")
inMaskData = unicode(r"D:\0000Desktop\AAAAAwork\研究方向\矢量图\鹿邑县\luyixian.shp", "utf-8") #shp文件可以换成tif文件
arcpy.CheckOutExtension("Spatial")
env.workspace = filepath
os.chdir(filepath)
names = glob.glob("*.tif")
for name in names:
arcpy.gp.ExtractByMask_sa(name, inMaskData, outfile + "\\" + name.split(".")[0] + ".tif")
print(name + '-----掩膜成功')
#删除多余的文件
for file_i in glob.glob(os.path.join(outfile, '*.xml')):
os.remove(file_i)
for file_i in glob.glob(os.path.join(outfile, '*.tfw')):
os.remove(file_i)
1.2 单个掩膜代码
# coding=UTF-8
import arcpy
if __name__ == '__main__':
filepath = unicode(r"D:\0000Desktop\AAAAAwork\研究方向\data\河南NDVI\modis_ndvi_2019_10_08.tif", "utf-8")
outfile = unicode(r"D:\0000Desktop\AAAAAwork\研究方向\data\鹿邑NDVI\new_modis_ndvi_2019_10_08.tif", "utf-8")
inMaskData = unicode(r"D:\0000Desktop\AAAAAwork\研究方向\矢量图\鹿邑县\luyixian.shp", "utf-8") #shp文件可以换成tif文件
arcpy.CheckOutExtension("Spatial")
arcpy.gp.ExtractByMask_sa(filepath, inMaskData, outfile)
2.使用GDAL进行批量掩膜
注意:GDAL是基于Python3.x系列的,需要自己安装GDAL,使用的时候arcpy使用的时候一定要区分开。
from osgeo import gdal
import os
import glob
def mask_gdal(inMaskData,filepath,outfile):
"""
GDAL掩膜提取
:param inMaskData: 圈选范围的路径 shp文件可以换成tif文件
:param filepath: 要掩膜的tif文件
:param outfile: 掩膜好的tif文件
:return:
"""
dataset = gdal.Open(filepath) # 打开遥感影像
gdal.Warp(outfile, dataset, cutlineDSName=inMaskData, cropToCutline=True) # 按掩膜提取
print(filepath + '-----掩膜成功')
if __name__ == '__main__':
filepath = r"D:\0000Desktop\AAAAAwork\研究方向\data\河南NDVI" # 要裁剪的tif文件所在的文件夹
outfile = r"D:\0000Desktop\AAAAAwork\研究方向\data\鹿邑NDVI"
inMaskData = r"D:\0000Desktop\AAAAAwork\研究方向\矢量图\鹿邑县\luyixian.shp" # 圈选范围的路径 shp文件可以换成tif文件
os.chdir(filepath)
names = glob.glob("*.tif") # 读取文件
for name in names:
out = outfile + "\\" + name.split(".")[0] + ".tif" # 按照圈选范围提取出的影像所存放的路径
mask_gdal(inMaskData,name,out)
3.使用rasterio进行批量裁剪
注意:rasterio也是基于python3.x系列的,需要自己安装rasterio包
import os
import glob
import rasterio as rio
import rasterio.mask
import numpy as np
from geopandas import GeoDataFrame
def clip_raster(inMaskData,filepath,outfile):
"""
raster裁剪提取
:param inMaskData:圈选范围的路径 shp文件可以换成tif文件
:param filepath:要裁剪的tif文件
:param outfile:裁剪好的tif文件
:return:
"""
#shp转geojson
shpdata = GeoDataFrame.from_file(inMaskData)
geo = shpdata.geometry[0]
feature = [geo.__geo_interface__]
#裁剪
src = rio.open(filepath)
out_image, out_transform = rio.mask.mask(src, feature, crop=True, nodata=np.nan)
out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff","height": out_image.shape[1],"width": out_image.shape[2],"transform": out_transform})
band_mask = rasterio.open(outfile, "w", **out_meta)
band_mask.write(out_image)
print(filepath + '-----裁剪成功')
if __name__ == '__main__':
filepath = r"D:\0000Desktop\AAAAAwork\研究方向\data\河南NDVI" # 要裁剪的tif文件所在的文件夹
outfile = r"D:\0000Desktop\AAAAAwork\研究方向\data\鹿邑NDVI"
inMaskData = r"D:\0000Desktop\AAAAAwork\研究方向\矢量图\鹿邑县\luyixian.shp" # 圈选范围的路径 shp文件可以换成tif文件
os.chdir(filepath)
names = glob.glob("*.tif") # 读取文件
for name in names:
out = outfile + "\\" + name.split(".")[0] + ".tif" # 按照圈选范围提取出的影像所存放的路径
clip_raster(inMaskData,name,out)