Welcome to iraf.net Sunday, September 24 2017 @ 04:52 AM GMT

New task to group data on the sky

  • Friday, February 10 2006 @ 09:36 PM GMT
  • Contributed by:
  • Views: 2,059
Code Snippets A new task allows taking a list of RA, Dec, and data, such as image names, that one wants to group into sublists. The grouping is based on a minimum separation in RA and Dec. The SKYGROUP help page is available from this site already, the code itself is in the body of the message and in the upcoming release in the NPROTO package.
# SKYGROUP -- Group coordinate list on the sky.
#
# The input is list of ra (0h-24h or 0d-360d) and dec (-90d to 90d) in the
# first two columns followed by arbitrary data (usually a filename).
# This is complicated by the periodicities at 0h.

procedure skygroup (input, output)
    
file    input                   {prompt="Input list"}
string  output                  {prompt="Output rootname"}
string  extn = ""               {prompt="Optional output extension"}
real    sep = 60                {prompt="Separation between groups (arcsec)"}
string  raunit = "hr"           {prompt="RA unit (hr|deg)"}
bool    keepcoords = yes        {prompt="Keep coordinates in output?"}
string  raformat = "%.2h"       {prompt="Output RA format"}
string  decformat = "%.1h"      {prompt="Output DEC format"}

struct  *fd1, *fd2

begin
        file    in, out, fname, temp1, temp2, temp3
        int     i, j, n, n1
        real    dra, r1, d1, r2, d2, r3, d3, r4
        struct  fmtstr, data1, data2, data3

        # Temporary files.
        fname = mktemp ("tmp")
        in = fname // "in"
        temp1 = fname // "1"
        temp2 = fname // "2"
        temp3 = fname // "3"

        # Set parameters.
        fname = input
        out = output

        # Check for existing output files.
        files (out//"_[0-9][0-9][0-9][0-9]*"//extn, > temp1)
        count (temp1) | scan (n); delete (temp1, verify-)
        if (access(out) || n > 0)
            error (1, "Output files already exist")

        if (raunit == "hr")
            dra = 24
        else
            dra = 360

        if (keepcoords)
            fmtstr = "%d " // raformat // " " // decformat // " %s"
        else
            fmtstr = "%d %s"
        fmtstr += "n"

        # We start by sorting in dec.
        sort (fname, col=2, num+, > in)

        # Find jumps in dec bigger than the separation and then sort
        # in ra and find jumps in ra bigger than separation.  Handle
        # the wrap around at 0h by duplicating to extend beyond 24h.
        # The duplicates will be eliminated during the merging process.

        n = 1
        fd1 = in
        if (fscan (fd1, r1, d1, data1) == EOF)
            error (1, "No data or badly formated data")
        while (fscan (fd1, r2, d2, data2) != EOF) {
            print (r1, d1, data1, >> temp1)
            if (r1 / max (0.001, dcos(d1)) * 3600 > temp1)
            }
            if (abs(d2-d1)  temp2)
            delete (temp1, verify-)

            fd2 = temp2
            if (fscan (fd2, r1, d1, data1) == EOF);
            while (fscan (fd2, r2, d2, data2) != EOF) {
                if (keepcoords)
                    printf (fmtstr, n, r1, d1, data1, >> temp3)
                else
                    printf (fmtstr, n, data1, >> temp3)
                skysep (r1, d1, r2, d2, raunit=raunit, verb-)
                if (skysep.sep > temp3)
            else
                printf (fmtstr, n, data1, >> temp3)
            n += 1

            r1 = r3; d1 = d3; data1 = data3
        }         fd1 = ""; delete (in, verify-)

        print (r1, d1, data1, >> temp1)
        if (r1 / max (0.001, dcos(d1)) * 3600 > temp1)
        }

        sort (temp1, col=1, num+, > temp2)
        delete (temp1, verify-)

        fd2 = temp2
        if (fscan (fd2, r1, d1, data1) == EOF);
        while (fscan (fd2, r2, d2, data2) != EOF) {
            if (keepcoords)
                printf (fmtstr, n, r1, d1, data1, >> temp3)
            else 
                printf (fmtstr, n, data1, >> temp3)
            skysep (r1, d1, r2, d2, raunit=raunit, verb-)
            if (skysep.sep > temp3)
        else
            printf (fmtstr, n, data1, >> temp3)

        # Now write out the lists and check for duplicate which must be
        # merged.

        sort (temp3, col=1, num+, > temp1); delete (temp3, verify-)
        touch (temp2)
        fd1 = temp1
        if (fscan (fd1, i, data1) == EOF);
        while (fscan (fd1, j, data2) != EOF) {
            if (data1 == data2) {
                print (j, i, >> temp2)
                next
            }
            printf ("%s_%03d%sn", out, i, extn) | scan (fname)
            print (data1, >> fname)
            i = j; data1 = data2
        }
        fd1 = ""; delete (temp1, verify-)
        printf ("%s_%03d%sn", out, i, extn) | scan (fname)
        print (data1, >> fname)
        sort (temp2, col=1, num+, rev+) | unique (> temp1)
        delete (temp2, verify-)

        # Merge the lists.
        n1 = n
        fd1 = temp1
        while (fscan (fd1, j, i) != EOF) {
            printf ("%s_%03d%sn", out, j, extn) | scan (fname)
            if (access (fname)) {
                printf ("%s_%03d%sn", out, i, extn) | scan (in)
                concat (in, fname, append+)
                delete (in, verify-)
                n1 -= 1
            }
        }
        fd1 = ""; delete (temp1, verify-)

        # Renumber if needed.
        if (n1 != n) {
            i = 1
            for (j=1; j out//extn)
end
New task to group data on the sky | 0 comments | Create New Account

The following comments are owned by whomever posted them. This site is not responsible for what they say.



Privacy Policy
Terms of Use

User Functions

Login