用Python实现矢量对栅格数据的裁剪

矢量对栅格数据的裁剪在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])