# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# # File Size Investigation
# Deleting items from bottom to top to try to get under 100 MB
#
# Item | After Removal | Size (kB)
# ---|---|---
# Original | 119424 | -
# Single Feature Choropleth | 111034 | 8390
# Final Cluster Maps (two rows) | 102255 | 8779
# P-value DF | 102227 | 28
# Feature DF Sorted | 102183 | 44
# Feature DF Unsorted | 102139 | 44
# Distance Matrix Image | 102125 | 14
# Dendrogram | 102110 | 15
# Distance Matrix Image | 102096 | 14
# Sorted Cluster Index Plot | 102068 | 28
# Distance Matrix Image | 102054 | 14
# Kopt P Values | 102051 | 3
# Scaled DF | 102007 | 44
# 3D Scatter Plot | 101610 | 397
# Pairplot 3-wide Figure | 101353 | 257
# GMM Separation Statistic | 101333 | 20
# Gap Statistic 3 Image | 101271 | 52
# Gap Statistic 2 Image | 101210 | 61
# Gap Statistic 1 Image | 101134 | 74
# Silhouette | 101102 | 32
# Silhouette | 101072 | 30
# Silhouette | 101049 | 23
# Gaussian Weights Grid 26 | 100870 | 179
# 3D Scatter Plot | 96883 | 3987
# Pairplot 3-wide Figure | 96649 | 234
# Gaussian Weights Grid 20 | 96511 | 138
# Pairplot 6-wide Image | 95042 | 1607
# Pairplot 6-wide Image | 93517 | 1525
# Gaussian Weights 30 | 93335 | 282
# Pairplot 11-wide Image | 91716 | 1619
# Pairplot 11-wide Image | 89995 | 1711
# Gaussian Weights 30 | 89828 | 167
# Pairplot 11-Wide Image One Color | 88001 |1827
# PCA Weights DF | 87972 | 29
# List - Pairplot K-Means 5 Features 11 Normal | |
# List - Pairplot K-Means 5 Features 11 Uniform | |
# List - Pairplot GMM 10/11 Uniform | |
# List - Pairplot GMM 8/11 Normal | |
# List - Pairplot Mean Shift 10/11 Uniform | |
# Single Feature Choropleth | 79665 | 8307
# K-Means Cluster Maps (one row) | 71354 | 8311
# Map: Census Area Difference | 62681 | 8673
# Area Difference Plots | 62625 | 56
# Map: Montreal Zero Pop | 62539 | 86
# Map: FSA Image | 61419 | 1120
# Map: FSA Interactive | 54289 | 7130 *to image*
# Map: Invalid FSAs | 39946 | 14343 *to single interactive map*
# Map: CT Image | 38816 | 1130
# Map: ADA Image | 37685 | 1131
# Map: DA Image | 36518 | 1167
# Map: Halifax Final Image | 35715 | 1203
# Map: Halifax Interactive | 26765 | 8950 *to image*
# Map: Vancouver Image | 25956 | 809
# Map: Vancouver Select Image | 25143 | 813
# Map: Vancouver Interactive | 14049 | 11094 *to image*
# Map: DA Image | 12883 | 1166
#

# ### Original df_loc (Canada) Description
#
# [Statistics Canada](https://www.statcan.gc.ca/eng/start) provides several geographic divisions for reporting of census statistics:
#
# Dissemination Area: DAs are small, relatively stable geographic unit composed of one or more adjacent dissemination blocks with an average population of 400 to 700 persons based on data from the previous Census of Population Program. It is the smallest standard geographic area for which all census data are disseminated.
#
# Aggregate Dissemination Area: ADAs cover the entire country and, where possible, have a population count between 5,000 and 15,000 people, and respect provincial, territorial, census division (CD), census metropolitan area (CMA) and census agglomeration (CA) with census tract (CT) boundaries in effect for the 2016 Census.
#
# Census Tract: CTs are small, relatively stable geographic areas that usually have a population of less than 10,000 persons, based on data from the previous Census of Population Program. They are located in census metropolitan areas and in census agglomerations that had a core population of 50,000 or more in the previous census.
#
# Census Metropolitan Area/ Census Agglomeration: A CMA or CA is formed by one or more adjacent municipalities centred on a population centre (known as the core). A CMA must have a total population of at least 100,000 of which 50,000 or more must live in the core based on adjusted data from the previous Census of Population Program. A CA must have a core population of at least 10,000 also based on data from the previous Census of Population Program. To be included in the CMA or CA, other adjacent municipalities must have a high degree of integration with the core, as measured by commuting flows derived from data on place of work from the previous Census Program.
#
# Census Subdivision: CSD is the general term for municipalities (as determined by provincial/territorial legislation) or areas treated as municipal equivalents for statistical purposes (e.g., Indian reserves, Indian settlements and unorganized territories).
#
# Census Consolidated Subdivision:A CCS is a group of adjacent census subdivisions within the same census division. Generally, the smaller, more densely-populated census subdivisions (towns, villages, etc.) are combined with the surrounding, larger, more rural census subdivision, in order to create a geographic level between the census subdivision and the census division.
#
# Census Division: CDs are a group of neighbouring municipalities joined together for the purposes of regional planning and managing common services (such as police or ambulance services).  Census divisions are intermediate geographic areas between the province/territory level and the municipality (census subdivision).
#
# Economic Region: An ER is a grouping of complete census divisions (CDs), with one exception in Ontario, created as a standard geographic unit for analysis of regional economic activity.
#
# Note: CSDNAME = 'Toronto' gives the same area as the FSAs, which makes sense as this is probably the definition of 'Toronto' 
#
# The meaning of prefixes for NAME (N), UID (U), PUID (PU), TYPE (T), and CODE (C):
# * DA U Dissemination Area unique identifier (composed of the 2-digit province/territory unique identifier followed by the 2-digit census division code and the 4-digit dissemination area code)
# * PR U,N Province or territory
# * CD U,N,T Census Division
# * CCS U,N Census Consolidated Subdivision
# * CSD U,N,T Census Subdivision
# * ER U,N Economic Region
# * SAC T,C Statistical Area Classification: Part of are a component of a census metropolitan area, a census agglomeration, a census metropolitan influenced zone or the territories?
# * CMA U,PU,N,T Census Metropolitan Area or Census Agglomeration name, PU Uniquely identifies the provincial or territorial part of a census metropolitan area or census agglomeration (composed of the 2-digit province or territory unique identifier followed by the 3-digit census metropolitan area or census agglomeration unique identifier)
# * CT U,N Census Tract within census metropolitan area/census agglomeration
# * ADA U Aggregate dissemination area unique identifier
#
# We will use census Distribution Areas as proxies for neighborhoods for cities in Canada.  In previous work where the Forward Sortation Areas (first three characters of the postal code) were used as neighborhood proxies, the sizes of many areas were quite large (several kilometers across) and therefore are likely internally non-homogeneous from a features perspective at the walking-distance (500 m) length scale.  To convert to neighborhood names we can look up the associated census tract as seen on [this](https://en.wikipedia.org/wiki/Demographics_of_Toronto_neighbourhoods) Wikipedia page.
#
# File lda_000b16g_e.gml was downloaded from the [Statistics Canada: Boundary Files](https://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2016-eng.cfm) site.
#
# Exploring the gml file and computing the area and centroid of the distribution areas can be done using the [geopandas module](https://geopandas.org/).  Geopandas builds upon [osgeo](https://gdal.org/python/index.html) which can also be used to explore and compute with the gml file, but in testing was 46 times faster due to vectorization of many calculations compared to a naive approach (see Appendix).
#
# Latitude and Longitude need to be obtained from a projection with those units, e.g. [EPSG-4326](https://epsg.io/4326).  Area can be calculated from an equal-area projection, e.g. [EPSG-6931](https://epsg.io/6931) (though a geodesic area calculation would be more accurate for larger regions, that would have to come from an additional package such as [proj](https://proj.org/), but since all regions here are small such the curvature of the earth is negligible, and altitude indtroduces an additional error likely comparable or larger to the curvature error, we will proceed with a simpler equal-area projection calculation).  The geometry is saved as text in the [Well-Known Text (WKT)](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry) representation.

# # Original Cities_Compared

# + [markdown] tags=[]
# Cities compared:
# * Canada
#     * Toronto
#     * Montréal
#     * Vancouver
#     * Halifax
# * United States of America
#     * New York City
#     * Boston
#     * Chicago
#     * San Fransisco
# * France
#     * Paris
# * England
#     * London
#
# First, get a list of Forward Sortation Areas as proxies for neighborhoods for cities in Canada.
# -

# # Original parsing code

# ### Osgeo package
#
# Before finding the geopandas library, I was using [gdal](https://gdal.org/python/index.html) (which is one of many dependencies of geopandas).  Using gdal's ogr and osr modules in the osgeo package I read in the gml file, converted coordinates to accurately compute area, and converted coordinates again to latitude and longitude.  As a guide to learning the library I adapted code from some [examples](https://pcjericks.github.io/py-gdalogr-cookbook/geometry.html#quarter-polygon-and-create-centroids).  Conversion to geojson may be done explicitly as demonstrated [here](https://gis.stackexchange.com/questions/77974/converting-gml-to-geojson-using-python-and-ogr-with-geometry-transformation), but by the time I was going to convert to geojson I found it much easier with geopandas, and this approach was abandoned.

# #### Parsing GML file

#     from osgeo import ogr
#     from osgeo import osr
#     
#     def parseGMLtoDF(fn_gml,limit=None):
#         '''Reads in a GML file and outputs a DataFrame
#
#         Adds columns for Geometry (in WKT format), Latitude (of centroid), Longitude (of centroid) and Area (in square meters)
#         '''
#         # Add lightweight progress bar, source copied from GitHub
#         import ipypb
#
#         # Read in the file
#         source = ogr.Open(fn_gml)
#         layer = source.GetLayer(0) # there is only one Layer (the FeatureCollection)
#
#         # Get a list of field names to extract (also could be gotten from schema e.g. lda_000b16g_e.xsd)
#         #   Fields are in order (sequence) and none are required (minOccurs=0)
#         layerDefinition = layer.GetLayerDefn()
#         layerFields = [layerDefinition.GetFieldDefn(i).GetName() for i in range(layerDefinition.GetFieldCount())]
#
#         # Initialize the output dataframe
#         df = pd.DataFrame(columns=[*layerFields, 'Geometry', 'Latitude', 'Longitude', 'Area'])
#
#         # Extract data from the features
#         parse_limit = layer.GetFeatureCount()
#         if not limit is None:
#             parse_limit = limit
#
#         for i in ipypb.track(range(parse_limit)):
#             # Get the feature to be processed
#             feature = layer.GetNextFeature()
#
#             # Copy all fields into an empty dataframe
#             df_tmp = pd.DataFrame(columns=df.columns)
#
#             # Copy all fields individually, in case some are missing
#             for i, fieldname in enumerate(layerFields):
#                 if feature.IsFieldSet(i):
#                     df_tmp.loc[0,fieldname] = feature.GetFieldAsString(fieldname)
#
#             # Get the latitude and longitude
#             inref = feature.GetGeometryRef().GetSpatialReference() # EPSG 3347
#             llref = osr.SpatialReference()
#             llref.ImportFromEPSG(4326) # coordinates for this geometry are latitude, longitude
#             lltransform = osr.CoordinateTransformation(inref, llref)
#             geom = feature.GetGeometryRef().Clone()
#             geom.Transform(lltransform)
#             centroid = geom.Centroid()
#             df_tmp['Latitude'] = centroid.GetX()
#             df_tmp['Longitude'] = centroid.GetY()
#             df_tmp['Geometry'] = geom.ExportToWkt()
#
#             # Get the area by converting to a locally-appropriate coordinate system
#             inref = feature.GetGeometryRef().GetSpatialReference()
#             arref = osr.SpatialReference()
#             arref.ImportFromEPSG(6931) # TODO: let this choose appropriate projection based on centroid for full-globe applicability, or alter to be a geodesic calculation.  Lambert Cylindrical Equal Area 6933, U.S. National Atlas Equal Area Projection 2163, NSIDC EASE-Grid North 3408, WGS 84 / NSIDC EASE-Grid 2.0 North 6931, WGS 84 / NSIDC EASE-Grid 2.0 South 6932, WGS 84 / NSIDC EASE-Grid 2.0 Global/Temperate 6933
#             artransform = osr.CoordinateTransformation(inref, arref)
#             geom = feature.GetGeometryRef().Clone() # Fresh instance to prevent error accumulation from multiple transforms
#             geom.Transform(artransform)
#             df_tmp['Area'] = geom.Area() # TODO: investigate proj module, related to osgeo, for geodesic area
#
#             # Add the feature dataframe to the output
#             df = df.append(df_tmp, ignore_index=True)
#
#         return df

