一、前言

前面一篇文章（ 使用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()