The code below is what I came up with to solve the problem:

1. We're given geographic coordinates for two opposite vertices of a

rectangle.

2. If any of the sides of the rectangle is larger than some number of

degrees, then subdivide the rectangle into squares/rectangles such

that all sub-units have sides smaller than the specified amount.

3. Define each sub-unit by providing a list of all vertices, where the

first and last vertices are identical so as to close the polygon.

I've ignored the "if" part of the problem to simplify.

The key to my solution was to use numpy's meshgrid to generate the

coordinates for defining the sub-units. However, it seems awfully

complex and contrived, and am wondering if there's a simpler solution,

or perhaps some package offers this functionality. I couldn't find any,

so any tips appreciated.

---<--------------------cut here---------------start------------------->---

# Let's say we have these coordinates

lons = [-96, -51.4]

lats = [60, 72]

# And we want to create 10x10 (max) degree polygons covering the rectangle

# defined by these coordinates

step = 10

# Calculate how many samples we need for linspace. The ceiling is required

# to cover the last step, and then add two to accomodate for the inclusion

# of the end points... what a pain.

xn = np.ceil((lons[1] - lons[0]) / step) + 2

yn = np.ceil((lats[1] - lats[0]) / step) + 2

xgrd = np.linspace(lons[0], lons[1], xn)

ygrd = np.linspace(lats[0], lats[1], yn)

# Create grids of longitudes and latitudes with dimension (yn, xn). The

# elements of the longitude grid are the longitude coordinates along the

# rows, where rows are identical. The elements of the latitude grid are

# the latitude coordinates along the columns, where columns are identical.

longrd, latgrd = np.meshgrid(xgrd, ygrd, sparse=False)

for i in range(int(xn) - 1):

for j in range(int(yn) - 1):

print [(longrd[j, i], latgrd[j, i]), # lower left

(longrd[j, i + 1], latgrd[j, i]), # lower right

(longrd[j, i + 1], latgrd[j + 1, i]), # upper right

(longrd[j, i], latgrd[j + 1, i]), # upper left

(longrd[j, i], latgrd[j, i])] # close at lower left

---<--------------------cut here---------------end--------------------->---

--

Seb