From “Geometric operations”#
Last updated: 2023-05-27 16:39:06
Exercise 09-b#
import geopandas as gpd
towns = gpd.read_file('data/muni_il.shp', encoding='utf-8')
towns = towns[['Muni_Heb', 'geometry']]
towns['area'] = towns.area
towns1 = towns.sort_values('area')
towns1['Muni_Heb'].iloc[0]
'ללא שיפוט - אזור נחל ארבל'
towns1['Muni_Heb'].iloc[towns1.shape[0]-1]
'רמת נגב'
Exercise 09-c#
import geopandas as gpd
towns = gpd.read_file('data/muni_il.shp', encoding='utf-8')
towns = towns[['Muni_Heb', 'Shape_Area', 'geometry']]
towns['area'] = towns.area
towns.plot.scatter('area', 'Shape_Area');
Exercise 09-g#
import geopandas as gpd
# Read data & subset columns
towns = gpd.read_file('data/muni_il.shp', encoding='utf-8')
towns = towns[['Machoz', 'geometry']]
towns
| Machoz | geometry | |
|---|---|---|
| 0 | צפון | POLYGON Z ((225369.655 770523.6... |
| 1 | דרום | POLYGON Z ((206899.135 552967.8... |
| 2 | מרכז | POLYGON Z ((194329.250 655299.1... |
| 3 | מרכז | POLYGON Z ((196236.573 657835.0... |
| 4 | צפון | POLYGON Z ((255150.135 764764.6... |
| ... | ... | ... |
| 405 | דרום | POLYGON Z ((164806.146 577898.8... |
| 406 | מרכז | POLYGON Z ((189803.359 684152.9... |
| 407 | צפון | POLYGON Z ((212294.953 763168.8... |
| 408 | ירושלים | POLYGON Z ((209066.649 635655.2... |
| 409 | דרום | POLYGON Z ((162082.027 592043.1... |
410 rows × 2 columns
# Dissolve
districts = towns.dissolve(by='Machoz').reset_index()
districts
| Machoz | geometry | |
|---|---|---|
| 0 | דרום | POLYGON Z ((214342.276 467536.8... |
| 1 | חיפה | MULTIPOLYGON Z (((205862.732 70... |
| 2 | ירושלים | POLYGON Z ((192314.625 638338.4... |
| 3 | מרכז | POLYGON Z ((171159.359 631505.0... |
| 4 | צפון | MULTIPOLYGON Z (((204805.266 71... |
| 5 | תל אביב | POLYGON Z ((177311.878 655799.7... |
# Add 'count' variable
districts['count'] = 1
districts
| Machoz | geometry | count | |
|---|---|---|---|
| 0 | דרום | POLYGON Z ((214342.276 467536.8... | 1 |
| 1 | חיפה | MULTIPOLYGON Z (((205862.732 70... | 1 |
| 2 | ירושלים | POLYGON Z ((192314.625 638338.4... | 1 |
| 3 | מרכז | POLYGON Z ((171159.359 631505.0... | 1 |
| 4 | צפון | MULTIPOLYGON Z (((204805.266 71... | 1 |
| 5 | תל אביב | POLYGON Z ((177311.878 655799.7... | 1 |
# Spatial join with self
districts = gpd.sjoin(
districts[['Machoz', 'geometry']],
districts[['count', 'geometry']],
how='left',
predicate='intersects'
)
districts = districts.drop('index_right', axis=1)
districts
| Machoz | geometry | count | |
|---|---|---|---|
| 0 | דרום | POLYGON Z ((214342.276 467536.8... | 1 |
| 0 | דרום | POLYGON Z ((214342.276 467536.8... | 1 |
| 0 | דרום | POLYGON Z ((214342.276 467536.8... | 1 |
| 1 | חיפה | MULTIPOLYGON Z (((205862.732 70... | 1 |
| 1 | חיפה | MULTIPOLYGON Z (((205862.732 70... | 1 |
| ... | ... | ... | ... |
| 3 | מרכז | POLYGON Z ((171159.359 631505.0... | 1 |
| 4 | צפון | MULTIPOLYGON Z (((204805.266 71... | 1 |
| 4 | צפון | MULTIPOLYGON Z (((204805.266 71... | 1 |
| 5 | תל אביב | POLYGON Z ((177311.878 655799.7... | 1 |
| 5 | תל אביב | POLYGON Z ((177311.878 655799.7... | 1 |
18 rows × 3 columns
# Aggregate
districts = districts.dissolve(by='Machoz', aggfunc='sum')
districts
| geometry | count | |
|---|---|---|
| Machoz | ||
| דרום | POLYGON Z ((214339.262 467527.3... | 3 |
| חיפה | MULTIPOLYGON Z (((205861.104 70... | 3 |
| ירושלים | POLYGON Z ((192338.062 638359.2... | 3 |
| מרכז | POLYGON Z ((171124.422 631538.1... | 5 |
| צפון | MULTIPOLYGON Z (((204779.953 71... | 2 |
| תל אביב | POLYGON Z ((177310.813 655789.8... | 2 |
# Subtract 1 (self)
districts['count'] = districts['count'] - 1
districts
| geometry | count | |
|---|---|---|
| Machoz | ||
| דרום | POLYGON Z ((214339.262 467527.3... | 2 |
| חיפה | MULTIPOLYGON Z (((205861.104 70... | 2 |
| ירושלים | POLYGON Z ((192338.062 638359.2... | 2 |
| מרכז | POLYGON Z ((171124.422 631538.1... | 4 |
| צפון | MULTIPOLYGON Z (((204779.953 71... | 1 |
| תל אביב | POLYGON Z ((177310.813 655789.8... | 1 |
# Plot
districts.plot(column='count', cmap='Spectral', edgecolor='grey', legend=True);
Exercise 09-h#
import pandas as pd
import geopandas as gpd
# Read & subset columns
stations = gpd.read_file('data/RAIL_STAT_ONOFF_MONTH.shp')
stat = gpd.read_file('data/statisticalareas_demography2019.gdb')
stations = stations[['STAT_NAMEH', 'geometry']]
stat = stat[['SHEM_YISHUV', 'Pop_Total', 'geometry']]
# Calculate area
stat['area'] = stat.area
# Replace "No Data" population with 0
stat['Pop_Total'] = stat['Pop_Total'].fillna(0)
# Buffer stations
stations['geometry'] = stations.buffer(1000)
# Intersect
stations = gpd.overlay(stations, stat, how="intersection")
# Calculate population per polygon
stations['area1'] = stations.area
stations['prop'] = stations['area1'] / stations['area']
stations['pop1'] = stations['Pop_Total'] * stations['prop']
# Dissolve
stations = stations[['STAT_NAMEH', 'pop1', 'geometry']].dissolve(by='STAT_NAMEH', aggfunc='sum').reset_index()
# Sort
stations = stations.sort_values(by='pop1', ascending=False)
stations
| STAT_NAMEH | geometry | pop1 | |
|---|---|---|---|
| 11 | בת ים יוספטל | POLYGON ((177152.750 657153.746... | 55761.058635 |
| 24 | ירושלים יצחק נבון | POLYGON ((218650.757 632087.423... | 46441.448009 |
| 62 | ת"א סבידור מרכז | POLYGON ((180306.847 665089.717... | 39493.940946 |
| 16 | חולון וולפסון | POLYGON ((176806.558 659551.427... | 35295.776253 |
| 48 | קרית מוצקין | POLYGON ((206654.795 747716.120... | 35148.940155 |
| ... | ... | ... | ... |
| 32 | מגדל העמק כפר ברוך | POLYGON ((218894.991 727874.541... | 29.732200 |
| 49 | קרית מלאכי יואב | POLYGON ((184266.430 628244.915... | 0.000000 |
| 17 | חוצות המפרץ | POLYGON ((205604.630 745103.145... | 0.000000 |
| 42 | פאתי מודיעין | POLYGON ((197394.885 644447.483... | 0.000000 |
| 23 | יקנעם כפר יהושוע | POLYGON ((212300.180 731084.998... | 0.000000 |
63 rows × 3 columns
Exercise 09-i#
# Read & subset columns
towns = gpd.read_file('data/muni_il.shp', encoding='utf-8')
rail = gpd.read_file('data/RAIL_STRATEGIC.shp')
towns = towns[['Muni_Heb', 'geometry']]
rail = rail[['ISACTIVE', 'geometry']]
# Subset railways
sel = rail['ISACTIVE'] == 'פעיל'
rail = rail[sel]
# Dissolve towns
towns = towns.dissolve(by='Muni_Heb').reset_index()
# Union railway lines
rail1 = rail.unary_union
# Subset towns
sel = towns.intersects(rail1)
towns = towns[sel]
# Plot
base = towns.plot(color='lightgrey', edgecolor='grey')
rail.plot(ax=base, color='black');
