1

Some background: I am trying to make region files using FITS data. I used astropy and matplotlib.pyplot to plot 2D histograms of the sections I'm interested in, and wrote a program that bins that data and identifies which bins belong in the region of interest. I then used matplotlib.patches to "highlight" the bins of a 2D histogram that qualified as part of the region.

Here's the plotting part of my code:

indices = range(1, 1295, 1)
NBINS=(360,360)
fig,ax = plt.subplots(1)
ax.hist2d(data['DET1X'], data['DET1Y'], NBINS, range=[[1, 360], [1,360]], cmap='gray')
ax.set_xlabel('DET1X')
ax.set_ylabel('DET1Y')
for i in region_bins:
    cornerx = (indices[i]// 36)*10
    cornery = (indices[i] % 36)*10
    SL = patches.Rectangle((cornerx,cornery), 10, 10, edgecolor='none', facecolor='yellow', alpha = 0.4)
    ax.add_patch(SL)
plt.show()

Here's the result I have: result

Now, I want to draw a perimeter around the region, so basically a polygon where all the yellow edges of bins meet the gray background. The polygon needs to be one shape that is a union of all the squares, so it will be jagged on the left and mostly straight on the bottom and right. It also needs to include all of the little holes, not have holes itself. How can I do this? Previous questions I've looked at deal with situations in languages I'm not familiar with that are very specific to particular projects, as is my own question. Ultimately, I want to be able to make a file that just lists the polygon coordinates of the region with respect to the original 2D histogram. Thoughts?

EllenT
  • 13
  • 4
  • Please explain how exactly that polygon (polygons?) should look like. Does it need to be jagged as a simple union of the yellow squares? The little holes will be holes in the polygon? – JohanC Nov 23 '19 at 00:47
  • Maybe `shapely` can be useful to create a union of rectangles? See e.g. https://stackoverflow.com/questions/25374459/find-holes-in-a-union-of-rectangles – JohanC Nov 23 '19 at 01:19
  • You've outline a somewhat tricky problem, but you might find matplotlib's contouring tools are what you need. e.g., see https://stackoverflow.com/questions/19418901/get-coordinates-from-the-contour-in-matplotlib – keflavich Nov 25 '19 at 17:59

0 Answers0