CODEWITHSUNDEEP
pythonarraysalgorithmpython-3.xnumpy
Geoff
Geoff

Reputation: 2491

How can I pull the values from a NumPy 2-D array that are inside a polygon laid on top of the array?

My goal is the following: given a 2-D array and a polygon that can be laid on top of the array, I want to pull out the values within the polygon (boundaries included) and sum the values.

Suppose I have a 2-D array like this:

1 6 8 4 2 3
5 4 1 3 7 9
1 0 2 3 4 8
6 7 0 5 7 4
0 1 2 5 4 2

And suppose I have a polygon with the following vertices:

(1, 2)
(2, 3)
(3, 3)
(4, 2)
(4, 1)
(3, 1)
(2, 1)
(1, 2)

Laying the polygon on top of the array would result in the following:

_ _ _ _ _ _
_ _ 1 _ _ _
_ 0 2 3 _ _
_ 7 0 5 _ _
_ 1 2 _ _ _

Finally, I want to sum the values. So the final output would be 21.

What's the most efficient way to do this? I've seen references to matplotlib.path and Shapely, but I haven't found anything that quite works for what I need unless I'm just missing something (which is certainly possible). It seems like this functionality should be built into NumPy, but I haven't come across it if it is.

Why I need this: I have an ASCII grid representing world-wide population, and I have political boundaries. I'm wanting to overlay the political boundaries on the population grid in order to determine that location's total population count.

Upvotes: 1

Views: 1791

Answers (2)

Divakar
Divakar

Reputation: 221574

Steps involved -

  • Create a boolean array of the same shape as the image array.

  • Fill the outline using the given outline points.

  • Create the polygon blob by filling that hole created by the outlined points.

  • Boolean-index into image array with the filled mask and get the total sum.

Implementation -

from scipy.ndimage.morphology import binary_fill_holes as imfill

mask = np.zeros(img.shape,dtype=bool)
mask[idx[:,0], idx[:,1]] = 1
out = img[imfill(mask)].sum()

Sample run -

1) Input image and polygon outline pts :

In [107]: img
Out[107]: 
array([[1, 6, 8, 4, 2, 3],
       [5, 4, 1, 3, 7, 9],
       [1, 0, 2, 3, 4, 8],
       [6, 7, 0, 5, 7, 4],
       [0, 1, 2, 5, 4, 2]])

In [108]: idx
Out[108]: 
array([[1, 2],
       [2, 3],
       [3, 3],
       [4, 2],
       [4, 1],
       [3, 1],
       [2, 1],
       [1, 2]])

2) Proposed codes :

In [109]: mask = np.zeros(img.shape,dtype=bool)
     ...: mask[idx[:,0], idx[:,1]] = 1
     ...: out = img[imfill(mask)].sum()
     ...:

3) Analyze the results :

In [134]: imfill(mask) # Polygon mask
Out[134]: 
array([[False, False, False, False, False, False],
       [False, False,  True, False, False, False],
       [False,  True,  True,  True, False, False],
       [False,  True,  True,  True, False, False],
       [False,  True,  True, False, False, False]], dtype=bool)

In [135]: img[imfill(mask)]
Out[135]: array([1, 0, 2, 3, 7, 0, 5, 1, 2])

In [136]: img[imfill(mask)].sum()
Out[136]: 21

Upvotes: 2

Paul Panzer
Paul Panzer

Reputation: 53029

I don't know it's the most efficient (most probably not), but it's vectorised and it works. a is the 'picture', b is the polygon.

i, k = np.indices(a.shape)
z = i+k*1j
p = b[:,0]+b[:,1]*1j
m = np.unwrap(np.angle(z[..., None]-p), axis=-1)
m = ~np.isclose (m[...,-1], m[...,0]) | np.any(np.isclose(z[..., None], p), axis=-1)
res = np.sum(a[m])
# 21 -> res

Upvotes: 1

Related Questions