A Netcdf file were generate from Google Earth Engine with the aid of packages Xee: Xarray + GEE and I can't open correctly it in QGIS

Question

I am trying to generated Netcdf files from a imageCollection. The code is working fine. I can save as nc file and reopened it in Colab using Xarray package. I also can open it in Panoply, however it is not open correct in QGis (QGIS 3.22.4-Białowieża). What am I missing? Thanks

The code in colab is here

The code:

# capture
if 'google.colab' in str(get_ipython()):
    !pip install xee rioxarray netCDF4

import ee
import geemap
import xee
import xarray as xr

ee.Authenticate()

ee.Initialize(project='ee-my-xavierac', opt_url='https://earthengine-highvolume.googleapis.com')

"""ROI"""

# Create a geodesic polygon.
# https://developers.google.com/earth-engine/guides/geometry_visualization_info
roi = ee.Geometry.Polygon([
  [[-46.757813, -23.337212],
      [-46.757813, -1.158979],
      [-33.398438, -1.158979],
      [-33.398438, -23.337212],
      [-46.757813, -23.337212]],
]);

"""Data from: https://gee-community-catalog.org/projects/br_dwgd/"""

name_var = 'PR'
# Define a scaling factor
offset = 225
scale  = 0.006866665

# Function to apply the scaling factor to a specific band
def scaleBand(image):
    scaledImage = image.select("b1").rename([name_var]) \
                       .multiply(scale).add(offset) \
                       .set({'date': image.date().format('yyyy-MM-dd')}) #\
                       # .clip(roi)
    return scaledImage.copyProperties(image, image.propertyNames())

# cmap
color_pal = ['000066', '001199', '0044BB', '0077DD', '33AAEE', '66CCFF', 'FFDDCC', 'FFBB99', 'FF9966', 'FF6644']

collection = ee.ImageCollection(f"projects/sat-io/open-datasets/BR-DWGD/{name_var}") \
                         .filterDate('2019-01-01', '2019-01-15') \
                         .map(scaleBand)

"""from: https://www.geopythontutorials.com/notebooks/xee_ic_to_netcdf.html"""

ds = xr.open_dataset(
    collection,
    engine='ee',
    scale=.1,
    projection=collection.first().select(0).projection()
)

ds

# https://www.geopythontutorials.com/notebooks/xee_ic_to_netcdf.html
clipped_ds = ds \
  .rename({'lon': 'x', 'lat':'y'}) \
  .rio.clip([roi.getInfo()])
clipped_ds

export_ds = clipped_ds.transpose()
export_ds[name_var].isel(time=0).plot()

"""Export"""

# Enable compression
encoding = {name_var: {'zlib': True}}
export_ds.to_netcdf('test.nc', encoding=encoding)

"""Open file to check"""

data_test = xr.open_mfdataset('test.nc')
data_test

data_test[name_var].isel(time=0).plot()

Answer 1

This code is able to save NetCDF file that was open correct in QGis.

# capture
if 'google.colab' in str(get_ipython()):
    !pip install xee rioxarray netCDF4

import ee
import geemap
import xee
import xarray as xr

ee.Authenticate()

ee.Initialize(project='ee-my-xavierac', opt_url='https://earthengine-highvolume.googleapis.com')

"""ROI"""

# Create a geodesic polygon.
# https://developers.google.com/earth-engine/guides/geometry_visualization_info
roi = ee.Geometry.Polygon([
  [[-46.757813, -23.337212],
      [-46.757813, -1.158979],
      [-33.398438, -1.158979],
      [-33.398438, -23.337212],
      [-46.757813, -23.337212]],
]);

"""Data from: https://gee-community-catalog.org/projects/br_dwgd/"""

name_var = 'PR'
# Define a scaling factor
offset = 225
scale  = 0.006866665

# Function to apply the scaling factor to a specific band
def scaleBand(image):
    scaledImage = image.select("b1").rename([name_var]) \
                       .multiply(scale).add(offset) \
                       .set({'date': image.date().format('yyyy-MM-dd')}) #\
                       # .clip(roi)
    return scaledImage.copyProperties(image, image.propertyNames())

# cmap
color_pal = ['000066', '001199', '0044BB', '0077DD', '33AAEE', '66CCFF', 'FFDDCC', 'FFBB99', 'FF9966', 'FF6644']

collection = ee.ImageCollection(f"projects/sat-io/open-datasets/BR-DWGD/{name_var}") \
                         .filterDate('2019-01-01', '2019-01-15') \
                         .map(scaleBand)

"""from: https://www.geopythontutorials.com/notebooks/xee_ic_to_netcdf.html"""

ds = xr.open_dataset(
    collection,
    engine='ee',
    scale=.1,
    projection=collection.first().select(0).projection()
)
ds

ds[name_var].isel(time=0).plot()

# https://www.geopythontutorials.com/notebooks/xee_ic_to_netcdf.html
clipped_ds = ds \
  .rename({'lon': 'x', 'lat':'y'}) \
  .rio.clip([roi.getInfo()])
clipped_ds

clipped_ds[name_var].isel(time=0).plot()

export_ds = clipped_ds.transpose('time', 'y', 'x')
export_ds

export_ds[name_var].isel(time=0).plot()

"""Export"""

# Enable compression
encoding = {name_var: {'zlib': True}}
export_ds.to_netcdf('test.nc', encoding=encoding)

"""Open file to check"""

data_test = xr.open_mfdataset('test.nc')
data_test

data_test['PR'].isel(time=0).plot()

Answer 2

Try:

export_ds.rio.set_spatial_dims(x_dim='x', y_dim='y', inplace=True)
export_ds.rio.write_grid_mapping(inplace=True)
export_ds.rio.write_crs(4326, inplace=True)
export_ds.to_netcdf('test.nc')

Normally, QGIS will have a better understanding of the structure of your netcdf, and the question mark displayed in QGIS next to your raster layer will disappear.