search

sunpy / sunpy-soar

A sunpy plugin for accessing data in the Solar Orbiter Archive (SOAR).
https://docs.sunpy.org/projects/soar/
BSD 2-Clause "Simplified" License
17 stars 11 forks source link

Query SOAR FoV Metadata for L2 Remote sensing instruments EUI, PHI and SPICE #141

Open JonCook opened 1 month ago

JonCook commented 1 month ago

Describe the feature

Recently the solar orbiter archive now has FoV metadata information available for all files via TAP for several remote sensing instruments and sensors:

Detailed information about the columns including descriptions is available via the TAP tables page published at: https://soar.esac.esa.int/soar-sl-tap/tap/tables

Using standard metadata requests it is possible to retrieve FoV information about a given file or set of files: For example for a given file:

But of course it is also possible to query on any column in the table. There are a few indexes in the columns

With this information it makes it 'quite easy' to make beautiful plots without the knowing the details and internal metadata of each instruments FITS files like:

image

As well as earth perspective: image

And positional: image

It could be nice to expose some of this functionality from the soar-sunpy package and also create some nice examples in the galleries

Proposed solution

Consider exposing some of this functionality via the soar-sunpy package and even creating a nice gallery

hayesla commented 1 month ago

This is great @JonCook ! thanks for advertising that these tables are now available.

nabobalis commented 1 month ago

We should make this a gallery example. @NucleonGodX can we add that to the gsoc list?

JonCook commented 1 month ago

If it helps @nabobalis - I have some code which makes the first plot using TAP requests. Perhaps it could be useful to get started. I'm sure it can be cleaned up a lot

nabobalis commented 1 month ago

That would be great! Thank you @JonCook

NucleonGodX commented 1 month ago

We should make this a gallery example. @NucleonGodX can we add that to the gsoc list?

Alright, lets do that.

JonCook commented 1 month ago

That would be great! Thank you @JonCook

Here it is, it is a bit 'rough', but this is the python to produce the first image in the example using TAP requests.. It also assumes you have the FITS files in question downloaded.

#!/usr/bin/env python
# coding: utf-8

# In[2]:

#Set up our dependencies

import sunpy.map
from sunpy.coordinates import frames, get_earth, get_horizons_coord, get_body_heliographic_stonyhurst
from sunpy.net import Fido, attrs as a 
import sunpy_soar
import sunraster.instr.spice
from astropy import units as u 
#from astroquery.utils.tap.core import TapPlus

from astropy.io import fits
from astropy.coordinates import SkyCoord, SkyOffsetFrame
from astropy.time import Time

import pyvo as vo

import matplotlib.pyplot as plt
import math

from dateutil import parser as date_parser
import warnings

warnings.filterwarnings('ignore')

eui_fsi_filename = '/Users/jcook/Downloads/fov-test-files-notebook/solo_L2_eui-fsi174-image_20221024T190050195_V01.fits'
eui_hri_filename = '/Users/jcook/Downloads/fov-test-files-notebook/solo_L2_eui-hrieuv174-image_20221024T190000171_V01.fits'
phi_hrt_filename = '/Users/jcook/Downloads/fov-test-files-notebook/solo_L2_phi-hrt-blos_20221024T191503_V01.fits'
spice_filename   = '/Users/jcook/Downloads/fov-test-files-notebook/solo_L2_spice-n-ras_20221024T190134_V05_150995395-058.fits'

# Lets get some stuff from TAP
service = vo.dal.TAPService("https://soar.esac.esa.int/soar-sl-tap/tap")
fits_item_resulset = service.search("SELECT * FROM soar.v_eui_sc_fits WHERE filename='solo_L2_eui-fsi174-image_20221024T190050195_V01.fits'")
#fits_item_resulset = service.search("SELECT * FROM soar.v_eui_sc_fits")

#print(fits_item_resulset)
print(fits_item_resulset.table.columns)

first_item = fits_item_resulset[0]
fsi_obs_date = first_item["date_average"]
print(first_item["soop_name"])
print(type(first_item["soop_name"]))
print(type(first_item["wavemax"]))

fov_eui_resultset = service.search("SELECT * FROM soar.v_eui_hri_fov WHERE filename='solo_L2_eui-hrieuv174-image_20221024T190000171_V01.fits'")
fov_phi_resultset = service.search("SELECT * FROM soar.v_phi_hrt_fov WHERE filename='solo_L2_phi-hrt-blos_20221024T191503_V01.fits'")
fov_spice_resultset = service.search("SELECT * FROM soar.v_spice_fov WHERE filename='solo_L2_spice-n-ras_20221024T190134_V05_150995395-058.fits'")

# Make a query to get FoV information for a given file
service = vo.dal.TAPService("https://soar.esac.esa.int/soar-sl-tap/tap")
fov_eui_resultset = service.search("SELECT * FROM soar.v_eui_hri_fov WHERE filename='solo_L2_eui-hrieuv174-image_20221024T190000171_V01.fits'")

# Print out the DALResultsTable object (it includes the correct type already)
#print(fov_eui_resultset)
# Print out the TableColumns
#print(fov_eui_resultset.table.columns)

