Using Rioxarray / QGIS Projection

Question

For the past couple weeks I have been trying to view GOES17 data from netCDF files (converted to geotif) in QGIS 3.10 but cannot get the projection to work correctly.

I have attempted numerous methods but most recently tried the procedure of the top answer in this post: Converting NetCDF dataset array to GeoTiff using rasterio Python

When loading the .tif into QGIS, it appears in the wrong place relative the OpenStreetMap. I have tried several projections including EPSG:3857, which is what appears in the bottom right in QGIS when the OSM is loaded. They are all wrong.

I have also tried this answer: How do I add projection to this NetCDF file? (Satellite)

When attempting the reproject function I get an error.

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

Error:

DimensionError: x dimension not found. 'set_spatial_dims()' can address this.

xds:

<xarray.Dataset>
Dimensions:                                 (number_of_LZA_bounds: 2, number_of_SZA_bounds: 2, number_of_image_bounds: 2, number_of_time_bounds: 2, x: 1086, y: 1086)
Coordinates:
    t                                       datetime64[ns] 2020-02-03T19:05:05.476645888
  * y                                       (y) float32 0.1519 ... -0.15190002
  * x                                       (x) float32 -0.1519 ... 0.15190002
    goes_imager_projection                  int32 -2147483647
    y_image                                 float32 0.0
    x_image                                 float32 0.0
    retrieval_local_zenith_angle            float32 85.0
    quantitative_local_zenith_angle         float32 70.0
    solar_zenith_angle                      float32 180.0
    time                                    int32 -2147483647
    spatial_ref                             int64 0

The issue continues to persist after doing the suggestion.

xds.rio.set_spatial_dims("x","y",inplace=True)

Answer 1

Upon inspecting the dataset, I realized that the units of the data are in radians.

import xarray
import rioxarray
from pyproj import CRS

xds = xarray.open_dataset("OR_ABI-L2-LSTF-M6_G17_s20200341900321_e20200341909388_c20200341910038.nc")

Inside the x variable the attributes say the data is in radians:

xds.x.attrs

{'units': 'rad',
 'axis': 'X',
 'long_name': 'GOES Projection x-Coordinate',
 'standard_name': 'projection_x_coordinate'}

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.

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

Next, set the CRS of the dataset:

cc = CRS.from_cf(xds.goes_imager_projection.attrs)
xds.rio.write_crs(cc, inplace=True)

Also, there are only two variables with the x and y dimensions:

Data variables:
    LST                                     (y, x) float32 ...
    DQF                                     (y, x) float32 ...

As such, only those ones will work, so you need to pull those out:

xds = xds[["LST", "DQF"]]

Then, you can reproject and export to raster:

xds3857 = xds.rio.reproject("EPSG:3857")
xds3857.rio.to_raster("geos17.tif")

It seemed to work:

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

Answer 2

The magic sauce is indeed the conversion of the coordinates to meters.

Xarray syntax changed slightly since the answer above, so here is a currently working example:

import s3fs 
import xarray as xr
import rioxarray
from pyproj import CRS

fs = s3fs.S3FileSystem(anon=True)

path = "s3://noaa-goes16/ABI-L1b-RadF/2017/129/11/OR_ABI-L1b-RadF-M3C03_G16_s20171291115392_e20171291126159_c20171291126196.nc"

with fs.open(path) as fileObj:
    with xr.open_dataset(fileObj, engine='h5netcdf') as ds:
        sat_height = ds["goes_imager_projection"].attrs["perspective_point_height"]
        ds = ds.assign_coords({
            "x": ds["x"].values * sat_height,
            "y": ds["y"].values * sat_height,
        })
    
        crs = CRS.from_cf(ds["goes_imager_projection"].attrs)
        ds.rio.write_crs(crs.to_string(), inplace=True)
    
        da = ds["Rad"]

        da_web = da.rio.reproject("epsg:3857")