# SPHERESORT is the jMap module of jLab.

``` SPHERESORT  Sorted great circle distances to nearby points on the earth.

Computing great circle distances on the earth between two sets of
points --- data points and a fixed grid, say --- is computationally
expensize.  SPHERESORT speeds this up substantially.

[DS,XS,YS]=SPHERESORT(LAT,LON,LATO,LONO,CUTOFF) returns the great
circle distances DS between data points at locations LAT, LON and
nearby grid points located at LATO, LONO, sorting in order of distance.

XS and YS are the corresponding positions of the sorted data points
in a local tangent plane about each grid point.

SPHERESORT only computes distances for nearby points.  CUTOFF is the
maximum distance (in kilometers) for which we wish the great circle
distance to be computed.

CUTOFF should be less that 1/4 of the circumference of the earth,
RADEARTH * PI/2, so that (LAT,LON) points lie in the same hemisphere
as the corresponding grid points.

LAT and LON are arrays of the same size into data point locations.
Any NaN values in LAT/LON are ignored.

LATO and LONO are arrays of length M and N, say, specifying the
latitudes and longitudes of an M x N matrix of grid points, i.e.

LATO =  [LATO_1;    LONO= [LONO_1 LONO_2 ... LONO_N].
LATO_2;
...
LATO_M]

The output arrays are then each M x N x P arrays of column vectors,
where P is the maximum number of points in the CUTOFF neighborhood at
any grid point.  Points farther away are filled with NANs.

DS gives the distances of data points less than CUTOFF kilometers
from the (m,n)th grid point, sorted in order of increasing distance.

XS and YS are also M x N x P, and give the coordinates of the
corresponding data points in a Cartesian expansion about the (m,n)th
grid point, in kilometers.  See LATLON2XY for details.
_________________________________________________________________

Let's say some additional variables Z1, Z2,...,ZK are given at
the data locations LAT, LON.  Then

[DS,XS,YS,Z1S,Z2S,...,ZKS]=
SPHERESORT(LAT,LON,Z1,Z2,...,ZK,LATO,LONO,CUTOFF);

also returns the sorted values of these variables.

Z1S, Z2S,...,ZKS are the same size as the other output arguments.
Z1S then gives the value of Z1 at data points no more than CUTOFF
kilometers from the (m,n)th grid point, etc.
_________________________________________________________________

One grid, many data fields

It is often the case that the field to be mapped, Z, consists of many
copies of observations at the same LAT/LON points.  For example, X and
Y could be 1-D arrays, and Z a matrix with LENGTH(Z) rows and many(K)
columns.  In this case, the format SPHERESORT(LAT,LON,Z1,Z2,...,ZK...)
described above is cumbersome, and may lead to memory problems.

Instead, one can pre-compute the DS, XS, and YS arrays once, and use
these for all K copies of Z. To do this, call SPHERESORT as

[DS,XS,YS,INDEX]=
SPHERESORT(LAT,LON,1:LENGTH(LAT(:)),LATO,LONO,CUTOFF)

which outputs an index INDEX into the data point locations, where INDEX
is the same size as the other output arrays.

Then for each of the K copies of Z at the same LAT/LON locations, form

ZS=NAN*ZEROS(SIZE(XS))
ZS(~ISNAN(INDEX))=ZK(INDEX(~ISNAN(INDEX)))

and then use this ZS together with DS, XS, and YS to call POLYSMOOTH.
In this way, SPHERESORT only needs to be called once and not K times.

When inputting LAT/LON values for which the corresponding data is
always undefined (e.g., altimeter tracks over land), once should set
the LAT/LON coordinates of these points to NaNs such that they will be
omitted from DS, XY, YS, and INDEX.
_________________________________________________________________

Parellel algorithm

With Matlab's Parallel Computing Toolbox installed, SPHERESORT can
take advantage of multiple cores to speed up operations.

SPHERESORT(...'parallel') parallelizes the computation with the
maximum possible number of workers for your system.  To set the number
of workers explicitly, use SPHERESORT(...,'parallel',NWORKERS).

As the additional efficiency is not dramatic, this is typically only an
advantage for very large datasets.  Due to the nature of the
calculation we have to use SPMD, which is less efficient than PARFOR.

As an example, for a dataset with 1 million points, a 12 core Mac Pro
takes about 44 seconds to complete the sort on a 1x1 degree grid,
versus 184 seconds for the standard algorithm, a factor of 4 faster.
_________________________________________________________________