metno / pyaerocom

Python tools for climate and air quality model evaluation
https://pyaerocom.readthedocs.io/
GNU General Public License v3.0
25 stars 13 forks source link

Colocation after running GriddedData.interpolate() does not work #1361

Closed thorbjoernl closed 1 week ago

thorbjoernl commented 2 weeks ago

Describe the bug Please provide a clear and concise description of what the bug is.

Traceback (most recent call last):
  File "/home/thlun8736/Documents/work/aerotools-scripts/bilinear.py", line 36, in <module>
    colocated_data = pya.colocation.colocation_utils.colocate_gridded_ungridded(
  File "/home/thlun8736/Documents/work/aerotools-scripts/venv/lib/python3.10/site-packages/pyaerocom/colocation/colocation_utils.py", line 818, in colocate_gridded_ungridded
    grid_stat_data = data.to_time_series(longitude=ungridded_lons, latitude=ungridded_lats)
  File "/home/thlun8736/Documents/work/aerotools-scripts/venv/lib/python3.10/site-packages/pyaerocom/griddeddata.py", line 1174, in to_time_series
    result = self._to_time_series_xarray(scheme=scheme, add_meta=add_meta, **coords)
  File "/home/thlun8736/Documents/work/aerotools-scripts/venv/lib/python3.10/site-packages/pyaerocom/griddeddata.py", line 1217, in _to_time_series_xarray
    subset = extract_latlon_dataarray(
  File "/home/thlun8736/Documents/work/aerotools-scripts/venv/lib/python3.10/site-packages/pyaerocom/helpers.py", line 189, in extract_latlon_dataarray
    subset = arr.sel(where, method=method)
  File "/home/thlun8736/Documents/work/aerotools-scripts/venv/lib/python3.10/site-packages/xarray/core/dataarray.py", line 1670, in sel
    ds = self._to_temp_dataset().sel(
  File "/home/thlun8736/Documents/work/aerotools-scripts/venv/lib/python3.10/site-packages/xarray/core/dataset.py", line 3184, in sel
    query_results = map_index_queries(
  File "/home/thlun8736/Documents/work/aerotools-scripts/venv/lib/python3.10/site-packages/xarray/core/indexing.py", line 185, in map_index_queries
    grouped_indexers = group_indexers_by_index(obj, indexers, options)
  File "/home/thlun8736/Documents/work/aerotools-scripts/venv/lib/python3.10/site-packages/xarray/core/indexing.py", line 144, in group_indexers_by_index
    raise KeyError(f"no index found for coordinate {key!r}")
KeyError: "no index found for coordinate 'lat'"

To Reproduce

if name == "main": mod_reader = pya.io.mscw_ctm.reader.ReadMscwCtm( "EMEP", data_dir="/lustre/storeB/project/fou/kl/emep/ModelRuns/2024_REPORTING/TRENDS/2020", )

mod_data: pya.griddeddata.GriddedData = mod_reader.read_var("concnh3", "daily")
obs_reader = pya.io.ReadUngridded("EBASMC")

obs_data: pya.ungriddeddata.UngriddedData = obs_reader.read(
    vars_to_retrieve="concnh3"
)

coord_pairs = obs_data._get_stat_coords()[1]
lat, lon = list(zip(*coord_pairs))

#mod_data = mod_data.interpolate(scheme=iris.analysis.Linear(), latitude=lat, longitude=lon)

colocated_data = pya.colocation.colocation_utils.colocate_gridded_ungridded(
    mod_data, obs_data
)

print(colocated_data)


**Expected behavior**
Both cases should work without error.

**Screenshots

![Pasted image 20241007092439](https://github.com/user-attachments/assets/d2f1b872-d8f0-440c-9710-ddb04b58058c)

![Pasted image 20241007092619](https://github.com/user-attachments/assets/f36448b9-49b6-4aa2-bcc2-57f6fdff70f9)

![Pasted image 20241007092714](https://github.com/user-attachments/assets/a4d6bedf-8e6e-4304-9ff2-51f11f73b8fa)

**Additional context**
Add any other context about the problem here.
thorbjoernl commented 1 week ago

The issue was on my end. lat and lon passed to sample_coords needs to be a ascending ordered list with duplicates removed. Otherwise interpolation will turn the dimcoords int auxcoords during interpolation (since dimcoords must be monotonic)

I'll make a small PR to clarify this requirement to GriddedData.interpolate()