NetCDF Projection – How to Add Projection to NetCDF File (Satellite)?

gdalgoes-16rioxarraysatellite

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?

Best Answer

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

Related Solutions

[GIS] Geotransformation for polar stereographic

Too long for a comment, this is to accompany @Matej's answer.

Add the “.grd” data into ArcGIS.
Use the “Raster to Other Format” function and convert your .grd file into an ESRII GRID format. This is important because most of the raster functions in ArcGIS are accessible for this format only, either that or it usually gets too slow when you use it on other formats.
As it already has the projection file associated to the file. Rather than projecting the new converted data, define its projection. ArcToolbox > Data Management Tools > Projections and Transformations > Define Projection. You can navigate to the pre-defined polar stereo graphic ESRII projection and see if its parameters match the one given in the metadata (it doesn’t), so you can modify it as per @Matej. Only here - rather than modifying, create a new one based on the NPS projection with the central meridian and Latitude of origin changed and save it on disk, then navigate to the new projection and use that when defining your projection. This is because your on the fly modification won't be available later on when you want to use it to set the coordinate system for your data frame, which is necessary in order to align your other boundary files to the same modified NPS system.

GDAL – How to Reproject a NetCDF File?

gdalinfo on your data returned below:

Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Driver: netCDF/Network Common Data Format
Files: C:\Users\userr\Downloads\OR_ABI-L1b-RadC-M3C01_G16_s20190321802143_e20190321804516_c20190321804558.nc
Size is 5000, 3000
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unknown",
        DATUM["unknown",
            SPHEROID["Spheroid",6378137,298.2572221]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Geostationary_Satellite"],
    PARAMETER["central_meridian",-75],
    PARAMETER["satellite_height",35786023],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    EXTENSION["PROJ4","+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs  +sweep=x"]]
Origin = (-3627271.328111211769283,1583173.512170388596132)
Pixel Size = (1002.008657743770527,1002.008657743770527)

Geostationary_Satellite... I honestly don't know anything about this projection... So I tried gdalwarp giving proj4 string explicitly to your nc (Rad layer) as below:

gdalwarp -of netCDF 
-s_srs "+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0
+ellps=GRS80 +units=m +no_defs +sweep=x" 
-t_srs EPSG:3857  -tr 1002.008657743770527 1002.008657743770527 
NETCDF:"OR_ABI-L1b-RadC-M3C01_G16_s20190321802143_e20190321804516_c20190321804558.nc":Rad 
"OR_ABI-L1b-RadC-M3C01_G16_s20190321802143_e20190321804516_c20190321804558_3857.nc"

[Below] Both images are underlain by Bing image (in QGIS...Left image is displayed on On-the-fly projection).

[Update]

Following this post we started working on a new scope to create True color image from three GOES-16 images (Band1, 2, and 3).
Initial idea was to reproject the provided NetCDF files along with all variables (Rad for radiance and DQF as QC value). Then with three reprojected bands David's software should be able to convert them into color enhanced Jpeg image... Using GDAL tool I could reproject one variable (Rad) but not two variables together.
Instead, I extracted Rad variable from each band and created RGB composite by SAGA. There was a nice article Bah et al. (2018) how to create "Green" band by mixing red, blue and NIR bands. Following this paper I mixed 45% of Band2 (red), 45% of Band1 (blue) and 10% of Band3 (Nir) by raster calculator. Output composite image is as below:

Best Answer

Related Solutions

[GIS] Geotransformation for polar stereographic

GDAL – How to Reproject a NetCDF File?

Related Question