Subversion Repositories colinrmitchell.com

Rev

Rev 673 | Blame | Compare with Previous | Last modification | View Log | RSS feed


import glob
from skimage import io
from skimage.color import rgb2gray, rgb2hsv, hsv2rgb
import numpy as np
import cv2
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
import datetime
import math
import pickle
from os.path import exists

# No X needed for servers.
mpl.use('Agg')

# Horizontal band of the image to analyze
# x in pixels
#  Can be used to remove nautical dusk/dawn
x_margin = 65
# y in percents of image height
#  Can be used to remove moon path or other obstructions
#  The top of the image is 0%
min_y_p, max_y_p = 50, 85

# The HSV hue value (0-1) that colors must differ by
# It will reduce the number of K-means values until
# the hues all differ by more than this value.
min_color_distance = 0.10

# The HSV hue value that is the threshold for clear skies.
hue_clear_limit = 0.35

# Need a link to the keograms image directory
files = glob.glob("keograms/*.jpg")

class Night:
   freqs = []
   palette = []
   date = datetime.datetime.now()
   filename = ""

nights = []

# Stores a cache file containing the processing
# of previous pictures.  Delete this file to
# reprocess all keograms.
if exists("cache.pkl"):
   with open("cache.pkl", mode='rb') as hand:
      nights = pickle.load(hand)

for file in files:
   # Don't process if we have cached data
   cont = False
   
   for night in nights:
      if night.filename == file:
         cont = True
         
   if cont:
      continue
     
   fimg = io.imread(file)
   
   # Crop image
   min_y = int(len(fimg) * (min_y_p / 100))
   max_y = int(len(fimg) * (max_y_p / 100))
   
   img = fimg[min_y:max_y, x_margin:len(fimg[0]) - x_margin]
   
   # Put all pixels in a 1x array  
   pixels = np.float32(img.reshape(-1, 3))

   # Remove the grid lines
   pixels = pixels[[t[0] for t in ~np.isin(pixels, [181, 184, 189])]]

   hues = rgb2hsv(pixels)[:, 0]
   
   # Number of distinct colors to solve for
   n_colors = 4

   while n_colors > 0:
      criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 500, 0.001)
      flags = cv2.KMEANS_RANDOM_CENTERS

      _, labels, palette = cv2.kmeans(hues, n_colors, None, criteria, 10, flags)
      _, counts = np.unique(labels, return_counts=True)

      # How often each color occurs
      freqs = [c / len(hues) for c in counts]

      pairs = [(a, b) for idx, a in enumerate(palette) for b in palette[idx + 1:]]

      b = True

      # We check the distance between all pairs of the solved colors.
      # If any distance is less than the threshold, reduce the number
      # of colors to solve for and try again.
      for pair in pairs:
         dist = np.linalg.norm(pair[0] - pair[1])
         
         if dist < min_color_distance:
            b = False

      n_colors -= 1
     
      if b:
         break

   night = Night()

   night.freqs = freqs
   # Convert the hues into "normalized" RGB values.
   night.palette = [hsv2rgb([h[0], 0.3, 128]) for h in palette]
   night.date = datetime.datetime.strptime(file[17:25], "%Y%m%d")
   night.filename = file

   nights.append(night)

# Save the cache file.
with open("cache.pkl", mode='wb') as hand:
   hand.truncate(0)
   pickle.dump(nights, hand)

# Clear night chart
fig, ax1 = plt.subplots(1, figsize=(12,6))

clears = [0, 0]

for n, night in enumerate(nights):
   rows = list(zip(night.freqs, night.palette))
   
   blue_count = 0
   
   # Considered blue if hue > hue_clear_limit
   for row in rows:  
      if rgb2hsv(row[1])[0] > hue_clear_limit:
         blue_count += row[0]
     
   # A clear night is at least 75% blue sky
   if blue_count > 0.75:
      clears[0] += 1
   else:
      clears[1] += 1

plt.bar("Clear", clears[0], 1, bottom=0, color=np.array([59, 67, 80]) / 255)
plt.bar("Cloudy", clears[1], 1, bottom=0, color=np.array([108, 108, 82]) / 255)

ax1.set_title("Nights over 75% clear")
ax1.set_ylabel('Count')
plt.grid(True, axis="y", linestyle='--')

plt.savefig("bar.png")

# Night clearness percentage
fig, ax1 = plt.subplots(1, figsize=(12,6))

# Quantized to 10%
percents = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

for n, night in enumerate(nights):
   rows = list(zip(night.freqs, night.palette))
   
   blue_count = 0
   
   # Considered blue if hue > hue_clear_limit
   for row in rows:  
      if rgb2hsv(row[1])[0] > hue_clear_limit:
         blue_count += row[0]
   
   percents[round(round(blue_count * 100) / 10)] += 1

for i, percent in enumerate(percents):
   plt.bar(int(i * 10), percent, 10, bottom=0, color=hsv2rgb([0.2 + i * 0.4 / 10, 0.3, 128]) / 255)

ax1.set_title("Night Clearness Percentage")
ax1.set_xlabel('Percentage')
ax1.set_ylabel('Count')
plt.grid(True, axis="y", linestyle='--')

plt.savefig("percents.png")

# Color histogram
fig, ax1 = plt.subplots(1, figsize=(12,6))

hist_divisions = 10

hues_count = []

for n, night in enumerate(nights):
   for t in zip(night.freqs, night.palette):
      for x in range(1, int(t[0] * 100)):
         hues_count.append(round(hist_divisions * rgb2hsv(t[1])[0]) / hist_divisions)

counts = 0
hues = np.unique(hues_count)

colors = [hsv2rgb([h, 0.3, 128]) / 255 for h in hues]

for hue in hues:
   counts += hues_count.count(hue)

for hue in hues:
   plt.bar(hue, 100 * hues_count.count(hue) / counts, 1 / hist_divisions, bottom=0, color=hsv2rgb([hue, 0.3, 128]) / 255)

ax1.set_title("Histogram of Sky Color by Time")
ax1.set_xlabel('Color (HSV Hue)')
ax1.set_ylabel('Time (Percent)')
plt.grid(True, axis="y", linestyle='--')

plt.savefig("histo.png")

# Trend graph
fig, ax1 = plt.subplots(1, figsize=(12,6))

for n, night in enumerate(nights):
   # Sort by blue channel of color
   sorted_rows = list(zip(night.freqs, night.palette))  
   sorted_rows.sort(key = lambda x: (x[1] * (128 / rgb2gray(x[1])))[2])
   
   # Turn absolute ratios into stacked bar values.
   sorted_rows = [(sum([t1[0] for t1 in sorted_rows[0:i+1]]), t[1]) for i, t in enumerate(sorted_rows)]
   
   pal = [t[1] for t in sorted_rows]
   rows = [100 * t[0] for t in sorted_rows]
   pal.reverse()
   rows.reverse()
   
   color_map = np.vectorize(lambda c: min(c, 255) / 255)
   
   plt.bar(night.date, rows, 1, bottom=0, color=[color_map(c) for c in pal])

plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%m/%d/%Y'))
plt.gca().xaxis.set_major_locator(mdates.DayLocator())
plt.gcf().autofmt_xdate()
# Only show every 7th day label
ax1.xaxis.set_major_locator(ticker.MultipleLocator(7))

ax1.set_title('Nightly Sky Condition')
ax1.set_xlabel('Date')
ax1.set_ylabel('Sky Color in Percent Time')
plt.grid(True, axis="y", linestyle='--')

plt.savefig("trends.png")