{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import copy\nimport json\nimport re\nimport time\nimport unicodedata\nimport warnings\nfrom functools import partial\n\nfrom tqdm.autonotebook import tqdm\nfrom urllib.parse import urlencode\nimport contextily\nimport folium\nimport geopandas as gpd\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pyproj\nimport rasterio\nimport rasterio.warp\nimport rasterio.features\nimport rasterio.plot\nimport requests\nimport requests_cache\nimport seaborn as sns\nimport shapely.wkt\n\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"code","source":"mpl.rcParams['figure.dpi'] = 200\npd.set_option('display.float_format', lambda x: f'{x:,.3f}')\npd.set_option(\"display.max_rows\", 100)\npd.set_option('display.width', 2000)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Disclosing\ndf_disclosing = pd.read_csv('../input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/2020_Cities_Disclosing_to_CDP.csv')\ndf_disclosing = df_disclosing.drop(['Access', 'Last update', 'Year Reported to CDP'], axis=1)\ndf_disclosing = df_disclosing.drop(['Reporting Authority', 'First Time Discloser'], axis=1)\n\n# Use None instead of NaN for strings.\ndf_disclosing['City'] = df_disclosing['City'].fillna('')\ndf_disclosing['city_location_raw'] = df_disclosing['City Location'].fillna('').apply(lambda x: shapely.wkt.loads(x) if x else None)\ndf_disclosing['Organization'] = df_disclosing['Organization'].fillna('')\n\n# Fix points.\ndf_disclosing.loc[df_disclosing['Account Number'] == 68383, 'City Location'] = 'POINT(-46.5 -23)'  # Missing -ve sign.\ndf_disclosing.loc[df_disclosing['Account Number'] == 826212, 'City Location'] = 'POINT(-103.6167 18.820)'  # No location given.\ndf_disclosing.loc[df_disclosing['Account Number'] == 826380, 'City Location'] = 'POINT(-104.361 19.7713)'  # No location given.\n\n# Assign point geometry.\ndf_disclosing['City Location'] = df_disclosing['City Location'].fillna('').apply(lambda x: shapely.wkt.loads(x) if x else None)\ndf_disclosing = gpd.GeoDataFrame(df_disclosing, geometry='City Location', crs='EPSG:4326')\n\n# Tidy population year.\nfixed_population = {19: 2019, 214: 2014, 216: 2016, 217: 2017}\ndf_disclosing['Population Year'] = df_disclosing['Population Year'].apply(lambda x: fixed_population.get(x, x))\n\n# Fix incorrect populations.\ndf_disclosing.loc[(df_disclosing['City'] == 'Montevideo') & (df_disclosing['Population'] == 1_383_432_134), 'Population'] = 1_383_432\ndf_disclosing.loc[(df_disclosing['Account Number'] == 826211) & (df_disclosing['Population'] == 399_724_000), 'Population'] = 399_724\ndf_disclosing.loc[(df_disclosing['Account Number'] == 54652) & (df_disclosing['Population'] == 676), 'Population'] = 676_000\ndf_disclosing.loc[(df_disclosing['Account Number'] == 841003) & (df_disclosing['Population'] == 5_247), 'Population'] = 524_700\ndf_disclosing.loc[(df_disclosing['Account Number'] == 839980) & (df_disclosing['Population'] == 26_492), 'Population'] = 264_920\n\ndf_disclosing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def log_print(*args, **kwargs):\n    # print(*args, **kwargs)\n    return None\n          \n\ndef osm_query(**kwargs):\n    \n    params = copy.deepcopy(kwargs)\n    params.update({\n        'format': 'jsonv2',\n        'polygon_geojson': '1',\n        'email': 'andrew@ajnisbet.com',\n    })\n    url = 'https://nominatim.openstreetmap.org/search'\n    \n    log_print('    https://nominatim.openstreetmap.org/ui/search.html?' + urlencode(kwargs))\n\n    response = requests.get(url, params=params)\n    df_response = pd.DataFrame.from_records(response.json())\n    if len(df_response) == 0:\n        return None\n        \n    def parse_geojson(x):\n        try:\n            return shapely.geometry.shape(x)\n        except AttributeError:\n            return None\n            \n    df_response['geometry'] = df_response.geojson.apply(parse_geojson)\n    df_response = df_response.drop([x for x in ['licence', 'boundingbox', 'place_id', 'osm_id', 'icon'] if x in df_response.columns], axis=1)\n    df_response['q'] = json.dumps(kwargs)\n    df_response['q_rank'] = df_response.index.values\n    df_response = pd.concat([df_response[df_response.category == 'boundary'], df_response[df_response.category != 'boundary']])\n    \n    return df_response\n    \n    \n    \ndef validate_osm_response(df, disclosure_point=None, country=None):\n    \n    # If there's no dataframe, nothing we can do about that.\n    if df is None or len(df) == 0:\n        log_print('    Empty response.')\n        return None\n    \n    # Try only to get polygons.\n    df = df[df.category == 'boundary']\n    if len(df) == 0:\n        log_print('    No polygons in response.')\n        return None\n    \n    # If country is provided, try to validate that.\n    if country:\n        display_name = df.display_name.iloc[0]\n        df = df[df.display_name.str.lower().str.contains(country.lower(), regex=False)]\n        if len(df) == 0:\n            log_print(f'    Country mismatch {country} vs {display_name}.')\n            return None\n    \n    # Check that polygon contains point.\n    if not disclosure_point:\n        log_print('    Success! No point to validate.')\n        return df.iloc[0]\n    \n    \n    df = gpd.GeoDataFrame(df, geometry='geometry', crs='EPSG:4326')\n    df['contains_point'] = df.intersects(disclosure_point)\n    df = df[df.contains_point]\n    if len(df) == 0:\n        log_print('    Point not in any polygon.')\n        return None\n    else:\n        log_print('    Success! Point in polygon.')\n        return df.iloc[0]\n    \n    \ndef parse_city_from_org(org):\n    org = org.split('Municipal Government of ')[-1].strip()\n    org = org.replace('Junta Municipal de Medio Ambiente ', '').strip()\n    org = org.replace('Junta Intermunicipal de Medio Ambiente ', '').strip()\n    org = org.replace(\"Municipal People's Government\", '').strip()\n    org = org.replace('County Council', '').strip()\n    org = org.replace('City Council', '').strip()\n    org = org.replace('Municipalidad Distrital de ', '').strip()\n    org = org.replace('Municipalidad de ', '').strip()\n    org = org.replace('Alcaldía Municipal de ', '').strip()\n    org = org.replace('Ayuntamiento de ', '').strip()\n    org = org.replace('Prefeitura de ', '').strip()\n    org = org.replace('Municipalidad de Provincial de ', '').strip()\n    org = org.replace('Junta Intermunicipal de la ', '').strip()\n    org = org.replace('Town of ', '').strip()\n    org = org.replace('City Government', '').strip()\n    org = org.replace('Gobierno Municipal de', '').strip()\n    org = org.replace('City of ', '').strip()\n    org = org.replace('Concejo Municipal de Distrito de ', '').strip()\n    org = org.replace('Municipalidad Provincial de ', '').strip()\n    org = org.replace('Municipal Government', '').strip()\n    org = org.replace('Municipio de ', '').strip()\n    org = org.replace('Municipality of', '').strip()\n    org = org.replace(' Council', '').strip()\n    org = org.replace('Gemeente', '').strip()\n    org = org.split('Estado de')[-1].strip()\n    org = org.split('(')[0].strip()\n    return org\n\ndef cdp_to_osm_country(c):\n    return {\n        'Taiwan, Greater China': 'Taiwan',\n        'United Kingdom of Great Britain and Northern Ireland': 'United Kingdom',\n        'China, Hong Kong Special Administrative Region': 'Hong Kong',\n        'Republic of Korea': 'South Korea',\n        'Bolivia (Plurinational State of)': 'Bolivia',\n        'Brazil': 'Brasil',\n        'Mexico': 'México',\n    }.get(c, c)\n    \ndef load_boundary_from_osm(i):\n    city_name = df_disclosing.City.iloc[i]\n    country_name = df_disclosing.Country.iloc[i]\n    osm_country = cdp_to_osm_country(country_name)\n    organization_name = df_disclosing.Organization.iloc[i]\n    disclosure_point = df_disclosing.geometry.iloc[i]\n    city_from_org = parse_city_from_org(organization_name)\n    city_filled = city_name or city_from_org\n    \n    dfs = []\n    \n    log_print()\n    log_print(f'[{i}] {city_name or None} ({organization_name or None}), {country_name} [{\"point\" if disclosure_point else \"no_point\"}]')\n    \n    if country_name in {'United States of America', 'Canada'} and city_name and re.search(', [A-Z]{2}$', city_name):\n        log_print('Checking structured City, State')\n        city, state = city_name.rsplit(',', 1)\n        df_osm = osm_query(city=city, state=state, country=country_name)\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point)\n        if val is not None:\n            return val\n        \n    # County and state provided.   \n    if country_name == 'United States of America' and organization_name and re.search('County, [A-Z]{2}$', organization_name):\n        log_print('Checking structured US County, State')\n        res = re.match('(.*?) County, ([A-Z]{2})$', organization_name)\n        county, state = res.groups()\n        df_osm = osm_query(county=county, state=state, country=country_name)\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point)\n        if val is not None:\n            return val\n        \n    # Just county.\n    if country_name == 'United States of America' and organization_name and organization_name.endswith(' County'):\n        log_print('Checking structured US County')\n        county = organization_name.split(' County')[0]\n        county = county.split('of')[-1].strip()\n        df_osm = osm_query(county=county, country=country_name)\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point)\n        if val is not None:\n            return val\n        \n    # Structured city, county.\n    if city_name:\n        log_print('Checking structured City, County')\n        df_osm = osm_query(city=city_name, country=country_name)\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point)\n        if val is not None:\n            return val\n        \n    # Check base organisation name with point safety.\n    if organization_name and disclosure_point:\n        log_print('Checking origanization only with point')\n        df_osm = osm_query(q=organization_name)\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point)\n        if val is not None:\n            return val\n        \n    # Check city with point safety.\n    if city_filled and disclosure_point:\n        log_print('Checking city only with point')\n        df_osm = osm_query(city=city_filled)\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point)\n        if val is not None:\n            return val\n        \n    # Filled city, structured.\n    if city_filled:\n        log_print('Structured City, Country')\n        df_osm = osm_query(city=city_filled, country=country_name)\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point)\n        if val is not None:\n            return val\n        \n    # Filled city, unstructured.\n    if city_filled:\n        log_print('Unstructured City, Country')\n        df_osm = osm_query(q=f'{city_filled}, {osm_country}')\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point)\n        if val is not None:\n            return val\n        \n    # City only, check country.\n    if city_filled:\n        log_print('City only, check country')\n        df_osm = osm_query(q=f'{city_filled}')\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point, country=country_name)\n        if val is not None:\n            return val\n        \n    # OSM country.\n    if city_filled:\n        log_print('City only, check country')\n        df_osm = osm_query(q=f'{city_filled}')\n        dfs.append(df_osm)\n        val = validate_osm_response(df_osm, disclosure_point=disclosure_point, country=osm_country)\n        if val is not None:\n            return val\n    \n    # If OSM didn't find anything, return the disclosure point.\n    dfs = [x for x in dfs if x is not None]\n    if len(dfs) > 0:\n        df = dfs[0]\n        return df.iloc[0]\n    elif disclosure_point:\n        print(f'{i}: returning dp')\n        return disclosure_point\n    else:\n        raise ValueError('No polygon found')\n    \n\n    \ndef buffer_meters(geom, buffer):\n    local_azimuthal_projection = f\"+proj=aeqd +R=6371000 +units=m +lat_0={geom.centroid.y} +lon_0={geom.centroid.x}\"\n\n    wgs84_to_aeqd = partial(\n        pyproj.transform,\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n        pyproj.Proj(local_azimuthal_projection),\n    )\n\n    aeqd_to_wgs84 = partial(\n        pyproj.transform,\n        pyproj.Proj(local_azimuthal_projection),\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n    )\n\n    geom_transformed = shapely.ops.transform(wgs84_to_aeqd, geom)\n\n    buffer = geom_transformed.buffer(buffer)\n\n    buffer_wgs84 = shapely.ops.transform(aeqd_to_wgs84, buffer)\n    \n    return buffer_wgs84        \n        \n        \n    \n    \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"geoms = []\ngeom_sources = []\n\nfor i in tqdm(range(len(df_disclosing))):\n    geom = load_boundary_from_osm(i)\n    raw_geom = geom\n    geom_source = 'CDP point'\n    \n    if isinstance(geom, pd.Series) or isinstance(geom, gpd.GeoSeries):\n        geom = geom.geometry\n        geom_source = 'OSM polygon'\n        \n    if isinstance(geom, shapely.geometry.MultiPolygon) or isinstance(geom, shapely.geometry.LineString):\n        geom = geom.convex_hull\n        geom_source = 'OSM polygon'\n\n    if isinstance(geom, shapely.geometry.Point):\n        if isinstance(raw_geom, shapely.geometry.Point):\n            geom_source = 'CDP point'\n        else:\n            geom_source = 'OSM point'\n            \n        # Exact radius doesn't matter so much, it will be adjusted to match population.\n        geom = buffer_meters(geom, 1000)\n        \n    if geom is None:\n        raise ValueError('No boundary found.')\n        \n    geoms.append(geom)\n    geom_sources.append(geom_source)\n\ndf_disclosing['city_geom'] = gpd.GeoSeries(geoms)\ndf_disclosing['city_geom_source'] = geom_sources\ndf_disclosing = df_disclosing.set_geometry('city_geom')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def point_to_latlon(point):\n    return point.y, point.x\n\ndef point_to_lonlat(point):\n    return point.x, point.y\n\ndef point_to_xy(point):\n    return point_to_lonlat(point)\n\ndef total_population_from_geom(f, a, window, g):\n    mask = rasterio.features.geometry_mask([g], a.shape, f.window_transform(window), invert=True).astype(int)\n    a_geom = a * mask\n    return a_geom.sum()\n\ndef shrink_around_point(point, poly, meters_to_shrink):\n    local_azimuthal_projection = f\"+proj=aeqd +R=6371000 +units=m +lat_0={point.y} +lon_0={point.x}\"\n\n    wgs84_to_aeqd = partial(\n        pyproj.transform,\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n        pyproj.Proj(local_azimuthal_projection),\n    )\n\n    poly_t = shapely.ops.transform(wgs84_to_aeqd, poly)\n    point_t = shapely.ops.transform(wgs84_to_aeqd, point)\n\n    max_radius = poly_t.hausdorff_distance(point_t)\n    new_radius = max_radius - meters_to_shrink\n    circle_to_cut = buffer_meters(point, new_radius)\n\n    new_poly = poly.intersection(circle_to_cut)\n    if not isinstance(new_poly, shapely.geometry.Polygon):\n        new_poly = new_poly.convex_hull\n    \n    return new_poly\n\n\ndef max_distance(point, poly):\n    local_azimuthal_projection = f\"+proj=aeqd +R=6371000 +units=m +lat_0={point.y} +lon_0={point.x}\"\n\n    wgs84_to_aeqd = partial(\n        pyproj.transform,\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n        pyproj.Proj(local_azimuthal_projection),\n    )\n\n    poly_t = shapely.ops.transform(wgs84_to_aeqd, poly)\n    point_t = shapely.ops.transform(wgs84_to_aeqd, point)\n\n    return poly_t.hausdorff_distance(point_t)\n\nPOP_PATH_1k = '../input/ppp_2020_1km_Aggregated.tif'\nPOP_PATH_100m = '../input/MOSAIC_ppp_prj_2020/MOSAIC_ppp_prj_2020.vrt'\ndef load_pop_raster(geom, target_population, buffer=10_000, pop_path=POP_PATH_1k):\n\n    \n\n    # Load population raster.\n    geom_buffer = buffer_meters(geom, buffer)\n    with rasterio.open(pop_path) as f:\n        window = rasterio.windows.from_bounds(*geom_buffer.bounds, f.transform)\n        a = f.read(1, window=window, masked=True)\n        a = np.ma.filled(a, 0)\n        \n    # Check population in raster is enough.\n    \n    raster_population = a.sum()\n    geom_population = total_population_from_geom(f, a, window, geom)\n    if raster_population < target_population*2:\n        if target_population > geom_population:\n            if (target_population - geom_population) / target_population > 0.25:\n                print(f'Target population {target_population} more than half raster {raster_population}')\n                return load_pop_raster(geom, target_population, buffer=buffer*3, pop_path=pop_path)\n\n    return f, a, window\n\ndef shrink_with_point(f, a, window, geom, point, target_population):\n    min_to_remove = 0\n    max_to_remove = int(max_distance(point, geom)) - 100\n\n    dx = (max_to_remove - min_to_remove) / 2 + min_to_remove\n    x_test = max_to_remove - dx\n\n    for i in range(10):\n        x_test = int(x_test)\n        geom_test = shrink_around_point(point, geom, x_test)\n        pop_test = total_population_from_geom(f, a, window, geom_test)\n        dx /= 2\n        if pop_test < target_population:\n            x_test -= dx\n        else:\n            x_test += dx\n    return geom_test, pop_test\n\ndef shrink_without_point(f, a, window, geom, target_population):\n    min_to_remove = 0\n    max_to_remove = int(max_distance(geom.centroid, geom)) - 100\n\n    dx = (max_to_remove - min_to_remove) / 2 + min_to_remove\n    x_test = max_to_remove - dx\n\n    for i in range(10):\n        x_test = int(x_test)\n        geom_test = buffer_meters(geom, -x_test)\n        pop_test = total_population_from_geom(f, a, window, geom_test)\n        dx /= 2\n        if pop_test < target_population:\n            x_test -= dx\n        else:\n            x_test += dx\n    return geom_test, pop_test\n\ndef grow_without_point(f, a, window, geom, target_population):\n    max_radius = max_distance(geom.centroid, geom)\n    to_add = max_radius / 10\n    \n    geom_test = buffer_meters(geom, to_add)\n    while total_population_from_geom(f, a, window, geom_test) < target_population:\n        to_add *= 2\n        geom_test = buffer_meters(geom, to_add)\n        \n    min_to_add = to_add / 2\n    max_to_add = to_add\n    dx = (max_to_add - min_to_add) / 2 + min_to_add\n    x_test = max_to_add - dx\n\n    for i in range(10):\n        x_test = int(x_test)\n        geom_test = buffer_meters(geom, x_test)\n        pop_test = total_population_from_geom(f, a, window, geom_test)\n        dx /= 2\n        if pop_test < target_population:\n            x_test += dx\n        else:\n            x_test -= dx\n    return geom_test, pop_test\n        \n    \ndef point_to_polygon_meters(point, poly):\n    if not point:\n        return None\n    \n    if point.intersects(poly):\n        return 0\n    \n    local_azimuthal_projection = f\"+proj=aeqd +R=6371000 +units=m +lat_0={point.y} +lon_0={point.x}\"\n\n    wgs84_to_aeqd = partial(\n        pyproj.transform,\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n        pyproj.Proj(local_azimuthal_projection),\n    )\n\n    poly_t = shapely.ops.transform(wgs84_to_aeqd, poly)\n    point_t = shapely.ops.transform(wgs84_to_aeqd, point)\n    \n    return poly.exterior.distance(point)\n    \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Rescale each geometry to match population.\nnew_geoms = []\nraw_geom_pops = []\n\nfor i in tqdm(range(len(df_disclosing))):\n    print()\n    print(f'i: {i}')\n    # Get data from df.\n    geom = df_disclosing.city_geom.iloc[i]\n    point = df_disclosing['City Location'].iloc[i]\n    if point and not geom.intersects(point):\n        point = None\n    target_population = df_disclosing.Population.iloc[i]\n    \n    if not target_population or pd.isnull(target_population):\n        new_geoms.append(geom)\n        raw_geom_pops.append(None)\n        continue\n\n    \n    # Load raster.\n    f, a, window = load_pop_raster(geom, target_population)\n    if max(a.shape) < 100:\n        f, a, window = load_pop_raster(geom, target_population, pop_path=POP_PATH_100m)\n    geom_population = total_population_from_geom(f, a, window, geom)\n    raw_geom_pops.append(geom_population)\n    if abs(geom_population - target_population) / target_population < 0.25:\n        new_geoms.append(geom)\n        continue\n    \n    # Do the transformation.\n    if geom_population <= target_population:\n        resized_geom, resized_population = grow_without_point(f, a, window, geom, target_population)\n    elif geom_population > target_population and point:\n        resized_geom, resized_population = shrink_with_point(f, a, window, geom, point, target_population)\n    elif geom_population > target_population and not point:\n        resized_geom, resized_population = shrink_with_point(f, a, window, geom, geom.centroid, target_population)\n\n    new_geoms.append(resized_geom)\n    del a","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_disclosing['city_geom_resized'] = new_geoms\ndf_disclosing = df_disclosing.set_geometry('city_geom_resized', crs='epsg:4326')\ndf_disclosing['area_km2'] = df_disclosing.to_crs('epsg:6933').area / 1000**2\ndf_disclosing['city_location_error'] = [point_to_polygon_meters(point, poly) for point, poly in zip(df_disclosing['city_location_raw'], new_geoms)]\ndf_disclosing['city_name_from_org'] = df_disclosing.Organization.fillna('').apply(parse_city_from_org)\ndf_disclosing['raw_geom_population'] = raw_geom_pops","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Save dataset.\ndf_dataset = df_disclosing[[\n    'Account Number',\n    'Population',\n    'city_geom_source',\n    'area_km2',\n    'city_geom_resized',\n    'city_location_error',\n    'city_name_from_org',\n    'raw_geom_population',\n]]\ndf_dataset = df_dataset.rename({'Account Number': 'account_number'}, axis=1)\ndf_dataset = df_dataset.rename({'Population': 'population_corrected'}, axis=1)\ndf_dataset.to_file('../data/city_boundaries_100m.geojson', driver='GeoJSON')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot a single result.\n\ni = 4\n\nprint()\nprint(f'i: {i}')\n# Get data from df.\ngeom = df_disclosing.city_geom.iloc[i]\npoint = df_disclosing['City Location'].iloc[i]\nif point and not geom.intersects(point):\n    point = None\ntarget_population = df_disclosing.Population.iloc[i]\n\nif not target_population or pd.isnull(target_population):\n    pass\n\n\n# Load raster.\nf, a, window = load_pop_raster(geom, target_population)\nif max(a.shape) < 100:\n    f, a, window = load_pop_raster(geom, target_population, pop_path=POP_PATH_100m)\ngeom_population = total_population_from_geom(f, a, window, geom)\nraw_geom_pops.append(geom_population)\n# if abs(geom_population - target_population) / target_population < 0.25:\n#     new_geoms.append(geom)\n#     continue\n\n# Do the transformation.\nif geom_population <= target_population:\n    resized_geom, resized_population = grow_without_point(f, a, window, geom, target_population)\nelif geom_population > target_population and point:\n    resized_geom, resized_population = shrink_with_point(f, a, window, geom, point, target_population)\nelif geom_population > target_population and not point:\n    resized_geom, resized_population = shrink_with_point(f, a, window, geom, geom.centroid, target_population)\n\n# new_geoms.append(resized_geom)\n# del a\n\n\n# Plot results.\nfig, axes = plt.subplots(1, 3, figsize=(20, 8))\nax = axes[0]\na_re, t_re = rasterio.warp.reproject(\n    source=a,\n    src_transform=f.window_transform(window),\n    src_crs='epsg:4326',\n    dst_crs='epsg:3857',\n)\nrasterio.plot.show(np.log10(a_re + 0.0001), transform=t_re, ax=ax)\nax.set_title(f'Bounds population: {int(a.sum()):,}')\nax.axis('off')\n\nax = axes[1]\ngpd.GeoSeries([geom], crs='EPSG:4326').to_crs('epsg:3857').plot(ax=ax, alpha=0.5, edgecolor='k')\nif point:\n    gpd.GeoSeries([point], crs='EPSG:4326').to_crs('epsg:3857').plot(ax=ax, color='red', alpha=0.5, edgecolor='k', markersize=200)\nax.set_title(f'Starting population: {int(geom_population):,}')\naxes[1].set_xlim(axes[0].get_xlim())\naxes[1].set_ylim(axes[0].get_ylim())\ncontextily.add_basemap(ax, source=contextily.providers.OpenStreetMap.Mapnik)\nax.axis('off')\n\nax = axes[2]\ngpd.GeoSeries([resized_geom], crs='EPSG:4326').to_crs('epsg:3857').plot(ax=ax, alpha=0.5, edgecolor='k')\nif point:\n    gpd.GeoSeries([point], crs='EPSG:4326').to_crs('epsg:3857').plot(ax=ax, color='red', alpha=0.5, edgecolor='k', markersize=200)\naxes[2].set_xlim(axes[0].get_xlim())\naxes[2].set_ylim(axes[0].get_ylim())\ncontextily.add_basemap(ax, source=contextily.providers.OpenStreetMap.Mapnik)\nax.set_title(f'Final population: {int(resized_population):,} (target {int(target_population):,})')\nax.axis('off')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}