I've been trying to generate an evenly spaced grid centred at a given point on a sphere, such that the angular separation between any neighbouring pair of points is the same (e.g., 1 degree). The grid should be oriented with the central column of points along the same meridian as the given point. I have little experience with spherical trigonometry (aside from knowing about the Haversine formula and the spherical law of cosines) but after some reading I'm fairly certain that each cell will require four unique great circle arcs. The previous approach I followed using vectors only found great circles perpendicular to the central column (like the coloured grids shown in the first plot) and so the points were not evenly spaced, especially at the edges. Each point in the second plot (showing the result of my latest attempt) is plotted with a 1 degree-radius circle surrounding it, which ideally should coincide with the surrounding four points. The updated method described below finally generates a regular grid. However, for a 20-by-20 degree grid the angular separation of the corner point to the central column is ~4.5 arcseconds smaller than the ideal 10 degree separation. Does anyone know of a method that can determine the required grid points to greater precision?
Current method:
The procedure I've followed so far to generate the n-by-n degree grid centred at a given point (e.g. longitude=0, latitude=0; let's call this point A) involves three stages:
Note: a 10-by-10 degree grid would consist of 11-by-11 points each separated from the nearest four points in the adjacent rows and columns by 1 degree. All calculations assume a sphere of radius 1.
-
Find the points along the `central' column:
Simply increment the latitude (i.e. following the meridian running through the centre point) and wrap the latitude and longitude if a pole is crossed (i.e. if the current point's latitude is less than -90 or greater than +90). -
Find the points along the `central' row (through the grid centre):
Since the row of points through A follow a great circle that intersects the meridian at a right angle, we need to find the plane that this second great circle lies on. Consider two points at the given central longitude; let point A be the grid centre at 0 degrees longitude and 0 degrees latitude, and point B at 0 degrees longitude and 1 degree north. The two points can be expressed in Cartesian coordinates, and so the unit vector (N) normal to the great circle that contains A and B is found using the cross product, i.e. N = BxA. The cross product of the unit vector N with a unit vector in the direction of A gives another normal vector (M). The vector M can then be used as an axis to rotate N by (90-x) degrees^, where x is the offset (in this case 1 deg) from the central column. This is repeated outwards in steps of 1 degree for one half of the central row then the other half using a unit vector in the opposite direction to N.
^http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle -
Fill in the remaining points in each quadrant:
Each point in the grid forms a rectangle lying on the same plane as three other points in the grid – let A be the point at the grid centre, B is a point 'above' or 'below' A in the same column and C is a point on the central row, and the unknown point D is in the same row as B and same column as C. For the upper left quadrant, the rectangle (sometimes a square) would have corners clockwise from the upper-left D, B, A, and C.
First compute vectors AB and AC then a unit vector normal to the triangle BAC can be found from their cross-product, i.e. t=ABxAC. To find D, the approach I've used is to rotate the plane on which A/B/C/D lie to align it with the xy-plane. The rotation matrix requires a rotation axis and rotation angle that I found here: https://gamedev.stackexchange.com/a/48102
The dot product of the rotation matrix and each of vectors A, B and C allows vector D' (the equivalent of vector D in the xy-plane) to be found from their components,
i.e., $D' = (b_x-a_x+c_x, b_y-a_y+c_y, b_z-a_z+c_z)$
The dot product of the rotation matrix using the same rotation axis but rotation angle of opposite sign to that given above with vector D' gives the coordinates of vector D, which can be converted to the latitude and longitude for point D. This is repeated for every remaining point in the four quadrants.
I've spent a few days looking for a better method but without any luck.
Fig 1. Sphere with grids:
Fig 2. Grid points – note circles shown are 1 degree in radius (i.e. the perimeter of a 1 degree spherical cap centred at the given grid point). I had used a gnomonic projection to present the grid in figure 2 (using the python plotting package APLpy) but replaced it with an orthographic projection since the former shows great circles as straight lines but in this case the rows of points are joined by different great circle arcs, except the (n/2+1)th row.
Best Answer
You can't do what you are hoping because you can't cover a sphere with squares. In fact, you can't even put one square on the sphere. If you look closely, the angles of your squares add up to more than $360^\circ$, by an amount that is determined by the solid angle taken up by the square.
A few choices that maintain some of what you want. 1) The black grid in your first figure comes from the classic lines of latitude and longitude. You get quadrilaterals except near the poles. You can choose what latitude to make the cells almost square, but they will be much larger near the equator than near the poles. 2) Use the latitude/longitude lines as in 1, but have fewer cells as you get near the poles. You can keep the cells close to square, but won't have a nice grid. 3) Use the square grid as you are, but recognize that you can't have one grid that covers the sphere. Maybe a grid that covers about half the sphere centered on a point of current interest meets your needs.