GeoDjango Tutorial

Introduction

GeoDjango is an add-on for Django that includes support for geometry fields and extends the ORM to allow spatial queries. This tutorial assumes a familiarity with Django; thus, if you’re brand new to Django please read through the regular tutorial to introduce yourself with Django concepts. GeoDjango has special prerequisites over what is required by Django – please consult the GeoDjango installation documentation for more details.

This tutorial is going to guide you through guide the user through the creation of a geographic web application for viewing the world borders. [1] Some of the code used in this tutorial is taken from and/or inspired by the GeoDjango basic apps project. [2]

Setting Up

Create a Spatial Database

First, a spatial database needs to be created for our project. If using PostgreSQL and PostGIS, then the following commands will create the database from a spatial database template:

$ sudo su - postgres
$ createdb -T template_postgis -O geo geodjango
$ exit

Here are becoming the postgres database super user to issue the createdb command because creating a PostGIS database requires elevated privileges. The -O geo option tells createdb to make the geo user the owner of the database. Replace geo with the username of your choice.

MySQL and Oracle

Create GeoDjango Project

Use the django-admin.py script like normal to create a geodjango project:

$ django-admin.py startproject geodjango

With the project initialized, now create a world Django application within the geodjango project:

$ cd geodjango
$ python manage.py startapp world

Configure settings.py

The geodjango project settings are stored in the settings.py file. Edit the database connection settings appropriately:

DATABASE_ENGINE = 'postgresql_psycopg2'
DATABASE_NAME = 'geodjango'
DATABASE_USER = 'geo'

In addition, modify the INSTALLED_APPS setting to include django.contrib.admin, django.contrib.gis, and geodjango.world (our newly created application):

INSTALLED_APPS = (
    'django.contrib.auth',
    'django.contrib.contenttypes',
    'django.contrib.sessions',
    'django.contrib.sites',
    'django.contrib.admin',
    'django.contrib.gis',
    'geodjango.world'
)

Geographic Data

World Borders

The world borders data is available in this zip file. Create a data directory in the world application, download the world borders data, and unzip:

$ mkdir world/data
$ cd world/data
$ wget http://thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip
$ unzip TM_WORLD_BORDERS-0.3.zip
$ cd ../..

The world borders ZIP file contains a set of data files collectively known as an ESRI Shapefile, one of the most popular geospatial data formats. When unzipped the world borders dataset includes files with the following extensions:

• .shp: Holds the vector data for the world borders geometries.
• .shx: Spatial index file for geometries stored in the .shp.
• .dbf: Database file for holding non-geometric attribute data (e.g., integer and character fields).
• .prj: Contains the spatial reference information for the geographic data stored in the shapefile.

Use ogrinfo to examine spatial data

The GDAL ogrinfo utility is excellent for examining metadata about shapefiles (or other vector data sources):

$ ogrinfo world/data/TM_WORLD_BORDERS-0.3.shp
INFO: Open of `world/data/TM_WORLD_BORDERS-0.3.shp'
      using driver `ESRI Shapefile' successful.
1: TM_WORLD_BORDERS-0.3 (Polygon)

Here ogrinfo is telling us that the shapefile has one layer, and that layer contains polygon data. To find out more we’ll specify the layer name and use the -so option to get only important summary information:

$ ogrinfo -so world/data/TM_WORLD_BORDERS-0.3.shp TM_WORLD_BORDERS-0.3
INFO: Open of `world/data/TM_WORLD_BORDERS-0.3.shp'
      using driver `ESRI Shapefile' successful.

Layer name: TM_WORLD_BORDERS-0.3
Geometry: Polygon
Feature Count: 246
Extent: (-180.000000, -90.000000) - (180.000000, 83.623596)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_1984",6378137.0,298.257223563]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]
FIPS: String (2.0)
ISO2: String (2.0)
ISO3: String (3.0)
UN: Integer (3.0)
NAME: String (50.0)
AREA: Integer (7.0)
POP2005: Integer (10.0)
REGION: Integer (3.0)
SUBREGION: Integer (3.0)
LON: Real (8.3)
LAT: Real (7.3)

This detailed summary information tells us the number of features in the layer (246), the geographical extent, the spatial reference system (“SRS WKT”), as well as detailed information for each attribute field. For example, FIPS: String (2.0) indicates that there’s a FIPS character field with a maximum length of 2; similarly, LON: Real (8.3) is a floating-point field that holds a maximum of 8 digits up to three decimal places. Although this information may be found right on the world borders website, this shows you how to determine this information yourself when such metadata is not provided.

