Generate regular hex and square grids covering area of interest

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

# the projected CRS
crs_vic = {'init': 'epsg:28355'}

# maybe try this one: epsg:3395

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

# north east most point of melb
p_upper = [146.742224, -36.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)

# 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 = 1000
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)

# 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_hex = gpd.GeoDataFrame(geometry=hexes)
gdf_hex = gdf_hex.to_crs(crs_abs)

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

# 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] + a*m*r, coords0[1] + b*n*r]).reshape(1,-1)
        polys.append(Polygon(d+poly_square*r))
        
squares = gpd.GeoSeries(polys, crs=crs_vic)

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

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

# Greater Melbourne
gdf_gm = gdf_vic[gdf_vic['GCC_NAME16']=='Greater Melbourne']

# The SA4's that make up Greater Melbourne and Geelong
sa4_aoi = ['Melbourne - Inner', 'Melbourne - Inner East',
       'Melbourne - Inner South', 'Melbourne - North East',
       'Melbourne - North West', 'Melbourne - Outer East',
       'Melbourne - South East', 'Melbourne - West',
       'Mornington Peninsula']

gdf_aoi = gdf_vic[gdf_vic['SA4_NAME16'].isin(sa4_aoi)]

# join hexagon grid with AOI
gdf_new = gpd.sjoin(gdf_aoi, gdf_hex, op='intersects', how='inner')

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

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

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

# just get the hexagons that cover melbourne
sqr_gm = gdf_sqr[gdf_sqr['SquareLabel'].isin(gdf_new['SquareLabel'].values)]

# write to a shapefile
sqr_gm.to_file('sqrgrid_aoi_1km.shp', index=False)

# 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

# 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(sqr_gm.geometry).add_to(m)

# display all of this
m