创建用1填充的的数组python_Python:在数组中填充多边形区域

我有一个光栅图像(Tiff格式)和一个多边形区域，在数组中转换为shapefile格式。我希望找到一种优雅的方法来创建一个数组，其中多边形边界内的所有元素都有1个值，而多边形外部的所有元素都有0个值。我的最终目标是用从shapefile派生的数组屏蔽从图像派生的数组。

我有以下问题，谢谢你的帮助：

在使用np.zeros((ds.RasterYSize，ds.RasterXSize))和多边形边界地理空间坐标的像素位置创建空数组之后，用1填充数组中的多边形的最佳解决方案是什么？from osgeo import gdal, gdalnumeric, ogr, osr

import osgeo.gdal

import math

import numpy

import numpy as np

def world2Pixel(geoMatrix, x, y):

"""

Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate

the pixel location of a geospatial coordinate

(source http://www2.geog.ucl.ac.uk/~plewis/geogg122/vectorMask.html)

geoMatrix

[0] = top left x (x Origin)

[1] = w-e pixel resolution (pixel Width)

[2] = rotation, 0 if image is "north up"

[3] = top left y (y Origin)

[4] = rotation, 0 if image is "north up"

[5] = n-s pixel resolution (pixel Height)

"""

ulX = geoMatrix[0]

ulY = geoMatrix[3]

xDist = geoMatrix[1]

yDist = geoMatrix[5]

rtnX = geoMatrix[2]

rtnY = geoMatrix[4]

pixel = np.round((x - ulX) / xDist).astype(np.int)

line = np.round((ulY - y) / xDist).astype(np.int)

return (pixel, line)

# Open the image as a read only image

ds = osgeo.gdal.Open(inFile,gdal.GA_ReadOnly)

# Get image georeferencing information.

geoMatrix = ds.GetGeoTransform()

ulX = geoMatrix[0] # top left x (x Origin)

ulY = geoMatrix[3] # top left y (y Origin)

xDist = geoMatrix[1] # w-e pixel resolution (pixel Width)

yDist = geoMatrix[5] # n-s pixel resolution (pixel Height)

rtnX = geoMatrix[2] # rotation, 0 if image is "north up"

rtnY = geoMatrix[4] #rotation, 0 if image is "north up"

# open shapefile (= border of are of interest)

shp = osgeo.ogr.Open(poly)

source_shp = ogr.GetDriverByName("Memory").CopyDataSource(shp, "")

# get the coordinates of the points from the boundary of the shapefile

source_layer = source_shp.GetLayer(0)

feature = source_layer.GetNextFeature()

geometry = feature.GetGeometryRef()

pts = geometry.GetGeometryRef(0)

points = []

for p in range(pts.GetPointCount()):

points.append((pts.GetX(p), pts.GetY(p)))

pnts = np.array(points).transpose()

print pnts

pnts

array([[ 558470.28969598, 559495.31976318, 559548.50931402,

559362.85560495, 559493.99688721, 558958.22572622,

558529.58862305, 558575.0174293 , 558470.28969598],

[ 6362598.63707171, 6362629.15167236, 6362295.16466266,

6362022.63453845, 6361763.96246338, 6361635.8559779 ,

6361707.07684326, 6362279.69352024, 6362598.63707171]])

# calculate the pixel location of a geospatial coordinate (= define the border of my polygon)

pixels, line = world2Pixel(geoMatrix,pnts[0],pnts[1])

pixels

array([17963, 20013, 20119, 19748, 20010, 18939, 18081, 18172, 17963])

line

array([35796, 35734, 36402, 36948, 37465, 37721, 37579, 36433, 35796])

#create an empty array with value zero using

data = np.zeros((ds.RasterYSize, ds.RasterXSize))