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