Before we can show the selected country on a map, we need to calculate the bounding box for that country, that is, the minimum and maximum latitude and longitude values. Knowing the bounding box allows us to draw a map centered over the desired country.

Add the following to the end of your selectArea.py script:

cursor.execute("SELECT name," +
               "ST_YMin(ST_Envelope(outline))," +
               "ST_XMin(ST_Envelope(outline))," +
               "ST_YMax(ST_Envelope(outline))," +
               "ST_XMax(ST_Envelope(outline)) " +
               "FROM countries WHERE id=%s",
               (countryID,))

row = cursor.fetchone()
if row != None:
    name     = row[0]
    min_lat  = row[1]
    min_long = row[2]
    max_lat  = row[3]
    max_long = row[4]
else:
    print(HEADER)
    print('<b>Missing country</b>')
    print(FOOTER)
    sys.exit(0)

minLong = minLong - 0.2
minLat = minLat - 0.2
maxLong = maxLong + 0.2
maxLat = maxLat + 0.2

As you can see, we extract the country's bounding box from the database and then increase the bounds slightly so that the country won't butt up against the edge of the map. We also grab the name of the country at the same time, as we'll need this shortly.

Because countries aren't square, we can't simply draw the country onto a map of a fixed size. Instead, we have to calculate the aspect ratio of the country's bounding box (its width as a proportion of its height) and then calculate the size of the map image based on this aspect ratio, while limiting the overall size of the image so that it can fit within our web page. Here's the necessary code, which you should add to the end of your selectArea.py script:

width = float(maxLong - minLong)
height = float(maxLat - minLat)
aspectRatio = width/height

mapWidth = MAX_WIDTH
mapHeight = int(mapWidth / aspectRatio)

if mapHeight > MAX_HEIGHT:
    # Scale the map to fit.
    scaleFactor = float(MAX_HEIGHT) / float(mapHeight)
    mapWidth = int(mapWidth * scaleFactor)
    mapHeight = int(mapHeight * scaleFactor)

Doing this means that the map is correctly sized to reflect the dimensions of the country we are displaying.

With the bounding box, the map's dimensions, and the data source all set up, we are finally ready to render the map into an image file. This is done using a single function call, as follows:

hilite = "[id] = " + str(countryID)
imgFile = mapGenerator.generateMap("countries",
                                   minLong, minLat,
                                   maxLong, maxLat,
                                   mapWidth, mapHeight,
                                   hiliteExpr=hilite)

Notice that we set the hiliteExpr parameter to the value [id] = "+str(countryID). This visually highlights the country with the given ID.

The mapGenerator.generateMap() function returns a reference to a PNG-format image file containing the generated map. This image file is stored in a temporary directory, and the file's relative pathname is returned to the caller. This allows us to use the returned imgFile directly within our CGI script, like this:

print(HEADER)
print('<b>' + name + '</b>')
print('<p>')
print('<form method="POST" action="showResults.py">')
print('Select all features within')
print('<input type="text" name="radius" value="10" size="2">')
print('miles of a point.')
print('<p>')
print('Click on the map to identify your starting point:')
print('<br>')
print('<input type="image" src="' + imgFile + '" ismap>')
print(HIDDEN_FIELD.format("countryID", countryID))
print(HIDDEN_FIELD.format("countryName", name))
print(HIDDEN_FIELD.format("mapWidth", mapWidth))
print(HIDDEN_FIELD.format("mapHeight", mapHeight))
print(HIDDEN_FIELD.format("minLong", minLong))
print(HIDDEN_FIELD.format("minLat", minLat))
print(HIDDEN_FIELD.format("maxLong", maxLong))
print(HIDDEN_FIELD.format("maxLat", maxLat))
print('</form>')
print(FOOTER)

<input type="hidden" name="..." value="...">

The use of <input type="image" src="..." ismap> in this CGI script has the interesting effect of making the map clickable: when the user clicks on the image, the enclosing HTML form will be submitted with two extra parameters named x and y. These contain the coordinate within the image that the user clicked on.

This completes the selectArea.py CGI script. If you are running this on Mac OS X or a Linux machine, you will need to make sure you added an appropriate shebang line to the start of your program and made it executable, as described earlier, so that it can run as a CGI script.

All going well, you should be able to point your web browser to http://127.0.0.1:8000/cgi-bin/selectCountry.py.

When you select a country, you should see a map of that country displayed within your web browser. If you click within the map, you'll get a 404 error, indicating that the final CGI script hasn't been written yet.

The selectArea.py script generates an HTML form that is submitted when the user clicks on the low-resolution country map. The showResults.py script receives the form parameters, including the x and y coordinates of the point the user clicked on.

