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. _________________________________________________________________ Additional input parameters 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. _________________________________________________________________ See also TWODSORT, POLYSMOOTH. 'spheresort --t' runs a test. Usage: ds=spheresort(lat,lon,lato,lono,cutoff); ds=spheresort(lat,lon,lato,lono,cutoff,'parallel',Nworkers); [ds,xs,ys]=spheresort(lat,lon,lato,lono,cutoff); [ds,xs,ys,z1s,z2s]=spheresort(lat,lon,z1,z2,lato,lono,cutoff); __________________________________________________________________ This is part of JLAB --- type 'help jlab' for more information (C) 2008--2017 J.M. Lilly --- type 'help jlab_license' for details