3

ラスター画像（Tiff形式）とシェープファイル形式のポリゴン領域が配列に変換されています。ポリゴンの境界の内側のすべての要素の値が1で、ポリゴンの外側のすべての要素の値が0である配列を作成するエレガントな方法を見つけたいと思います。最後の目標は、画像から派生した配列をシェープファイルから派生した配列でマスクすることです。 。

私は次の質問があり、助けてくれてありがとう：

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

2 に答える 2

2

これは本質的にポリゴンの点の問題です。

この問題を解決するための小さなライブラリがあります。このページから、読みやすくするためにいくつかの変更を加えています。

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)
于 2012-11-16T19:17:11.890 に答える
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)
于 2015-07-09T18:41:39.607 に答える