6

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

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

https://www.dropbox.com/s/f53fkqqyamd1rb9/glm.nc?dl=0

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

Driver: netCDF/Network Common Data Format
Files: glm.nc
Size is 1499, 2499
Metadata:
goes_imager_projection#grid_mapping_name=geostationary
goes_imager_projection#inverse_flattening=298.2572221
goes_imager_projection#latitude_of_projection_origin=0
goes_imager_projection#longitude_of_projection_origin=-75
goes_imager_projection#long_name=GOES-R ABI fixed grid projection
goes_imager_projection#perspective_point_height=35786023
goes_imager_projection#semi_major_axis=6378137
goes_imager_projection#semi_minor_axis=6356752.31414
goes_imager_projection#sweep_angle_axis=x
NETCDF_DIM_EXTRA={ntimes}
NETCDF_DIM_ntimes_DEF={1,5}
NETCDF_DIM_ntimes_VALUES=36920
total_energy#grid_mapping=goes_imager_projection
total_energy#long_name=Total radiant energy
total_energy#missing_value=-9999
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
Metadata:
grid_mapping=goes_imager_projection
long_name=Total radiant energy
missing_value=-9999
NETCDF_DIM_ntimes=36920
NETCDF_VARNAME=total_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 is 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?

2
  • 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
    Dec 20, 2019 at 18:06
  • for gdalwarp use -s_srs, similarly in rasterio
    – mdsumner
    Apr 14 at 1:23

1 Answer 1

10

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:

<xarray.Dataset>
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:

<xarray.Dataset>
Dimensions:                 (x: 2499, y: 1499)
Coordinates:
    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.to_string(), inplace=True)
<xarray.Dataset>
Dimensions:                 (x: 2499, y: 1499)
Coordinates:
    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 ...
Attributes:
    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:

<xarray.Dataset>
Dimensions:                 (x: 2495, y: 1506)
Coordinates:
  * 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
Attributes:
    grid_mapping:  spatial_ref

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

xds3857.total_energy.rio.to_raster("total_energy.tif")

You can install everything with conda:

conda install -c conda-forge rioxarray pyproj

And my ~/.condarc file looks like:

channels:
  - conda-forge
  - defaults
channel_priority: strict

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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