SEARCH_AREA_ID = 3600365331 # Italia
#SEARCH_AREA_ID = 3600042611 # Emilia-Romagna
#SEARCH_AREA_ID = 3600040218 # Campania
# ID calculated with https://wiki.openstreetmap.org/wiki/Overpass_API/Overpass_QL#By_element_id
PROVINCES_FILE_PATH = f"./provinces_{SEARCH_AREA_ID}.4326.parquet"
MUNICIPALITIES_FILE_PATH = f"./municipalities_{SEARCH_AREA_ID}.4326.parquet"
TOWNHALLS_FILE_PATH = f"./townhalls_{SEARCH_AREA_ID}.4326.parquet"
BAD_TOWNHALLS_FILE_PATH = f"./bad_townhalls_{SEARCH_AREA_ID}.4326.geojson"
WITHOUT_TOWNHALL_FILE_PATH = f"./without_townhall_{SEARCH_AREA_ID}.4326.geojson"
DBSN_FILE_PATH = f"./dbsn_townhalls.geojson"
DBSN_CONFLICTS_FILE_PATH = f"./dbsn_conflicts_{SEARCH_AREA_ID}.4326.geojson"
DBSN_MISSING_FILE_PATH = f"./dbsn_missing_{SEARCH_AREA_ID}.4326.geojson"
UNTAGGED_FILE_PATH = f"./untagged_{SEARCH_AREA_ID}.4326.geojson"
UNTAGGED_MISSING_FILE_PATH = f"./untagged_missing_{SEARCH_AREA_ID}.4326.geojson"
%pip install geopandas contextily pyproj rtree shapely mapclassify branca folium
Requirement already satisfied: geopandas in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (0.13.2) Requirement already satisfied: contextily in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (1.3.0) Requirement already satisfied: pyproj in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (3.6.0) Requirement already satisfied: rtree in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (1.0.1) Requirement already satisfied: shapely in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (2.0.1) Requirement already satisfied: mapclassify in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (2.6.0) Requirement already satisfied: branca in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (0.6.0) Requirement already satisfied: folium in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (0.14.0) Requirement already satisfied: fiona>=1.8.19 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from geopandas) (1.9.4.post1) Requirement already satisfied: packaging in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from geopandas) (23.0) Requirement already satisfied: pandas>=1.1.0 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from geopandas) (1.5.3) Requirement already satisfied: geopy in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from contextily) (2.3.0) Requirement already satisfied: matplotlib in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from contextily) (3.7.1) Requirement already satisfied: mercantile in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from contextily) (1.2.1) Requirement already satisfied: pillow in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from contextily) (9.4.0) Requirement already satisfied: rasterio in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from contextily) (1.3.8) Requirement already satisfied: requests in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from contextily) (2.31.0) Requirement already satisfied: joblib in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from contextily) (1.2.0) Requirement already satisfied: xyzservices in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from contextily) (2022.9.0) Requirement already satisfied: certifi in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from pyproj) (2023.7.22) Requirement already satisfied: numpy>=1.14 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from shapely) (1.24.3) Requirement already satisfied: scipy>=1.0 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from mapclassify) (1.10.1) Requirement already satisfied: scikit-learn in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from mapclassify) (1.3.0) Requirement already satisfied: networkx in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from mapclassify) (3.1) Requirement already satisfied: jinja2 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from branca) (3.1.2) Requirement already satisfied: attrs>=19.2.0 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (22.1.0) Requirement already satisfied: click~=8.0 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (8.0.4) Requirement already satisfied: click-plugins>=1.0 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (1.1.1) Requirement already satisfied: cligj>=0.5 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (0.7.2) Requirement already satisfied: six in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (1.16.0) Requirement already satisfied: MarkupSafe>=2.0 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from jinja2->branca) (2.1.1) Requirement already satisfied: python-dateutil>=2.8.1 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (2.8.2) Requirement already satisfied: pytz>=2020.1 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (2022.7) Requirement already satisfied: geographiclib<3,>=1.52 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from geopy->contextily) (2.0) Requirement already satisfied: contourpy>=1.0.1 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from matplotlib->contextily) (1.0.5) Requirement already satisfied: cycler>=0.10 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from matplotlib->contextily) (0.11.0) Requirement already satisfied: fonttools>=4.22.0 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from matplotlib->contextily) (4.25.0) Requirement already satisfied: kiwisolver>=1.0.1 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from matplotlib->contextily) (1.4.4) Requirement already satisfied: pyparsing>=2.3.1 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from matplotlib->contextily) (3.0.9) Requirement already satisfied: affine in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from rasterio->contextily) (2.4.0) Requirement already satisfied: snuggs>=1.4.1 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from rasterio->contextily) (1.4.7) Requirement already satisfied: setuptools in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from rasterio->contextily) (68.0.0) Requirement already satisfied: charset-normalizer<4,>=2 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from requests->contextily) (2.0.4) Requirement already satisfied: idna<4,>=2.5 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from requests->contextily) (3.4) Requirement already satisfied: urllib3<3,>=1.21.1 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from requests->contextily) (1.26.16) Requirement already satisfied: threadpoolctl>=2.0.0 in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from scikit-learn->mapclassify) (2.2.0) Requirement already satisfied: colorama in c:\users\daniele.santini\appdata\local\anaconda3\lib\site-packages (from click~=8.0->fiona>=1.8.19->geopandas) (0.4.6) Note: you may need to restart the kernel to use updated packages.
from pandas import merge, Series
from geopandas import GeoDataFrame, read_file, read_parquet
from shapely.geometry import GeometryCollection, shape, Point, MultiPolygon, LineString
from shapely.ops import polygonize
from urllib.request import urlopen, urlretrieve
from urllib.error import HTTPError
from urllib.parse import quote_plus
from os.path import exists
import contextily as cx
import json
from numpy import array
from os import system
PROVINCE_OVERPASS_QUERY=f"""
[out:json][timeout:90];
area({SEARCH_AREA_ID})->.searchArea;
(
relation["boundary"="administrative"]["admin_level"="6"]["ISO3166-2"!="FR-74"](area.searchArea);
relation["boundary"="administrative"]["admin_level"="4"]["ISO3166-2"="IT-23"](area.searchArea);
);
convert item ::=::,::geom=geom(),_osm_type=type(),_osm_id=id();
out geom;
"""
MUNICIPALITY_OVERPASS_QUERY=f"""
[out:json][timeout:90];
area({SEARCH_AREA_ID})->.searchArea;
relation["boundary"="administrative"]["admin_level"="8"](area.searchArea);
convert item ::=::,::geom=geom(),_osm_type=type(),_osm_id=id();
out geom;
"""
TOWNHALL_OVERPASS_QUERY=f"""
[out:json][timeout:90];
area({SEARCH_AREA_ID})->.searchArea;
nwr["amenity"="townhall"](area.searchArea);
convert item ::=::,::geom=geom(),_osm_type=type(),_osm_id=id();
out geom;
"""
def fetchOverpassGeoDataFrame(overpass_query:str, geometry_transform=shape):
url = "https://overpass-api.de/api/interpreter?data="+quote_plus(overpass_query)
try:
with urlopen(url) as response:
data = response.read()
encoding = response.info().get_content_charset('utf-8')
json_content = data.decode(encoding)
if "Query timed out" in json_content:
raise Exception("Query timed out")
#print(json_content)
json_object = json.loads(json_content)
#print(json_object['elements'][0] if json_object['elements'] else "No elments")
elements = [{
"id": element["id"],
"osm_id": element["tags"]["_osm_id"],
"osm_type": element["tags"]["_osm_type"],
"osm_url": 'https://www.openstreetmap.org/'+element["tags"]["_osm_type"]+'/'+element["tags"]["_osm_id"],
"name": element["tags"]["name"] if "name" in element["tags"] else None,
"geometry": geometry_transform(element['geometry'])
} for element in json_object['elements']]
#print(elements[0])
# OSM uses WGS 84: https://wiki.openstreetmap.org/wiki/Converting_to_WGS84
crs = 'EPSG:4326' # Use the SRID for WGS 84 - https://epsg.io/4326
gdf = GeoDataFrame(elements, crs=crs)
except HTTPError as err:
print("Failed downloading data from Overpass, retry later")
raise err
except json.JSONDecodeError as err:
print("Failed interpreting JSON data from Overpass")
raise err
return gdf
def readOrFetchOverpassGeoDataFrame(file_path:str, overpass_query:str, geometry_transform=shape):
if exists(file_path):
if file_path.endswith(".parquet"):
gdf = read_parquet(file_path)
else:
gdf = read_file(file_path, driver='GeoJSON')
else:
gdf = fetchOverpassGeoDataFrame(overpass_query, geometry_transform)
if file_path.endswith(".parquet"):
gdf.to_parquet(file_path)
else:
gdf.to_file(file_path, driver='GeoJSON')
return gdf
# Convert Overpass geometries into MultiPolygons - https://stackoverflow.com/a/72677231/2347196
convert_geom_to_multipolygon = lambda x: MultiPolygon(polygonize(shape(x)))
province_gdf = readOrFetchOverpassGeoDataFrame(PROVINCES_FILE_PATH, PROVINCE_OVERPASS_QUERY, convert_geom_to_multipolygon)
province_gdf.count()
id 110 osm_id 110 osm_type 110 osm_url 110 name 110 geometry 110 dtype: int64
province_gdf.head()
| id | osm_id | osm_type | osm_url | name | geometry | |
|---|---|---|---|---|---|---|
| 0 | 1 | 39151 | relation | https://www.openstreetmap.org/relation/39151 | Agrigento | MULTIPOLYGON (((13.98172 37.19299, 13.98217 37... |
| 1 | 2 | 39979 | relation | https://www.openstreetmap.org/relation/39979 | Nuoro | MULTIPOLYGON (((9.62535 40.25649, 9.62514 40.2... |
| 2 | 3 | 40021 | relation | https://www.openstreetmap.org/relation/40021 | Aristanis/Oristano | MULTIPOLYGON (((8.38245 40.33860, 8.38324 40.3... |
| 3 | 4 | 276369 | relation | https://www.openstreetmap.org/relation/276369 | Cagliari | MULTIPOLYGON (((9.16386 39.83254, 9.16450 39.8... |
| 4 | 5 | 12998776 | relation | https://www.openstreetmap.org/relation/12998776 | Sulcis Iglesiente | MULTIPOLYGON (((8.37422 39.02318, 8.37403 39.0... |
municipality_gdf = readOrFetchOverpassGeoDataFrame(MUNICIPALITIES_FILE_PATH, MUNICIPALITY_OVERPASS_QUERY, convert_geom_to_multipolygon)
municipality_gdf.count()
id 7903 osm_id 7903 osm_type 7903 osm_url 7903 name 7903 geometry 7903 dtype: int64
municipality_gdf.head()
| id | osm_id | osm_type | osm_url | name | geometry | |
|---|---|---|---|---|---|---|
| 0 | 1 | 39150 | relation | https://www.openstreetmap.org/relation/39150 | Lampedusa e Linosa | MULTIPOLYGON (((12.87805 35.85517, 12.87809 35... |
| 1 | 2 | 39777 | relation | https://www.openstreetmap.org/relation/39777 | Santu Antiogu/Sant'Antioco | MULTIPOLYGON (((8.38404 39.00591, 8.38379 39.0... |
| 2 | 3 | 39809 | relation | https://www.openstreetmap.org/relation/39809 | Câdesédda/Calasetta | MULTIPOLYGON (((8.37490 39.10915, 8.37532 39.1... |
| 3 | 4 | 39853 | relation | https://www.openstreetmap.org/relation/39853 | Igrèsias/Iglesias | MULTIPOLYGON (((8.43347 39.30784, 8.43347 39.3... |
| 4 | 5 | 39915 | relation | https://www.openstreetmap.org/relation/39915 | Bugerru/Buggerru | MULTIPOLYGON (((8.41075 39.44130, 8.41059 39.4... |
townhall_gdf = readOrFetchOverpassGeoDataFrame(TOWNHALLS_FILE_PATH, TOWNHALL_OVERPASS_QUERY)
townhall_gdf.count()
id 7313 osm_id 7313 osm_type 7313 osm_url 7313 name 5760 geometry 7313 dtype: int64
townhall_gdf.head()
| id | osm_id | osm_type | osm_url | name | geometry | |
|---|---|---|---|---|---|---|
| 0 | 1 | 4492704609 | node | https://www.openstreetmap.org/node/4492704609 | Comune di Carloforte | POINT (8.30562 39.14578) |
| 1 | 2 | 2440099045 | node | https://www.openstreetmap.org/node/2440099045 | None | POINT (8.37907 39.20365) |
| 2 | 3 | 1853454108 | node | https://www.openstreetmap.org/node/1853454108 | Comune di Teulada | POINT (8.77381 38.96790) |
| 3 | 4 | 5358907076 | node | https://www.openstreetmap.org/node/5358907076 | Comune di Sant'Antioco | POINT (8.45543 39.06645) |
| 4 | 5 | 2126087185 | node | https://www.openstreetmap.org/node/2126087185 | Comune di San Giovanni Suergiu | POINT (8.52207 39.11028) |
from matplotlib import pyplot as plt
def show_map(geo_df:GeoDataFrame, background_gdf:GeoDataFrame=None, color_column:str=None, cmap:str=None):
df_wm = geo_df.to_crs(epsg=3857)
figsize=(20,10)
fig,ax = plt.subplots(1, 1, figsize=figsize)
legend = False
if background_gdf is not None:
background_df_wm = background_gdf.to_crs(epsg=3857)
ax = background_df_wm.plot(ax=ax, figsize=figsize, alpha=0.3, edgecolor='k')
if color_column is not None:
legend = True
ax = df_wm.plot(ax=ax, figsize=figsize, edgecolor='k', column=color_column, cmap=cmap, legend=legend)
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLite)
cx.add_basemap(ax, source=cx.providers.Stamen.TonerLabels)
Map of municipalities available on OSM:
show_map(municipality_gdf)
Map of town halls available on OSM:
show_map(townhall_gdf)
#show_map(townhall_gdf, municipality_gdf)
#municipality_gdf.explore()
#townhall_gdf.explore()
with_townhall_gdf = townhall_gdf.sjoin(
municipality_gdf,
how="inner",
predicate="within",
lsuffix="hall",
rsuffix="town"
)
with_townhall_gdf.count()
id_hall 7319 osm_id_hall 7319 osm_type_hall 7319 osm_url_hall 7319 name_hall 5765 geometry 7319 index_town 7319 id_town 7319 osm_id_town 7319 osm_type_town 7319 osm_url_town 7319 name_town 7319 dtype: int64
without_townhall_gdf = municipality_gdf[ # Anti-join
~municipality_gdf["id"].isin(with_townhall_gdf["id_town"])
]
without_townhall_gdf.to_file(WITHOUT_TOWNHALL_FILE_PATH, driver="GeoJSON")
without_townhall_gdf.count()
id 1231 osm_id 1231 osm_type 1231 osm_url 1231 name 1231 geometry 1231 dtype: int64
without_townhall_gdf.head()
| id | osm_id | osm_type | osm_url | name | geometry | |
|---|---|---|---|---|---|---|
| 27 | 28 | 39855 | relation | https://www.openstreetmap.org/relation/39855 | Bidda Matzràxia/Villamassargia | MULTIPOLYGON (((8.73100 39.28096, 8.73147 39.2... |
| 44 | 45 | 39803 | relation | https://www.openstreetmap.org/relation/39803 | Sa Baronia/Villaperuccio | MULTIPOLYGON (((8.61680 39.12812, 8.61717 39.1... |
| 185 | 186 | 12338682 | relation | https://www.openstreetmap.org/relation/12338682 | Misiliscemi | MULTIPOLYGON (((12.49852 37.96577, 12.50396 37... |
| 426 | 427 | 40827 | relation | https://www.openstreetmap.org/relation/40827 | Caragnani/Calangianus | MULTIPOLYGON (((9.30121 40.84693, 9.30056 40.8... |
| 453 | 454 | 40950 | relation | https://www.openstreetmap.org/relation/40950 | Sant'Antoni di Gaddura/Sant'Antonio di Gallura | MULTIPOLYGON (((9.31750 40.92267, 9.31720 40.9... |
Map of municipalities without town hall:
show_map(without_townhall_gdf)
#without_townhall_gdf.explore()
province_municipality_df = province_gdf[["id","geometry"]].sjoin(
municipality_gdf[["geometry"]],
how="left",
predicate="contains",
lsuffix="pro",
rsuffix="mun"
)
province_gdf["num_municipalities"] = province_municipality_df.groupby(province_municipality_df.index).count()["index_mun"]
province_without_townhall_df = province_gdf[["id","geometry"]].sjoin(
without_townhall_gdf[["geometry"]],
how="left",
predicate="contains",
lsuffix="pro",
rsuffix="mun"
)
province_gdf["num_without_townhall"] = province_without_townhall_df.groupby(province_without_townhall_df.index).count()["index_mun"]
province_gdf["osm_availability"] = 100 * (1 - (province_gdf["num_without_townhall"] / province_gdf["num_municipalities"]))
province_gdf["osm_availability"].hist()
<Axes: >
show_map(province_gdf, None, "osm_availability", "RdYlGn")
UNTAGGED_OVERPASS_QUERY=f"""
[out:json][timeout:300];
area({SEARCH_AREA_ID})->.searchArea;
nwr["building"]["amenity"!="townhall"]["name"~"^(\s|palazzo|del|nuovo|comune|-)*municipio",i](area.searchArea);
convert item ::=::,::geom=geom(),_osm_type=type(),_osm_id=id();
out geom;
"""
untagged_gdf = readOrFetchOverpassGeoDataFrame(UNTAGGED_FILE_PATH, UNTAGGED_OVERPASS_QUERY)
untagged_gdf.count()
id 30 osm_id 30 osm_type 30 osm_url 30 name 30 geometry 30 dtype: int64
untagged_gdf.head()
| id | osm_id | osm_type | osm_url | name | geometry | |
|---|---|---|---|---|---|---|
| 0 | 1 | 277834707 | way | https://www.openstreetmap.org/way/277834707 | Municipio di Monserrato | LINESTRING (9.14363 39.25416, 9.14336 39.25394... |
| 1 | 2 | 933053288 | way | https://www.openstreetmap.org/way/933053288 | Municipio | LINESTRING (13.09184 37.95253, 13.09212 37.952... |
| 2 | 3 | 332108931 | way | https://www.openstreetmap.org/way/332108931 | Municipio Roma V - Sala del Consiglio | LINESTRING (12.54286 41.88062, 12.54275 41.880... |
| 3 | 4 | 60145522 | way | https://www.openstreetmap.org/way/60145522 | Municipio | LINESTRING (12.95500 41.58178, 12.95523 41.581... |
| 4 | 5 | 1141570499 | way | https://www.openstreetmap.org/way/1141570499 | Municipio di Montecastrilli | LINESTRING (12.48795 42.64981, 12.48778 42.649... |
untagged_missing_gdf = untagged_gdf.sjoin(
without_townhall_gdf,
how="inner",
predicate="within",
lsuffix="hall",
rsuffix="town"
)
untagged_missing_gdf.count()
id_hall 2 osm_id_hall 2 osm_type_hall 2 osm_url_hall 2 name_hall 2 geometry 2 index_town 2 id_town 2 osm_id_town 2 osm_type_town 2 osm_url_town 2 name_town 2 dtype: int64
untagged_missing_gdf.head()
| id_hall | osm_id_hall | osm_type_hall | osm_url_hall | name_hall | geometry | index_town | id_town | osm_id_town | osm_type_town | osm_url_town | name_town | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 4 | 5 | 1141570499 | way | https://www.openstreetmap.org/way/1141570499 | Municipio di Montecastrilli | LINESTRING (12.48795 42.64981, 12.48778 42.649... | 874 | 875 | 42101 | relation | https://www.openstreetmap.org/relation/42101 | Montecastrilli |
| 10 | 11 | 672589063 | way | https://www.openstreetmap.org/way/672589063 | Municipio di Atella di Napoli (1928-46) | LINESTRING (14.26055 40.96097, 14.26056 40.961... | 1834 | 1835 | 40939 | relation | https://www.openstreetmap.org/relation/40939 | Sant'Arpino |
untagged_missing_gdf.to_file(UNTAGGED_MISSING_FILE_PATH, driver='GeoJSON')
#show_map(untagged_missing_gdf)
untagged_missing_gdf.explore()
bad_townhall_mask = Series(False, index=townhall_gdf.index)
for word in ["consorzio","uffici"]:
bad_townhall_mask |= townhall_gdf["name"].str.lower().str.contains(word, na=False)
bad_townhall_mask.value_counts()
False 7274 True 39 dtype: int64
bad_townhall_gdf = townhall_gdf[bad_townhall_mask]
bad_townhall_gdf.to_file(BAD_TOWNHALLS_FILE_PATH, driver='GeoJSON')
bad_townhall_gdf.explore()
Download the data from IGM DBSN ( https://www.igmi.org/it/dbsn-database-di-sintesi-nazionale ).
Then extract the catageory of elements you are interested in and save it as GeoJSON with EPSG 4326 SRID. To accomplish you can use the script "filtra_dbsn.sh" in the project repository root (or https://www.dsantini.it/dbsn/filtra_dbsn.sh ). The file downloaded below from dbsn_url has been generated with this script.
def download_file_if_not_exists(file_path, url):
if not exists(file_path):
try:
urlretrieve(url, file_path)
except HTTPError as err:
print("Failed downloading data from Overpass, retry later")
raise err
# https://www.dsantini.it/dbsn/
dbsn_url = "https://www.dsantini.it/dbsn/notebooks/dbsn_townhalls.geojson"
download_file_if_not_exists(DBSN_FILE_PATH, dbsn_url)
dbsn_gdf = read_file(DBSN_FILE_PATH)
dbsn_gdf.count()
edifc_uso 8382 edifc_ty 8382 edifc_sot 8382 classid 8354 edifc_nome 8381 edifc_stat 8382 edifc_at 8321 scril 8055 meta_ist 8367 edifc_mon 8382 shape_Length 8382 shape_Area 8382 geometry 8382 dtype: int64
dbsn_gdf.head()
| edifc_uso | edifc_ty | edifc_sot | classid | edifc_nome | edifc_stat | edifc_at | scril | meta_ist | edifc_mon | shape_Length | shape_Area | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0201 | 93 | 93 | 65D8ED81-5047-4CEA-9103-AA2034B665DD | Municipio di Lampedusa e Linosa | 91 | -9999.0 | 10k | 04 | 02 | 66.022667 | 272.435503 | MULTIPOLYGON Z (((12.60941 35.50254 0.00000, 1... |
| 1 | 0201 | 93 | 93 | D659C260-EDD5-4540-BC6C-EACB33EC447B | Municipio di Licata | 91 | -9999.0 | 10k | 04 | 02 | 99.677476 | 358.678321 | MULTIPOLYGON Z (((13.93720 37.10145 0.00000, 1... |
| 2 | 0201 | 93 | 93 | b4949764-f626-442c-89ad-200e56ea653d | Municipio di Cammarata | 91 | -9999.0 | 10k | 04 | 02 | 133.927947 | 670.469227 | MULTIPOLYGON Z (((13.63779 37.63324 0.00000, 1... |
| 3 | 0201 | 93 | 93 | ece5bd6f-834b-4230-a100-89232be50ff0 | Municipio di San Biagio Platani | 91 | -9999.0 | 10k | 04 | 02 | 134.259563 | 628.250694 | MULTIPOLYGON Z (((13.52434 37.50919 0.00000, 1... |
| 4 | 0201 | 93 | 93 | 7c14edd9-e910-4de2-9cb9-2cc44be6a4fc | Municipio di Cianciana | 91 | -9999.0 | 10k | 04 | 02 | 126.859157 | 633.197564 | MULTIPOLYGON Z (((13.43422 37.51843 0.00000, 1... |
municipality_dbsn_gdf = municipality_gdf[["id","geometry"]].sjoin(
dbsn_gdf[["geometry"]],
how="left",
predicate="contains",
lsuffix="mun",
rsuffix="dbsn"
)
municipality_gdf["dbsn_townhall_num"] = municipality_dbsn_gdf.groupby(municipality_dbsn_gdf.index).count()["index_dbsn"]
municipality_gdf["dbsn_townhall_num"].hist()
<Axes: >
show_map(municipality_gdf[municipality_gdf["dbsn_townhall_num"] != 1], None, "dbsn_townhall_num", "RdYlGn")
province_dbsn_df = province_gdf[["id","geometry"]].sjoin(
dbsn_gdf[["geometry"]],
how="left",
predicate="contains",
lsuffix="pro",
rsuffix="dbsn"
)
province_gdf["dbsn_townhall_num"] = province_dbsn_df.groupby(province_dbsn_df.index).count()["index_dbsn"]
province_gdf["dbsn_availability"] = (100 * (province_gdf["dbsn_townhall_num"] / province_gdf["num_municipalities"])).astype(int)
province_gdf["dbsn_availability"].hist()
<Axes: >
show_map(province_gdf[province_gdf["dbsn_availability"] != 100], None, "dbsn_availability", "RdYlGn")
province_gdf[["id","name","dbsn_townhall_num","num_municipalities","dbsn_availability"]].to_csv("province_stats.csv")
dbsn_existing_df = dbsn_gdf.sjoin(
municipality_gdf,
how="inner",
predicate="within",
lsuffix="dbsn",
rsuffix="town"
)
dbsn_existing_df.count()
edifc_uso 8374 edifc_ty 8374 edifc_sot 8374 classid 8360 edifc_nome 8373 edifc_stat 8374 edifc_at 8327 scril 8060 meta_ist 8373 edifc_mon 8374 shape_Length 8374 shape_Area 8374 geometry 8374 index_town 8374 id 8374 osm_id 8374 osm_type 8374 osm_url 8374 name 8374 dbsn_townhall_num 8374 dtype: int64
compare_gdf = merge(
dbsn_existing_df,
with_townhall_gdf,
how="inner",
left_on="index_town",
right_on="index_town"
)
compare_gdf.count()
edifc_uso 7838 edifc_ty 7838 edifc_sot 7838 classid 7829 edifc_nome 7837 edifc_stat 7838 edifc_at 7794 scril 7496 meta_ist 7833 edifc_mon 7838 shape_Length 7838 shape_Area 7838 geometry_x 7838 index_town 7838 id 7838 osm_id 7838 osm_type 7838 osm_url 7838 name 7838 dbsn_townhall_num 7838 id_hall 7838 osm_id_hall 7838 osm_type_hall 7838 osm_url_hall 7838 name_hall 6279 geometry_y 7838 id_town 7838 osm_id_town 7838 osm_type_town 7838 osm_url_town 7838 name_town 7838 dtype: int64
compare_gdf["distance"] = compare_gdf["geometry_x"].to_crs(epsg=3857).distance(
compare_gdf["geometry_y"].to_crs(epsg=3857)
)
compare_gdf["geometry"] = compare_gdf.apply(
lambda row: LineString([row['geometry_x'].centroid, row['geometry_y'].centroid]),
axis=1
).set_crs(epsg=4326)
compare_gdf.hist(log=True, column="distance")
array([[<Axes: title={'center': 'distance'}>]], dtype=object)
conflict_mask=compare_gdf["distance"] > 1000
for word in [
"delegazione",
"circoscrizione",
"quartiere",
"frazione",
"municipalit",
"municipio roma",
"consorzio",
"° municipio",
"uffici",
"sede di",
"distaccat"
]:
conflict_mask &= ~(compare_gdf["name_hall"].str.lower().str.contains(word, na=False))
conflict_mask.value_counts()
False 7489 True 349 Name: distance, dtype: int64
dbsn_conflict_gdf = compare_gdf[conflict_mask]
dbsn_conflict_gdf["name_hall"].dropna().drop_duplicates()
87 Municipio di Alluvioni Piovera
88 Municipio di Piovera
131 Municipio di Lu e Cuccaro Monferrato
408 Comune dell'Aquila Ricostruzione Beni Pubblici
410 Comune dell'Aquila Equità tributaria
...
7582 Sede Municipale di Colceresa
7606 Ex municipio di Campolongo sul Brenta
7617 Municipio di Lusiana Conco
7669 Comune di Verona - Edilizia Scolastica
7784 Comune di Ionadi
Name: name_hall, Length: 177, dtype: object
dbsn_conflict_gdf.explore()
dbsn_conflict_gdf.drop(["geometry_x", "geometry_y"], axis=1).to_file(DBSN_CONFLICTS_FILE_PATH, driver="GeoJSON")
dbsn_missing_df = dbsn_gdf.sjoin(
without_townhall_gdf,
how="inner",
predicate="within"
)
dbsn_missing_df.count()
edifc_uso 1278 edifc_ty 1278 edifc_sot 1278 classid 1273 edifc_nome 1278 edifc_stat 1278 edifc_at 1273 scril 1255 meta_ist 1278 edifc_mon 1278 shape_Length 1278 shape_Area 1278 geometry 1278 index_right 1278 id 1278 osm_id 1278 osm_type 1278 osm_url 1278 name 1278 dtype: int64
dbsn_missing_df.head()
| edifc_uso | edifc_ty | edifc_sot | classid | edifc_nome | edifc_stat | edifc_at | scril | meta_ist | edifc_mon | shape_Length | shape_Area | geometry | index_right | id | osm_id | osm_type | osm_url | name | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 48 | 0201 | 01 | 01 | e6635754-4083-4030-8ceb-77084f2c9ba3 | Municipio di PADERNA | 03 | -9999.0 | 1:2000 ... | 04 | 02 | 76.283735 | 210.957406 | MULTIPOLYGON Z (((8.89020 44.82071 0.00000, 8.... | 3852 | 3853 | 43685 | relation | https://www.openstreetmap.org/relation/43685 | Paderna |
| 49 | 0201 | 01 | 01 | 93480b43-8fd4-49be-a7f0-c85487b62546 | Municipio di MONTEGIOCO | 03 | -9999.0 | 1:2000 ... | 04 | 02 | 45.669440 | 129.941911 | MULTIPOLYGON Z (((8.95418 44.83776 0.00000, 8.... | 3839 | 3840 | 43687 | relation | https://www.openstreetmap.org/relation/43687 | Montegioco |
| 50 | 0201 | 01 | 01 | 5d77b0ec-1cc8-491e-8a92-425b38cb3111 | Municipio di CASTELLAR GUIDOBONO | 03 | -9999.0 | 1:2000 ... | 04 | 02 | 46.584363 | 123.354444 | MULTIPOLYGON Z (((8.94677 44.90624 0.00000, 8.... | 3857 | 3858 | 43791 | relation | https://www.openstreetmap.org/relation/43791 | Castellar Guidobono |
| 52 | 0201 | 01 | 01 | 94106d54-3dd9-422e-9199-93b079e01f6c | Municipio di VILLALVERNIA | 03 | -9999.0 | 1:2000 ... | 04 | 02 | 92.197220 | 387.012599 | MULTIPOLYGON Z (((8.85671 44.81768 0.00000, 8.... | 3837 | 3838 | 43637 | relation | https://www.openstreetmap.org/relation/43637 | Villalvernia |
| 53 | 0201 | 01 | 01 | 287da1e2-505b-4295-b484-9e3ec07571c5 | Municipio di DERNICE | 03 | -9999.0 | 1:2000 ... | 04 | 02 | 48.885326 | 84.226670 | MULTIPOLYGON Z (((9.04985 44.76721 0.00000, 9.... | 3830 | 3831 | 43572 | relation | https://www.openstreetmap.org/relation/43572 | Dernice |
dbsn_missing_df.to_file(DBSN_MISSING_FILE_PATH, driver="GeoJSON")
#%pip install folium matplotlib mapclassify
dbsn_missing_df.explore()