diff --git a/examples/tracing_loops.py b/examples/tracing_loops.py index 3cdf9286..493dfad2 100644 --- a/examples/tracing_loops.py +++ b/examples/tracing_loops.py @@ -1,18 +1,18 @@ """ -===================== -Tracing Coronal Loops -===================== +================================================ +Tracing Coronal Loops and Extracting Intensities +================================================ This example traces out the coronal loops in a FITS image -using `~sunkit_image.trace.occult2`. +using `~sunkit_image.trace.occult2` and then extracts the intensity +along one traced loop. In this example we will use the settings and the data from Markus Aschwanden's tutorial on his IDL implementation of the ``OCCULT2`` algorithm, which can be found `here `__. -The aim of this example is to demonstrate that `~sunkit_image.trace.occult2` provides similar -results compared to the IDL version. """ +# sphinx_gallery_thumbnail_number = 1 # NOQA: ERA001 import matplotlib.pyplot as plt import numpy as np @@ -70,4 +70,30 @@ fig.tight_layout() +############################################################################### +# Finally, we can use the traced loops location information to extract the intensity values. + +# Since we only currently get pixel locations, we need to get the word coordinates of the first loop. +first_loop = np.array(loops[0]) +loop_coords = trace_map.pixel_to_world(first_loop[:, 0] * u.pixel, first_loop[:, 1] * u.pixel) + +# Now we can extract the intensity along the loop +intensity = sunpy.map.sample_at_coords(trace_map, loop_coords) + +# Finally, we can calculate the angular separation along the loop +angular_separation = loop_coords.separation(loop_coords[0]).to(u.arcsec) + +# Plot the loop location and its intensity profile +fig = plt.figure(figsize=(10, 4)) +ax = fig.add_subplot(121, projection=trace_map) +trace_map.plot(axes=ax, clip_interval=(1, 99.99) * u.percent) +ax.plot_coord(loop_coords, color="r") + +ax = fig.add_subplot(122) +ax.plot(angular_separation, intensity) +ax.set_xlabel("Distance along loop [Arcsec]") +ax.set_ylabel("Intensity") + +fig.tight_layout() + plt.show()