PYTHON SCRIPT FOR GOES SATELLITES

Here is my current Python script for GOES-16 and GOES-17 satellites. If you followed the GOES Satellite Data Download tutorial and saved the file as suggested you need only change the path variable in the script. For example, if your path looks like this:

E:/goes-16-data/full-disk/multi-band/2020/0415/1430/

The path variable will be this:

path = 'E:/goes-16-data/full-disk/multi-band/'

Do not include the YYYY/MMDD/hhmm/ part of your path; the script assumes you used this directory structure and will handle it automatically. Also don't forget the quotes and forward slash at the end!

You should already have Python installed as explained in Installing Python on Windows 10. Run Spyder from the Start menu and copy and paste the text below directly into the script pane and change the path.

###############################
#    EDIT THESE VARIABLES     #
###############################

# Tutorial: https://loneskyimages.blogspot.com/2020/04/goes-satellite-data-download.html
# IMPORTANT! Change the path variable
# to where you saved the GOES file
# If you are following the tutorial you
# don't need to make any further changes

# Define variables
path = 'E:/goes-16-data/full-disk/multi-band/'
year = '2020' # YYYY
date = '0415' # MMDD
time = '1430' # hhmm

###############################
# DO NOT EDIT BELOW THIS LINE #
###############################

# Load Python libraries
from glob import glob
from netCDF4 import Dataset
import numpy as np
import cv2

# Variables
full_path = path + year + '/' + date + '/' + time + '/'
date_name = year + date + time
goes_file = str(glob(full_path + '*.nc')).split('\\')[-1][:-2]
satellite = goes_file.split('/')[-1][21:23]
save_file = 'GOES-' + satellite + '_' + date_name + '.png'

# GOES file
nc_file = full_path + goes_file

# Read the file using NetCDF4 with the 'r' flag (read)
nc = Dataset(nc_file, 'r')

# Load the array variables
B = nc.variables['CMI_C01'][:] # Band 1 is blue (0.47 um)
R = nc.variables['CMI_C02'][:] # Band 2 is red (0.64 um)
I = nc.variables['CMI_C03'][:] # Band 3 is "veggie" infrared (0.865 um)

# Close the file
nc.close()

# Turn empty values into NaNs
B[B == -1] = np.nan
R[R == -1] = np.nan
I[I == -1] = np.nan

# Blend the B channel and R channel
G_base = cv2.addWeighted(B, 0.38, R, 0.62, 0)

# Use the I channel to lighten G_base using np.maximum
# This will increase the vegetation signal
G_veg = cv2.addWeighted(G_base, 0.97, np.maximum(G_base, I), 0.03, 0)

# Blend the B channel to decrease the signal
# in predominantly red/orange soil regions
G_soil = new_G = cv2.addWeighted(G_veg, 0.41, np.minimum(G_veg, B), 0.59, 0)

# Create a "true color" BGR image
img = cv2.merge((B, new_G, R))

# Convert BGR to HLS
hls = cv2.cvtColor(img, cv2.COLOR_BGR2HLS)

# Split image to HLS components
h, l, s = cv2.split(hls)

# Adjust the tone curve
# function: y = (x^b(x+1-xa))^1/c
# https://www.desmos.com/calculator/wf7oex5edi
a = 1.85
b = 0.90
c = 1.90

new_l = np.power(np.power(l, b * (l + 1 - (l * a))), 1/c)

# Merge to HLS
new_img = cv2.merge((h, new_l, s))

# Convert to BGR
img = cv2.cvtColor(new_img, cv2.COLOR_HLS2BGR)

# Change the background from white to black
img[img == -1] = 0

# Final image
img = img

# Display the image in a window
cv2.namedWindow('Image', cv2.WINDOW_NORMAL)
cv2.resizeWindow('Image', 1000, 1000)
cv2.imshow('Image', img)

# Save the image
cv2.imwrite(full_path + save_file, img * 255)

# Close the window in the specified milliseconds e.g. (1000)
# or by pressing any key
cv2.waitKey(0)
cv2.destroyAllWindows()

Click the green Run file button near the top or press F5 to run the script.

Paste the code into the Spyder script pane as shown and run.

If all goes well a GOES-16 satellite image from April 15, 2020, will appear in a separate window on your Task bar and the file will be saved to your data folder.

GOES-16 using a new experimental formula to generate a pseudo green band.

While likely an improvement over previous simple blending formulas for generating a first-order synthetic green band, this will never be able to fully reproduce the subtle, sublime depth of color visible in true color imagery.

Eventually you may want to explore different dates, resolutions, or modes. You will need to make the appropriate changes to the script variables. I will be expanding this tutorial in the near future to cover at the very least individual bands and higher resolutions.

CAVEATS
When you examine the image you may see there is quite a lot of noise beyond the terminator. I'm still working on a way to gradually darken those pixels. In fact the entire script is still a work in progress.

Lastly, coding and math are not my strong points so I welcome any suggestions to improve and optimize the script. Any future updates and tweaks will be noted below so check back occasionally.

[UPDATE: May 4, 2020] Changed the method for making a black background (greatly improves the appearance of the mask); removed the # on the cv2.imwrite() line that was erroneously left in.

CREDIT
Parts of this script are based on examples by Brian Blaylock and others.

NASA
NOAA

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.