Extract coordinates enclosed by matplotlib patch.
I have created an ellipse using matplotlib.patches.ellipse as shown below:
patch = mpatches.Ellipse(center, major_ax, minor_ax, angle_deg, fc='none', ls='solid', ec='g', lw='3.')
I need a list of all integer coordinates enclosed within this patch. That is, if I were going to plot this ellipse along with each integer point in the same grid, how many of those points are enclosed in the ellipse?
I tried to see if I can extract the equation of the ellipse so that I can loop through each point and see if it hits the line, but I can't find an obvious way to do this, it gets more complicated since the main axis of the ellipse can be oriented to any angle. The information for this should be stored in patches somewhere, but I can't find it.
Any advice on this would be much appreciated.
source to share
Ellipse
objects have a method contains_point
that will return 1 if the point is in an ellipse, 0 is another wise.
Stealing from @DrV's answer:
import matplotlib.pyplot as plt
import matplotlib.patches
import numpy as np
# create an ellipse
el = matplotlib.patches.Ellipse((50,-23), 10, 13.7, 30, facecolor=(1,0,0,.2), edgecolor='none')
# calculate the x and y points possibly within the ellipse
y_int = np.arange(-30, -15)
x_int = np.arange(40, 60)
# create a list of possible coordinates
g = np.meshgrid(x_int, y_int)
coords = list(zip(*(c.flat for c in g)))
# create the list of valid coordinates (from untransformed)
ellipsepoints = np.vstack([p for p in coords if el.contains_point(p, radius=0)])
# just to see if this works
fig = plt.figure()
ax = fig.add_subplot(111)
ax.add_artist(el)
ep = np.array(ellipsepoints)
ax.plot(ellipsepoints[:,0], ellipsepoints[:,1], 'ko')
plt.show()
This will give you the output as shown below:
source to share
If you really want to use the methods suggested matplotlib
then:
import matplotlib.pyplot as plt
import matplotlib.patches
import numpy as np
# create an ellipse
el = matplotlib.patches.Ellipse((50,-23), 10, 13.7, 30, facecolor=(1,0,0,.2), edgecolor='none')
# find the bounding box of the ellipse
bb = el.get_window_extent()
# calculate the x and y points possibly within the ellipse
x_int = np.arange(np.ceil(bb.x0), np.floor(bb.x1) + 1, dtype='int')
y_int = np.arange(np.ceil(bb.y0), np.floor(bb.y1) + 1, dtype='int')
# create a list of possible coordinates
g = np.meshgrid(x_int, y_int)
coords = np.array(zip(*(c.flat for c in g)))
# create a list of transformed points (transformed so that the ellipse is a unit circle)
transcoords = el.get_transform().inverted().transform(coords)
# find the transformed coordinates which are within a unit circle
validcoords = transcoords[:,0]**2 + transcoords[:,1]**2 < 1.0
# create the list of valid coordinates (from untransformed)
ellipsepoints = coords[validcoords]
# just to see if this works
fig = plt.figure()
ax = fig.add_subplot(111)
ax.add_artist(el)
ep = np.array(ellipsepoints)
ax.plot(ellipsepoints[:,0], ellipsepoints[:,1], 'ko')
Seems to work:
(Scaling shows that even the dots hanging on the edge are inside.)
The fact is that it matplotlib
processes ellipses in the form of transformed circles (translation, rotation, scaling, something affine). If the transformation is applied in reverse order, the result is a unit circle at the origin, and it is very easy to check if a point is in it.
Just a warning: get_window_extent
May not be very reliable as it seems to use a spline approximation of the circle. Also, see the comment tcaswell
for rendering dependencies.
To find a more reliable bounding box, you can:
-
create a horizontal and vertical vector in the coordinates of the plot (their position is not important, [[0,0], [1,0]) and ([0,0], [0,1]) will do)
-
convert these vectors to ellipse coordinates (
get_transform
etc.) -
find in the ellipse coordinate system (i.e. the system where the ellipse is the unit circle around the origin), four tangent circles that are parallel to these two vectors
-
find the intersection points of vectors (4 intersections, but 2 diagonals will be enough)
-
convert the intersection points back to plot coordinates
This will give a precise (but of course limited numerical precision) square bounding box.
However, you can use a simple approximation:
- all possible points are inside a circle whose center coincides with the center of the ellipse and whose diameter is the same as that of the main axis of the ellipse
In other words, all possible points are in a square bounding box that is between x0 + -m / 2, y0 + -m / 2, where (x0, y0) is the center of the ellipse and m is the major axis.
source to share
I would like to suggest another solution that uses the Path
object method contains_points()
instead contains_point()
:
First, we get the coordinates of the ellipse and turn it into an object Path
:
elpath=Path(el.get_verts())
(NOTE that el.get_paths()
won't work for some reason.)
Then call the path contains_points()
:
validcoords=elpath.contains_points(coords)
Below, I compare @ tacaswell's solution (method 1), @Drv (method 2) and my own (method 3) (I enlarged the ellipse by about 5 times):
import numpy
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from matplotlib.path import Path
import time
#----------------Create an ellipse----------------
el=Ellipse((50,-23),50,70,30,facecolor=(1,0,0,.2), edgecolor='none')
#---------------------Method 1---------------------
t1=time.time()
for ii in range(50):
y=numpy.arange(-100,50)
x=numpy.arange(-30,130)
g=numpy.meshgrid(x,y)
coords=numpy.array(zip(*(c.flat for c in g)))
ellipsepoints = numpy.vstack([p for p in coords if el.contains_point(p, radius=0)])
t2=time.time()
print 'time of method 1',t2-t1
#---------------------Method 2---------------------
t2=time.time()
for ii in range(50):
y=numpy.arange(-100,50)
x=numpy.arange(-30,130)
g=numpy.meshgrid(x,y)
coords=numpy.array(zip(*(c.flat for c in g)))
invtrans=el.get_transform().inverted()
transcoords=invtrans.transform(coords)
validcoords=transcoords[:,0]**2+transcoords[:,1]**2<=1.0
ellipsepoints=coords[validcoords]
t3=time.time()
print 'time of method 2',t3-t2
#---------------------Method 3---------------------
t3=time.time()
for ii in range(50):
y=numpy.arange(-100,50)
x=numpy.arange(-30,130)
g=numpy.meshgrid(x,y)
coords=numpy.array(zip(*(c.flat for c in g)))
#------Create a path from ellipse vertices------
elpath=Path(el.get_verts())
# call contains_points()
validcoords=elpath.contains_points(coords)
ellipsepoints=coords[validcoords]
t4=time.time()
print 'time of method 3',t4-t3
#---------------------Plot it ---------------------
fig,ax=plt.subplots()
ax.add_artist(el)
ep=numpy.array(ellipsepoints)
ax.plot(ellipsepoints[:,0],ellipsepoints[:,1],'ko')
plt.show(block=False)
I got this runtime:
time of method 1 62.2502269745
time of method 2 0.488734006882
time of method 3 0.588987112045
So the approach is contains_point()
slower. The coordinate transformation method is faster than mine, but when you get irregularly shaped contours / polygons, this method will still work.
Finally, the result is:
source to share