reproject_match() returns: ValueError: IndexVariable objects must be 1-dimensional

[dkczk]

Hey folks,
I'm working with unequal spaced data and want to regularise using rioxarrays reproject_match(). For that I extract lan/lot, the desired variable and create a new DataArray with that. Afterwards I create a dummy DataArray with the target grid. In the last step I'm using reproject_match() to project the irregular grid to the target grid. But when executing reproject_match() I get: "ValueError: IndexVariable objects must be 1-dimensional".
I did:

import xarray as xr
import numpy as np
import rioxarray as rxr
from rasterio.enums import Resampling


file_path = "/path/to/irregular_gridded_data.nc"
ds = xr.open_dataset(file_path, group="PRODUCT")

ds_var = ds["carbonmonoxide_total_column"]
ds_lat = ds.coords["latitude"]
ds_lon = ds.coords["longitude"]

ds_var_val = ds_var.values
ds_lat_val = ds_lat.values
ds_lon_val = ds_lon.values

res = 0.25

da = xr.DataArray(
    data=ds_var_val[0],
    coords={"latitude": (("y", "x"), ds_lat_val[0]),
            "longitude": (("y", "x"), ds_lon_val[0])
            },
    dims=["y", "x"]
)

da.attrs["_FillValue"] = np.nan
da = da.rio.write_crs("EPSG:4326")


new_lon = np.arange(ds_lon.min() + res / 2, ds_lon.max() + res / 2, res)
new_lat = np.arange(ds_lat.min() + res / 2, ds_lat.max() + res / 2, res)

lon2d, lat2d = np.meshgrid(new_lon, new_lat)

dummy_vals = np.ones_like(lat2d)
dummy_vals[dummy_vals == 1] = np.nan
dummy = xr.DataArray(
    data=np.ones_like(lat2d),
    coords={"latitude": (("y", "x"), lat2d),
            "longitude": (("y", "x"), lon2d)
            },
    dims=["y", "x"]
)
dummy = dummy.rio.write_crs("EPSG:4326")

da_resampled = da.rio.reproject_match(dummy, Resampling.average)

I would like to understand why the error happens and also (obviously) how to fix it.

[snowman2]

Related: #209

[dkczk]

Thank you very much for your reply. Your suggested solution only works at first glance. It produces an array with a lot of artefacts. I guess its a result of the interpolation method. From the look of it, it's linear interpolated which doesn't work well if you have a lot of NaN's in your data. In addition to that it is very slow. It takes nearly sixty times the time compared to R's terra approach.

[snowman2]

It does look like it is trying to fill in data.

How does the default rasterization work?

from geocube.api.core import make_geocube

geo_grid = make_geocube(
    vector_data=gdf,
    resolution=(-0.25, 0.25),
)

[dkczk]

That works, thank you! But still the solution takes a lot of time (1 m 4 s vs. 1.2 s) and is therefore unsuitable for bulk processing. I adapted your solution as follows:

import xarray as xr
import numpy as np
from geocube.api.core import make_geocube
from geocube.rasterize import rasterize_points_griddata
import geopandas as gpd


file_path = "/path/to/irregular_gridded_data.nc"
ds = xr.open_dataset(file_path, group="PRODUCT")
var = "carbonmonoxide_total_column"
res = 0.25

ds_var = ds[var]
ds_lat = ds.coords["latitude"]
ds_lon = ds.coords["longitude"]

ds_var_val = ds_var.values
ds_lat_val = ds_lat.values
ds_lon_val = ds_lon.values

da = xr.DataArray(
    data=ds_var_val[0],
    coords={"latitude": (("lat", "lon"), ds_lat_val[0]),
            "longitude": (("lat", "lon"), ds_lon_val[0])
            },
    dims=["lat", "lon"]
)

da.attrs["_FillValue"] = np.nan
da.name = var

df = da.squeeze().to_dataframe().reset_index()
geometry = gpd.points_from_xy(df.longitude, df.latitude)
gdf = gpd.GeoDataFrame(df, crs="EPSG:4326", geometry=geometry)
geo_grid = make_geocube(
    vector_data=gdf,
    resolution=(-res, res)
)

I don't know if there is something to improve the processing speed using geocube.