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');
_images/8f7e93fede451defcbc599fd57410cceca956394a1c83729415c9ec011b308ac.png

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);
_images/153dbada486a72d9b51bae5affe423997768a310be964732a54de50329c56d07.png

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');
_images/0d05b933134af25ce135dcf42eb9f84009583ab7adc438339cc666227ff304a5.png