# #### Viewing fields
#
#     layerDefinition = s2.GetLayerDefn()
#
#     for i in range(layerDefinition.GetFieldCount()):
#         print(layerDefinition.GetFieldDefn(i).GetName(),f" \t",s2.GetFeature(0).GetFieldAsString(s2.GetFeature(0).GetFieldIndex(layerDefinition.GetFieldDefn(i).GetName())))

df_CA_DA = parseGMLtoDF('lda_000b16g_e.gml')

display(df_CA_DA.head(2))
print(df_CA_DA.shape)

# #### Layered parsing code
#
# From when I thought dataframe copying was the limiting step
#
#     # Read in the file
#     source = ogr.Open('lda_000b16g_e.gml')
#     layer = source.GetLayer(0) # there is only one Layer (the FeatureCollection)
#
#     # Get a list of field names to extract (see lda_000b16g_e.xsd)
#     #   Fields are in order (sequence) and none are required (minOccurs=0)
#     layerDefinition = layer.GetLayerDefn()
#     layerFields = [layerDefinition.GetFieldDefn(i).GetName() for i in range(layerDefinition.GetFieldCount())]
#
#     # %%time
#     import ipypb # Lightweight progress bar, source copied from GitHub
#
#     # Initialize the output dataframe - REUSE FOR ALL OF CANADA
#     df_CA_DA = pd.DataFrame(columns=[*layerFields, 'Geometry', 'Latitude', 'Longitude', 'Area'])
#
#     # Only process a few entries
#     parse_limit = 5
#
#     # Setup for 2-step dataframe appending, should be ~O(n^(1+e)) instead of ~O(n^2)
#     n_features = layer.GetFeatureCount()
#     #max_n_merge = 1000
#     n_merge = n_features**0.5//1
#     #n_merge = min(n_features**0.5//1,max_n_merge) # How many rows to process before a merge
#     df_tmp_collect = pd.DataFrame(columns=df_CA_DA.columns)
#
#     # Extract data from the features
#     layer.ResetReading()
#     for i in ipypb.track(range(layer.GetFeatureCount())):
#         # Get the feature to be processed
#         feature = layer.GetNextFeature()
#
#         # Copy all fields into an empty dataframe
#         df_tmp = pd.DataFrame(columns=df_CA_DA.columns)
#
#         # Copy all fields individually, in case some are missing
#         for i, fieldname in enumerate(layerFields):
#             if feature.IsFieldSet(i):
#                 df_tmp.loc[0,fieldname] = feature.GetFieldAsString(fieldname)
#
#         # Get the latitude and longitude
#         inref = feature.GetGeometryRef().GetSpatialReference() # EPSG 3347
#         llref = osr.SpatialReference()
#         llref.ImportFromEPSG(4326) # coordinates for this geometry are latitude, longitude
#         lltransform = osr.CoordinateTransformation(inref, llref)
#         geom = feature.GetGeometryRef().Clone()
#         geom.Transform(lltransform)
#         centroid = geom.Centroid()
#         df_tmp['Latitude'] = centroid.GetX()
#         df_tmp['Longitude'] = centroid.GetY()
#         df_tmp['Geometry'] = geom.ExportToWkt()
#
#         # Get the area by converting to a locally-appropriate coordinate system
#         inref = feature.GetGeometryRef().GetSpatialReference()
#         arref = osr.SpatialReference()
#         arref.ImportFromEPSG(6931) # Lambert Cylindrical Equal Area 6933, U.S. National Atlas Equal Area Projection 2163, NSIDC EASE-Grid North 3408, WGS 84 / NSIDC EASE-Grid 2.0 North 6931, WGS 84 / NSIDC EASE-Grid 2.0 South 6932, WGS 84 / NSIDC EASE-Grid 2.0 Global/Temperate 6933
#         artransform = osr.CoordinateTransformation(inref, arref)
#         geom = feature.GetGeometryRef().Clone() # Fresh instance to prevent error accumulation from multiple transforms
#         geom.Transform(artransform)
#         df_tmp['Area'] = geom.Area()
#
#         # Add the feature dataframe to the output
#         #df_tmp_collect = df_tmp_collect.append(df_tmp, ignore_index=True)
#         df_tmp_collect = pd.concat([df_tmp_collect, df_tmp], ignore_index=True)
#
#         # Check if 2nd-step dataframe needs to be reset
#         if df_tmp_collect.shape[0] >= n_merge:
#             #df_CA_DA = df_CA_DA.append(df_tmp_collect, ignore_index=True)
#             df_CA_DA = pd.concat([df_CA_DA,df_tmp_collect], ignore_index=True)
#             df_tmp_collect = pd.DataFrame(columns=df_CA_DA.columns)
#
#         if not parse_limit is None:
#             parse_limit -= 1
#             if parse_limit<=0:
#                 break
#     #df_CA_DA = df_CA_DA.append(df_tmp_collect, ignore_index=True)
#     df_CA_DA = pd.concat([df_CA_DA,df_tmp_collect], ignore_index=True)
#
# ##### Timing Notes
#
# limit | n_merge | s/it | s total
# ---|---|---|---
# 500 | 1 | 0.13 | 1:06
# 500 | 2 | 0.13 | 1:02
# 500 | 10 | 0.12 | 1:01
# 500 | 50 | 0.12 | 0:59.2
# 500 | 237 | 0.11 | 0:57.5
# 500 | 237 | 0.13 | 1:03
# 1000 | 1 | 0.25 | 4:05
# 1000 | 237 | 0.23 | 3:52 # maybe this was max instead of min?
# 1000 | 1000 | 0.23 | 3:54
# 1000 | 1000 | 0.24 | 3:55  # pd.concat instead of dataframe.append
# 1000 | 237 | 0.22 | 3:38  # pd.concat instead of dataframe.append
# 1000 | - | 0.27 | 4:25  # original, no 2-levels
# 500  | - | 0.15 | 1:15  # original, no 2-levels
# 500 | 237 | 0.07 | 0:34.3  # pd.concat instead of dataframe.append
# 2000 | 237 | 0.07 | 2:13  # pd.concat instead of dataframe.append # Have not been resetting df_CA_DA!
# 4000 | 237 | 0.06 | 4:15  # pd.concat instead of dataframe.append # Resetting df from here on, and concat too, and GetNextFeature instead of GetFeature(int) which may have been the bottleneck then...
# 1000 | 237 | 0.06 | 1:01
# 500 | 237 | 0.06 | 0:31.3
# 1000 | 50 | 0.06 | 1:00
# 1000 | 1 | 0.06 | 1:04
# 1000 | 1 | 0.23 | 3:52  # Try once with GetFeature(i), yes, this was the culprit, it must parse anew each time.
# 4000 | 1 | 0.09 | 5:41
# 4000 | 5 | 0.06 | 4:11
# 4000 | 20 | 0.06 | 3:41
# 4000 | 50 | 0.06 | 3:54
# 4000 | 50 | 0.06 | 3:57
# 4000 | 100 | 0.06 | 3:58
# 4000 | 175 | 0.06 | 4:03
# 4000 | 237 | 0.06 | 4:08
# 4000 | 500 | 0.06 | 4:18
# 4000 | 500 | 0.06 | 3:50
# 4000 | 1000 | 0.06 | 4:16
# 4000 | 4000 | 0.07 | 4:58 # Oh, also hadn't reset the counter... adding that before last 500 and 30, that's been introducing variability too
# 4000 | 1 | 0.06 | 4:16
# 4000 | 30 | 0.06 | 4:17, 4:09 # Might change when tabs are changed (first with changing, second staying on same tab)
# 4000 | 60 | 0.06 | 3:53, 4:03
# 4000 | 200 | 0.06 | 3:56
# 4000 | 500 | 0.06 | 3:44
# 4000 | 1000 | 0.06 | 3:55
# 4000 | 4000 | 0.06 | 4:03
# 4000 | - | 0.06 | 4:04 # original, no 2-levels, with GetNextFeature
#
# Gain is marginal with the double layering... might be important for very large sets, but at 4000 features results in at most a 10% speedup.  Note that geopandas speedup is ~4,600%.

# #### Getting an appropriate EPSG projection
#
# Using UTM projections (of which there are 30, each subtending 6 degrees of longitude) was a commonly proposed method for projection for area calculation.  I went with the simpler EPSG-6931/2/3, which is easier due to one zone encompassing all of Canada (though with unknown accuracy tradeoff).  This approach is included just for reference.
#
#     srcSR = osr.SpatialReference()
#     srcSR.ImportFromEPSG(4326) # WGS84 Geographic
#     destSR   = osr.SpatialReference()
#
#     lyr.ResetReading()
#     for feat in lyr:
#         geom = feat.GetGeometryRef()
#         if not geom.IsEmpty():                 # make sure the geometry isn't empty
#             geom.AssignSpatialReference(srcSR) # you only need to do this if the shapefile isn't set or is set wrong
#             env = geom.GetEnvelope()           # get the Xmin, Ymin, Xmax, Ymax bounds
#             CentX = ( env[0] + env[2] ) / 2    # calculate the centre X of the whole geometry
#             Zone  = int((CentX + 180)/6) + 1   # see http://gis.stackexchange.com/questions/13291/computing-utm-zone-from-lat-long-point
#             EPSG  = 32600 + Zone               # get the EPSG code from the zone and the constant 32600 (all WGS84 UTM North start with 326)
#             destSR.ImportFromEPSG(EPSG)        # create the 'to' spatial reference
#             geom.TransformTo(destSR)           # project the geometry
#             print geom.GetArea()               # get the area in square metres

