From ca2596a99e0798a845ca354a38ffbaac56973b2f Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Thu, 24 Oct 2024 14:36:57 -0400 Subject: [PATCH 01/25] initial commit --- city_metrix/layers/__init__.py | 1 + city_metrix/layers/isochrone.py | 13 +++++++ .../percent_population_within_isochrone.py | 38 +++++++++++++++++++ tests/test_layers.py | 7 ++++ 4 files changed, 59 insertions(+) create mode 100644 city_metrix/layers/isochrone.py create mode 100644 city_metrix/metrics/percent_population_within_isochrone.py diff --git a/city_metrix/layers/__init__.py b/city_metrix/layers/__init__.py index e701e52..1c42593 100644 --- a/city_metrix/layers/__init__.py +++ b/city_metrix/layers/__init__.py @@ -21,3 +21,4 @@ from .nasa_dem import NasaDEM from .era_5_hottest_day import Era5HottestDay from .impervious_surface import ImperviousSurface +from .isochrone import Isochrone diff --git a/city_metrix/layers/isochrone.py b/city_metrix/layers/isochrone.py new file mode 100644 index 0000000..9b3d8a3 --- /dev/null +++ b/city_metrix/layers/isochrone.py @@ -0,0 +1,13 @@ +import geopandas as gpd +import shapely + +from .layer import Layer, get_utm_zone_epsg + +class Isochrone(Layer): + def __init__(self, geojson_filepath, **kwargs): + super().__init__(**kwargs) + gdf = gpd.read_file(geojson_filepath) + self.gdf = gpd.GeoDataFrame({'is_accessible': [1] * len(gdf), 'geometry': gdf.to_crs('EPSG:4326')['geometry']}).to_crs('EPSG:4326').set_crs('EPSG:4326').set_geometry('geometry') + + def get_data(self, bbox): + return self.gdf.clip(shapely.box(*bbox)) diff --git a/city_metrix/metrics/percent_population_within_isochrone.py b/city_metrix/metrics/percent_population_within_isochrone.py new file mode 100644 index 0000000..313c490 --- /dev/null +++ b/city_metrix/metrics/percent_population_within_isochrone.py @@ -0,0 +1,38 @@ +from geopandas import GeoDataFrame, GeoSeries +import numpy as np + +from city_metrix.layers.isochrone import Isochrone +from city_metrix.layers.world_pop import WorldPop + + +def get_accessible_population(access_features_layer, popraster_layer, zones): + if len(access_features_layer.gdf): + result_series = popraster_layer.mask(access_features_layer).groupby(zones).sum() + else: + result_series = pd.Series([0] * len(zones)) + return result_series + +def percent_population_within_isochrone(zones: GeoDataFrame, isochrone_filename, agesex_classes=[], worldpop_year=2020) -> GeoSeries: + population_layer = WorldPop(agesex_classes=agesex_classes, worldpop_year=worldpop_year) + accesszone_layer = Isochrone(isochrone_filename) + + try: + access_pop = get_accessible_population(accesszone_layer, population_layer, zones) + total_pop = population_layer.groupby(zones).sum() + result = (access_pop / total_pop) * 100 + + except: + # Sometimes doing entire zones gdf causes groupby to throw empty-GDF error -- workaraound is to go district-by-district + result_gdf = GeoDataFrame({'geometry': zones['geometry']}).set_geometry('geometry').set_crs('EPSG:4326') + access_fraction = [] + for idx in zones.index: + access_pop = get_accessible_population(accesszone_layer, population_layer, zones.loc[[zones.index[idx]]])[0] + total_pop = population_layer.groupby(zones.loc[[zones.index[idx]]]).sum()[0] + if total_pop != 0: + access_fraction.append(access_pop / total_pop) + else: + access_fraction.append(np.nan) + + result = access_fraction.replace([np.inf,], np.nan) * 100 + + return result \ No newline at end of file diff --git a/tests/test_layers.py b/tests/test_layers.py index f8e1fe2..6677252 100644 --- a/tests/test_layers.py +++ b/tests/test_layers.py @@ -10,6 +10,7 @@ EsaWorldCoverClass, HighLandSurfaceTemperature, ImperviousSurface, + Isochrone, LandsatCollection2, LandSurfaceTemperature, NasaDEM, @@ -66,6 +67,12 @@ def test_impervious_surface(): data = ImperviousSurface().get_data(BBOX) assert np.size(data) > 0 +def test_isochrone(): + layer = Isochrone('https://wri-cities-indicators.s3.us-east-1.amazonaws.com/data/isochrones/nairobi_schools.geojson') + nairobi_bbox = (36.66446402, -1.44560888, 37.10497899, -1.16058296) + data = layer.get_data(nairobi_bbox) + assert np.size(data) > 0 + def test_land_surface_temperature(): data = LandSurfaceTemperature().get_data(BBOX) assert np.size(data) > 0 From 03aa10cef742cd10e33cd42b6fc78d2652f39b61 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Wed, 20 Nov 2024 17:21:08 -0500 Subject: [PATCH 02/25] change isochrone to isoline, address filename naming convention --- city_metrix/layers/__init__.py | 2 +- city_metrix/layers/isochrone.py | 13 -------- city_metrix/layers/isoline.py | 32 +++++++++++++++++++ city_metrix/metrics/__init__.py | 1 + ...y => percent_population_within_isoline.py} | 13 +++++--- tests/test_layers.py | 6 ++-- 6 files changed, 46 insertions(+), 21 deletions(-) delete mode 100644 city_metrix/layers/isochrone.py create mode 100644 city_metrix/layers/isoline.py rename city_metrix/metrics/{percent_population_within_isochrone.py => percent_population_within_isoline.py} (66%) diff --git a/city_metrix/layers/__init__.py b/city_metrix/layers/__init__.py index 1c42593..694bd1c 100644 --- a/city_metrix/layers/__init__.py +++ b/city_metrix/layers/__init__.py @@ -21,4 +21,4 @@ from .nasa_dem import NasaDEM from .era_5_hottest_day import Era5HottestDay from .impervious_surface import ImperviousSurface -from .isochrone import Isochrone +from .isoline import Isoline diff --git a/city_metrix/layers/isochrone.py b/city_metrix/layers/isochrone.py deleted file mode 100644 index 9b3d8a3..0000000 --- a/city_metrix/layers/isochrone.py +++ /dev/null @@ -1,13 +0,0 @@ -import geopandas as gpd -import shapely - -from .layer import Layer, get_utm_zone_epsg - -class Isochrone(Layer): - def __init__(self, geojson_filepath, **kwargs): - super().__init__(**kwargs) - gdf = gpd.read_file(geojson_filepath) - self.gdf = gpd.GeoDataFrame({'is_accessible': [1] * len(gdf), 'geometry': gdf.to_crs('EPSG:4326')['geometry']}).to_crs('EPSG:4326').set_crs('EPSG:4326').set_geometry('geometry') - - def get_data(self, bbox): - return self.gdf.clip(shapely.box(*bbox)) diff --git a/city_metrix/layers/isoline.py b/city_metrix/layers/isoline.py new file mode 100644 index 0000000..f799b0c --- /dev/null +++ b/city_metrix/layers/isoline.py @@ -0,0 +1,32 @@ +import geopandas as gpd +import shapely +import boto3 + +from .layer import Layer, get_utm_zone_epsg + +BUCKET_NAME = 'wri-cities-indicators' +PREFIX = 'data/isolines/' + +class Isochrone(Layer): + def __init__(self, cityname, amenityname, travelmode, threshold_type, threshold_value, year=None, **kwargs): + super().__init__(**kwargs) + + # Get list of isoline files from S3 and find the most recent file with correct cityname, amenity, etc + s3 = boto3.client('s3') + obj_list = s3.list_objects(Bucket=BUCKET_NAME, Prefix=PREFIX)['Contents'] + objnames = [i['Key'].replace(PREFIX, '') for i in obj_list if len(i['Key'].replace(PREFIX, '')) > 0] + matches = [oname for oname in objnames if oname.split('_')[:-1] == [cityname, amenityname, travelmode, threshold_type, threshold_value]] + if year is not None: + matches = [oname for oname in matches if oname.split('.')[0].split('_')[-1][:4] == str(year)] + if len(matches) == 0: + raise Exception('No isoline file found.') + matches.sort(key=lambda x: int(x.split('.')[0].split('_')[-1])) + objname = matches[-1] + + # Get data from S3 + url = f'https://{BUCKET_NAME}.s3.us-east-1.amazonaws.com/{PREFIX}{objname}' + gdf = gpd.read_file(url) + self.gdf = gpd.GeoDataFrame({'is_accessible': [1] * len(gdf), 'geometry': gdf.to_crs('EPSG:4326')['geometry']}).to_crs('EPSG:4326').set_crs('EPSG:4326').set_geometry('geometry') + + def get_data(self, bbox): + return self.gdf.clip(shapely.box(*bbox)) diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index d95cfa5..f2a01e3 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -5,3 +5,4 @@ from .urban_open_space import urban_open_space from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing +from .percent_population_within_isoline import percent_population_within_isoline \ No newline at end of file diff --git a/city_metrix/metrics/percent_population_within_isochrone.py b/city_metrix/metrics/percent_population_within_isoline.py similarity index 66% rename from city_metrix/metrics/percent_population_within_isochrone.py rename to city_metrix/metrics/percent_population_within_isoline.py index 313c490..79dcf07 100644 --- a/city_metrix/metrics/percent_population_within_isochrone.py +++ b/city_metrix/metrics/percent_population_within_isoline.py @@ -1,7 +1,7 @@ from geopandas import GeoDataFrame, GeoSeries import numpy as np -from city_metrix.layers.isochrone import Isochrone +from city_metrix.layers.isochrone import Isoline from city_metrix.layers.world_pop import WorldPop @@ -12,9 +12,14 @@ def get_accessible_population(access_features_layer, popraster_layer, zones): result_series = pd.Series([0] * len(zones)) return result_series -def percent_population_within_isochrone(zones: GeoDataFrame, isochrone_filename, agesex_classes=[], worldpop_year=2020) -> GeoSeries: - population_layer = WorldPop(agesex_classes=agesex_classes, worldpop_year=worldpop_year) - accesszone_layer = Isochrone(isochrone_filename) +def percent_population_within_isoline(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], worldpop_year=2020) -> GeoSeries: + # cityname example: ARG-Buenos-Aires + # amenityname is OSMclass names, in lowercase + # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) + # threshold_type is distance or time + # threshold_value is integer, in minutes or meters + population_layer = WorldPop(agesex_classes=agesex_classes, worldpop_year=worldpop_year) + accesszone_layer = Isoline(cityname, amenityname, travelmode, threshold_type, threshold_value, year) try: access_pop = get_accessible_population(accesszone_layer, population_layer, zones) diff --git a/tests/test_layers.py b/tests/test_layers.py index 6677252..8ca615f 100644 --- a/tests/test_layers.py +++ b/tests/test_layers.py @@ -10,7 +10,7 @@ EsaWorldCoverClass, HighLandSurfaceTemperature, ImperviousSurface, - Isochrone, + Isoline, LandsatCollection2, LandSurfaceTemperature, NasaDEM, @@ -67,8 +67,8 @@ def test_impervious_surface(): data = ImperviousSurface().get_data(BBOX) assert np.size(data) > 0 -def test_isochrone(): - layer = Isochrone('https://wri-cities-indicators.s3.us-east-1.amazonaws.com/data/isochrones/nairobi_schools.geojson') +def test_isoline(): + layer = Isoline('KEN-Nairobi', 'schools', 'walk', 'time', '15') nairobi_bbox = (36.66446402, -1.44560888, 37.10497899, -1.16058296) data = layer.get_data(nairobi_bbox) assert np.size(data) > 0 From 1996fa6f23d13d202018e6e9d20af4532c572437 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Thu, 21 Nov 2024 21:24:39 -0500 Subject: [PATCH 03/25] interact with AWS --- city_metrix/layers/isoline.py | 44 ++++++++------ city_metrix/layers/layer.py | 1 + city_metrix/metrics/__init__.py | 2 +- .../metrics/percent_population_access.py | 60 +++++++++++++++++++ .../percent_population_within_isoline.py | 43 ------------- 5 files changed, 89 insertions(+), 61 deletions(-) create mode 100644 city_metrix/metrics/percent_population_access.py delete mode 100644 city_metrix/metrics/percent_population_within_isoline.py diff --git a/city_metrix/layers/isoline.py b/city_metrix/layers/isoline.py index f799b0c..523a489 100644 --- a/city_metrix/layers/isoline.py +++ b/city_metrix/layers/isoline.py @@ -7,24 +7,34 @@ BUCKET_NAME = 'wri-cities-indicators' PREFIX = 'data/isolines/' -class Isochrone(Layer): - def __init__(self, cityname, amenityname, travelmode, threshold_type, threshold_value, year=None, **kwargs): +class Isoline(Layer): + def __init__(self, params, aws_profilename=None, **kwargs): super().__init__(**kwargs) - - # Get list of isoline files from S3 and find the most recent file with correct cityname, amenity, etc - s3 = boto3.client('s3') - obj_list = s3.list_objects(Bucket=BUCKET_NAME, Prefix=PREFIX)['Contents'] - objnames = [i['Key'].replace(PREFIX, '') for i in obj_list if len(i['Key'].replace(PREFIX, '')) > 0] - matches = [oname for oname in objnames if oname.split('_')[:-1] == [cityname, amenityname, travelmode, threshold_type, threshold_value]] - if year is not None: - matches = [oname for oname in matches if oname.split('.')[0].split('_')[-1][:4] == str(year)] - if len(matches) == 0: - raise Exception('No isoline file found.') - matches.sort(key=lambda x: int(x.split('.')[0].split('_')[-1])) - objname = matches[-1] - - # Get data from S3 - url = f'https://{BUCKET_NAME}.s3.us-east-1.amazonaws.com/{PREFIX}{objname}' + # params is dict with keys cityname, amenityname, travelmode, threshold_type, threshold_value, year + cityname = params['cityname'] + amenityname = params['amenityname'] + travelmode = params['travelmode'] + threshold_type = params['threshold_type'] + threshold_value = params['threshold_value'] + year = params['year'] + # Get list of isoline files from S3 and find the most recent file with correct cityname, amenity, etc + if aws_profilename is None: + session = boto3.Session() + else: + session = boto3.Session(profile_name=aws_profilename) + s3 = session.client('s3') + obj_list = s3.list_objects(Bucket=BUCKET_NAME, Prefix=PREFIX)['Contents'] + objnames = [i['Key'].replace(PREFIX, '') for i in obj_list if len(i['Key'].replace(PREFIX, '')) > 0] + matches = [oname for oname in objnames if oname.split('.')[0].split('_')[:-1] == [cityname, amenityname, travelmode, threshold_type, str(threshold_value)]] + if year is not None: + matches = [oname for oname in matches if oname.split('.')[0].split('_')[-1][:4] == str(year)] + if len(matches) == 0: + raise Exception('No isoline file found.') + matches.sort(key=lambda x: int(x.split('.')[0].split('_')[-1])) + objname = matches[-1] + + # Get data from S3 + url = f'https://{BUCKET_NAME}.s3.us-east-1.amazonaws.com/{PREFIX}{objname}' gdf = gpd.read_file(url) self.gdf = gpd.GeoDataFrame({'is_accessible': [1] * len(gdf), 'geometry': gdf.to_crs('EPSG:4326')['geometry']}).to_crs('EPSG:4326').set_crs('EPSG:4326').set_geometry('geometry') diff --git a/city_metrix/layers/layer.py b/city_metrix/layers/layer.py index dce9603..7b3c42c 100644 --- a/city_metrix/layers/layer.py +++ b/city_metrix/layers/layer.py @@ -353,6 +353,7 @@ def write_layer(path, data): data.to_file(path, driver="GeoJSON") else: raise NotImplementedError("Can only write DataArray, Dataset, or GeoDataFrame") + raise NotImplementedError("Can only write DataArray, Dataset, or GeoDataFrame") def write_dataarray(path, data): # for rasters, need to write to locally first then copy to cloud storage diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index f2a01e3..3d23fe4 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -5,4 +5,4 @@ from .urban_open_space import urban_open_space from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing -from .percent_population_within_isoline import percent_population_within_isoline \ No newline at end of file +from .percent_population_access import percent_population_access \ No newline at end of file diff --git a/city_metrix/metrics/percent_population_access.py b/city_metrix/metrics/percent_population_access.py new file mode 100644 index 0000000..53a4872 --- /dev/null +++ b/city_metrix/metrics/percent_population_access.py @@ -0,0 +1,60 @@ +from geopandas import GeoDataFrame, GeoSeries +import numpy as np + +from city_metrix.layers.isoline import Isoline +from city_metrix.layers.world_pop import WorldPop + + + + +def percent_population_access(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], worldpop_year=2020, aws_profilename=None) -> GeoSeries: + + def get_accessible_population(access_features_layer, popraster_layer, zones): + if len(access_features_layer.gdf): + result_series = popraster_layer.mask(access_features_layer).groupby(zones).mean() * popraster_layer.mask(access_features_layer).groupby(zones).count() + else: + result_series = pd.Series([0] * len(zones)) + return result_series + + # cityname example: ARG-Buenos-Aires + # amenityname is OSMclass names, in lowercase + # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) + # threshold_type is distance or time + # threshold_value is integer, in minutes or meters + population_layer = WorldPop(agesex_classes=agesex_classes, year=worldpop_year) + params = { + 'cityname': cityname, + 'amenityname': amenityname, + 'travelmode': travelmode, + 'threshold_type': threshold_type, + 'threshold_value': threshold_value, + 'year': isoline_year + } + accesszone_layer = Isoline(params, aws_profilename=aws_profilename) + + try: + access_pop = get_accessible_population(accesszone_layer, population_layer, zones) + total_pop = population_layer.groupby(zones).mean() * population_layer.groupby(zones).count() + result = (access_pop / total_pop) * 100 + + except: + # Sometimes doing entire zones gdf causes groupby to throw empty-GDF error -- workaraound is to go district-by-district + print('Calculating district-by-district') + result_gdf = GeoDataFrame({'geometry': zones['geometry']}).set_geometry('geometry').set_crs('EPSG:4326') + access_fraction = [] + for idx in zones.index: + try: # Sometimes there's still an empty-gdf error + access_pop = get_accessible_population(accesszone_layer, population_layer, zones.loc[[zones.index[idx]]])[0] + total_pop = (population_layer.groupby(zones.loc[[zones.index[idx]]]).mean() * population_layer.groupby(zones.loc[[zones.index[idx]]]).count())[0] + if total_pop != 0: + access_fraction.append(access_pop / total_pop) + else: + access_fraction.append(np.nan) + except: + print('Empty-GDF error for index {0}'.format(idx)) + access_fraction.append(np.nan) + result_gdf['access_fraction'] = access_fraction + result_gdf['access_fraction'].replace([np.inf, -np.inf], np.nan, inplace=True) + result_gdf['access_fraction'] = result_gdf['access_fraction'] * 100 + + return result_gdf \ No newline at end of file diff --git a/city_metrix/metrics/percent_population_within_isoline.py b/city_metrix/metrics/percent_population_within_isoline.py deleted file mode 100644 index 79dcf07..0000000 --- a/city_metrix/metrics/percent_population_within_isoline.py +++ /dev/null @@ -1,43 +0,0 @@ -from geopandas import GeoDataFrame, GeoSeries -import numpy as np - -from city_metrix.layers.isochrone import Isoline -from city_metrix.layers.world_pop import WorldPop - - -def get_accessible_population(access_features_layer, popraster_layer, zones): - if len(access_features_layer.gdf): - result_series = popraster_layer.mask(access_features_layer).groupby(zones).sum() - else: - result_series = pd.Series([0] * len(zones)) - return result_series - -def percent_population_within_isoline(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], worldpop_year=2020) -> GeoSeries: - # cityname example: ARG-Buenos-Aires - # amenityname is OSMclass names, in lowercase - # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) - # threshold_type is distance or time - # threshold_value is integer, in minutes or meters - population_layer = WorldPop(agesex_classes=agesex_classes, worldpop_year=worldpop_year) - accesszone_layer = Isoline(cityname, amenityname, travelmode, threshold_type, threshold_value, year) - - try: - access_pop = get_accessible_population(accesszone_layer, population_layer, zones) - total_pop = population_layer.groupby(zones).sum() - result = (access_pop / total_pop) * 100 - - except: - # Sometimes doing entire zones gdf causes groupby to throw empty-GDF error -- workaraound is to go district-by-district - result_gdf = GeoDataFrame({'geometry': zones['geometry']}).set_geometry('geometry').set_crs('EPSG:4326') - access_fraction = [] - for idx in zones.index: - access_pop = get_accessible_population(accesszone_layer, population_layer, zones.loc[[zones.index[idx]]])[0] - total_pop = population_layer.groupby(zones.loc[[zones.index[idx]]]).sum()[0] - if total_pop != 0: - access_fraction.append(access_pop / total_pop) - else: - access_fraction.append(np.nan) - - result = access_fraction.replace([np.inf,], np.nan) * 100 - - return result \ No newline at end of file From 697c2bc4a9ebbf9a69a620f9fdbc4c0449a4d8d1 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Thu, 21 Nov 2024 21:39:44 -0500 Subject: [PATCH 04/25] fix bug when district-by-district not necessary --- city_metrix/metrics/percent_population_access.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/city_metrix/metrics/percent_population_access.py b/city_metrix/metrics/percent_population_access.py index 53a4872..edddd3f 100644 --- a/city_metrix/metrics/percent_population_access.py +++ b/city_metrix/metrics/percent_population_access.py @@ -32,15 +32,18 @@ def get_accessible_population(access_features_layer, popraster_layer, zones): } accesszone_layer = Isoline(params, aws_profilename=aws_profilename) + result_gdf = GeoDataFrame({'geometry': zones['geometry']}).set_geometry('geometry').set_crs('EPSG:4326') + try: access_pop = get_accessible_population(accesszone_layer, population_layer, zones) total_pop = population_layer.groupby(zones).mean() * population_layer.groupby(zones).count() result = (access_pop / total_pop) * 100 + result_gdf['access_fraction'] = result except: # Sometimes doing entire zones gdf causes groupby to throw empty-GDF error -- workaraound is to go district-by-district print('Calculating district-by-district') - result_gdf = GeoDataFrame({'geometry': zones['geometry']}).set_geometry('geometry').set_crs('EPSG:4326') + access_fraction = [] for idx in zones.index: try: # Sometimes there's still an empty-gdf error From 300f8668be2df6033266db07c6f2e7897142a28c Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Wed, 8 Jan 2025 11:03:07 -0500 Subject: [PATCH 05/25] fix test to pass dict instead of positional args --- tests/test_layers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_layers.py b/tests/test_layers.py index 8ca615f..3383d51 100644 --- a/tests/test_layers.py +++ b/tests/test_layers.py @@ -68,7 +68,7 @@ def test_impervious_surface(): assert np.size(data) > 0 def test_isoline(): - layer = Isoline('KEN-Nairobi', 'schools', 'walk', 'time', '15') + layer = Isoline({'cityname': 'KEN-Nairobi', 'amenityname': 'schools', 'travelmode': 'walk', 'threshold_type': 'time', 'threshold_value': '15', 'year': 2023}) nairobi_bbox = (36.66446402, -1.44560888, 37.10497899, -1.16058296) data = layer.get_data(nairobi_bbox) assert np.size(data) > 0 From 2ee7cb6cbf0014b01e1ec0e08f0d08397b576a6d Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Thu, 9 Jan 2025 14:24:43 -0500 Subject: [PATCH 06/25] initial commit --- .../layers/count_accessible_amenities.py | 42 +++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 city_metrix/layers/count_accessible_amenities.py diff --git a/city_metrix/layers/count_accessible_amenities.py b/city_metrix/layers/count_accessible_amenities.py new file mode 100644 index 0000000..523a489 --- /dev/null +++ b/city_metrix/layers/count_accessible_amenities.py @@ -0,0 +1,42 @@ +import geopandas as gpd +import shapely +import boto3 + +from .layer import Layer, get_utm_zone_epsg + +BUCKET_NAME = 'wri-cities-indicators' +PREFIX = 'data/isolines/' + +class Isoline(Layer): + def __init__(self, params, aws_profilename=None, **kwargs): + super().__init__(**kwargs) + # params is dict with keys cityname, amenityname, travelmode, threshold_type, threshold_value, year + cityname = params['cityname'] + amenityname = params['amenityname'] + travelmode = params['travelmode'] + threshold_type = params['threshold_type'] + threshold_value = params['threshold_value'] + year = params['year'] + # Get list of isoline files from S3 and find the most recent file with correct cityname, amenity, etc + if aws_profilename is None: + session = boto3.Session() + else: + session = boto3.Session(profile_name=aws_profilename) + s3 = session.client('s3') + obj_list = s3.list_objects(Bucket=BUCKET_NAME, Prefix=PREFIX)['Contents'] + objnames = [i['Key'].replace(PREFIX, '') for i in obj_list if len(i['Key'].replace(PREFIX, '')) > 0] + matches = [oname for oname in objnames if oname.split('.')[0].split('_')[:-1] == [cityname, amenityname, travelmode, threshold_type, str(threshold_value)]] + if year is not None: + matches = [oname for oname in matches if oname.split('.')[0].split('_')[-1][:4] == str(year)] + if len(matches) == 0: + raise Exception('No isoline file found.') + matches.sort(key=lambda x: int(x.split('.')[0].split('_')[-1])) + objname = matches[-1] + + # Get data from S3 + url = f'https://{BUCKET_NAME}.s3.us-east-1.amazonaws.com/{PREFIX}{objname}' + gdf = gpd.read_file(url) + self.gdf = gpd.GeoDataFrame({'is_accessible': [1] * len(gdf), 'geometry': gdf.to_crs('EPSG:4326')['geometry']}).to_crs('EPSG:4326').set_crs('EPSG:4326').set_geometry('geometry') + + def get_data(self, bbox): + return self.gdf.clip(shapely.box(*bbox)) From 690f0b136570dbb427e3f5c69e1b92123fcbf33d Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Thu, 9 Jan 2025 14:33:24 -0500 Subject: [PATCH 07/25] initial commit --- .../layers/count_accessible_amenities.py | 42 ----------- .../metrics/count_accessible_amenities.py | 69 +++++++++++++++++++ 2 files changed, 69 insertions(+), 42 deletions(-) delete mode 100644 city_metrix/layers/count_accessible_amenities.py create mode 100644 city_metrix/metrics/count_accessible_amenities.py diff --git a/city_metrix/layers/count_accessible_amenities.py b/city_metrix/layers/count_accessible_amenities.py deleted file mode 100644 index 523a489..0000000 --- a/city_metrix/layers/count_accessible_amenities.py +++ /dev/null @@ -1,42 +0,0 @@ -import geopandas as gpd -import shapely -import boto3 - -from .layer import Layer, get_utm_zone_epsg - -BUCKET_NAME = 'wri-cities-indicators' -PREFIX = 'data/isolines/' - -class Isoline(Layer): - def __init__(self, params, aws_profilename=None, **kwargs): - super().__init__(**kwargs) - # params is dict with keys cityname, amenityname, travelmode, threshold_type, threshold_value, year - cityname = params['cityname'] - amenityname = params['amenityname'] - travelmode = params['travelmode'] - threshold_type = params['threshold_type'] - threshold_value = params['threshold_value'] - year = params['year'] - # Get list of isoline files from S3 and find the most recent file with correct cityname, amenity, etc - if aws_profilename is None: - session = boto3.Session() - else: - session = boto3.Session(profile_name=aws_profilename) - s3 = session.client('s3') - obj_list = s3.list_objects(Bucket=BUCKET_NAME, Prefix=PREFIX)['Contents'] - objnames = [i['Key'].replace(PREFIX, '') for i in obj_list if len(i['Key'].replace(PREFIX, '')) > 0] - matches = [oname for oname in objnames if oname.split('.')[0].split('_')[:-1] == [cityname, amenityname, travelmode, threshold_type, str(threshold_value)]] - if year is not None: - matches = [oname for oname in matches if oname.split('.')[0].split('_')[-1][:4] == str(year)] - if len(matches) == 0: - raise Exception('No isoline file found.') - matches.sort(key=lambda x: int(x.split('.')[0].split('_')[-1])) - objname = matches[-1] - - # Get data from S3 - url = f'https://{BUCKET_NAME}.s3.us-east-1.amazonaws.com/{PREFIX}{objname}' - gdf = gpd.read_file(url) - self.gdf = gpd.GeoDataFrame({'is_accessible': [1] * len(gdf), 'geometry': gdf.to_crs('EPSG:4326')['geometry']}).to_crs('EPSG:4326').set_crs('EPSG:4326').set_geometry('geometry') - - def get_data(self, bbox): - return self.gdf.clip(shapely.box(*bbox)) diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py new file mode 100644 index 0000000..fb97d99 --- /dev/null +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -0,0 +1,69 @@ +from geopandas import GeoDataFrame, GeoSeries +import numpy as np + +from city_metrix.layers.isoline import Isoline +from city_metrix.layers.world_pop import WorldPop + +def count_accessible_amenities(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], worldpop_year=2020, aws_profilename=None) -> GeoSeries: + + class AccessCountTmp(Layer): + def __init__(self, access_gdf, **kwargs): + super().__init__(**kwargs) + self.gdf = access_gdf + def get_data(self, bbox): + return self.gdf.clip(shapely.box(*bbox)) + + def get_accessible_population(access_features_layer, popraster_layer, zones): + + # cityname example: ARG-Buenos-Aires + # amenityname is OSMclass names, in lowercase + # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) + # threshold_type is distance or time + # threshold_value is integer, in minutes or meters + population_layer = WorldPop(agesex_classes=agesex_classes, year=worldpop_year) + params = { + 'cityname': cityname, + 'amenityname': amenityname, + 'travelmode': travelmode, + 'threshold_type': threshold_type, + 'threshold_value': threshold_value, + 'year': isoline_year + } + iso_layer = Isoline(params, aws_profilename=aws_profilename) + iso_gdf = accesszone_layer.get_data(zones.total_bounds) + + # Collect individual boundary linestrings of each isoline polygon + iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] + iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] + iso_boundarylines = [] + for i in iso_boundarylinesmulti: + if type(i) == shapely.geometry.linestring.LineString: + iso_boundarylines.append(i) + else: + for j in list(i.geoms): + iso_boundarylines.append(j) + iso_bl_gdf = gpd.GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) + + # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" + iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) + accesscount_gdf = gpd.GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}) + + # For each dissected polygon, find how many of the original isoline polys contain the centroid + # This is the number of amenity points are accessible within the dissected polygon + count_amenities = iso_dissected_polys.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 + for iso_idx in range(1, len(iso_gdf)): + count_amenities = count_amenities + (iso_dissected_polys.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) + + accesscount_gdf['count_amenities'] = count_amenities + + count_layers = {count: AccessCountTmp(gpd.GeoDataFrame({'count_amenities': polys_gdf['count_amenities']==count, 'geometry': polys_gdf['geometry']})) for count in range(1, polys_gdf['count_amenities'].max() + 1)} + + # For each zone, find average number of accessible amenities, and store in result_gdf + max_count = polys_gdf['count_amenities'].max() + result = population_layer.mask(count_layers[1]).count() / max_count + for count in range(1, max_count+1): + result += population_layer.mask(count_layers[count]).count() / max_count + + result_gdf['count_accessible_amenities'] = result + + return result_gdf \ No newline at end of file From d68398823252c68ee0b20bc87f1623f016169b38 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Thu, 9 Jan 2025 18:33:27 -0500 Subject: [PATCH 08/25] fixes after initial testing --- .../metrics/count_accessible_amenities.py | 111 +++++++++--------- 1 file changed, 54 insertions(+), 57 deletions(-) diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index fb97d99..9081101 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -1,69 +1,66 @@ from geopandas import GeoDataFrame, GeoSeries -import numpy as np +from city_metrix.layers import Layer +import shapely from city_metrix.layers.isoline import Isoline from city_metrix.layers.world_pop import WorldPop -def count_accessible_amenities(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], worldpop_year=2020, aws_profilename=None) -> GeoSeries: - - class AccessCountTmp(Layer): - def __init__(self, access_gdf, **kwargs): - super().__init__(**kwargs) - self.gdf = access_gdf - def get_data(self, bbox): - return self.gdf.clip(shapely.box(*bbox)) +class AccessCountTmp(Layer): + def __init__(self, access_gdf, **kwargs): + super().__init__(**kwargs) + self.gdf = access_gdf + def get_data(self, bbox): + return self.gdf.clip(shapely.box(*bbox)) - def get_accessible_population(access_features_layer, popraster_layer, zones): - +def count_accessible_amenities(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, aws_profilename=None) -> GeoSeries: # cityname example: ARG-Buenos-Aires # amenityname is OSMclass names, in lowercase # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) # threshold_type is distance or time # threshold_value is integer, in minutes or meters - population_layer = WorldPop(agesex_classes=agesex_classes, year=worldpop_year) - params = { - 'cityname': cityname, - 'amenityname': amenityname, - 'travelmode': travelmode, - 'threshold_type': threshold_type, - 'threshold_value': threshold_value, - 'year': isoline_year - } - iso_layer = Isoline(params, aws_profilename=aws_profilename) - iso_gdf = accesszone_layer.get_data(zones.total_bounds) - - # Collect individual boundary linestrings of each isoline polygon - iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] - iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] - iso_boundarylines = [] - for i in iso_boundarylinesmulti: - if type(i) == shapely.geometry.linestring.LineString: - iso_boundarylines.append(i) - else: - for j in list(i.geoms): - iso_boundarylines.append(j) - iso_bl_gdf = gpd.GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) - - # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" - iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) - accesscount_gdf = gpd.GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}) - - # For each dissected polygon, find how many of the original isoline polys contain the centroid - # This is the number of amenity points are accessible within the dissected polygon - count_amenities = iso_dissected_polys.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 - for iso_idx in range(1, len(iso_gdf)): - count_amenities = count_amenities + (iso_dissected_polys.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) - - accesscount_gdf['count_amenities'] = count_amenities - - count_layers = {count: AccessCountTmp(gpd.GeoDataFrame({'count_amenities': polys_gdf['count_amenities']==count, 'geometry': polys_gdf['geometry']})) for count in range(1, polys_gdf['count_amenities'].max() + 1)} + population_layer = WorldPop(agesex_classes=[], year=2020) + params = { + 'cityname': cityname, + 'amenityname': amenityname, + 'travelmode': travelmode, + 'threshold_type': threshold_type, + 'threshold_value': threshold_value, + 'year': isoline_year + } + iso_layer = Isoline(params, aws_profilename=aws_profilename) + iso_gdf = iso_layer.get_data(zones.total_bounds) + + # Collect individual boundary linestrings of each isoline polygon + iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] + iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] + iso_boundarylines = [] + for i in iso_boundarylinesmulti: + if type(i) == shapely.geometry.linestring.LineString: + iso_boundarylines.append(i) + else: + for j in list(i.geoms): + iso_boundarylines.append(j) + iso_bl_gdf = GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) + + # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" + iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) + accesscount_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') + + # For each dissected polygon, find how many of the original isoline polys contain the centroid + # This is the number of amenity points are accessible within the dissected polygon + count_amenities = iso_dissected_polys.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 + for iso_idx in range(1, len(iso_gdf)): + count_amenities = count_amenities + (iso_dissected_polys.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) - # For each zone, find average number of accessible amenities, and store in result_gdf - max_count = polys_gdf['count_amenities'].max() - result = population_layer.mask(count_layers[1]).count() / max_count - for count in range(1, max_count+1): - result += population_layer.mask(count_layers[count]).count() / max_count - - result_gdf['count_accessible_amenities'] = result - - return result_gdf \ No newline at end of file + accesscount_gdf['count_amenities'] = count_amenities + + count_layers = {count: AccessCountTmp(GeoDataFrame({'count_amenities': accesscount_gdf['count_amenities']==count, 'geometry': accesscount_gdf['geometry']}).set_crs('EPSG:4326')) for count in range(1, accesscount_gdf['count_amenities'].max() + 1)} + + # For each zone, find average number of accessible amenities, and store in result_gdf + max_count = accesscount_gdf['count_amenities'].max() + result = population_layer.mask(count_layers[1]).groupby(zones).count() * 1 + for count in range(2, max_count+1): + result += population_layer.mask(count_layers[count]).groupby(zones).count() * count + result = result / population_layer.groupby(zones).count() + result_gdf = GeoDataFrame({'count_accessible_amenities': result, 'geometry': zones['geometry']}).set_crs('EPSG:4326') + return result_gdf \ No newline at end of file From 39711687bb93052024f630679b2d76a92425d933 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Fri, 10 Jan 2025 10:33:43 -0500 Subject: [PATCH 09/25] add tests --- city_metrix/metrics/__init__.py | 3 ++- tests/test_metrics.py | 15 +++++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index 3d23fe4..964896c 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -5,4 +5,5 @@ from .urban_open_space import urban_open_space from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing -from .percent_population_access import percent_population_access \ No newline at end of file +from .percent_population_access import percent_population_access +from .count_accessible_amenities import count_accessible_amenities \ No newline at end of file diff --git a/tests/test_metrics.py b/tests/test_metrics.py index 8fd42cc..2b02ced 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -3,6 +3,7 @@ import pytest + def test_built_land_with_high_lst(): indicator = built_land_with_high_land_surface_temperature(ZONES) expected_zone_size = ZONES.geometry.size @@ -23,6 +24,13 @@ def test_built_land_without_tree_cover(): actual_indicator_size = indicator.size assert expected_zone_size == actual_indicator_size +def test_count_accessible_amenities(): + from .conftest import create_fishnet_grid + NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296).reset_index() + indicator = count_accessible_amenities(NAIROBI_BBOX, 'KEN-Nairobi', 'schools', 'walk', 'time', '15', 2024, aws_profilename=None) + expected_zone_size = NAIROBI_BBOX.geometry.size + actual_indicator_size = indicator.size + assert expected_zone_size == actual_indicator_size @pytest.mark.skip(reason="CDS API needs personal access token file to run") def test_era_5_met_preprocess(): @@ -43,6 +51,13 @@ def test_natural_areas(): actual_indicator_size = indicator.size assert expected_zone_size == actual_indicator_size +def test_percent_population_access(): + from .conftest import create_fishnet_grid + NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296).reset_index() + indicator = percent_population_access(NAIROBI_BBOX, 'KEN-Nairobi', 'schools', 'walk', 'time', '15', 2024, aws_profilename=None) + expected_zone_size = NAIROBI_BBOX.geometry.size + actual_indicator_size = indicator.size + assert expected_zone_size == actual_indicator_size def test_urban_open_space(): indicator = urban_open_space(ZONES) From 7b4aace8bce1d77918d2c858a9f182aebf79d6cd Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Fri, 10 Jan 2025 16:37:19 -0500 Subject: [PATCH 10/25] fix bug in test --- tests/test_metrics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_metrics.py b/tests/test_metrics.py index 2b64c63..7f59135 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -26,7 +26,7 @@ def test_built_land_without_tree_cover(): def test_count_accessible_amenities(): from .conftest import create_fishnet_grid - NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296).reset_index() + NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296, 0.01).reset_index() indicator = count_accessible_amenities(NAIROBI_BBOX, 'KEN-Nairobi', 'schools', 'walk', 'time', '15', 2024, aws_profilename=None) expected_zone_size = NAIROBI_BBOX.geometry.size actual_indicator_size = indicator.size From c09fd1ab9ab7e33cdf2fd40595aa8524f5423235 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Fri, 10 Jan 2025 16:48:46 -0500 Subject: [PATCH 11/25] bug fix again --- tests/test_metrics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_metrics.py b/tests/test_metrics.py index 7f59135..fe714e6 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -53,7 +53,7 @@ def test_natural_areas(): def test_percent_population_access(): from .conftest import create_fishnet_grid - NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296).reset_index() + NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296, 0.01).reset_index() indicator = percent_population_access(NAIROBI_BBOX, 'KEN-Nairobi', 'schools', 'walk', 'time', '15', 2024, aws_profilename=None) expected_zone_size = NAIROBI_BBOX.geometry.size actual_indicator_size = indicator.size From e094141f5645ca54e5b4d8e077effc225968b332 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Sun, 12 Jan 2025 18:02:56 -0500 Subject: [PATCH 12/25] add demographic-specific metrics --- city_metrix/metrics/__init__.py | 2 +- .../metrics/percent_population_access.py | 35 ++++++++++++++++--- 2 files changed, 31 insertions(+), 6 deletions(-) diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index da75299..0c9524e 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -5,6 +5,6 @@ from .urban_open_space import urban_open_space from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing -from .percent_population_access import percent_population_access +from .percent_population_access import percent_population_access, percent_population_access_elderly, percent_population_access_children, percent_population_access_female, percent_population_access_informal from .count_accessible_amenities import count_accessible_amenities from .recreational_space_per_capita import recreational_space_per_capita diff --git a/city_metrix/metrics/percent_population_access.py b/city_metrix/metrics/percent_population_access.py index edddd3f..2e2e5bd 100644 --- a/city_metrix/metrics/percent_population_access.py +++ b/city_metrix/metrics/percent_population_access.py @@ -2,12 +2,17 @@ import numpy as np from city_metrix.layers.isoline import Isoline -from city_metrix.layers.world_pop import WorldPop +from city_metrix.layers.world_pop import WorldPop, UrbanLandUse +def percent_population_access(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], informal_only=False, worldpop_year=2020, aws_profilename=None) -> GeoSeries: - -def percent_population_access(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], worldpop_year=2020, aws_profilename=None) -> GeoSeries: + class InformalTemp(UrbanLandUse): + def get_data(self, bbox): + data = super().get_data(bbox) + result = xr.where(data==3, 1, np.nan).rio.write_crs(data.crs)) + result = result.assign_attrs(**data.attrs) + return result def get_accessible_population(access_features_layer, popraster_layer, zones): if len(access_features_layer.gdf): @@ -21,7 +26,11 @@ def get_accessible_population(access_features_layer, popraster_layer, zones): # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) # threshold_type is distance or time # threshold_value is integer, in minutes or meters - population_layer = WorldPop(agesex_classes=agesex_classes, year=worldpop_year) + if informal_only: + informal_mask = InformalTemp() + population_layer = WorldPop(agesex_classes=agesex_classes, year=worldpop_year, masks=[informal_layer,]) + else: + population_layer = WorldPop(agesex_classes=agesex_classes, year=worldpop_year) params = { 'cityname': cityname, 'amenityname': amenityname, @@ -60,4 +69,20 @@ def get_accessible_population(access_features_layer, popraster_layer, zones): result_gdf['access_fraction'].replace([np.inf, -np.inf], np.nan, inplace=True) result_gdf['access_fraction'] = result_gdf['access_fraction'] * 100 - return result_gdf \ No newline at end of file + return result_gdf + +def percent_population_access_elderly(zones, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, worldpop_year=2020, aws_profilename=None): + agesex_classes=['F_60', 'F_65', 'F_70', 'F_75', 'F_80', 'M_60', 'M_65', 'M_70', 'M_75', 'M_80'] + return percent_population_access(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, agesex_classes=agesex_classes, informal_only=False, worldpop_year=worldpop_year, aws_profilename=aws_profilename) + +def percent_population_access_children(zones, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, worldpop_year=2020, aws_profilename=None): + agesex_classes=['F_0', 'F_1', 'F_5', 'F_10', 'M_0', 'M_1', 'M_5', 'M_10'] + return percent_population_access(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, agesex_classes=agesex_classes, informal_only=False, worldpop_year=worldpop_year, aws_profilename=aws_profilename) + +def percent_population_access_female(zones, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, worldpop_year=2020, aws_profilename=None): + agesex_classes=['F_0', 'F_1', 'F_5', 'F_10', 'F_15', 'F_20', 'F_25', 'F_30', 'F_35', 'F_40', 'F_45', 'F_50', 'F_55', 'F_60', 'F_65', 'F_70', 'F_75', 'F_80'] + return percent_population_access(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, agesex_classes=agesex_classes, informal_only=False, worldpop_year=worldpop_year, aws_profilename=aws_profilename) + +def percent_population_access_informal(zones, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, worldpop_year=2020, aws_profilename=None): + agesex_classes=[] + return percent_population_access(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, agesex_classes=agesex_classes, informal_only=True, worldpop_year=worldpop_year, aws_profilename=aws_profilename) \ No newline at end of file From 3e2615b47f2e029b472234aaa88f862ebf3911e8 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Sun, 12 Jan 2025 18:12:28 -0500 Subject: [PATCH 13/25] test bypass --- tests/test_layers.py | 1 + tests/test_metrics.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/tests/test_layers.py b/tests/test_layers.py index 9256496..7b692da 100644 --- a/tests/test_layers.py +++ b/tests/test_layers.py @@ -78,6 +78,7 @@ def test_impervious_surface(): data = ImperviousSurface().get_data(BBOX) assert np.size(data) > 0 +@pytest.mark.skipif(EXECUTE_IGNORED_TESTS == False, reason="AWS redentials needed") def test_isoline(): layer = Isoline({'cityname': 'KEN-Nairobi', 'amenityname': 'schools', 'travelmode': 'walk', 'threshold_type': 'time', 'threshold_value': '15', 'year': 2023}) nairobi_bbox = (36.66446402, -1.44560888, 37.10497899, -1.16058296) diff --git a/tests/test_metrics.py b/tests/test_metrics.py index fe714e6..4cc9d72 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -24,6 +24,7 @@ def test_built_land_without_tree_cover(): actual_indicator_size = indicator.size assert expected_zone_size == actual_indicator_size +@pytest.mark.skipif(EXECUTE_IGNORED_TESTS == False, reason="AWS credentials needed") def test_count_accessible_amenities(): from .conftest import create_fishnet_grid NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296, 0.01).reset_index() @@ -51,6 +52,7 @@ def test_natural_areas(): actual_indicator_size = indicator.size assert expected_zone_size == actual_indicator_size +@pytest.mark.skipif(EXECUTE_IGNORED_TESTS == False, reason="AWS credentials needed") def test_percent_population_access(): from .conftest import create_fishnet_grid NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296, 0.01).reset_index() From df71d8203ca214447ca85861f450edf6eed98d03 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Sun, 12 Jan 2025 18:19:59 -0500 Subject: [PATCH 14/25] bug fix --- city_metrix/metrics/percent_population_access.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/city_metrix/metrics/percent_population_access.py b/city_metrix/metrics/percent_population_access.py index 2e2e5bd..4f04b79 100644 --- a/city_metrix/metrics/percent_population_access.py +++ b/city_metrix/metrics/percent_population_access.py @@ -10,7 +10,7 @@ def percent_population_access(zones: GeoDataFrame, cityname, amenityname, travel class InformalTemp(UrbanLandUse): def get_data(self, bbox): data = super().get_data(bbox) - result = xr.where(data==3, 1, np.nan).rio.write_crs(data.crs)) + result = xr.where(data==3, 1, np.nan).rio.write_crs(data.crs) result = result.assign_attrs(**data.attrs) return result From 96d556cdd1f209f37db3d390b7e4fcf69b89bb3b Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Sun, 12 Jan 2025 18:57:15 -0500 Subject: [PATCH 15/25] bug fix --- city_metrix/metrics/percent_population_access.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/city_metrix/metrics/percent_population_access.py b/city_metrix/metrics/percent_population_access.py index 4f04b79..3b0a713 100644 --- a/city_metrix/metrics/percent_population_access.py +++ b/city_metrix/metrics/percent_population_access.py @@ -2,7 +2,8 @@ import numpy as np from city_metrix.layers.isoline import Isoline -from city_metrix.layers.world_pop import WorldPop, UrbanLandUse +from city_metrix.layers.world_pop import WorldPop +from city_metrix.layers.urban_land_use import UrbanLandUse def percent_population_access(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], informal_only=False, worldpop_year=2020, aws_profilename=None) -> GeoSeries: From 58a64c8c9743360703bbe5bdcfaaad9d398a1114 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Mon, 13 Jan 2025 07:17:54 -0500 Subject: [PATCH 16/25] distinguish populationmean and areamean metrics --- city_metrix/metrics/__init__.py | 2 +- .../metrics/count_accessible_amenities.py | 22 ++++++++++++++----- tests/test_metrics.py | 2 +- 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index 0c9524e..968c4be 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -6,5 +6,5 @@ from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing from .percent_population_access import percent_population_access, percent_population_access_elderly, percent_population_access_children, percent_population_access_female, percent_population_access_informal -from .count_accessible_amenities import count_accessible_amenities +from .count_accessible_amenities import count_accessible_amenities_areamean, count_accessible_amenities_populationmean from .recreational_space_per_capita import recreational_space_per_capita diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index 9081101..0ce09c7 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -12,12 +12,13 @@ def __init__(self, access_gdf, **kwargs): def get_data(self, bbox): return self.gdf.clip(shapely.box(*bbox)) -def count_accessible_amenities(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, aws_profilename=None) -> GeoSeries: +def _count_accessible_amenities(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, weighting='population', aws_profilename=None) -> GeoSeries: # cityname example: ARG-Buenos-Aires # amenityname is OSMclass names, in lowercase # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) # threshold_type is distance or time # threshold_value is integer, in minutes or meters + population_layer = WorldPop(agesex_classes=[], year=2020) params = { 'cityname': cityname, @@ -58,9 +59,20 @@ def count_accessible_amenities(zones: GeoDataFrame, cityname, amenityname, trave # For each zone, find average number of accessible amenities, and store in result_gdf max_count = accesscount_gdf['count_amenities'].max() - result = population_layer.mask(count_layers[1]).groupby(zones).count() * 1 - for count in range(2, max_count+1): - result += population_layer.mask(count_layers[count]).groupby(zones).count() * count + if weighting == 'area': + result = population_layer.mask(count_layers[1]).groupby(zones).count() + for count in range(2, max_count+1): + result += population_layer.mask(count_layers[count]).groupby(zones).count() * count + else: # weighting == 'population' + result = population_layer.mask(count_layers[1]).groupby(zones).sum() + for count in range(2, max_count+1): + result += population_layer.mask(count_layers[count]).groupby(zones).sum() result = result / population_layer.groupby(zones).count() result_gdf = GeoDataFrame({'count_accessible_amenities': result, 'geometry': zones['geometry']}).set_crs('EPSG:4326') - return result_gdf \ No newline at end of file + return result_gdf.count_accessible_amenities + +def count_accessible_amenities_areamean(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, aws_profilename=None) -> GeoSeries: + return _count_accessible_amenities(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, weighting='area', aws_profilename=aws_profilename) + +def count_accessible_amenities_populationmean(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, aws_profilename=None) -> GeoSeries: + return _count_accessible_amenities(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, weighting='population', aws_profilename=aws_profilename) \ No newline at end of file diff --git a/tests/test_metrics.py b/tests/test_metrics.py index 4cc9d72..87f6bfb 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -25,7 +25,7 @@ def test_built_land_without_tree_cover(): assert expected_zone_size == actual_indicator_size @pytest.mark.skipif(EXECUTE_IGNORED_TESTS == False, reason="AWS credentials needed") -def test_count_accessible_amenities(): +def test_count_accessible_amenities_areamean(): from .conftest import create_fishnet_grid NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296, 0.01).reset_index() indicator = count_accessible_amenities(NAIROBI_BBOX, 'KEN-Nairobi', 'schools', 'walk', 'time', '15', 2024, aws_profilename=None) From ffc15eda7bd79e58e156fdb3861f923bd4c83171 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Mon, 13 Jan 2025 08:42:07 -0500 Subject: [PATCH 17/25] correct pop-weighted indicator calculation --- city_metrix/metrics/count_accessible_amenities.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index 0ce09c7..48ef588 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -62,12 +62,13 @@ def _count_accessible_amenities(zones: GeoDataFrame, cityname, amenityname, trav if weighting == 'area': result = population_layer.mask(count_layers[1]).groupby(zones).count() for count in range(2, max_count+1): - result += population_layer.mask(count_layers[count]).groupby(zones).count() * count + result += population_layer.mask(count_layers[count]).groupby(zones).count() + result = result / population_layer.groupby(zones).count() else: # weighting == 'population' result = population_layer.mask(count_layers[1]).groupby(zones).sum() for count in range(2, max_count+1): result += population_layer.mask(count_layers[count]).groupby(zones).sum() - result = result / population_layer.groupby(zones).count() + result = result / population_layer.groupby(zones).sum() result_gdf = GeoDataFrame({'count_accessible_amenities': result, 'geometry': zones['geometry']}).set_crs('EPSG:4326') return result_gdf.count_accessible_amenities From 6a7f86ca2a26ff74f018a4395c7590ef78f245d5 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Wed, 15 Jan 2025 17:37:36 -0500 Subject: [PATCH 18/25] clean up percent_population_access --- city_metrix/layers/isoline.py | 32 +---- city_metrix/layers/urban_land_use.py | 5 +- city_metrix/metrics/__init__.py | 2 +- .../metrics/percent_population_access.py | 118 ++++++------------ 4 files changed, 48 insertions(+), 109 deletions(-) diff --git a/city_metrix/layers/isoline.py b/city_metrix/layers/isoline.py index 523a489..d5b82a9 100644 --- a/city_metrix/layers/isoline.py +++ b/city_metrix/layers/isoline.py @@ -2,40 +2,20 @@ import shapely import boto3 -from .layer import Layer, get_utm_zone_epsg +from .layer import Layer +from city_metrix.layers.layer import get_utm_zone_epsg BUCKET_NAME = 'wri-cities-indicators' PREFIX = 'data/isolines/' class Isoline(Layer): - def __init__(self, params, aws_profilename=None, **kwargs): + def __init__(self, filename, **kwargs): super().__init__(**kwargs) - # params is dict with keys cityname, amenityname, travelmode, threshold_type, threshold_value, year - cityname = params['cityname'] - amenityname = params['amenityname'] - travelmode = params['travelmode'] - threshold_type = params['threshold_type'] - threshold_value = params['threshold_value'] - year = params['year'] - # Get list of isoline files from S3 and find the most recent file with correct cityname, amenity, etc - if aws_profilename is None: - session = boto3.Session() - else: - session = boto3.Session(profile_name=aws_profilename) - s3 = session.client('s3') - obj_list = s3.list_objects(Bucket=BUCKET_NAME, Prefix=PREFIX)['Contents'] - objnames = [i['Key'].replace(PREFIX, '') for i in obj_list if len(i['Key'].replace(PREFIX, '')) > 0] - matches = [oname for oname in objnames if oname.split('.')[0].split('_')[:-1] == [cityname, amenityname, travelmode, threshold_type, str(threshold_value)]] - if year is not None: - matches = [oname for oname in matches if oname.split('.')[0].split('_')[-1][:4] == str(year)] - if len(matches) == 0: - raise Exception('No isoline file found.') - matches.sort(key=lambda x: int(x.split('.')[0].split('_')[-1])) - objname = matches[-1] - # Get data from S3 - url = f'https://{BUCKET_NAME}.s3.us-east-1.amazonaws.com/{PREFIX}{objname}' + url = f'https://{BUCKET_NAME}.s3.us-east-1.amazonaws.com/{PREFIX}{filename}' + print(f'Attempting to retrieve isoline file from {url}', end=' ') gdf = gpd.read_file(url) + print('(Succeeded)') self.gdf = gpd.GeoDataFrame({'is_accessible': [1] * len(gdf), 'geometry': gdf.to_crs('EPSG:4326')['geometry']}).to_crs('EPSG:4326').set_crs('EPSG:4326').set_geometry('geometry') def get_data(self, bbox): diff --git a/city_metrix/layers/urban_land_use.py b/city_metrix/layers/urban_land_use.py index 1c81d39..fea653f 100644 --- a/city_metrix/layers/urban_land_use.py +++ b/city_metrix/layers/urban_land_use.py @@ -13,10 +13,11 @@ class UrbanLandUse(Layer): spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster) """ - def __init__(self, band='lulc', spatial_resolution=5, **kwargs): + def __init__(self, band='lulc', return_value=None, spatial_resolution=5, **kwargs): super().__init__(**kwargs) self.band = band self.spatial_resolution = spatial_resolution + self.return_value = return_value def get_data(self, bbox): ulu = ee.ImageCollection("projects/wri-datalab/cities/urban_land_use/V1") @@ -42,5 +43,7 @@ def get_data(self, bbox): self.spatial_resolution, "urban land use" ).lulc + if self.return_value is not None: + data = data.where(data==return_value) return data diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index 76f40ca..da9473d 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -5,7 +5,7 @@ from .urban_open_space import urban_open_space from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing -from .percent_population_access import percent_population_access, percent_population_access_elderly, percent_population_access_children, percent_population_access_female, percent_population_access_informal +from .percent_population_access import percent_population_access#, percent_population_access_elderly, percent_population_access_children, percent_population_access_female, percent_population_access_informal from .count_accessible_amenities import count_accessible_amenities_areamean, count_accessible_amenities_populationmean from .recreational_space_per_capita import recreational_space_per_capita from .vegetation_water_change import vegetation_water_change_gain_area diff --git a/city_metrix/metrics/percent_population_access.py b/city_metrix/metrics/percent_population_access.py index 3b0a713..7624582 100644 --- a/city_metrix/metrics/percent_population_access.py +++ b/city_metrix/metrics/percent_population_access.py @@ -1,89 +1,45 @@ from geopandas import GeoDataFrame, GeoSeries +import shapely import numpy as np from city_metrix.layers.isoline import Isoline from city_metrix.layers.world_pop import WorldPop from city_metrix.layers.urban_land_use import UrbanLandUse - - -def percent_population_access(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, agesex_classes=[], informal_only=False, worldpop_year=2020, aws_profilename=None) -> GeoSeries: - - class InformalTemp(UrbanLandUse): +from city_metrix.layers.layer import Layer, get_utm_zone_epsg + + +def percent_population_access(zones, city_name, amenity_name, travel_mode, threshold_type, threshold_value, retrieval_date, worldpop_agesex_classes=[], worldpop_year=2020, informal_only=False): + # Example: city_name='MEX-Mexico-City', amenity_name='schools', travel_mode='walk', threshold_type='time', threshold_value=30, retrieval_date='20241105' + + class IsolineSimplified(Isoline): + # Stores and returns isochrone or isodistance polygons with simplified geometry, and dissolved into fewer non-overlapping polygons + def __init__(self, filename, **kwargs): + super().__init__(filename=filename) + # params is dict with keys cityname, amenityname, travelmode, threshold_type, threshold_value, year + iso_gdf = self.gdf + if iso_gdf.crs in ('EPSG:4326', 'epsg:4326'): + utm_epsg = get_utm_zone_epsg(iso_gdf.total_bounds) + iso_gdf = iso_gdf.to_crs(utm_epsg).set_crs(utm_epsg) + poly_list = [shapely.simplify(p.buffer(0.1), tolerance=10) for p in iso_gdf.geometry] # Buffer and simplify + uu = shapely.unary_union(shapely.MultiPolygon(poly_list)) # Dissolve + shorter_gdf = GeoDataFrame({'accessible': 1, 'geometry': list(uu.geoms)}).set_crs(utm_epsg) + self.gdf = shorter_gdf.to_crs('EPSG:4326').set_crs('EPSG:4326') + def get_data(self, bbox): - data = super().get_data(bbox) - result = xr.where(data==3, 1, np.nan).rio.write_crs(data.crs) - result = result.assign_attrs(**data.attrs) - return result - - def get_accessible_population(access_features_layer, popraster_layer, zones): - if len(access_features_layer.gdf): - result_series = popraster_layer.mask(access_features_layer).groupby(zones).mean() * popraster_layer.mask(access_features_layer).groupby(zones).count() - else: - result_series = pd.Series([0] * len(zones)) - return result_series - - # cityname example: ARG-Buenos-Aires - # amenityname is OSMclass names, in lowercase - # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) - # threshold_type is distance or time - # threshold_value is integer, in minutes or meters + return self.gdf.clip(shapely.box(*bbox)) + filename = f"{city_name}_{amenity_name}_{travel_mode}_{threshold_type}_{threshold_value}_{retrieval_date}.geojson" + iso_layer = IsolineSimplified(filename) + accesspop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year, masks=[iso_layer,]) + totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) if informal_only: - informal_mask = InformalTemp() - population_layer = WorldPop(agesex_classes=agesex_classes, year=worldpop_year, masks=[informal_layer,]) - else: - population_layer = WorldPop(agesex_classes=agesex_classes, year=worldpop_year) - params = { - 'cityname': cityname, - 'amenityname': amenityname, - 'travelmode': travelmode, - 'threshold_type': threshold_type, - 'threshold_value': threshold_value, - 'year': isoline_year - } - accesszone_layer = Isoline(params, aws_profilename=aws_profilename) - - result_gdf = GeoDataFrame({'geometry': zones['geometry']}).set_geometry('geometry').set_crs('EPSG:4326') - - try: - access_pop = get_accessible_population(accesszone_layer, population_layer, zones) - total_pop = population_layer.groupby(zones).mean() * population_layer.groupby(zones).count() - result = (access_pop / total_pop) * 100 - result_gdf['access_fraction'] = result - - except: - # Sometimes doing entire zones gdf causes groupby to throw empty-GDF error -- workaraound is to go district-by-district - print('Calculating district-by-district') - - access_fraction = [] - for idx in zones.index: - try: # Sometimes there's still an empty-gdf error - access_pop = get_accessible_population(accesszone_layer, population_layer, zones.loc[[zones.index[idx]]])[0] - total_pop = (population_layer.groupby(zones.loc[[zones.index[idx]]]).mean() * population_layer.groupby(zones.loc[[zones.index[idx]]]).count())[0] - if total_pop != 0: - access_fraction.append(access_pop / total_pop) - else: - access_fraction.append(np.nan) - except: - print('Empty-GDF error for index {0}'.format(idx)) - access_fraction.append(np.nan) - result_gdf['access_fraction'] = access_fraction - result_gdf['access_fraction'].replace([np.inf, -np.inf], np.nan, inplace=True) - result_gdf['access_fraction'] = result_gdf['access_fraction'] * 100 - - return result_gdf - -def percent_population_access_elderly(zones, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, worldpop_year=2020, aws_profilename=None): - agesex_classes=['F_60', 'F_65', 'F_70', 'F_75', 'F_80', 'M_60', 'M_65', 'M_70', 'M_75', 'M_80'] - return percent_population_access(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, agesex_classes=agesex_classes, informal_only=False, worldpop_year=worldpop_year, aws_profilename=aws_profilename) - -def percent_population_access_children(zones, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, worldpop_year=2020, aws_profilename=None): - agesex_classes=['F_0', 'F_1', 'F_5', 'F_10', 'M_0', 'M_1', 'M_5', 'M_10'] - return percent_population_access(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, agesex_classes=agesex_classes, informal_only=False, worldpop_year=worldpop_year, aws_profilename=aws_profilename) - -def percent_population_access_female(zones, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, worldpop_year=2020, aws_profilename=None): - agesex_classes=['F_0', 'F_1', 'F_5', 'F_10', 'F_15', 'F_20', 'F_25', 'F_30', 'F_35', 'F_40', 'F_45', 'F_50', 'F_55', 'F_60', 'F_65', 'F_70', 'F_75', 'F_80'] - return percent_population_access(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, agesex_classes=agesex_classes, informal_only=False, worldpop_year=worldpop_year, aws_profilename=aws_profilename) - -def percent_population_access_informal(zones, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, worldpop_year=2020, aws_profilename=None): - agesex_classes=[] - return percent_population_access(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, agesex_classes=agesex_classes, informal_only=True, worldpop_year=worldpop_year, aws_profilename=aws_profilename) \ No newline at end of file + informal_layer = UrbanLandUse(return_value=3) + accesspop_layer.masks.append(informal_layer) + totalpop_layer.masks.append(informal_layer) + res = [] + for rownum in range(len(zones)): # Doing it district-by-district to avoid empty-GDF error + try: + res.append(100 * accesspop_layer.groupby(zones.iloc[[rownum]]).sum()[0] / totalpop_layer.groupby(zones.iloc[[rownum]]).sum()[0]) + except: + res.append(0) + result_gdf = GeoDataFrame({'access_percent': res, 'geometry': zones.geometry}).set_crs(zones.crs) + return result_gdf \ No newline at end of file From c5889be0866cf3e156511bbae017357e841499e6 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Wed, 15 Jan 2025 18:03:12 -0500 Subject: [PATCH 19/25] bug fix --- city_metrix/layers/urban_land_use.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/city_metrix/layers/urban_land_use.py b/city_metrix/layers/urban_land_use.py index fea653f..0716a92 100644 --- a/city_metrix/layers/urban_land_use.py +++ b/city_metrix/layers/urban_land_use.py @@ -43,7 +43,8 @@ def get_data(self, bbox): self.spatial_resolution, "urban land use" ).lulc + if self.return_value is not None: - data = data.where(data==return_value) + data = data.where(data==self.return_value) return data From e625c390ef54cb14d4f728742a773964d578c2d2 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Fri, 17 Jan 2025 15:34:13 -0500 Subject: [PATCH 20/25] overhaul count_accessible_amenities metric --- city_metrix/metrics/__init__.py | 2 +- .../metrics/count_accessible_amenities.py | 184 +++++++++++------- tests/test_metrics.py | 9 +- 3 files changed, 123 insertions(+), 72 deletions(-) diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index da9473d..c40d230 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -6,7 +6,7 @@ from .natural_areas import natural_areas from .era_5_met_preprocessing import era_5_met_preprocessing from .percent_population_access import percent_population_access#, percent_population_access_elderly, percent_population_access_children, percent_population_access_female, percent_population_access_informal -from .count_accessible_amenities import count_accessible_amenities_areamean, count_accessible_amenities_populationmean +from .count_accessible_amenities import count_accessible_amenities#_areamean, count_accessible_amenities_populationmean from .recreational_space_per_capita import recreational_space_per_capita from .vegetation_water_change import vegetation_water_change_gain_area from .vegetation_water_change import vegetation_water_change_loss_area diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index 48ef588..5cd60d7 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -1,79 +1,131 @@ from geopandas import GeoDataFrame, GeoSeries +from pandas import Series from city_metrix.layers import Layer +from math import floor import shapely -from city_metrix.layers.isoline import Isoline -from city_metrix.layers.world_pop import WorldPop +from city_metrix.layers import Isoline, UrbanLandUse, WorldPop +from city_metrix.layers.layer import get_utm_zone_epsg -class AccessCountTmp(Layer): - def __init__(self, access_gdf, **kwargs): - super().__init__(**kwargs) - self.gdf = access_gdf - def get_data(self, bbox): - return self.gdf.clip(shapely.box(*bbox)) - -def _count_accessible_amenities(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, weighting='population', aws_profilename=None) -> GeoSeries: + +def count_accessible_amenities(zones: GeoDataFrame, city_name, amenity_name, travel_mode, threshold_type, threshold_value, retrieval_date, weighting='population', worldpop_agesex_classes=[], worldpop_year=2020, informal_only=False) -> GeoSeries: # cityname example: ARG-Buenos-Aires # amenityname is OSMclass names, in lowercase # travelmode is walk, bike, automobile, publictransit (only walk implemented for now) # threshold_type is distance or time # threshold_value is integer, in minutes or meters - population_layer = WorldPop(agesex_classes=[], year=2020) - params = { - 'cityname': cityname, - 'amenityname': amenityname, - 'travelmode': travelmode, - 'threshold_type': threshold_type, - 'threshold_value': threshold_value, - 'year': isoline_year - } - iso_layer = Isoline(params, aws_profilename=aws_profilename) - iso_gdf = iso_layer.get_data(zones.total_bounds) - - # Collect individual boundary linestrings of each isoline polygon - iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] - iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] - iso_boundarylines = [] - for i in iso_boundarylinesmulti: - if type(i) == shapely.geometry.linestring.LineString: - iso_boundarylines.append(i) - else: - for j in list(i.geoms): - iso_boundarylines.append(j) - iso_bl_gdf = GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) - - # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" - iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) - accesscount_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') - - # For each dissected polygon, find how many of the original isoline polys contain the centroid - # This is the number of amenity points are accessible within the dissected polygon - count_amenities = iso_dissected_polys.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 - for iso_idx in range(1, len(iso_gdf)): - count_amenities = count_amenities + (iso_dissected_polys.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) - - accesscount_gdf['count_amenities'] = count_amenities - - count_layers = {count: AccessCountTmp(GeoDataFrame({'count_amenities': accesscount_gdf['count_amenities']==count, 'geometry': accesscount_gdf['geometry']}).set_crs('EPSG:4326')) for count in range(1, accesscount_gdf['count_amenities'].max() + 1)} - - # For each zone, find average number of accessible amenities, and store in result_gdf - max_count = accesscount_gdf['count_amenities'].max() - if weighting == 'area': - result = population_layer.mask(count_layers[1]).groupby(zones).count() - for count in range(2, max_count+1): - result += population_layer.mask(count_layers[count]).groupby(zones).count() - result = result / population_layer.groupby(zones).count() - else: # weighting == 'population' - result = population_layer.mask(count_layers[1]).groupby(zones).sum() - for count in range(2, max_count+1): - result += population_layer.mask(count_layers[count]).groupby(zones).sum() - result = result / population_layer.groupby(zones).sum() - result_gdf = GeoDataFrame({'count_accessible_amenities': result, 'geometry': zones['geometry']}).set_crs('EPSG:4326') - return result_gdf.count_accessible_amenities + class AccessCountTmp(Layer): + def __init__(self, access_gdf, return_value, **kwargs): + super().__init__(**kwargs) + self.gdf = access_gdf[access_gdf['count_amenities']==return_value] + def get_data(self, bbox): + return self.gdf.clip(shapely.box(*bbox)) + + class IsolinesBuffered(Isoline): + # Stores and returns isochrone or isodistance polygons with polygons slightly buffered and simplified to avoid invalid-geometry errors + def __init__(self, filename, **kwargs): + super().__init__(filename=filename) + # params is dict with keys cityname, amenityname, travelmode, threshold_type, threshold_value, year + iso_gdf = self.gdf + if iso_gdf.crs in ('EPSG:4326', 'epsg:4326'): + utm_epsg = get_utm_zone_epsg(iso_gdf.total_bounds) + iso_gdf = iso_gdf.to_crs(utm_epsg).set_crs(utm_epsg) + poly_list = [shapely.simplify(p.buffer(0.1), tolerance=10) for p in iso_gdf.geometry] # Buffer and simplify + buffered_gdf = GeoDataFrame({'accessible': 1, 'geometry': poly_list}).set_crs(utm_epsg) + self.gdf = buffered_gdf.to_crs('EPSG:4326').set_crs('EPSG:4326') + + def get_data(self, bbox): + return self.gdf.clip(shapely.box(*bbox)) + + filename = f"{city_name}_{amenity_name}_{travel_mode}_{threshold_type}_{threshold_value}_{retrieval_date}.geojson" + + population_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + if informal_only: + informal_layer = UrbanLandUse(return_value=3) + population_layer.masks.append(informal_layer) + + filename = f"{city_name}_{amenity_name}_{travel_mode}_{threshold_type}_{threshold_value}_{retrieval_date}.geojson" + iso_layer = IsolinesBuffered(filename=filename) + + results = [] + for rownum in range(len(zones)): + print("\n{0} (geo {1} of {2})".format(zones['geo_name'][rownum], rownum+1, len(zones))) + zone = zones.iloc[[rownum]] + iso_gdf = iso_layer.get_data(zone.total_bounds) + # Collect individual boundary linestrings of each isoline polygon + iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] + iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] + iso_boundarylines = [] + for i in iso_boundarylinesmulti: + if type(i) == shapely.geometry.linestring.LineString: + iso_boundarylines.append(i) + else: + for j in list(i.geoms): + iso_boundarylines.append(j) + iso_bl_gdf = GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) + print("Converted isoline polygons into boundary multilines") + # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" + iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) + print("Creating gdf of dissected polygons") + dissected_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') + # For each dissected polygon, find how many of the original isoline polys contain the centroid + # This is the number of amenity points are accessible within the dissected polygon + print("Counting how many original isoline polygons intersect with each dissected polygon") + count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, with num rows == num of dissected polys + print("|==" + "PROGRESS=OF=INCLUSION=TESTS" + ("=" * 71) + "|") + print(" ", end="") + progress = 0 + for iso_idx in range(1, len(iso_gdf)): # Iterate through all original isoline polys: test dissected polys to see whether centroids are included in isoline poly, add to running sum Series + if floor(100 * iso_idx / len(iso_gdf)) > progress: + progress = floor(100 * iso_idx / len(iso_gdf)) + print("X", end="") + count_amenities = count_amenities + (dissected_gdf.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) + print("X") + dissected_gdf['count_amenities'] = count_amenities + + # Create dict of polygons each with a single asset-count + max_count = dissected_gdf['count_amenities'].max() + count_layers = {count: AccessCountTmp(access_gdf=dissected_gdf, return_value=count) for count in range(1, max_count + 1)} -def count_accessible_amenities_areamean(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, aws_profilename=None) -> GeoSeries: - return _count_accessible_amenities(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, weighting='area', aws_profilename=aws_profilename) + # For each zone, find average number of accessible amenities, and store in result_gdf + rowresults = [] + for rownum in range(len(zone)): + running_sum = Series([0] * len(zone)) + for count in range(1, max_count+1): + try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() + pop_layer + except NameError: + pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + pop_layer.masks = [] + try: + totalpop_layer + except NameError: + totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + totalpop_layer.masks = [] + if informal_only: + pop_layer.masks.append(informal_layer) + totalpop_layer.masks.append(informal_layer) + pop_layer.masks.append(count_layers[count]) + groupby = pop_layer.groupby(zone.iloc[[rownum]]) + if weighting == 'area': + try: + running_sum += count * groupby.count().fillna(0) + except: + running_sum += 0 + else: # weighting == 'population' + try: + running_sum += count * groupby.sum().fillna(0) + except: + running_sum += 0 + if weighting == 'area': + rowresults.append(running_sum / totalpop_layer.groupby(zone).count()) + else: # weighting == 'population' + rowresults.append(running_sum / totalpop_layer.groupby(zone).sum()) -def count_accessible_amenities_populationmean(zones: GeoDataFrame, cityname, amenityname, travelmode, threshold_type, threshold_value, isoline_year=2024, aws_profilename=None) -> GeoSeries: - return _count_accessible_amenities(zones=zones, cityname=cityname, amenityname=amenityname, travelmode=travelmode, threshold_type=threshold_type, threshold_value=threshold_value, isoline_year=isoline_year, weighting='population', aws_profilename=aws_profilename) \ No newline at end of file + rowresult_gdf = GeoDataFrame({f'{weighting}_averaged_num_accessible_amenities': rowresults, 'geometry': zone['geometry']}).set_crs('EPSG:4326') + results.append(rowresult_gdf[f'{weighting}_averaged_num_accessible_amenities'][rownum]) + result_gdf = GeoDataFrame({f'{weighting}_averaged_num_accessible_amenities': results, 'geometry': zones['geometry']}).set_crs('EPSG:4326') + return result_gdf \ No newline at end of file diff --git a/tests/test_metrics.py b/tests/test_metrics.py index eeaca06..d4f32f5 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -24,11 +24,10 @@ def test_built_land_without_tree_cover(): actual_indicator_size = indicator.size assert expected_zone_size == actual_indicator_size -@pytest.mark.skipif(EXECUTE_IGNORED_TESTS == False, reason="AWS credentials needed") -def test_count_accessible_amenities_areamean(): - from .conftest import create_fishnet_grid - NAIROBI_BBOX = create_fishnet_grid(36.66446402, -1.44560888, 37.10497899, -1.16058296, 0.01).reset_index() - indicator = count_accessible_amenities(NAIROBI_BBOX, 'KEN-Nairobi', 'schools', 'walk', 'time', '15', 2024, aws_profilename=None) +@pytest.mark.skipif(EXECUTE_IGNORED_TESTS == False, reason="Needs specific zone file to run") +def test_count_accessible_amenities(): + nairobi_gdf = None # Need link to a GDF for which the isoline file has been calculated + indicator = count_accessible_amenities(nairobi_gdf, 'KEN-Nairobi', 'schools', 'walk', 'time', 15, '20241105', weighting='population') expected_zone_size = NAIROBI_BBOX.geometry.size actual_indicator_size = indicator.size assert expected_zone_size == actual_indicator_size From 9be9c619a3c517ae57784a3653aee6bd5c55c96f Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Fri, 17 Jan 2025 17:12:20 -0500 Subject: [PATCH 21/25] bug fixes --- .../metrics/count_accessible_amenities.py | 69 +++++++++---------- 1 file changed, 34 insertions(+), 35 deletions(-) diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index 5cd60d7..7300a97 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -52,6 +52,7 @@ def get_data(self, bbox): for rownum in range(len(zones)): print("\n{0} (geo {1} of {2})".format(zones['geo_name'][rownum], rownum+1, len(zones))) zone = zones.iloc[[rownum]] + iso_gdf = iso_layer.get_data(zone.total_bounds) # Collect individual boundary linestrings of each isoline polygon iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] @@ -89,43 +90,41 @@ def get_data(self, bbox): count_layers = {count: AccessCountTmp(access_gdf=dissected_gdf, return_value=count) for count in range(1, max_count + 1)} # For each zone, find average number of accessible amenities, and store in result_gdf - rowresults = [] - for rownum in range(len(zone)): - running_sum = Series([0] * len(zone)) - for count in range(1, max_count+1): - try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() - pop_layer - except NameError: - pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) - else: - pop_layer.masks = [] - try: - totalpop_layer - except NameError: - totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) - else: - totalpop_layer.masks = [] - if informal_only: - pop_layer.masks.append(informal_layer) - totalpop_layer.masks.append(informal_layer) - pop_layer.masks.append(count_layers[count]) - groupby = pop_layer.groupby(zone.iloc[[rownum]]) - if weighting == 'area': - try: - running_sum += count * groupby.count().fillna(0) - except: - running_sum += 0 - else: # weighting == 'population' - try: - running_sum += count * groupby.sum().fillna(0) - except: - running_sum += 0 + + running_sum = Series([0]) + for count in range(1, max_count+1): + try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() + pop_layer + except NameError: + pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + pop_layer.masks = [] + try: + totalpop_layer + except NameError: + totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + totalpop_layer.masks = [] + if informal_only: + pop_layer.masks.append(informal_layer) + totalpop_layer.masks.append(informal_layer) + pop_layer.masks.append(count_layers[count]) + groupby = pop_layer.groupby(zone) if weighting == 'area': - rowresults.append(running_sum / totalpop_layer.groupby(zone).count()) + try: + running_sum += count * groupby.count().fillna(0) + except: + running_sum += Series([0]) else: # weighting == 'population' - rowresults.append(running_sum / totalpop_layer.groupby(zone).sum()) + try: + running_sum += count * groupby.sum().fillna(0) + except: + running_sum += Series([0]) + if weighting == 'area': + rowresult = running_sum / totalpop_layer.groupby(zone).count() + else: # weighting == 'population' + rowresult = running_sum / totalpop_layer.groupby(zone).sum() - rowresult_gdf = GeoDataFrame({f'{weighting}_averaged_num_accessible_amenities': rowresults, 'geometry': zone['geometry']}).set_crs('EPSG:4326') - results.append(rowresult_gdf[f'{weighting}_averaged_num_accessible_amenities'][rownum]) + results.append(rowresult[0]) result_gdf = GeoDataFrame({f'{weighting}_averaged_num_accessible_amenities': results, 'geometry': zones['geometry']}).set_crs('EPSG:4326') return result_gdf \ No newline at end of file From a65984e5672aa8061c8223ed63491aaccc88251f Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Fri, 17 Jan 2025 19:58:59 -0500 Subject: [PATCH 22/25] bug fix for district with zero isoline polys --- .../metrics/count_accessible_amenities.py | 135 +++++++++--------- 1 file changed, 69 insertions(+), 66 deletions(-) diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index 7300a97..7f44be1 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -54,77 +54,80 @@ def get_data(self, bbox): zone = zones.iloc[[rownum]] iso_gdf = iso_layer.get_data(zone.total_bounds) - # Collect individual boundary linestrings of each isoline polygon - iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] - iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] - iso_boundarylines = [] - for i in iso_boundarylinesmulti: - if type(i) == shapely.geometry.linestring.LineString: - iso_boundarylines.append(i) - else: - for j in list(i.geoms): - iso_boundarylines.append(j) - iso_bl_gdf = GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) - print("Converted isoline polygons into boundary multilines") - # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" - iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) - print("Creating gdf of dissected polygons") - dissected_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') - # For each dissected polygon, find how many of the original isoline polys contain the centroid - # This is the number of amenity points are accessible within the dissected polygon - print("Counting how many original isoline polygons intersect with each dissected polygon") - count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, with num rows == num of dissected polys - print("|==" + "PROGRESS=OF=INCLUSION=TESTS" + ("=" * 71) + "|") - print(" ", end="") - progress = 0 - for iso_idx in range(1, len(iso_gdf)): # Iterate through all original isoline polys: test dissected polys to see whether centroids are included in isoline poly, add to running sum Series - if floor(100 * iso_idx / len(iso_gdf)) > progress: - progress = floor(100 * iso_idx / len(iso_gdf)) - print("X", end="") - count_amenities = count_amenities + (dissected_gdf.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) - print("X") - dissected_gdf['count_amenities'] = count_amenities + if len(iso_gdf) == 0: + results.append(0) + else: + # Collect individual boundary linestrings of each isoline polygon + iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] + iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] + iso_boundarylines = [] + for i in iso_boundarylinesmulti: + if type(i) == shapely.geometry.linestring.LineString: + iso_boundarylines.append(i) + else: + for j in list(i.geoms): + iso_boundarylines.append(j) + iso_bl_gdf = GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) + print("Converted isoline polygons into boundary multilines") + # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" + iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) + print("Creating gdf of dissected polygons") + dissected_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') + # For each dissected polygon, find how many of the original isoline polys contain the centroid + # This is the number of amenity points are accessible within the dissected polygon + print("Counting how many original isoline polygons intersect with each dissected polygon") + count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, with num rows == num of dissected polys + print("|==" + "PROGRESS=OF=INCLUSION=TESTS" + ("=" * 71) + "|") + print(" ", end="") + progress = 0 + for iso_idx in range(1, len(iso_gdf)): # Iterate through all original isoline polys: test dissected polys to see whether centroids are included in isoline poly, add to running sum Series + if floor(100 * iso_idx / len(iso_gdf)) > progress: + progress = floor(100 * iso_idx / len(iso_gdf)) + print("X", end="") + count_amenities = count_amenities + (dissected_gdf.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) + print("X") + dissected_gdf['count_amenities'] = count_amenities - # Create dict of polygons each with a single asset-count - max_count = dissected_gdf['count_amenities'].max() - count_layers = {count: AccessCountTmp(access_gdf=dissected_gdf, return_value=count) for count in range(1, max_count + 1)} + # Create dict of polygons each with a single asset-count + max_count = dissected_gdf['count_amenities'].max() + count_layers = {count: AccessCountTmp(access_gdf=dissected_gdf, return_value=count) for count in range(1, max_count + 1)} - # For each zone, find average number of accessible amenities, and store in result_gdf + # For each zone, find average number of accessible amenities, and store in result_gdf - running_sum = Series([0]) - for count in range(1, max_count+1): - try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() - pop_layer - except NameError: - pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) - else: - pop_layer.masks = [] - try: - totalpop_layer - except NameError: - totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) - else: - totalpop_layer.masks = [] - if informal_only: - pop_layer.masks.append(informal_layer) - totalpop_layer.masks.append(informal_layer) - pop_layer.masks.append(count_layers[count]) - groupby = pop_layer.groupby(zone) - if weighting == 'area': + running_sum = Series([0]) + for count in range(1, max_count+1): + try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() + pop_layer + except NameError: + pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + pop_layer.masks = [] try: - running_sum += count * groupby.count().fillna(0) - except: - running_sum += Series([0]) + totalpop_layer + except NameError: + totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + totalpop_layer.masks = [] + if informal_only: + pop_layer.masks.append(informal_layer) + totalpop_layer.masks.append(informal_layer) + pop_layer.masks.append(count_layers[count]) + groupby = pop_layer.groupby(zone) + if weighting == 'area': + try: + running_sum += count * groupby.count().fillna(0) + except: + running_sum += Series([0]) + else: # weighting == 'population' + try: + running_sum += count * groupby.sum().fillna(0) + except: + running_sum += Series([0]) + if weighting == 'area': + rowresult = running_sum / totalpop_layer.groupby(zone).count() else: # weighting == 'population' - try: - running_sum += count * groupby.sum().fillna(0) - except: - running_sum += Series([0]) - if weighting == 'area': - rowresult = running_sum / totalpop_layer.groupby(zone).count() - else: # weighting == 'population' - rowresult = running_sum / totalpop_layer.groupby(zone).sum() + rowresult = running_sum / totalpop_layer.groupby(zone).sum() - results.append(rowresult[0]) + results.append(rowresult[0]) result_gdf = GeoDataFrame({f'{weighting}_averaged_num_accessible_amenities': results, 'geometry': zones['geometry']}).set_crs('EPSG:4326') return result_gdf \ No newline at end of file From 63536a440f26678425f8ef0ae4d2efeaa75a37b8 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Sat, 18 Jan 2025 09:41:14 -0500 Subject: [PATCH 23/25] Bug fix for single simple isoline polygon --- .../metrics/count_accessible_amenities.py | 119 +++++++++--------- 1 file changed, 59 insertions(+), 60 deletions(-) diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index 7f44be1..c438c06 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -46,16 +46,12 @@ def get_data(self, bbox): population_layer.masks.append(informal_layer) filename = f"{city_name}_{amenity_name}_{travel_mode}_{threshold_type}_{threshold_value}_{retrieval_date}.geojson" - iso_layer = IsolinesBuffered(filename=filename) - - results = [] - for rownum in range(len(zones)): - print("\n{0} (geo {1} of {2})".format(zones['geo_name'][rownum], rownum+1, len(zones))) - zone = zones.iloc[[rownum]] - - iso_gdf = iso_layer.get_data(zone.total_bounds) - if len(iso_gdf) == 0: - results.append(0) + iso_gdf = iso_layer.get_data(zone.total_bounds) + if len(iso_gdf) == 0: + result.append(0) + else: + if len(iso_gdf) == 1 and type(iso_gdf.iloc[0].geometry) == shapely.geometry.polygon.Polygon: # only one simple-polygon isoline + dissected_gdf = GeoDataFrame({'poly_id': [0], 'geometry': [iso_gdf.iloc[0].geometry]}).set_crs('EPSG:4326') else: # Collect individual boundary linestrings of each isoline polygon iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] @@ -70,64 +66,67 @@ def get_data(self, bbox): iso_bl_gdf = GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) print("Converted isoline polygons into boundary multilines") # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" + iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) print("Creating gdf of dissected polygons") dissected_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') - # For each dissected polygon, find how many of the original isoline polys contain the centroid - # This is the number of amenity points are accessible within the dissected polygon - print("Counting how many original isoline polygons intersect with each dissected polygon") - count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, with num rows == num of dissected polys - print("|==" + "PROGRESS=OF=INCLUSION=TESTS" + ("=" * 71) + "|") - print(" ", end="") - progress = 0 - for iso_idx in range(1, len(iso_gdf)): # Iterate through all original isoline polys: test dissected polys to see whether centroids are included in isoline poly, add to running sum Series - if floor(100 * iso_idx / len(iso_gdf)) > progress: - progress = floor(100 * iso_idx / len(iso_gdf)) - print("X", end="") - count_amenities = count_amenities + (dissected_gdf.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) - print("X") - dissected_gdf['count_amenities'] = count_amenities + # For each dissected polygon, find how many of the original isoline polys contain the centroid + # This is the number of amenity points are accessible within the dissected polygon + print("Counting how many original isoline polygons intersect with each dissected polygon") + count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, with num rows == num of dissected polys + print("|==" + "PROGRESS=OF=INCLUSION=TESTS" + ("=" * 71) + "|") + print(" ", end="") + progress = 0 - # Create dict of polygons each with a single asset-count - max_count = dissected_gdf['count_amenities'].max() - count_layers = {count: AccessCountTmp(access_gdf=dissected_gdf, return_value=count) for count in range(1, max_count + 1)} + for iso_idx in range(1, len(iso_gdf)): # Iterate through all original isoline polys: test dissected polys to see whether centroids are included in isoline poly, add to running sum Series + if floor(100 * iso_idx / len(iso_gdf)) > progress: + progress = floor(100 * iso_idx / len(iso_gdf)) + print("X", end="") + count_amenities = count_amenities + (dissected_gdf.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) + print("X") + dissected_gdf['count_amenities'] = count_amenities - # For each zone, find average number of accessible amenities, and store in result_gdf + # Create dict of polygons each with a single asset-count + max_count = dissected_gdf['count_amenities'].max() + count_layers = {count: AccessCountTmp(access_gdf=dissected_gdf, return_value=count) for count in range(1, max_count + 1)} - running_sum = Series([0]) - for count in range(1, max_count+1): - try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() - pop_layer - except NameError: - pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) - else: - pop_layer.masks = [] - try: - totalpop_layer - except NameError: - totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) - else: - totalpop_layer.masks = [] - if informal_only: - pop_layer.masks.append(informal_layer) - totalpop_layer.masks.append(informal_layer) - pop_layer.masks.append(count_layers[count]) - groupby = pop_layer.groupby(zone) - if weighting == 'area': - try: - running_sum += count * groupby.count().fillna(0) - except: - running_sum += Series([0]) - else: # weighting == 'population' - try: - running_sum += count * groupby.sum().fillna(0) - except: - running_sum += Series([0]) + # For each zone, find average number of accessible amenities, and store in result_gdf + + running_sum = Series([0]) + for count in range(1, max_count+1): + try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() + pop_layer + except NameError: + pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + pop_layer.masks = [] + try: + totalpop_layer + except NameError: + totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + totalpop_layer.masks = [] + if informal_only: + pop_layer.masks.append(informal_layer) + totalpop_layer.masks.append(informal_layer) + pop_layer.masks.append(count_layers[count]) + groupby = pop_layer.groupby(zone) if weighting == 'area': - rowresult = running_sum / totalpop_layer.groupby(zone).count() + try: + running_sum += count * groupby.count().fillna(0) + except: + running_sum += Series([0]) else: # weighting == 'population' - rowresult = running_sum / totalpop_layer.groupby(zone).sum() + try: + running_sum += count * groupby.sum().fillna(0) + except: + running_sum += Series([0]) + if weighting == 'area': + rowresult = running_sum / totalpop_layer.groupby(zone).count() + else: # weighting == 'population' + rowresult = running_sum / totalpop_layer.groupby(zone).sum() + + results.append(rowresult[0]) - results.append(rowresult[0]) result_gdf = GeoDataFrame({f'{weighting}_averaged_num_accessible_amenities': results, 'geometry': zones['geometry']}).set_crs('EPSG:4326') return result_gdf \ No newline at end of file From 161c71a8f7dffef9f104a3379c71c4e21ba46f84 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Sat, 18 Jan 2025 15:22:21 -0500 Subject: [PATCH 24/25] typo fix --- .../metrics/count_accessible_amenities.py | 149 +++++++++--------- 1 file changed, 77 insertions(+), 72 deletions(-) diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index c438c06..5bb5839 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -46,87 +46,92 @@ def get_data(self, bbox): population_layer.masks.append(informal_layer) filename = f"{city_name}_{amenity_name}_{travel_mode}_{threshold_type}_{threshold_value}_{retrieval_date}.geojson" - iso_gdf = iso_layer.get_data(zone.total_bounds) - if len(iso_gdf) == 0: - result.append(0) - else: - if len(iso_gdf) == 1 and type(iso_gdf.iloc[0].geometry) == shapely.geometry.polygon.Polygon: # only one simple-polygon isoline - dissected_gdf = GeoDataFrame({'poly_id': [0], 'geometry': [iso_gdf.iloc[0].geometry]}).set_crs('EPSG:4326') + iso_layer = IsolinesBuffered(filename=filename) + results = [] + for rownum in range(len(zones)): + zone = zones.iloc[[rownum]] + print('\n{0} (geo {1} of {2})'.format(zone.geo_name[rownum], rownum+1, len(zones))) + iso_gdf = iso_layer.get_data(zone.total_bounds) + if len(iso_gdf) == 0: + results.append(0) else: - # Collect individual boundary linestrings of each isoline polygon - iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] - iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] - iso_boundarylines = [] - for i in iso_boundarylinesmulti: - if type(i) == shapely.geometry.linestring.LineString: - iso_boundarylines.append(i) - else: - for j in list(i.geoms): - iso_boundarylines.append(j) - iso_bl_gdf = GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) - print("Converted isoline polygons into boundary multilines") - # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" + if len(iso_gdf) == 1 and type(iso_gdf.iloc[0].geometry) == shapely.geometry.polygon.Polygon: # only one simple-polygon isoline + dissected_gdf = GeoDataFrame({'poly_id': [0], 'geometry': [iso_gdf.iloc[0].geometry]}).set_crs('EPSG:4326') + else: + # Collect individual boundary linestrings of each isoline polygon + iso_eachboundary = [iso_gdf.iloc[rownum]['geometry'].boundary for rownum in range(len(iso_gdf))] + iso_boundarylinesmulti = [i for i in iso_eachboundary if i is not None] + iso_boundarylines = [] + for i in iso_boundarylinesmulti: + if type(i) == shapely.geometry.linestring.LineString: + iso_boundarylines.append(i) + else: + for j in list(i.geoms): + iso_boundarylines.append(j) + iso_bl_gdf = GeoDataFrame({'idx': range(len(iso_boundarylines)), 'geometry': iso_boundarylines}) + print("Converted isoline polygons into boundary multilines") + # Dissolve all linestrings into large multilinestring, and polygonize into "dissected polygons" - iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) - print("Creating gdf of dissected polygons") - dissected_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') - # For each dissected polygon, find how many of the original isoline polys contain the centroid - # This is the number of amenity points are accessible within the dissected polygon - print("Counting how many original isoline polygons intersect with each dissected polygon") - count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, with num rows == num of dissected polys - print("|==" + "PROGRESS=OF=INCLUSION=TESTS" + ("=" * 71) + "|") - print(" ", end="") - progress = 0 + iso_dissected_polys = shapely.polygonize(list(iso_bl_gdf.dissolve().geometry[0].geoms)) + print("Creating gdf of dissected polygons") + dissected_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') + # For each dissected polygon, find how many of the original isoline polys contain the centroid + # This is the number of amenity points are accessible within the dissected polygon + print("Counting how many original isoline polygons intersect with each dissected polygon") + count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, with num rows == num of dissected polys + print("|==" + "PROGRESS=OF=INCLUSION=TESTS" + ("=" * 71) + "|") + print(" ", end="") + progress = 0 - for iso_idx in range(1, len(iso_gdf)): # Iterate through all original isoline polys: test dissected polys to see whether centroids are included in isoline poly, add to running sum Series - if floor(100 * iso_idx / len(iso_gdf)) > progress: - progress = floor(100 * iso_idx / len(iso_gdf)) - print("X", end="") - count_amenities = count_amenities + (dissected_gdf.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) - print("X") - dissected_gdf['count_amenities'] = count_amenities + for iso_idx in range(1, len(iso_gdf)): # Iterate through all original isoline polys: test dissected polys to see whether centroids are included in isoline poly, add to running sum Series + if floor(100 * iso_idx / len(iso_gdf)) > progress: + progress = floor(100 * iso_idx / len(iso_gdf)) + print("X", end="") + count_amenities = count_amenities + (dissected_gdf.centroid.within(iso_gdf.iloc[[iso_idx]].geometry[iso_gdf.index[iso_idx]]) * 1) + print("X") + dissected_gdf['count_amenities'] = count_amenities - # Create dict of polygons each with a single asset-count - max_count = dissected_gdf['count_amenities'].max() - count_layers = {count: AccessCountTmp(access_gdf=dissected_gdf, return_value=count) for count in range(1, max_count + 1)} + # Create dict of polygons each with a single asset-count + max_count = dissected_gdf['count_amenities'].max() + count_layers = {count: AccessCountTmp(access_gdf=dissected_gdf, return_value=count) for count in range(1, max_count + 1)} - # For each zone, find average number of accessible amenities, and store in result_gdf + # For each zone, find average number of accessible amenities, and store in result_gdf - running_sum = Series([0]) - for count in range(1, max_count+1): - try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() - pop_layer - except NameError: - pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) - else: - pop_layer.masks = [] - try: - totalpop_layer - except NameError: - totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) - else: - totalpop_layer.masks = [] - if informal_only: - pop_layer.masks.append(informal_layer) - totalpop_layer.masks.append(informal_layer) - pop_layer.masks.append(count_layers[count]) - groupby = pop_layer.groupby(zone) - if weighting == 'area': + running_sum = Series([0]) + for count in range(1, max_count+1): + try: # Because adding masks to pop_layer adds them to WorldPop(), and they cannot be removed from WorldPop() + pop_layer + except NameError: + pop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + pop_layer.masks = [] try: - running_sum += count * groupby.count().fillna(0) - except: - running_sum += Series([0]) + totalpop_layer + except NameError: + totalpop_layer = WorldPop(agesex_classes=worldpop_agesex_classes, year=worldpop_year) + else: + totalpop_layer.masks = [] + if informal_only: + pop_layer.masks.append(informal_layer) + totalpop_layer.masks.append(informal_layer) + pop_layer.masks.append(count_layers[count]) + groupby = pop_layer.groupby(zone) + if weighting == 'area': + try: + running_sum += count * groupby.count().fillna(0) + except: + running_sum += Series([0]) + else: # weighting == 'population' + try: + running_sum += count * groupby.sum().fillna(0) + except: + running_sum += Series([0]) + if weighting == 'area': + rowresult = running_sum / totalpop_layer.groupby(zone).count() else: # weighting == 'population' - try: - running_sum += count * groupby.sum().fillna(0) - except: - running_sum += Series([0]) - if weighting == 'area': - rowresult = running_sum / totalpop_layer.groupby(zone).count() - else: # weighting == 'population' - rowresult = running_sum / totalpop_layer.groupby(zone).sum() + rowresult = running_sum / totalpop_layer.groupby(zone).sum() - results.append(rowresult[0]) + results.append(rowresult[0]) result_gdf = GeoDataFrame({f'{weighting}_averaged_num_accessible_amenities': results, 'geometry': zones['geometry']}).set_crs('EPSG:4326') return result_gdf \ No newline at end of file From 7b6a0ad1e44e6eaf78d8b699a0ebc9afbb4bce18 Mon Sep 17 00:00:00 2001 From: Ted Wong Date: Sun, 19 Jan 2025 09:57:14 -0500 Subject: [PATCH 25/25] clarifying comments --- city_metrix/metrics/count_accessible_amenities.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/city_metrix/metrics/count_accessible_amenities.py b/city_metrix/metrics/count_accessible_amenities.py index 5bb5839..ef9cf32 100644 --- a/city_metrix/metrics/count_accessible_amenities.py +++ b/city_metrix/metrics/count_accessible_amenities.py @@ -77,8 +77,8 @@ def get_data(self, bbox): dissected_gdf = GeoDataFrame({'poly_id': range(len(list(iso_dissected_polys.geoms))), 'geometry': list(iso_dissected_polys.geoms)}).set_crs('EPSG:4326') # For each dissected polygon, find how many of the original isoline polys contain the centroid # This is the number of amenity points are accessible within the dissected polygon - print("Counting how many original isoline polygons intersect with each dissected polygon") - count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, with num rows == num of dissected polys + print("Counting original isoline polygons that intersect with each dissected polygon") + count_amenities = dissected_gdf.centroid.within(iso_gdf.iloc[[0]].geometry[iso_gdf.index[0]]) * 1 # This is a Series storing running sum, a vector with num rows == num of dissected polys print("|==" + "PROGRESS=OF=INCLUSION=TESTS" + ("=" * 71) + "|") print(" ", end="") progress = 0