From “Geometric operations”
Contents
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");
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);
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");
