cd ~/IDL rm -rf Demos rm -rf Data wget -r -nc -nH --cut-dirs=1 http://ross.iasfbo.inaf.it/IDL/Demos wget -r -nc -nH --cut-dirs=1 http://ross.iasfbo.inaf.it/IDL/Data wget http://ross.iasfbo.inaf.it/IDL/healpix201.tar.gz wget http://ross.iasfbo.inaf.it/IDL/jhuapl_sep07.tar.gz wget http://ross.iasfbo.inaf.it/IDL/solarpath.tar.gz tar zxvf healpix201.tar.gz tar zxvf jhuapl_sep07.tar.gz tar zxvf solarpath.tar.gz cd Demos
Reading a columns formatted ASCII file:
We assume the data file demodat.asc is in the dir ~/IDL/Data/
; Let's build the full path to the file
file = FILEPATH('demodat.asc',SUBDIR='IDL/Data',ROOT=GETENV('HOME'))
a = READ_ASCII(file)
; Let's see what's in a (it is a "structure" so I use the STRUCTURE kwd)
HELP, a, /STRU
; READ_ASCII is very flexible but you need to get familiar with it.
; Can get help on READ_ASCII using "doc_library"
DOC_LIBRARY, 'read_ascii'
; Let's use some keywords:
; we ask for data from line "23",
; put the first 22 "header" lines into h
; return the number of data lines read into n
a = READ_ASCII(file, DATA_START=22, HEADER=h, COUNT=n)
PRINT, 'There are '+ STRTRIM(n,1) +' data records.'
PRINT, h ; all the header as a single string
PRINT, h, FORMAT='(A)' ; formatted output, 1 line per read row
PRINT, a.field1[*,0:3] ; the first 3 rows (all columns) of "tag" field1
; To set the template structure of a file, use ASCII_TEMPLATE.
tmpl = ASCII_TEMPLATE(file)
; Data start at line 23. Can then simply click "Next", "Next", "Finish" or
; change the column names as you like. Data types are what matters.
; See how the returned structure looks like
HELP, tmpl, /STRU
PRINT, tmpl.fieldnames
; Now we can read the file using this template and see what is changed
a = READ_ASCII(file, TEMPLATE=tmpl, HEADER=h)
HELP, a, /STRU
; You'll see 7 distinct tags corresponding to 7 FLOAT arrays with the columns
; data values!
; Let's print column 3
PRINT, a.field3
For "simple" cases, e.g. want to map all the columns into a 2-d FLOAT or
DOUBLE array, it would be handy to use my personal routine nxyreadf:
; I use the keyword "ncol" to force reading all the data. The def. is 2 columns a = nxyreadf(file,ncol=7) HELP, a ; Plot 2nd column vs first column PLOT, a[0,*], a[1,*] ; Plot 3rd column vs first column PLOT, a[0,*], a[2,*]
; Read a SAX-PDS light curve of the Crab pulsar. Bins are "pulse phase".
a = lcqdpread('crab_pdspulse.qdp', /phase)
a[2,*] = a[2,*]/MAX(a[2,*]) ; normalise in the range 0 - 1
plot, [0], /NODATA, /NOERASE, xsty=1, ysty=1, chars=1.5, $
xrange=[-0.30,0.80], yrange=[0.5,1.05], $
xtitle='!17Pulse Phase!6', ytitle='!17Normalized cts/s!6'
binplot, a[0,*]-1, a[2,*], thick=2, /connect
; Read the spectra indices files. Data also have errors, so use the
; personal routine erplot to plot them.
a = lcqdpread('crab_pdsphin.qdp', /phase)
plot, [0], /NODATA, /NOERASE, xsty=1, ysty=1, chars=1.5, $
xrange=[-0.10,0.50], yrange=[1.5,2.2], $
xtitle='!17Pulse Phase!6', ytitle='!17Photon index'
binplot, a[0,*], a[2,*], a[1,*], thick=2
erplot, a[0,*], a[2,*], a[3,*]
; Run the demo program crabpulse which plots light curves and
; phase spectral indices for the four SAX instruments
crabpulse
; To get a PS version simply add /PS
crabpulse, /ps
x = FINDGEN(16)*10 & ex = SQRT(x) & y=FINDGEN(16)*100 & ey=SQRT(y) PLOT, x, y, PSYM=4, SYMSIZE=1.5 erplot, x, y, ey & binplot, x, y, exIf not done yet, check the Wed21 section on how to read and write FITS files.
An example of FITS image manipulation:
Here we display 2 images around the stars WR45 and Feige56 taken with an optical camera. Sky coordinates grid is overplotted. The images are smoothed using FFTs and re-displayed in shaded surface. As usual, images are in the "IDL/Data" directory.
fits_read, '../Data/IMG122R126.fits.gz', image, h
EXTAST, h, astrom
HELP, astrom, /STRU
; Let's get the image size: we use SIZE: first element is Ndim, 2nd, 3rd Dims
s = SIZE(image)
PRINT, s
nx = s[1] & ny = s[2]
; Let's make an empty plot area where to show the image.
; Note that we leave 2 extra pixels for the axis
PLOT, [0], XRANGE=[0,nx+1], YRANGE=[0,ny+1], CHARS=1.3, xsty=5, ysty=5
; Now get the plot area in DEVICE coordinates (from DATA, the default)
cv = CONVERT_COORD([0,nx], [0,ny], /TO_DEV)
; Byte scale using a reasonable Z range (can try other ranges)
imin = MIN(image, MAX=imax)
bima = BYTSCL(image, MIN=imin, MAX=imin*1.2)
; Display as 2-d pseudo image using color table 3
loadct, 3
TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
; Now overplot the sky frame. Note the keywords used.
skyframe, '../Data/IMG122R126.fits.gz', /dogrid, gridcol=255, CHARS=1.3, $
XTITLE='RA', YTITLE='Dec', /noera
; The image is "square" but the plot area is not? Use sqplset.
; The second parameter is the CHARacter Size we'll use. The returned
; parameters are passed to PLOT via XMARGIN and YMARGIN
sqplset, 1.3, xm, ym
PLOT, [0], XRANGE=[0,nx+1],YRANGE=[0,ny+1], CHARS=1.3, xsty=5, ysty=5, $
XMARGIN=xm, YMARGIN=ym
; Now just re-convert the DATA coordinates to DEVICE coordinates.
; Note how we use XMARGIN and YMARGIN also for "skyframe"
cv = CONVERT_COORD([0,nx],[0,ny],/TO_DEV)
TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
skyframe, '../Data/IMG122R126.fits.gz', /dogrid, gridcol=255, CHARS=1.3, $
XTITLE='RA', YTITLE='Dec', /noera, XMARGIN=xm, YMARGIN=ym
; You'll probably see that the "Dec" label is not visible. Try to understand
; why and find a solution!
; Let's have a better view in a zoom window. Simply type:
zoom
; Pay attention to the message! Use the correct mouse buttons!
; If you want a "continuous" zoom use the CONTINUOUS keyword. Type:
zoom, /CONT
; Have a look to the other options with the help
doc_library, 'zoom'
; Let's get the cordinates at a given position. First get position using CURSOR
CURSOR, x, y, /DATA
; Then convert into sky coordinates (fractional dregrees) using the astrometry
; info in the "astrom" structure
XY2AD, x, y, astrom, xd, yd
PRINT, 'RA=', xd, ' Dec=', yd
; Then convert into string format
coord_s = ADSTRING(xd, yd, 1)
PRINT, 'Coords:', coord_s
; Smooth with a 5 pixels FWHM gaussian using FFT
simage = filter_image(image, FWHM_GAUSSIAN=5)
; Rescale image and replot
imin = MIN(simage, MAX=imax)
bima = BYTSCL(simage, MIN=imin, MAX=imin*1.2)
TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
skyframe, '../Data/IMG122R126.fits.gz', /dogrid, gridcol=255, CHARS=1.3, $
XTITLE='RA', YTITLE='Dec', /noera, XMARGIN=xm, YMARGIN=ym
; As a shaded surface (giving some ranges)
SHADE_SURF, simage, XSTY=1, YSTY=1, XRANGE=[0,1024], YRANGE=[0,1024]
; The base of the main star
SHADE_SURF, simage[441:530,461:530]
; Surface plot in "LEGO" format
SURFACE, simage[441:530,461:530], /LEGO
; Colored mesh using the "Rainbow + white" table (exclude index 255 = white!)
simage_c = simage[441:530,461:530]
imin = MIN(simage_c, MAX=imax)
shad = (simage_c - imin)/(imax - imin) * 254
loadct, 39
SURFACE, simage_c, SHADES=shad
; or
SURFACE, simage_c, SHADES=shad, /LEGO
Let's do a similar work with a slitless optical image, but within a
PROgram.
Edit a file in the Demos dir. with your preferred text editor, for example:
$ xemacs image_test.proThen type in the following code (cut-and-paste):
PRO image_test, filename
IF N_PARAMS() NE 1 THEN $
file = '../Data/IMG039A119.fits.gz' $
ELSE $
file = '../Data/' + filename
fits_read, file, image, h
EXTAST, h, astrom, ptype
IF (ptype GT 0) THEN $
PRINT, 'The image has valid astrometry keywords.' $
ELSE BEGIN
MESSAGE, 'The image has NOT valid astrometry keywords.'
ENDELSE
s = SIZE(image)
nx = s[1] & ny = s[2]
PRINT, 'This is a '+ STRTRIM(nx,2) +'x'+ STRTRIM(ny,2) +' image.'
loadct, 5
WINDOW, 1, XSIZE=700, YSIZE=600
sqplset, 1.3, xm, ym
PLOT, [0], XRANGE=[0,nx+1],YRANGE=[0,ny+1], CHARS=1.3, xsty=5, ysty=5, $
XMARGIN=xm, YMARGIN=ym
cv = CONVERT_COORD([0,nx],[0,ny],/TO_DEV)
imin = MIN(image, MAX=imax)
bima = BYTSCL(image, MIN=imin, MAX=imin*1.2)
TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
skyframe, file, /dogrid, gridcol=255, CHARS=1.3, $
XTITLE='RA', YTITLE='Dec', /noera, XMARGIN=xm, YMARGIN=ym
simage = filter_image(image, FWHM_GAUSSIAN=5)
imin = MIN(simage, MAX=imax)
bima = BYTSCL(simage, MIN=imin, MAX=imin*1.2)
WINDOW, 2, XSIZE=700, YSIZE=600
TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
skyframe, file, /dogrid, gridcol=255, CHARS=1.3, $
XTITLE='RA', YTITLE='Dec', /noera, XMARGIN=xm, YMARGIN=ym
WINDOW, 3, XSIZE=700, YSIZE=600
SHADE_SURF, simage, XSTY=1, YSTY=1, XRANGE=[0,1024], YRANGE=[0,1024], $
ZTITLE = 'WR45: Smoothed image'
END
; Then can run it image_test ; Note that you could first try to compile it .r image_test ; As you note, the file name can be given as parameter. Try with the previous ; image name (omit "../Data")
Aitoff (galactic) projection of celestial sources:
Se the aitoff routines used below.
; AGNs and unidentified EGRET (3EG) sources
cieloagn
See cieloagn.pro source code (in the Demos dir.).
; Produce a PS file
cieloagn, /PS
Have a look to the used functions and routines!
HEALPix maps:
Plot the astronomical catalogue UCAC-2 as objects density in orthographic and mollview all-sky maps + local orthographic view. Here the sky resolution is about 14 arcmin. The "ucac2_cover_k8.fits" file was produceded reading directly from a MySQL database.
orthview, '../Data/ucac2_cover_k8.fits', rot=[45,45], /grat, col=5, $
title='UCAC 2 - objects density'
mollview, '../Data/ucac2_cover_k8.fits', coord=['C','G'], /grat, col=5, $
title='UCAC 2 - objects density'
; 12 degrees around the galactic center in gnomic view (graticule 1 deg.):
gnomview, '../Data/ucac2_cover_k8.fits', coord=['C','G'], res=10, grat=[6,6], $
title='UCAC 2 - Galactic Center'
; To produce a PNG output file use the PNG keyword. For example:
orthview, '../Data/ucac2_cover_k8.fits', rot=[45,45], /grat, col=5, $
title='UCAC 2 - objects density', png='ucac2_density_orth.png'
Graphics AstroTools:
The Astronomy Users' library and other libraries (e.g. the JHU/APL/S1R Johns Hopkins Applied Physics Laboratory S1R group library) have plenty of Astronomical routines and tools. Please check routines description (e.g. at idlastro.gsfc.nasa.gov and fermi.jhuapl.edu). Here we only show a few graphics tools.
The Javier Corripio's tool to compute and plot solar path for a given day, year, longitude and latitude. Source filea are in the Solarpath dir.
; Run the widget based user interface
plotsun
; Input local latitude and longitude: for Bologna it is lat.= 44.5 lon.= 11.6
; The showed data is also written into an output file named "solpos.dat"
Position of the the Sun projected on the Earth map with day, twilight, and night
zones. Uses routines in the Jhuapl directory.
; Simply run the program sunclock ; or in Orthographic projection sunclock, /ORTH ; or in Mercator projection sunclock, /MERC
; Make a square area in a 600x600 graphics window
WINDOW, xs=600, ys=600
sqplset, 1.0, xm, ym
PLOT, [-130,130], [-130,130], XSTY=1, YSTY=1, XMARGIN=xm, YMARGIN=ym, /NODATA
; Make the grid, centered on the South pole
eqpole_grid, /label, /south
; If you don't want to see the axes, use "XSTY=5, YSTY=5" in the PLOT command
; Longitude/Latidude or RA/Dec are in degrees
glon = [15, 30, 45, 60]
glat = [0, -20, -30, -60]
; Convert into cartesian coordinates for plotting.
; The keyword "south" tells we give -ve Lat
eqpole, glon, glat, x, y, /south
c = loadct19()
plots, x, y, psym=2, syms=2, col=c[2]
Put some objects onto an Aitoff projection map (Galactic coordinates).
WINDOW, xs=700, ys=400 ; Let's assume we have 5 objects with Euqtorial coordinates (RA/Dec) ra = [2.56, 50.90, 80.14, 259.30, 340.26] dec = [73.17, 51.37, 25.75, -27.63, -67.60] : Let's convert them into Galactic coordinates glactc, ra, dec, 2000., gl, gb, 1, /DEG
; Convert into cartesian coordinates for plotting. Use a filled circle marker
aitoff, gl, gb, x, y
pp = FINDGEN(32)*(!pi*2/31.)
USERSYM, COS(pp), SIN(pp), /FILL
PLOT, x, y, XRANGE=[-180,180], YRANGE=[-90,90], XSTY=5, YSTY=5, $
PSYM=8, COLOR=c[2]
; The Aitoff grid with marked corner values
aitoff_grid2, 30, 15, label=1
Routines involving a Sphere. Use the spherical plot (sph*) routines
in the Jhuapl directory. See the S1R pages.
; Plot the great circle between two points on Earth. From Bologna to La Silla ; Observatory (Chile) lat1 = 44.5 lng1 = 11.6 lat2 = -29.25 lng2 = -70.73 MAP_SET, /CONT ; Plot Earth map ll2rb, lng1, lat1, lng2, lat2, r, b ; Find r and b of pt2 rb2ll, lng1, lat1, maken(0,r,100), b, lng, lat ; Step along track PLOTS, lng, lat, col=c[3] PLOTS, lng1, lat1, psym=2, syms=2, col=c[2] ; Plot point 1 PLOTS, lng2, lat2, psym=2, syms=2, col=c[4] ; Plot point 2
; To view the plot in Orthographic projection use: ERASE, 'ffffff'x MAP_SET, 0, -20, /ORTH, /HORIZON, /ISOTROPIC, /CONT PLOTS, lng, lat, col=c[3] PLOTS, lng1, lat1, psym=2, syms=2, col=c[2] PLOTS, lng2, lat2, psym=2, syms=2, col=c[4] ; Note the center of the projection is Lat= 0, Long= -20 ; If you want a projection "without" the continents coastlines, simply ; omit the /CONT keyword
; Spirals around a sphere window, xs=500, ys=500 psav = !P ; save current !P structure !P.POSITION = [0,0,1,1] ; the whole graph window set_isoxy, -1.1, 1.1, -1.1, 1.1 LOADCT, 4 ERASE, 80 ; use color index 80 for the background sphinit, long=-30, lat=40, rad=.4, fill=140, pa=30 lng = makex(0,90,2) lat = lng*0 rad = maken(0.9, 0.4, n_elements(lat)) FOR a = 0, 360, 10 DO sphplot, lng+a, lat, rad, maxrad=.4 !P = psav ; restore !P structure
; Simply run the program demo_sat ...
Plot 2-D data in 3-D space (to be continued...):
See spaceplot.pro source code (in the Demos dir.).
The relevant keywords are /SAVE given to SURFACE
and /T3D used
in all the PLOT, OPLOT, PLOTS commands. The command T3D
allows you to change the plane in the space where to plot your data.
Note that here you could have problems viewing correctly the plot due to the
transformation done by demo_sat. Perhaps it's better restarting
IDL (this is temporary, I'll suggest what to do later).
; Simply execute the program spaceplot ; Use the personal widgets routine "setaspect" to change viewing angle setaspect ; Replot spaceplot
3-D data in 3-D space (to be continued...):
To plot a function of two variables (e.g., Z=f(X, Y)) inside a 3-D box can either make custom programs using commands allowing the T3D keyword, or use the IDL library procedure PLOT_3DBOX. Here is an example copied form the on-line help pages:
; Create some data to be plotted:
X = REPLICATE(5., 10.)
X1 = COS(FINDGEN(36)*10.*!DTOR)*2.+5.
X = [X, X1, X]
Y = FINDGEN(56)
Z = REPLICATE(5., 10)
Z1 = SIN(FINDGEN(36)*10.*!DTOR)*2.+5.
Z = [Z, Z1, Z]
; Create the box plot with data projected on all of the walls. The
; PSYM value of -4 plots the data as diamonds connected by lines:
PLOT_3DBOX, X, Y, Z, /XY_PLANE, /YZ_PLANE, /XZ_PLANE, $
/SOLID_WALLS, GRIDSTYLE=1, XYSTYLE=3, XZSTYLE=4, $
YZSTYLE=5, AZ=40, TITLE='Example Plot Box', $
XTITLE='X Coordinate', YTITLE='Y Coodinate', $
ZTITLE='Z Coordinate', SUBTITLE='Sub Title', $
/YSTYLE, ZRANGE=[0,10], XRANGE=[0,10], $
PSYM=-4, CHARSIZE=1.6