Lon = [-180:0.25:180];
Lat = [-90:0.25:90];
X, Y = meshgrid(Lon, Lat);

I have many pairs of station data with corresponding longitude and latitude values and my goal is to find the nearest grid point for each station. For example, the results for a station with longitude of 101.02, and latitude of 24.97, would be [101, 25].

I did a google search and it seems that in Python, people can do this:

from scipy.interpolate import RegularGridInterpolator
itp = RegularGridInterpolator( (lat, lon), data, method='nearest')
res = itp(some_new_point)

I wonder what would be the Julia solution? Many thanks!

In that example yes, but in the original post, the grid points vary by 0.25 degrees, not integers. Could be easy to solve if the mesh is regular, of course.

Thanks, Steven! Does this mean I would have to interpolate my lon lat pairs as if it is on a smooth line? In reality, my data points can be quite random and unrelated.

Is there a function that allows me to find the last value in [-180:0.25:180] that is smaller than my longitude? For example, if my longitude is 101.02, the value I’m looking for would be 101. Then I can find out whether my value is closer to 101 or closer to the next value 101.25, and decide which one to use.

Absolutely. However, in global models, it is common to use very complex grids to highlight regions with large changes. That said, I’m able to use the function findfirst for my needs.

using GMT
# Have to add *z* too because this program expects *x y z*
blockmean([11.4 11.6 0], region=(10,15,10,15), spacing=1, center=true)
Vector{GMTdataset} with 1 segments
First segment DATA
Global BoundingBox: [11.0, 11.0, 12.0, 12.0, 0.0, 0.0]
First seg BoundingBox: [11.0, 11.0, 12.0, 12.0, 0.0, 0.0]
1×3 Matrix{Float64}:
11.0 12.0 0.0

I posted the Interpolations.jl code based on the Python code you originally posted, but upon further inspection it looks like you are not interpolating but discretizing. So what I posted seems irrelevant to what you’re looking for.