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

Code Access

A follow-along Jupyter notebook can be found here:

Let’s go!

First, import the necessary libraries. We’ll be opening GOES data on Amazon Web Services’ S3 storage, so we need s3fs as well as xarray to open the data. We also need numpy for some extra array operations and matplotlib.pyplot for plotting.

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

This can be accomplished using a bit of math, which I got from this blog post. The math is dependent on some projection parameters, which are stored in the variable goes_imager_projection . The function below uses those parameters, along with the transform math, to return a dataset with lat and lon as xarray coordinates.

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

Adding lat/lon as a coordinate allows for easy taking of a subset of the data. First, define a bounding box in latitude and longitude:

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

We can use cartopy along with xarray’s built-in plotting tools (built on matplotlib) to plot the imagery on a map. I chose the blue channel (“CMI_C01”) to plot.

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)

--

--

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

19 Followers

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