How to Plot a Map Using Healpix and Real Data?

[fineTuner]

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 which strategy i should use.

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

Thanks in advance.

 

[Chalnoth]

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?

 

[fineTuner]

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?

 

[Chalnoth]

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?

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)

 

[fineTuner]

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!

 

[Chalnoth]

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!

Yes. How are the coordinates stored in the .txt file? What is their format?

 

[fineTuner]

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 :)

 

[Chalnoth]

To convert from RA/DEC, you can make use of NASA's Astro library, which includes the routine glactc.

 

[fineTuner]

Is coordinate change necessary? Can i produce the map with celestial one?

Here's my code:

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?

 

[Chalnoth]

Is coordinate change necessary? Can i produce the map with celestial one?

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.

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?

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.

 

[fineTuner]

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.

Ok, i'll make a coordinate transformation using that routine!

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

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

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.

 

[Chalnoth]

(28121 is the number of lines in the .txt file)

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.

 

[fineTuner]

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:

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

 

[Chalnoth]

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?

 

[fineTuner]

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?

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!).

;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...

 

[Chalnoth]

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!).

;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'

Oh, I see. That makes sense.

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...

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.

 

[fineTuner]

You might want to check the latitude coordinate.

You're right!

I used this formula:

THETA = (!pi/2.)*(1.-THETA/90.)

I solved all the problems, the map is correct!

Thank you Chalnoth :)