"""A module providing facilities to analyze the moon in an image.
Methods in this module predominantly deal with determining the relationship
between the phase of the moon and the apparent size of the moon in an all-sky
image. The phase of the moon is provided on a scale from 0.0 (a new moon) to
1.0 (a full moon). Methods are provided to find the position of the moon and
sun. The phase of the moon in an eclipse and outside of an eclipse are
separated into their own methods. There is also a method provided to plot this
data and the model linking the apparent size of the image to the moon phase.
"""
import os
import math
import warnings
import ephem
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from astropy.modeling import models, fitting
import coordinates
import image
# Sets up a pyephem object for the camera.
camera = ephem.Observer()
camera.lat = "31.959417"
camera.lon = "-111.598583"
camera.elevation = 2120
[docs]def eclipse_phase(d):
"""Calculate the proportion of the moon that is lit up during an eclipse.
Parameters
----------
d : float
The distance between the center of Earth"s shadow and the center of the
moon in kilometers.
Returns
-------
float
The phase of the moon, ranging between 0.0 and 1.0, where 0.0 is a new
moon and 1.0 is a full moon.
"""
# In kilometers.
r = 1737 # R_moon
R = 4500 #R_earth
# This makes addition work as addition and not concatenates
d = np.asarray(d)
d = np.abs(d) # Required as for after totality times d < 0
r2 = r * r
R2 = R * R
d2 = np.square(d)
# Part 1 of the shaded area equation
a = (d2 + r2 - R2)
b = d * 2 * r
p1 = r2 * np.arccos(a / b)
# Part 2 of the shaded area equation
a = (d2 + R2 - r2)
b = d * 2 * R
p2 = R2 * np.arccos(a / b)
# Part 3 of the shaded area equation
a1 = (r + R - d)
a2 = (d + r - R)
a3 = (d - r + R)
a4 = (d + r + R)
p3 = (0.5) * np.sqrt(a1 * a2 * a3 * a4)
# Add them together to get the shaded area
A = p1 + p2 - p3
# Get the shaded proportion by divding the shaded area by the total area
# Assumes r is the radius of the moon being shaded.
P = A / (np.pi * r2)
# P is the shaded area, so 1-P is the lit up area.
P = 1 - P
return P
# 1.0 = Full moon, 0.0 = New Moon
[docs]def moon_phase(img):
"""Calculate the proportion of the moon that is lit up for non-eclipse
nights.
Parameters
----------
img : image.AllSkyImage
The image.
Returns
-------
float
The phase of the moon, ranging between 0.0 and 1.0, where 0.0 is a new
moon and 1.0 is a full moon.
"""
# Sets the calculation date.
camera.date = img.formatdate
# Makes a moon object and calculates it for the observation location/time
moon = ephem.Moon()
moon.compute(camera)
return moon.moon_phase
[docs]def moon_size(img):
"""Calculate the area of the moon in pixels in a given image.
Parameters
----------
img : image.AllSkyImage
The image.
Returns
-------
int
The size of the moon in pixels.
Notes
-----
This method first converts the image to a black and white two-tone image
where pixels with greyscale values above or equal to 250 is set to
white and everything below is set to black. A binary closing is performed
on the image, which smooths over any small black pixel regions within the
moon. These black regions are created when the pixels are brighter than the
maximum of 255 for white pixels and overflow back to 0.
The white regions are labeled and their sizes are found using
ndimage.label. Then, the approximate position of the center of the moon
is found using find_moon.
If the pixel at the moon"s center is black (due to aforementioned pixel
value overflow), the nearest white region along the x axis is found and
the size of this region is returned.
If this pixel is white, the size of this white region is returned.
"""
thresh = 5
dat = np.where(img.data >= 255 - thresh, 1, 0)
# Runs a closing to smooth over local minimums (which are mainly caused by
# a rogue antenna). Then labels the connected white regions. Structure s
# Makes it so that regions connected diagonally are counted as 1 region.
s = [[1, 1, 1], [1, 1, 1], [1, 1, 1]]
dat = ndimage.morphology.binary_closing(dat, structure=s)
labeled, nums = ndimage.label(dat, structure=s)
# Want to find the size of each labeled region.
sizes = [0] * (nums + 1)
for row in labeled:
for val in row:
sizes[val] = sizes[val] + 1
# We want to exclude the background from the "biggest region" calculation.
# It"s quicker to just set the background (0) region to 0 than to check
# every single value above before adding it to the array.
sizes[0] = 0
# Following code calculates d, the distance between the center of
# Earth"s shadow and the center of the moon. Basically just d = v*t.
# Use pyephem to find the labeled region that the moon is in.
posx, posy, _ = find_moon(img)
posx = math.floor(posx)
posy = math.floor(posy)
reg = labeled[posy, posx]
# Very large and bright moons have a dark center (region 0) but I want the
# region of the moon.
while reg == 0 and posx < 511:
posx = posx + 1
reg = labeled[posy, posx]
biggest = sizes[reg]
return biggest
[docs]def find_moon(img):
"""Find the (x, y, alt) coordinate of the moon"s center in a given image.
Parameters
----------
img : image.AllSkyImage
The image.
Returns
-------
x : float
The x coordinate of the moon"s center.
y : float
The y coordinate of the moon"s center.
alt : float
The altitude angle of the moon"s center.
Notes
-----
The x and y coordinates are corrected for irregularities in the lens using
coordinates.galactic_conv.
"""
# Sets the date of calculation.
camera.date = img.formatdate
# Calculates the moon position.
moon = ephem.Moon()
moon.compute(camera)
# Conversion to x,y positions on the image.
alt = np.degrees(moon.alt)
az = np.degrees(moon.az)
x, y = coordinates.altaz_to_xy(alt, az)
x, y = coordinates.galactic_conv(x, y, az)
return (x, y, alt)
[docs]def find_sun(img):
"""Find the (alt, az) coordinate of the sun"s center in a given image.
Parameters
----------
img : image.AllSkyImage
The image.
Returns
-------
alt : float
The altitude coordinate of the sun"s center.
az : float
The azimuth coordinate of the sun"s center.
"""
# Sets the date of calculation.
camera.date = img.formatdate
# Calculates the sun position.
sun = ephem.Sun()
sun.compute(camera)
# Conversion to x,y positions on the image.
alt = np.degrees(sun.alt)
az = np.degrees(sun.az)
return (alt, az)
# Fits a Moffat fit to the moon and returns the estimated radius of the moon.
# Radius of the moon is the FWHM of the fitting function.
[docs]def fit_moon(img, x, y):
"""Fit a Moffat function to the moon in a given image.
Parameters
---------
img : image.AllSkyImage
The image.
x : float
The x coordinate of the moon"s center.
y : float
The y coordinate of the moon"s center.
Returns
-------
float
The Full Width at Half Maximum of the Moffat function.
"""
# This block of code runs straight vertical from the center of the moon
# It gives a predicted rough radius of the moon, it starts counting at the
# first white pixel it encounters (the center may be black)
# and stops at the last white pixel. White here defined as > 250 greyscale.
yfloor = math.floor(y)
count = False
size = 0
xfloor = math.floor(x)
start = xfloor
# The only reason we have this if block is to ensure we don"t run for
# moon radii greater than 35 in this case.
for i in range(0, 35):
start += 1
# Breaks if it reaches the edge of the image.
if start == img.data.shape[1]:
break
if not count and img.data[yfloor, start] >= 250:
count = True
elif count and img.data[yfloor, start] >= 250:
size += 1
elif count and img.data[yfloor, start] < 250:
break
# Add some buffer pixels in case the center is black and the edges of the
# moon are fuzzed and then convert radius to diameter.
size = (size + 10) * 2
# Makes sure the lower/upper slices don"t out of bounds error.
lowerx = xfloor - size if (xfloor - size > 0) else 0
lowery = yfloor - size if (yfloor - size > 0) else 0
upperx = xfloor + size if (xfloor + size < 511) else 511
uppery = yfloor + size if (yfloor + size < 511) else 511
# Size of the moon enclosing square.
deltax = (upperx - lowerx)
deltay = (uppery - lowery)
# Creates two arrays, with the array values being the x or y coordinate of
# that location in the array.
y, x = np.mgrid[0:deltay, 0:deltax]
# Slices out the moon square and finds center coords.
z = img.data[lowery:uppery, lowerx:upperx]
midy = deltay / 2
midx = deltax / 2
# Moffat fit, centered in square, stdev of 20 as a start.
stddev = 20
model_init = models.Moffat2D(amplitude=200, x_0=midx, y_0=midy,
gamma=stddev)
fit = fitting.LevMarLSQFitter()
with warnings.catch_warnings():
# Ignore model linearity warning from the fitter
warnings.simplefilter("ignore")
model = fit(model_init, x, y, z)
# /2 is average FWHM but FWHM = diameter, so divide by two again.
#fwhm = (model.x_fwhm + model.y_fwhm) / 4
fwhm = model.fwhm / 2
return fwhm
# Generates the size vs illuminated fraction for the two eclipse nights.
[docs]def generate_eclipse_data(regen=False):
"""Generate moon phase data for eclipses.
Parameters
----------
regen : bool, optional
If True, regen from scratch. Otherwise, load it from a file.
Defaults to False.
Returns
-------
truevis : list of lists
A list containing one or more lists where each list contains values
corresponding to an eclipse. Each value is the phase of the moon
ranging from 0.0 to 1.0, in an image taken during that night.
0.0 corresponds to a new moon, while 1.0 corresponds to a full moon.
The list is ordered in chronological order; that is, the first value
within the list is calculated from an image that was taken at the
beginning of the eclipse, and the last value is taken from the end of
the eclipse. Currently, this is hardcoded such that the first list
represents the eclipse on 2018/01/31, and the second list represents
the eclipse on 2015/04/04.
imvis : list of lists
A list containing one or more lists where each list contains values
corresponding to an eclipse. Each value is the area of the moon,
in pixels, in an image taken during that night. The list is ordered in
chronological order; that is, the first value within the list is
calculated from an image that was taken at the beginning of the eclipse,
and the last value is taken from the end of the eclipse.
Currently, this is hardcoded such that the first list represents
the eclipse on 2018/01/31, and the second list represents the eclipse
on 2015/04/04.
"""
dates = ["20180131", "20150404"]
# Function within a function to avoid code duplication.
def data(date):
# Necessary lists
distances = []
imvis = []
truevis = []
# Check to see if the data has been generated already.
# If it has then read it from the file.
save = os.path.join(os.path.dirname(__file__), *["data", "eclipse-" + date + ".txt"])
if os.path.isfile(save) and not regen:
f = open(save)
for line in f:
line = line.rstrip().split(",")
truevis.append(float(line[0]))
imvis.append(float(line[1]))
f.close()
return (truevis, imvis)
# If we're regenerating the data we do it here.
directory = os.path.join(os.path.dirname(__file__), *["Images", "Original", "KPNO", date])
images = sorted(os.listdir(directory))
# I need a better way to check this.
if ".DS_Store" in images:
images.remove(".DS_Store")
# Finds the size of the moon in each image.
for name in images:
print(name)
print(date)
print(directory)
img = image.load_image(name, date, "KPNO")
# This basically hacks us to use the center of the earth as our
# observation point.
camera.elevation = - ephem.earth_radius
camera.date = img.formatdate
# Calculates the sun and moon positions.
moon = ephem.Moon()
sun = ephem.Sun()
moon.compute(camera)
sun.compute(camera)
# Finds the angular separation between the sun and the moon.
sep = ephem.separation((sun.az, sun.alt), (moon.az, moon.alt))
# Radius of moon orbit to convert angular separation -> distance
R = 385000
# For angles this small theta ~ sin(theta), so I dropped the sine
# to save computation time.
# Angle between moon and earth"s shadow + angle between moon and sun
# should ad d to pi, i.e. the earth"s shadow is across from the sun.
d = R * (np.pi - sep)
size = moon_size(img)
imvis.append(size)
distances.append(d)
print("Processed: " + date + "/" + name)
# Calculates the proportion of visible moon for the given distance
# between the centers.
truevis = eclipse_phase(distances)
imvis = np.asarray(imvis)
# If the moon is greater than 40,000 pixels then I know that the moon
# has merged with the light that comes from the sun and washes out the
# horizon.
imvis = np.where(imvis < 80000, imvis, float("NaN"))
f = open(save, "w")
# Writes the data to a file so we can read it later for speed.
for i in range(0, len(truevis)):
f.write(str(truevis[i]) + "," + str(imvis[i]) + "\n")
f.close()
return (truevis, imvis)
trues = []
ims = []
for date in dates:
true, im = data(date)
trues.append(true)
ims.append(im)
return (trues, ims)
[docs]def moon_circle(frac):
"""Calculate the estimated pixel radius of the moon based on the
fraction of the moon that is illuminated.
Parameters
----------
frac : float
The proportion of the moon that is illuminated by sunlight.
Returns
-------
float
The estimated radius of the moon.
Notes
-----
The model used in this method to convert between the fraction of the moon
that is illuminated to the estimated pixel area was found by plotting the
moon pixel area versus the moon phase and picking representative points
to model the relation. The
model is designed to always overestimate the size of the moon. The model
is defined by interpolating from the following table of representative
points:
========================= ==================
Moon fraction illuminated Moon area (pixels)
------------------------- ------------------
0 650
0.345 4000
0.71 10500
0.88 18000
0.97 30000
1.0 35000
========================= ==================
"""
illuminated = [0, 0.345, 0.71, 0.88, 0.97, 1.0]
size = [650, 4000, 10500, 18000, 30000, 35000]
A = np.interp(frac, illuminated, size)
return np.sqrt(A/np.pi)
[docs]def moon_mask(img):
"""Generate a masking array that covers the moon in a given image.
Parameters
---------
img : image.AllSkyImage
The image.
Returns
-------
numpy.ndarray
An array where pixels inside the moon are marked with False and those
outside the moon are marked with True.
"""
# Get the fraction visible for interpolation and find the
# location of the moon.
vis = moon_phase(img)
x, y, _ = find_moon(img)
# Creates the circle patch we use.
r = moon_circle(vis)
circ = Circle((x, y), r, fill=False)
# The following code converts the patch to a 512x512 mask array, with True
# for values outside the circle and False for those inside.
# This is the same syntax as np.ma.make_mask returns.
# This section of code generates an 262144x2 array of the
# 512x512 pixel locations. 262144 = 512^2
points = np.zeros((512**2, 2))
index = 0
for i in range(0, 512):
for j in range(0, 512):
# These are backwards as expected due to how reshape works later.
# Points is in x,y format, but reshape reshapes such that
# it needs to be in y,x format.
points[index, 0] = j
points[index, 1] = i
index += 1
# Checks all the points are inside the circle, then reshapes it to the
# 512x512 size.
mask = circ.contains_points(points)
mask = mask.reshape(512, 512)
return mask
[docs]def generate_plots():
"""Generate a plot of illuminated fraction versus apparent moon size.
Notes
-----
The eclipse dataset is loaded using generate_eclipse_data and then plotted.
The file images.txt contains the illuminated fraction of the moon
and the pixel area of the moon in that image. This data is plotted on top
of the eclipse data. The illuminated fraction of the moon is found using
moon_phase, and the size of the moon in the image is found using moon_size.
On top of this data the theoretical moon size model used in moon_circle is
plotted. Once all of these are plotted, two versions of the plot are saved,
one with a standard y and x axis, and one with a logarithmic y axis.
These plots are saved directly to Images/ under the names "moon-size.png"
and "moon-size-log.png."
"""
# Loads the eclipse data
vis, found = generate_eclipse_data()
print("Eclipse data loaded!")
# Eclipse normalization code.
#found[0] = np.asarray(found[0]) / np.nanmax(found[0])
#found[1] = np.asarray(found[1]) / np.nanmax(found[1])
# Plots the two eclipses, the first in blue (default), the second in green
plt.scatter(vis[0], found[0], label="2018/01/31 Eclipse", s=7)
#plt.scatter(vis[1], found[1], label="2015/04/04 Eclipse", s=7, c="g")
plt.ylabel("Approx Moon Size (pixels)")
plt.xlabel("Illuminated Fraction")
# Vis is the portion of the moon illuminated by the sun that night
# Found is the approximate size of the moon in the image
vis = []
found = []
loc = os.path.join(os.path.dirname(__file__), *["data", "images.txt"])
with open(loc, "r") as f:
for line in f:
line = line.rstrip()
info = line.split(",")
img = image.load_image(info[1], info[0], "KPNO")
vis.append(moon_phase(img))
found.append((moon_size(img)))
print("Processed: " + info[0] + "/" + info[1] + ".png")
# Removes out any moons that appear too large in the images to be
# considered valid.
found = np.asarray(found)
found = np.where(found < 40000, found, float("NaN"))
# Normalizes the non eclipse data.
#found1 = found / np.nanmax(found)
# Adds the noneclipse data to the plot.
plt.scatter(vis, found, label="Regular", s=7)
# This plots the estimated model of moon size on top of the graph.
vis2 = [0, 0.345, 0.71, 0.88, 0.97, 1.0]
found2 = [650, 4000, 10500, 18000, 30000, 35000]
plt.plot(vis2, found2, label="Model", c="r")
# Interpolation estimate for the moon size in the image based on the
# illuminated fractions.
found3 = np.interp(vis, vis2, found2)
plt.scatter(vis, found3, label="Interpolated", s=7)
plt.legend()
# Saves the figure, and then saves the same figure with a log scale.
plt.savefig("Images/moon-size.png", dpi=256)
ax = plt.gca()
ax.set_yscale("log")
plt.savefig("Images/moon-size-log.png", dpi=256)
plt.close()
if __name__ == "__main__":
# 20160326/r_ut071020s70020
generate_plots()