Nuventure Blog

Knowledge.transmit!

Python

GeoDjango

Steve

GeoDjango is an included contrib module to make easy GIS (Geographic Information System) web apps with location-based services. GeoDjango provides a toolbox of utilities for building GIS web applications and also bindings to popular spatial libraries such as GEOS, GDAL, and GeoIP, which can be used separately without Django in any Python application or interactively in the shell.

Let’s get started by installing Python 3.

$ python3 --version

Now let’s install GeoDjango dependencies by running the following commands.

$ sudo aptitude install gdal-bin libgdal-dev 
$ sudo aptitude install python3-gdal

Now, we have to create spatial databases. PostGIS 2 ships with PostgreSQL as an extension for spatial functionality. Geodjango also supports Sqlite (via Spatialite), MySql and Oracle backends.

$ createdb  <db name> 
$ psql <db name> 
> CREATE EXTENSION postgis;

Let’s now get started with our project. Let’s call our project geodjango by running:

$ django-admin startproject geodjango

This will initialize the project. Now, we can create a django project within our project.

$ cd geodjango 
$ python manage.py startapp world

Let’s edit the database connection in geodjango/settings.py to match our setup.

DATABASES = {
     'default': {
          'ENGINE': 'django.contrib.gis.db.backends.postgis',
          'NAME': 'geodjango',
          'USER': 'geo',
     },
 }

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

INSTALLED_APPS = [
     'django.contrib.admin',
     'django.contrib.auth',
     'django.contrib.contenttypes',
     'django.contrib.sessions',
     'django.contrib.messages',
     'django.contrib.staticfiles',
     'django.contrib.gis',
     'world',
 ]

Geographic Data

World Borders

We can avail the data of the world borders from Here. Let’s create a data directory in the world application, download the world borders data, and unzip. Use the following commands:

$ mkdir world/data
$ cd world/data
$ wget https://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 which is one of the most popular geospatial data formats known as an ESRI Shapefile. 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.

The GDAL ogrinfo utility allows examining the metadata of 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)

ogrinfo tells us that the shapefile has one layer and that this layer contains polygon data. Let’s specify the layer name and use the -so option to get the 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, the geographic bounds of the data, the spatial reference system (“SRS WKT”), as well as type information for each attribute field. For example, FIPS: String (2.0) indicates that the FIPS character field has 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.

Geographic Models

Defining a Geographic Model

Let’s create a GeoDjango model to represent the data:

from django.contrib.gis.db import models

class WorldBorder(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)
    mpoly = models.MultiPolygonField()

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

The default spatial reference system for geometry fields is WGS84 . In other words, the field coordinates are in longitude, latitude pairs in units of degrees. To use a different coordinate system, let’s set the SRID of the geometry field with the srid argument. Use an integer representing the coordinate system’s EPSG code.

Run migrate

After defining your model, we need to sync it with the database. First, we need to create a database migration. Run the following command:

$ python manage.py makemigrations
Migrations for 'world':
  world/migrations/0001_initial.py:
    - Create model WorldBorder

You may inspect the raw code generated by the above migration:

$ python manage.py sqlmigrate world 0001

The output:

BEGIN;
--
-- Create model WorldBorder
--
CREATE TABLE "world_worldborder" (
    "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
    "mpoly" geometry(MULTIPOLYGON,4326) NOT NULL
)
;
CREATE INDEX "world_worldborder_mpoly_id" ON "world_worldborder" USING GIST ( "mpoly" );
COMMIT;

Now, let’s run migrate to execute the sql:

$ python manage.py migrate
Operations to perform:
  Apply all migrations: admin, auth, contenttypes, sessions, world
Running migrations:
  ...
  Applying world.0001_initial... OK

GDAL Interface

GeoDjango also includes a Pythonic interface to GDAL’s powerful OGR library that can work with all the vector data sources that OGR supports.

First, let’s invoke the Django shell:

$ python manage.py shell

We can determine the path of World Border data using Python’s built-in os module:

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

Now, open the world borders shapefile using GeoDjango’s DataSource interface:

>>> from django.contrib.gis.gdal import DataSource
>>> 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

We can see the layer’s geometry type the features it contains:

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

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 '

This shapefile is in the popular WGS84 spatial reference system.

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']

The following code will let you examine the OGR types (e.g. integer or 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']

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

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

Boundary geometries may be exported as WKT and GeoJSON. Both are plaintext representations of shapes which are originally encoded in binary.

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

To import the data, use a LayerMapping in a Python script. Create a file called load.py inside the world application, with the following code:

import os
from django.contrib.gis.utils import LayerMapping
from .models import WorldBorder

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):
    lm = LayerMapping(WorldBorder, world_shp, world_mapping, transform=False)
    lm.save(strict=True, verbose=verbose)

Let’s invoke the Django shell from the geodjango project directory:

$ python manage.py shell

Next, let’s import the load module, call the run routine:

>>> from world import load
>>> load.run()

grinspect

The ogrinspect command introspects a GDAL-supported vector data source (e.g., a shapefile) and generates a model definition and LayerMapping dictionary automatically.

The general usage of the command goes as follows:

$ python manage.py ogrinspect [options] <data_source> <model_name> [options]

data_source is the path to the GDAL-supported data source and model_name is the name to use for the model. Command-line options may be used to further define how the model is generated.

For example, the following command nearly reproduces the WorldBorder model and mapping dictionary created above, automatically:

$ python manage.py ogrinspect world/data/TM_WORLD_BORDERS-0.3.shp WorldBorder \
    --srid=4326 --mapping --multi

The command produces the following output, which may be copied directly into the models.py of a GeoDjango application:

# This is an auto-generated Django model module created by ogrinspect.
from django.contrib.gis.db import models

