How to find the nearest point on a grid?

General Usage
1

I have a grid defined as

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!

2

I think the Julia equivalent to the Python code would use Interpolations.jl:

using Interpolations

itp = interpolate((lat, lon), data, Gridded(Constant()))
res = itp(some_new_point)
1 Like
3

Maybe this is a simple rounding problem?

  • round(x) “round x to the nearest integer”
julia> Int8(round(101.02)) , Int8(round(24.97))
(101, 25)

EDIT:

  • my proposal is not working for the Antimeridian ; so need some post-checking …

1 Like
4

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.

2 Likes
5

Many thanks all for the replies!

Actually, my 0.25 example is a simplified one. I have grids like the below:

  -53.6250
  -53.3750
  -53.1250
  -52.8750
  -52.6250
  -52.3750
  -52.1250
  -51.8750
  -51.6250
  -51.3750
6

so complex …

julia> (round( (  -52.1251  - 0.1250) *4*1000 ) / (4*1000) ) + 0.125
-52.125

julia> (round( (  -52.1249  - 0.1250) *4*1000 ) / (4*1000) ) + 0.125
-52.125
2 Likes
7

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.

Something like the below?

Ind = Lon .< 101.02;
Lon1 = Lon[Ind][end];
Lon2 = Lon1+0.25
8

If the spacing is regular, the mapping can always be expressed in integer multiples of the grid spacing, wouldn’t you agree?

9

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.

Here is my code:

Ind0 = findfirst(Lon_grid .> Lon1[k]);
L2 = Lon_grid[Ind0]; 
L1 = Lon_grid[Ind0-1];
if abs(L2 - Lon1[k]) > abs(Lon1[k]-L1)
    Lon_id = Ind0-1;
else
    Lon_id = Ind0;
end
11

Example

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
1 Like
12

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.

1 Like