Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a certain geographic region defined by the bottom left and top right coordinates. How can I divide this region into areas of 20x20km. I mean in practial the shape of the earth is not flat it's round. The bounding box is just an approximation. It's not even rectangular in actual sense. It's just an assumption. Lets say the bottomleft coordinate is given by x1,y1 and the topright coordinate is given by x2,y2, the length of x1 to x2 at y1 is different than that of the length between x1 to x2 at y2. How can I overcome this issue

Actually, I have to create a spatial meshgrid for this region using matlab's meshgrid function. So that the grids are of area 20x20km.


As you can see I can have only one deltaX and one deltaY. I want to choose deltaX and deltaY such that the increments create grid of size 20x20km. However this deltaX and deltaY are supposed to vary based upon the location. Any suggestions?

I mean lets say deltaX=del1. Then distance between points (x1,y1) to (x1,y1+del1) is 20km. BUt when I measure the distance between points (x2,y1) to (x2, y1_del1) the distance is < 20km. The meshgrid function above does creates mesh. But the distances are not consistent. Any ideas how to overcome this issue?

For finding point p(1,0)

deg = km2deg(20)

deg =


>> reckon(25,-120,0.1799,-18)

ans =

   25.1711 -120.0614

Function to calculate distance between two points

function [ distance ] = calculateDistance( latitude1,longitude1,latitude2,longitude2 )
radius = 6371;
dLat = degtorad(latitude2-latitude1);
dLon = degtorad(longitude2-longitude1);
a = sin(dLat/2) * sin(dLat/2) + cos(degtorad(latitude1)) * cos(degtorad(latitude2)) * sin(dLon/2) * sin(dLon/2);
c = 2 * atan2(sqrt(a), sqrt(1-a));
distance = radius * c;

Referring to whuber's suggestions

syms p(1,1).lat p(1,1).lon
solve ( calculateDistance(p(0,1).lat,p(0,1).lon,p(1,1).lat,p(1,1).lon), calculateDistance(p(1,0).lat,p(1,0).lon,p(1,1).lat,p(1,1).lon))
share|improve this question
up vote 3 down vote accepted

This can be accomplished in the same way we would lay out a grid of squares in the plane:

  1. Lay off a baseline in any direction. Mark it off in equal increments (of 20 km). Call these points p(0,0), p(0,1), ..., p(0,n), in order. It does not have to be geodesically straight, but it should be close to straight.

  2. Starting with the first marked point p(0,0), choose a direction beta not parallel to the baseline (perhaps perpendicular to the baseline, for instance). Move 20 km in that direction and stop. Call this point p(1,0).

  3. Begin with the first two points on the baseline, p(0,0) and p(0,1). There will be two locations that are 20 km from each of p(0,1) and p(1,0): these are the intersections of circles of radius 20 km centered at p(0,1) and p(1,0). One of these points, by construction, is p(0,0). Define p(1,1) to be the other point.

    Iteratively determine p(1,j), j = 2, ..., n, in this fashion.

  4. Repeating steps 2 and 3, iteratively determine as many rows as desired, using the most recent row p(i,0), p(i,1), ..., p(i,n) as the baseline and keeping beta constant (at least approximately).

Retain only those grid points falling within one's region of interest.

The earth's curvature begins having a pronounced effect on grids more than a thousand kilometers on a side (or so). To illustrate this, here is a grid with 400 km spacing constructed to cover the conterminous United States.


This is a Mercator projection. The baseline of the grid is the bottom row of points, starting at (lat, lon) = (25, -120) and proceeding along the geodesic oriented 89.5 degrees east of north. The angle 'beta' is 18 degrees west of north: that determines the orientation of the left side of the grid.

(There's a limit to how far one can go with this approach: eventually such grids will wrap around the globe and return to cover the same areas. They will never match up unless they are extremely large (and effectively create a cube, octahedron, or prism).)

Note that apart from the aforementioned exceptions, it is impossible to create a grid of congruent squares (or any other shape) on the sphere (or any ellipsoid). This is because the congruencies would generate a discrete subgroup of the group of rotations of space and those subgroups correspond to the five Platonic solids, prisms, and pyramids. So, something has to give: if you want the grid cells to have equal areas or equal angles, the lengths of their sides must vary; if you want the side lengths to be the same, then the areas and angles must vary. This is why the cell shapes appear so distorted and non-square in the figure: such distortion is unavoidable.

share|improve this answer
Hi Whuber, Thanks for the explanation. However, I tried what you have showed here before. My question is how come the point p(1,1) is at distance 20km from each of the points p(0,1) and p(1,0). I tried with the base row as you said with a distance of 20km. Then I tried to get the next row. But point p(1,1) is not at a distance of 20km. I mean lets say for the base row I define delta_lon so that point p(0,0) and point (0,1) have same lat but lon different by delta_lon which makes it 20km. Now if I use the same delta_lon in the next row, the distance is not equal to 20km – rajan sthapit Jun 4 '12 at 17:49
My question is how can I get the point p(1,1). I mean as you said the point is the intersection of the circles at point p(1,0) and p(0,1). How can I get that in matlab. Are there any specific functions. Also how did you get the angle beta. I mean I tried in the same longitude as point p(0,0) to get the point p(1,0) at a distance of 20km. I need some more details – rajan sthapit Jun 4 '12 at 17:55
Rajan, As you point out, the crux of the matter is finding p(1,1) given p(1,0), p(0,1), and p(0,0). You can solve this numerically: that is, solve the simultaneous (nonlinear) equations distance(p(1,0), p(0,0)) = 20 km = distance(p(0,1), p(0,0)) and remove the solution equal to (or extremely close to) p(0,0). This requires (a) solving simultaneous nonlinear equations and (b) finding distances on the sphere or ellipsoid. – whuber Jun 4 '12 at 18:51
Hi Whuber, Thanks I will try that. But my other question is how you found out beta. Can I choose any beta I want. I mean in my case for finding point p(1,0) at a distance of 20km from point(0,0), can I use the direction of the same longitude as point(0,0) has. I mean how come you got 18 degrees west of north in your case. Did you follow a specific procedure? – rajan sthapit Jun 4 '12 at 18:57
Whuber, I have made some changes in the question itself. Could you suggest me if it is the right way to go? I am pretty confused how to exactly use the and In the distance function, I cannot actually get the equations to pass the solve function. So I thought that I should write my own distance function. – rajan sthapit Jun 4 '12 at 20:10

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.