Solution 04

Setup

try:
  import pygeos
except ModuleNotFoundError as e:
  !pip install pygeos==0.10.2
  import pygeos

try:
  import geopandas as gpd
except ModuleNotFoundError as e:
  !pip install geopandas==0.10.0
  import geopandas as gpd

if gpd.__version__ != "0.10.0":
  !pip install -U geopandas==0.10.0
  import geopandas as gpd

try:
  import pyrosm
except ModuleNotFoundError as e:
  !pip install pyrosm==0.6.1
  import pyrosm

Exercise

download from OpenStreetMap all supermarkets inside the bounding box of the city in this point
latitude: 45.4395
longitude: 11.09351
identify the longest road of the city (state roads, walking routes, motorways are excluded). Please use “unclassified”
How many drinking water are in this city?
How many benches in this city have the backrest?

download from OpenStreetMap all supermarkets inside the bounding box of the city in this point and

latitude: 45.4395
longitude: 11.09351

import geopandas as gpd

latitude = 45.4395
longitude = 12.3478

find the city

method 1 - reverse geocoding

from geopy.geocoders import ArcGIS

latlon = str(latitude) + "," + str(longitude);

geolocator = ArcGIS(user_agent="Mozilla/5.0 (Linux; Android 10) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/86.0.4240.75 Mobile Safari/537.36")

location = geolocator.reverse(latlon)

location.raw

{'Match_addr': '30122',
 'LongLabel': '30122, ITA',
 'ShortLabel': '30122',
 'Addr_type': 'Postal',
 'Type': '',
 'PlaceName': '30122',
 'AddNum': '',
 'Address': '',
 'Block': '',
 'Sector': '',
 'Neighborhood': '',
 'District': '',
 'City': 'Venezia',
 'MetroArea': '',
 'Subregion': 'Venezia',
 'Region': 'Veneto',
 'RegionAbbr': '',
 'Territory': '',
 'Postal': '30122',
 'PostalExt': '',
 'CntryName': 'Italia',
 'CountryCode': 'ITA'}

city = location.raw['City']

city

'Venezia'

method 2 - spatial relation

from shapely.geometry import Point

point = Point(longitude, latitude)

urlmunicipalities2022 = 'https://github.com/napo/geospatial_course_unitn/raw/master/data/istat/shapefile_istat_municipalities_2022.zip'
municipalities2022 = gpd.read_file(urlmunicipalities2022)

muncipality = municipalities2022[municipalities2022.to_crs(epsg=4326).geometry.contains(point)]

muncipality

	COD_RIP	COD_REG	COD_PROV	COD_CM	COD_UTS	PRO_COM	PRO_COM_T	COMUNE	COMUNE_A	CC_UTS	Shape_Area	Shape_Leng	geometry
