I’ve just recently been looking into a similar problem.
I do amateur solar imaging with my Lunt50 telescope and Saturn-M camera. After initial processing using common tools for stacking and sharpening, I’ve been working on a python script to do additional processing to enhance contrast beyond simple image-wide gamma curves and histogram equalization.
A key part of that is the need to identify the center and radius of the solar disk in the image, so that’s what I first worked on using opencv. Take a look at my script:
Look at the findValidCircle and findCircle functions, using the EDCircles algorithm. It does a good job of finding the center and radius.
Recently, I’ve been looking at SunPy and thinking how I might use it. What strikes me is that it seems very focused on data from professional instruments where all of the relevant metadata is available. For someone like me with my personal Lunt50 and you with your PST, there’s little meta data available. I was looking at this documentation for how to create a “map” from an image file:
And it seems like what I need to do is call make_fitswcs_header(), and give it info like date/time, center reference pixel, and scale of my pixels.
My circle-finding script can find the center. The scale in arcseconds can be estimated using the radius in pixels and the known size of the sun in arcseconds for the day of the image capture. Or can be determined directly if you know the telescope’s optical parameters and the camera’s pixel size. Though that might be complicated if capturing through an eyepiece.
You can use a tool to convert image file formats, but I think the important thing to understand is that your tiff files are missing important meta data that sunpy requires, so the fits files will necessarily be missing it as well. If the meta data was all there in the fits file, SunPy would automatically extract it and use it, but that’s not really an option. So instead, I think best is to piece together what info you can and feed it into make_fitswcs_header().
If I understand this bit from the example, it essentially says that the reference pixel is the center of the sun, and it’s captured from earth on a given date/time - so aside from date/time, this would be the same for every image:
coord = SkyCoord(0u.arcsec, 0u.arcsec, obstime=‘2013-10-28 08:24’,
observer=‘earth’, frame=frames.Helioprojective)
This piece then would get the specific center found and calculated pixel scale:
header = sunpy.map.make_fitswcs_header(data, coord,
reference_pixel=[0, 0]*u.pixel,
scale=[2, 2]u.arcsec/u.pixel,
telescope=‘Fake Telescope’, instrument=‘UV detector’,
wavelength=1000u.angstrom)
Probably unnecessary, but you could set telescope to the PST, instrument to your camera, and wavelength to 6562.8 angstroms for h-alpha.
Anyway, this is the sort of thing I was looking into… having a somewhat automated function to fill in this data. If there was a SunPy function that gave the size of the sun in arcseconds for a given day, that would make things simpler.