SOHO/MDI

Example of Basic MDI Data Manipulation Using IDL

This example demonstrates basic IDL manipulation of an MDI intensitygram. It illustrates reading in a file; rotating, cropping, scaling and rebinning it; writing output to FITS and gif files; examining a single row of pixels; and computing the average of a subset of pixels. Some of the IDL procedures used are from SolarSoft. The commands are also combined in one procedure file with comments: MDI_example.pro.

Read in the FITS image and header of an MDI intensitygram. The /scale option converts to physical units using image header keywords BSCALE and BZERO:
```
    data = rfits('fd_Ic_6h_01d.4599.0000.fits',head=hd,/scale)
```

From the image header, read in the values for observation time, p angle, solar radius, and image center. Shorten the observations timestamp to leave off seconds and '_TAI':


    obs_time = strtrim(sxpar(hd, 'T_OBS'),2)
    pangle = sxpar(hd, 'P_ANGLE')
    radius = sxpar(hd,'R_SUN')
    x0 = sxpar(hd, 'X0')
    y0 = sxpar(hd, 'Y0')
    dattim = strmid(obs_time,0,16)
    print, '  TIME:  ', dattim, '     P_ANGLE: ', pangle

  TIME:  2005.08.05_00:00     P_ANGLE:  -179.99999

Crop the image by setting values not on the solar disk to 0 and rotate the image so north is pointing up. Write the result to a fits file:


    cropped = circle_mask(data, x0, y0, 'GE', radius, mask=0)
    rotated = rot(cropped, pangle, 1, x0, y0)
    
    writefits, 'fd_Ic_cr_rot.4599.0000.fits', rotated

Display the result to find the position of a sunspot (or use a FITS reader to examine the image):
```
IDL> xtv, rotated(*,*)
```

Look at a sunspot profile by reading the pixels along the x-axis that goes through the lower sunspot in the image:


    sunspot_profile = intarr(2,936)
    for x = 45, 980 do begin
       i = x - 45
       sunspot_profile(0,i) = x
       sunspot_profile(1,i) = rotated(x,366)
    end

Plot the result and write it to a gif file:


    pix = sunspot_profile(0,*)
    int = sunspot_profile(1,*)
    
    title = 'Sunspot Profile'
    ylabel = 'Intensity'
    xlabel = 'Pixel'
    plot, pix, int, title=title, xtitle=xlabel, ytitle = ylabel
    write_gif, 'sunspot_profile.gif', tvrd()

Rescale and rebin the cropped, rotated image and write it to a gif file for easier displaying:


    scaled = bytscl(rotated, min=1000, max=14000)
    im = rebin(scaled,512,512)
    tv2, 512, 512, /init
    tv2, im(*,*)
    title = 'MDI Intensitygram:  ' + dattim
    xyouts2, 100, 500, title, charsize=1.5
    write_gif, 'im.gif', tvrd()

Calculate the average of the central 100 pixels of the rebinned, rescaled image:


    total = 0
    for x = 251, 260  do begin
       for y = 251, 260  do begin
         total = total + im(x,y)
       endfor
    endfor
    avg = total / 100
    print, 'AVERAGE:  ', avg

    end

AVERAGE:       228

Oct 2005