PYTHON SCRIPT FOR HIMAWARI-8

Here is my current Python script for Himawari-8. If you followed the Himawari-8 Data Download tutorial and saved the files as suggested you need only change the path variable in the script. For example, if your path looks like this:

D:/himawari-8-data/2020/0415/0500/

The path variable will be this:

path = 'D:/himawari-8-data/'

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.

[IMPORTANT UPDATE: March 22, 2024] The Australian Bureau of Meteorology now serves only full-resolution data, which means the red band is at 500-meter resolution while the blue, green, and infrared bands are at 1000-meter resolution. The script has been updated and now resamples the red band down to 1000-meters / 11000 x 11000 pixels or all bands down to 2000-meters / 5500 x 5500 pixels depending on the size variable. Please note that full-resolution data will be much more memory intensive.

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

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

# Define variables
path = 'D:/himawari-8-data/'
year = '2020' # YYYY
date = '0415' # MMDD
time = '0500' # hhmm
mode = 'OBS'  # Either OBS or BRF
size = '2000' # Either 2000 or 1000

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

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

# Variables
if mode == 'OBS':
    v_mode = 'scaled_radiance'

if mode == 'BRF':
    v_mode = 'brf'

full_path = path + year + '/' + date + '/' + time + '/'
date_name = year + date + time
save_file = 'Himawari-8_' + date_name + '_' + mode + '_' + size + '.png'

# Himawari 8 files
nc1_file = date_name + '00-P1S-ABOM_' + mode + '_B01-PRJ_GEOS141_' + '1000' + '-HIMAWARI8-AHI.nc'
nc2_file = date_name + '00-P1S-ABOM_' + mode + '_B02-PRJ_GEOS141_' + '1000' + '-HIMAWARI8-AHI.nc'
nc3_file = date_name + '00-P1S-ABOM_' + mode + '_B03-PRJ_GEOS141_' + '500'  + '-HIMAWARI8-AHI.nc'
nc4_file = date_name + '00-P1S-ABOM_' + mode + '_B04-PRJ_GEOS141_' + '1000' + '-HIMAWARI8-AHI.nc'

# Read the files using NetCDF4 with the 'r' flag (read)
nc1 = Dataset(full_path + nc1_file, 'r')
nc2 = Dataset(full_path + nc2_file, 'r')
nc3 = Dataset(full_path + nc3_file, 'r')
nc4 = Dataset(full_path + nc4_file, 'r')

# Load the array variables
B = nc1.variables['channel_0001_' + v_mode][0,:,:] # Band 1 is blue (0.47 um)
G = nc2.variables['channel_0002_' + v_mode][0,:,:] # Band 2 is green (0.51 um)
R = nc3.variables['channel_0003_' + v_mode][0,:,:] # Band 3 is red (0.64 um)
I = nc4.variables['channel_0004_' + v_mode][0,:,:] # Band 4 is "veggie" infrared (0.865 um)

# Close the files
nc1.close()
nc2.close()
nc3.close()
nc4.close()

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

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

# Change all 1 values to 0.999999
# This will allow the white background to be changed to
# black later without affecting any pixels in the image
if mode == 'BRF':
    B[B == 1] = 0.999999
    G[G == 1] = 0.999999
    R[R == 1] = 0.999999
    I[I == 1] = 0.999999

# Resample the channel(s) to the desired resolution
def rebin(a, shape):
    sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]
    return a.reshape(sh).mean(-1).mean(1)

if size == '2000':
    B = rebin(B, (5500, 5500)).astype(np.float32)
    G = rebin(G, (5500, 5500)).astype(np.float32)
    R = rebin(R, (5500, 5500)).astype(np.float32)
    I = rebin(I, (5500, 5500)).astype(np.float32)
elif size == '1000':
    R = rebin(R, (11000, 11000)).astype(np.float32)
else:
    print(f"The size variable {size} is invalid. It must be either '1000' or '2000'.")

# Blend the G channel and R channel to shift the hue
# This will compensate for the non-ideal G wavelength
G_hue = cv2.addWeighted(G, 0.66, R, 0.34, 0)

# Use the I channel to lighten G_hue using np.maximum
# This will increase the vegetation signal
G_veg = new_G = cv2.addWeighted(G_hue, 0.97, np.maximum(G_hue, I), 0.03, 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
if mode == 'OBS':
    a = 1.85
    b = 0.90
    c = 1.90

if mode == 'BRF':
    a = 1.25
    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
# Useful for BRF images
if mode == 'BRF':
    img[img == 1] = 0

# Create a mask to hide overscan regions
# Useful for OBS images
if mode == 'OBS':
    # Return image width and height
    w, h = img.shape[:2]

    # Return image center
    center = (int(w / 2), int(h / 2))

    # Ratio of Earth's diameter at poles to the equator
    ratio = 0.9966472

    # Overscan margin
    margin = w * 0.002363636

    # Radius at equator
    radius_w = int((w - margin) / 2)

    # Radius at poles
    radius_h = int((w - margin) / 2 * ratio)

    # Axes
    axes = (radius_w, radius_h)

    # Create an elliptical mask
    mask = np.zeros(img.shape, dtype = np.uint8)
    mask_ellipse = cv2.ellipse(mask, center, axes, 0, 0, 360, (255, 255, 255), -1, cv2.LINE_AA)
 
    # Apply the mask
    mask_img = img * (mask_ellipse / 255)

    # Final image
    img = mask_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 Himawari-8 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.

Himawari-8: Satellite image originally processed by the Bureau of Meteorology from the geostationary satellite Himawari-8 operated by the Japan Meteorological Agency.

Not bad! But it needs more work and an update is coming soon.

Please note that this is not the same technique used in my earlier post, Himawari-8 Color Correction. That was done using pre-processed PNG images. This is the actual satellite data.

Eventually you may want to explore different dates, resolutions, or modes. You will need to make the appropriate changes to the script variables.

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.

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

Australian Bureau of Meteorology http://www.bom.gov.au/metadata/19115/ANZCW0503900400

Japan Meteorological Agency

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