Take a photo of earth from space

README.md

Get a current photo of earth from space!

All you have to do is run

make

download.sh

PRODUCT=s3://noaa-goes16/ABI-L2-MCMIPF
NOW=$(date -v-20M -u '+%Y/%j/%H')
DIRECTORY="${PRODUCT}/${NOW}/"
LATEST=$(aws s3 ls $DIRECTORY --no-sign-request \
  | tail -1 \
  | awk '{print $NF}')
aws s3 cp ${DIRECTORY}${LATEST} data/full_disk.nc --no-sign-request

fulldisk.py

from datetime import datetime

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import metpy  # noqa: F401
import numpy as np
import xarray

FILE = ('data/full_disk.nc')
C = xarray.open_dataset(FILE)

# Scan's start time, converted to datetime object
scan_start = datetime.strptime(C.time_coverage_start, '%Y-%m-%dT%H:%M:%S.%fZ')
# Scan's end time, converted to datetime object
scan_end = datetime.strptime(C.time_coverage_end, '%Y-%m-%dT%H:%M:%S.%fZ')
# File creation time, convert to datetime object
file_created = datetime.strptime(C.date_created, '%Y-%m-%dT%H:%M:%S.%fZ')
# The 't' variable is the scan's midpoint time
midpoint = str(C['t'].data)[:-8]
scan_mid = datetime.strptime(midpoint, '%Y-%m-%dT%H:%M:%S.%f')

print('Scan Start    : {}'.format(scan_start))
print('Scan midpoint : {}'.format(scan_mid))
print('Scan End      : {}'.format(scan_end))
print('File Created  : {}'.format(file_created))
print('Scan Duration : {:.2f} minutes'.format((scan_end-scan_start).seconds/60))


# True Color RGB Recipe
# ---------------------
# 
# Color images are a Red-Green-Blue (RGB) composite of three different
# channels. To make a "Natural True Color" image we assign the following
# channels as our R, G, and B values:
# 
# | --                   | RED         | GREEN          | BLUE         |
# |----------------------|-------------|----------------|--------------|
# | **Name**             | Red Visible | Near-IR Veggie | Blue Visible |
# | **Wavelength**       | 0.64 µm     | 0.86 µm        | 0.47 µm      |
# | **Channel**          | 2           | 3              | 1            |
# | **Units**            | Reflectance | Reflectance    | Reflectance  |
# | **Range of Values**  | 0-1         | 0-1            | 0-1          |
# | **Gamma Correction** | 2.2         | 2.2            | 2.2          |
# 
# 
# Some important details to know about...
# 
# **Value Range**: The data units of channel 1, 2, and 3 are in reflectance and
# have a range of values between 0 and 1. RGB values must also be between 0 and
# 1.
# 
# **Gamma Correction**: A gamma correction is applied to control the brightness
# and make the image not look too dark.
# `corrected_value = value^(1/gamma)`.
# Most displays have a decoding gamma of 2.2. Read more about gamma correction
# at the following links...
# [source1](https://en.wikipedia.org/wiki/Gamma_correction) and
# [source2](https://www.cambridgeincolour.com/tutorials/gamma-correction.htm)).
# 
# **True Green**: The GREEN "veggie" channel on GOES-16 does not measure
# visible green light. Instead, it measures a near-infrared band sensitive to
# chlorophyll. We could use that channel in place of green, but it would make
# the green in our image appear too vibrant. Instead, we will tone-down the
# green channel by interpolating the value to simulate a natural green color.
# 
#       `TrueGreen = (0.45*RED) + (0.1*GREEN) + (0.45*BLUE)`
# 
# Now we can begin putting the pieces together...
# 
# 

# In[5]:


# Confirm that each band is the wavelength we are interested in
for band in [2, 3, 1]:
    print('{} is {:.2f} {}'.format(
        C['band_wavelength_C{:02d}'.format(band)].long_name,
        float(C['band_wavelength_C{:02d}'.format(band)][0]),
        C['band_wavelength_C{:02d}'.format(band)].units))


# In[6]:


# Load the three channels into appropriate R, G, and B variables
R = C['CMI_C02'].data
G = C['CMI_C03'].data
B = C['CMI_C01'].data


# In[7]:


# Apply range limits for each channel. RGB values must be between 0 and 1
R = np.clip(R, 0, 1)
G = np.clip(G, 0, 1)
B = np.clip(B, 0, 1)


# In[8]:


