6

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

1 Answer 1

9

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

6
  • 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
    Commented Feb 10, 2020 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
    Commented Feb 11, 2020 at 14:10
  • 1
    What method did you use for installation?
    – snowman2
    Commented Feb 11, 2020 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
    Commented Feb 11, 2020 at 15:33
  • @snowman2 , what is CRS ?‌ How should I import it?
    – Muser
    Commented Jul 15, 2020 at 14:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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