Creating a filled numpy array of ones from a polygon points Python

Question

I would be very gratefull if anyone could help me find a faster solution to my problem.

Heres the scenario:- I have a polygon of float points which I want to map to a grid. The grid cells can be different width and height not uniformed like my image shows. ie rectangular.

I have tried using Image draw but it only uses ints. Converting floats to ints means I have to scale the floats up and remove decimal to keep some precision but image draw will not work with the larger polygon of points.

Is there a more eloquent and fast way of achieving a numpy array of ones (blue) for the filled area of the polygon and zeros (red) for the rest. I have read a little on mesh grid but cant see how it can be of use for this scenario.

Many thanks

The results from the code are

cols = 4
rows = 4
points = [[1535116L, 1725047L], [1535116L, 2125046L], [-464884L, 2125046L], [-464884L, 125046L]]    

bbCut = getPythonBoundBox(points)
cutWidth =  bbCut[1][0]-bbCut[0][0]
scale = float(cutWidth) / float(rows)
###Center data to origin
for p in range(len(points)):
    points[p][0] -= (bbCut[1][0] - bbCut[0][0])/2
    points[p][1] -= (bbCut[1][1] - bbCut[0][1])/2
    points[p][0] /= scale
    points[p][1] /= scale


##move points to Zero
bbCut = getPythonBoundBox(points)

for p in range(len(points)):
    points[p][0] -=bbCut[0][0]
    points[p][1] -=bbCut[0][1]

pointToTuple= []
for p in range(len(points)):
  pointToTuple.append((points[p][0], points[p][1]))

imgWidth = float(rows)
imgHeight = float(cols)

img = Image.new('L', (int(imgWidth), int(imgHeight)), 0) 
draw = ImageDraw.Draw(img)

draw.polygon(pointToTuple, fill=1)

array =  np.reshape(list(img.getdata()), (cols, rows))

############This is the result from the array############
##If you compare this array to the coloured scaled image ive have drawn
##its missing a 1 on the second value in the first row
##and another 1 on the second row 3rd value
##I'm assuming there is some parsing happening here with float to int?                       
array([1, 0, 0, 0])
array([1, 1, 0, 0])
array([1, 1, 1, 1])
array([1, 1, 1, 1])
#########################################################

def getPythonBoundBox(points):
    bigNumber = 10e10
    xmin = bigNumber
    xmax = -bigNumber
    ymin = bigNumber
    ymax = -bigNumber
    g = []
    a = len(points)
    for i in xrange(a):
        if points[i][0] < xmin: xmin = points[i][0]
        if points[i][0] > xmax: xmax = points[i][0]
        if points[i][1] < ymin: ymin = points[i][1]
        if points[i][1] > ymax: ymax = points[i][1]
    p1 = [xmin,ymin]
    g.append(p1)
    p2 = [xmax,ymax]
    g.append(p2)
    return (g)

Answer 1

matplotlib.path.Path has a method contains_points. Hence simply instantiate a path with your polygon points and then check your grid coordinates if they fall within that path. Your grid can have any resolution you want. This is controlled by nx and ny (or alternatively dx and dy) in the code below.

Code:

import numpy as np
import matplotlib.pyplot as plt

from matplotlib.patches import PathPatch
from matplotlib.path import Path

# create a matplotlib path
points = [[1535116L, 1725047L],
          [1535116L, 2125046L],
          [-464884L, 2125046L],
          [-464884L, 125046L],
          [1535116L, 1725047L]]
codes = [Path.MOVETO,
         Path.LINETO,
         Path.LINETO,
         Path.LINETO,
         Path.CLOSEPOLY,
         ]
path = Path(points, codes)

# check the path
fig, (ax1, ax2, ax3) = plt.subplots(1,3)
patch = PathPatch(path, facecolor='k')
ax1.add_patch(patch)
xmin, ymin = np.min(points, axis=0)
xmax, ymax = np.max(points, axis=0)
ax1.set_ylim(ymin,ymax)
ax1.set_xlim(xmin,xmax)
ax1.set_aspect('equal')

# create a grid
nx, ny = 1000, 1000
x = np.linspace(xmin, xmax, nx)
y = np.linspace(ymin, ymax, ny)
xgrid, ygrid = np.meshgrid(x, y)
pixel_coordinates = np.c_[xgrid.ravel(), ygrid.ravel()]

# find points within path
img = path.contains_points(pixel_coordinates).reshape(nx,ny)

# plot
ax2.imshow(img, cmap='gray_r', interpolation='none', origin='lower')

# repeat, but this time specify pixel widths explicitly
dx, dy = 2000, 2000
x = np.arange(xmin, xmax, dx)
y = np.arange(ymin, ymax, dy)
xgrid, ygrid = np.meshgrid(x, y)
pixel_coordinates = np.c_[xgrid.ravel(), ygrid.ravel()]
img = path.contains_points(pixel_coordinates).reshape(len(x), len(y))
ax3.imshow(img, cmap='gray_r', interpolation='none', origin='lower')

UPDATE:

Ok, so this now tests the if any of the corners of each tile are within the path. For some reason, I still get another answer than the picture suggests. How sure are you, that the points that you provide are exact?

Code + image:

import numpy as np
import matplotlib.pyplot as plt

from matplotlib.patches import PathPatch
from matplotlib.path import Path

# create a matplotlib path
points = [[1535116L, 1725047L],
          [1535116L, 2125046L],
          [-464884L, 2125046L],
          [-464884L, 125046L],
          [1535116L, 1725047L]]

codes = [Path.MOVETO,
         Path.LINETO,
         Path.LINETO,
         Path.LINETO,
         Path.CLOSEPOLY,
         ]
path = Path(points, codes)

fig, (ax1, ax2) = plt.subplots(1,2)
patch = PathPatch(path, facecolor='k')
ax1.add_patch(patch)
xmin, ymin = np.min(points, axis=0)
xmax, ymax = np.max(points, axis=0)
ax1.set_ylim(ymin,ymax)
ax1.set_xlim(xmin,xmax)
ax1.set_aspect('equal')

nx, ny = 4, 4
x = np.linspace(xmin, xmax, nx)
y = np.linspace(ymin, ymax, ny)
xgrid, ygrid = np.meshgrid(x, y)
pixel_centers = np.c_[xgrid.ravel(), ygrid.ravel()]

def pixel_center_to_corners(x, y, dx, dy, precision=0.):
    """
    Returns array indexed by (pixel, corner, (x,y))
    """

    # make dx and dy ever so slightly smaller,
    # such that the points fall **inside** the path (not **on** the path)
    dx -= precision
    dy -= precision

    return np.array([(x - dx/2., y - dy/2.), # lower left
                     (x + dx/2., y - dy/2.), # lower right
                     (x + dx/2., y + dy/2.), # upper right
                     (x - dx/2., y + dy/2.), # upper left
    ]).transpose([2,0,1])

# get pixel corners
dx = (xmax - xmin) / float(nx)
dy = (ymax - ymin) / float(ny)
pixel_corners = pixel_center_to_corners(pixel_centers[:,0], pixel_centers[:,1], dx, dy)

# test corners of each pixel;
# set img to True, iff any corners within path;
img = np.zeros((len(pixel_corners)))
for ii, pixel in enumerate(pixel_corners):
    is_inside_path = path.contains_points(pixel)
    img[ii] = np.any(is_inside_path)

img = img.reshape(len(x), len(y))
ax2.imshow(img, cmap='gray_r', interpolation='none', origin='lower')