# Apply a gamma correction to the image to correct ABI detector brightness
gamma = 2.2
# gamma = 1.8
R = np.power(R, 1/gamma)
G = np.power(G, 1/gamma)
B = np.power(B, 1/gamma)


# In[9]:


# Calculate the "True" Green
G_true = 0.45 * R + 0.1 * G + 0.45 * B
G_true = np.clip(G_true, 0, 1)  # apply limits again, just in case.


# Simple Image
# -----------------
# 
# Use `plt.imshow` to get a quick look at the channels and RGB composite we
# created.
# 
# First, plot each channel individually. The deeper the color means the
# satellite is observing more light in that channel. Clouds appear white because
# they reflect lots of red, green, and blue light. Notice that the land reflects
# a lot of "green" in the veggie channel because this channel is sensitive to
# the chlorophyll.
# 
# 

# In[10]:


fig, ([ax1, ax2, ax3, ax4]) = plt.subplots(1, 4, figsize=(16, 3))

ax1.imshow(R, cmap='Reds', vmax=1, vmin=0)
ax1.set_title('Red', fontweight='bold')
ax1.axis('off')

ax2.imshow(G, cmap='Greens', vmax=1, vmin=0)
ax2.set_title('Veggie', fontweight='bold')
ax2.axis('off')

ax3.imshow(G_true, cmap='Greens', vmax=1, vmin=0)
ax3.set_title('"True" Green', fontweight='bold')
ax3.axis('off')

ax4.imshow(B, cmap='Blues', vmax=1, vmin=0)
ax4.set_title('Blue', fontweight='bold')
ax4.axis('off')

plt.subplots_adjust(wspace=.02)


# The addition of the three channels results in a color image. Combine the three
# channels with a stacked array and display the image with `imshow`.
# 
# 

# In[11]:


# The RGB array with the raw veggie band
RGB_veggie = np.dstack([R, G, B])

# The RGB array for the true color image
RGB = np.dstack([R, G_true, B])

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

# The RGB using the raw veggie band
ax1.imshow(RGB_veggie)
ax1.set_title('GOES-16 RGB Raw Veggie', fontweight='bold', loc='left',
              fontsize=12)
ax1.set_title('{}'.format(scan_start.strftime('%d %B %Y %H:%M UTC ')),
              loc='right')
ax1.axis('off')

# The RGB for the true color image
ax2.imshow(RGB)
ax2.set_title('GOES-16 RGB True Color', fontweight='bold', loc='left',
              fontsize=12)
ax2.set_title('{}'.format(scan_start.strftime('%d %B %Y %H:%M UTC ')),
              loc='right')
ax2.axis('off')


# Adjust Image Contrast
# ---------------------
# 
# I think the color looks a little dull. We could get complicated and make a
# Rayleigh correction to the data to fix the blue light scattering, but that can
# be intense. More simply, we can make the colors pop out by adjusting the image
# contrast. Adjusting image contrast is easy to do in Photoshop, and also easy
# to do in Python.
# 
# We are still using the RGB values from the day/night GOES-16 ABI scan.
# 
# Note: you should adjust the contrast _before_ you add in the Clean IR channel.
# 
# 

# In[12]:


def contrast_correction(color, contrast):
    """
    Modify the contrast of an RGB
    See:
    https://www.dfstudios.co.uk/articles/programming/image-programming-algorithms/image-processing-algorithms-part-5-contrast-adjustment/

    Input:
        color    - an array representing the R, G, and/or B channel
        contrast - contrast correction level
    """
    F = (259*(contrast + 255))/(255.*259-contrast)
    COLOR = F*(color-.5)+.5
    COLOR = np.clip(COLOR, 0, 1)  # Force value limits 0 through 1.
    return COLOR


# Amount of contrast
contrast_amount = 200

# Apply contrast correction
RGB_contrast = contrast_correction(RGB, contrast_amount)


# In[16]:


# Plot on map with Cartopy

fig = plt.figure(figsize=(15, 12))
plt.imshow(RGB_contrast, interpolation="antialiased")
plt.imsave("out/color.png", RGB_contrast)

Makefile

all: earth.png

earth.png: out/color.png 
	convert $< \
		-resize 2972x1820\> \
		-size 3072x1920 xc:black +swap \
		-gravity center \
	  -composite $@

out/color.png: full_disk.py data/full_disk.nc
	python $<

data/full_disk.nc: download.sh
	bash -e $<