# GOES-R Satellite Latitude and Longitude Grid Projection Algorithm

In this tutorial, I will discuss the methods used for calculation of latitude and longitude from a GOES-R L1b data file. The GOES-R L1b radiance files contain radiance data and geometry scan information in radians. This information is not enough to plot geographic radiance data right from the file, however, after some geometric manipulation harnessing satellite position and ellipsoid parameters, we can derive latitude and longitude values from the one-dimensional scan angles and plot our data in projected formats familiar to many geographic information tools.

#### Blog Post Index:

__Theory: Satellite and Earth Geometry__

Understanding Scan Angles (Radians)

Reprojecting Scan Angles (Radians) to Geographic Coordinates (Lat/Lon)

Exploring GOES-R netCDF Data Files

Algorithm Implementation: GOES-R Scan Angles (Radians) to Geographic Coordinates (Lat/Lon)

Visualizing GOES-R Data with The Derived Geographic Coordinate Grid Projection

Conclusion

Understanding Scan Angles (Radians)

Reprojecting Scan Angles (Radians) to Geographic Coordinates (Lat/Lon)

Exploring GOES-R netCDF Data Files

Algorithm Implementation: GOES-R Scan Angles (Radians) to Geographic Coordinates (Lat/Lon)

Visualizing GOES-R Data with The Derived Geographic Coordinate Grid Projection

Conclusion

Satellite pixels can be approximated as an area on the surface of the earth. Each pixel can be identified by an angle at a certain height above Earth. The satellite produces angles in the x and y direction, while also giving us information about the height of the instrument and information about the earth. These parameters will allows us to approximate the latitude and longitude of each pixel after extrapolation onto the earth’s surface.

Each pixel needs information about the follow:

Geometry of the earth ellipsoid

Reference with respect to the earth’s center

Satellite properties from the point and center of Earth

Below is a depiction of the complex geometry involved in projecting the scan information from the satellite onto the earth’s geodetic coordinate system:

Now, we will derive the equations needed to go from the scan angles (in radians) to the latitude and longitude coordinates that we often associate with geographic information.

The GOES-R satellite produces netCDF4 data files that contain data and scan variables. The scan variables are called ‘x’ and ‘y’ GOES fixed grid projection coordinates in radians:

The ‘x’ and ‘y’ projection coordinates merely give information about the angle of each pixel in reference to the imager. This looks like the following:

Using the scan angles and the geometry above, we can derive the latitude and longitude grid to relate the data variables to ground-truth behavior. This method creates latitude and longitude arrays, which also aid in the calibration of satellite products, as well as integration with weather models based on geographic locations.

It is easy to think of the relationship above as one-to-one, but in reality, we need both scan angles to create latitude and longitude.

The scan angles can be gridded using the simple square array above, but then need to be converted using the complex geometry shown in the first figure above. In the next section, I outline the equations needed to transform the radian angles taken from the netCDF file and convert them to meaningful geographic latitude and longitude values. The reprojection will create a 2-D grid that we can then use to plot the data from each data file.

Now that we have a rough understanding of satellite-earth geometry and the gridded scheme used in the GOES-R data files, we can introduce the complex series of equations used to reproject the scan angles down to the earth’s surface in terms of latitude and longitude coordinates.

From here, *x* represents the horizontal scan value at a given point, and *y* represents a vertical scan value at a given point. *φ* is the geodetic latitude of a corresponding point on the earth's surface, and *λ* is the geodetic longitude of the corresponding point. Below is the equation used to convert scan angles (*x, y*) to latitude and longitude (*φ, λ*):

where:

I will not cover the derivation of each coordinate rotation (*s _{x,y,z}*) here, however, the methods can be found in nearly any orbital mechanics textbook. The coordinate rotations are defined as follows:

Again, the *x, y* are the scan angles given by the netCDF data file. From here, we further define:

Which was found using the quadratic equation and the geometry above. The variables are defined as:

