Xarray returns duplicate Longitude/Latitude values with sel method

[rogomichael]

What is your issue?

I am trying to select air temperature from Netcdf file and polgon regions from NUTS3 shape file with the aim of masking both to extract the values. However, for some regions, xarray returns a data frame with the same latitude or longitude values which cannot be masked by geopandas. The shape file and the netcdf I am using are zipped here

Here is the code:

import os
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import seaborn as sn
import geopandas as gpd
import earthpy as et
import xarray as xr
import regionmask
from sultan.api import Sultan
#Plotting options
sn.set(font_scale=1.3)
sn.set_style("white")

nc_file = xr.open_mfdataset('year2019.nc') #, chunks = {'time': 10})
shapefile = "NUTS_RG_60M_2021_4326/NUTS_RG_60M_2021_4326.shp"
countries = gpd.read_file(shapefile)
querry = 3
countries = countries[countries['LEVL_CODE'] == querry]

```#Select the longitude and latitude coordinates of the region of interest in the shape file.```
#These values will be passed to xarray
german_state = countries[countries.NUTS_ID.isin(["DE255"])]
german_state.plot()
boundaries = german_state.total_bounds
print(boundaries)

ger_mask = regionmask.mask_3D_geopandas(german_state, nc_file.longitude, nc_file.latitude, method="pygeos")
print(ger_mask)
state_lat = [float(boundaries[1]), float(boundaries[3])]
print("Boundary Latitude", state_lat)
state_lon = [float(boundaries[0]), float(boundaries[2])]
print(" Boundary Longitude", state_lon)
print(state_lat, state_lon)

```#Slice the data by time and spatial extent by passing the values selected from the shape file```
two_months_255 = nc_file["tg"].sel(
    longitude=[state_lon[0], state_lon[1]],
    latitude=[state_lat[0], state_lat[1]], method="nearest")

two_months_255.where(ger_mask)

Error:
File /opt/homebrew/lib/python3.9/site-packages/xarray/core/indexes.py:449, in PandasIndex.reindex_like(self, other, method, tolerance)
    445 def reindex_like(
    446     self, other: PandasIndex, method=None, tolerance=None
    447 ) -> dict[Hashable, Any]:
    448     if not self.index.is_unique:
--> 449         raise ValueError(
    450             f"cannot reindex or align along dimension {self.dim!r} because the "
    451             "(pandas) index has duplicate values"
    452         )
    454     return {self.dim: get_indexer_nd(self.index, other.index, method, tolerance)}

ValueError: cannot reindex or align along dimension 'latitude' because the (pandas) index has duplicate values

ask.zip

[TomNicholas]

Hi @rogomichael, thanks for raising your issue.

Could you possibly turn your problem into a minimal reproducible example? That would make it much easier for us to help you.

You've included a lot of code, but from looking at it it seems to me that the only lines which matter are

nc_file = xr.open_mfdataset('year2019.nc')

nc_file["tg"].sel(
    longitude=[state_lon[0], state_lon[1]],
    latitude=[state_lat[0], state_lat[1]],
    method="nearest"
)

This suggests to me that your original file has duplicate values, as all xarray has done here is open it.

(However I still am not going to be able to reproduce this error without access to the file year2019.nc.)

[rogomichael]

@TomNicholas ,
I have reduced the code to a smaller portion. The reason I include a long piece of code is because the longitude and latitude values are read from the shape file and passed to xarray to read from year2019.nc. I previously attached the file in the issue above but you can still access it here.

https://1drv.ms/u/s!AvYOTGjStibEmGGsJwRlpoBfLvHm?e=LBLm95

I hope it is accessible.

import xarray as xr

nc_file = xr.open_mfdataset('year2019.nc') #, chunks = {'time': 10})
two_months_255 = nc_file["tg"].sel(
    longitude = [49.294, 49.376],
    latitude = [10.971, 11.087])

[benbovy]

Hi @rogomichael,

It is expected in your example that sel with your given values and method='nearest' returns coordinates with duplicate longitude & latitude values.

The resolution of the grid is too coarse compared to your ROI, which is why the two provided values for latitude (longitude) both select the same grid point.

If you want to use sel for extracting a range over the latitude and longitude dimensions, it is better to use slices. It better illustrates that your ROI is smaller than the grid resolution:

ds.sel(
    longitude=slice(49.294, 49.376),
    latitude=slice(10.971, 11.087)
)
# <xarray.Dataset>
# Dimensions:    (time: 365, longitude: 1, latitude: 0)
# Coordinates:
#   * time       (time) datetime64[ns] 2019-01-01 2019-01-02 ... 2019-12-31
#   * longitude  (longitude) float64 49.38
#   * latitude   (latitude) float64 
# Data variables:
#     tg         (time, latitude, longitude) float32 ...

[rogomichael]

Hi @benbovy ,
As you can see, the slice option returns an empty latitude. From your explanation above, does it mean that it is not possible to select the latitude values and in some cases longitude values? I needed both the longitude and latitude values.

[benbovy]

You'd need to adjust the latitude / longitude values passed to sel so that your ROI contains (or is nearest to) more than one gird point in latitude or longitude. I'm afraid there's not automatic way in Xarray to do that.