Demo 1: GeoDjango Exploration

Django has nice geospatial support built-in. This is particularly good when used with a Postgres db. This page presents an exploration of basic retrieval and manipulation of geospatial data. To ease presentation the folio package is used which allows maps to be configured in the backend and displayed without the need for any javascript.

The folium package can be found here: https://python-visualization.github.io/folium/latest/index.html.
The install instructions for GeoDjango are available here. A docker image for postgres/postgis is available here.

I've found a setup with Docker to be the best approach. I use a standard python docker image (slim-buster, which uses Debian) to which I add in the necesssary additional geospatial libs. This together with the postgis image works pretty much straight out of the box and is pretty simple if you are already familliar with Docker.

GeoDjango Documentation Links

Not quite sure why, but I found navigating the documentation a little painful (wierd as Django docs are fab). Anyway, here are the main pages for quick reference.

Main Contents Page
Queryset API Ref (including spatial queries)
Database functions - Area, distance, centroid, i/o ops
GEOS API ref


Features can be made using the django (python-wrapped GEOS classes) such as Point, Line, Polygon.

Below a point object is created by passing in the coordinates and the spatial reference. This point object is then passed into the geom field when creating a new venue object.

    
    from django.contrib.gis.geos import Point

    coords = [-1.568, 53.511]
    point_obj = Point(coords, srid=4326)

    new_obj = Venue.objects.create(
        name="A New Venue", 
        lat=coords[0], 
        lon=coords[1], 
        geom=point_obj)

       

Points can also be created using well-known text. In this instance, creating a well-known-text string from the coordinates using an f-string interpolation.

    
    from django.contrib.gis.geos import GEOSGeometry, Point

    coords = [-1.568, 53.511]
    point_wkt = f"POINT({coords[0]} {coords[1]})"   
   
    point_obj = GEOSGeometry(point_wkt)           
    point_coordinates = Point(point_obj.coords)

    new_obj = Venue.objects.create(
        name="A New Venue", 
        lat=coords[0], 
        lon=coords[1], 
        geom=point_obj)

       

Transforms and Functions

In this example the city boundary is in a different spatial reference (EPSG 28992) to the default map (EPSG 4326 / WGS 84). The queryset is annotated with a transformed geometry using the Transform() function.

The centroid of the boundary polygon is also calucated and annotated using the Centroid() function.

    
    en_boundary = (
        EnschedeBuiltBoundary.objects
        .annotate(transformed=Transform("geom", srid=4326))
        .annotate(centroid=Centroid("transformed"))
        .first()
    )

    # passed in to get area and length 
    # note that using originaal projected srs.
    en_boundary.geom.area
    en_boundary.geom.length
       

Serialization

The output from a GeoDjango query will be a single object (.get) or queryset of objects (.filter). The geometry is an GEOSGeometry object which can be manipulated. These can be output as geojson or wkt. See the documentation for examples here.

Where multiple features are returned, these can be combined into a single GeoJson by serialising the queryset. Discussed in the documentation here.

    
    from django.core.serializers import serialize
    import folium

    query = CountryBorder.objects.filter(
        name__in=[
        "United Kingdom", 
        "France", 
        "Germany",
        ])

    serialized_data = serialize(
        "geojson", 
        query, 
        geometry_field="geom", 
        fields=["id", "pop2005", "name",
        ])
    
    m = folium.Map(location=(53, -1.43), zoom_start=2)
    folium.GeoJson(serialized_data).add_to(m)
    context["map"] = m._repr_html_()

       

Geospatial Queries

Geospatial filteringis available in geodjango in a similar vein to postgis. Below an example intersect of point with a polygon (country).

NOTE: that the definition of the coordinates used by folium are X, Y cartesian whereas postgis and osgeom, gdal etc., depending on the projection, are returned as lon/lat. folium.GeoJson accepts lon/lat.

A few hours of my life I'll never get back

    
    from django.contrib.gis.geos import GEOSGeometry, Point

    pnt = Point([-1.48, 53.369], srid=4326)
    country = CountryBorder.objects.filter(geom__intersects=pnt)
       

Spatial Query - Touches

In this example the "touches" spatial query is used to identify adjoining buildings. This is the equivalent of PostGIS ST_Touches.

The result of the spatial query can be found as a layer.

Note: If you zoom in you'll notice that the data does not align with underlying basemap. In some instances the issue is related to buildings having changed (the dataset is over 12 years old), whereas in others it would appear to be a projection/data issue. I have tried in QGIS and ArcGIS-pro to reproject from the original RD New(28992)/Amersfort to WGS84 but none of the available transforms appear to work. Requires further investigation.

    

   # Select a building from id (standard query)
   target_building = (
        EnschedeBuildings.objects
        .filter(id=22371)
        .annotate(transformed=Transform("geom", srid=4326))
        .first()
    )
    
    # Get buildings touching target
    touching_buildings = EnschedeBuildings.objects
                            .filter(geom__touches=target_building.geom)
    
    #serialize the response
    touching_buildings_serialized = serialize(
        "geojson",
        touching_buildings,
        geometry_field="geom",
        id_field="id",
    )
       

Spatial Query : Distance Within

In this example a distance witin spatial query is used to identify all the buildings within 20m of the railway. As the railway is multi-line Union of the features is taken to pass into the spatial query. dwithin is the equivalent of postGIS ST_DWithin().



    # Get the boundary, transform for display and centroid to center the map.
    en_boundary = (
        EnschedeBuiltBoundary.objects
        .annotate(transformed=Transform("geom", srid=4326))
        .annotate(centroid=Centroid("transformed"))
        .first()
    )

    # get a Union of the railways queryset
    railways = EnschedeRailwayLines.objects.aggregate(Union("geom"))

    # dwithin spatial query betweeen railway geometry and the buildings
    buildings = EnschedeBuildings.objects
        .filter(geom__dwithin=(railways["geom__union"], D(m=30)))

    # Serialise the returned buildings for display in folium.GeoJson()
    serialized_buildings = serialize(
        "geojson",
        buildings,
        geometry_field="geom",
        fields=["id"], 
    )

   

Spatial Function: Distance

In this example we take two points (home centroid and the central station). We then calculate the distance between the two using the Geo Distance function (not be confused with a Distance object!).

Separately a simple line is created (in folium) for presentation on the map.

    

    # Get home and its centroid
    target_building = (
        EnschedeBuildings.objects
        .filter(id=22371)
        .annotate(transformed=Transform("geom", srid=4326))
        .first()
    )
    home_centroid = target_building.transformed.centroid

    # Get all train stations and serialize them
    train_stns = EnschedeRailwayStations.objects.all()
                    .annotate(transformed=Transform("geom", srid=4326))

    ts_serialized = serialize(
        "geojson",
        train_stns,
        geometry_field="geom",
        id_field="id",
    )
    
    # get the main trainstation and annotate with distance to home.
    main_station = (
        EnschedeRailwayStations.objects
        .filter(stationnam="Enschede")
        .annotate(transformed=Transform("geom", srid=4326))
        .annotate(dist=Distance("transformed", home_centroid))
        .first()
    )