From “Geometric operations”
Contents
From “Geometric operations”#
Last updated: 2023-02-25 13:37:27
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
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
/tmp/ipykernel_12589/1772616394.py in <module>
1 # Spatial join with self
----> 2 districts = gpd.sjoin(
3 districts[['Machoz', 'geometry']],
4 districts[['count', 'geometry']],
5 how='left',
~/.local/lib/python3.10/site-packages/geopandas/tools/sjoin.py in sjoin(left_df, right_df, how, predicate, lsuffix, rsuffix, **kwargs)
122 _basic_checks(left_df, right_df, how, lsuffix, rsuffix)
123
--> 124 indices = _geom_predicate_query(left_df, right_df, predicate)
125
126 joined = _frame_join(indices, left_df, right_df, how, lsuffix, rsuffix)
~/.local/lib/python3.10/site-packages/geopandas/tools/sjoin.py in _geom_predicate_query(left_df, right_df, predicate)
214 # all other predicates are symmetric
215 # keep them the same
--> 216 sindex = right_df.sindex
217 input_geoms = left_df.geometry
218
~/.local/lib/python3.10/site-packages/geopandas/base.py in sindex(self)
2704 [2]])
2705 """
-> 2706 return self.geometry.values.sindex
2707
2708 @property
~/.local/lib/python3.10/site-packages/geopandas/array.py in sindex(self)
289 def sindex(self):
290 if self._sindex is None:
--> 291 self._sindex = _get_sindex_class()(self.data)
292 return self._sindex
293
~/.local/lib/python3.10/site-packages/geopandas/sindex.py in _get_sindex_class()
19 if compat.HAS_RTREE:
20 return RTreeIndex
---> 21 raise ImportError(
22 "Spatial indexes require either `rtree` or `pygeos`. "
23 "See installation instructions at https://geopandas.org/install.html"
ImportError: Spatial indexes require either `rtree` or `pygeos`. See installation instructions at https://geopandas.org/install.html
# 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
| SHEM_YISHUV | Pop_Total | geometry | area | |
|---|---|---|---|---|
| 0 | המעפיל | 869.0 | MULTIPOLYGON (((198790.704 698580.560, 198811.... | 4.195856e+05 |
| 1 | משגב עם | 339.0 | MULTIPOLYGON (((251692.761 794858.649, 251719.... | 2.482362e+05 |
| 2 | גאולים | 968.0 | MULTIPOLYGON (((194649.280 690123.970, 194803.... | 2.274099e+06 |
| 3 | להבות הבשן | 915.0 | MULTIPOLYGON (((260535.621 783379.248, 260519.... | 7.285277e+05 |
| 4 | מכמורת | 1403.0 | MULTIPOLYGON (((188690.124 701820.428, 188707.... | 6.281249e+05 |
| ... | ... | ... | ... | ... |
| 3190 | קריית נטפים | 958.0 | MULTIPOLYGON (((210706.963 669359.635, 210710.... | 7.847603e+03 |
| 3191 | דולב | 1448.0 | MULTIPOLYGON (((212744.459 648214.855, 212747.... | 7.847603e+03 |
| 3192 | עתניאל | 1044.0 | MULTIPOLYGON (((202798.088 594213.306, 202801.... | 7.847603e+03 |
| 3193 | יצהר | 1726.0 | MULTIPOLYGON (((222570.662 675002.822, 222574.... | 7.847603e+03 |
| 3194 | מבואות יריחו | NaN | MULTIPOLYGON (((239670.270 645948.661, 239673.... | 7.847603e+03 |
3195 rows × 4 columns
# Replace "No Data" population with 0
stat['Pop_Total'] = stat['Pop_Total'].fillna(0)
| SHEM_YISHUV | Pop_Total | geometry | area | |
|---|---|---|---|---|
| 0 | המעפיל | 869.0 | MULTIPOLYGON (((198790.704 698580.560, 198811.... | 4.195856e+05 |
| 1 | משגב עם | 339.0 | MULTIPOLYGON (((251692.761 794858.649, 251719.... | 2.482362e+05 |
| 2 | גאולים | 968.0 | MULTIPOLYGON (((194649.280 690123.970, 194803.... | 2.274099e+06 |
| 3 | להבות הבשן | 915.0 | MULTIPOLYGON (((260535.621 783379.248, 260519.... | 7.285277e+05 |
| 4 | מכמורת | 1403.0 | MULTIPOLYGON (((188690.124 701820.428, 188707.... | 6.281249e+05 |
| ... | ... | ... | ... | ... |
| 3190 | קריית נטפים | 958.0 | MULTIPOLYGON (((210706.963 669359.635, 210710.... | 7.847603e+03 |
| 3191 | דולב | 1448.0 | MULTIPOLYGON (((212744.459 648214.855, 212747.... | 7.847603e+03 |
| 3192 | עתניאל | 1044.0 | MULTIPOLYGON (((202798.088 594213.306, 202801.... | 7.847603e+03 |
| 3193 | יצהר | 1726.0 | MULTIPOLYGON (((222570.662 675002.822, 222574.... | 7.847603e+03 |
| 3194 | מבואות יריחו | 0.0 | MULTIPOLYGON (((239670.270 645948.661, 239673.... | 7.847603e+03 |
3195 rows × 4 columns
# Buffer stations
stations['geometry'] = stations.buffer(1000)
| STAT_NAMEH | geometry | |
|---|---|---|
| 0 | ראשון לציון משה דיין | POLYGON ((178206.094 655059.936, 178201.279 65... |
| 1 | קרית ספיר נתניה | POLYGON ((188592.123 687587.598, 188587.308 68... |
| 2 | ת"א השלום | POLYGON ((181621.940 664537.210, 181617.125 66... |
| 3 | שדרות | POLYGON ((161722.849 602798.889, 161718.033 60... |
| 4 | רמלה | POLYGON ((189508.910 648466.870, 189504.095 64... |
| ... | ... | ... |
| 58 | מגדל העמק כפר ברוך | POLYGON ((220851.931 728164.825, 220847.116 72... |
| 59 | באר שבע אוניברסיטה | POLYGON ((182813.614 574550.411, 182808.798 57... |
| 60 | חיפה בת גלים | POLYGON ((199610.010 748443.790, 199605.195 74... |
| 61 | לוד | POLYGON ((189315.977 650289.373, 189311.162 65... |
| 62 | נהריה | POLYGON ((210570.590 767769.220, 210565.775 76... |
63 rows × 2 columns
# Intersect
stations = gpd.overlay(stations, stat, how="intersection")
| STAT_NAMEH | SHEM_YISHUV | Pop_Total | area | geometry | |
|---|---|---|---|---|---|
| 0 | ראשון לציון משה דיין | חולון | 0.0 | 4.161271e+05 | POLYGON ((177401.184 656040.721, 177496.379 65... |
| 1 | בת ים קוממיות | חולון | 0.0 | 4.161271e+05 | POLYGON ((177911.240 656463.090, 177888.670 65... |
| 2 | ראשון לציון משה דיין | בת ים | 1549.0 | 1.667471e+06 | POLYGON ((176734.697 655941.857, 176823.411 65... |
| 3 | בת ים יוספטל | בת ים | 1549.0 | 1.667471e+06 | POLYGON ((176702.235 657724.386, 176678.391 65... |
| 4 | בת ים קוממיות | בת ים | 1549.0 | 1.667471e+06 | POLYGON ((176480.670 656118.579, 176447.609 65... |
| ... | ... | ... | ... | ... | ... |
| 584 | נהריה | נהרייה | 4513.0 | 4.696705e+05 | POLYGON ((208797.580 768403.613, 208863.483 76... |
| 585 | נהריה | נהרייה | 3947.0 | 3.870092e+05 | POLYGON ((208646.710 768151.903, 208688.669 76... |
| 586 | נהריה | נהרייה | 3228.0 | 4.423230e+05 | POLYGON ((208575.405 767867.237, 208589.805 76... |
| 587 | נהריה | נהרייה | 5745.0 | 4.866801e+05 | POLYGON ((208739.120 767213.650, 208688.669 76... |
| 588 | נהריה | נהרייה | 7076.0 | 1.308312e+06 | POLYGON ((209375.500 766788.435, 209280.305 76... |
589 rows × 5 columns
# Calculate population per polygon
stations['area1'] = stations.area
stations['prop'] = stations['area1'] / stations['area']
stations['pop1'] = stations['Pop_Total'] * stations['prop']
| STAT_NAMEH | SHEM_YISHUV | Pop_Total | area | geometry | area1 | prop | pop1 | |
|---|---|---|---|---|---|---|---|---|
| 0 | ראשון לציון משה דיין | חולון | 0.0 | 4.161271e+05 | POLYGON ((177401.184 656040.721, 177496.379 65... | 119577.976361 | 0.287359 | 0.000000 |
| 1 | בת ים קוממיות | חולון | 0.0 | 4.161271e+05 | POLYGON ((177911.240 656463.090, 177888.670 65... | 416127.085983 | 1.000000 | 0.000000 |
| 2 | ראשון לציון משה דיין | בת ים | 1549.0 | 1.667471e+06 | POLYGON ((176734.697 655941.857, 176823.411 65... | 165230.066812 | 0.099090 | 153.490774 |
| 3 | בת ים יוספטל | בת ים | 1549.0 | 1.667471e+06 | POLYGON ((176702.235 657724.386, 176678.391 65... | 45069.519276 | 0.027029 | 41.867413 |
| 4 | בת ים קוממיות | בת ים | 1549.0 | 1.667471e+06 | POLYGON ((176480.670 656118.579, 176447.609 65... | 765448.317038 | 0.459048 | 711.064620 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 584 | נהריה | נהרייה | 4513.0 | 4.696705e+05 | POLYGON ((208797.580 768403.613, 208863.483 76... | 352598.187026 | 0.750735 | 3388.068019 |
| 585 | נהריה | נהרייה | 3947.0 | 3.870092e+05 | POLYGON ((208646.710 768151.903, 208688.669 76... | 324775.549789 | 0.839193 | 3312.296416 |
| 586 | נהריה | נהרייה | 3228.0 | 4.423230e+05 | POLYGON ((208575.405 767867.237, 208589.805 76... | 427565.966690 | 0.966637 | 3120.305430 |
| 587 | נהריה | נהרייה | 5745.0 | 4.866801e+05 | POLYGON ((208739.120 767213.650, 208688.669 76... | 387469.851531 | 0.796149 | 4573.875376 |
| 588 | נהריה | נהרייה | 7076.0 | 1.308312e+06 | POLYGON ((209375.500 766788.435, 209280.305 76... | 160615.106359 | 0.122765 | 868.686321 |
589 rows × 8 columns
# Dissolve
stations = stations[['STAT_NAMEH', 'pop1', 'geometry']].dissolve(by='STAT_NAMEH', aggfunc='sum').reset_index()
| STAT_NAMEH | geometry | pop1 | |
|---|---|---|---|
| 0 | אופקים | POLYGON ((165365.894 580214.219, 165268.821 58... | 87.255153 |
| 1 | אחיהוד | POLYGON ((216050.708 756566.100, 216041.383 75... | 499.515455 |
| 2 | אשדוד עד הלום | POLYGON ((168257.270 630397.885, 168165.251 63... | 1980.780315 |
| 3 | אשקלון | POLYGON ((163065.155 619740.598, 162976.442 61... | 688.438355 |
| 4 | באר יעקב | POLYGON ((183715.588 647952.959, 183620.394 64... | 8294.583098 |
| ... | ... | ... | ... |
| 58 | שער צומת חולון | POLYGON ((178470.478 659675.076, 178391.655 65... | 14188.762474 |
| 59 | ת"א אוניברסיטה | POLYGON ((181155.300 667044.540, 181076.477 66... | 9160.396912 |
| 60 | ת"א ההגנה | POLYGON ((179349.547 661510.335, 179347.178 66... | 23976.189521 |
| 61 | ת"א השלום | POLYGON ((180523.923 663542.025, 180426.850 66... | 25821.530059 |
| 62 | ת"א סבידור מרכז | POLYGON ((180306.847 665089.717, 180277.411 66... | 39493.940946 |
63 rows × 3 columns
# Sort
stations = stations.sort_values(by='pop1', ascending=False)
stations
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');
