@kmuehlbauer

The code I provided above (with the fn_NASA_Earthdata_download file) works for you too, but with fn_NSIDC_output it does not?

You are right. I've gone made updates to both the issue description (i.e., the first post) and the Jupyter notebook based on our discussion. So, there might be an issue with the data or the rasterio engine, but for now, we'll use the NetCDF engine for NSIDC data, and hope that someone figures out the problem. Also, thanks for your suggestion on SO; it's much appreciated!

@RY4GIT Glad you can use the code. And good there's a way to use the data as is without tampering.

Regarding the different rasterio-approaches, just to be on the same page. The code I provided above (with the fn_NASA_Earthdata_download file) works for you too, but with fn_NSIDC_output it does not? Out of the box, the second also doesn't work for me. Didn't check this.

I've no immediate idea, what's going on in that case. The only thing I can think of is, that in the conversion process something has gone wrong with the georeferencing and/or the image origin. I hope you can figure that out, eventually.

Here for reference my package versions:

The post on SO is actually what I wrote. I didn't feel right about messing around with y coordinates, so I came here to figure out what was happening.

Thank you so much for taking the time out of your busy schedule to download the data and test it.

The Xarray Netcdf engine solution worked on my end as well. It was really helpful in exploring what works best for this specific data, where latitude and longitude are stored in different locations. I appreciate your code and will use it for my future processing.

As for the rasterio solution, the output is still flipped along the Y axis. I updated both the rasterio and xarray packages to the latest versions, but the results were still the same. I used the same lines of code that you provided, but with a different file (instead of direct NASA output, I used NSIDC output, as the filename was specified as a variable, fn_NSIDC_output).

@RY4GIT I forgot to add my package versions. I'll add them the next day. But I'm pretty sure to have the most recent one's installed in a Python 3.10 conda-forge environment.

For the rasterio-approach:

python ds_NASA_download_rasterio = xr.open_dataset(os.path.join(input_path, fn_NASA_Earthdata_download), 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")

First, I've found this on SO: https://gis.stackexchange.com/questions/454543/fixing-the-flipped-inverted-y-axis-in-the-xarray-with-rasterio

For your data reading method 1, it works for me like this:

```python

load root group with coordinates

ds_NSIDC_root = xr.open_dataset(os.path.join(input_path, fn_NSIDC_output), group="/", engine='netcdf4')

load data from Geophysical_Data group

ds_NSIDC_precip = xr.open_dataset(os.path.join(input_path, fn_NSIDC_output), 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")

```

This can be expanded to use the cell_lon/cell_lat:

```

set lat/lon as coords

ds_NSIDC = ds_NSIDC.set_coords(["cell_lon", "cell_lat"]) ds_NSIDC.precipitation_total_surface_flux.plot(x="cell_lon", y="cell_lat") ```

Update: This should also work like above when directly using NASA Earth Data.

Thank you for checking the code, @kmuehlbauer. However, I am still experiencing issues even after fixing it.

Regarding the coordinates: - The y coordinate starts from the maximum values and decreases as you move down the array. - The latitude values increase with increasing y-coordinate values. - i.e., The latitude corresponding to the maximum y coordinate is around 85, and the one corresponding to the minimum y coordinate is around -85.

Regarding the precipitation data: - The structure is the same as the cell_lat/cell_lon array. - The y coordinate starts from the maximum values and decreases as you move down the array. - If I plot the precipitation with the xy coordinate, it looks like this, indicating that the data is still inverted. Please refer to these images for further clarification: Note that the unit is read incorrectly as degrees_north and degrees_east; please ignore it.

Here is the updated code:

``` lons = ds_NSIDC_output_rasterio.cell_lon.load() lons_array = lons[0][0]

lats = ds_NSIDC_output_rasterio.cell_lat.load() lats_array = np.arange(lats[0][0][0], lats[0][-1][0], -1*(lats.max().values-lats.min().values)/ds_NSIDC_output.y.size)

ds_NSIDC_output_rasterio_xymanual = xr.DataArray( data=ds_NSIDC_output_rasterio.precipitation_total_surface_flux.sel(band=1).values, dims = ["y","x"], coords = dict( y = lats_array, x = lons_array.values ) ) ds_NSIDC_output_rasterio_xymanual.plot() ```

@RY4GIT You would need to assign cell_lat/cell_lon as coordinates of the dataset and then use x="cell_lon" and y="cell_lat" in the call to .plot().

Plotting without that, will use the dimension coordinates (x,y) and there might be a difference in origin (upper vs. lower).