From “Geometric operations”

Last updated: 2022-06-22 18:51:36

Exercise 09-b

import geopandas as gpd
towns = gpd.read_file("data/muni_il.shp")
towns = towns[["Muni_Heb", "Muni_Eng", "Machoz", "AreaSQM", "geometry"]]
towns["area"] = towns.area
towns1 = towns.sort_values("area")
towns1["Muni_Heb"].iloc[0]
/home/michael/.local/lib/python3.8/site-packages/geopandas/_compat.py:111: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.10.1-CAPI-1.16.0). Conversions between both will be slow.
  warnings.warn(
'ללא שיפוט - אזור נחל ארבל'
towns1["Muni_Heb"].iloc[towns1.shape[0]-1]
'רמת נגב'

Exercise 09-c

import geopandas as gpd
towns = gpd.read_file("data/muni_il.shp")
towns = towns[["Muni_Heb", "Muni_Eng", "Machoz", "AreaSQM", "geometry"]]
towns["area"] = towns.area
towns.plot.scatter("area", "AreaSQM");
_images/exercise_solutions_geopandas2_9_0.png

Exercise 09-g

import geopandas as gpd
# Read data & subset columns
towns = gpd.read_file("data/muni_il.shp")
towns = towns[["Machoz", "geometry"]]
towns
Machoz geometry
0 מרכז POLYGON ((196372.646 684381.632...
1 מרכז POLYGON ((199366.832 687791.637...
2 דרום POLYGON ((223717.609 523705.241...
3 דרום POLYGON ((179909.061 607123.385...
4 דרום POLYGON ((169289.462 601344.604...
... ... ...
422 צפון POLYGON ((250468.035 779727.733...
423 צפון POLYGON ((251434.056 771376.157...
424 צפון POLYGON ((241900.500 776693.188...
425 צפון POLYGON ((254712.930 767150.100...
426 ירושלים POLYGON ((209273.654 635045.173...

427 rows × 2 columns

# Dissolve
districts = towns.dissolve(by="Machoz").reset_index()
districts
Machoz geometry
0 דרום POLYGON ((214086.053 466729.572...
1 חיפה MULTIPOLYGON (((188127.359 7026...
2 ירושלים POLYGON ((209724.750 637761.125...
3 מרכז POLYGON ((172850.921 635138.036...
4 צפון MULTIPOLYGON (((222775.875 7154...
5 תל אביב POLYGON ((183862.073 659328.831...
# Add 'count' variable
districts['count'] = 1
districts
Machoz geometry count
0 דרום POLYGON ((214086.053 466729.572... 1
1 חיפה MULTIPOLYGON (((188127.359 7026... 1
2 ירושלים POLYGON ((209724.750 637761.125... 1
3 מרכז POLYGON ((172850.921 635138.036... 1
4 צפון MULTIPOLYGON (((222775.875 7154... 1
5 תל אביב POLYGON ((183862.073 659328.831... 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 ((214086.053 466729.572... 1
0 דרום POLYGON ((214086.053 466729.572... 1
0 דרום POLYGON ((214086.053 466729.572... 1
1 חיפה MULTIPOLYGON (((188127.359 7026... 1
1 חיפה MULTIPOLYGON (((188127.359 7026... 1
... ... ... ...
3 מרכז POLYGON ((172850.921 635138.036... 1
4 צפון MULTIPOLYGON (((222775.875 7154... 1
4 צפון MULTIPOLYGON (((222775.875 7154... 1
5 תל אביב POLYGON ((183862.073 659328.831... 1
5 תל אביב POLYGON ((183862.073 659328.831... 1

18 rows × 3 columns

# Aggregate
districts = districts.dissolve(by='Machoz', aggfunc='sum')
districts
geometry count
Machoz
דרום POLYGON ((213816.718 465853.995... 3
חיפה MULTIPOLYGON (((192270.440 7241... 3
ירושלים POLYGON ((209753.703 637767.375... 3
מרכז POLYGON ((172832.602 635073.264... 5
צפון MULTIPOLYGON (((222662.984 7155... 2
תל אביב POLYGON ((183823.783 659264.058... 2
# Subtract 1 (self)
districts['count'] = districts['count'] - 1
districts
geometry count
Machoz
דרום POLYGON ((213816.718 465853.995... 2
חיפה MULTIPOLYGON (((192270.440 7241... 2
ירושלים POLYGON ((209753.703 637767.375... 2
מרכז POLYGON ((172832.602 635073.264... 4
צפון MULTIPOLYGON (((222662.984 7155... 1
תל אביב POLYGON ((183823.783 659264.058... 1
# Plot
districts.plot(column='count', cmap='Spectral', edgecolor='grey', legend=True);
_images/exercise_solutions_geopandas2_19_0.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")
rail = gpd.read_file("data/RAIL_STRATEGIC.shp")
towns = towns[["Muni_Heb", "geometry"]]
rail = rail[["geometry", "ISACTIVE"]]
# 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/exercise_solutions_geopandas2_39_0.png