class WorldBorder(models.Model):
    fips = models.CharField(max_length=2)
    iso2 = models.CharField(max_length=2)
    iso3 = models.CharField(max_length=3)
    un = models.IntegerField()
    name = models.CharField(max_length=50)
    area = models.IntegerField()
    pop2005 = models.IntegerField()
    region = models.IntegerField()
    subregion = models.IntegerField()
    lon = models.FloatField()
    lat = models.FloatField()
    geom = models.MultiPolygonField(srid=4326)

# Auto-generated `LayerMapping` dictionary for WorldBorder model
worldborders_mapping = {
    'fips' : 'FIPS',
    'iso2' : 'ISO2',
    'iso3' : 'ISO3',
    'un' : 'UN',
    'name' : 'NAME',
    'area' : 'AREA',
    'pop2005' : 'POP2005',
    'region' : 'REGION',
    'subregion' : 'SUBREGION',
    'lon' : 'LON',
    'lat' : 'LAT',
    'geom' : 'MULTIPOLYGON',
}

Spatial Queries

Spatial Lookups

GeoDjango adds spatial lookups to the Django ORM. We can find the country in the WorldBorder table that a point belongs to. Let’s fire up the management shell:

$ python manage.py shell

Define the point of interest:

>>> pnt_wkt = 'POINT(-95.3385 29.7245)'

The pnt_wkt string represents the point at -95.3385 degrees longitude, 29.7245 degrees latitude. The geometry is in a format known as Well Known Text (WKT), a standard issued by the Open Geospatial Consortium (OGC). [ Import the WorldBorder model, and perform a contains lookup using the pnt_wkt as the parameter:

>>> from world.models import WorldBorder
>>> WorldBorder.objects.filter(mpoly__contains=pnt_wkt)
<QuerySet [<WorldBorder: United States>]>

We can also use a GEOS geometry object. Here, we can combine the intersects spatial lookup with the get method to retrieve only the WorldBorder instance for San Marino instead of a queryset:

>>> from django.contrib.gis.geos import Point
>>> pnt = Point(12.4604, 43.9420)
>>> WorldBorder.objects.get(mpoly__intersects=pnt)
<WorldBorder: San Marino>

Automatic Spatial Transformations

When doing spatial queries, GeoDjango automatically transforms geometries if they’re in a different coordinate system.

>>> from django.contrib.gis.geos import GEOSGeometry, Point
>>> pnt = Point(954158.1, 4215137.1, srid=32140)

pnt may also be constructed with EWKT, an “extended” form of WKT that includes the SRID:

>>> pnt = GEOSGeometry('SRID=32140;POINT(954158.1 4215137.1)')

GeoDjango’s ORM will automatically wrap geometry values in transformation SQL, allowing the developer to work at a higher level of abstraction:

>>> qs = WorldBorder.objects.filter(mpoly__intersects=pnt)
>>> print(qs.query) # Generating the SQL
SELECT "world_worldborder"."id", "world_worldborder"."name", "world_worldborder"."area",
"world_worldborder"."pop2005", "world_worldborder"."fips", "world_worldborder"."iso2",
"world_worldborder"."iso3", "world_worldborder"."un", "world_worldborder"."region",
"world_worldborder"."subregion", "world_worldborder"."lon", "world_worldborder"."lat",
"world_worldborder"."mpoly" FROM "world_worldborder"
WHERE ST_Intersects("world_worldborder"."mpoly", ST_Transform(%s, 4326))
>>> qs # printing evaluates the queryset
<QuerySet [<WorldBorder: United States>]>

Lazy Geometries

GeoDjango loads geometries in a standardized textual representation. When the geometry field is first accessed, GeoDjango creates a GEOSGeometry object, exposing powerful functionality, such as serialization properties for popular geospatial formats:

>>> sm = WorldBorder.objects.get(name='San Marino')
>>> sm.mpoly
<MultiPolygon object at 0x24c6798>
>>> sm.mpoly.wkt # WKT
MULTIPOLYGON (((12.4157980000000006 43.9579540000000009, 12.4505540000000003 43.9797209999999978, ...
>>> sm.mpoly.wkb # WKB (as Python binary buffer)
<read-only buffer for 0x1fe2c70, size -1, offset 0 at 0x2564c40>
>>> sm.mpoly.geojson # GeoJSON
'{ "type": "MultiPolygon", "coordinates": [ [ [ [ 12.415798, 43.957954 ], [ 12.450554, 43.979721 ], ...

This includes access to all of the advanced geometric operations provided by the GEOS library:

>>> pnt = Point(12.4604, 43.9420)
>>> sm.mpoly.contains(pnt)
True
>>> pnt.contains(sm.mpoly)
False

Putting your data on the map

Geographic Admin

GeoDjango extends Django’s admin application with support for editing geometry fields.

Basics

GeoDjango also supplements the Django admin by allowing users to create and modify geometries on a JavaScript slippy map.

Let’s create a file called admin.py inside the world application with the following code:

from django.contrib.gis import admin
from .models import WorldBorder

admin.site.register(WorldBorder, admin.GeoModelAdmin)

Next, edit the urls.py in the geodjango application folder as follows:

from django.contrib.gis import admin
from django.urls import include, path

urlpatterns = [
    path('admin/', admin.site.urls),
]

Create an admin user:

$ python manage.py createsuperuser

Next, start up the Django development server:

$ python manage.py runserver

Browse in to http://localhost:8000/admin/, and log in with the user you just created. Browse to any of the WorldBorder entries. The borders may be edited by clicking on a polygon and dragging the vertices to the desired position.

GeoDjango turns Django into an effective geographic Web framework. GeoDjango eases the effort of making GIS (Geographic information system) web apps with location-based services.

Related Posts

Begin typing your search term above and press enter to search. Press ESC to cancel.

Back To Top