矢量对栅格数据的裁剪在GIS软件中是基本功能，本文使用Python来实现该功能。其中，矢量数据是面（Polygon)类型，且矢量数据和栅格数据的坐标系一致。在这里，用到的矢量数据为geojson格式，栅格数据为tif格式。

数据读取

import geopandas as gpd
from osgeo import gdal, gdalnumeric

ht = gdal.Open(raster_file_path)
data_df = gpd.read_file(vector_file_path)
poly = data_df.ix[0]['geometry']

接下来，将根据poly这一面状要素，对栅格数据中相应范围进行裁剪。

地理坐标转像素坐标

利用GDAL中的GeoTransform参数来实现。

def geoToPixel(record,GeoTransform):
    """
    http://www.gdal.org/gdal_tutorial.html

    GeoTransform[0] /* top left x */
    GeoTransform[1] /* w-e pixel resolution */
    GeoTransform[2] /* 0 */
    GeoTransform[3] /* top left y */
    GeoTransform[4] /* 0 */
    """
    uper_left_x = float(GeoTransform[0])
    uper_left_y = float(GeoTransform[3])
    pixel_width = float(GeoTransform[1])
    pixel_height = float(GeoTransform[5])

    (mx,my) =(record[0],record[1])
    px = int((mx - uper_left_x) / pixel_width) #x pixel
    py = int((my - uper_left_y) / pixel_height) #y pixel

    return (px,py)

然后，获取面的外接矩形

minX,maxX = min(poly.envelope.exterior.coords.xy[0]),max(poly.envelope.exterior.coords.xy[0])
minY,maxY = min(poly.envelope.exterior.coords.xy[1]),max(poly.envelope.exterior.coords.xy[1])

uper_left_X, uper_left_Y = geoToPixel((minX,maxY),gt)
lower_right_X, lower_right_Y = geoToPixel((maxX,minY),gt)
uper_left_X, uper_left_Y,lower_right_X, lower_right_Y

输出类似

(18450, 5358, 18645, 5466)

根据外接矩形裁剪

pxWidth = int(lower_right_X - uper_left_X)
pxHeight = int(lower_right_Y - uper_left_Y)

clip = ht.ReadAsArray(int(uper_left_X),int(uper_left_Y),int(pxWidth),int(pxHeight))
plt.imshow(clip[0])

现在只是得到了规则矩形裁剪的结果，面状要素往往是不规则的，接下来用掩模处理。

掩模Mask生成

## 计算当前裁剪后栅格的GeoTransform
gt2 = list(gt)
gt2[0] = minX
gt2[3] = maxY

利用opencv的fillPoly实现掩模效果

import opencv
ext_pix = np.array([geoToPixel(points,gt2) for points in poly.exterior.coords])
ext_pix = np.expand_dims(ext_pix,axis=0)
mask = np.ones([pxHeight,pxWidth], np.uint8)

mask = cv2.fillPoly(mask,ext_pix, 0)
plt.imshow(mask)

最终结果

clip = gdalnumeric.choose(mask, (clip, -9999))
plt.imshow(clip[0])