使用python的geopands和rasterio批量裁剪遥感数据

---

前言

这几天在跑模型的时候，需要对遥感影像批量裁剪和镶嵌，本来想使用arcpy完成，但是由于arcpy裁剪工具Spatial Analyst-Extraction-Extract by Mask工具裁剪出来是栅格数据集的形式，不能直接得到tif文件，于是就想直接用python完成。想着以后可能会经常用到这个工具，于是把裁剪和镶嵌代码一起实现了。

---

一、geopandas和rasterio的安装

这里建议在用conda新环境里安装，因为很多库之间有依赖性，新环境下可以保证安装的都是最新版本或者比较新的版本。我电脑上首先安装了Anaconda库，可以直接使用conda命令创建新环境，创建、删除、激活环境命令可以百度。我电脑上新建的是python3.8环境，刚开始想使用pip在线安装，但是各种报错，于是直接下载离线包。下面给出rasterio包的下载地址：https://www.lfd.uci.edu/~gohlke/pythonlibs/#rasterio
rasterio包依赖于GDAL，可以先在线或者离线安装GDAL（推荐离线安装）。
geopandas依赖库有四个，参考https://geopandas.org/en/latest/getting_started/install.html
四个包在线安装成功了三个，fiona包好像是关联GDAL还是啥安装报错，于是直接离线安装。
安装完四个依赖包后，最后直接使用pip install geopandas完成第一步。

二、批量裁剪和镶嵌

1.数据准备与批量裁剪

我使用的的是中国dem90米分辨率数据，下载完后是hgt格式，可以使用arcmap批量格式转换工具将数据转成tif。


转换完成后使用一下代码进行批量裁剪，这里只简单交代下代码的用法，输入需要裁剪的tif文件夹pathDir，shpfile为裁剪矢量，save_dir为保存路径，nodata=32767是表示背景的值，由于这里数据格式是int16，所以用最大值表示背景。

代码如下（示例）：

import rasterio as rio
from rasterio.mask import mask
from tqdm import tqdm
import glob
import os
import geopandas as gpd
from geopandas import GeoSeries
from osgeo import gdal
import math
import numpy as np

def data_cut():
    pathDir = r'*.tif'
    shpfile = r'Yangtze.shp'
    save_dir = r'' #保存tif
    tif_list = glob.glob(pathDir)
    for i in tqdm(range(len(tif_list))):
        # 读入栅格文件
        rasterfile = tif_list[i]
        rasterdata = rio.open(rasterfile)
        #获取栅格信息
        # profile = rasterdata.profile
        out_meta = rasterdata.meta.copy()
        shpdata = gpd.read_file(shpfile)
        # 投影变换，使矢量数据与栅格数据投影参数一致
        shpdata = shpdata.to_crs(rasterdata.crs)
        
        # 按照所有矢量进行循环裁剪
        for j in range(0, len(shpdata)):
            try:
                # 获取矢量数据的features
                geo = shpdata.geometry[j