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)

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; — ba0a2794, Feb 07 '20 at 16:03
Looks like the first step is to convert it to geodetic lat/lon. https://makersportal.com/blog/2018/11/25/goes-r-satellite-latitude-and-longitude-grid-projection-algorithm — snowman2, Feb 08 '20 at 02:32
This also looks promising: http://meteothink.org/examples/meteoinfolab/satellite/geos-16.html — snowman2, Feb 08 '20 at 04:36

snowman2 · Accepted Answer · 2020-07-15T15:11:24.000

8

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.to_string(), 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:

edited Jul 15 '20 at 15:11

answered Feb 08 '20 at 04:56

snowman2

7,321
12
29
54

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:
– ba0a2794 Feb 10 '20 at 15:58
1

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 Feb 11 '20 at 14:10
1

What method did you use for installation? – snowman2 Feb 11 '20 at 14:11
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! – ba0a2794 Feb 11 '20 at 15:33
@snowman2 , what is CRS ?‌ How should I import it? – Muser Jul 15 '20 at 14:40
from pyproj import CRS – snowman2 Jul 15 '20 at 15:11

Using Rioxarray / QGIS Projection

1 Answers1

Linked