I have been working with satellite data from the GOES-16 satellite for quite some time now. A sample .nc file for this data is located at this NOAA AWS link

I use gdalwarp to transform it for Mapbox use (GeoTIFF) and it works great.

I recently started working on some new data from the same satellite, but the projection data seems to not be applied. This is because the data is cutting-edge and is generated by a Python script that does not project the data the same way as the AWS files.

Below is a link to Dropbox to the .nc file I am working on now :


When I run gdalinfo on the new data, I see this :

Driver: netCDF/Network Common Data Format
Files: glm.nc
Size is 1499, 2499
goes_imager_projection#long_name=GOES-R ABI fixed grid projection
total_energy#long_name=Total radiant energy
total_energy#units=J per flash
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 2499.0)
Upper Right ( 1499.0,    0.0)
Lower Right ( 1499.0, 2499.0)
Center      (  749.5, 1249.5)
Band 1 Block=1499x1 Type=Float32, ColorInterp=Undefined
NoData Value=-9999
Unit Type: J per flash
long_name=Total radiant energy
units=J per flash

How do I convert the projection into the same one from the AWS data? (first link).

I found a Python script online which I have put up on Pastebin (https://pastebin.com/P9edec4H) that seems to have all of the projection data needed. However, after trying to get GDAL working in Python for 2 hours I had to give up due to package conflicts. (I use the command line executables normally)

I have tried all sorts of commands (trying to add the proj string, etc) but honestly, it would be a bit embarrassing if I had to post them all here because I don't exactly know what I am doing.

Does anyone have any idea?

  • Your pastebin link isn't working for me, so not sure what you've tried, but I did see an SO post with a similar problem with links to how to do this with cartopy
    – Ryan
    Commented Dec 20, 2019 at 18:06
  • for gdalwarp use -s_srs, similarly in rasterio
    – mdsumner
    Commented Apr 14, 2022 at 1:23

1 Answer 1


It appears that the projection information is stored using CF conventions.

import rioxarray # for 'rio' accessor
import xarray

xds = xarray.open_dataset("glm.nc")

Here is what is in xds:

Dimensions:                 (ntimes: 1, nx: 2499, ny: 1499)
Dimensions without coordinates: ntimes, nx, ny
Data variables:
    goes_imager_projection  int32 ...
    x                       (nx) float32 ...
    y                       (ny) float32 ...
    time                    (ntimes) datetime64[ns] ...
    total_energy            (ntimes, nx, ny) float32 ...

Here is what is in xds.goes.imager_projection.attrs:

{'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'}

So, I would recommend building the CRS using pyproj.CRS.from_cf.

from pyproj import CRS

cc = CRS.from_cf(xds.goes_imager_projection.attrs)

Here is what cc looks like:

<Projected CRS: +proj=geos +h=35786023.0 +a=6378137.0 +b=6356752.3 ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Geostationary Satellite (Sweep X)
Datum: unknown
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

The next step is to reorganize the netCDF file into standard names/locations:

xds = xds.squeeze().rename_dims({"nx": "x", "ny": "y"}).transpose('y', 'x')
xds.coords["x"] = xds.x
xds.coords["y"] = xds.y
xds.coords["goes_imager_projection"] = xds.goes_imager_projection
xds.coords["time"] = xds.time

Here is what xds looks like now:

Dimensions:                 (x: 2499, y: 1499)
    goes_imager_projection  int32 ...
  * x                       (x) float32 -0.101304 -0.101248 ... 0.038584
  * y                       (y) float32 0.044296 0.044352 ... 0.128128 0.128184
    time                    int32 ...
Data variables:
    total_energy            (y, x) float32 ...

After that, write the CRS to the dataset using rioxarray's rio.write_crs:

xds.rio.write_crs(cc, inplace=True)
Dimensions:                 (x: 2499, y: 1499)
    goes_imager_projection  int32 ...
  * x                       (x) float32 -0.101304 -0.101248 ... 0.038584
  * y                       (y) float32 0.044296 0.044352 ... 0.128128 0.128184
    time                    int32 ...
    spatial_ref             int64 0
Data variables:
    total_energy            (y, x) float32 ...
    grid_mapping:  spatial_ref

According to this post http://meteothink.org/examples/meteoinfolab/satellite/geos-16.html, you just need to multiply by the perspective_point_height to convert to meters from radians.

sat_height = xds.goes_imager_projection.attrs["perspective_point_height"]
xds.x.values *= sat_height
xds.y.values *= sat_height

Then, you can reporoject the netCDF file using rioxarray's reproject functionality:

xds3857 = xds.rio.reproject("epsg:3857")

Here is what xds3857 looks like:

Dimensions:                 (x: 2495, y: 1506)
  * x                       (x) float64 -8.349e+06 -8.349e+06 ... -8.349e+06
  * y                       (y) float64 0.129 0.129 0.1289 ... 0.04467 0.04462
    time                    int32 -2147483647
    goes_imager_projection  int32 -2147483647
    spatial_ref             int64 0
Data variables:
    total_energy            (y, x) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    grid_mapping:  spatial_ref

And finally, you can write it to a geotiff with rioxarray using rio.to_raster.


You can install everything with conda:

conda install -c conda-forge rioxarray pyproj

And my ~/.condarc file looks like:

  - conda-forge
  - defaults
channel_priority: strict

Note: alternate approach using CRS definition: https://github.com/cf-convention/cf-conventions/issues/248#issuecomment-586350202

  • Hey @snowman2, after I followed these steps to reproject GOES data to wgs84 (epsg:4326), I lose the resolution and some pixels are lost. I there anyway to maintain the 2km resolution and not lose any pixels?
    – Clouseau
    Commented Oct 19, 2022 at 19:32
  • There is a resolution parameter you can pass in when reprojecting.
    – snowman2
    Commented Oct 19, 2022 at 20:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.