Add Lat/Lon Coordinates to GOES-16/GOES-17 L2 Data and Plot with Cartopy

Code Access

Let’s go!

import xarray as xr 
import s3fs
import matplotlib.pyplot as plt
import numpy as np

import cartopy.crs as ccrs
import cartopy.feature as cfeature
fs = s3fs.S3FileSystem(anon=True)f = fs.open("s3://noaa-goes17/ABI-L2-MCMIPC/2021/050/18/OR_ABI-L2-MCMIPC-M6_G17_s20210501801176_e20210501803549_c20210501804089.nc")ds = xr.open_dataset(f) 

Do the conversion

def calc_latlon(ds):
# The math for this function was taken from
# https://makersportal.com/blog/2018/11/25/goes-r-satellite-latitude-and-longitude-grid-projection-algorithm
x = ds.x
y = ds.y
goes_imager_projection = ds.goes_imager_projection

x,y = np.meshgrid(x,y)

r_eq = goes_imager_projection.attrs["semi_major_axis"]
r_pol = goes_imager_projection.attrs["semi_minor_axis"]
l_0 = goes_imager_projection.attrs["longitude_of_projection_origin"] * (np.pi/180)
h_sat = goes_imager_projection.attrs["perspective_point_height"]
H = r_eq + h_sat

a = np.sin(x)**2 + (np.cos(x)**2 * (np.cos(y)**2 + (r_eq**2 / r_pol**2) * np.sin(y)**2))
b = -2 * H * np.cos(x) * np.cos(y)
c = H**2 - r_eq**2

r_s = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)

s_x = r_s * np.cos(x) * np.cos(y)
s_y = -r_s * np.sin(x)
s_z = r_s * np.cos(x) * np.sin(y)

lat = np.arctan((r_eq**2 / r_pol**2) * (s_z / np.sqrt((H-s_x)**2 +s_y**2))) * (180/np.pi)
lon = (l_0 - np.arctan(s_y / (H-s_x))) * (180/np.pi)

ds = ds.assign_coords({
"lat":(["y","x"],lat),
"lon":(["y","x"],lon)
})
ds.lat.attrs["units"] = "degrees_north"
ds.lon.attrs["units"] = "degrees_east"
return ds
def get_xy_from_latlon(ds, lats, lons):
lat1, lat2 = lats
lon1, lon2 = lons

lat = ds.lat.data
lon = ds.lon.data

x = ds.x.data
y = ds.y.data

x,y = np.meshgrid(x,y)

x = x[(lat >= lat1) & (lat <= lat2) & (lon >= lon1) & (lon <= lon2)]
y = y[(lat >= lat1) & (lat <= lat2) & (lon >= lon1) & (lon <= lon2)]

return ((min(x), max(x)), (min(y), max(y)))
ds = calc_latlon(ds)
lat                  (y, x) float64 53.5 53.49 53.49 ... 14.8 14.8 14.81lon                  (y, x) float64 -184.4 -184.3 -184.2 ... -112.5 -112.4

Take a lat/lon subset

lats = (30, 55)
lons = (-152, -112)
((x1,x2), (y1, y2)) = get_xy_from_latlon(ds, lats, lons)
subset = ds.sel(x=slice(x1, x2), y=slice(y2, y1))

Plot with Cartopy

p = subset.CMI_C01.plot(
x='lon', y='lat',
subplot_kws={'projection' : ccrs.Orthographic(-130, 40)},
transform = ccrs.PlateCarree()
)
p.axes.add_feature(cfeature.COASTLINE)
p.axes.add_feature(cfeature.STATES)
GOES blue channel plotted over the west coast of the United States
# Get channels 
r = subset['CMI_C02'].data; r = np.clip(r, 0, 1)
g = subset['CMI_C03'].data; g = np.clip(g, 0, 1)
b = subset['CMI_C01'].data; b = np.clip(b, 0, 1)
# Apply a gamma of 2.5
gamma = 2.5; r = np.power(r, 1/gamma); g = np.power(g, 1/gamma); b = np.power(b, 1/gamma)
# Calculate "true green" from the veggie channel
g_true = 0.45 * r + 0.1 * g + 0.45 * b
g_true = np.clip(g_true, 0, 1)
rgb = np.dstack((r, g_true, b))
fig = plt.figure(figsize=(8,5))
plt.imshow(rgb)

--

--

--

I'm a PhD student at UC Davis studying clouds, aerosols, and numerical modeling. Follow me on twitter @lucassterzinger

Love podcasts or audiobooks? Learn on the go with our new app.

Recommended from Medium

My Transition Story

4 Less-Known Yet Very Functional Pandas Operations

Knowledge Graph

Python vs Microsoft Excel: The Better Option for Data Analysis

Frazzled This Holiday Season? Check Out This List Of Calming Essential Oils

Altair: Statistical Visualization Library for Python

Predict the Gender and Age range of an individual using Python

What is EDA?

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store
Lucas Sterzinger

Lucas Sterzinger

I'm a PhD student at UC Davis studying clouds, aerosols, and numerical modeling. Follow me on twitter @lucassterzinger

More from Medium

Make Pivot Tables in Python with Mito

GIS project with GeoPandas

Mapping Data Using Folium