1

The data I read with xarray with rasterio engine is inverted along Y axis (This is SMAPL4 data). The correct dimension values are y=ds.y*(-1).

ds = xr.open_mfdataset(file_paths, engine="rasterio", chunks=chunks, combine='nested', concat_dim='time', parallel=True)

enter image description here

How can I fix this, without messing up the chunks?

Debug attempt 1

ds = ds.reindex(y=ds.y*(-1))

Result: The chunking over the y-axis disappears. Plus, numeric values in the variable data field somehow disappears.

Debug attempt 2:

ds= ds.isel(y=slice(None, None, -1))

Result: This code just re-selects the grids and does not change the dimension values themselves

Debug attempt 3:

ds= ds.reindex_like(ds_template)

Result: Again, numeric values within precipitation data disappear

Debug attempt 4: Change engine from rasterio to default (auto-detect).

ds = xr.open_mfdataset(file_paths, chunks=chunks, combine='nested', concat_dim='time', parallel=True)

Result: My file has some issues with the default engine; the data is read as "attributes."

2
  • In discussion with @raraki we've found a solution at least partly using xarray netcdf4-engine. Using rasterio-engine works too, but only for the original Earth Data dataset. Please follow up in github.com/pydata/xarray/issues/7621. It would be great if raraki could wrap this up and present as answer here. Mar 15 at 7:23
  • Hi @kmuehlbauer, thanks a lot for your help. I will post it here
    – raraki
    Mar 15 at 17:49

2 Answers 2

2

There is a better way to set the coordinates. For the SMAP dataset, you have to assign "cell_lat" and "cell_lon" as the coordinates of the dataset, and use them to plot. I am posting on behalf of kmuehlbauer based on our discussion in the Github issue. Thank you kmuehlbauer!

Import xarray as xr

Use xarray with netcdf engine

# load root group with coordinates
ds_NSIDC_root = xr.open_dataset(os.path.join(input_path, fn), group="/", engine='netcdf4')

# load data from Geophysical_Data group
ds_NSIDC_precip = xr.open_dataset(os.path.join(input_path, fn), group="Geophysical_Data", engine='netcdf4')

# merge groups
ds_NSIDC = xr.merge([ds_NSIDC_root, ds_NSIDC_precip])

# plot
ds_NSIDC.precipitation_total_surface_flux.plot()
# the above is essentially the same as 
# ds_NSIDC.precipitation_total_surface_flux.plot(x="x", y="y")

enter image description here

Use xarray with rasterio engine

ds_NASA_download_rasterio = xr.open_dataset(os.path.join(input_path, fn), engine='rasterio')
ds_NASA_download_rasterio = ds_NASA_download_rasterio.set_coords(["cell_lat", "cell_lon"])
ds_NASA_download_rasterio.Geophysical_Data_precipitation_total_surface_flux[0].plot(y="cell_lat", x="cell_lon")

enter image description here

Note that the rasterio method only works for files directly downloaded from NASA Earth Data, and not preprocessed data on the NSIDC server.

1

You can change the sign of the coordinates like this:

# Create an example dataset
ds = xr.Dataset({"data": ("x", np.arange(11))}, 
                coords = {"x": np.arange(-5,6,1)})

print(f"[x =-4]: {ds['data'].sel(x=-4).values}")

# Make the negative coordinates positive and vice versa.
ds = ds.assign_coords({"x": ds.x * -1})

print(f"[x =-4]: {ds['data'].sel(x=-4).values}")

Gives:

>>> [x = -4]: 1
>>> [x = -4]: 9
3
  • I am coming to think this is a bug in the dataset, and cannot be fixed. Even if you sort the data according to Y coordinate, the original coordinates dont change.
    – raraki
    Mar 8 at 19:29
  • Yeah you're right, think my updated answer will fix your problem. Mar 9 at 9:17
  • Agreed, that's the only way left and I incorporated it to my code. Thank you, Bert!
    – raraki
    Mar 10 at 20:04

Your Answer

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

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