By itself, this coordinate isn't very useful, as it's simply the x and y offset, measured in pixels, of the point the user clicked on. We need to translate the submitted (x, y) pixel coordinate into a latitude and longitude value corresponding to the clicked-on point on the earth's surface.

To do this, we need to have the following information:

These variables were all calculated in the previous section and passed to us using hidden form variables, along with the country name and ID, the desired search radius, and the (x,y) coordinate of the clicked-on point. We can retrieve all of these using the cgi module. To do this, add the following to the end of your showResults.py file:

form = cgi.FieldStorage()

countryID   = int(form['countryID'].value)
radius      = int(form['radius'].value)
x           = int(form['x'].value)
y           = int(form['y'].value)
mapWidth    = int(form['mapWidth'].value)
mapHeight   = int(form['mapHeight'].value)
countryName = form['countryName'].value
minLong     = float(form['minLong'].value)
minLat      = float(form['minLat'].value)
maxLong     = float(form['maxLong'].value)
maxLat      = float(form['maxLat'].value)

With this information, we can now calculate the latitude and longitude that the user clicked on. We start by calculating how far across the image the user clicked, as a number in the range from 0 to 1:

xFract = float(x)/float(mapWidth)

An xFract value of 0.0 corresponds to the left side of the image, while an xFract value of 1.0 corresponds to the right side of the image. We then combine this with the minimum and maximum longitude values to calculate the longitude of the clicked-on point:

longitude = minLong + xFract * (maxLong-minLong)

We then do the same to convert the Y coordinate into a latitude value:

yFract = float(y)/float(mapHeight)
latitude = minLat + (1-yFract) * (maxLat-minLat)

Notice that we are using (1-yFract) rather than yFract in the preceding calculation. This is because the minLat value refers to the latitude of the bottom of the image, while a yFract value of 0.0 corresponds to the top of the image. By using (1-yFract), we flip the values vertically so that the latitude is calculated correctly.

Let's review what we have achieved so far: the user has selected a country, viewed a simple map of the country's outline, entered the desired search radius, and clicked on a point on the map to identify the origin for the search. We then converted this clicked-on point into a latitude and longitude value.

All of this provides us with three numbers: the desired search radius and the latitude and longitude of the point at which to start the search. Our task now is to identify which place names are within the given search radius of the clicked-on point:

Fortunately, because our places table uses the GEOGRAPHY column type, we can ask PostGIS to do all the hard work for us. Here is the necessary code to find the matching place names:

cursor.execute("SELECT ST_X(position::geometry)," +
               "ST_Y(position::geometry),name " +
               "FROM places WHERE " +
               "ST_DWithin(position, ST_MakePoint(%s, %s), %s)",
               (longitude, latitude, radius * METERS_PER_MILE))

points = []
for name,long,lat in cursor:
    points.append([long, lat, name])

The matching place names will be stored in the points variable, as a list of (longitude, latitude, name) tuples.

Now that we have calculated the list of place names within the desired search radius, we can use the mapGenerator.py module to display them. To do this, we're going to have to calculate the bounding box for the area to display, based on the clicked-on point, the supplied search radius, and a margin to ensure that we don't draw place names right on the edge of the map:

To do this, we'll use pyproj to calculate the position of four coordinates directly north, south, east, and west of the starting point:

How far north, south, east, or west we go depends on the search radius, but we'll also add a 10% margin to this distance. Using pyproj, this calculation is quite straightforward:

geod = pyproj.Geod(ellps="WGS84")
distance = radius * METERS_PER_MILE * 1.1

x,y,angle = geod.fwd(longitude, latitude, 0, distance)
maxLat = y

x,y,angle = geod.fwd(longitude, latitude, 90, distance)
maxLong = x

x,y,angle = geod.fwd(longitude, latitude, 180, distance)
minLat = y

x,y,angle = geod.fwd(longitude, latitude, 270, distance)
minLong = x

Once this has been done, we can use our mapGenerator.py module to generate a map showing the results we have found:

imgFile = mapGenerator.generateMap("shorelines",
                                   minLong, minLat,
                                   maxLong, maxLat,
                                   600, 600,
                                   points=points)

When we called the map generator previously, we used a filter expression to highlight particular features. In this case, we don't need to highlight anything. Instead, we pass it the list of place names to display on the map in the keyword parameter named points.

The map generator creates a PNG-format file and returns a reference to that file. Let's finish our program by displaying the map so the user can see it:

print(HEADER)
print('<b>' + countryName + '</b>')
print('<p>')
print('<img src="' + imgFile + '">')
print(FOOTER)

This completes our first attempt at the entire DISTAL application. Let's now see how it runs.