def parameter_estimate(gdf1, gdf2, area_field, param_field, label2_field):
    
    # calculate the overlay
    gdf_overlay = gpd.overlay(gdf1, gdf2, how='intersection')
    
    # calculate the area of each of the new polygons
    gdf_overlay['area_new'] = gdf_overlay['geometry'].apply(lambda p: p.area)
    
    # calculate the fraction of the new polygons intersecting with the old ones
    gdf_overlay['fraction_overlap'] = gdf_overlay['area_new'] / gdf_overlay[area_field]
    
    # estimate the value of the parameter for each of the new polygons in their overlaps
    gdf_overlay[param_field + '_new'] = gdf_overlay['fraction_overlap'] * gdf_overlay[param_field]
    
    # calculate the weighted sum to estimate the new value for each of the new polygons
    df_ws = gdf_overlay.groupby(label2_field)[param_field + '_new'].sum().reset_index()
    
    return df_ws

# define a CRS
crs_abs = {'init': 'epsg:4283'}

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

# load up the Vic meshblock codes and corresponding dwelling and person counts
df_mb_vic = pd.read_csv('data/refined/MB_VIC_2016.csv')

# let's have a look at this
df_mb_vic.head(5)

gdf_1km_sqr = gpd.read_file('data/refined/sqrgrid_aoi_1km.shp')
gdf_1km_sqr.head(5)

gdf_1km_hex = gpd.read_file('data/refined/hexgrid_aoi_1km.shp')
gdf_1km_hex.head(5)

# load up the meshblocks shapefiles
gdf_mb = gpd.read_file('data/opendata/Census/1270055001_mb_2016_vic_shape/MB_2016_VIC.shp')

# rename the column to be consistent with the other data sources
gdf_mb = gdf_mb.rename(columns={'MB_CODE16': 'MB_CODE_2016'})
gdf_mb = gdf_mb[['MB_CODE_2016', 'geometry']]

# convert the MB_CODE_2016 to a integer
gdf_mb['MB_CODE_2016'] = gdf_mb['MB_CODE_2016'].astype(np.int)

gdf_mb.head(5)

# look at just one suburb
#SA2_NAME_2016 = 'Carlton'
#SA3_NAME_2016 = 'Melbourne City'
#df_smaller = df_mb_vic[df_mb_vic['SA3_NAME_2016']==SA3_NAME_2016]
#df_smaller.head(5)

GCCSA_NAME_2016 = 'Greater Melbourne'
df_aoi = df_mb_vic[df_mb_vic['GCCSA_NAME_2016']==GCCSA_NAME_2016]
df_aoi.head(5)

df_aoi['Person'].sum()

# join the shapefiles with the Vic. meshblock data
gdf_mb_sa2 = df_aoi.merge(gdf_mb, on='MB_CODE_2016')

nice_cols = ['MB_CODE_2016', 'MB_CATEGORY_NAME_2016_y', 'SA1_MAINCODE_2016', 'SA2_NAME_2016', 'SA3_NAME_2016',
            'Dwelling', 'Person']

meshblocks_sa2 = gpd.GeoDataFrame(data=gdf_mb_sa2[nice_cols], geometry=gdf_mb_sa2['geometry'])
meshblocks_sa2.crs = crs_abs

gdf_sa2 = gpd.sjoin(gdf_1km_sqr, meshblocks_sa2, op='intersects')
sqr_sa2 = gdf_sa2.dissolve(by='SquareLabe', aggfunc='first')

gdf_sa2 = gpd.sjoin(gdf_1km_hex, meshblocks_sa2, op='intersects')
hex_sa2 = gdf_sa2.dissolve(by='HexLabel', aggfunc='first')

meshblocks_sa2_new = meshblocks_sa2.to_crs(crs_vic)
sqr_sa2_new = sqr_sa2.to_crs(crs_vic) # convert the squares
hex_sa2_new = hex_sa2.to_crs(crs_vic) # convert the hexagons

meshblocks_sa2_new.head(5)

squares = sqr_sa2_new[['MB_CODE_2016', 'geometry']].reset_index()
squares.head(5)

hexes = hex_sa2_new[['MB_CODE_2016', 'geometry']].reset_index()
hexes.head(5)

mesh = meshblocks_sa2_new[['MB_CODE_2016', 'Person', 'geometry']]

# calculate the area of each of the meshblocks in the projected CRS
mesh['area'] = mesh['geometry'].apply(lambda p: p.area)

mesh.head(5)

square_data = np.round(parameter_estimate(mesh, squares, 'area', 'Person', 'SquareLabe'))

square_data.head(5)

hex_data = np.round(parameter_estimate(mesh, hexes, 'area', 'Person', 'HexLabel'))

hex_data.head(5)

square_data.to_csv('data/refined/population_sqr_1km.csv', index=False)

hex_data.to_csv('data/refined/population_hex_1km.csv', index=False)

square_data.head(5)

hex_data.head(5)

square_data['Person_new'].sum()

hex_data['Person_new'].sum()