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

NOAA’s GOES-16/17 satellite data have been made accessible on Amazon AWS’s S3 object storage. These data are stored on the commonly-used NetCDF4 data format and are easy to access in Python thanks to great tools such as xarray and s3fs.

If you’ve worked with GOES data before though, you’ll know that the data coordinates are on a x/y scan grid, not latitudes and longitudes.

But there’s a (relatively) easy way to convert between the GOES x/y coordinates and latitude/longitude as long as you know a few key parameters about the GOES projection. The full mathematical derivation of the coordinate transform can be found here:

While the guide above goes (pun slightly intended) into great detail on the projections and coordinates of the GOES satellite, I wanted to combine that information with information on how to access the GOES data on AWS, plotting with xarray and cartopy, and put it all in one place.

# Code Access

A follow-along Jupyter notebook can be found here:

A Binder environment is also available to run this code interactively on the cloud without needing to set up anything on your local machine. Either click on the badge in the repo’s README file or this link.

To see the code with all the outputs, see here: https://nbviewer.jupyter.org/github/lsterzinger/goes-l2-latlon-tutorial/blob/main/tutorial-with-outputs.ipynb

# 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.

We’ll also import `cartopy.crs`

and `cartopy.feature`

for mapping and border options on our plot

`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

Let’s open the file. We’ll use s3fs and xarray to open the remote file on S3 and convert it into a dataset that we can use for plotting.

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)

The `ds`

dataframe has coordinate dimensions `x`

and `y`

in radians, which are part of the GOES projection. Converting these to latitude and longitude is useful for selecting subsets of the data over certain geographical areas.

## 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

Now, `ds`

has 2D coordinates `lat`

and `lon`

. Unfortunately, xarray doesn’t support using `ds.sel()`

on non-dimension coordinates, and 2D arrays cannot be added as dimension coordinates.

Instead, we’ll write a small function to return the x/y maxima and minima for a given latitude/longitude box:

`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)))

Now when we pass `ds`

to this function, we get a new dataset with `lat`

and `lon`

as coordinates.

`ds = calc_latlon(ds)`

Now, if we print `ds.coords`

we see the following coordinates appended to the list:

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

These coordinates are 2D arrays of latitude and longitude, with dimensions `(y,x)`

## 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)

Next, use the `get_xy_from_latlon()`

function we defined above to get the `x/y`

max/min for the given latitude and longitude pairs.

`((x1,x2), (y1, y2)) = get_xy_from_latlon(ds, lats, lons)`

Now that we have `x1, x2, y1, y2`

defined, we can use `ds.sel()`

to take a subset of the data

`subset = ds.sel(x=slice(x1, x2), y=slice(y2, y1))`

Note: the y-slice must be done in order `y2, y1`

because the y-coordinate is decreasing from left to right so the larger value must be given first.

## 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)

In this example, we define the projection (Orthographic) meaning that the data will be plotted on a globe. The transform (PlateCarree) tells cartopy that the data is in lat/lon coordinates. This yields the following plot:

We can also build a “fake” true-color image by combining the red, blue, and “veggie” channels (GOES does not capture a green band; the vegetation channel can be used instead, but must be corrected else our image would appear too green. For more information, see here)

# 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))

Now, we can plot the RGB image with:

`fig = plt.figure(figsize=(8,5))`

plt.imshow(rgb)

While constructing an RGB image in matplotlib does not allow for plotting on different coordinate grids (so we won’t be able to get coastlines/state borders), we still plotted the following image based on the lat/lon box defined above. For more advanced satellite plotting, see PyTroll/SatPy.

I hope this guide was useful! If you have questions, comments, or suggestions, please feel free reach out to me on twitter @lucassterzinger or via email (see git profile). If you’d like to add something to the tutorial, please feel free to open a PR on the notebook repository.