Dividing a geographic region

Question

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.

meshgrid(x1:deltaY:x2,y1:deltaX:y2)

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 =

    0.1799

>> 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;
end

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))

Answer 1

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.