Skip to content

Commit

Permalink
Animate gallery example (#323)
Browse files Browse the repository at this point in the history 
Fixes #113

---------

Co-authored-by: Nabil Freij <[email protected]>
  • Loading branch information
@wtbarnes @nabobalis
wtbarnes and nabobalis authored Dec 18, 2023
1 parent 86b59e2 commit 5e50365
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 104 deletions.
1 change: 1 addition & 0 deletions changelog/323.doc.rst
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Improved the "Requesting specific AIA images from the JSOC" Example to animate and use Fido instead of DRMS.
144 changes: 40 additions & 104 deletions examples/download_specific_data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,128 +5,54 @@
This example shows how to request a specific series of AIA images from the JSOC.
We will be filtering the data we require by keywords and requesting short exposure images from a recent flare.
Unfortunately, this can not be done using the sunpy downloader `~sunpy.net.Fido`
and instead we will use the `drms <https://docs.sunpy.org/projects/drms/en/stable/>`__ Python library directly.
We will be filtering the data we require by keywords and requesting short exposure images from a recent (at time of writing) flare.
"""

import os
from pathlib import Path

import astropy.units as u
import drms
import matplotlib.pyplot as plt
import sunpy.map
from astropy.coordinates import SkyCoord
from astropy.visualization import AsinhStretch, ImageNormalize
from sunpy.net import Fido, attrs

from aiapy.calibrate import correct_degradation, register, update_pointing
from aiapy.calibrate.util import get_correction_table, get_pointing_table

#####################################################
# Exporting data from the JSOC requires registering your
# email first. Please replace this with your email
# address once you have registered.
# See `this page <http://jsoc.stanford.edu/ajax/register_email.html>`__
# for more details.
# email first. Please replace the text after the ``=``
# with your email address once you have registered.
# `See this page for more details. <http://jsoc.stanford.edu/ajax/register_email.html>`__

jsoc_email = os.environ.get("JSOC_EMAIL")

#####################################################
# Our goal is to request data of a recent (of time of writing)
# X-class flare. However, we will request the explanation of
# the keywords we want from the JSOC.

client = drms.Client(email=jsoc_email)
keys = ["EXPTIME", "QUALITY", "T_OBS", "T_REC", "WAVELNTH"]

print("Querying series info")
# We plan to only use the EUV 12s data for this example.
series_info = client.info("aia.lev1_euv_12s")
for key in keys:
linkinfo = series_info.keywords.loc[key].linkinfo
note_str = series_info.keywords.loc[key].note
print(f"{key:>10} : {note_str}")

#####################################################
# We will construct the query. The X-class flare occurred
# on the 2021/07/03 at 14:30:00 UTC. We will focus on the 5 minutes
# before and after this time.

qstr = "aia.lev1_euv_12s[2021-07-03T14:25:00Z-2021-07-03T14:35:00Z]"
print(f"Querying data -> {qstr}")
results = client.query(qstr, key=keys)
print(f"{len(results)} records retrieved.")

#####################################################
# As you can see from the output, we have received a
# a list of AIA images that were taken during the flare.
# What we want to do now is to filter the list of images
# to only include shorter expsoures.
# However, before we do this, let us check what the exposure times are.

# Filter out entries with EXPTIME > 2 seconds
results = results[results.EXPTIME < 2]
print(results)

#####################################################
# This style of filtering can be done to any column
# in the results. For example, we can filter the WAVELNTH
# column to only include 171 data with short expsoures.

# Only use entries with WAVELNTH == 211
results = results[results.WAVELNTH == 211]
print(results)

#####################################################
# .. note::
# **Only complete searches can be downloaded from JSOC**,
# this means that no slicing operations performed on the results object
# will affect the number of files downloaded.
#
# We can filter and do analysis on the metadata that was returned.
# The issue is is that if we only want this data, you can not use
# this "filtered results" to download only the data we want.
# To do this, we will have to do a second query to the JSOC,
# this time using the query string syntax the
# `lookdata <http://jsoc.stanford.edu/ajax/lookdata.html>`__ web page.
# You can use the website to validate the string before you export the query.

updated_qstr = "aia.lev1_euv_12s[2021-07-03T14:25:00Z-2021-07-03T14:35:00Z][? EXPTIME<2.0 AND WAVELNTH=211 ?]{image}"
print(f"Querying data -> {updated_qstr}")
# The trick here is to use the "image" keyword for ``seg`` to only download the
# image data only and this gives us direct filenames as well.
records, filenames = client.query(updated_qstr, key=keys, seg="image")
print(f"{len(records)} records retrieved. \n")

# We do a quick comparision to ensure the final results are the same.
# For this to work, we just need to deal with the different indexes.
print("Quick Comparison")
print(results.reset_index(drop=True) == records.reset_index(drop=True))

#####################################################
# From here you can now request (export) the data.
# This will download this specific subset of data to your
# local machine when the export request has been completed.
# Depending on the status of the JSOC, this might take a while.
#
# Please be aware the script will hold until the export is complete.

export = client.export(updated_qstr, method="url", protocol="fits")
files = export.download(Path("~/sunpy/").expanduser().as_posix())
# Our goal is to request data of a recent X-class flare.
# The X-class flare occurred on the 2021/07/03 at 14:30:00 UTC.
# We will focus on the 5 minutes before and after this time.
# What we want to do is only get the shorter exposures,
# which are going to be the flare related.

query = Fido.search(
attrs.Time("2021-07-03 14:25:00", "2021-07-03 14:35:00"),
attrs.jsoc.Series("aia.lev1_euv_12s"),
attrs.Wavelength(211 * u.AA),
attrs.jsoc.Notify(jsoc_email),
attrs.jsoc.Keyword("EXPTIME") <= 2,
attrs.jsoc.Segment("image"),
)

print(query)

#####################################################
# With AIA files, it is possible to bypass the export stage.
# We can manually construct the URLS of the the data.
# Be aware that each file will have the same filename based on the URL.
# You will have to then use your preferred downloader to download the files.
# Now we will download the data and "aia prep" the
# data with every feature of `aiapy` and plot the
# data sequence using `sunpy`.

urls = [f"http://jsoc.stanford.edu{filename}" for filename in filenames.image]
files = Fido.fetch(query)

#####################################################
# Now we will "prep" the data with every feature of
# `aiapy` and plot the data sequence using `sunpy`.

level_1_maps = sunpy.map.Map(files.download.to_list())
level_1_maps = sunpy.map.Map(files)
# We get the pointing table outside of the loop for the relevant time range.
# Otherwise you're making a call to the JSOC every single time.
pointing_table = get_pointing_table(level_1_maps[0].date - 3 * u.h, level_1_maps[-1].date + 3 * u.h)
Expand All @@ -139,8 +65,18 @@
map_registered = register(map_updated_pointing)
map_degradation = correct_degradation(map_registered, correction_table=correction_table)
map_normalized = map_degradation / map_degradation.exposure_time
level_15_maps.append(map_normalized)
bottom_left = SkyCoord(500 * u.arcsec, 100 * u.arcsec, frame=map_normalized.coordinate_frame)
top_right = SkyCoord(1500 * u.arcsec, 700 * u.arcsec, frame=map_normalized.coordinate_frame)
map_cropped = map_normalized.submap(bottom_left, top_right=top_right)
level_15_maps.append(map_cropped)

#####################################################
# Finally, we create a sequence of maps and animate it

sequence = sunpy.map.Map(level_15_maps, sequence=True)
sequence.peek()

fig = plt.figure()
ax = fig.add_subplot(projection=sequence.maps[0])
ani = sequence.plot(axes=ax, norm=ImageNormalize(vmin=0, vmax=1e3, stretch=AsinhStretch()))

plt.show()

0 comments on commit 5e50365

Please sign in to comment.