Geographic Models

Defining a Geographic Model

Now that we’ve examined our world borders data set using ogrinfo, we can create a GeoDjango model to represent this data:

from django.contrib.gis.db import models

class WorldBorders(models.Model):
    # Regular Django fields corresponding to the attributes in the
    # world borders shapefile.
    name = models.CharField(max_length=50)
    area = models.IntegerField()
    pop2005 = models.IntegerField('Population 2005')
    fips = models.CharField('FIPS Code', max_length=2)
    iso2 = models.CharField('2 Digit ISO', max_length=2)
    iso3 = models.CharField('3 Digit ISO', max_length=3)
    un = models.IntegerField('United Nations Code')
    region = models.IntegerField('Region Code')
    subregion = models.IntegerField('Sub-Region Code')
    lon = models.FloatField()
    lat = models.FloatField()

    # GeoDjango-specific: a geometry field (MultiPolygonField), and
    # overriding the default manager with a GeoManager instance.
    mpoly = models.MultiPolygonField()
    objects = models.GeoManager()

    # So the model is pluralized correctly in the admin.
    class Meta:
        verbose_name_plural = "World Borders"

    # Returns the string representation of the model.
    def __unicode__(self):
        return self.name

Two important things to note:

1. The models module is imported from django.contrib.gis.db.
2. The model overrides its default manager with GeoManager; this is required to perform spatial queries.

When declaring a geometry field on your model the default spatial reference system is WGS84 (meaning the SRID is 4326) – in other words, the field coordinates are in longitude/latitude pairs in units of degrees. If you want the coordinate system to be different, then SRID of the geometry field may be customized by setting the srid with an integer corresponding to the coordinate system of your choice.

Run syncdb

After you’ve defined your model, it needs to be synched with the spatial database. First, let’s look at the SQL that will generate the table for the WorldBorders model:

$ python manage.py sqlall world

This management command should produce the following output:

BEGIN;
CREATE TABLE "world_worldborders" (
    "id" serial NOT NULL PRIMARY KEY,
    "name" varchar(50) NOT NULL,
    "area" integer NOT NULL,
    "pop2005" integer NOT NULL,
    "fips" varchar(2) NOT NULL,
    "iso2" varchar(2) NOT NULL,
    "iso3" varchar(3) NOT NULL,
    "un" integer NOT NULL,
    "region" integer NOT NULL,
    "subregion" integer NOT NULL,
    "lon" double precision NOT NULL,
    "lat" double precision NOT NULL
)
;
SELECT AddGeometryColumn('world_worldborders', 'mpoly', 4326, 'MULTIPOLYGON', 2);
ALTER TABLE "world_worldborders" ALTER "mpoly" SET NOT NULL;
CREATE INDEX "world_worldborders_mpoly_id" ON "world_worldborders" USING GIST ( "mpoly" GIST_GEOMETRY_OPS );
COMMIT;

If satisfied, you may then create this table in the database by running the syncdb management command:

$ python manage.py syncdb
Creating table world_worldborders
Installing custom SQL for world.WorldBorders model

Importing Spatial Data

This section will show you how to take the data from the world borders shapefile and import it into GeoDjango models using the LayerMapping utility. There are many different different ways to import data in to a spatial database – besides the tools included within GeoDjango, you may also use the following to populate your spatial database:

• ogr2ogr: Command-line utility, included with GDAL, that supports loading a multitude of vector data formats into the PostGIS, MySQL, and Oracle spatial databases.
• shp2pgsql: This utility is included with PostGIS and only supports ESRI shapefiles.

GDAL Interface

Earlier we used the the ogrinfo to explore the contents fo the world borders shapefile. Included within GeoDjango is an interface to GDAL’s powerful OGR library – in other words, you’ll be able explore all the vector data sources that OGR supports via a Pythonic API.

First, invoke the Django shell:

$ python manage.py shell

If the World Borders data was downloaded like earlier in the tutorial, then we can determine the path using Python’s built-in os module:

>>> import os
>>> from geodjango import world
>>> world_shp = os.path.abspath(os.path.join(os.path.dirname(world.__file__),
...                             'data/TM_WORLD_BORDERS-0.3.shp'))

Now, the world borders shapefile may be opened using GeoDjango’s DataSource interface:

