Python遥感开发之批量掩膜和裁剪

前言：主要介绍了使用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)