Python：填充一個數組中的多邊形區域

Question

我在一個數組中轉換了一個格式文件格式的光柵圖像（Tiff格式）和多邊形區域。我希望找到一種優雅的方式來創建一個數組，其中多邊形的邊界內的所有元素都有1個值，多邊形外的所有元素的值都爲0.我的最終目標是掩蓋從圖像派生的數組，其中數組從派生自shapefile 。

我有以下問題，並感謝幫助：

之後創建使用np.zeros（（ds.RasterYSize，ds.RasterXSize））空數組和地理空間的像素位置的座標我的多邊形的邊框，填充數組中的多邊形的最佳解決方案是什麼？

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)) 

Answer 1

這實質上是一個point-in-polygon問題。

這裏有一個小庫來解決這個問題。它從this頁面進行了一些修改，使其更具可讀性。

pip.py

#From http://www.ariel.com.au/a/python-point-int-poly.html 
# Modified by Nick ODell 
from collections import namedtuple 

def point_in_polygon(target, poly): 
    """x,y is the point to test. poly is a list of tuples comprising the polygon.""" 
    point = namedtuple("Point", ("x", "y")) 
    line = namedtuple("Line", ("p1", "p2")) 
    target = point(*target) 

    inside = False 
    # Build list of coordinate pairs 
    # First, turn it into named tuples 

    poly = map(lambda p: point(*p), poly) 

    # Make two lists, with list2 shifted forward by one and wrapped around 
    list1 = poly 
    list2 = poly[1:] + [poly[0]] 
    poly = map(line, list1, list2) 

    for l in poly: 
     p1 = l.p1 
     p2 = l.p2 

     if p1.y == p2.y: 
      # This line is horizontal and thus not relevant. 
      continue 
     if max(p1.y, p2.y) < target.y <= min(p1.y, p2.y): 
      # This line is too high or low 
      continue 
     if target.x < max(p1.x, p2.x): 
      # Ignore this line because it's to the right of our point 
      continue 
     # Now, the line still might be to the right of our target point, but 
     # still to the right of one of the line endpoints. 
     rise = p1.y - p2.y 
     run = p1.x - p2.x 
     try: 
      slope = rise/float(run) 
     except ZeroDivisionError: 
      slope = float('inf') 

     # Find the x-intercept, that is, the place where the line we are 
     # testing equals the y value of our target point. 

     # Pick one of the line points, and figure out what the run between it 
     # and the target point is. 
     run_to_intercept = target.x - p1.x 
     x_intercept = p1.x + run_to_intercept/slope 
     if target.x < x_intercept: 
      # We almost crossed the line. 
      continue 

     inside = not inside 

    return inside 

if __name__ == "__main__": 
    poly = [(2,2), (1,-1), (-1,-1), (-1, 1)] 
    print point_in_polygon((1.5, 0), poly) 

Answer 2

接受的答案不適合我的工作。

我結束了使用勻稱的庫。

sudo pip install shapely 

代碼：

import shapely.geometry 

poly = shapely.geometry.Polygon([(2,2), (1,-1), (-1,-1), (-1, 1)]) 
point = shapely.geometry.Point(1.5, 0) 

point.intersects(poly)