Hosted with nbsanity. See source notebook on GitHub.

Geopandas, Shapely, Folium learning: Generate a regular hexagonal grid to cover a region

27 May 2019

As part of one of my work projects, I’ve been trying to create a regular grid of hexagons that tile the Greater Melbourne region. This is my attempt at this problem. It was also a nice way for me to start learning some geospatial libraries in Python: GeoPandas (spatial joins, filtering, co-ordinate reference systems); Shapely; and Folium.

The basic approach is: - Define a regular hexagon - Define a flat CRS and a projected CRS. So x, y and lat, lng values respectively. - Write the code the create a Shapely Polygon that defines a regular hexagon with a particular radius in metres. - Define an area (in lat, lng values) over which we would like to create a tiling. - Intersect this area with the geographic area of interest, in this case the Greater Melbourne region.

1.0. Setup

Import all the necessary libraries and then define the Coordinate Reference Systems (CRS) that we will be working with.

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np

import folium

from shapely.geometry import Polygon, Point, MultiPolygon

%matplotlib inline

1.1. Coordinate Reference Systems

# The CRS that the ABS SA2 shapefile uses
crs_abs = {'init': 'epsg:4283'}

# the 'flattened' CRS used to define the grid before projecting it.
crs_vic = {'init': 'epsg:28355'}

# south west most point of melb
p_lower = [143.338683, -39.185049]

# north east most point of melb
p_upper = [146.742224, -37.559378]

# put these in a Geoseries
gs1 = gpd.GeoSeries([Point(p_lower), Point(p_upper)], crs=crs_abs)

# convert each of these points to the 'flat coordinate system'
gs1b = gs1.to_crs(crs_vic)

2.0. Define the regular hexagonal grid

2.1. Define a regular hexagon in the flat geometry

Firstly we define a regular hexagon and then tile this over the given rectangle

# startpoint coords
coords0 = gs1b[0].x, gs1b[0].y

# the x and y extent of the grid
Xdiff = gs1b[1].x - gs1b[0].x
Ydiff = gs1b[1].y - gs1b[0].y

# the cell width and height (in metres)
celllength = 2000
r = celllength

# fundamental hexagon
a, b = np.sqrt(3) / 2, 1 / 2
poly_hex = np.array([(0, r), (-r*a, r*b), (-r*a, -r*b), (0, -r), (r*a, -r*b), (r*a, r*b)])

# calculate the number of cells in the X and Y directions
num_cells_x = int(Xdiff / celllength / 2)
num_cells_y = int(Ydiff / celllength / 2)

2.2. Tile this hexagon over the rectangle covering the specified area.

Clean this up to void using the for loops.

# centres of a regular grid of hexagons of unit radius
polys = []

# using a for loop
for m in range(0, num_cells_x):
    for n in range(0, num_cells_y):
        # generate the two sublattices
        d = np.array([coords0[0] + 2*m*a*r, coords0[1] + 3*2*n*b*r]).reshape(1,-1) # displacement vector
        polys.append(Polygon(d+poly_hex))
        
        d = np.array([coords0[0] + (2*m+1)*a*r, coords0[1] + 3*(2*n+1)*b*r]).reshape(1,-1)
        polys.append(Polygon(d+poly_hex))
        
hexes = gpd.GeoSeries(polys, crs=crs_vic)

# convert to dataframe and change crs
gdf = gpd.GeoDataFrame(geometry=hexes)
gdf = gdf.to_crs(crs_abs)

# add a label for each hexagon for when we do a spatial join
gdf['HexLabel'] = np.arange(gdf.shape[0])

3.0. Load up the Vic SA2 values and shapefiles

The SA2 files are here: https://www.abs.gov.au/AUSSTATS/[email protected]/DetailsPage/1270.0.55.003July%202016?OpenDocument

# load up the shape files for Vic and filter down to Greater Melbourne ones
gdf_vic = gpd.read_file('data/shapefiles/SA2_2016_AUST.shp')

#sa2_gm = gdf_vic[gdf_vic['GCC_NAME16']=='Greater Melbourne'].geometry
gdf_vic = gdf_vic[gdf_vic['geometry'].isnull() == False] # remove the null values
gdf_gm = gdf_vic[gdf_vic['GCC_NAME16']=='Greater Melbourne'] # only interested in Greater Melbourne

# join hexagon grid with Greater Melbourne
gdf_new = gpd.sjoin(gdf_gm, gdf, op='intersects', how='inner')

# just get the hexagons that cover melbourne
hex_gm = gdf[gdf['HexLabel'].isin(gdf_new['HexLabel'].values)]

# write to a shapefile
hex_gm.to_file('hexgrid_1km.shp', index=False)

4.0. Plot the regular hexgrid covering Greater Melbourne

# initialize the map
m = folium.Map(location=[-37.8, 145], zoom_start=8.5, tiles='cartodbpositron')

# add the hexagon grid to the map
folium.GeoJson(hex_gm.geometry).add_to(m)

# display all of this
m