3692	2	5	27	227	227	27042	027042	Venezia	None	1	4.158927e+08	170837.245485	POLYGON ((780145.784 5049170.898, 780115.597 5...

city = muncipality.COMUNE.values[0]

city

'Venezia'

download from OpenStreetMap

find the boundig box of Venice

municipality = municipalities2022[municipalities2022.COMUNE==city]

Donwload the PBF of Venice from OSM Estratti

url_download_venice_pbf = 'https://osmit-estratti-test.wmcloud.org/dati/poly/comuni/pbf/027042_Venezia_poly.osm.pbf'
import urllib.request
urllib.request.urlretrieve(url_download_venice_pbf ,"venice_osm.pbf")

('venice_osm.pbf', <http.client.HTTPMessage at 0x7efc48639960>)

find all the supermarkets in the area

import pyrosm

osm = pyrosm.OSM("venice_osm.pbf") 

custom_filter = {'shop': ['supermarket']}

supermarkets = osm.get_pois(custom_filter=custom_filter)

/home/napo/.local/lib/python3.10/site-packages/pyrosm/pyrosm.py:576: ShapelyDeprecationWarning: The array interface is deprecated and will no longer work in Shapely 2.0. Convert the '.coords' to a numpy array instead.
  gdf = get_poi_data(
/home/napo/.local/lib/python3.10/site-packages/geopandas/array.py:1406: UserWarning: CRS not set for some of the concatenation inputs. Setting output's CRS as WGS 84 (the single non-null crs provided).

supermarkets.shape

(102, 23)

supermarkets.head(5)

	lon	version	tags	id	lat	timestamp	addr:city	addr:country	addr:housenumber	...	email	name	opening_hours	operator	phone	website	organic	shop	geometry	osm_type
0	12.327674	8	{"addr:neighbourhood":"Santa Croce","brand":"C...	443122709	45.440384	1628972013	None	None	1491a	...	None	Coop	Mo-Sa 08:30-20:00; Su 09:00-20:00	None	+39 041 2750218	None	None	supermarket	POINT (12.32767 45.44038)	node
1	12.200390	2	None	566696229	45.482265	1258587411	None	None	None	...	None	Alì	None	None	None	None	None	supermarket	POINT (12.20039 45.48227)	node
2	12.208943	1	None	566696230	45.485161	1258568581	None	None	None	...	None	Cadoro	None	None	None	None	None	supermarket	POINT (12.20894 45.48516)	node
3	12.253382	4	None	617964249	45.492580	1492417346	None	None	43a	...	None	Spak	None	None	+39 041 5351288	None	None	supermarket	POINT (12.25338 45.49258)	node
4	12.324074	10	{"brand":"Conad City","brand:wikidata":"Q57543...	626923862	45.433987	1603553460	Venezia	None	3017	...	None	Conad City	Mo-Sa 07:30-20:30; Su 09:00-14:00,16:00-19:30	Conad	None	https://www.conad.it	None	supermarket	POINT (12.32407 45.43399)	node

5 rows × 23 columns

print("there are %s supermarkets in Venice" % supermarkets.shape[0])

there are 102 supermarkets in Venice

supermarkets.name.unique()

array(['Coop', 'Alì', 'Cadoro', 'Spak', 'Conad City', 'Simply', 'Prix',
       'Spazio Conad', 'Despar', 'inCoop', 'Supermercato SISA',
       'Minimarket', None, 'Pam', 'Supermarket "Cadoro"', 'Lidl',
       'Camping Fusina', 'Eurospesa', 'Pam Local Corso del Popolo',
       'Supermercati Maxi', 'Punto Simply', 'Ali Supermercato',
       'Cuore Bio', 'Crai', 'NaturaSì', 'InCoop', 'Alì via Sforza',
       'Alì Supermercati', "In's Mercato", 'A&O', 'L’angolo Giallo',
       'Ballarin Paolo E C. S.N.C.', 'Dolce Bianchi Venezia', 'Interspar',
       'Conad', 'Tigota', 'Aldi', 'Aliper', 'Mega Marghera', "iN's",
       'Eurospar', 'EuroSpin', 'Famila Superstore', 'IperLando'],
      dtype=object)

identify the longest road of the city

roads = osm.get_network(network_type='all')

roads.columns

Index(['access', 'area', 'bicycle', 'bridge', 'cycleway', 'est_width', 'foot',
       'footway', 'highway', 'int_ref', 'junction', 'lanes', 'lit', 'maxspeed',
       'motorcar', 'motor_vehicle', 'name', 'oneway', 'overtaking', 'path',
       'psv', 'ref', 'service', 'segregated', 'sidewalk', 'smoothness',
       'surface', 'tracktype', 'tunnel', 'width', 'id', 'timestamp', 'version',
       'tags', 'osm_type', 'geometry', 'length'],
      dtype='object')

roads.highway.unique()

array(['pedestrian', 'service', 'unclassified', 'residential', 'footway',
       'tertiary', 'primary', 'cycleway', 'steps', 'secondary',
       'motorway_link', 'track', 'secondary_link', 'primary_link',
       'motorway', 'path', 'tertiary_link', 'living_street',
       'construction', 'rest_area', 'services', 'elevator', 'corridor',
       'bridleway', 'proposed'], dtype=object)

lenght of unclassified and residential togheter

roads_unclassified_residential = roads[(roads.highway=='unclassified') |  (roads.highway == 'residential')]

roads_unclassified_residential

	access	area	bicycle	bridge	cycleway	est_width	foot	footway	highway	int_ref	...	tracktype	tunnel	width	id	timestamp	version	tags	osm_type	geometry	length
2	None	None	None	None	None	None	None	None	unclassified	None	...	None	None	None	4438335	1611572910	9	None	way	MULTILINESTRING ((12.38501 45.42321, 12.38511 ...	538.0
3	None	None	None	None	None	None	None	None	unclassified	None	...	None	None	None	4438345	1611572458	6	None	way	MULTILINESTRING ((12.38398 45.42382, 12.38346 ...	51.0
4	None	None	None	None	None	None	None	None	residential	None	...	None	None	None	4448298	1399909744	9	None	way	MULTILINESTRING ((12.36842 45.40751, 12.36896 ...	168.0
8	None	None	None	None	None	None	None	None	unclassified	None	...	None	None	None	4597767	1643983716	17	{"cycleway:right":"no","layer":"1","noname":"y...	way	MULTILINESTRING ((12.31524 45.44175, 12.31523 ...	39.0
13	destination	None	None	None	None	None	None	None	residential	None	...	None	None	None	4710516	1467135189	5	{"note":"durante i giorni di scuola accesso ri...	way	MULTILINESTRING ((12.27803 45.50600, 12.27798 ...	67.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
21153	None	None	None	None	None	None	None	None	residential	None	...	None	None	None	1054635814	1650871543	1	None	way	MULTILINESTRING ((12.24553 45.50209, 12.24560 ...	32.0
21154	None	None	None	None	None	None	None	None	residential	None	...	None	None	None	1054635815	1650871543	1	None	way	MULTILINESTRING ((12.24561 45.50221, 12.24556 ...	17.0
21164	None	None	None	None	None	None	None	None	unclassified	None	...	None	None	None	1055712580	1651150935	3	{"cycleway:both":"no"}	way	MULTILINESTRING ((12.29539 45.53549, 12.29539 ...	8.0
21165	None	None	None	None	None	None	None	None	unclassified	None	...	None	None	None	1055712581	1651150947	2	{"cycleway:both":"no"}	way	MULTILINESTRING ((12.29430 45.53416, 12.29319 ...	87.0
21177	None	None	None	None	None	None	None	None	residential	None	...	None	None	None	1058284862	1652007383	1	{"vehicle:conditional":"private @ (22:00-04:00)"}	way	MULTILINESTRING ((12.23248 45.48717, 12.23248 ...	137.0

4022 rows × 37 columns

roads_unclassified_residential['name'].value_counts()

Via Eugenio Carlo Pertini    35
Via Ca' Solaro               23
Via Vallon                   23
Via Colombara                23
Via Evaristo Scaramuzza      21
                             ..
Via Prete Ermanno             1
Via dei Gradenighi            1
Via dei Mandarini             1
Via Egidio Marcon             1
Santa Maria dei Battuti       1
Name: name, Length: 1461, dtype: int64

names = roads_unclassified_residential.name.unique()

names

array(['Viale Umberto Klinger', 'Via Giannantonio Selva',
       'Via Dardanelli', ..., 'via Caravaggio', 'Via Trento',
       'Santa Maria dei Battuti'], dtype=object)

The best projection for Venice is WGS84 UTM 33N - epsg:32633

roads_lenght = {}
rodas_in_meters = roads_unclassified_residential.to_crs(epsg=32633)
for name in names:
  road = rodas_in_meters[rodas_in_meters.name==name]
  road_lenght = road.length.sum()
  roads_lenght[name] = road_lenght

#roads_lenght

max(roads_lenght, key=roads_lenght.get)

'Via Litomarino'

longest_road = max(roads_lenght, key=roads_lenght.get)

roads_lenght[longest_road]

4882.770029371328

print("the longest road (unclassified + residential) in Venice is %s and is long %.1f km" % (str(longest_road),roads_lenght[longest_road]/1000))

the longest road (unclassified + residential) in Venice is Via Litomarino and is long 4.9 km

the case only with unclassified

roads_unclassified = roads[(roads.highway=='unclassified')] 

names = roads_unclassified.name.unique()

roads_lenght = {}
rodas_in_meters = roads_unclassified.to_crs(epsg=32632)
for name in names:
  road = rodas_in_meters[rodas_in_meters.name==name]
  road_lenght = road.length.sum()
  roads_lenght[name] = road_lenght

longest_road = max(roads_lenght, key=roads_lenght.get)

print("the longest road (unclassified) is %s and is long %.1f km" % (str(longest_road),roads_lenght[longest_road]/1000))

the longest road (unclassified) is Via Litomarino and is long 4.9 km

where is this road?

… we can use overpass-turbo.eu with this wizard query

name='Via Litomarino in Venezia'

https://overpass-turbo.eu/s/1nf9

How many drinking water are in this city?

custom_filter = {'amenity': ['drinking_water']}

drinking_water = osm.get_pois(custom_filter=custom_filter)

drinking_water.amenity.unique()

array(['drinking_water'], dtype=object)

drinking_water.shape

(181, 15)

drinking_water.head(4)

	lon	version	tags	id	lat	timestamp	addr:street	name	amenity	drinking_water	fountain	source	geometry	osm_type
0	12.234090	2	{"check_date":"2022-03-26"}	315040103	45.534821	1648284842	None	None	drinking_water	None	None	None	POINT (12.23409 45.53482)	node
1	12.324445	5	None	639811824	45.437145	1493048539	None	None	drinking_water	None	None	None	POINT (12.32444 45.43715)	node
2	12.286287	3	{"access":"private"}	699811494	45.520275	1612044310	None	None	drinking_water	None	None	None	POINT (12.28629 45.52028)	node
3	12.315189	2	None	1023518847	45.427383	1291582192	None	None	drinking_water	None	None	None	POINT (12.31519 45.42738)	node

drinking_water.explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

How many benches in this city have the backrest?

custom_filter = {'amenity': ['bench']}

benchs = osm.get_pois(custom_filter=custom_filter)

benchs.shape

(597, 11)

benchs.columns

Index(['lon', 'version', 'tags', 'id', 'changeset', 'lat', 'timestamp',
       'amenity', 'source', 'geometry', 'osm_type'],
      dtype='object')

benchs.tags.unique()

array([None, '{"backrest":"yes"}',
       '{"backrest":"yes","check_date":"2021-02-15"}',
       '{"backrest":"yes","note":"still exists?"}',
       '{"backrest":"yes","check_date":"2021-05-31","material":"wood","natural":"tree"}',
       '{"backrest":"yes","check_date":"2021-05-31","material":"wood"}',
       '{"backrest":"yes","direction":"110","material":"wood"}',
       '{"backrest":"yes","check_date":"2022-04-06"}',
       '{"backrest":"no","material":"stone"}',
       '{"backrest":"yes","check_date":"2021-05-27"}',
       '{"backrest":"no"}',
       '{"backrest":"yes","direction":"290","material":"wood"}',
       '{"backrest":"yes","direction":"200","material":"wood"}',
       '{"backrest":"yes","direction":"20","material":"wood"}',
       '{"backrest":"yes","colour":"red","material":"wood","seats":"3","survey:date":"2021-08-25"}',
       '{"backrest":"yes","colour":"red","seats":"3","survey:date":"2021-08-25"}',
       '{"backrest":"yes","direction":"130"}',
       '{"backrest":"yes","direction":"35"}',
       '{"backrest":"yes","direction":"215"}',
       '{"backrest":"yes","material":"metal","seats":"1"}',
       '{"backrest":"yes","material":"wood","seats":"3"}',
       '{"backrest":"yes","colour":"green","material":"metal","seats":"4"}',
       '{"backrest":"yes","colour":"red","material":"wood"}',
       '{"backrest":"no","survey:date":"2021-12-17"}',
       '{"backrest":"yes","survey:date":"2021-12-17"}',
       '{"backrest":"yes","colour":"brown","material":"wood"}',
       '{"backrest":"yes","colour":"red","material":"wood","shop":"newsagent"}',
       '{"backrest":"yes","colour":"Rosso","material":"wood","seats":"3"}'],
      dtype=object)

quick&dirty solution

benchs_with_backrest = benchs[benchs.tags.str.find('"backrest":"yes"') > -1]

benchs_with_backrest.tags.unique()

array(['{"backrest":"yes"}',
       '{"backrest":"yes","check_date":"2021-02-15"}',
       '{"backrest":"yes","note":"still exists?"}',
       '{"backrest":"yes","check_date":"2021-05-31","material":"wood","natural":"tree"}',
       '{"backrest":"yes","check_date":"2021-05-31","material":"wood"}',
       '{"backrest":"yes","direction":"110","material":"wood"}',
       '{"backrest":"yes","check_date":"2022-04-06"}',
       '{"backrest":"yes","check_date":"2021-05-27"}',
       '{"backrest":"yes","direction":"290","material":"wood"}',
       '{"backrest":"yes","direction":"200","material":"wood"}',
       '{"backrest":"yes","direction":"20","material":"wood"}',
       '{"backrest":"yes","colour":"red","material":"wood","seats":"3","survey:date":"2021-08-25"}',
       '{"backrest":"yes","colour":"red","seats":"3","survey:date":"2021-08-25"}',
       '{"backrest":"yes","direction":"130"}',
       '{"backrest":"yes","direction":"35"}',
       '{"backrest":"yes","direction":"215"}',
       '{"backrest":"yes","material":"metal","seats":"1"}',
       '{"backrest":"yes","material":"wood","seats":"3"}',
       '{"backrest":"yes","colour":"green","material":"metal","seats":"4"}',
       '{"backrest":"yes","colour":"red","material":"wood"}',
       '{"backrest":"yes","survey:date":"2021-12-17"}',
       '{"backrest":"yes","colour":"brown","material":"wood"}',
       '{"backrest":"yes","colour":"red","material":"wood","shop":"newsagent"}',
       '{"backrest":"yes","colour":"Rosso","material":"wood","seats":"3"}'],
      dtype=object)

benchs_with_backrest.shape

(279, 11)

print("the total of benchs tagged with the presence of a backrest is %s" % (str(benchs_with_backrest.shape[0])))

the total of benchs tagged with the presence of a backrest is 279

Share on

Twitter Facebook LinkedIn

Geospatial Lab