Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help Healpix: from data to map

  1. Mar 26, 2014 #1
    Hi,

    i've just installed Healpix on IDL, and i'm starting to try all the subroutines.
    First of all, i learned how to make a dipole map with random numbers (it works!:) ). Now i'd like to plot a map using real data.

    The .txt file is organized like this:

    DEC(°) AR(°) DATA(uK) ERROR(uK)

    (i don't need the ERROR column)
    First of all, i need a to-do list, to draw a fast block diagram. The Healpix manual is just a list of subroutines, and i don't know wich strategy i should use.

    Can someone help me? I'm new either in Healpix or IDL.

    Thanks in advance.
     
  2. jcsd
  3. Mar 27, 2014 #2

    Chalnoth

    User Avatar
    Science Advisor

    Does the data you got the source from offer their data in FITS format? Most full-sky astrophysical data sets do. Examples would be the Planck and WMAP data sets.

    Note that most astrophysical images of individual objects really aren't very good for Healpix format. What data are you trying to visualize?
     
  4. Mar 27, 2014 #3
    Hi Chalnoth and thank you,

    unfortunately i don't have the data in FITS format, only .txt file. The data comes from a balloon experiment for CMB measurements, so it covers only a limited sky region (no objects on particular).
    Should i convert the .txt file in FITS?
     
  5. Mar 27, 2014 #4

    Chalnoth

    User Avatar
    Science Advisor

    That's not an easy task.

    Easier, I think, to read it into an array using the flat sky approximation and display it in IDL using the TV procedure.

    If you *really* want to use Healpix, convert the coordinates on the sky to Healpix pixel indices using the following procedures:
    http://healpix.jpl.nasa.gov/html/idlnode41.htm

    (Specifically you'd want to use ang2pix_nest or ang2pix_ring, depending upon which pixelization you want to use)
     
  6. Mar 29, 2014 #5
    Unfotunately i have to use healpix, even if i haven't a full-sky map!
    So, i need to use the ring scheme. The manual (http://healpix.jpl.nasa.gov/html/intronode4.htm) says that with this scheme the Fourier transform is easy to implement (after the map i'll have to plot the angular power spectrum..).

    The routine you mentioned is:

    ang2pix_ring, nside, ipring, theta, phi

    Where:
    nside: resolution parameter.
    ipring: the pixel number whose coordinates are theta and phi.

    Now, how can i put the coordinates from the .txt file inside the theta and phi variables? Should i use two arrays called, for example, theta and phi, to store coordinates?

    Thanks!
     
  7. Mar 29, 2014 #6

    Chalnoth

    User Avatar
    Science Advisor

    Yes. How are the coordinates stored in the .txt file? What is their format?
     
  8. Mar 29, 2014 #7
    Here's an example:

    78.002930 -40.033035 84.2312 60.1127

    The order is AR, DEC, DATA (uK), dummy.

    The 4th column isn't necessary.
    I found this tutorial https://www.idlcoyote.com/tips/ascii_column_data.html
    Now i'm going to work on this. I think it's an easy task (i'm starting to understand IDL logic).

    Let me write some code and then i'll show you what i've done :)
     
  9. Mar 29, 2014 #8

    Chalnoth

    User Avatar
    Science Advisor

    To convert from RA/DEC, you can make use of NASA's Astro library, which includes the routine glactc.
     
  10. Mar 29, 2014 #9
    Is coordinate change necessary? Can i produce the map with celestial one?

    Here's my code:

    Code (Text):
    pro makeMap

    init_healpix                                

    ;From data file to 3 arrays
    OPENR, lun, 'dati.txt', /GET_LUN            
    nLines = FILE_LINES('dati.txt')            
    data = FLTARR(4, nLines)                    
    READF, lun, data                            
    AR = data(0,*)                              
    DEC = data(1,*)                            
    TEMPERATURE = data(2,*)                            

    ;Degrees -> Radians
    AR[*] = (AR[*]*!pi)/180.
    DEC[*] = (DEC[*]*!pi)/180.

    nPixel = nLines
    nSide = npix2nside(nPixel)

    ;Map creation
    ang2pix_ring, nSide, AR, DEC, map              ;map is the output (?)
    map[*] = TEMPERATURE[*]                       ;associate pixel to data

    WRITE_FITS_map,'Mappa.fits', map, /ring
    mollview,'Mappa.fits',titleplot='Mappa'

    end
    And the error: % ANG2PIX_RING: Invalid Nside: -1

    The problem is nPixel. I put nPixel = nLines, in this way each pixel has a data point...but it seems to be wrong!
    Maybe because the number of lines/data points (28121) is not a power of 2?
     
    Last edited: Mar 29, 2014
  11. Mar 29, 2014 #10

    Chalnoth

    User Avatar
    Science Advisor

    Yes. Right ascension isn't measured in degrees, but in hours. Your map will be very distorted with your code as it stands. It's really easiest to use a library somebody else has written for this, so you don't make simple errors.

    Nside does need to be a power of two.
    Npix = 12 * Nside^2

    I'd suggest starting with nside=1024. This is a reasonably high resolution. If you see a lot of empty spots in the resulting map, you can reduce the resolution.
     
  12. Mar 30, 2014 #11
    Ok, i'll make a coordinate transformation using that routine!

    Using this formula (writing before that nside = 1024), it reports this:

    Code (Text):

    WRITE_FITS_MAP: Non-Healpix data set
     npix =        28121 nside =       -1
     *** file NOT written !! ***
    Mappa.fits not found
     
    (28121 is the number of lines in the .txt file)

    Maybe i'm using the wrong routine to draw the map? I used this for the dipole, quadrupole, etc.. maps and it worked well.
     
  13. Mar 30, 2014 #12

    Chalnoth

    User Avatar
    Science Advisor

    That's wrong. Either set nside, or set npix to be 12 * nside^2. Npix should not be the number of lines in the .txt file.
     
  14. Mar 31, 2014 #13
    I set nside = 1024. Now i don't calculate npix because ang2pix_ring needs only nside to work (i also took a look to the source code).

    So, here's the new code (now i'm using degrees for coordinates), but the error remains the same:

    Code (Text):
    pro makemap

    init_healpix                                

    ;From data file to 3 arrays
    OPENR, lun, 'dati.txt', /GET_LUN            
    nLines = FILE_LINES('dati.txt')            
    data = FLTARR(4, nLines)
    READF, lun, data                            
    LAT = data(0,*)                              
    LON = data(1,*)                            
    TEMPERATURE = data(2,*)                            

    ;Degrees -> Radians
    LAT[*] = (LAT[*]*!pi)/180.
    LON[*] = (LON[*]*!pi)/180.

    ;Map creation
    nSide = 1024
    map = FLTARR(nLines)
    ang2pix_ring, nSide, LAT, LON, map              ;map is the output
    map[*] = TEMPERATURE[*]                         ;associate pixel to data
    WRITE_FITS_map,'Mappa.fits', map, /ring
    mollview,'Mappa.fits',titleplot='Mappa'

    end
     
  15. Mar 31, 2014 #14

    Chalnoth

    User Avatar
    Science Advisor

    You're never creating the map.

    The following line:

    ang2pix_ring, nSide, LAT, LON, map

    ...fills the array "map" with the indices of the pixels corresponding to the latitude and longitude. A better name for it would be "pixel indices".

    You then need to create a new map of the appropriate nSide, and fill its values with the temperature at the corresponding pixels. Does that make sense?
     
  16. Apr 2, 2014 #15
    No, that's why i decided to use the write_fits_cut4 routine, that uses the pixIndex array and the original TEMPERATURE array. Then i plotted the map with mollview (it works!!).

    Code (Text):
    ;Map creation
    nside=512
    ang2pix_ring, nside, AR, DEC, pixIndex

    WRITE_FITS_CUT4,'Mappa.fits', pixIndex, TEMPERATURE, Coordsys = 'Q', nside=nside, /ring
    MOLLVIEW,'Mappa.fits', COORD='Q'
    The problem is that the data square is located in the wrong place on the sphere (i use mollcursor to check it). It is located in [214,-53] instead of [80,-45]...i know that the coordinates are equatorial (letter Q). It's really strange, i didn't convert the coordinates, i only made a transformation from degrees to radians...
     
  17. Apr 2, 2014 #16

    Chalnoth

    User Avatar
    Science Advisor

    Oh, I see. That makes sense.

    You might want to check the latitude coordinate. Note that in degrees, latitude goes from 90 degrees at the north pole to -90 degrees at the south pole, while in radians is goes from 0 at the north pole to pi at the south pole.
     
  18. Apr 3, 2014 #17
    You're right!

    I used this formula:
    Code (Text):
    THETA = (!pi/2.)*(1.-THETA/90.)
    I solved all the problems, the map is correct!

    Thank you Chalnoth :)
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook