Point in Polygon Using Geopandas
Tuseday 06 - October - 2020 (trip to Kaolack, Madina Baye)
### Point in Polygon Using Geopandas
""" Previously, we sow how to find point in polygon or polygon in point with shapely by using
- .within()
- .contains()
Now, we will see how to use that techniques with Geopandas. We have a file called addresses,
our file contains a data saved in .shp format, the addresses are places located in Southern district of Helsinki.
"""
# Upload the file, read the addresses before going to other step
import geopandas as gpd
source_file = 'D:/Programm files/python/Scripts/Automatic Gis precess/data/addresses.shp'
addresses = gpd.read_file(source_file)
"""
>>> addresses.head()
address ... geometry
0 Kampinkuja 1, 00100 Helsinki, Finland ... POINT (24.93017 60.16837)
1 Kaivokatu 8, 00101 Helsinki, Finland ... POINT (24.94189 60.16987)
2 Hermanstads strandsväg 1, 00580 Helsingfors, F... ... POINT (24.97740 60.18736)
3 Itäväylä, 00900 Helsinki, Finland ... POINT (25.09196 60.21448)
4 Tyynenmerenkatu 9, 00220 Helsinki, Finland ... POINT (24.92148 60.15658)
[5 rows x 3 columns]
>>> addresses.columns.values
array(['address', 'id', 'geometry'], dtype=object)
>>> addresses.info()
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 35 entries, 0 to 34
Data columns (total 3 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 address 35 non-null object
1 id 35 non-null int64
2 geometry 35 non-null geometry
dtypes: geometry(1), int64(1), object(1)
memory usage: 624.0+ bytes
>>> addresses.plot()
<AxesSubplot:>
>>> import matplotlib.pyplot as plt
>>> plt.show()
"""
# Reading a KML file in geopandas, first we need to activate it.
gpd.io.file.fiona.drvsupport.supported_drivers
gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'
""" Output:
>>> gpd.io.file.fiona.drvsupport.supported_drivers
{'AeronavFAA': 'r', 'ARCGEN': 'r', 'BNA': 'rw', 'DXF': 'rw', 'CSV': 'raw', 'OpenFileGDB': 'r', 'FlatGeobuf': 'r', 'ESRIJSON': 'r', 'ESRI Shapefile': 'raw', 'GeoJSON': 'raw', 'GeoJSONSeq': 'rw', 'GPKG': 'raw', 'GML': 'rw', 'OGR_GMT': 'rw', 'GPX': 'rw', 'GPSTrackMaker': 'rw', 'Idrisi': 'r', 'MapInfo File': 'raw', 'DGN': 'raw', 'PCIDSK': 'rw', 'OGR_PDS': 'r', 'S57': 'r', 'SEGY': 'r', 'SUA': 'r', 'TopoJSON': 'r', 'KML': 'rw'}
"""
source_kml_file = 'D:/Studying/GEOSPATIAL TUTORIALS/site-master/data/PKS_suuralue.kml'
polys = gpd.read_file(source_kml_file, driver='KML')
""" Output
# The names of columns of our data that contains 23 districts of Southern Helsinki
>>> polys.columns.values
array(['Name', 'Description', 'geometry'], dtype=object)
>>> polys.info() # A lot of information about our 23 districts in Southern Helsinki
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 23 entries, 0 to 22
Data columns (total 3 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 Name 23 non-null object
1 Description 23 non-null object
2 geometry 23 non-null geometry
dtypes: geometry(1), object(2)
memory usage: 340.0+ bytes
>>> import matplotlib.pyplot as plt
>>> polys.plot() # Visualizing the 23 districts layers
<AxesSubplot:>
>>> plt.show()
"""
# Select the Etelainen district and see its location on map
import matplotlib.pyplot as plt
southern = polys.loc[polys['Name']=='Eteläinen']
# If we print southern we get its index 10 in polys data
"""
>>> southern
Name Description geometry
10 Eteläinen POLYGON Z ((24.78277 60.09997 0.00000, 24.8197...
"""
# Here, we can reset the index and keep it normal
southern.reset_index(drop=True, inplace=True)
"""
>>> southern.head()
Name Description geometry
0 Eteläinen POLYGON Z ((24.78277 60.09997 0.00000, 24.8197...
"""
# Displaying our south district
"""
>>> fig, ax = plt.subplots()
>>> polys.plot(ax=ax, facecolor='gray')
<AxesSubplot:>
>>> addresses.plot(ax=ax, color='blue', markersize=5)
>>> southern.plot(ax=ax, facecolor='red')
>>> southern
Name Description geometry
10 Eteläinen POLYGON Z ((24.78277 60.09997 0.00000, 24.8197...
>>> polys.plot()
>>> plt.show()
We can absolutely see certain points (blue color) are within the selected red polygon.
These blues points are the addresses of southern locations in Helsinki.
Note: Our main location here is Eteläinen
"""
# As we see some blue points inside the red district, let's see the PIP or Point in Polygon to find those points
from shapely import speedups
speedups.enabled
speedups.enabled
""" output: True """
# Let's check which Points (within blue points) are within southern Polygon.
blue_points_in_red = addresses.within(southern.at[0, 'geometry'])
""" output
>>> blue_points_in_red.head()
0 True
1 True
2 False
3 False
4 True
dtype: bool
>>> blue_points_in_red.value_counts()
False 27
True 8
dtype: int64
If the Point was inside the Polygon means True (then we have 8 True Points), else if the Point was
False means tt's not within the Polygon (then we have 27 False)
"""
# Retrieve those 8 True Points inside the Polygon
points_in_polygon = addresses.loc[blue_points_in_red]
""" Output
>>> points_in_polygon
address ... geometry
0 Kampinkuja 1, 00100 Helsinki, Finland ... POINT (24.93017 60.16837)
1 Kaivokatu 8, 00101 Helsinki, Finland ... POINT (24.94189 60.16987)
4 Tyynenmerenkatu 9, 00220 Helsinki, Finland ... POINT (24.92148 60.15658)
10 Rautatientori 1, 00100 Helsinki, Finland ... POINT (24.94251 60.17119)
30 Urho Kekkosen katu 1, 00100 Helsinki, Finland ... POINT (24.93376 60.16948)
31 Gräsviksgatan 17, 00101 Helsingfors, Finland ... POINT (24.92501 60.16500)
32 Stillahavsgatan 3, 00220 Helsingfors, Finland ... POINT (24.92140 60.15907)
33 Vilhelmsgatan 4, 00101 Helsingfors, Finland ... POINT (24.94685 60.17191)
[8 rows x 3 columns]
>>> points_in_polygon.plot()
<AxesSubplot:>
>>> plt.show()
>>> points_in_polygon.to_excel('8 Points in addresses point in polygon.xlsx')
"""
# Display only the 8 points in polygon
fig, ax = plt.subplots()
polys.plot(ax=ax, facecolor='gray')
southern.plot(ax=ax, facecolor='red')
points_in_polygon.plot(ax=ax, color='gold', markersize=2)
plt.tight_layout()
plt.show()
Comments
Post a Comment