A program to identify neighboring countries

For our first real geospatial analysis program, we are going to write a Python script that identifies neighboring countries. The basic concept is to extract the polygon or multipolygon for each country and see which other countries each polygon or multipolygon touches. For each country, we will display a list of other countries that border that country.

Let's start by creating the Python script. Create a new file named borderingCountries.py and place it in the same directory as the TM_WORLD_BORDERS-0.3.shp shapefile you downloaded earlier. Then enter the following into this file:

import osgeo.ogr
import shapely.wkt

def main():
    shapefile = osgeo.ogr.Open("TM_WORLD_BORDERS-0.3.shp")
    layer = shapefile.GetLayer(0)

    countries = {} # Maps country name to Shapely geometry.

    for i in range(layer.GetFeatureCount()):
        feature = layer.GetFeature(i)
        country = feature.GetField("NAME")
        outline = shapely.wkt.loads(feature.GetGeometryRef().ExportToWkt())

        countries[country] = outline

    print "Loaded %d countries" % len(countries)

if __name__ == "__main__":
    main()

So far, this is pretty straightforward. We are using the techniques we learned earlier to read the contents of the shapefile into memory and converting each country's geometry into a Shapely object. The results are stored in the countries dictionary. Finally, notice that we've placed the program logic into a function called main()—this is good practice as it lets us use a return statement to handle errors.

Now run your program just to make sure it works:

$ python borderingCountries.py
Loaded 246 countries

Our next task is to identify the bordering countries. Our basic logic will be to iterate through each country and then find the other countries that border this one. Here is the relevant code, which you should add to the end of your main() function:

    for country in sorted(countries.keys()):
        outline = countries[country]

        for other_country in sorted(countries.keys()):

            if country == other_country: continue

            other_outline = countries[other_country]

            if outline.touches(other_outline):

                print "%s borders %s" % (country, other_country)

As you can see, we use the touches() method to check if the two countries' geometries are touching.

Running this program will now show you the countries that border each other:

$ python borderingCountries.py
Loaded 246 countries
Afghanistan borders Tajikistan
Afghanistan borders Uzbekistan
Albania borders Montenegro
Albania borders Serbia
Albania borders The former Yugoslav Republic of Macedonia
Algeria borders Libyan Arab Jamahiriya
Algeria borders Mali
Algeria borders Morocco
Algeria borders Niger
Algeria borders Western Sahara
Angola borders Democratic Republic of the Congo
Argentina borders Bolivia
...

Congratulations! You have written a simple Python program to analyze country outlines. Of course, there is a lot that could be done to improve and extend this program. For example:

• You could add command-line arguments to let the user specify the name of the shapefile and which attribute to use to display the country name.
• You could add error checking to handle invalid and non-existent shapefiles.
• You could add error checking to handle invalid geometries.
• You could use a spatial database to speed up the process. The program currently takes about a minute to complete, but using a spatial database would speed that up dramatically. If you are dealing with a large amount of spatial data, properly indexed databases are absolutely critical or your program might take weeks to run.

We will look at all these things later in the book.