# Print out the soop name and some column types
row = fov_eui_resultset[0]
print(row["soop_name"])
print(type(row["fov_earth_bot_left_arcsec_tx"]))
print(row["fov_earth_bot_left_arcsec_tx"])

eui_fsi_map = sunpy.map.Map(eui_fsi_filename)
eui_hri_map = sunpy.map.Map(eui_hri_filename)
phi_hrt_map = sunpy.map.Map(phi_hrt_filename)

# Generation of SunPy Map for a zoomed version of the FSI image
m_fsi_zoom = eui_fsi_map.submap(SkyCoord(-3000,-3000, unit='arcsec', frame=eui_fsi_map.coordinate_frame),\
                                     top_right=SkyCoord(3000, 3000, unit='arcsec', \
                                     frame=eui_fsi_map.coordinate_frame))

fig = plt.figure(figsize=(20, 10))
ax2 = fig.add_subplot(1, 3, 2, projection=m_fsi_zoom)
## ORBITAL plot
#fsi_obs_date=eui_fsi_info[1].header['DATE-AVG']
m_fsi_zoom.plot(axes=ax2, title='Solar Orbiter perspective (EUI 174 $\mathrm{\AA}$)\n''%s' % fsi_obs_date)

original_bl = SkyCoord(
    fov_eui_resultset["fov_solo_bot_left_arcsec_tx"]*u.arcsec, fov_eui_resultset["fov_solo_bot_left_arcsec_ty"]*u.arcsec, frame=eui_hri_map.coordinate_frame)
original_tr = SkyCoord(
    fov_eui_resultset["fov_solo_top_right_arcsec_tx"]*u.arcsec, fov_eui_resultset["fov_solo_top_right_arcsec_ty"]*u.arcsec, frame=eui_hri_map.coordinate_frame)

rotated_frame = SkyOffsetFrame(origin=eui_hri_map.reference_coordinate,
                                rotation=-fov_eui_resultset["crota"]*u.deg,
                                observer=eui_hri_map.observer_coordinate)

rotated_bl = original_bl.transform_to(rotated_frame)
rotated_tr = original_tr.transform_to(rotated_frame)

# Draw the quadrangle
m_fsi_zoom.draw_quadrangle(
    rotated_bl[0],
    top_right=rotated_tr[0],
    label='EUI HRI',
    edgecolor='C0',
    lw=2)

original_bl = SkyCoord(
    fov_phi_resultset["fov_solo_bot_left_arcsec_tx"]*u.arcsec, fov_phi_resultset["fov_solo_bot_left_arcsec_ty"]*u.arcsec, frame=phi_hrt_map.coordinate_frame)
original_tr = SkyCoord(
    fov_phi_resultset["fov_solo_top_right_arcsec_tx"]*u.arcsec, fov_phi_resultset["fov_solo_top_right_arcsec_ty"]*u.arcsec, frame=phi_hrt_map.coordinate_frame)
rotated_frame = SkyOffsetFrame(origin=phi_hrt_map.reference_coordinate,
                                rotation=-fov_phi_resultset["crota"]*u.deg,
                                observer=phi_hrt_map.observer_coordinate)
rotated_bl = original_bl.transform_to(rotated_frame)
rotated_tr = original_tr.transform_to(rotated_frame)
m_fsi_zoom.draw_quadrangle(
    rotated_bl[0],
    top_right=rotated_tr[0],
    label='PHI HRT',
    edgecolor='C1',
    lw=2)

raster_spice = sunraster.instr.spice.read_spice_l2_fits(spice_filename)
spice_window = raster_spice['O III 703 / Mg IX 706 (Merged)']
# The spice window has data in t, lambda, y, x
spice_map = sunpy.map.Map(spice_window[0, 20, :, :].data, \
                          spice_window[0, 20, :, :].meta)

original_bl = SkyCoord(
    fov_spice_resultset["fov_solo_bot_left_arcsec_tx"]*u.arcsec, fov_spice_resultset["fov_solo_bot_left_arcsec_ty"]*u.arcsec, frame=spice_map.coordinate_frame)
original_tr = SkyCoord(
    fov_spice_resultset["fov_solo_top_right_arcsec_tx"]*u.arcsec, fov_spice_resultset["fov_solo_top_right_arcsec_ty"]*u.arcsec, frame=spice_map.coordinate_frame)
rotated_frame = SkyOffsetFrame(origin=spice_map.reference_coordinate,
                                rotation=-fov_spice_resultset["crota"]*u.deg,
                                observer=spice_map.observer_coordinate)
rotated_bl = original_bl.transform_to(rotated_frame)
rotated_tr = original_tr.transform_to(rotated_frame)
m_fsi_zoom.draw_quadrangle(
    rotated_bl[0],
    top_right=rotated_tr[0],
    label='SPICE',
    edgecolor='C2',
    lw=2)

ax2.legend()

Thanks here to @hayesla and others. I just changed a bit here where we get the values from. We also have the code for generating the earth view and positional images, but I did not update it yet to work with TAP

Thanks