File Exchange

## Grid Sphere

version 1.15 (10 KB) by Kurt von Laven

### Kurt von Laven (view profile)

Produces a nearly even grid over the surface of a sphere.

Updated 13 Mar 2015

The grid may have a total of exactly 12, 42, 162, 642, ... points. Mathematically speaking, the grid may have either 12 or 2 + ( 10 * (4 ^ k) ) points, where k is a positive integer. The user may request any number of points, and the closest attainable value will be produced. All code is compatible with GNU Octave. The algorithm was developed by Nick A. Teanby of Oxford University. Refer to his website for the publication describing the approach and more elaborate geodesic grid software written in IDL: http://www.atm.ox.ac.uk/user/teanby/software.html#icos. Use the FindNearestNeighbors function, available on the MATLAB file exchange at http://www.mathworks.com/matlabcentral/fileexchange/28844-find-nearest-neighbors-on-sphere, to find the grid points closest to arbitrary query points.
GridSphere and FindNearestNeighbors share some functions in common. Each package contains a copy of these functions so that both can stand alone. To eliminate duplicates, simply move all the files into a single folder and replace the shared files when prompted.
Sample usage:
[latGridInDegrees, longGridInDegrees] = GridSphere(12)
latGridInDegrees =
-58.28253
-58.28253
-31.71747
-31.71747
0.00000
0.00000
0.00000
0.00000
31.71747
31.71747
58.28253
58.28253
longGridInDegrees =

0.00000
180.00000
-90.00000
90.00000
-148.28253
-31.71747
31.71747
148.28253
-90.00000
90.00000
0.00000
180.00000

### Cite As

