使用Python以优雅的方式实现根据shp数据对栅格影像进行切割

一、前言

       前面一篇文章（使用Python实现子区域数据分类统计）讲述了通过geopandas库实现对子区域数据的分类统计，说白了也就是如何根据一个shp数据对另一个shp数据进行切割。本篇作为上一篇内容的姊妹篇讲述如何采用优雅的方式根据一个shp数据对一个栅格影像数据进行切割。废话不多说，直接进入主题。

二、涉及到的技术

本方案涉及以下技术点：

1. geopandas：已经在上一篇文章中简单介绍。
2. numpy：这是一个开源的数据分析处理库，非常高效、简洁。
3. rasterio：这是一个开源的影像处理库，非常好用，基本涵盖了常用的功能。也可以配合numpy进行数据计算。
4. datashader：这是一个开源的大数据可视化库，可以进行遥感影像、矢量数据的可视化。其基于bokeh，bokeh是一个通用的可视化工具，有兴趣的可以参考github，我之前采用Scala语言对其进行了简单的封装，请参考使用bokeh-scala进行数据可视化以及使用bokeh-scala进行数据可视化(2)。

本文基本涉及以上技术。另，最近Github貌似被墙了，所以你懂的。推荐使用Lantern，请自行百度之。

三、优雅切割

       为什么叫优雅的切割，其实我这里倒不是卖弄文字，主要是为了与Gdal的方式相区别。传统的方式可以采用Gdal命令行进行一点点的手动处理，稍微智能化一点可以在python程序中发送控制台语句的方式调用gdal命令。作为程序员我们都是想采用最简单、最不需要手工操作、看上去最舒服的方式。所以我这里称其为优雅的方式。

       我们大致需要经历读取影像、投影转换、读取shp、切割、显示等几个步骤。下面逐一介绍。

3.1 读取影像

       采用rasterio进行影像读取。代码如下：

import rasterio as rio

band = rio.open(path)

       非常简单，只要传入影像的路径即可。rasterio支持tif、hdf格式（亲测）。

3.2 投影转换

       这是比较难的一块，需要注意很多细节。代码如下：

from rasterio.warp import (reproject,RESAMPLING, transform_bounds,calculate_default_transform as calcdt)

affine, width, height = calcdt(src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
    'crs': dst_crs,
    'transform': affine,
    'affine': affine,
    'width': width,
    'height': height,
    'geotransform':(0,1,0,0,0,-1) ,
    'driver': 'GTiff'
})

dst = rio.open(newtiffname, 'w', **kwargs)

for i in range(1, src.count + 1):
    reproject(
        source = rio.band(src, i), 
        destination = rio.band(dst, i), 
        src_transform = src.affine,
        src_crs = src.crs,
        dst_transform = affine,
        dst_crs = dst_crs,
        dst_nodata = src.nodata,
        resampling = RESAMPLING.bilinear)

       首先使用calculate_default_transform函数计算投影转换后的影像尺寸以及仿射变换参数。其中src表示原始影像，src.crs可以取到原始投影，dst_crs位定义的目标投影，与上一篇中介绍shp投影变换时基本一致。

       然后计算投影后的tiff元数据信息。src.meta.copy()读出原始元数据信息并进行拷贝，kwargs.update将原始元数据更新为目标元数据。

       dst = rio.open(newtiffname, 'w', **kwargs)打开一个新的影像其模式w表示写入。

       最后循环原始影像的所有波段，逐一进行投影变换并写入新的影像。其参数一目了然，不再赘述。

       上一个影像的整体截图，以与下述切割后的效果进行对比。

3.3 读取shp

       这在上一篇文章中也已经做了详细描述，不再赘述，需要强调的时此处也需要将shp进行投影转换，使其与我们要处理的影像一致，所以简单的方式就是直接读取影像的投影信息，将shp数据转换到此投影，详情请参考使用Python实现子区域数据分类统计。

3.4 切割

       我们要对一个完整的影像进行切割，可以分为两步。首先将shp数据转换为geojson，然后使用rasterio进行切割。

3.4.1 shp数据转换为geojson

       rasterio进行切割时需要传入的时geojson对象，而不是普通的GeoSeries对象，所以我们需要进行一步转换。代码如下：

from geopandas import GeoSeries
features = [shpdata.geometry.__geo_interface__]

       其中shpdata为读出的shp数据对象。如果我们想要获取shp中的某条空间数据而不是全部，可以采用如下方式：

from geopandas import GeoSeries
features = [GeoSeries(shpdata.geometry[i]).__geo_interface__]

       其中i表示的是取出元素的序号，最后都要采用[]将结果变成数组，因为rasterio最后需要传入一个数组参数。

3.4.2 使用rasterio进行切割

       其实有了前面的准备这一步也就变的简单了，直接调用rio.mask.mask函数，该函数返回该栅格数据与features相交部分的数组结果以及变换信息。代码如下：

import rasterio.mask
out_image, out_transform = rio.mask.mask(src, features, crop=True, nodata=src.nodata)
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(newtiffname, "w", **out_meta)
band_mask.write(out_image)

       其中src为切割前的影像数据，features为上一步得到的shp数据转换后的geojson，crop表示是否对原始影像进行切割，如果为True表示将该geojson的外界框以外的数据全部删除，既缩小原始影像的大小，只保留外界框以内部分，nodata表示无值数据，凡是geojson外部的数据都会转换成此值。

       后面的基本与投影转换后的一致，根据切割的结果生成一个新的影像数据。这样我们就实现了根据shp数据对遥感影像进行切割。效果如下：

四、总结

       本文所介绍的技术可以用于对全国的影像数据进行分省切割，或者省的影像数据进行县市切割等。同理与上一篇文章一致的是凡是这种处理子区域的方式都可以采用此技术。当然本文没有介绍如何对遥感影像进行处理，其实非常简单，当我们读出影像数据之后，其就是一个numpy的array对象，已经变成了纯数学问题，处理完之后只需要附加投影等信息写入新的tiff文件即可。