sl.findnn.curvilin: Find Nearest Neighbour in a Gridded Field

In FESOM/spheRlab: Spherical Geometry, Analysis, and Plotting of Geoscientific Data on Arbitrary Grids

Find Nearest Neighbour in a Gridded Field

Description

Find the nearest neighbour of a gridded field relative to a specific point on a sphere.

Usage

sl.findnn.curvilin(lon = NULL, lat = NULL, x = NULL, y = NULL, z = NULL, lon.grid = NULL, lat.grid = NULL, x.grid = NULL, y.grid = NULL, z.grid = NULL, bruteforce = FALSE, Rsphere = 1)

Arguments

lon	a scalar or numeric vector of length L with the longitude(s) of the point(s) for which to find the nearest neighbour(s). Not required if (x.grid,y.grid,z.grid) are provided.
lat	a scalar or numeric vector of length L with the latitude(s) of the point(s) for which to find the nearest neighbour(s). Not required if (x.grid,y.grid,z.grid) are provided.
x	a scalar or numeric vector of length L with the x-coordinates of the point(s) for which to find the nearest neighbour(s). Note that (x,y,z) must be on a unit sphere, that is, sqrt(x^2+y^2+z^2)=1 must hold. (x,y,z) are not required if (lon,lat) are provided.
y	a scalar or numeric vector of length L with the y-coordinates of the point(s) for which to find the nearest neighbour(s). Note that (x,y,z) must be on a unit sphere, that is, sqrt(x^2+y^2+z^2)=1 must hold. (x,y,z) are not required if (lon,lat) are provided.
z	a scalar or numeric vector of length L with the z-coordinates of the point(s) for which to find the nearest neighbour(s). Note that (x,y,z) must be on a unit sphere, that is, sqrt(x^2+y^2+z^2)=1 must hold. (x,y,z) are not required if (lon,lat) are provided.
lon.grid	an MxN matrix with the longitudes of the points of a curvilinear grid. Not required if (x.grid,y.grid,z.grid) are provided.
lat.grid	an MxN matrix with the latitudes of the points of a curvilinear grid. Not required if (x.grid,y.grid,z.grid) are provided.
x.grid	an MxN matrix with the x-coordinates of the points of a curvilinear grid. Note that (x.grid,y.grid,z.grid) must be on a unit sphere, that is, sqrt(x.grid^2+y.grid^2+z.grid^2)=1 must hold. (x.grid,y.grid,z.grid) are not required if (lon.grid,lat.grid) are provided.
y.grid	an MxN matrix with the y-coordinates of the points of a curvilinear grid. Note that (x.grid,y.grid,z.grid) must be on a unit sphere, that is, sqrt(x.grid^2+y.grid^2+z.grid^2)=1 must hold. (x.grid,y.grid,z.grid) are not required if (lon.grid,lat.grid) are provided.
z.grid	an MxN matrix with the z-coordinates of the points of a curvilinear grid. Note that (x.grid,y.grid,z.grid) must be on a unit sphere, that is, sqrt(x.grid^2+y.grid^2+z.grid^2)=1 must hold. (x.grid,y.grid,z.grid) are not required if (lon.grid,lat.grid) are provided.
bruteforce	a logical value specifying whether to find the nearest neighbour(s) by computing the distance(s) from all grid points. If FALSE (default) and if max(M,N)>5, the nearest neighbour(s) is (are) found in an iterative manner (see Details), which can be much faster for large grids.
Rsphere	a scalar specifying the radius of the sphere to be used when computing the distance; default is Rsphere=1. This affects only the return value(s); if the point(s) or the grid are provided in Cartesian coordinates, these must correspond to a unit sphere in any case.

Details

If FALSE (default) and if max(M,N)>5, the nearest neighbour(s) is (are) found in an iterative manner by subsetting (coarsening) the grid such that distances are first computed only for 5x5 quasi-evenly distributed points (using the indices [round(seq(0,4)*((M-1)/4))+1),round(seq(0,4)*((N-1)/4))+1)]). In the next step, another 5x5 subset refined by a factor 2, centered on the nearest neighbour of the previous subset, is used. This procedure is continued until the subgrid resolution reaches the original grid resolution (in both dimensions), which will have 'zoomed' into the location of the overall nearest neighbour. This iterative procedure can be much faster for large grids because distances don't need to be computed from all grid points (scaling much better).

Value

a list with three elements:

nn.index	an Lx2 data frame with columns named 'row' and 'col' for the indices of the nearest neighbour(s). If L=1, it might be convenient in some cases to convert the data frame into a vector, e.g., using unlist().
nn.dist	a scalar or numeric vector of length L providing the distance from the respective nearest neighbour.

Author(s)

Helge Goessling

See Also

sl.grid2path

Examples

# generate an example curvilinear grid; first just lon-lat, then rotate (to around the North Pole) to make it more general
lon.seq = seq(-60,60,5)
lat.seq = seq(-30,30,5)
Nx = length(lon.seq)
Ny = length(lat.seq)
lon.0 = matrix(rep(lon.seq,Ny),nrow=Nx)
lat.0 = matrix(rep(lat.seq,each=Nx),nrow=Nx)
abg = sl.lonlatrot2abg(c(0,0,30))
rot.lonlat = sl.rot(lon = lon.0, lat = lat.0, alpha = abg[1], beta = abg[2], gamma = abg[3])
lon.grid = matrix(rot.lonlat$lon,nrow=Nx)
lat.grid = matrix(rot.lonlat$lat,nrow=Nx)

# define a few points, 'spiralling' to the North Pole
points.lon = seq(-135,180,45)
points.lat = seq(50,85,5)

# compute the nearest neighbours
nn = sl.findnn.curvilin(lon = points.lon, lat = points.lat, lon.grid = lon.grid, lat.grid = lat.grid)

# visualize the grid, points, and nearest neighbours in a map
pir = sl.plot.init(projection = "polar", do.init.device = FALSE)
sl.plot.naturalearth(pir, what = "land", resolution = "coarse", fill.col = "grey")
sl.plot.lonlatgrid(pir, col = "black", pole.hole = TRUE)
sl.plot.points(pir, lon = lon.grid, lat = lat.grid, pch = 16, cex = 0.3, col = "black")
sl.plot.points(pir, lon = points.lon, lat = points.lat, col = "red")
for (i in 1:length(points.lon)) {
  sl.plot.points(pir, lon = lon.grid[nn$nn.index$row[i],nn$nn.index$col[i]], lat = lat.grid[nn$nn.index$row[i],nn$nn.index$col[i]], col = "blue")
}
sl.plot.end(pir, do.close.device = FALSE)

FESOM/spheRlab documentation built on April 6, 2024, 6:52 p.m.