DataCube.md

I am working with pollution data from ECMWF using the Python modules xarray and xvec.

ds = xarray.open_dataset("/content/drive/MyDrive/UNaM/FCF/TUSIGyT/Geomática_Analisis_Territorial/2024/Datos/2022_pm2p5.nc")
ds

Plotting it:

ds.pm2p5.isel(time=slice(0,3)).plot(x="longitude", y="latitude", col="time", col_wrap=3)

So, using a polygon GeoDataFrame, I ran the zonal_stats analysis:

dep = gpd.read_file("/content/drive/MyDrive/UNaM/FCF/TUSIGyT/Geomática_Analisis_Territorial/2024/Datos/departamento_continente.shp")
dep.plot(facecolor="None")

aggregated = ds.xvec.zonal_stats(
    dep.geometry.values,
    x_coords="longitude",
    y_coords="latitude
)
aggregated

What I am trying to do is plot the zonal_stats results, believing that it is possible.

I expect a plot like a coropleth, like the following, but not needing to convert back to geodataframe for eac visual inspection on time index, like the following:

dep_ts = aggregated.xvec.to_geodataframe()
dep_ts

dep_ts.loc['2022-01-01 06:00:00	'].plot(column='pm2p5')

Would exist another way of doing it? How could I represent an xarray DataSet's geometry in a plot?