lat/lon of grid points

[juanpablosimarro]

Hello,

I wonder if there is a method to get the longitudes and latitudes of all grid points. From this code:

importer = io.get_method(importer_name, "importer")
Z, _, metadata = io.read_timeseries(fns, importer, **importer_kwargs)

Given Z and metadata from the file (HDF5 or other format), is there a method to get the longitude and latitude of the grid points of Z given the information kept in metadata?

Thank you

Juan Simarro

[pulkkins]

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)

[pulkkins]

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.

[dnerini]

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?

[juanpablosimarro]

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

[juanpablosimarro]

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!