Howdy,
I’m working with the below script. I got to work truecolor, but I want it to do truecolor and Shortwave IR for the night. It only seems to be doing day and not IR. So at night it is black. Does anyone have experience with this who might be able to assist? Thank you
import numpy as np
import os
import s3fs
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from datetime import datetime, timedelta
# Set up the connection to the AWS S3 bucket
fs = s3fs.S3FileSystem(anon=True)
# Define the save path and the maximum number of images to keep
save_path = '/var/www/site/static/images/floaters/GM/truecolor'
max_images = 30
# Define the path to the GOES-16 ABI bands (true color requires Red, Green, Blue bands)
year_path = 'noaa-goes16/ABI-L2-MCMIPC/2024/'
days = fs.ls(year_path)
latest_day = sorted(days)[-1]
time_segments = fs.ls(latest_day)
latest_time_segment = sorted(time_segments)[-1]
nc_files = fs.ls(latest_time_segment)
latest_nc_file = sorted([f for f in nc_files if f.endswith('.nc')])[-1]
# Extract the timestamp from the latest filename (in the format "sYYYYJJJHHMMSS")
timestamp_str = latest_nc_file.split("_")[4][1:14] # Extract only the first 14 characters after 's'
dt = datetime.strptime(timestamp_str, "%Y%j%H%M%S") # Parse the datetime (YYYY + Julian Day + HHMMSS)
# Round the timestamp to the nearest 5 minutes
minute = (dt.minute + 4) // 5 * 5
dt = dt.replace(minute=minute, second=0)
timestamp_display = dt.strftime("%H:%MZ %b %d, %Y") # Format as "HH:MMZ %b %d, %Y"
# Set the filename format and check if the image already exists
filename = f"goes16-gm-truecolor-{dt.strftime('%H%MZ-%b-%d-%Y')}.jpeg"
output_path = os.path.join(save_path, filename)
# Skip if the image already exists
if os.path.exists(output_path):
print(f"Image {filename} already exists. Skipping this run.")
exit()
# Open the latest file directly from the S3 bucket
with fs.open(latest_nc_file) as file:
ds = xr.open_dataset(file)
# Load the Red, Green, and Blue channels (C01 = Blue, C02 = Red, C03 = Veggie band)
R = ds['CMI_C02'].data # Red channel
G = ds['CMI_C03'].data # Green channel (veggie band)
B = ds['CMI_C01'].data # Blue channel
# Ensure values are between 0 and 1
R = np.clip(R, 0, 1)
G = np.clip(G, 0, 1)
B = np.clip(B, 0, 1)
# Apply gamma correction
gamma = 2.2
R = np.power(R, 1/gamma)
G = np.power(G, 1/gamma)
B = np.power(B, 1/gamma)
# Calculate the "True" Green channel
G_true = 0.45 * R + 0.1 * G + 0.45 * B
G_true = np.clip(G_true, 0, 1)
# Stack the corrected RGB channels into a single image
RGB = np.dstack([R, G_true, B])
# Set up satellite projection parameters
sat_height = ds['goes_imager_projection'].perspective_point_height
sat_lon = ds['goes_imager_projection'].longitude_of_projection_origin
# Define the extent using x and y coordinates
x = ds['x'] * sat_height
y = ds['y'] * sat_height
extent = [x.min(), x.max(), y.min(), y.max()] # Set extent based on satellite coordinates
# Create a figure with a map projection
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
# Apply the extent for the region specified
ax.set_extent([-100, -79, 17, 32], crs=ccrs.PlateCarree())
# Plot the RGB image with the correct transformation and extent
ax.imshow(RGB, origin='upper', extent=extent,
transform=ccrs.Geostationary(central_longitude=sat_lon, satellite_height=sat_height))
# Add coastlines and state boundaries for reference
ax.coastlines(resolution='10m', color='black', linewidth=1)
ax.add_feature(cfeature.BORDERS, edgecolor='black')
ax.add_feature(cfeature.STATES, linewidth=0.25)
# Customize gridlines: Show only left and bottom labels
gl = ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False # Disable top labels
gl.right_labels = False # Disable right labels
gl.left_labels = True # Enable left labels
gl.bottom_labels = True # Enable bottom labels
# Add the title and timestamp on a single line, slightly moving the text down
title_text = f"GOES-16 True Color (Day) and Shortwave IR (Night) - {timestamp_display}"
ax.text(0.02, 1.03, title_text, transform=ax.transAxes, fontsize=9, fontweight='bold', ha='left', va='top', color='black')
# Add "NAME.COM" in the top-right corner
ax.text(0.98, 1.03, 'NAME.COM', transform=ax.transAxes, fontsize=10, fontweight='bold', ha='right', va='top', color='#A9A9A9')
# Save the image as JPEG
plt.savefig(output_path, format='jpeg', dpi=300, bbox_inches='tight')
# Close the figure
plt.close(fig)
print(f"True color image saved at {output_path}")
# Step 6: Keep only the latest 30 images in the directory
images = sorted([f for f in os.listdir(save_path) if f.endswith('.jpeg')])
if len(images) > max_images:
# Delete the oldest images
oldest_images = images[:len(images) - max_images]
for img in oldest_images:
os.remove(os.path.join(save_path, img))
print(f"Deleted old image: {img}")