Finding contours of a map

This example shows how to find and plot contours on a map.

Reference: corresponding sunpy webpage

using SunPy, PythonCall
using SunPy: Map
@py import sunpy.data.sample: AIA_193_IMAGE
@py import astropy.units as u

aiamap = Map(AIA_193_IMAGE)

contours =  SunPy.find_contours_py(aiamap, 50000.0)
contours_py =  aiamap.find_contours(50000 * u.DN)

Python:
[<SkyCoord (Helioprojective: obstime=2011-06-07T06:33:08.840, rsun=696000.0 km,…
    (-0.00405795, 0.04788327, 1.51846027e+11)>): (Tx, Ty) in arcsec
    [(683.27903798, -439.85353936), (684.50944624, -440.83467707),
     (685.24686413, -439.85158592), (685.59350099, -437.44823271),
     (684.50539784, -436.72726943), (682.10283552, -437.17844294),
     (681.93892109, -437.4518604 ), (682.10383581, -438.1933212 ),
     (683.27903798, -439.85353936)]>, <SkyCoord (Helioprojective: obstime=2011-…
    (-0.00405795, 0.04788327, 1.51846027e+11)>): (Tx, Ty) in arcsec
    [(692.9544155 , -437.44092553), (694.11892783, -438.25296287),
     (694.88704773, -437.43900692), (696.51900478, -435.28112493),
                        ... 998 more lines ...
     (-863.51063935, 387.66138148), (-863.85341055, 388.14142254),
     (-866.04479465, 390.06186699), (-866.2585235 , 390.29463432),
     (-868.36589515, 392.4625638 ), (-868.66404582, 392.8634119 ),
     (-871.06865301, 394.50397712), (-873.46982828, 392.66287077),
     (-873.57156947, 392.45739089), (-875.87051858, 390.32996917),
     (-876.2812805 , 390.05169542), (-877.82731085, 387.64715637),
     (-877.53399293, 385.24444497), (-876.74779008, 382.84222314),
     (-876.04252634, 380.43992072), (-875.86042381, 380.08816137),
     (-873.52281716, 378.03942066)]>]

@py import matplotlib.pyplot as plt

let fig = plt.figure(figsize=(12, 5))
    # First subplot: Julia implementation
    ax = fig.add_subplot(1, 2, 1, projection=aiamap)
    # aiamap.plot(axes=ax)
    ax.plot_coord.(contours)
    ax.set_title("Julia find_contours")

    # Second subplot: Python implementation
    ax2 = fig.add_subplot(1, 2, 2, projection=aiamap)
    # aiamap.plot(axes=ax2)
    ax2.plot_coord.(PyList(contours_py))
    ax2.set_title("Python find_contours")

    fig
end

Acceleration

julia> using Chairmarks
julia> # Python
       @b aiamap.find_contours(50000 * u.DN)175.753 ms (16 allocs: 352 bytes)
julia> # vs Julia
       @b SunPy.find_contours_py(aiamap, 50000.0)7.007 ms (9477 allocs: 363.367 KiB)