Code covered by the BSD License

Highlights from Grid Sphere

4.66667
4.7 | 3 ratings Rate this file 24 Downloads (last 30 days) File Size: 10 KB File ID: #28842 Version: 1.15

Grid Sphere

Kurt von Laven (view profile)

26 Sep 2010 (Updated )

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

File Information
Description

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

Acknowledgements

Find Nearest Neighbors On Sphere, Golden Ratio, Spherical To Azimuthal Equidistant, and Geodesic Midpoints inspired this file.

Required Products MATLAB
MATLAB release MATLAB 7.9 (R2009b)
MATLAB Search Path
```/
/GridSphere```
17 Feb 2015 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.

Comment only
09 Feb 2015 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 :-)

Comment only
02 Jul 2013 Andrew

Andrew (view profile)

Fantastic file.

10 Jun 2013 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.

Comment only
06 Jun 2013 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?

14 Feb 2013 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.

Comment only
14 Feb 2013 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.

Comment only
25 Jan 2013 Kurt von Laven

Kurt von Laven (view profile)

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

Comment only
25 Jan 2013 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.

Comment only
24 Jan 2013 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.)

Comment only
24 Jan 2013 Kurt von Laven

Kurt von Laven (view profile)

Comment only
23 Jan 2013 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.

Comment only
19 May 2012 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.

Comment only
16 May 2012 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!

08 Feb 2012 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.

Comment only
07 Feb 2012 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;

Comment only
27 Sep 2010 1.1

27 Sep 2010 1.2

10 Oct 2010 1.3

Replaced the ElementWiseMax function with an equivalent built-in function.

08 Feb 2012 1.4

Added usage example per user request.

21 May 2012 1.5

Replaced tabs with spaces so that the source code displays consistently in all text editors.

14 Feb 2013 1.7

n/a

14 Feb 2013 1.8

Fixes a serious bug in RemoveEqualPairs that resulted in a small number of latitude, longitude pairs being omitted from the final output of GridSphere. Big thanks to Francisco for the bug report.

01 Apr 2013 1.9

A MATLAB user reported a crash on line 23 of RemoveEqualPairs.m. This change fixes the type mismatch error by converting a double to a uint32. GNU Octave users will be unaffected by this change as GNU Octave handles the type conversion automatically.

10 Mar 2015 1.10

Made code available as a toolbox and removed some unnecessary hidden files.

10 Mar 2015 1.11

13 Mar 2015 1.12

Removed trailing whitespace and switched from 2 spaces after each period to 1 in the interest of readability.

13 Mar 2015 1.13

Removed write and execute permissions from group/others on all *.m files.

13 Mar 2015 1.14

Moved some words in comments onto previous lines where they fit without exceeding 80 columns.

13 Mar 2015 1.15

Added a link to another FX submission that this submission uses code from.