additional pydl Mangle bugs

[mdipompe]

I've found a few more apparent issues in the mangle codes, when using is_in_polygon. I've attached some test files - a rectangular lune polygon (with 4 caps) and a single circular polygon (1 cap), both in both polygon and fits formats, along with random points in the circular region (fits file with tags ra and dec). These were generated with idlutils and/or mangle itself to illustrate.

test_files.zip

import numpy as np
import pydl.pydlutils.mangle as mangle
from astropy.table import Table

poly1=mangle.read_fits_polygons('lune.fits')
poly2=mangle.read_fits_polygons('circ.fits')
rand=Table.read('rand_circ.fits')
pos=np.transpose(np.array([rand['RA'],rand['DEC']]))

isin1=mangle.is_in_polygon(poly1,pos)

First, note the need to transpose the position data table - it might be nice to adjust the mangle code to do this transpose internally, or avoid it by rearranging some dimensions.

Second, the is_in_polygon call will crash because of a dimension problem:

in is_in_polygon(polygon, points, caps)

    512         if is_cap_used(polygon.use_caps, cap):

    513             in_polygon &= is_in_cap(polygon.x[icap, :],
--> 514                                     polygon.cm[icap], points)
    515     return in_polygon
    516 

/Users/Mike/computing/python/anaconda/lib/python2.7/site-packages/pydl/pydlutils/mangle.py in is_in_cap(x, cm, points)
    482         A boolean vector giving the result for each point.
    483     """
--> 484     return cap_distance(x, cm, points) >= 0.0
    485 
    486 

/Users/Mike/computing/python/anaconda/lib/python2.7/site-packages/pydl/pydlutils/mangle.py in cap_distance(x, cm, points)
    402     else:
    403         raise ValueError("Inappropriate shape for point!")
--> 404     dotprod = np.dot(xyz, x)
    405     cdist = np.degrees(np.arccos(1.0 - np.abs(cm)) - np.arccos(dotard))I
    406     if cm < 0:

ValueError: shapes (100,3) and (4,3) not aligned: 3 (dim 1) != 4 (dim 0)

I tracked this down to an extra dimension of length 1 in the polygon cap parameters (x and cm). I was able to fix this (though there might be other ways) adjusting line 513 from:

in_polygon &= is_in_cap(polygon.x[icap, :],
                                    polygon.cm[icap], points)

to

in_polygon &= is_in_cap(np.squeeze(polygon.x)[icap, :],
                                    np.squeeze(polygon.cm)[icap], points)

This works, and gives correct results (compared to idlutils) for any random set.

However, if instead I'm using the circular polygon with a single cap, I get an indexing error

isin1=mangle.is_in_polygon(poly2,pos)
is_in_polygon(polygon, points, ncaps)
    511     for icap in range(usencaps):
    512         if is_cap_used(polygon.use_caps, icap):
--> 513             in_polygon &= is_in_cap(np.squeeze(polygon.x)[icap, :],
    514                                     np.squeeze(polygon.cm)[icap], points)
    515     return in_polygon

IndexError: too many indices for array

I was able to circumvent this with a conditional in is_in_polygon based on the number of caps, but was unable to get it to return correct results.

Finally, I noticed that if I try reading in these polygons in Mangle format using read_mangle_polygons, none of this works because the polygons don't have the appropriate attributes (e.g. .x or .cm). This seems to be a bug in read_mangle_polygons.

[weaverba137]

First, note the need to transpose the position data table - it might be
nice to adjust the mangle code to do this transpose internally, or avoid it
by rearranging some dimensions.

If you go back to the original idlutils functions, you will notice that they accept two types of inputs: RA, Dec and unit 3-vector. The number of columns (2 versus 3) is used to distinguish between these types in the PyDL version. This is documented, see http://pydl.readthedocs.io/en/latest/api/pydl.pydlutils.mangle.is_in_polygon.html#pydl.pydlutils.mangle.is_in_polygon

Second, the is_in_polygon call will crash because of a dimension problem:

However, if instead I'm using the circular polygon with a single cap, I get an indexing error

Finally, I noticed that if I try reading in these polygons in Mangle format using read_mangle_polygons, none of this works because the polygons don't have the appropriate attributes (e.g. .x or .cm). This seems to be a bug in read_mangle_polygons.

is_in_polygon() takes a single polygon.

import numpy as np
from astropy.io import fits
from pydl.pydlutils.mangle import read_mangle_polygons, is_in_polygon
with fits.open('rand_circ.fits') as hdulist:
    data = hdulist[1].data

points = np.vstack((data.RA, data.DEC)).T
lune = read_mangle_polygons('lune.ply')
circ = read_mangle_polygons('circ.ply')
is_in_polygon(lune[0], points)
is_in_polygon(circ[0], points)

In IDL you can get away with passing, e.g. lune instead of lune[0], because a single IDL structure is basically indistinguishable from an array of structures (of the same type) of length 1. In Python, you have to be more careful about this distinction, & that's almost always good thing. A low priority tweak would be to demote PolygonList of length 1 to a single polygon.

There is a genuine bug in how polygons are read from FITS files, but that has to do with the object used to store the polygon data rather than the is_in_polygon() code itself.

Finally, I'm somewhat confused, because the polygon contained in lune.fits or lune.ply contains four caps, but the formal definition of a lune is defined by two caps. See https://en.wikipedia.org/wiki/Lune_%28geometry%29.

Let me know if I've missed anything here.

[mdipompe]

Ok, thanks for clarifying. One reason for my confusion is I'm actually used to using is_in_window.pro in idlutils, which takes one or several polygons, rather than is_in_polygon.

Now I'm able to get things to run reading in using read_mangle_polygons and specifying the index of the polygon. This doesn't work using read_fits_polygons, but it seems you've identified a bug there.

However, it seems that when using the single cap circular polygon, I don't get the right answer. By definition, all of the points in rand_circ.fits are inside this polygon, but is_in_polygon returns "False" for all of them. Using the "lune", I get agreement with the idlutils implementation (some of the rand_circ.fits points fall inside that polygon, others don't).

As for the definition of a lune, apologies for the confusion. I had a professor (Adam Myers, you probably know him), who defined a lune in his course as a rectangle on the sphere bounded by 4 caps - two for latitude (or dec) and two for longitude (or RA). Even though this isn't officially correct, for better or worse that definition has stuck in my head - I'll be more careful in the future!

[weaverba137]

I'm not able to reproduce the behavior you describe for the single cap polygon:

import numpy as np
from astropy.io import fits
from pydl.pydlutils.mangle import read_mangle_polygons, is_in_polygon
with fits.open('rand_circ.fits') as hdulist:
    data = hdulist[1].data

points = np.vstack((data.RA, data.DEC)).T
circ = read_mangle_polygons('circ.ply')
assert is_in_polygon(circ[0], points).all()

But make sure you are using the latest tagged version, 0.5.1, or any later commit.

[mdipompe]

Working now. Had a previous version loaded in the kernel that wasn't reset properly. Thanks for the clarification on all this, and for the work you've put into porting the idlutils!

[weaverba137]

Closing in preparation for next tag.