Skip to content

Commit

Permalink
Fix conflating the lunar disk and the solar disk
Browse files Browse the repository at this point in the history
  • Loading branch information
@ayshih
ayshih committed Apr 12, 2024
1 parent da6e9d6 commit b1911d4
Showing 1 changed file with 29 additions and 24 deletions.
53 changes: 29 additions & 24 deletions posts/2024/2024-04-03-eclipse.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@
"import astropy.units as u\n",
"import scipy.ndimage as ndimage\n",
"\n",
"from astropy.constants import R_earth\n",
"from astropy.coordinates import CartesianRepresentation, EarthLocation, SkyCoord\n",
"from matplotlib.patches import Circle\n",
"from skimage.color import rgb2gray\n",
Expand Down Expand Up @@ -163,10 +164,10 @@
"id": "4127a8f9-3230-4645-958a-0a919afea3cf",
"metadata": {},
"source": [
"## Find the Sun\n",
"## Find the Moon\n",
"\n",
"Now that we've loaded our image and accompanying metadata, the next step is to locate the edge of the Sun in our image.\n",
"This allows us to find the center of the Sun, which is needed when aligning our data, as well as allowing us to determine the scale of the Sun in the image.\n",
"Now that we've loaded our image and accompanying metadata, the next step is to locate the edge of the Moon in our image.\n",
"This allows us to find the center of the Moon, which is needed when aligning our data, as well as allowing us to determine the scale of the Moon in the image.\n",
"In order to do this we are going to use several different image manipulation techniques."
]
},
Expand All @@ -176,7 +177,7 @@
"metadata": {},
"source": [
"\n",
"We start with applying a [Gaussian blur](https://en.wikipedia.org/wiki/Gaussian_filter) to help segment the solar disk from the corona and mask all parts of the image above a given threshold."
"We start with applying a [Gaussian blur](https://en.wikipedia.org/wiki/Gaussian_filter) to help segment the lunar disk from the corona and mask all parts of the image above a given threshold."
]
},
{
Expand Down Expand Up @@ -239,7 +240,7 @@
"metadata": {},
"source": [
"Finally, we use scikit-image to apply the [Hough Transform](https://en.wikipedia.org/wiki/Hough_transform) to identify circles in the image.\n",
"We then use this to extract the size in pixels of the solar disk and its center."
"We then use this to extract the size in pixels of the lunar disk and its center."
]
},
{
Expand Down Expand Up @@ -292,7 +293,7 @@
"id": "c18ed50e-1ad3-4df2-a182-5d86f39ade26",
"metadata": {},
"source": [
"Lastly, let's add units to our circle that we fit the solar limb."
"Lastly, let's add units to our circle that we fit the lunar limb."
]
},
{
Expand Down Expand Up @@ -320,7 +321,7 @@
"id": "a61bd22f-57f1-4379-82e2-45ae1a343fbb",
"metadata": {},
"source": [
"At this point, we have our image data, it's metadata and we've located the Sun.\n",
"At this point, we have our image data, it's metadata and we've located the Moon.\n",
"Now we are ready to load our eclipse photograph into a SunPy `Map`!\n",
"A `Map` allows us to carry around both the data and metadata of our image together and, importantly, makes it easy to combine solar observations from multiple observatories."
]
Expand All @@ -330,8 +331,8 @@
"id": "c7fbab37",
"metadata": {},
"source": [
"First, we need to know the [angular size](https://en.wikipedia.org/wiki/Angular_diameter) of the Sun.\n",
"Thankfully, we can use sunpy to calculate the angular size of the solar disk as seen from Earth at the time we took our photograph."
"First, we define the observer location (i.e., where on Earth did we take our picture?) and the orientation of the Sun in the sky.\n",
"For the observer location, we can use the GPS coordinates from our image metadata."
]
},
{
Expand All @@ -341,17 +342,16 @@
"metadata": {},
"outputs": [],
"source": [
"rsun_obs = sunpy.coordinates.sun.angular_radius(camera_metadata[\"time\"])\n",
"print(rsun_obs)"
"loc = EarthLocation(lat=camera_metadata[\"gps\"][0], lon=camera_metadata[\"gps\"][1], height=camera_metadata[\"gps\"][2])\n",
"observer = loc.get_itrs(camera_metadata[\"time\"])"
]
},
{
"cell_type": "markdown",
"id": "d365d554-96bd-4fc9-bc84-60f0771a6c66",
"metadata": {},
"source": [
"Combining this angular radius with the radius of the solar disc in pixels gives us the angular size of the pixels in the image.\n",
"In the parlance of astronomical imaging, this is often referred to as the *plate scale*."
"Second, we determine the [angular size](https://en.wikipedia.org/wiki/Angular_diameter) of the Moon using its radius and its distance from the observer."
]
},
{
Expand All @@ -361,17 +361,20 @@
"metadata": {},
"outputs": [],
"source": [
"plate_scale = rsun_obs / im_radius[0]\n",
"print(plate_scale)"
"moon = SkyCoord(sunpy.coordinates.get_body_heliographic_stonyhurst('moon', camera_metadata[\"time\"], observer=observer, quiet=True))\n",
"R_moon = 0.2725076 * R_earth # IAU mean radius\n",
"dist_moon = SkyCoord(observer).separation_3d(moon)\n",
"moon_obs = np.arcsin(R_moon / dist_moon).to('arcsec')\n",
"print(moon_obs)"
]
},
{
"cell_type": "markdown",
"id": "bd49ea6d-14a1-4591-bcda-c2b1949b0e8f",
"metadata": {},
"source": [
"The last two pieces of information we need, are the location of our observer (i.e., where on Earth did we take our picture?) and the orientation of the Sun in the sky.\n",
"For the observer location, we can use the GPS coordinates from our image metadata."
"Combining this angular radius with the radius of the lunar disk in pixels gives us the angular size of the pixels in the image.\n",
"In the parlance of astronomical imaging, this is often referred to as the *plate scale*."
]
},
{
Expand All @@ -381,15 +384,16 @@
"metadata": {},
"outputs": [],
"source": [
"loc = EarthLocation(lat=camera_metadata[\"gps\"][0], lon=camera_metadata[\"gps\"][1], height=camera_metadata[\"gps\"][2])"
"plate_scale = moon_obs / im_radius[0]\n",
"print(plate_scale)"
]
},
{
"cell_type": "markdown",
"id": "4044ce0d",
"metadata": {},
"source": [
"We can then use that location to calculate the orientation of the Sun as seen from that location on Earth.\n",
"We also use the observer location to calculate the orientation of the Sun as seen from that location on Earth.\n",
"This gives us a rotation angle between our image coordinate system and the solar north pole.\n",
"If your camera is not perfectly horizontal then you may need to fudge this value slightly."
]
Expand Down Expand Up @@ -420,11 +424,12 @@
"metadata": {},
"outputs": [],
"source": [
"frame = sunpy.coordinates.Helioprojective(observer=loc.get_itrs(camera_metadata[\"time\"]),\n",
" obstime=camera_metadata[\"time\"])\n",
"frame = sunpy.coordinates.Helioprojective(observer=observer, obstime=camera_metadata[\"time\"])\n",
"moon_hpc = moon.transform_to(frame)\n",
"\n",
"header = make_fitswcs_header(\n",
" im,\n",
" SkyCoord(0, 0, unit=u.arcsec, frame=frame),\n",
" moon_hpc,\n",
" reference_pixel=u.Quantity([im_cx, im_cy]),\n",
" scale=u.Quantity([plate_scale, plate_scale]),\n",
" rotation_angle=solar_rotation_angle,\n",
Expand Down Expand Up @@ -460,7 +465,7 @@
"id": "3bb5a923-0c08-4a6a-8766-f6caef7b147a",
"metadata": {},
"source": [
"## Find the Stars\n",
"## Find a Star\n",
"\n",
"Though we already have all of the metadata we need to make a SunPy `Map`, we can further improve the accuracy by \n",
"locating a known star in the image and using that to determine the rotation angle.\n",
Expand Down Expand Up @@ -639,7 +644,7 @@
"source": [
"header = make_fitswcs_header(\n",
" im,\n",
" SkyCoord(0, 0, unit=u.arcsec, frame=frame),\n",
" moon_hpc,\n",
" reference_pixel=u.Quantity([im_cx, im_cy]),\n",
" scale=u.Quantity([plate_scale, plate_scale]),\n",
" rotation_angle=solar_rotation_angle + fudge_angle,\n",
Expand Down

0 comments on commit b1911d4

Please sign in to comment.