dendrograms / astrodendro

Generate a dendrogram from a dataset
https://dendrograms.readthedocs.io/
Other
37 stars 36 forks source link

(literal) edge-case problem with periodic neighbours #122

Closed tomr-stargazer closed 9 years ago

tomr-stargazer commented 9 years ago

Hi, For the most part, creating dendrograms with periodic_neighbours (#108) has been working correctly (with the help of #121). However, I recently created a wraparound dendrogram containing 2,177 structures, and discovered that five of those structures had nan values for their positional values, even though their sizes seemed to compute properly:

In [8]: catalog[np.isnan(catalog['x_cen'])].pprint()
_idx  area_ellipse  area_exact      flux      major_sigma    minor_sigma   position_angle     radius     v_cen     v_rms     x_cen y_cen flux_kelvin_kms_deg2
          deg2         deg2          Jy           deg            deg            deg            deg                 km / s                   deg2 K km / s
---- -------------- ---------- ------------- -------------- -------------- -------------- -------------- ----- ------------- ----- ----- --------------------
1260  6.72880856597    17.1875 2402344.57932  1.47463656379   1.0477441366  58.4699370791  1.24299710914   nan  4.3322404787   nan   nan        50.4645520866
1321  5.47662014899     13.125 2012821.87577  1.15082683443  1.09270946325  64.2470089952  1.12139171236   nan 4.45440055903   nan   nan        42.2820919469
1451  1.73235213726     4.8125 996618.077528   1.1081083067 0.358968243592  168.418661404 0.630694611175   nan 3.82949792643   nan   nan        20.9353334725
1517 0.218526355915        1.0 258177.577029 0.320772473391 0.156425947805  145.561964693 0.224002540566   nan 1.29269923017   nan   nan        5.42337510437
1524 0.927216048493        2.5 546511.560744 0.791179090111 0.269096595562  178.762590177 0.461414780462   nan 1.90060983097   nan   nan        11.4802270085

In the viewer, these structures look like this: stupid_viewer_view

In case you missed that, here's some zoomed in-views: image image

(image axes are in coordinate space.)

Structure 1517 (the smallest one) looks more like this -- note that part of its footprint is a single-pixel-wide strip (on this velocity channel) at the far right edge of the data : image image

I'm currently adding a few more tests to test_compute and test_analysis that may explore what's going on here. It's also possible that the fault is in the WCS conversion -- I'll play around with putting a warning if the catalog gets passed a nan.

ChrisBeaumont commented 9 years ago

Structure 1517 looks small enough that you should be able to directly extract a mask of pixel locations, and from that build a simple script to reproduce the bug (ideally passing the dendrogram and catalog steps, since I bet the logic bug is somewhere in the structure or statistic classes)

tomr-stargazer commented 9 years ago

Looks like the error may come from here, in the self.wcs.all_pix2world call (and there's a cryptic todo comment) :

https://github.com/dendrograms/astrodendro/blob/master/astrodendro/analysis.py#L293

class SpatialBase(object):
    ... 
    def _world_pos(self):
        xyz = self.stat.mom1()[::-1]
        if self.wcs is None:
            return xyz[::-1] * u.pixel
        else:
            # TODO: set units correctly following WCS
            # We use origin=0 since the indices come from Numpy indexing
            return self.wcs.all_pix2world([xyz], 0).ravel()[::-1]

Structure 1517's (ScalarStatistic) mom1 looks like this

In [103]: stat1517.mom1()
Out[103]: [123.44708399481173, 9.5652454276340659, 1440.223699313673]

when the data's shape tuple is (246, 40, 1440). I'm still playing with this -- maybe wcs.all_pix2world always throws up when the nominal pixel value is bigger than the data, in which case maybe I'll throw a modulo 1440 operation somewhere before passing mom1 to _world_pos.

tomr-stargazer commented 9 years ago

I believe I've solved this problem with #124 , but let's see if others think that code is ready for merging.

tomr-stargazer commented 9 years ago

wcs_object.wcs.bounds_check(False) fixes this as noted in #124. Closing.