Kurt von Laven (2019). Grid Sphere (https://www.mathworks.com/matlabcentral/fileexchange/28842-grid-sphere), MATLAB Central File Exchange. Retrieved .

Bobbie WU

### Bobbie WU (view profile)

W Ryan Williamson

### W Ryan Williamson (view profile)

Elapsed time: 0.0035 seconds for 642 points. Seems plenty fast. Suited my needs perfectly, along with Find Nearest Neighbors (http://www.mathworks.com/matlabcentral/fileexchange/28844-find-nearest-neighbors-on-sphere)

Peter Gagarinov

### Peter Gagarinov (view profile)

spheretri is much faster than GridSphere, see http://www.mathworks.com/matlabcentral/fileexchange/58453-spheretri

sumana

Kurt von Laven

### Kurt von Laven (view profile)

Hello, Peter! Sorry for the delayed response. Yes, it is intentional. Also, the formula you linked to assume the following coordinate system: http://en.wikipedia.org/wiki/File:3D_Spherical.svg, where r corresponds to rho. The theta computed by arccos(z / ((x ^ 2 + y ^ 2 + z ^ 2) ^ (1 / 2))) is not the latitude, but rather the angle between the positive z-axis and the point in question. Also, acosd (http://www.mathworks.com/help/matlab/ref/acosd.html) returns results in the range [0, 180] rather than the range [-90, 90] returned by http://www.mathworks.com/help/matlab/ref/asind.html that we want for longitude.

Peter Markowsky

### Peter Markowsky (view profile)

Hi Kurt,

Is the use of the arcus sinus instead of the tradition arcus cosinus (see
http://en.wikipedia.org/wiki/List_of_common_coordinate_transformations#From_Cartesian_coordinates)
for sphere coordinates inentional? It created some confusion for me :-)

Andrew

Fantastic file.

Kurt von Laven

### Kurt von Laven (view profile)

Hi, Ben. Thanks for rating. I haven't tested it, but see the first comment on this file (from Paul on February 7th, 2012). I believe this will work if you're using MATLAB but possibly not in GNU Octave.

Ben

### Ben (view profile)

Hi Kurt, thanks a lot for the sharing! I have a naive question. How to visualize the result? Could you also kindly provide some examples on using this code?

Kurt von Laven

### Kurt von Laven (view profile)

The fix has been approved and is now available for download. Thanks again for catching the bug, Francisco.

Kurt von Laven

### Kurt von Laven (view profile)

Finally had time to come up with a more elegant fix and submitted it for review. Will post again once the submission is approved.

Kurt von Laven

### Kurt von Laven (view profile)

Sweet; glad to hear it. Yeah, my fix is pretty brutish too.

Francisco

### Francisco (view profile)

Dont worry !, y already fixed it (in a quite brute way, but it's ok) just wanted to report the bug.

Kurt von Laven

### Kurt von Laven (view profile)

Thank you for the clear bug report, Francisco. I don't have MATLAB on hand at the moment, but I have an untested fix ready. If you like I can email it straight to you. Either way, I'll make sure I didn't introduce any additional bugs and then go through the normal MATLAB File Exchange upload process. (It usually takes them a few days to review submissions.)

Kurt von Laven

Francisco

### Francisco (view profile)

Hi Kurt ! awesome comments in the code, i could understand it entirely.

just one thing, i found something i believe is an error, if i calculate de size of the results (i.e. length(GridSphere(n))

im getting this:
>> length(GridSphere(42)) = 41

>> length(GridSphere(162)) = 161

>> length(GridSphere(600)) = 641

the ammount of points are not correct, so i entered the code and found this at the 'RemoveEqualPairs' function:

uniquePairIndices = DistinctFromPrevious(sortedTwos);
% Assuming the greatest two in each cluster does not equal the least two in
% the next cluster, the pairs for which the two does not equal the previous
% two are the unique pairs.

the asumption there is wrong i believe, for instance for the 42 point grid, when i get the sorted points (with repeated points), the first ones are:

>> [sortedOnes, sortedTwos]

ans =

-90.0000 0
-90.0000 0
-58.2825 0
-58.2825 0
-58.2825 180.0000
-54.0000 -121.7175
... ....

so the 'cleaning' function 'DistinctFromPrevious' should be changed, in particular it must consider the 'sortedOnes' vector.

That's all,
Regards,
Francisco.

Kurt von Laven

### Kurt von Laven (view profile)

Thank you, LU! That's a great question. How concerned are you concerned with performance? If you can stomach an O(n ^ 2) algorithm, where n is the number of points output by GridSphere, then I think the following approach will be easiest.

1. Find the distance between all the pairs of points output by GridSphere.
2. Sort by these distances to find all the pairs of points that are closest to one another. There are at least three gotchas here: (1) you'll need to make sure that either you're not computing the distance between a point and itself or that you're filtering out the pairs with distance 0, (2) you'll probably want to avoid computing both the distance from point a to point b and the distance from point b to point a, and (3) you'll need to account for the imperfect precision of floating point arithmetic, which will result in clumps of distances that are very close but not quite equal. Some of the functions that come with GridSphere may be helpful to you, such as DistinctFromPrevious, which will find the boundaries between the aforementioned clumps.
3. Draw great circles (or just straight lines) connecting the pairs of points you found in step 2. Each of these pairs corresponds to one edge of a triangle. As Paul mentioned, you may need to call conversion functions like deg2rad and sph2cart on your great circle (or line) vectors to project them appropriately.

LU

### LU (view profile)

Hi Kurt, Wonderful program. I was just wondering if it is possible to connect the vertices with lines? Using Paul's example plots the vertiecs but i would like to display it as a mesh grid showing the edges of the triangles.

Thanks!

Kurt von Laven

### Kurt von Laven (view profile)

Thank you for the feedback, Paul. I went ahead and put in a request for the description to be modified to include a usage example. It should appear within five business days.

Paul

### Paul (view profile)

Blazingly fast and well documented. Only an example would not have been unbecoming, e.g. to illustrate the fact that returned values are in degrees (whereas Matlab's sph2cart requires radians).

[lat,long] = GridSphere(1000);
scatter3(x, y, z, 22);
axis equal vis3d;