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月03日T19: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)

• snowman2: Can you post a link to the input file?
• ba0a2794: Here is the aws command: aws s3 cp s3://noaa-goes17/ABI-L2-LSTF/2020/034/19/OR_ABI-L2-LSTF-M6_G17_s20200341900321_e20200341909388_c20200341910038.nc ./ --no-sign-request;
• snowman2: Looks like the first step is to convert it to geodetic lat/lon. makersportal.com/blog/2018/11/25/…
• snowman2: This also looks promising: meteothink.org/examples/meteoinfolab/satellite/geos-16.html

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:

enter image description here

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

• ba0a2794: I am getting this error at the reprojection line: CRSError: The WKT could not be parsed. OGR Error code 5. I have GDAL 2.2.3, rasterio 1.1.2, and rioxarray 0.0.21. Also, not sure if it is relevant but I have to import rioxarray after opening the dataset, otherwise I get OSError: [Errno -101] NetCDF: HDF error:
• snowman2: Looks like you are using a version of GDAL that does not support WKT2. I updated the script so that it should work now.
• snowman2: What method did you use for installation?
• ba0a2794: rioxarray, xarray, and rasterio were through pip. Not sure, but I think GDAL was through apt. Thanks for all the help and for making rioxarray!
• Muser: @snowman2 , what is CRS ?‌ How should I import it?

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")