Solution 02
Exercise 02: Spatial Relationships and Operations
learning objectives
- repeat the concepts on the previous lesson
- errors with the simplified boundaries
- convex hull / concave hull / alphashape
- nearest_points
Exercise
1 - create the geodataframe of the gas&oil stations of Italy
- data from the italian Ministry of Economic Development
- count the total of the gas&oil stations for each municipality of Trentino
2 - identify the difference of municipalities in Trentino in the year 2019 with the year 2022
- identify which municipalities are created from aggregation to others
- find the biggest new municipality of Trentino and show all the italian municipalities with bordering it
- create the macroarea of all the municipalities bordering with it
- for each gas&oil station in the macro-area, calculate how many monumental trees have been within a 500m radius
3 - creates a polygon that contains all the monumental trees inside the area
- identify all the gas&oil stations in this area which are within 2km of each other
- save the polygon in geopackage with the attribute “description” with the name of the gas&oil station
4 - create the polygon of the Island of Elba from the layer of municipalities with functions of overlay
Setup
try:
import pygeos
except ModuleNotFoundError as e:
!pip install pygeos==0.13
import pygeos
# needed to use the method .explore in geopandas in colab
try:
import mapclassify
except ModuleNotFoundError as e:
!pip install mapclassify==2.4.3
import mapclassify
if mapclassify.__version__ != "2.4.3":
!pip install mapclassify==2.4.3
import mapclassify
try:
import geopandas as gpd
except ModuleNotFoundError as e:
!pip install geopandas==0.10.2
import geopandas as gpd
if gpd.__version__ != "0.10.2":
!pip install -U geopandas==0.10.2
import geopandas as gpd
Import of the packages
import geopandas as gpd
import requests
import matplotlib.pyplot as plt
import pandas as pd
pd.options.mode.chained_assignment = None
1 - create the geodataframe of the gas&oil stations of Italy
urlfile = "https://www.mise.gov.it/images/exportCSV/anagrafica_impianti_attivi.csv"
stations = None
try:
stations = pd.read_csv(urlfile,skiprows=1,sep=";",encoding="ISO-8859-1")
except Exception as e:
print(e)
Error tokenizing data. C error: Expected 10 fields in line 1498, saw 11
There several erros
stations = pd.read_csv(urlfile,skiprows=1,on_bad_lines='warn',sep=";",encoding="ISO-8859-1")
b'Skipping line 1498: expected 10 fields, saw 11\nSkipping line 8633: expected 10 fields, saw 11\nSkipping line 8634: expected 10 fields, saw 11\nSkipping line 8763: expected 10 fields, saw 11\nSkipping line 13173: expected 10 fields, saw 11\n'
this message show you the rows with problems to solve
stations = pd.read_csv(urlfile,skiprows=1,on_bad_lines='skip',sep=";",encoding="ISO-8859-1")
stations.head(5)
| idImpianto | Gestore | Bandiera | Tipo Impianto | Nome Impianto | Indirizzo | Comune | Provincia | Latitudine | Longitudine | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 53637 | ODELLI LUCIA | Pompe Bianche | Stradale | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 52829 | ERRA NICOLA | Api-Ip | Stradale | NaN | NaN | NaN | NaN | 40.716039 | 14.941328 |
| 2 | 52614 | FONZI PAOLO | Api-Ip | Stradale | stazione di servizio IP di Fonzi Paolo | via s francesco d'assisi snc - - | NaN | NaN | NaN | NaN |
| 3 | 53552 | KHALIL NOMAN | Agip Eni | Stradale | NaN | VIA LINCOLN 69 20092 | CINISELLO BALSAMO | NaN | NaN | NaN |
| 4 | 53906 | GE.CA.R SOCIETA' A RESPONSABILITA' LIMITATA SE... | Q8 | Stradale | NaN | VIA L.CAVALLARO SNC 84018 | SCAFATI | NaN | NaN | NaN |
stations.columns
Index(['idImpianto', 'Gestore', 'Bandiera', 'Tipo Impianto', 'Nome Impianto',
'Indirizzo', 'Comune', 'Provincia', 'Latitudine', 'Longitudine'],
dtype='object')
columns = {
'idImpianto': 'id',
'Gestore': 'manager',
'Bandiera':'company',
'Tipo Impianto':'type',
'Nome Impianto':'name',
'Indirizzo':'address',
'Comune':'city',
'Provincia':'province',
'Latitudine':'latitude',
'Longitudine':'longitude'
}
stations.rename(columns=columns,inplace=True)
stations.head(3)
| id | manager | company | type | name | address | city | province | latitude | longitude | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 53637 | ODELLI LUCIA | Pompe Bianche | Stradale | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 52829 | ERRA NICOLA | Api-Ip | Stradale | NaN | NaN | NaN | NaN | 40.716039 | 14.941328 |
| 2 | 52614 | FONZI PAOLO | Api-Ip | Stradale | stazione di servizio IP di Fonzi Paolo | via s francesco d'assisi snc - - | NaN | NaN | NaN | NaN |
create the geodataframe
geo_stations = gpd.GeoDataFrame(
stations,
crs='EPSG:4326',
geometry=gpd.points_from_xy(stations.longitude, stations.latitude))
geo_stations[geo_stations.geometry.is_empty].shape[0]
8
Error:
the values should be zero: the geodataframe should contains points.
Maybe there are some rows where the values of latitude and longitude aren’t present
stations.latitude.isnull().sum()
8
stations.longitude.isnull().sum()
8
8 … the same value for the invalid geometries
stations = stations[~stations.latitude.isnull()]
geo_stations = gpd.GeoDataFrame(
stations,
crs='EPSG:4326',
geometry=gpd.points_from_xy(stations.longitude, stations.latitude))
geo_stations[geo_stations.geometry.is_empty].shape[0]
0
Now it’s ZERO :)
geo_stations.plot()
plt.show()
Ehm … there is a point outside Italy
fix the problem
We need to identify the points outside the boundary of Italy
- create the polygon with the administrative boundaries of Italy
- check each point if is inside the polygon
macroregions22 = gpd.read_file("https://github.com/napo/geospatial_course_unitn/raw/master/data/istat/macroregions_2022_not_generalized.zip")
italy = macroregions22.to_crs(epsg=4326).dissolve()
italy.plot()
plt.show()
italy.geometry.values[0]
The fast way approach
USE RTree Index
using str-packed-r-tree
Shapely provides an interface to the query-only GEOS R-tree packed using the Sort-Tile-Recursive algorithm. Pass a list of geometry objects to the STRtree constructor to create a spatial index that you can query with another geometric object. Query-only means that once created, the STRtree is immutable. You cannot add or remove geometries.
R-trees are tree data structures used for spatial access methods, i.e., for indexing multi-dimensional information such as geographical coordinates, rectangles or polygons
R-Tree
Guttman, A. (1984). “R-Trees: A Dynamic Index Structure for Spatial Searching” (PDF). Proceedings of the 1984 ACM SIGMOD international conference on Management of data – SIGMOD ‘84. p. 47. doi:10.1145/602259.602266. ISBN 978-0897911283. S2CID 876601
Sort-Tile-Recursive
Leutenegger, Scott T.; Edgington, Jeffrey M.; Lopez, Mario A. (February 1997). “STR: A Simple and Efficient Algorithm for R-Tree Packing”.
from shapely.strtree import STRtree
def checkInsideItaly(p,i):
touch = len(STRtree([p]).query(i))
if touch == 0:
return False
else:
return True
geo_stations['in_italy'] = geo_stations.geometry.apply(lambda g: checkInsideItaly(g,italy.geometry.values[0]))
%time
CPU times: user 4 µs, sys: 0 ns, total: 4 µs
Wall time: 7.39 µs
wrong_positions = geo_stations[geo_stations['in_italy'] == False]
wrong_positions
| id | manager | company | type | name | address | city | province | latitude | longitude | geometry | in_italy | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 11410 | 52311 | PALMIGIANO CARBURANTI SAS DI PALMIGIANO MARIA ... | Agip Eni | Stradale | eni Palmigiano Carburanti s.a.s. | via nuova nola 80036 | PALMA CAMPANIA | NaN | 2.24347 | 16.409572 | POINT (16.40957 2.24347) | False |
The fastest ways
points = list(geo_stations.geometry.values)
strtree = STRtree(points)
points_cross = strtree.query(italy.geometry.values[0])
d = {'geometry': points}
gdf_points = gpd.GeoDataFrame(d, crs="EPSG:4326")
d = {'geometry': points_cross}
gdf_points_cross = gpd.GeoDataFrame(d, crs="EPSG:4326")
point_to_exclude = pd.concat([gdf_points,gdf_points_cross]).drop_duplicates(keep=False)
%time
CPU times: user 4 µs, sys: 0 ns, total: 4 µs
Wall time: 6.44 µs
wrong_position = geo_stations[geo_stations.geometry == point_to_exclude.geometry.values[0]]
wrong_position
| id | manager | company | type | name | address | city | province | latitude | longitude | geometry | in_italy | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 11410 | 52311 | PALMIGIANO CARBURANTI SAS DI PALMIGIANO MARIA ... | Agip Eni | Stradale | eni Palmigiano Carburanti s.a.s. | via nuova nola 80036 | PALMA CAMPANIA | NaN | 2.24347 | 16.409572 | POINT (16.40957 2.24347) | False |
.. the slow version …
uncomment this code, run … and drink a big cup of coffee
#geo_stations['in_italy'] = geo_stations.geometry.apply(lambda g: g.within(italy.geometry.values[0]))
#wrong_positions = geo_stations[geo_stations['in_italy'] == False]
#%time
#wrong_position
#stations_province_trento = geo_stations[geo_stations.within(italy.geometry.values[0])]
#%time
geo_stations = geo_stations[geo_stations.geometry != point_to_exclude.geometry.values[0]]
geo_stations.plot()
plt.show()
count the total of the gas&oil stations for each muncipality of Trentino
On the GitHub repository of the course there are the geopackage files with the administrative limits of ISTAT 2019 and 2022 with generalized geometries
We download the data of both due the second issue of the exercise.
url2019 = 'https://github.com/napo/geospatial_course_unitn/raw/master/data/istat/istat_administrative_units_2019.gpkg'
url2022 = 'https://github.com/napo/geospatial_course_unitn/raw/master/data/istat/istat_administrative_units_2022.gpkg'
istat2019 = "istat_administrative_units_2019.gpkg"
istat2022 = "istat_administrative_units_2022.gpkg"
r = requests.get(url2019, allow_redirects=True)
open(istat2019, 'wb').write(r.content)
117116928
r = requests.get(url2022, allow_redirects=True)
open(istat2022, 'wb').write(r.content)
116932608
import fiona
fiona.listlayers(istat2019)
['provinces', 'regions', 'macroregions', 'municipalities']
fiona.listlayers(istat2022)
['provinces', 'regions', 'municipalities', 'macroregions']
provinces2022 = gpd.read_file(istat2022,layer="provinces")
provinces2022.head(3)
| COD_RIP | COD_REG | COD_PROV | COD_CM | COD_UTS | DEN_PROV | DEN_CM | DEN_UTS | SIGLA | TIPO_UTS | SHAPE_AREA | Shape_Leng | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1 | 201 | 201 | - | Torino | Torino | TO | Citta metropolitana | 6.826908e+09 | 593389.667001 | MULTIPOLYGON (((411015.006 5049970.983, 411070... |
| 1 | 1 | 1 | 2 | 0 | 2 | Vercelli | - | Vercelli | VC | Provincia | 2.081602e+09 | 458754.449021 | MULTIPOLYGON (((437900.552 5088796.204, 437915... |
| 2 | 1 | 1 | 3 | 0 | 3 | Novara | - | Novara | NO | Provincia | 1.340250e+09 | 276722.284585 | MULTIPOLYGON (((459146.367 5079451.275, 459180... |
provinces2022.DEN_PROV.unique()
array(['-', 'Vercelli', 'Novara', 'Cuneo', 'Asti', 'Alessandria', 'Aosta',
'Imperia', 'Savona', 'La Spezia', 'Varese', 'Como', 'Sondrio',
'Bergamo', 'Brescia', 'Pavia', 'Cremona', 'Mantova', 'Bolzano',
'Trento', 'Verona', 'Vicenza', 'Belluno', 'Treviso', 'Padova',
'Rovigo', 'Udine', 'Gorizia', 'Trieste', 'Piacenza', 'Parma',
"Reggio nell'Emilia", 'Modena', 'Ferrara', 'Ravenna',
"Forli'-Cesena", 'Pesaro e Urbino', 'Ancona', 'Macerata',
'Ascoli Piceno', 'Massa Carrara', 'Lucca', 'Pistoia', 'Livorno',
'Pisa', 'Arezzo', 'Siena', 'Grosseto', 'Perugia', 'Terni',
'Viterbo', 'Rieti', 'Latina', 'Frosinone', 'Caserta', 'Benevento',
'Avellino', 'Salerno', "L'Aquila", 'Teramo', 'Pescara', 'Chieti',
'Campobasso', 'Foggia', 'Taranto', 'Brindisi', 'Lecce', 'Potenza',
'Matera', 'Cosenza', 'Catanzaro', 'Trapani', 'Agrigento',
'Caltanissetta', 'Enna', 'Ragusa', 'Siracusa', 'Sassari', 'Nuoro',
'Pordenone', 'Isernia', 'Oristano', 'Biella', 'Lecco', 'Lodi',
'Rimini', 'Prato', 'Crotone', 'Vibo Valentia',
'Verbano-Cusio-Ossola', 'Monza e della Brianza', 'Fermo',
'Barletta-Andria-Trani', 'Sud Sardegna'], dtype=object)
choose the province of Trento
province_of_trento = provinces2022[provinces2022['DEN_PROV']=='Trento']
province_of_trento.crs
<Derived Projected CRS: EPSG:32632>
Name: WGS 84 / UTM zone 32N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 6°E and 12°E, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Austria. Cameroon. Denmark. Equatorial Guinea. France. Gabon. Germany. Italy. Libya. Liechtenstein. Monaco. Netherlands. Niger. Nigeria. Norway. Sao Tome and Principe. Svalbard. Sweden. Switzerland. Tunisia. Vatican City State.
- bounds: (6.0, 0.0, 12.0, 84.0)
Coordinate Operation:
- name: UTM zone 32N
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
boundary_province_of_trento = province_of_trento.to_crs(epsg=4326).geometry.values[0]
plot it
boundary_province_of_trento
stations_province_trento = geo_stations[geo_stations.within(boundary_province_of_trento)]
stations_province_trento.plot()
plt.show()
stations_province_trento.shape[0]
211
without spatial relationship
stations.province.unique()
array([nan, 'AG', 'AL', 'AN', 'AO', 'AP', 'AQ', 'AR', 'AT', 'AV', 'BA',
'BG', 'BI', 'BL', 'BN', 'BO', 'BR', 'BS', 'BT', 'BZ', 'CA', 'CB',
'CE', 'CH', 'CI', 'CL', 'CN', 'CO', 'CR', 'CS', 'CT', 'CZ', 'EN',
'FC', 'FE', 'FG', 'FI', 'FM', 'FR', 'GE', 'GO', 'GR', 'IM', 'IS',
'KR', 'LC', 'LE', 'LI', 'LO', 'LT', 'LU', 'MB', 'MC', 'ME', 'MI',
'MN', 'MO', 'MS', 'MT', 'NO', 'NU', 'OG', 'OR', 'OT', 'PA', 'PC',
'PD', 'PE', 'PG', 'PI', 'PN', 'PO', 'PR', 'PT', 'PU', 'PV', 'PZ',
'RA', 'RC', 'RE', 'RG', 'RI', 'RM', 'RN', 'RO', 'SA', 'SI', 'SO',
'SP', 'SR', 'SS', 'SU', 'SV', 'TA', 'TE', 'TN', 'TO', 'TP', 'TR',
'TS', 'TV', 'UD', 'VA', 'VB', 'VC', 'VE', 'VI', 'VR', 'VS', 'VT',
'VV'], dtype=object)
provinces2022[provinces2022['DEN_PROV']=='Trento']['SIGLA'].unique()
array(['TN'], dtype=object)
stations[stations['province']=='TN'].shape[0]
211
#municipalities2022 = gpd.read_file(urlmunicipalities2022, encoding='iso-8859-15')
municipalities2022 = gpd.read_file(istat2022,layer="municipalities")
municipalities2022.head(3)
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1 | 201 | 201 | 1077 | 001077 | Chiaverano | None | 0 | 1.202212e+07 | 18164.236621 | MULTIPOLYGON (((414358.390 5042001.044, 414381... |
| 1 | 1 | 1 | 1 | 201 | 201 | 1079 | 001079 | Chiesanuova | None | 0 | 4.118911e+06 | 10777.318814 | MULTIPOLYGON (((394621.039 5031581.116, 394716... |
| 2 | 1 | 1 | 1 | 201 | 201 | 1089 | 001089 | Coazze | None | 0 | 5.657268e+07 | 41591.122092 | MULTIPOLYGON (((364914.897 4993224.894, 364929... |
cod_prov_trento = provinces2022[provinces2022.DEN_PROV == 'Trento'].COD_PROV.values[0]
municipalities_trentino_2022 = municipalities2022[municipalities2022.COD_PROV == cod_prov_trento]
municipalities_trentino_2022.head(3)
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 113 | 2 | 4 | 22 | 0 | 22 | 22172 | 022172 | Segonzano | None | 0 | 2.071212e+07 | 21551.885438 | MULTIPOLYGON (((677853.445 5121202.433, 677866... |
| 114 | 2 | 4 | 22 | 0 | 22 | 22177 | 022177 | Sover | None | 0 | 1.482203e+07 | 18613.838987 | MULTIPOLYGON (((680005.555 5123944.912, 680049... |
| 115 | 2 | 4 | 22 | 0 | 22 | 22195 | 022195 | Terzolas | None | 0 | 5.588756e+06 | 15881.333074 | MULTIPOLYGON (((648742.270 5136116.640, 648747... |
province_of_trento = italy = municipalities_trentino_2022.dissolve(by='COD_PROV')
%time
CPU times: user 4 µs, sys: 0 ns, total: 4 µs
Wall time: 7.15 µs
boundary_province_of_trento = province_of_trento.to_crs(epsg=4326).geometry.values[0]
boundary_province_of_trento
stations_province_trento = geo_stations[geo_stations.within(boundary_province_of_trento)]
stations_province_trento.shape[0]
211
the total is right ;)
and also the spatial relationship
now we can count the number of gas&oil stations for each municipality of Trentino
with the dataframe way
stations_by_municipalities = stations_province_trento.groupby(['city']).size().reset_index().rename(columns={0:'total'}).sort_values(['total','city'],ascending=[False,True])
%time
CPU times: user 7 µs, sys: 0 ns, total: 7 µs
Wall time: 14.3 µs
stations_by_municipalities = stations_province_trento.groupby(['city']).size().reset_index().rename(columns={0:'total'}).sort_values(['total','city'],ascending=[False,True])
stations_by_municipalities
| city | total | |
|---|---|---|
| 91 | TRENTO | 33 |
| 75 | ROVERETO | 13 |
| 65 | PERGINE VALSUGANA | 7 |
| 3 | ARCO | 6 |
| 47 | LAVIS | 6 |
| ... | ... | ... |
| 96 | VIGOLO VATTARO | 1 |
| 97 | VILLA AGNEDO | 1 |
| 98 | VILLA LAGARINA | 1 |
| 99 | VOLANO | 1 |
| 100 | ZUCLO | 1 |
101 rows × 2 columns
but … if the columns “city” is not present?
del stations_province_trento['city'] #delete the column city
stations_province_trento.head(3)
| id | manager | company | type | name | address | province | latitude | longitude | geometry | in_italy | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 18538 | 52336 | CR.VU. DI CRACIUN CRISTINA | Agip Eni | Stradale | CR.VU | CORSO PASSO BUOLE 1 38061 | TN | 45.763267 | 11.006058 | POINT (11.00606 45.76327) | True |
| 18539 | 5169 | ESSO | Esso | Stradale | ESSO | STRADA STAT. 12 KM 342 SNC 38061, ALA (TN) 38061 | TN | 45.803566 | 11.019186 | POINT (11.01919 45.80357) | True |
| 18540 | 7317 | FERRARI ATTILIO | Pompe Bianche | Stradale | FERRARI ATTILIO | CORSO VERONA 22 38061 | TN | 45.757343 | 10.999531 | POINT (10.99953 45.75734) | True |
reconstruct the name of the city associated for each location
def getNameCity(point,cities):
name = cities[cities.to_crs(epsg=4326).contains(point)].COMUNE.values[0]
return name
stations_province_trento['city'] = stations_province_trento.geometry.apply(lambda point: getNameCity(point,municipalities_trentino_2022))
%time
CPU times: user 3 µs, sys: 0 ns, total: 3 µs
Wall time: 6.44 µs
stations_province_trento.head(3)
| id | manager | company | type | name | address | province | latitude | longitude | geometry | in_italy | city | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 18538 | 52336 | CR.VU. DI CRACIUN CRISTINA | Agip Eni | Stradale | CR.VU | CORSO PASSO BUOLE 1 38061 | TN | 45.763267 | 11.006058 | POINT (11.00606 45.76327) | True | Ala |
| 18539 | 5169 | ESSO | Esso | Stradale | ESSO | STRADA STAT. 12 KM 342 SNC 38061, ALA (TN) 38061 | TN | 45.803566 | 11.019186 | POINT (11.01919 45.80357) | True | Ala |
| 18540 | 7317 | FERRARI ATTILIO | Pompe Bianche | Stradale | FERRARI ATTILIO | CORSO VERONA 22 38061 | TN | 45.757343 | 10.999531 | POINT (10.99953 45.75734) | True | Ala |
stations_by_municipalities = stations_province_trento.groupby(['city']).size().reset_index().rename(columns={0:'total'}).sort_values(['total','city'],ascending=[False,True])
stations_by_municipalities
| city | total | |
|---|---|---|
| 87 | Trento | 32 |
| 71 | Rovereto | 13 |
| 43 | Lavis | 7 |
| 61 | Pergine Valsugana | 7 |
| 4 | Arco | 6 |
| ... | ... | ... |
| 88 | Valdaone | 1 |
| 89 | Vallarsa | 1 |
| 91 | Vermiglio | 1 |
| 92 | Villa Lagarina | 1 |
| 94 | Volano | 1 |
95 rows × 2 columns
2 - identify the difference of municipalities in Trentino of year 2019 with year 2022
municipalities2019 = gpd.read_file(istat2019,layer="municipalities")
municipalities_trentino_2019 = municipalities2019[municipalities2019['COD_PROV'] == cod_prov_trento]
names2019 = list(municipalities_trentino_2019.COMUNE.unique())
names2022 = list(municipalities_trentino_2022.COMUNE.unique())
notpresentin2022 = list(set(names2019) - set(names2022))
list of municipalities that were present in 2019
notpresentin2022
['Malosco',
'Romallo',
'Carano',
'Varena',
'Cagnò',
'Daiano',
'Fondo',
'Cloz',
'Revò',
'Faedo',
'Brez',
'Castelfondo']
notpresentin2019 = list(set(names2022) - set(names2019))
list of the new municipalities in 2022
notpresentin2019
["Borgo d'Anaunia", 'Ville di Fiemme', 'Novella']
old_municipalities_2019 = municipalities_trentino_2019[municipalities_trentino_2019.COMUNE.isin(notpresentin2022)]
old_municipalities_2019.plot()
plt.show()
old_municipalities_2019
| COD_RIP | COD_REG | COD_PROV | COD_CM | COD_UTS | PRO_COM | PRO_COM_T | COMUNE | COMUNE_A | CC_UTS | SHAPE_LENG | SHAPE_AREA | SHAPE_LEN | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 724 | 2 | 4 | 22 | 0 | 22 | 22027 | 022027 | Brez | None | 0 | 26280.588774 | 1.916712e+07 | 26280.400804 | MULTIPOLYGON (((658244.658 5151579.353, 658250... |
| 2168 | 2 | 4 | 22 | 0 | 22 | 22041 | 022041 | Carano | None | 0 | 22704.690339 | 1.357006e+07 | 22704.527934 | MULTIPOLYGON (((685476.521 5132926.845, 685524... |
| 3247 | 2 | 4 | 22 | 0 | 22 | 22030 | 022030 | Cagnò | None | 0 | 16888.878461 | 3.407854e+06 | 16888.757603 | MULTIPOLYGON (((657219.683 5144541.638, 657205... |
| 3343 | 2 | 4 | 22 | 0 | 22 | 22111 | 022111 | Malosco | None | 0 | 23242.468818 | 6.728356e+06 | 23242.303282 | MULTIPOLYGON (((667424.906 5146341.256, 667423... |
| 3346 | 2 | 4 | 22 | 0 | 22 | 22070 | 022070 | Daiano | None | 0 | 14916.892170 | 9.510673e+06 | 14916.785768 | MULTIPOLYGON (((688996.496 5133085.842, 688983... |
| 4269 | 2 | 4 | 22 | 0 | 22 | 22063 | 022063 | Cloz | None | 0 | 13642.729963 | 8.211168e+06 | 13642.632958 | MULTIPOLYGON (((659208.764 5145192.728, 659220... |
| 4869 | 2 | 4 | 22 | 0 | 22 | 22154 | 022154 | Romallo | None | 0 | 7687.685830 | 2.444897e+06 | 7687.631108 | MULTIPOLYGON (((659180.319 5141366.583, 659191... |
| 5142 | 2 | 4 | 22 | 0 | 22 | 22046 | 022046 | Castelfondo | None | 0 | 31016.737999 | 2.587157e+07 | 31016.516642 | MULTIPOLYGON (((660043.059 5153438.578, 660045... |
| 5540 | 2 | 4 | 22 | 0 | 22 | 22152 | 022152 | Revò | None | 0 | 24487.292355 | 1.335782e+07 | 24487.116841 | MULTIPOLYGON (((657885.646 5145988.024, 657936... |
| 5607 | 2 | 4 | 22 | 0 | 22 | 22211 | 022211 | Varena | None | 0 | 25565.413219 | 2.306859e+07 | 25565.229395 | MULTIPOLYGON (((691744.479 5137215.811, 691809... |
| 5622 | 2 | 4 | 22 | 0 | 22 | 22088 | 022088 | Fondo | None | 0 | 32412.749470 | 3.062881e+07 | 32412.517789 | MULTIPOLYGON (((668115.125 5152346.817, 668120... |
| 5652 | 2 | 4 | 22 | 0 | 22 | 22080 | 022080 | Faedo | None | 0 | 16440.165284 | 1.068038e+07 | 16440.047652 | MULTIPOLYGON (((667690.769 5121538.436, 667726... |
new_municipalities_2022 = municipalities_trentino_2022[municipalities_trentino_2022.COMUNE.isin(notpresentin2019)]
new_municipalities_2022.plot()
plt.show()
new_municipalities_2022
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7887 | 2 | 4 | 22 | 0 | 22 | 22252 | 022252 | Borgo d'Anaunia | None | 0 | 6.322874e+07 | 64924.029558 | MULTIPOLYGON (((660043.059 5153438.578, 660045... |
| 7889 | 2 | 4 | 22 | 0 | 22 | 22253 | 022253 | Novella | None | 0 | 4.658886e+07 | 44596.283777 | MULTIPOLYGON (((658244.658 5151579.353, 658250... |
| 7902 | 2 | 4 | 22 | 0 | 22 | 22254 | 022254 | Ville di Fiemme | None | 0 | 4.614933e+07 | 42004.727128 | MULTIPOLYGON (((692073.476 5137147.812, 692107... |
identify which municipalities are created from aggregation to others
def whereincluded(geometry, geometries_gdf):
name = "not included"
found = geometries_gdf[geometries_gdf.geometry.contains(geometry)]
if len(found) > 0:
name = found.COMUNE.values[0]
return(name)
old_municipalities_2019['included_in'] = old_municipalities_2019.geometry.apply(lambda g: whereincluded(g,new_municipalities_2022))
old_municipalities_2019[['COMUNE','included_in']]
| COMUNE | included_in | |
|---|---|---|
| 724 | Brez | Novella |
| 2168 | Carano | Ville di Fiemme |
| 3247 | Cagnò | Novella |
| 3343 | Malosco | Borgo d'Anaunia |
| 3346 | Daiano | Ville di Fiemme |
| 4269 | Cloz | Novella |
| 4869 | Romallo | Novella |
| 5142 | Castelfondo | Borgo d'Anaunia |
| 5540 | Revò | Novella |
| 5607 | Varena | Ville di Fiemme |
| 5622 | Fondo | Borgo d'Anaunia |
| 5652 | Faedo | not included |
Where is “Faedo” ?
faedo = old_municipalities_2019[old_municipalities_2019.COMUNE == 'Faedo']
faedo
| COD_RIP | COD_REG | COD_PROV | COD_CM | COD_UTS | PRO_COM | PRO_COM_T | COMUNE | COMUNE_A | CC_UTS | SHAPE_LENG | SHAPE_AREA | SHAPE_LEN | geometry | included_in | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 5652 | 2 | 4 | 22 | 0 | 22 | 22080 | 022080 | Faedo | None | 0 | 16440.165284 | 1.068038e+07 | 16440.047652 | MULTIPOLYGON (((667690.769 5121538.436, 667726... | not included |
faedo_geometry = faedo.geometry.values[0]
municipalities_trentino_2022[municipalities_trentino_2022.geometry.within(faedo_geometry)]
| 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 |
|---|
municipalities_trentino_2022[municipalities_trentino_2022.geometry.intersects(faedo_geometry)]
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 4002 | 2 | 4 | 22 | 0 | 22 | 22092 | 022092 | Giovo | None | 0 | 2.081293e+07 | 33094.134912 | MULTIPOLYGON (((668394.632 5121453.303, 668403... |
| 6310 | 2 | 4 | 22 | 0 | 22 | 22116 | 022116 | Mezzocorona | None | 0 | 2.535304e+07 | 29851.314722 | MULTIPOLYGON (((665540.410 5124920.912, 665610... |
| 7900 | 2 | 4 | 22 | 0 | 22 | 22167 | 022167 | San Michele all'Adige | None | 0 | 1.599623e+07 | 26368.864693 | MULTIPOLYGON (((667690.769 5121538.436, 667726... |
municipalities_trentino_2022[municipalities_trentino_2022.geometry.intersects(faedo_geometry.representative_point())]
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7900 | 2 | 4 | 22 | 0 | 22 | 22167 | 022167 | San Michele all'Adige | None | 0 | 1.599623e+07 | 26368.864693 | MULTIPOLYGON (((667690.769 5121538.436, 667726... |
faedo_is_in = municipalities_trentino_2022[municipalities_trentino_2022.geometry.intersects(faedo_geometry.representative_point())]
faedo_is_in
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7900 | 2 | 4 | 22 | 0 | 22 | 22167 | 022167 | San Michele all'Adige | None | 0 | 1.599623e+07 | 26368.864693 | MULTIPOLYGON (((667690.769 5121538.436, 667726... |
faedo_new_municipality = faedo_is_in.COMUNE.values[0]
faedo_new_municipality
"San Michele all'Adige"
list_changed_municipalities = old_municipalities_2019[old_municipalities_2019.included_in != 'not included']
list_changed_municipalities = list(list_changed_municipalities.included_in.unique())
list_changed_municipalities.append(faedo_new_municipality)
list_changed_municipalities
['Novella', 'Ville di Fiemme', "Borgo d'Anaunia", "San Michele all'Adige"]
and we can do the same with the polygons
new_municipalities_trentino_2022 = municipalities_trentino_2022[municipalities_trentino_2022.COMUNE.isin(list_changed_municipalities)]
new_municipalities_trentino_2022.plot()
plt.show()
find the biggest new municipality of Trentino and show all the italian municipalities with bordering it
biggest_new_municipality_trentino = new_municipalities_trentino_2022[new_municipalities_trentino_2022.geometry.area == new_municipalities_trentino_2022.geometry.area.max()]
biggest_new_municipality_trentino.plot()
plt.show()
biggest_new_municipality_trentino
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7887 | 2 | 4 | 22 | 0 | 22 | 22252 | 022252 | Borgo d'Anaunia | None | 0 | 6.322874e+07 | 64924.029558 | MULTIPOLYGON (((660043.059 5153438.578, 660045... |
boundary_borgo_anaunia = biggest_new_municipality_trentino.geometry.values[0]
around_borgo_anaunia = municipalities2022[municipalities2022.touches(boundary_borgo_anaunia)]
around_borgo_anaunia
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2617 | 2 | 4 | 21 | 0 | 21 | 21043 | 021043 | Lauregno | Laurein | 0 | 1.388482e+07 | 23321.072254 | MULTIPOLYGON (((655860.135 5153156.388, 655884... |
| 3437 | 2 | 4 | 21 | 0 | 21 | 21084 | 021084 | San Pancrazio | St. Pankraz | 0 | 6.317246e+07 | 42883.848255 | MULTIPOLYGON (((659549.008 5165260.582, 659593... |
| 3615 | 2 | 4 | 21 | 0 | 21 | 21118 | 021118 | Senale-San Felice | Unsere Liebe Frau im Walde-St. Felix | 0 | 2.763132e+07 | 25151.845343 | MULTIPOLYGON (((660357.881 5155829.884, 660405... |
| 4603 | 2 | 4 | 21 | 0 | 21 | 21004 | 021004 | Appiano sulla strada del vino | Eppan an der Weinstraße | 0 | 5.945015e+07 | 49491.808002 | MULTIPOLYGON (((669345.871 5154254.383, 669361... |
| 5181 | 2 | 4 | 22 | 0 | 22 | 22170 | 022170 | Sarnonico | None | 0 | 1.218983e+07 | 46246.416637 | MULTIPOLYGON (((668419.900 5147274.750, 668412... |
| 5283 | 2 | 4 | 22 | 0 | 22 | 22159 | 022159 | Ronzone | None | 0 | 5.297326e+06 | 17188.159009 | MULTIPOLYGON (((670031.369 5149451.313, 670162... |
| 7889 | 2 | 4 | 22 | 0 | 22 | 22253 | 022253 | Novella | None | 0 | 4.658886e+07 | 44596.283777 | MULTIPOLYGON (((658244.658 5151579.353, 658250... |
around_borgo_anaunia.plot()
plt.show()
create the macroarea of all the municipalities bordering with it
new_area = around_borgo_anaunia.append(biggest_new_municipality_trentino).dissolve()
/tmp/ipykernel_32044/1213046248.py:1: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
new_area = around_borgo_anaunia.append(biggest_new_municipality_trentino).dissolve()
new_area = new_area[['geometry']]
new_area['name'] = "area of borgo d'anaunia and bordering municipalities"
new_area.plot()
plt.show()
for each gas&oil station in the macro-area, calculate how many monumental trees have been within a 500m radius
the dataset in GeoJSON of the italian monumental trees is created with the code of the lesson 02
You can find the dataset here
macroarea_geometry = new_area.to_crs(epsg=4326).geometry.values[0]
stations_in_macroarea = geo_stations[geo_stations.within(macroarea_geometry)]
monumental_trees = gpd.read_file('https://github.com/napo/geospatial_course_unitn/raw/master/data/mipaaf/geo_monumental_trees_20020726.geojson')
monumental_trees_in_macroarea = monumental_trees[monumental_trees.within(macroarea_geometry)]
def fivehundredfrom(point,points):
present = False
found = stations_in_macroarea[stations_in_macroarea.within(point)]
if len(found) > 0:
present = True
return(present)
monumental_trees_in_macroarea.to_crs(epsg=3263).geometry.buffer(500).apply(lambda point: fivehundredfrom(point,stations_in_macroarea.to_crs(epsg=32632)))
297 False
326 False
3833 False
dtype: bool
it’s normal that a gas&oil station is far away from a monumental tree :)
3 - creates a polygon that contains all the monumental trees inside the area
convex hull
solution: create a convex hull
In geometry, the convex hull or convex envelope or convex closure of a shape is the smallest convex set that contains it. The convex hull may be defined either as the intersection of all convex sets containing a given subset of a Euclidean space, or equivalently as the set of all convex combinations of points in the subset. For a bounded subset of the plane, the convex hull may be visualized as the shape enclosed by a rubber band stretched around the subset. (source: wikipedia)
monumental_trees_in_macroarea.plot()
plt.show()
area_of_monumental_trees_in_macroarea = monumental_trees_in_macroarea.unary_union.convex_hull
area_of_monumental_trees_in_macroarea
Concave Hull
Contrary to a convex hull, a concave hull can describe the shape of a point cloud.
Convex hull
Concave hull
Alpha shapes
Alpha shapes are often used to generalize bounding polygons containing sets of points. The alpha parameter is defined as the value a, such that an edge of a disk of radius 1/a can be drawn between any two edge members of a set of points and still contain all the points. The convex hull, a shape resembling what you would see if you wrapped a rubber band around pegs at all the data points, is an alpha shape where the alpha parameter is equal to zero
https://alphashape.readthedocs.io/
try:
import alphashape
except ModuleNotFoundError as e:
!pip install alphashape==1.3.1
import alphashape
if alphashape.__version__ != "1.3.1":
!pip install -U alphashape==1.3.1
import alphashape
alpha_shape = alphashape.alphashape(monumental_trees_in_macroarea, 100)
alpha_shape.plot()
plt.show()
… we have only three points … but if you want try with more …
convex_hull_trento = gpd.GeoDataFrame(
geometry=[stations_province_trento.geometry.unary_union.convex_hull],
columns=['geometry'],
crs=stations_province_trento.crs)
convex_hull_trento.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook
stations_province_trento.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook
alpha_paramenter = 60
alphashape.alphashape(stations_province_trento, alpha_paramenter).explore()
Make this Notebook Trusted to load map: File -> Trust Notebook
Creating alpha shapes around sets of points usually requires a visually interactive step where the alpha parameter for a concave hull is determined by iterating over or bisecting values to approach a best fit.
More informations: https://alphashape.readthedocs.io/en/latest/readme.html#using-a-varying-alpha-parameter
identify all the gas&oil stations in this area which are within 2km of each other
stations_in_area_monumental_trees = stations_in_macroarea[stations_in_macroarea.within(area_of_monumental_trees_in_macroarea)]
len(stations_in_area_monumental_trees)
0
stations_out_area_monumental_trees = stations_in_macroarea[~stations_in_macroarea.within(area_of_monumental_trees_in_macroarea)]
len(stations_out_area_monumental_trees)
8
nearest points
from shapely.ops import nearest_points
shapely offers a method to identify the nearest points between two geometries
Documentation here
def get_nearest_id(id, points):
# Create a union of points (multipoint geometry
multipoints = points[points.id != id]["geometry"].unary_union
# identify the starting point
point = points[points.id == id]
# find the nearest points
nearest_geoms = nearest_points(point['geometry'].values[0], multipoints)
# get corresponding values of the nearest point
# note: in the position 0 there is the starting point
nearest_data = points[points["geometry"] == nearest_geoms[1]]
# extract the id of the nearest point
nearest_id = nearest_data['id'].values[0]
return (nearest_id)
stations_in_macroarea['id_nearest'] = stations_in_macroarea['id'].apply(lambda x :get_nearest_id(x,stations_in_macroarea))
stations_in_macroarea
| id | manager | company | type | name | address | city | province | latitude | longitude | geometry | in_italy | id_nearest | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3491 | 7866 | ENI | Agip Eni | Stradale | ENI | STRADA DI CIRCONVALLAZIONE 4 39057 | APPIANO SULLA STRADA DEL VINO | BZ | 46.456748 | 11.268874 | POINT (11.26887 46.45675) | True | 27263 |
| 3492 | 49587 | MEBORAST | Agip Eni | Stradale | MEBORAST | 238 delle Palade, Km. 221, SUD - 39057 | APPIANO SULLA STRADA DEL VINO | BZ | 46.494556 | 11.281309 | POINT (11.28131 46.49456) | True | 7866 |
| 3493 | 9080 | Q8 des Karl Pichler | Q8 | Stradale | Q8 des Karl Pichler | VIA CALDARO 8 39057 | APPIANO SULLA STRADA DEL VINO | BZ | 46.444788 | 11.260646 | POINT (11.26065 46.44479) | True | 27263 |
| 3494 | 27263 | TSCHIGG HELMUT | Esso | Stradale | TSCHIGG HELMUT | VIA BOLZANO 5 39057 | APPIANO SULLA STRADA DEL VINO | BZ | 46.458640 | 11.261110 | POINT (11.26111 46.45864) | True | 7866 |
| 18558 | 23679 | Flaim Carlo | Tamoil | Stradale | Flaim Carlo | LIBERTA' 4/B 38028 | BREZ | TN | 46.431569 | 11.107855 | POINT (11.10786 46.43157) | True | 23500 |
| 18600 | 23500 | ESSO RAINER DI ZUCOL PIETRO | Esso | Stradale | ESSO RAINER DI ZUCOL PIETRO | VIA PALADE 49 38013, FONDO (TN) 49 38013 | FONDO | TN | 46.436768 | 11.139896 | POINT (11.13990 46.43677) | True | 12750 |
| 18663 | 50275 | DISTRIBUTORE IP REVO' | Api-Ip | Stradale | DISTRIBUTORE IP REVO' | 42 del Tonale e della Mendola, Km. 192 + 370, ... | REVO' | TN | 46.393205 | 11.063137 | POINT (11.06314 46.39321) | True | 23679 |
| 18688 | 12750 | BONANI GIULIANO | Agip Eni | Stradale | bonani giuliano | C. BATTISTI 1 38011 | SARNONICO | TN | 46.427077 | 11.141900 | POINT (11.14190 46.42708) | True | 23500 |
def getdistance(id,points):
points = points.to_crs(epsg=32632)
point = points[points.id == id]
id_nearest = point.id_nearest.values[0]
point_nearest = points[points.id == id_nearest]
from_geometry = point.geometry.values[0]
to_geometry = point_nearest.geometry.values[0]
dist = from_geometry.distance(to_geometry)
return (dist)
stations_in_macroarea['distance_to_nearest'] = stations_in_macroarea['id'].apply(lambda x :getdistance(x,stations_in_macroarea))
stations_in_macroarea
| id | manager | company | type | name | address | city | province | latitude | longitude | geometry | in_italy | id_nearest | distance_to_nearest | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3491 | 7866 | ENI | Agip Eni | Stradale | ENI | STRADA DI CIRCONVALLAZIONE 4 39057 | APPIANO SULLA STRADA DEL VINO | BZ | 46.456748 | 11.268874 | POINT (11.26887 46.45675) | True | 27263 | 632.397005 |
| 3492 | 49587 | MEBORAST | Agip Eni | Stradale | MEBORAST | 238 delle Palade, Km. 221, SUD - 39057 | APPIANO SULLA STRADA DEL VINO | BZ | 46.494556 | 11.281309 | POINT (11.28131 46.49456) | True | 7866 | 4309.826911 |
| 3493 | 9080 | Q8 des Karl Pichler | Q8 | Stradale | Q8 des Karl Pichler | VIA CALDARO 8 39057 | APPIANO SULLA STRADA DEL VINO | BZ | 46.444788 | 11.260646 | POINT (11.26065 46.44479) | True | 27263 | 1540.195343 |
| 3494 | 27263 | TSCHIGG HELMUT | Esso | Stradale | TSCHIGG HELMUT | VIA BOLZANO 5 39057 | APPIANO SULLA STRADA DEL VINO | BZ | 46.458640 | 11.261110 | POINT (11.26111 46.45864) | True | 7866 | 632.397005 |
| 18558 | 23679 | Flaim Carlo | Tamoil | Stradale | Flaim Carlo | LIBERTA' 4/B 38028 | BREZ | TN | 46.431569 | 11.107855 | POINT (11.10786 46.43157) | True | 23500 | 2529.179916 |
| 18600 | 23500 | ESSO RAINER DI ZUCOL PIETRO | Esso | Stradale | ESSO RAINER DI ZUCOL PIETRO | VIA PALADE 49 38013, FONDO (TN) 49 38013 | FONDO | TN | 46.436768 | 11.139896 | POINT (11.13990 46.43677) | True | 12750 | 1088.202974 |
| 18663 | 50275 | DISTRIBUTORE IP REVO' | Api-Ip | Stradale | DISTRIBUTORE IP REVO' | 42 del Tonale e della Mendola, Km. 192 + 370, ... | REVO' | TN | 46.393205 | 11.063137 | POINT (11.06314 46.39321) | True | 23679 | 5477.457684 |
| 18688 | 12750 | BONANI GIULIANO | Agip Eni | Stradale | bonani giuliano | C. BATTISTI 1 38011 | SARNONICO | TN | 46.427077 | 11.141900 | POINT (11.14190 46.42708) | True | 23500 | 1088.202974 |
4 - create the polygon of the Island of Elba from the layer of municipalities with functions of overlay
- identify Elba Island
- generate an area around the island
- clip the geometry
open umap and draw an area around Elba Island and download it in geojson format
the geojson file is available also here https://raw.githubusercontent.com/napo/geospatial_course_unitn/master/data/self_made/elbaisland.geojson
url_area = "https://raw.githubusercontent.com/napo/geospatial_course_unitn/master/data/self_made/elbaisland.geojson"
area_around_elba = gpd.read_file(url_area)
area_around_elba.plot()
plt.show()
gdf_clipped_area = gpd.clip(macroregions22.to_crs(epsg=4326), area_around_elba, keep_geom_type=False)
gdf_clipped_area.plot()
plt.show()
elba_island = gdf_clipped_area.dissolve()
elba_island
| geometry | COD_RIP | DEN_RIP | Shape_Area | Shape_Leng | |
|---|---|---|---|---|---|
| 0 | MULTIPOLYGON (((10.36101 42.71325, 10.36105 42... | 3 | Centro | 5.802768e+10 | 2.405546e+06 |
remove not usefull columns
del elba_island['COD_RIP']
del elba_island['DEN_RIP']
del elba_island['Shape_Area']
del elba_island['Shape_Leng']
elba_island['name'] = "Elba"
elba_island
| geometry | name | |
|---|---|---|
| 0 | MULTIPOLYGON (((10.36101 42.71325, 10.36105 42... | Elba |
if you want to know the names of the municipalities on the Island of Elba
municipalities_in_elba = municipalities2022[municipalities2022.to_crs(epsg=4326).geometry.intersects(elba_island.geometry.values[0]) == True]
%time
CPU times: user 6 µs, sys: 0 ns, total: 6 µs
Wall time: 10 µs
municipalities_in_elba
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 691 | 3 | 9 | 49 | 0 | 49 | 49014 | 049014 | Portoferraio | None | 0 | 4.802935e+07 | 68845.371050 | MULTIPOLYGON (((603254.255 4742679.377, 603266... |
| 1426 | 3 | 9 | 49 | 0 | 49 | 49021 | 049021 | Rio | None | 0 | 3.651822e+07 | 43276.023829 | MULTIPOLYGON (((615535.906 4747526.833, 615534... |
| 2349 | 3 | 9 | 49 | 0 | 49 | 49004 | 049004 | Capoliveri | None | 0 | 3.956061e+07 | 62573.909927 | MULTIPOLYGON (((608048.964 4737102.921, 608560... |
| 3177 | 3 | 9 | 49 | 0 | 49 | 49011 | 049011 | Marciana Marina | None | 0 | 5.858972e+06 | 17093.632943 | MULTIPOLYGON (((596961.056 4740648.895, 596970... |
| 6250 | 3 | 9 | 49 | 0 | 49 | 49003 | 049003 | Campo nell'Elba | None | 0 | 5.578762e+07 | 67951.434125 | MULTIPOLYGON (((600235.028 4736821.926, 600304... |
| 6251 | 3 | 9 | 49 | 0 | 49 | 49010 | 049010 | Marciana | None | 0 | 4.545090e+07 | 46878.864895 | MULTIPOLYGON (((593200.087 4740268.899, 593205... |
| 6252 | 3 | 9 | 49 | 0 | 49 | 49013 | 049013 | Porto Azzurro | None | 0 | 1.333011e+07 | 17720.940797 | MULTIPOLYGON (((613705.918 4738969.904, 614452... |
or
municipalities_in_elba = gpd.clip(municipalities2022.to_crs(epsg=4326), area_around_elba, keep_geom_type=False)
%time
CPU times: user 4 µs, sys: 0 ns, total: 4 µs
Wall time: 7.63 µs
municipalities_in_elba
| 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 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 691 | 3 | 9 | 49 | 0 | 49 | 49014 | 049014 | Portoferraio | None | 0 | 4.802935e+07 | 68845.371050 | MULTIPOLYGON (((10.26345 42.82968, 10.26373 42... |
| 6250 | 3 | 9 | 49 | 0 | 49 | 49003 | 049003 | Campo nell'Elba | None | 0 | 5.578762e+07 | 67951.434125 | POLYGON ((10.22616 42.77691, 10.22701 42.77654... |
| 6251 | 3 | 9 | 49 | 0 | 49 | 49010 | 049010 | Marciana | None | 0 | 4.545090e+07 | 46878.864895 | POLYGON ((10.13998 42.80928, 10.14025 42.80929... |
| 3177 | 3 | 9 | 49 | 0 | 49 | 49011 | 049011 | Marciana Marina | None | 0 | 5.858972e+06 | 17093.632943 | POLYGON ((10.18608 42.81221, 10.18619 42.81221... |
| 2349 | 3 | 9 | 49 | 0 | 49 | 49004 | 049004 | Capoliveri | None | 0 | 3.956061e+07 | 62573.909927 | MULTIPOLYGON (((10.32705 42.77487, 10.33141 42... |
| 6252 | 3 | 9 | 49 | 0 | 49 | 49013 | 049013 | Porto Azzurro | None | 0 | 1.333011e+07 | 17720.940797 | POLYGON ((10.39946 42.79169, 10.39999 42.79184... |
| 1426 | 3 | 9 | 49 | 0 | 49 | 49021 | 049021 | Rio | None | 0 | 3.651822e+07 | 43276.023829 | MULTIPOLYGON (((10.41450 42.87120, 10.41447 42... |
PS:
If someone tells you that there are 7 municipalities in Elba, tell them they are wrong: The Municipality of Rio was established in 2018 following the merger of the municipalities of Rio d’Elba, Rio Marina and others.
(https://en.wikipedia.org/wiki/Rio,_Italy)[https://en.wikipedia.org/wiki/Rio,_Italy]
