SOCIS Final Status Report - Cleaning

After a great summer I would like to share with you what I worked on. 😊 But before that I would like to thank everyone who made that possible. Thank you Shane that you was my mentor this summer and guided me through every problem. I definitely learnt a lot and it was fun. 😃 Thanks for the SunPy community taking me in, and for their every help and I also want to thank ESA and Maxime Perrotin who organized this great program!

And now, we can check out the results! 😉

You can find the documentation in the project wiki  and the theoretical background in my previous posts. Now I want to concentrate on the actual usage of the library in case of an actual problem from A-Z.

Lets say that we have a fits file with several RHESSI visibilities and we would like to process it. First we would like to get the xrayvision module (I assuming that SunPy is already installed).

To get the module we have to check it out from GitHub:

 git clone https://github.com/sunpy/xrayvision.git 

After that we have to install it for example with pip (I am going to ask pip to keep the data in place):

 cd xrayvision
 pip install -e  .

Now we are ready to work on some data!

Ok, you have already downloaded some data. It is there because of the tests what were also implemented beside the features. (You can even find a precommit git hook in the hooks folder if you want to run the tests before every commit.)

I am going to work with the following file:

 xrayvision/data/hsi_20020220_110600_1time_4energies.fits

First we should import the RHESSIVisibility class, numpy, matplotlib and read in our data:

 import numpy as np
 import matplotlib.pyplot as plt
 from xrayvision.Visibility import RHESSIVisibility
 rhessi_objects = RHESSIVisibility.from_fits_file("xrayvision/data/hsi_20020220_110600_1time_4energies.fits")

We have 4 RHESSIVisibility objects now in a list. I am going to work with the first one.
I would like to check out the image first, what we can get from the data. To do that we have to transform it (I am going to create a 128x128 sized image from it):

 image = rhessi_objects[0].to_map_v2(np.zeros((128,128)))  

To view it we can use the matplotlib library.

 plt.imshow(image)
 plt.show()  

The result:

We would like to use the CLEAN algorithm on this image. To be able to do that we have to import it:

 from xrayvision.Clean import Hogbom
 from xrayvision.Clean import ReasonOfStop

We would need the dirty beam. I am going to create it from our visibility data. To do that I am creating a copy of the RHESSIVisibility object, and after that I simply replace all of the visibility values with one. I have to multiply this with the visibility data what I would get watching at a Dirac-delta with our instrument.

 import copy
 psf_vis = copy.deepcopy(rhessi_objects[0])
 psf_vis.vis = np.ones(len(psf_vis.vis))
 dirac_vis = copy.deepcopy(psf_vis)
 dirac_map = np.zeros((128,128))
 dirac_map[63:64,63:64] = 1.
 dirac_vis.from_map_v2(dirac_map, center=(0,0))
 
 psf_vis.vis = np.multiply(psf_vis.vis, dirac_vis.vis)
 psf = psf_vis.to_map_v2(np.zeros((128,128)))
 plt.imshow(psf)
 plt.show()




And now we got the dirty beam:

We have everything to use CLEAN! Because it is an iterative method we can set the maximal number of the iteration or a threshold where it should stop. We also have to set the gain.

 clean = Hogbom(rhessi_objects[0], psf, 0.01, (128,128), gain=0.01, niter=1000)
 while clean.iterate(gain=False) is ReasonOfStop.NOT_FINISHED:
     pass




We can check why it finished the iteration:

 print(clean.iterate(gain=False))




In our case it says that we reached the limit of the iteration count. Ok, how can we get our final CLEANed image?

 result_image = clean.finish(3)
 plt.imshow(result_image)
 plt.show()




When I called the finish function, I gave the standard deviation in pixels for a Gaussion which is convoluted with the sky map. Our result:




As you can see when you compare it to the original image, it really found some features on the image. To get good and valid result, you have to set right parameters. What you saw in this example were just random values to demonstrate the process! Have fun! 😄