# #### Read in the Canadian data file: Original
#     gdf = geopandas.read_file('lda_000b16g_e.gml')
#     gdf.rename_geometry('Geometry', inplace=True) # Default geometry column name is 'geometry'; changed for consistent capitalization of columns
#     gdf.set_geometry('Geometry') # Renaming is insufficient; this sets special variable gdf.geometry = gdf['Geometry']
#     gdf['Area'] = gdf['Geometry'].to_crs(epsg=6931).area
#     gdf['Centroid'] = gdf['Geometry'].centroid
#     gdf['Geometry'] = gdf['Geometry'].to_crs(epsg=4326)
#     gdf['Centroid'] = gdf['Centroid'].to_crs(epsg=4326) # Only the set geometry is converted with gdf.to_crs(); all other geometry-containing columns must be converted explicitly; here we convert all columns explicitly
#     gdf['Centroid Latitude'] = gdf['Centroid'].geometry.y
#     gdf['Centroid Longitude'] = gdf['Centroid'].geometry.x
#     gdf.drop(columns = 'Centroid', inplace=True) # Because WKT Point cannot be serialized to JSON, we drop the Centroid column and keep only its float components
#     gdf_CA_DA = gdf # Rename because we may be generating additional variables
#     gdf_CA_DA.head(2)

# # Old Vancouver Extended

# Old Vancouver Extended:

cityname = 'Vancouver'
gdf_select = selectRegion(gdf_CA_DA, 'CDNAME', method='contains', names='Greater Vancouver')
mp = folium.Map(location=getCityByName(cityname)['centroid'])
mp.fit_bounds(getGDFBounds(gdf_select.to_crs(epsg=4326)))
m = folium.GeoJson(
            gdf_select.to_crs(epsg=4326).to_json(),
            style_function=lambda feature: {
                'fillColor': 'yellow',
                'fillOpacity': 0.8,
                'color':'black',
                'weight': 1,
                'opacity': 0.2,
            },
        ).add_to(mp)
m.add_child(folium.features.GeoJsonTooltip(tooltip_DA))
mp

# # Unused Choropleth map development

# +
# Not used - this naive choropleth displays only one map and doesn't allow much scale customization

city_index = 1

# create a plain world map
city_map = folium.Map(location=cities[city_index]['centroid'], control_scale=True)
city_map.fit_bounds(cities[city_index]['bounds'])

folium.Choropleth(
    geo_data=cities[city_index]['geojson'],
    data=cities[city_index]['data'],
    columns=['DAUID', 'Area'],
    key_on='feature.properties.DAUID',
    fill_color='YlOrRd', 
    fill_opacity=0.7, 
    line_opacity=0.2,
    legend_name='Area [square meters]'
).add_to(city_map)

# display map
city_map

# +
# Not used - links choropleth to a layercontrol, but does not get the legend on the correct layer

from branca.element import MacroElement

from jinja2 import Template

class BindColormap(MacroElement):
    """Binds a colormap to a given layer.

    Parameters
    ----------
    colormap : branca.colormap.ColorMap
        The colormap to bind.
    """
    def __init__(self, layer, colormap):
        super(BindColormap, self).__init__()
        self.layer = layer
        self.colormap = colormap
        self._template = Template(u"""
        {% macro script(this, kwargs) %}
            {{this.colormap.get_name()}}.svg[0][0].style.display = 'block';
            {{this._parent.get_name()}}.on('overlayadd', function (eventLayer) {
                if (eventLayer.layer == {{this.layer.get_name()}}) {
                    {{this.colormap.get_name()}}.svg[0][0].style.display = 'block';
                }});
            {{this._parent.get_name()}}.on('overlayremove', function (eventLayer) {
                if (eventLayer.layer == {{this.layer.get_name()}}) {
                    {{this.colormap.get_name()}}.svg[0][0].style.display = 'none';
                }});
        {% endmacro %}
        """)  # noqa
        
# https://gitter.im/python-visualization/folium?at=5a36090a03838b2f2a04649d
print('Area Selection Criterion: CSDNAME contains the city name')
import branca.element as bre
from branca.colormap import LinearColormap

f = bre.Figure()
sf = [[],[],[]]
city_map = [[],[],[]]
colormap = [[],[],[]]
cp = [[],[],[]]
for i, city in enumerate(cities):
    sf[i] = f.add_subplot(1,len(cities),1+i)
    city_map[i] = folium.Map(location=city['centroid'], control_scale=True)
    sf[i].add_child(city_map[i])
    
    title_html = '''<h3 align="center" style="font-size:16px;charset=utf-8"><b>{}</b></h3>'''.format(city['name'])
    city_map[i].get_root().html.add_child(folium.Element(title_html))
    city_map[i].fit_bounds(city['bounds'])
    cp[i] = folium.Choropleth(geo_data=city['geojson'],
                                data=city['data'],
                                columns=['DAUID', 'Area'],
                                key_on='feature.properties.DAUID',
                                fill_color='YlOrRd',
                                fill_opacity=0.7, 
                                line_opacity=0.2
                               )
    city_map[i].add_child(cp[i])
    city_map[i].add_child(folium.map.LayerControl())
    city_map[i].add_child(BindColormap(cp[i],cp[i].color_scale))

    city['map'] = city_map[i]

display(f)
# -



# # To get a list of functions in a module:
# from inspect import getmembers, isfunction
#
# from somemodule import foo
# print(getmembers(foo, isfunction))

# Postal Codes
import requests
# from bs4 import BeautifulSoup # had to install to environment in Anaconda
import lxml # had to install to environment in Anaconda, backdated to 4.6.1 (4.6.2 current) for pandas read_html()
import html5lib # had to install to environment in Anaconda (1.1 current) for pandas read_html()

# List of coordinate tuples of nonzero elements of matrix r
list(zip(*np.nonzero(r)))

geopandas.show_versions()


# ### Header that won't display on GitHub (tags not supported)
#
#     # <center>Applied Data Science Capstone</center>
#
#     #### <center>In completion of requirements for the IBM Data Science Professional Certificate on Coursera</center>
#
#     <center><NAME></center>
#
#     <hr>
#

# ### Loading/Saving

# #### Loading/Saving Without Compression

# Standalone
#
#     import dill
#     with open('GDF_FSA-DA-D.db','wb') as file:
#         dill.dump(gdf_union,file)
#     with open('GDF_FSA-DA-D_times.db','wb') as file:
#         dill.dump(times,file)
#     with open('GDF_FSA-DA-D_areas.db','wb') as file:
#         dill.dump(areas,file)
#
#     import dill
#     with open('GDF_FSA-DA-D.db','r') as file:
#         gdf_union = dill.load(file)
#     with open('GDF_FSA-DA-D_times.db','r') as file:
#         times = dill.load(file)
#     with open('GDF_FSA-DA-D_areas.db','r') as file:
#         areas = dill.load(file)

# Option Block
#
#     try: # Results have been calculated previously, load them to save time
#         with open(DIR_RESULTS+'GDF_FSA-DA_D.db','rb') as file:
#             gdf_union = dill.load(file)
#         with open(DIR_RESULTS+'GDF_FSA-DA_D_times.db','rb') as file:
#             times = dill.load(file)
#         with open(DIR_RESULTS+'GDF_FSA-DA_D_areas.db','rb') as file:
#             areas = dill.load(file)
#         print('Results loaded from file')
#
#     except (FileNotFoundError, IOError):  # Results not found in file, regenerate them
#         gdf_union, times, areas = intersectGDF(gdf_CA_FSA_D,'CFSAUID',gdf_CA_DA_D,'DAUID',verbosity=1)
#
#         with open(DIR_RESULTS+'GDF_FSA-DA_D.db','wb+') as file:
#             dill.dump(gdf_union,file)
#         with open(DIR_RESULTS+'GDF_FSA-DA_D_times.db','wb+') as file:
#             dill.dump(times,file)
#         with open(DIR_RESULTS+'GDF_FSA-DA_D_areas.db','wb+') as file:
#             dill.dump(areas,file)
#         print('Results saved to file')

# #### Loading/saving with compression
#
# Original load/compute/save logic for overlap computation:
#
#     try: # Results have been calculated previously, load them to save time
#         with open(DIR_RESULTS+'GDF_FSA-DA_D_times.db.gz','rb') as file:
#             times = dill.loads(gzip.decompress(file.read()))
#         with open(DIR_RESULTS+'GDF_FSA-DA_D_areas.db.gz','rb') as file:
#             areas = dill.loads(gzip.decompress(file.read()))
#         try: # If not available (e.g. on github, due to size), reconstruct the intersection gdf
#             with open(DIR_RESULTS+'GDF_FSA-DA_D.db.gz','rb') as file:
#                 gdf_union = dill.loads(gzip.decompress(file.read()))
#         except (FileNotFoundError, IOError):
#             print('Recomputing gdf_union using areas loaded from file')
#             gdf_union_areas, times_areas, areas_areas= intersectGDFareas(gdf_CA_FSA_D,'CFSAUID',gdf_CA_DA_D,'DAUID',areas_in=areas,verbosity=1)
#             gdf_union = gdf_union_areas
#             with open(DIR_RESULTS+'GDF_FSA-DA_D.db.gz','wb+') as file:
#                 file.write(gzip.compress(dill.dumps(gdf_union)))
#         print('Results loaded from file')
#
#     except (FileNotFoundError, IOError):  # Results not found in file, regenerate them
#         gdf_union, times, areas = intersectGDF(gdf_CA_FSA_D,'CFSAUID',gdf_CA_DA_D,'DAUID',verbosity=1)
#
#         with open(DIR_RESULTS+'GDF_FSA-DA_D.db.gz','wb+') as file:
#             file.write(gzip.compress(dill.dumps(gdf_union)))
#         with open(DIR_RESULTS+'GDF_FSA-DA_D_times.db.gz','wb+') as file:
#             file.write(gzip.compress(dill.dumps(times)))
#         with open(DIR_RESULTS+'GDF_FSA-DA_D_areas.db.gz','wb+') as file:
#             file.write(gzip.compress(dill.dumps(areas)))
#         print('Results saved to file')

# ## ExtendBounds Development

# ### Original extendBounds Function and Testing

