This feature has not been implemented in pysteps. Implementing an utility function for doing this could be useful.

In the meanwhile, you can use the following example:

import numpy as np
import pyproj

xr = np.arange(geodata["x1"], geodata["x2"], geodata["xpixelsize"])
xr += 0.5 * (xr[1] - xr[0])
yr = np.arange(geodata["y1"], geodata["y2"], geodata["ypixelsize"])
yr += 0.5 * (yr[1] - yr[0])
grid_x, grid_y = np.meshgrid(xr, yr)

pr = pyproj.Proj(geodata["projection"])
grid_lon, grid_lat = pr(grid_x, grid_y, inverse=True)

Correction: a similar piece of code is used in io.exporters.initialize_forecast_exporter_netcdf. It's not too much effort to implement that in a separate utility function.

Hello @juanpablosimarro, @pulkkins already provided you with some good answers to your problem. Let us know should you need more support.

@pulkkins, maybe what you suggest could be flagged as a good first issue? @juanpablosimarro would you be interested in having a go at it?

Hello,
Thank you for your answers.
I will try today to use those solutions. If I encounter any problem, I will report here.
King regards
Juan

Hello again,
I had problems with the exporters and OPERA.
The code you provide above works well except that it is necessary to reverse the order of the y axis:
grid_x, grid_y = np.meshgrid(xr, np.flip(yr))
Thanks again!