利用Sentinel-2 A/B、Landsat8、高分等多光譜遙感數(shù)據(jù)，生成長時序植被指數(shù)數(shù)據(jù)，是監(jiān)測農(nóng)作物生長的常見方法。進行時序分析時，為了從高維時序圖像中的指定樣本地塊中快速提取數(shù)據(jù)，常采用編程的方式批量處理。gdal和rasterio都可以讀取柵格數(shù)據(jù)，但是，兩者進行讀取的過程有區(qū)別。

圖像進行分割后，生成大量的多邊形，保存在shp文件中。然后需要計算每個多邊形覆蓋范圍內(nèi)的圖像的紋理特征。gdal和rasterio都能讀取shp文件，并訪問指定多邊形后生成掩膜后，再計算紋理特征。但是兩者生成掩膜的過程和結(jié)果差別比較大，經(jīng)過比較和分析，發(fā)現(xiàn)使用rasterio能夠更快地生成合格的掩膜，相關(guān)代碼如下：

import gdal

from shapely.Geometry import Polygon

from affine import Affine

raster = gdal.Open('p.tif')

shp = gdal.Open('p.shp')

lyr = shp.GetLayer()

feat = lyr.GetNextFeature()

# TODO: 使用rasterio和shapely處理多邊形掩膜，速度快

# 將ogr.Geometry轉(zhuǎn)換為shapely.Geometry.Polygon

geom = [feat.GetGeometryRef()]

geom_dict =eval(geom[0].ExportToJson())

geom = [Polygon(geom_dict['coordinates'][0]).__geo_interface__]

# TODO: rasterio不支持gdal形式的transform，轉(zhuǎn)換為Affine

transform = raster.GetGeoTransform()

pixelWidth = transform[1]

pixelHeight = transform[5]

transform = Affine.from_gdal(xmin, pixelWidth, 0, ymax, 0, pixelHeight)

datamask = np.logical_not(features.geometry_mask(geom, [xcount, ycount], transform))

datamask就是多邊形生成的掩膜，用它就可以計算多邊形覆蓋范圍內(nèi)的圖像的各種特征，比如均值，方差，各種紋理特征，等等