目錄
一、前言
前面一篇文章（使用Python實現子區域數據分類統計）講述了通過geopandas庫實現對子區域數據的分類統計，說白了也就是如何根據一個shp數據對另一個shp數據進行切割。本篇作為上一篇內容的姊妹篇講述如何采用優雅的方式根據一個shp數據對一個柵格影像數據進行切割。廢話不多說，直接進入主題。
二、涉及到的技術
本方案涉及以下技術點：
- geopandas：已經在上一篇文章中簡單介紹。
- numpy：這是一個開源的數據分析處理庫，非常高效、簡潔。
- rasterio：這是一個開源的影像處理庫，非常好用，基本涵蓋了常用的功能。也可以配合numpy進行數據計算。
- 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文件即可。