def extendBounds(bounds,method='nearestLeadingDigit',scale=10):
    '''Extend bounds (low, high) to give round numbers for scalebars

    The low bound is decreased and the high bound is increased according to
        method to the first number satisfying the method conditions.
    Returns a new bound which includes the old bounds in its entirety

    Parameters
    ----------
    bounds: (low_bound, high_bound) list or tuple
    method: str describing the extension method
        'nearestLeadingDigit': Bounds are nearest numbers with leading digit followed by zeros
        'nearestPower': Bounds are nearest powers of scale (scale must be > 1).  For negative numbers, the sign and direction are reversed, the extension performed, then the sign of the result is reversed back.
        'nearestMultiple': Bounds are nearest multiple of scale (scale must be > 0)
        'round': Bounds are rounded
    scale: numeric as described in method options

    Returns
    -------
    2-element tuple of extended bounds e.g. (newlow,newhigh)
    '''
    if bounds[0]>bounds[1]:
        print('bounds must be ordered from least to greatest')
        return None    
    if method=='nearestLeadingDigit':
        iszero = np.array(bounds)==0
        isnegative = np.array(bounds) < 0
        offsets = [1 if isnegative[0] else 0, 0 if isnegative[1] else 1]
        power = [0 if z else np.floor(np.log10(abs(b))) for b, z in zip(bounds, iszero)]
        firstdigit = [abs(b)//np.power(10,p) for b, p in zip(bounds, power)]
        exceeds = [abs(b)>f*np.power(10,p) for b, f, p in zip(bounds, firstdigit, power)]
        newbounds = [abs(b) if not t else (f+o)*np.power(10,p) for b, t, n, f, o, p in zip(bounds, exceeds, isnegative, firstdigit, offsets, power)]
        newbounds = [-n if t else n for n, t in zip(newbounds,isnegative)]
    elif method=='nearestPower':
        try:
            scale = float(scale)
            if scale<=1:
                print('scale should be greater than 1')
                return None
        except ValueError:
            print('scale should be a number greater than 1')
            return None
        isnegative = np.array(bounds) < 0
        roundfuns = [np.ceil if isnegative[0] else np.floor, np.floor if isnegative[1] else np.ceil]
        newbounds = [0 if b==0 else np.power(scale, r(np.log10(abs(b))/np.log10(scale))) for b, r in zip(bounds,roundfuns)]
        newbounds = [-n if t else n for n, t in zip(newbounds,isnegative)]
    elif method=='nearestMultiple':
        try:
            scale = float(scale)
            if scale<=0:
                print('scale should be greater than 0')
                return None
        except ValueError:
            print('scale should be a number greater than 0')
            return None
        newbounds = [scale*(np.floor(bounds[0]/scale)), scale*(np.ceil(bounds[1]/scale))]
    elif method=='round':
        newbounds = [np.floor(bounds[0]), np.ceil(bounds[1])]
    else:
        print('Invalid method, see help(extendBounds)')
        return None
    return newbounds


#     print("Testing invalid method")
#     print("  Expect errors:")
#     print(f"{extendBounds([11,130],'invalid')}")
#     print()
#
#     print("Testing method 'nearestLeadingDigit'")
#     print("  Expect errors:")
#     print(f"{extendBounds([9,-930],'nearestLeadingDigit')}")
#     print(f"{extendBounds([-9,-930],'nearestLeadingDigit')}")
#     print("  Expect success:")
#     print(f"{extendBounds([11,130],'nearestLeadingDigit')}",)
#     print(f"{extendBounds([11,130],'nearestLeadingDigit',-1)}")
#     print(f"{extendBounds([9,930],'nearestLeadingDigit')}")
#     print(f"{extendBounds([-9,930],'nearestLeadingDigit')}")
#     print(f"{extendBounds([-990,-930],'nearestLeadingDigit')}")
#     print(f"{extendBounds([-990,0.05],'nearestLeadingDigit')}")
#     print(f"{extendBounds([0,0.052],'nearestLeadingDigit')}")
#     print()
#
#     print("Testing method 'nearestPower'")
#     print("  Expect errors:")
#     print(f"{extendBounds([11,130],'nearestPower',-2)}")
#     print(f"{extendBounds([11,130],'nearestPower',0)}")
#     print(f"{extendBounds([11,130],'nearestPower',1)}")
#     print(f"{extendBounds([-11,-130],'nearestPower',10)}")
#     print("  Expect success:")
#     print(f"{extendBounds([11,130],'nearestPower')}")
#     print(f"{extendBounds([10,100],'nearestPower')}")
#     print(f"{extendBounds([11,130],'nearestPower',1.1)}")
#     print(f"{extendBounds([11,130],'nearestPower',2)}")
#     print(f"{extendBounds([11,130],'nearestPower',10)}")
#     print(f"{extendBounds([11,130],'nearestPower',10.)}")
#     print(f"{extendBounds([-11,130],'nearestPower',10)}")
#     print(f"{extendBounds([-5100,-130],'nearestPower',10)}")
#     print(f"{extendBounds([-.0101,-0.00042],'nearestPower',10)}")
#     print(f"{extendBounds([0,0.00042],'nearestPower',10)}")
#     print()
#
#     print("Testing method 'nearestMultiple'")
#     print("  Expect errors:")
#     print(f"{extendBounds([11,132],'nearestMultiple',-2)}")
#     print(f"{extendBounds([11,132],'nearestMultiple',0)}")
#     print(f"{extendBounds([0,-10],'nearestMultiple',100)}")
#     print("  Expect success:")
#     print(f"{extendBounds([11,132],'nearestMultiple')}")
#     print(f"{extendBounds([10,130],'nearestMultiple')}")
#     print(f"{extendBounds([11.55,132.55],'nearestMultiple',0.1)}")
#     print(f"{extendBounds([11.55,132.55],'nearestMultiple',1)}")
#     print(f"{extendBounds([11.55,132.55],'nearestMultiple',100)}")
#     print(f"{extendBounds([-11,132],'nearestMultiple',10)}")
#     print(f"{extendBounds([-1121,-132],'nearestMultiple',10)}")
#     print(f"{extendBounds([-10,-10],'nearestMultiple',10)}")
#     print(f"{extendBounds([-10,-10],'nearestMultiple',100)}")
#     print()
#
#     print("Testing method 'round'")
#     print("  Expect errors:")
#     print(f"{extendBounds([-11.1,-132.1],'round')}")
#     print("  Expect success:")
#     print(f"{extendBounds([11.1,132.1],'round')}")
#     print(f"{extendBounds([10,130],'round')}")
#     print(f"{extendBounds([11.1,132.1],'round',-2)}")
#     print(f"{extendBounds([-11.1,132.1],'round')}")
#     print(f"{extendBounds([-1100.1,-132.1],'round')}")

#     Testing invalid method
#       Expect errors:
#     Invalid method, see help(extendBounds)
#     None
#
#     Testing method 'nearestLeadingDigit'
#       Expect errors:
#     bounds must be ordered from least to greatest
#     None
#     bounds must be ordered from least to greatest
#     None
#       Expect success:
#     [10.0, 200.0]
#     [10.0, 200.0]
#     [9, 1000.0]
#     [-9, 1000.0]
#     [-1000.0, -900.0]
#     [-1000.0, 0.05]
#     [0, 0.06]
#
#     Testing method 'nearestPower'
#       Expect errors:
#     scale should be greater than 1
#     None
#     scale should be greater than 1
#     None
#     scale should be greater than 1
#     None
#     bounds must be ordered from least to greatest
#     None
#       Expect success:
#     [10.0, 1000.0]
#     [10.0, 100.0]
#     [10.834705943388395, 142.04293198443193]
#     [8.0, 256.0]
#     [10.0, 1000.0]
#     [10.0, 1000.0]
#     [-100.0, 1000.0]
#     [-10000.0, -100.0]
#     [-0.1, -0.0001]
#     [0, 0.001]
#
#     Testing method 'nearestMultiple'
#       Expect errors:
#     scale should be greater than 0
#     None
#     scale should be greater than 0
#     None
#     bounds must be ordered from least to greatest
#     None
#       Expect success:
#     [10.0, 140.0]
#     [10.0, 130.0]
#     [11.5, 132.6]
#     [11.0, 133.0]
#     [0.0, 200.0]
#     [-20.0, 140.0]
#     [-1130.0, -130.0]
#     [-10.0, -10.0]
#     [-100.0, -0.0]
#
#     Testing method 'round'
#       Expect errors:
#     bounds must be ordered from least to greatest
#     None
#       Expect success:
#     [11.0, 133.0]
#     [10.0, 130.0]
#     [11.0, 133.0]
#     [-12.0, 133.0]
#     [-1101.0, -132.0]

# ### New extendBounds Function

# +
def extendBound(bound,direction='up',method='nearestLeadingDigit',scale=10):
    '''Extend bound to next 'round' number
    
    Parameters
    ----------
    bound: float or float castable number or a list thereof
    direction: {'up','down',nonzero number} or a list of these values indicating the direction to round in
    method: str describing the extension method
        'nearestLeadingDigit': Bound is nearest numbers with leading digit followed by zeros
        'nearestPower': Bound is nearest integer power of scale (scale must be > 1).  For negative numbers, the sign and direction are reversed, the extension performed, then the sign of the result is reversed back.
        'nearestMultiple': Bound is nearest multiple of scale (scale must be > 0)
        'round': Bound is rounded using the default method
    scale: numeric as described in method options or a list thereof
    
    Returns
    -------
    float: the extended bound
    
    Notes
    -----
    All inputs, if not single-valued, must be lists of length equal to input bound
    TODO: extend so that method may also be a list
    TODO: replace prints with raising errors
    '''
    import numpy as np
    
    # Check and adjust the length of inputs
    unlist_bound = False
    if not(type(bound) in {list,tuple,range}):
        bound = [bound]
        unlist_bound = True
    acceptable_len = set((1,len(bound)))
    if not(type(direction) in {list,tuple,range}):
        direction = [direction]
    if not(len(direction) in acceptable_len):
        print('"direction" must have length 1 or length equal to the length of "bound"')
        return None
    if (type(method) in {str}):
        method = [method]
    if not(len(method) in acceptable_len):
        print('"method" must have length 1 or length equal to the length of "bound"')
        return None
    if not(type(scale) in {list,tuple,range}):
        scale = [scale]
    if not(len(scale) in acceptable_len):
        print('"scale" must have length 1 or length equal to the length of "bound"')
        return None
    if len(bound)>1:
        if len(direction)==1: direction = [direction[0] for b in bound]
        if len(scale)==1: scale = [scale[0] for b in bound]
        
    # If multiple methods are specified, recursively call this function for each method and reassemble results
    if len(bound)>1 and len(method)>1:
        ret = np.array([None for b in bound])
        for m in list(set(method)):
            ind = np.where(np.array(method)==m)
            ret[ind] = extendBound(list(np.array(bound)[ind]),list(np.array(direction)[ind]),m,list(np.array(scale)[ind]))
        return list(ret)
    
    # Convert direction to a logical array roundup
    try:
        roundup = [True if d=='up' else False if d=='down' else True if float(d)>0 else False if float(d)<0 else None for d in direction]
    except:
        print('direction must be "up", "down", or a non-negative number')
        return None
    if any([r==None for r in roundup]):
        print('direction must be "up", "down", or a non-negative number')
        return None
    
    # Cases for multiple methods handled above, return to string method
    method = method[0]
    
    # Execute the conversions
    if method=='nearestLeadingDigit':
        iszero = np.array(bound)==0
        isnegative = np.array(bound) < 0
        offsets = np.logical_xor(roundup, isnegative)
        power = [0 if z else np.floor(np.log10(abs(b))) for b, z in zip(bound, iszero)]
        firstdigit = [abs(b)//np.power(10,p) for b, p in zip(bound, power)]
        exceeds = [abs(b)>f*np.power(10,p) for b, f, p in zip(bound, firstdigit, power)]
        newbound = [abs(b) if not t else (f+o)*np.power(10,p) for b, t, n, f, o, p in zip(bound, exceeds, isnegative, firstdigit, offsets, power)]
        newbound = [-n if t else n for n, t in zip(newbound,isnegative)]
    elif method=='nearestPower':
        try:
            scale = [float(s) for s in scale]
            if any([s<=1 for s in scale]):
                print('scale should be greater than 1')
                return None
        except ValueError:
            print('scale should be a number or list of numbers greater than 1')
            return None
        isnegative = np.array(bound) < 0
        offsets = np.logical_xor(roundup, isnegative)
        roundfuns = [np.ceil if o else np.floor for o in offsets]
        newbound = [0 if b==0 else np.power(s, r(np.log10(abs(b))/np.log10(s))) for b, r, s in zip(bound,roundfuns,scale)]
        newbound = [-n if t else n for n, t in zip(newbound,isnegative)]
    elif method=='nearestMultiple':
        try:
            scale = [float(s) for s in scale]
            if any([s<=0 for s in scale]):
                print('scale should be greater than 0')
                return None
        except ValueError:
            print('scale should be a number or list of numbers greater than 0')
            return None
        roundfuns = [np.ceil if r else np.floor for r in roundup]
        newbound = [s*(r(b/s)) for b, r, s in zip(bound,roundfuns,scale)]
    elif method=='round':
        roundfuns = [np.ceil if r else np.floor for r in roundup]
        newbound = [f(b) for b, f in zip(bound, roundfuns)]
    else:
        print('Invalid method, see help(extendBound)')
        return None
    return newbound[0] if unlist_bound else newbound

def extendBounds(bounds,method='nearestLeadingDigit',scale=10):
    if bounds[0]>bounds[1]:
        print('bounds must be ordered from least to greatest')
        return None    
    return extendBound(bounds,direction=['down','up'],method=method,scale=scale)


# -

# ### Newest extendBounds and testing: list-castable

# +
def extendBound(bound,direction='up',method='nearestLeadingDigit',scale=10):
    '''Extend bound to next 'round' number
    
    Parameters
    ----------
    bound: float or float castable number or a list thereof
    direction: {'up','down',nonzero number} or a list of these values indicating the direction to round in
    method: str describing the extension method
        'nearestLeadingDigit': Bound is nearest numbers with leading digit followed by zeros
        'nearestPower': Bound is nearest integer power of scale (scale must be > 1).  For negative numbers, the sign and direction are reversed, the extension performed, then the sign of the result is reversed back.
        'nearestMultiple': Bound is nearest multiple of scale (scale must be > 0)
        'round': Bound is rounded using the default method
    scale: numeric as described in method options or a list thereof
    
    Returns
    -------
    float: the extended bound
    
    Notes
    -----
    All inputs, if not single-valued, must be list-castable and of equal length
    If all inputs are single-valued, the output is a float, otherwise it is a list of floats
    '''
    import numpy as np
    
    # Check and adjust the length of inputs
    unlist = False
    try:
        bound = list(bound)
    except:
        try:
            bound = [bound]
            unlist = True
        except:
            print("Input 'bound' must be numeric or convertible to list type.")
            return None
    try:
        if type(direction)==str:
            direction = [direction]
        direction = list(direction)
    except:
        try:
            direction = [direction]
        except:
            print("Input 'direction' must be a string or nonzero number or convertible to list type.")
            return None
    try:
        if type(method)==str:
            method = [method]
        method = list(method)
    except:
        try:
            method = [method]
        except:
            print("Input 'method' must be a string or convertible to list type.")
            return None
    try:
        scale = list(scale)
    except:
        try:
            scale = [scale]
        except:
            print("Input 'scale' must be numeric or convertible to list type.")
            return None
    inputs = [bound, direction, method, scale]
    lengths = [len(i) for i in inputs]
    set_lengths = set(lengths)
    max_len = max(set_lengths)
    set_lengths.remove(1)
    if len(set_lengths)>1:
        print('Inputs must be of the same length or of length one.  See help(extendBound)')
        return None
    if max_len>1: # can this be converted to a looped statement?
        if len(bound)==1:
            bound = bound*max_len
        if len(direction)==1:
            direction = direction*max_len
        if len(method)==1:
            method = method*max_len
        if len(scale)==1:
            scale = scale*max_len
        unlist = False

    # If multiple methods are specified, recursively call this function for each method and reassemble results
    if len(bound)>1 and len(set(method))>1:
        ret = np.array([None for b in bound])
        for m in list(set(method)):
            ind = np.where(np.array(method)==m)
            ret[ind] = extendBound(list(np.array(bound)[ind]),list(np.array(direction)[ind]),m,list(np.array(scale)[ind]))
        return list(ret)
    
    # Convert direction to a logical array roundup
    try:
        roundup = [True if d=='up' else False if d=='down' else True if float(d)>0 else False if float(d)<0 else None for d in direction]
    except:
        print('direction must be "up", "down", or a non-negative number')
        return None
    if any([r==None for r in roundup]):
        print('direction must be "up", "down", or a non-negative number')
        return None
    
    # Cases for multiple methods handled above, return to string method
    method = method[0]
    
    # Execute the conversions
    if method=='nearestLeadingDigit':
        iszero = np.array(bound)==0
        isnegative = np.array(bound) < 0
        offsets = np.logical_xor(roundup, isnegative)
        power = [0 if z else np.floor(np.log10(abs(b))) for b, z in zip(bound, iszero)]
        firstdigit = [abs(b)//np.power(10,p) for b, p in zip(bound, power)]
        exceeds = [abs(b)>f*np.power(10,p) for b, f, p in zip(bound, firstdigit, power)]
        newbound = [abs(b) if not t else (f+o)*np.power(10,p) for b, t, n, f, o, p in zip(bound, exceeds, isnegative, firstdigit, offsets, power)]
        newbound = [-n if t else n for n, t in zip(newbound,isnegative)]
    elif method=='nearestPower':
        try:
            scale = [float(s) for s in scale]
            if any([s<=1 for s in scale]):
                print('scale should be greater than 1')
                return None
        except ValueError:
            print('scale should be a number or list of numbers greater than 1')
            return None
        isnegative = np.array(bound) < 0
        offsets = np.logical_xor(roundup, isnegative)
        roundfuns = [np.ceil if o else np.floor for o in offsets]
        newbound = [0 if b==0 else np.power(s, r(np.log10(abs(b))/np.log10(s))) for b, r, s in zip(bound,roundfuns,scale)]
        newbound = [-n if t else n for n, t in zip(newbound,isnegative)]
    elif method=='nearestMultiple':
        try:
            scale = [float(s) for s in scale]
            if any([s<=0 for s in scale]):
                print('scale should be greater than 0')
                return None
        except ValueError:
            print('scale should be a number or list of numbers greater than 0')
            return None
        roundfuns = [np.ceil if r else np.floor for r in roundup]
        newbound = [s*(r(b/s)) for b, r, s in zip(bound,roundfuns,scale)]
    elif method=='round':
        roundfuns = [np.ceil if r else np.floor for r in roundup]
        newbound = [f(b) for b, f in zip(bound, roundfuns)]
    else:
        print('Invalid method, see help(extendBound)')
        return None
    return newbound[0] if unlist else newbound

def extendBounds(bounds,method='nearestLeadingDigit',scale=10):
    if bounds[0]>bounds[1]:
        print('bounds must be ordered from least to greatest')
        return None    
    return extendBound(bounds,direction=['down','up'],method=method,scale=scale)


# +
# Unit test of extendBounds
#   TODO: check out the builtin unittest module and convert this code to use that testing structure

# Get default arguments from https://stackoverflow.com/questions/12627118/get-a-function-arguments-default-value
import inspect
def get_default_args(func):
    signature = inspect.signature(func)
    return {
        k: v.default
        for k, v in signature.parameters.items()
        if v.default is not inspect.Parameter.empty
    }

defaults = get_default_args(extendBound)
def test_extendBound(bounds,direction=defaults['direction'],method=defaults['method'],scale=defaults['scale'],expected=None): # default arguments taken from extendBound; not sure how to get defaults when not supplied
    output = extendBound(bounds,direction,method,scale)
    print('Input:',bounds,direction,method,scale,'  Output:',output,'  Expected:',expected,'  Passed:',output==expected)
    return output==expected

defaults = get_default_args(extendBounds)
def test_extendBounds(bounds,method=defaults['method'],scale=defaults['scale'],expected=None): # default arguments taken from extendBounds; not sure how to get defaults when not supplied
    output = extendBounds(bounds,method,scale)
    print('Input:',bounds,method,scale,'  Output:',output,'  Expected:',expected,'  Passed:',output==expected)
    return output==expected

print('Testing function extendBounds\n')

passed = True
print("Testing invalid method")
print("  Expect errors:")
passed = test_extendBounds([11,130],'invalid',expected=None) and passed
print()

print("Testing method 'nearestLeadingDigit'")
print("  Expect errors:")
passed = test_extendBounds([9,-930],'nearestLeadingDigit',expected=None) and passed
passed = test_extendBounds([-9,-930],'nearestLeadingDigit',expected=None) and passed
print("  Expect success:")
passed = test_extendBounds([11,130],'nearestLeadingDigit',expected=[10,200]) and passed
passed = test_extendBounds([11,130],'nearestLeadingDigit',-1,expected=[10,200]) and passed
passed = test_extendBounds([9,930],'nearestLeadingDigit',expected=[9,1000]) and passed
passed = test_extendBounds([-9,930],'nearestLeadingDigit',expected=[-9,1000]) and passed
passed = test_extendBounds([-990,-930],'nearestLeadingDigit',expected=[-1000,-900]) and passed
passed = test_extendBounds([-990,0.05],'nearestLeadingDigit',expected=[-1000,0.05]) and passed
passed = test_extendBounds([0,0.052],'nearestLeadingDigit',expected=[0,0.06]) and passed
print()

print("Testing method 'nearestPower'")
print("  Expect errors:")
passed = test_extendBounds([11,130],'nearestPower',-2,expected=None) and passed
passed = test_extendBounds([11,130],'nearestPower',0,expected=None) and passed
passed = test_extendBounds([11,130],'nearestPower',1,expected=None) and passed
passed = test_extendBounds([-11,-130],'nearestPower',10,expected=None) and passed
print("  Expect success:")
passed = test_extendBounds([11,130],'nearestPower',expected=[10,1000]) and passed
passed = test_extendBounds([10,100],'nearestPower',expected=[10,100]) and passed
passed = test_extendBounds([11,130],'nearestPower',1.1,expected=[10.834705943388395, 142.04293198443193]) and passed
passed = test_extendBounds([11,130],'nearestPower',2,expected=[8,256]) and passed
passed = test_extendBounds([11,130],'nearestPower',10,expected=[10,1000]) and passed
passed = test_extendBounds([11,130],'nearestPower',10.,expected=[10,1000]) and passed
passed = test_extendBounds([-11,130],'nearestPower',10,expected=[-100,1000]) and passed
passed = test_extendBounds([-5100,-130],'nearestPower',10,expected=[-10000,-100]) and passed
passed = test_extendBounds([-.0101,-0.00042],'nearestPower',10,expected=[-0.1,-0.0001]) and passed
passed = test_extendBounds([0,0.00042],'nearestPower',10,expected=[0,0.001]) and passed
print()

print("Testing method 'nearestMultiple'")
print("  Expect errors:")
passed = test_extendBounds([11,132],'nearestMultiple',-2,expected=None) and passed
passed = test_extendBounds([11,132],'nearestMultiple',0,expected=None) and passed
passed = test_extendBounds([0,-10],'nearestMultiple',100,expected=None) and passed
print("  Expect success:")
passed = test_extendBounds([11,132],'nearestMultiple',expected=[10,140]) and passed
passed = test_extendBounds([10,130],'nearestMultiple',expected=[10,130]) and passed
passed = test_extendBounds([11.55,132.55],'nearestMultiple',0.1,expected=[11.5,132.6]) and passed
passed = test_extendBounds([11.55,132.55],'nearestMultiple',1,expected=[11,133]) and passed
passed = test_extendBounds([11.55,132.55],'nearestMultiple',100,expected=[0,200]) and passed
passed = test_extendBounds([-11,132],'nearestMultiple',10,expected=[-20,140]) and passed
passed = test_extendBounds([-1121,-132],'nearestMultiple',10,expected=[-1130,-130]) and passed
passed = test_extendBounds([-10,-10],'nearestMultiple',10,expected=[-10,-10]) and passed
passed = test_extendBounds([-10,-10],'nearestMultiple',100,expected=[-100,0]) and passed
print()

print("Testing method 'round'")
print("  Expect errors:")
passed = test_extendBounds([-11.1,-132.1],'round',expected=None) and passed
print("  Expect success:")
passed = test_extendBounds([11.1,132.1],'round',expected=[11,133]) and passed
passed = test_extendBounds([10,130],'round',expected=[10,130]) and passed
passed = test_extendBounds([11.1,132.1],'round',-2,expected=[11,133]) and passed
passed = test_extendBounds([-11.1,132.1],'round',expected=[-12,133]) and passed
passed = test_extendBounds([-1100.1,-132.1],'round',expected=[-1101,-132]) and passed
print()

print("Testing array execution")
print("  Expect errors:")
print("  Expect success:")
print("    method")
passed = test_extendBound(1.5,'up','nearestLeadingDigit',4,expected=2) and passed
passed = test_extendBound(1.5,'up','nearestPower',4,expected=4) and passed
passed = test_extendBound(1.5,'up','nearestMultiple',4,expected=4) and passed
passed = test_extendBound(1.5,'up','round',4,expected=2) and passed
print("    direction")
passed = test_extendBound(1.5,'up','nearestMultiple',4,expected=4) and passed
passed = test_extendBound(1.5,0.1,'nearestMultiple',4,expected=4) and passed
passed = test_extendBound(1.5,'down','nearestMultiple',4,expected=0) and passed
passed = test_extendBound(1.5,-0.1,'nearestMultiple',4,expected=0) and passed
print("    broadcasting")
passed = test_extendBound([1.5],'up','nearestMultiple',1,expected=[2]) and passed
passed = test_extendBound([1.5,2.5],'up','nearestMultiple',1,expected=[2,3]) and passed
passed = test_extendBound(1.5,['up','down'],'nearestMultiple',1,expected=[2,1]) and passed
passed = test_extendBound(1.5,'up',['nearestLeadingDigit','nearestPower','nearestMultiple','round'],3,expected=[2,3,3,2]) and passed
passed = test_extendBound(1.5,'up','nearestMultiple',[1,5,10],expected=[2,5,10]) and passed
print()

print("All tests passed:",passed)


# -

# ## Import using GeoPandas
#
# Earlier version

def GMLtoGDF(filename):
    gdf = geopandas.read_file(filename)
    gdf.rename_geometry('Geometry', inplace=True) # Default geometry column name is 'geometry'; changed for consistent capitalization of columns
    gdf.set_geometry('Geometry') # Renaming is insufficient; this sets special variable gdf.geometry = gdf['Geometry']
    gdf = gdf.set_crs(epsg=3347) # Needed only for FSA file, the others are 3347 and parsed correctly by geopandas, and the pdf in the zip file has the same projection parameters (FSA vs. DA, ADA, CT)
    gdf['Area'] = gdf['Geometry'].to_crs(epsg=6931).area # Equal-area projection # MODIFY THIS to account for validity regions of each geometry
    gdf['Centroid'] = gdf['Geometry'].centroid
    gdf['Geometry'] = gdf['Geometry'].to_crs(epsg=4326) # Latitude/Longitude representation
    gdf['Centroid'] = gdf['Centroid'].to_crs(epsg=4326) # Only the set geometry is converted with gdf.to_crs(); all other geometry-containing columns must be converted explicitly; here we convert all columns explicitly
    gdf = gdf.set_crs(epsg=4326) # The series and geodataframe can have separate crs; this was found necessary for the geopandas.union function to operate easily
    gdf['Centroid Latitude'] = gdf['Centroid'].geometry.y
    gdf['Centroid Longitude'] = gdf['Centroid'].geometry.x
    gdf.drop(columns = 'Centroid', inplace=True) # Because WKT Point cannot be serialized to JSON, we drop the Centroid column and keep only its float components
    return gdf


# Modified with comments; use standard 'geometry' instead of 'Geometry', preserve original crs

def GMLtoGDF(filename):
    gdf = geopandas.read_file(filename)
    #gdf.rename_geometry('Geometry', inplace=True) # Removed to revert to standard geopandas naming 'geometry' # Default geometry column name is 'geometry'; changed for consistent capitalization of columns
    #gdf.set_geometry('Geometry') # Removed to revert to standard geopandas naming 'geometry' # Renaming is insufficient; this sets special variable gdf.geometry = gdf['Geometry']
    gdf = gdf.set_crs(epsg=3347) # Needed only for FSA file, the others are 3347 and parsed correctly by geopandas, and the pdf in the zip file has the same projection parameters (FSA vs. DA, ADA, CT)
    gdf['Area'] = gdf.geometry.to_crs(epsg=6931).area # Equal-area projection # MODIFY THIS to account for validity regions of each geometry
    gdf['Centroid'] = gdf.geometry.centroid
    #gdf['Geometry'] = gdf.geometry.to_crs(epsg=4326) # Removed to preserve original crs # Latitude/Longitude representation
    gdf['Centroid'] = gdf['Centroid'].to_crs(epsg=4326) # Only the set geometry is converted with gdf.to_crs(); all other geometry-containing columns must be converted explicitly; here we convert all columns explicitly
    #gdf = gdf.set_crs(epsg=4326) # Removed to preserve original crs # The series and geodataframe can have separate crs; this was found necessary for the geopandas.union function to operate easily
    gdf['Centroid Latitude'] = gdf['Centroid'].geometry.y
    gdf['Centroid Longitude'] = gdf['Centroid'].geometry.x
    gdf.drop(columns = 'Centroid', inplace=True) # Because WKT Point cannot be serialized to JSON, we drop the Centroid column and keep only its float components
    return gdf


# Final version without comments

def GMLtoGDF(filename):
    gdf = geopandas.read_file(filename)
    gdf = gdf.set_crs(epsg=3347) # Needed only for FSA file, the others are 3347 and parsed correctly by geopandas, and the pdf in the zip file has the same projection parameters (FSA vs. DA, ADA, CT)
    gdf['Area'] = gdf.geometry.to_crs(epsg=6931).area # Equal-area projection # MODIFY THIS to account for validity regions of each geometry
    gdf['Centroid'] = gdf.geometry.centroid
    gdf['Centroid'] = gdf['Centroid'].to_crs(epsg=4326) # Only the set geometry is converted with gdf.to_crs(); all other geometry-containing columns must be converted explicitly; here we convert all columns explicitly
    gdf['Centroid Latitude'] = gdf['Centroid'].geometry.y
    gdf['Centroid Longitude'] = gdf['Centroid'].geometry.x
    gdf.drop(columns = 'Centroid', inplace=True) # Because WKT Point cannot be serialized to JSON, we drop the Centroid column and keep only its float components
    return gdf


# ## Save/Load (original)

# +
# Function(s) to encapsulate loading and saving long calculations
def loadResults_(name,tuples,fileformat='db',compress=False):
    '''Loads variables from files
    
    Parameters
    ----------
    name: str, file name base (including directory if desired)
    tuples: list of tuples (varname, suffix),
        varname: str, the key of the output dict where the data will be stored
        suffix: str, the string appended to name to generate a full file name
    fileformat: str, suffix to save the file with (do not include period)
    compress: bool, True to zip results (appends '.gz' to filename)
    
    Returns
    -------
    None if an error was encountered, or
    Tuple the length of tuples containing for each element of tuples:
        None if there was an error, or
        the variable loaded from file at the same position from tuples
    
    Notes
    -----
    Files read in binary format with optional gzip encoding
    This function is the complement to saveResults_()
    
    TODO
    ----
    Add option to change save format (text vs. binary)
    Make fileformat select the save format
    '''
    if type(name)!=str:
        print('Error: name must be a string')
        return None
    if type(fileformat)!=str:
        print('Error: fileformat must be a string')
        return None
    
    ret = []
    for n, s in tuples:
        fn = name+s+'.'+fileformat+('.gz' if compress else '')
        try:
            with open(fn,'rb') as file:
                ret.append(dill.loads(gzip.decompress(file.read()) if compress else file.read()))
        except (FileNotFoundError, IOError) as e:
            ret.append(None)
            print(f'An error was encountered while reading from file {fn}: {e}')
    return tuple(ret)

def loadResults(name):
    '''Loads variables 'gdf_union', 'times', and 'areas' from zipped files
    
    Parameters
    ----------
    name: str containing the base name of the files
    
    Returns
    -------
    None if an error was encountered, or
    Tuple the length of tuples containing:
        None if there was an error, or
        the variable loaded from file at the same position from tuples
    
    Notes
    -----
    File names area <name>_<variable>.db.gz and are in gzip dill binary format
    Uses outside variable DIR_RESULTS if available, otherwise put path in name
    '''
    tuples = [('gdf_union',''),
              ('times','_times'),
              ('areas','_areas')]
    
    return loadResults_(name,tuples,fileformat='db',compress=True)

def saveResults_(name,tuples,fileformat='db',compress=False):
    '''Saves variables to files
    
    Parameters
    ----------
    name: str, file name base (including directory if desired)
    tuples: list of tuples (varname, suffix),
        var: <any>, the variable to be output to file
        suffix: str, the string appended to name to generate a full file name
    fileformat: str, suffix to save the file with (do not include period)
    compress: bool, True to zip results (appends '.gz' to filename)
    
    Returns
    -------
    None if an error was encountered, or
    Tuple the same length as tuples containing return codes:
        0 Failure
        1 Success
    
    Notes
    -----
    Files written in binary format
    Files are created if they do not already exist
    Files are overwritten if they already exist

    TODO
    ----
    Make fileformat determine save format

    '''
    if type(name)!=str:
        print('Error: name must be a string')
        return None
    if type(fileformat)!=str:
        print('Error: fileformat must be a string')
        return None
    
    ret = []
    for v, s in tuples:
        fn = name+s+'.'+fileformat+('.gz' if compress else '')
        try:
            with open(fn,'wb+') as file:
                file.write(gzip.compress(dill.dumps(v)) if compress else dill.dumps(v))
                ret.append(1)
        except IOError as e:
            ret.append(0)
            print(f'An error was encountered while writing to file {fn}: {e}')
    return tuple(ret)

def saveResults(name, gdf_union, times, areas):
    '''Saves variables 'times', 'areas', and 'gdf_union' to zipped files
    
    Parameters
    ----------
    name: str, file name base (including directory if desired)
    gdf_union: geodataframe of geographic areas, produced from intersectGDF()
    times: 1d array of computation times, produced from intersectGDF()
    areas: list of lists of overlap areas, produced from intersectGDF()
    
    Returns
    -------
    None if an error was encountered, or
    Tuple the same length as tuples containing return codes:
        0 Failure
        1 Success
    
    Notes
    -----
    File names area <name><variable>.db.gz and are in gzip dill binary format
    Use outside variable DIR_RESULTS in construction of name
    '''
    tuples = [(gdf_union,''),
              (times,'_times'),
              (areas,'_areas')]
    
    return saveResults_(name,tuples,fileformat='db',compress=True)

def loadComputeSave(gdf1, key1, gdf2, key2, loadname=None, savename=None, verbosity=1, area_epsg=6931, gdf1b=None, gdf2b=None):
    '''Returns the overlap of geometries, defaulting to file versions if possible
    
    Parameters
    ----------
    gdf_1: GeoDataFrame (must match crs of gdf2, will be utilized for vectorized overlap calculation)
    keyfield1: column name in gdf1 which uniquely identifies each row and will be used to label the results
    gdf2: GeoDataFrame (must match crs of gdf1, will be iterated over for overlap calculation)
    keyfield2: column name in gdf2 which uniquely identifies each row and will be used to label the results
    loadname: str or None, base name of files to load data from (None -> 'DEFAULT'), see saveResults()
    savename: str or None, base name of files to save data to (None -> loadname), see loadResults()
    verbosity: int, detail level of reporting during execution: 0=none, 1=10-100 updates, 2=update every loop and announce exceptions
    area_epsg: int, convert to this epsg for area calculation
    gdf1b: gdf1 with all geometries valid, to be used in case of failed overlap with gdf1, if None use gdf1.buffer(0)
    gdf2b: gdf2 with all geometries valid, to be used in case of failed overlap with gdf2, if None use gdf2.buffer(0)

    
    Returns
    -------
    gdf_union: Geodataframe containing columns of nonzero overlap geometries, corresponding gdf1[keyfield1], and corresponding gdf2[keyfield2], where only one value of gdf1[keyfield1] is selected which is the one with maximum overlap area
    times: List of execution times for each overlap calculation; len(times)=gdf2.shape[0]
    areas: List of pandas Series of overlap areas; len(areas)=gdf2.shape[0], len(areas[i])=gdf1.shape[0]
    
    Notes
    -----
    gdf1 and gdf2 must be set to the same crs
    Iterates over gdf2, which should have the larger number of rows of {gdf1,gdf2} in order to minimize required memory (assuming geometries are of roughly equal size)
    '''
    verbosity = 1
    
    if savename is None:
        savename = loadname if not loadname is None else 'DEFAULT'
    
    ret = None if loadname is None else loadResults(DIR_RESULTS+loadname)
    recompute = False
    saveresults = False
    if ret is None:
        recompute = True
        saveresults = True
    else:
        gdf_union, times, areas = ret
        if gdf_union is None:
            if areas is None: # Recompute
                recompute = True
                saveresults = True
            else:                # Reconstruct from areas
                print("Overlaps will be recomputed based on loaded variable 'areas'")
                gdf_union, times, areas = intersectGDF(gdf1,key1,gdf2,key2,areas_in=temp_areas,verbosity=1)
                saveresults = True
        else:
            print("Overlaps loaded from file")

    if recompute:
        print("Overlaps must be computed")
        gdf_union, times, areas = intersectGDF(gdf1,key1,gdf2,key2,verbosity=1)
    
    if saveresults:
        saveResults(DIR_RESULTS+savename, gdf_union, times, areas)
        print("Variables saved to file at "+DIR_RESULTS+savename)

    return gdf_union, times, areas
# -

# ### Save/Load new



# ## mapCitiesAdjacent Old:

def mapCitiesAdjacent(cities,propertyname,title='',tooltiplabels=None):
    '''Displays cities in separate adjacent maps
    
    Cities dict as above
    Displayed as choropleth keyed to propertyname
    Header displays title and colormap labeled with keyed propertyname
    Tooltips pop up on mouseover showing properties listed in tooltiplabels
    
    Inspired by https://gitter.im/python-visualization/folium?at=5a36090a03838b2f2a04649d
    
    Assumes propertyname is identical in gdf and geojson
    '''
    if (not type(cities)==list) and type(cities)==dict:
        cities = [cities]
    
    f = bre.Figure()
    div_header = bre.Div(position='absolute',height='10%',width='100%',left='0%',top='0%').add_to(f)

    map_header = folium.Map(location=[0,0],control_scale=False,zoom_control=False,tiles=None,attr=False).add_to(div_header)
    div_header2 = bre.Div(position='absolute',height='10%',width='97%',left='3%',top='0%').add_to(div_header)
    html_header = '''<h3 align="left" style="font-size:16px;charset=utf-8"><b>{}</b></h3>'''.format(title)
    div_header2.get_root().html.add_child(folium.Element(html_header))

    vbounds = getCityBounds(cities,propertyname)
    vbounds[0] = 0
    vbounds = extendBounds(vbounds,'nearestLeadingDigit')
    
    cm_header = LinearColormap(
        colors=['yellow', 'orange', 'red'],
        index=None,
        vmin=vbounds[0],
        vmax=vbounds[1],
        caption=propertyname
        ).add_to(map_header) # .to_step(method='log', n=5, data=?), log has log labels but linear color scale

    for i, city in enumerate(cities):
        div_map = bre.Div(position='absolute',height='80%',width=f'{(100/len(cities))}%',left=f'{(i*100/len(cities))}%',top='10%').add_to(f)

        city_map = folium.Map(location=city['centroid'], control_scale=True)
        div_map.add_child(city_map)
        title_html = '''<h3 align="center" style="font-size:16px;charset=utf-8"><b>{}</b></h3>'''.format(city['name'])
        city_map.get_root().html.add_child(folium.Element(title_html))

        city_map.fit_bounds(city['bounds'])

        m = folium.GeoJson(
                    city['geojson'],
                    style_function=lambda feature: {
                        'fillColor': cm_header.rgb_hex_str(feature['properties'][propertyname]),
                        'fillOpacity': 0.8,
                        'color':'black',
                        'weight': 1,
                        'opacity': 0.2,
                    },
                    name=f'Choropleth_{i}'
                ).add_to(city_map)
        if not tooltiplabels==None:
            m.add_child(folium.features.GeoJsonTooltip(tooltiplabels))
    
    return display(f)









scalerN = sklearn.preprocessing.Normalizer()
scalerN.fit_transform(cities[0]['df_venues'])

ind_max = df_tmp['Dwelling Occupancy'].sort_values(ascending=False).index
val_max = df_tmp.loc[ind_max,'Dwelling Occupancy'].values
df_tmp.loc[ind_max[0],'Dwelling Occupancy'] = val_max[1]

# Venues scaled globally with additional factor, Census data and Venue density scaled;
from sklearn.preprocessing import StandardScaler
scale_column_names = ['Dwelling Density','Log10 (Population Density)','Log10 (1 - Residentiality)','Log10 (Venue Density)']
for i in range(max(df_venues['Category Depth'])+1):
    features = locals()[f"toronto_grouped_{i}"]
    vectors = features.drop('Postal Code',1).values
    vectors = scaler.fit_transform(vectors.reshape(-1,1)).reshape(*list(vectors.shape))
    vectors = vectors*1.4 # This factor increases MSS by about the same amount as the census and venue density does
    features[features.columns.drop('Postal Code')] = vectors
    features = features.merge(df_final[['Postal Code','Dwelling Density','Log10 (Population Density)','Log10 (1 - Residentiality)']],on='Postal Code').merge(df_venue_density[['Postal Code','Log10 (Venue Density)']],on='Postal Code')
    add_dict = {'name':f"Super-Scaled Venues Depth {i}\n with scaled Census Data\n and Venue Density",
                'labels':features['Postal Code'],
                'features':features.drop('Postal Code',1)
               }
    add_dict['features'][scale_column_names] = scaler.fit_transform(add_dict['features'][scale_column_names])
    k_means_data.append(add_dict)

from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import RobustScaler
from sklearn.preprocessing import QuantileTransformer
scaler = MinMaxScaler() # StandardScaler(), RobustScaler(), MinMaxScaler(), QuantileTransformer(output_distribution='normal', random_state=0)

# +
m = folium.Map(location=totalCentroid(df_merged), control_scale=True)
vbounds = [min(df_merged['Area Ratio']), max(df_merged['Area Ratio'])]

cm = LinearColormap(
                    colors=['blue', 'white', 'yellow', 'orange', 'red'],
                    index=[vbounds[0],0,-vbounds[0],vbounds[1]/2,vbounds[1]],
                    vmin=vbounds[0],
                    vmax=vbounds[1],
                    caption='Fractional Area Change'
                    ).add_to(m) # .to_step(method='log', n=5, data=?), log has log labels but linear color scale

g = folium.GeoJson(
            df_merged.to_crs(epsg=4326).to_json(),
            style_function=lambda feature: {
                'fillColor': cm.rgb_hex_str(feature['properties']['Area Ratio']) if abs(feature['properties']['Area Ratio'])>0.01 else 'white',# 'red' if feature['properties']['Area Diff'] else 'yellow',
                'fillOpacity': 0.4,
                'color':'black',
                'weight': 1,
                'opacity': 0.2,
            },
            name=f'Choropleth_{i}'
        ).add_to(m)
g.add_child(folium.features.GeoJsonTooltip(['Calc Census Area','Census Area','Area Diff','Area Ratio','CTUID'],localize=True))
display(m)
# -

# To simplify the data collection, let's construct a single dataframe with all CTs labeled by their city:

gdf_cities = pd.concat([c['gdf'].assign(**{'City':c['name']}) for c in cities],axis=0)
print(gdf_cities.shape)
gdf_cities.head(2)

# Old NN Distances

# %%time
for city in cities:
    city['gdf']['City'] = city['name'] # For definitive sorting after dataframe concatenation later
    city['gdf'][['NN Index','NN CTUID','NN Distance']] = None
    city['gdf'][['NN Index','NN CTUID','NN Distance']] = city['gdf'][['NN Index','NN CTUID','NN Distance']].astype(object)
    for ii, row in city['gdf'].iterrows():
        df_tmp = city['gdf'][['CTUID','Centroid Latitude','Centroid Longitude']].copy()
        df_tmp['Distance'] = [geopy.distance.distance((row['Centroid Latitude'],row['Centroid Longitude']),(x['Centroid Latitude'],x['Centroid Longitude'])).meters for jj, x in df_tmp.iterrows()]
        df_tmp.sort_values('Distance',inplace=True)
        # df_tmp.drop(index=row.name,axis=0,inplace=True) # To drop self-distance of 0
        city['gdf']['NN Index'].loc[ii] = df_tmp.index.values
        city['gdf']['NN CTUID'].loc[ii] = df_tmp['CTUID'].index.values
        city['gdf']['NN Distance'].loc[ii] = df_tmp['Distance'].index.values

#     df_vensum:   A pandas dataframe, collects summary of venue information, keyed to DAUID
#

for city in cities:
    city['isComplete'], city['df_venues'] = resumableGetNearbyVenues(city['gdf'], None, radius=RADIUS, extend=True)
    print(city['name'],'venue lookup completed successfully' if city['isComplete'] else 'venue lookup interrupted')
    display(city['df_venues'].tail(2))

for city in cities:
    if not city['isComplete']:
        city['isComplete'], city['df_venues'] = resumableGetNearbyVenues(city['gdf'], city['df_venues'], radius=RADIUS, extend=True)
        print(city['name'],'venue lookup completed successfully' if isComplete else 'venue lookup interrupted')
        display(df_venues.tail(2))
    else:
        print(city['name'],'venue lookup is already finished')

for jj in range(len(categories_depth_list)):
    fig, axs = plt.subplots(1,len(cities), figsize=(20,4))
    df_venues = pd.DataFrame()
    for ii, city in enumerate(cities):
        colname = f"Venue Category, Depth {jj}"
        df_venues = df_venues.append(city['df_venues'][['CTUID',colname]])
    df_category_counts_all = df_venues[['CTUID',colname]].groupby(colname).count()
    df_category_counts_all.rename(columns={'CTUID':'Unique Categories'},inplace=True)
    for ii, city in enumerate(cities):
        colname = f"Venue Category, Depth {jj}"
        df_category_counts = city['df_venues'][['CTUID',colname]].groupby(colname).count()
        df_category_counts.rename(columns={'CTUID':'Unique Categories'},inplace=True)
        plt.sca(axs[ii])
        df_category_counts_all.hist(log=True,bins=np.logspace(0,5,20),ax=plt.gca(),color='orange')
        df_category_counts.hist(log=True,bins=np.logspace(0,5,20),ax=plt.gca(),color='blue',alpha=0.5)
        plt.gca().set_xscale('log')
        plt.title(f"{city['name']}, Depth {jj}")
        plt.xlabel('Number of venues in a venue category')
        plt.ylabel('Count of venue categories');

# Reset defaults after using yellowbrick
# matplotlib.rcParams.update(matplotlib.rcParamsDefault)

matplotlib.rcParams.update(matplotlib.rcParamsDefault)

matplotlib.rcParams.update(matplotlib.rcParamsDefault)



# Note: yellowbrick changed matplotlib.rcParams and matplotlib.rcParamsDefault
#
# Before yellowbrick: matplotlib.rcParams['figure.figsize']=[6.0,4.0]
#
# And something set dpi to 72 for saving figures of gap statistic; default is 100.  Could it be folium or import order?
#
# matplotlib.pyplot.rcParams==matplotlib.rcParams
#
# For now the only solution is to not do plots after yellowbrick, else restart kernel (inelegant).

# + jupyter={"outputs_hidden": true} tags=[]
MyrcParams = matplotlib.rcParams.copy()
MyrcParams
# -

matplotlib.rcParams.update(MyrcParams)

matplotlib.rcParams.update(matplotlib.pyplot.rcParams)

# + tags=[]
matplotlib.rcParamsDefault

# + tags=[]
matplotlib.pyplot.rcParamsDefault
# -

matplotlib.pyplot.rcParams==matplotlib.rcParams

# +
import mpld3
#mpld3.enable_notebook()
import mpl_toolkits

fig = plt.figure()
ax = mpl_toolkits.mplot3d.Axes3D(fig)
#colors = [matplotlib.colors.to_hex(matplotlib.cm.hsv(x/(df_labels.shape[0]))) for x in range(df_labels.shape[0])]
ax.scatter(xs=df_features.iloc[:,0],ys=df_features.iloc[:,1],zs=df_features.iloc[:,2],alpha=0.2)#,c=clusters)
mpld3.display(fig)
# -

import ipyvolume as ipv

# +
fig = ipv.figure()
#scatter = ipv.scatter(df_features.iloc[:,0].values, df_features.iloc[:,1].values, df_features.iloc[:,2].values)
#ipv.show()

ipv.quickscatter(df_features.iloc[:,0], df_features.iloc[:,1], df_features.iloc[:,2], size=1, marker="sphere")
# -





# ## Appendix

# + [markdown] tags=[]
# DONE: TODO: For overlaps, store all thresholded overlap areas in a geodataframe, do area extraction afterward?
#
# DONE: TODO: Convert overlapGDF, overlapGDFareas, and subsequent area examination to units of meters instead of degrees (e.g. add to_crs())
#
# DONE: (supersceded) TODO: In the FSA geometry and area section, simplify code to use the newly declared indices instead of repeating computation.
#
# DONE: TODO: Additional overlap calculations use dataframes to simplify index tracking
#
# DONE: TODO: Indexing of overlap etc. should go into a dataframe, not be dealt with by indices...
#         Separate summarization in dataframe completed, still need to convert analysis sections
#
# TODO: NN Distances should be from centroid projected to other geometry exterior?
#
# TODO: Refine Foursquare lookup so that it doesn't require editing the functions to adapt to a new dataset
#
# TODO: For finding venues, also return a summary of counts vs radius steps for visualization
#
# DONE: TODO: Add comparison with Gaussian Mixture Model, Spectral Clustering, and other Estimation-Maximization clustering approaches like Maximum A-Posteriori Dirichlet Process Mixtures (MAP-DP, which can concurrently optimize for k). (See Extensions section, though Spectral Clustering was left out)
#
# TODO: Correct LoadResults_ docstring: input tuple - name not used, returns list
#
# DONE: TODO: Make df_features contain the unnormalized data, and df_features_scaled contain the adjusted data (reversed: df_features_unscaled instead)
#
# TODO: Functions for map image display/generate/save: abstract to use any map function e.g. with dict mapargs
#
# TODO: Functions for map feature display: a) look into DualMap to sync layer controls, b) input functions to control scaling e.g. log or feature arithmetic
#
# DONE: TODO: Next time, just copy params into data dicts?  Check if this would be too unreasonable to write around for plotting... (Done BGM section in Extensions)
#
# TODO: Consider if PCA clusters could originate as an artifact somehow, possibly due to quantization of the venue vector.
#
# TODO: Add a cluster
# -

import sys
modules = list(set(sys.modules) & set(globals()))
print(modules)

for module_name in modules:
    module = sys.modules[module_name]
    print(module_name, getattr(module, '__version__', 'unknown'))

# +
# Realative directories
DIR_DATA = 'Data/'
DIR_RESULTS = 'Results/'

import config # login credentials (for Here (geocoding) and Foursquare (venue lookup))

import importlib # for reloading modules (e.g. after editing config)
import copy # for performing a deepcopy to preserve original variables
import inspect # for determining function signatures and default arguments
import os # for getting current directory for relative save/load, and filepath functions
import sys # for checking system recursion limit
import time # for timing and web query rate throttling
import dill # for saving/loading results to save on computation time
import gzip # for compressing results to fit on github
import collections # for testing if parameters are list-like
import requests # for web queries (Foursquare)

import line_profiler # line profiling  # https://mortada.net/easily-profile-python-code-in-jupyter.html
# %load_ext line_profiler

import pandas # data management tool including dataframes and series, built on numpy
import pandas as pd

import geopandas # for creating dataframes with geometries, built on pandas, shapely, fiona, and matplotlib

import numpy # for high performance array mathematics, used by pandas
import numpy as np

import geopy # for geocoding and geodesic distance calculation
from geopy import distance
import pyproj # for more efficient coordinate transformations and geodesic area
import shapely # lower-level geometry library used in pyproj, geopy, and geopandas; for specific geometry error identification
import ogr # lower-level geospatial library, paired with gdal wrapper library

import matplotlib # for plotting
import matplotlib.cm # colormaps
import matplotlib.cm as cm
import matplotlib.colors
import matplotlib.pyplot as plt # for convenient plotting functions
# %matplotlib inline

import folium # for map display
import branca # for HTML capability for folium extension
import branca.element as bre
from branca.colormap import LinearColormap

import selenium # for rendering folium HTML maps and saving as images
from selenium import webdriver # also requires geckodriver in the environment; available in anaconda

import sklearn # machine learning module
import sklearn.cluster
import sklearn.feature_selection # for p-value calculation
import sklearn.metrics # for silhouette scoring
import sklearn.mixture # for gaussian mixture models

import yellowbrick # for silhouette visualization
import yellowbrick.cluster

import seaborn
import seaborn as sns # for convenient plotting

import scipy # for dendrogram display

import plotly
import plotly.express as px # for interactive 3D scatter plot

# Suppress unsightly warningsthat occur during geometric calculation
import warnings
warnings.filterwarnings("ignore")
import logging
logging.getLogger('shapely.geos').setLevel(logging.CRITICAL)

# +
import pkg_resources
import types
def get_imports():
    for name, val in globals().items():
        if isinstance(val, types.ModuleType):
            # Split ensures you get root package, 
            # not just imported function
            name = val.__name__.split(".")[0]

        elif isinstance(val, type):
            name = val.__module__.split(".")[0]

        # Some packages are weird and have different
        # imported names vs. system/pip names. Unfortunately,
        # there is no systematic way to get pip names from
        # a package's imported name. You'll have to add
        # exceptions to this list manually!
        poorly_named_packages = {
            "PIL": "Pillow",
            "sklearn": "scikit-learn"
        }
        if name in poorly_named_packages.keys():
            name = poorly_named_packages[name]

        yield name
imports = list(set(get_imports()))

# The only way I found to get the version of the root package
# from only the name of the package is to cross-check the names 
# of installed packages vs. imported packages
requirements = []
for m in pkg_resources.working_set:
    if m.project_name in imports and m.project_name!="pip":
        requirements.append((m.project_name, m.version))
# -

[x.project_name for x in pkg_resources.working_set]





# +
# From https://stackoverflow.com/questions/2572582/return-a-list-of-imported-python-modules-used-in-a-script
import ast

modules = set()

def visit_Import(node):
    for name in node.names:
        modules.add(name.name.split(".")[0])

def visit_ImportFrom(node):
    # if node.module is missing it's a "from . import ..." statement
    # if level > 0 it's a "from .submodule import ..." statement
    if node.module is not None and node.level == 0:
        modules.add(node.module.split(".")[0])

node_iter = ast.NodeVisitor()
node_iter.visit_Import = visit_Import
node_iter.visit_ImportFrom = visit_ImportFrom
# -

with open("Appendix.ipynb") as f:
    fstr = f"{chr(10)}#------------------------------------------{chr(10)}".join([''.join([s for s in cell['source'] if not s[0] in {'%','!'}]) for cell in json.loads(f.read())['cells'] if cell['cell_type']=='code'])
    node_iter.visit(ast.parse(fstr))
print(modules)

print([globals()[m] for m in modules if m in globals()])

# ! jupyter notebook --version

# ! jupyter nbconvert --version

__file__

print(os.getcwd())

import pipdeptree
# deps = !pipdeptree
print(chr(10).join([s for s in deps[:-1] if not s[0] in {' ','-','*'}]))

import pipdeptree

pipdeptree.main()

# deps = !pipdeptree

deps = deps[:-1]



lens = [len(s) for s in deps]

inds = np.array(lens)==0

import gdal

sum(inds)

[s for s in deps if not s[0] in {' ','-','*'}]

# !python --version

folium.__version__



import collections


def isIterable(obj):
    '''Determines if an object is list-like'''
    if not isinstance(obj, (str, collections.abc.ByteString)):
        if isinstance(obj, collections.abc.Sequence):
            return True
        try:
            a = [x for x in obj] # Caution: this could potentially consume significant resources
            if len(a)<1:
                return False
        except TypeError:
            return False
        return True
    return False


isIterable([])

len([])



StarcoderdataJupyterScriptsDedupFiltered