And now we’re ready to calculate our latitude and longitude grid from the information give to us by the GOES-R satellite data files.

The GOES-R files contain information needed to use the equations above. GOES-R data can be found on the Google Cloud server located here:

https://console.cloud.google.com/storage/browser/gcp-public-data-goes-16

Download a test file from one of the folders (I recommend using ** ‘gcp-public-data-goes-16/ABI-L1b-RadC/2018’ **which is recent raw radiance data from 2018. Once the file is downloaded and put into a nearby directory (I put mine into ‘./rad_nc_files/’) we can start analzing the .nc file.

The following code reads a netCDF file and creates a projection variables called ‘proj_info’ which contains many of the variables needed in the equations above:

#!/usr/bin/python from netCDF4 import Dataset import os # navigate to directory with .nc data files (raw radiance files in my case) os.chdir('./rad_nc_files/') full_direc = os.listdir() nc_files = [ii for ii in full_direc if ii.endswith('.nc')] nc_indx = 5 g16_data_file = nc_files[nc_indx] # select .nc file print(nc_files[nc_indx]) # print file name # designate dataset g16nc = Dataset(g16_data_file, 'r') # GOES-R projection info proj_info = g16nc.variables['goes_imager_projection']

Printing out the ‘proj_info’ variable yields the following:

<class 'netCDF4._netCDF4.Variable'> int32 goes_imager_projection() long_name: GOES-R ABI fixed grid projection grid_mapping_name: geostationary perspective_point_height: 35786023.0 semi_major_axis: 6378137.0 semi_minor_axis: 6356752.31414 inverse_flattening: 298.2572221 latitude_of_projection_origin: 0.0 longitude_of_projection_origin: -75.0 sweep_angle_axis: x unlimited dimensions: current shape = () filling on, default _FillValue of -2147483647 used

There are several very important variables to be taken from this list:

With the variables above, we now have enough information to go from the data files that contain 1-D radian angles to latitude and longitude grids. In the next section I will introduce the actual algorithm for calculating lat/lon grids from radian angles.

Now that all the groundwork has been built, we can calculate the actual grid for a given netCDF data file. The Python program is introduced below:

#!/usr/bin/python from netCDF4 import Dataset import numpy as np import os # navigate to directory with .nc data files (raw radiance files in my case) os.chdir('./rad_nc_files/') full_direc = os.listdir() nc_files = [ii for ii in full_direc if ii.endswith('.nc')] nc_indx = 5 g16_data_file = nc_files[nc_indx] # select .nc file print(nc_files[nc_indx]) # print file name # designate dataset g16nc = Dataset(g16_data_file, 'r') # GOES-R projection info and retrieving relevant constants proj_info = g16nc.variables['goes_imager_projection'] lon_origin = proj_info.longitude_of_projection_origin H = proj_info.perspective_point_height+proj_info.semi_major_axis r_eq = proj_info.semi_major_axis r_pol = proj_info.semi_minor_axis # Data info lat_rad_1d = g16nc.variables['x'][:] lon_rad_1d = g16nc.variables['y'][:] # close file when finished g16nc.close() g16nc = None # create meshgrid filled with radian angles lat_rad,lon_rad = np.meshgrid(lat_rad_1d,lon_rad_1d) # lat/lon calc routine from satellite radian angle vectors lambda_0 = (lon_origin*np.pi)/180.0 a_var = np.power(np.sin(lat_rad),2.0) + (np.power(np.cos(lat_rad),2.0)*(np.power(np.cos(lon_rad),2.0)+(((r_eq*r_eq)/(r_pol*r_pol))*np.power(np.sin(lon_rad),2.0)))) b_var = -2.0*H*np.cos(lat_rad)*np.cos(lon_rad) c_var = (H**2.0)-(r_eq**2.0) r_s = (-1.0*b_var - np.sqrt((b_var**2)-(4.0*a_var*c_var)))/(2.0*a_var) s_x = r_s*np.cos(lat_rad)*np.cos(lon_rad) s_y = - r_s*np.sin(lat_rad) s_z = r_s*np.cos(lat_rad)*np.sin(lon_rad) lat = (180.0/np.pi)*(np.arctan(((r_eq*r_eq)/(r_pol*r_pol))*((s_z/np.sqrt(((H-s_x)*(H-s_x))+(s_y*s_y)))))) lon = (lambda_0 - np.arctan(s_y/(H-s_x)))*(180.0/np.pi) # print test coordinates print('{} N, {} W'.format(lat[318,1849],abs(lon[318,1849])))

If your algorithm is correct and your data file is 1500 x 2500, then your output above should look like the following:

The screenshot above tells the user a few things:

The name of the netCDF file read in

There were negative values under the square root in the calculation of r_s

There were negative values under the square root in the calculation of latitude

The values for latitude and longitude at a specific point were successful

If your screenshot looks similar to the one above - CONGRATULATIONS. The code worked. The values outputted by the program are the latitude and longitude of the pixel nearest to New York City. Now, you have a grid for latitude and longitude that is the same size as the data array in the netCDF file. In the next section we will use this to plot and visualize the data!

Below is a simple implementation of the grid projection algorithm above, which takes the grid and projects the data onto a basemap using Python’s Basemap Toolkit.

#!/usr/bin/python from netCDF4 import Dataset import matplotlib as mpl mpl.use('TkAgg') import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap, cm import numpy as np import os def lat_lon_reproj(nc_folder,nc_indx): os.chdir(nc_folder) full_direc = os.listdir() nc_files = [ii for ii in full_direc if ii.endswith('.nc')] g16_data_file = nc_files[nc_indx] # select .nc file print(nc_files[nc_indx]) # print file name # designate dataset g16nc = Dataset(g16_data_file, 'r') var_names = [ii for ii in g16nc.variables] var_name = var_names[0] try: band_id = g16nc.variables['band_id'][:] band_id = ' (Band: {},'.format(band_id[0]) band_wavelength = g16nc.variables['band_wavelength'] band_wavelength_units = band_wavelength.units band_wavelength_units = '{})'.format(band_wavelength_units) band_wavelength = ' {0:.2f} '.format(band_wavelength[:][0]) print('Band ID: {}'.format(band_id)) print('Band Wavelength: {} {}'.format(band_wavelength,band_wavelength_units)) except: band_id = '' band_wavelength = '' band_wavelength_units = '' # GOES-R projection info and retrieving relevant constants proj_info = g16nc.variables['goes_imager_projection'] lon_origin = proj_info.longitude_of_projection_origin H = proj_info.perspective_point_height+proj_info.semi_major_axis r_eq = proj_info.semi_major_axis r_pol = proj_info.semi_minor_axis # grid info lat_rad_1d = g16nc.variables['x'][:] lon_rad_1d = g16nc.variables['y'][:] # data info data = g16nc.variables[var_name][:] data_units = g16nc.variables[var_name].units data_time_grab = ((g16nc.time_coverage_end).replace('T',' ')).replace('Z','') data_long_name = g16nc.variables[var_name].long_name # close file when finished g16nc.close() g16nc = None # create meshgrid filled with radian angles lat_rad,lon_rad = np.meshgrid(lat_rad_1d,lon_rad_1d) # lat/lon calc routine from satellite radian angle vectors lambda_0 = (lon_origin*np.pi)/180.0 a_var = np.power(np.sin(lat_rad),2.0) + (np.power(np.cos(lat_rad),2.0)*(np.power(np.cos(lon_rad),2.0)+(((r_eq*r_eq)/(r_pol*r_pol))*np.power(np.sin(lon_rad),2.0)))) b_var = -2.0*H*np.cos(lat_rad)*np.cos(lon_rad) c_var = (H**2.0)-(r_eq**2.0) r_s = (-1.0*b_var - np.sqrt((b_var**2)-(4.0*a_var*c_var)))/(2.0*a_var) s_x = r_s*np.cos(lat_rad)*np.cos(lon_rad) s_y = - r_s*np.sin(lat_rad) s_z = r_s*np.cos(lat_rad)*np.sin(lon_rad) lat = (180.0/np.pi)*(np.arctan(((r_eq*r_eq)/(r_pol*r_pol))*((s_z/np.sqrt(((H-s_x)*(H-s_x))+(s_y*s_y)))))) lon = (lambda_0 - np.arctan(s_y/(H-s_x)))*(180.0/np.pi) # print test coordinates print('{} N, {} W'.format(lat[318,1849],abs(lon[318,1849]))) return lon,lat,data,data_units,data_time_grab,data_long_name,band_id,band_wavelength,band_wavelength_units,var_name nc_folder = './rad_nc_files/' # define folder where .nc files are located file_indx = 15 # be sure to pick the correct file. Make sure the file is not too big either, # some of the bands create large files (stick to band 7-16) # main data grab from function above lon,lat,data,data_units,data_time_grab,data_long_name,band_id,band_wavelength,band_units,var_name = lat_lon_reproj(nc_folder,file_indx) bbox = [np.min(lon),np.min(lat),np.max(lon),np.max(lat)] # set bounds for plotting # figure routine for visualization fig = plt.figure(figsize=(9,4),dpi=200) n_add = 0 m = Basemap(llcrnrlon=bbox[0]-n_add,llcrnrlat=bbox[1]-n_add,urcrnrlon=bbox[2]+n_add,urcrnrlat=bbox[3]+n_add,resolution='i', projection='cyl') m.drawcoastlines(linewidth=0.5) m.drawcountries(linewidth=0.25) m.pcolormesh(lon.data, lat.data, data, latlon=True) parallels = np.linspace(np.min(lat),np.max(lat),5.) m.drawparallels(parallels,labels=[True,False,False,False]) meridians = np.linspace(np.min(lon),np.max(lon),5.) m.drawmeridians(meridians,labels=[False,False,False,True]) cb = m.colorbar() data_units = ((data_units.replace('-','^{-')).replace('1','1}')).replace('2','2}') plt.rc('text', usetex=True) cb.set_label(r'%s $ \left[ \mathrm{%s} \right] $'% (var_name,data_units)) plt.title('{0}{2}{3}{4} on {1}'.format(data_long_name,data_time_grab,band_id,band_wavelength,band_units)) plt.tight_layout() #plt.savefig('goes_16_demo.png',dpi=200,transparent=True) # uncomment to save figure plt.show()

If analyzing raw radiance data, the figure below is the expected output. It has information about the data (ABI L1b Radiances), the band ID (8th band, the band wavelength (6.118 microns), and the date (July 7th, 2018).

In contrast to the figure above, the figure below is taken from the GOES-R processed land surface temperature (LST) product. The LST product is used for applications in weather modeling and prediction, real-time information on fires and other natural hazards, historical behavior in relation to crops, and much more. A snapshot taken from the 6th of June, 2018 is shown below:

Notice how the figure directly above masks the off-land data. This is how many algorithms handle product data - they filter out data that may be over water, typically using some sort of elevation mask. There are also other masks that filter cloudy data and other low quality data that may have interference between the satellite and the surface of the earth.

This concludes the tutorial on creating latitude and longitude grids from netCDF files that only contain 1-dimensional scan angles. This is a great practice in understanding satellites and orbital mechanics. In the tutorial I covered some basics on orbiting satellites and scanning methods used by satellites. While I didn’t derive the complex geometric equations, I demonstrated how to use them and implement the algorithms used for plotting satellite imagery in Python. Python was able to handle reading and plotting of data fairly quickly, which makes this a powerful method for plotting satellite imagery. One may want to create a grid file and save it, so that the algorithm can be sped up when reading geostationary satellite imagery. However, for an introduction into orbiting satellites and reading netCDF data, while visualizing imagery, it’s a great foray into visualizing radiances from satellites.

*See More in Engineering and GIS:*