>>> from django.contrib.gis.gdal import *
>>> ds = DataSource(world_shp)
>>> print ds
/ ... /geodjango/world/data/TM_WORLD_BORDERS-0.3.shp (ESRI Shapefile)

Data source objects can have different layers of geospatial features; however, shapefiles are only allowed to have one layer:

>>> print len(ds)
1
>>> lyr = ds[0]
>>> print lyr
TM_WORLD_BORDERS-0.3

You can see what the geometry type of the layer is and how many features it contains:

>>> print lyr.geom_type
Polygon
>>> print len(lyr)
246

Note

Unfortunately the shapefile data format does not allow for greater specificity with regards to geometry types. This shapefile, like many others, actually includes MultiPolygon geometries in its features. You need to watch out for this when creating your models as a GeoDjango PolygonField will not accept a MultiPolygon type geometry – thus a MultiPolygonField is used in our model’s definition instead.

The Layer may also have a spatial reference system associated with it – if it does, the srs attribute will return a SpatialReference object:

>>> srs = lyr.srs
>>> print srs
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_1984",6378137.0,298.257223563]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]
>>> srs.proj4 # PROJ.4 representation
'+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs '

Here we’ve noticed

In addition, shapefiles also support attribute fields that may contain additional data. Here are the fields on the World Borders layer:

>>> print lyr.fields
['FIPS', 'ISO2', 'ISO3', 'UN', 'NAME', 'AREA', 'POP2005', 'REGION', 'SUBREGION', 'LON', 'LAT']

Here we are examining the OGR types (e.g., whether a field is an integer or a string) associated with each of the fields:

>>> [fld.__name__ for fld in lyr.field_types]
['OFTString', 'OFTString', 'OFTString', 'OFTInteger', 'OFTString', 'OFTInteger', 'OFTInteger', 'OFTInteger', 'OFTInteger', 'OFTReal', 'OFTReal']

You can iterate over each feature in the layer and extract information from both the feature’s geometry (accessed via the geom attribute) as well as the feature’s attribute fields (whose values are accessed via get()):

>>> for feat in lyr:
...    print feat.get('NAME'), feat.geom.num_points
...
Guernsey 18
Jersey 26
South Georgia South Sandwich Islands 338
Taiwan 363

Layer objects may be sliced:

>>> lyr[0:2]
[<django.contrib.gis.gdal.feature.Feature object at 0x2f47690>, <django.contrib.gis.gdal.feature.Feature object at 0x2f47650>]

And individual features may be retrieved by their feature ID:

>>> feat = lyr[234]
>>> print feat.get('NAME')
San Marino

Here the boundary geometry for San Marino is extracted and looking exported to WKT and GeoJSON:

>>> geom = feat.geom
>>> print geom.wkt
POLYGON ((12.415798 43.957954,12.450554 ...
>>> print geom.json
{ "type": "Polygon", "coordinates": [ [ [ 12.415798, 43.957954 ], [ 12.450554, 43.979721 ], ...

LayerMapping

We’re going to dive right in – create a file called load.py inside the world application, and insert the following:

import os

world_mapping = {
    'fips' : 'FIPS',
    'iso2' : 'ISO2',
    'iso3' : 'ISO3',
    'un' : 'UN',
    'name' : 'NAME',
    'area' : 'AREA',
    'pop2005' : 'POP2005',
    'region' : 'REGION',
    'subregion' : 'SUBREGION',
    'lon' : 'LON',
    'lat' : 'LAT',
    'mpoly' : 'MULTIPOLYGON',
}

world_shp = os.path.abspath(os.path.join(os.path.dirname(__file__), 'data/TM_WORLD_BORDERS-0.3.shp'))

def run(verbose=True):
    from django.contrib.gis.utils import LayerMapping
    from models import WorldBorders

    lm = LayerMapping(WorldBorders, world_shp, world_mapping,
                      transform=False, encoding='iso-8859-1')

    lm.save(strict=True, verbose=verbose)

Try ogrinspect

Spatial Queries

Spatial Lookups

Lazy Geometries

GeoQuerySet Methods

Putting your data on the map

Google

Geographic Admin

Footnotes

[1]	Special thanks to Bjørn Sandvik of thematicmapping.org for providing and maintaining this dataset.

[2]	GeoDjango basic apps was written by Dane Springmeyer, Josh Livni, and Christopher Schmidt.