import geopandas as gpd
from pathlib import PathSpatial Analysis using Geopandas
Importing data with Geopandas
input_path = '/home/thulasiram/personal/going_deep_and_wide/togo/togo-targeting-replication/data/shapefiles/cantons.geojson'
data = gpd.read_file(input_path)type(data)geopandas.geodataframe.GeoDataFramedata.head()| canton | poverty | geometry | |
|---|---|---|---|
| 0 | 1 | 3.738084 | MULTIPOLYGON (((0.75228 6.83786, 0.75137 6.840... |
| 1 | 2 | 7.096286 | MULTIPOLYGON (((0.69026 6.80602, 0.69627 6.806... |
| 2 | 3 | 0.824586 | MULTIPOLYGON (((0.63102 6.74430, 0.63295 6.747... |
| 3 | 4 | 3.983729 | MULTIPOLYGON (((0.67259 6.85123, 0.67714 6.849... |
| 4 | 5 | 7.708810 | MULTIPOLYGON (((0.75269 6.84116, 0.75137 6.840... |
print("Number of rows",len(data["canton"]))
print("Number of classes",data["canton"].nunique())Number of rows 387
Number of classes 387Plotting the data on the Map
data.plot()<AxesSubplot: >
Checking the shape and area of the first Multipolygon in the data
data.at[0,"geometry"]
# Calculating area with lat and long is wrong, to calculate
# area correctly we need to change the co-ordinate reference system
round(data.at[0,"geometry"].area,3)0.014Creating an interactive map
data.plot("poverty",legend=True)<AxesSubplot: >
data.explore("poverty",legend=True)Make this Notebook Trusted to load map: File -> Trust Notebook
Geocoding
Geocoding is the process of transforming place names or addresses into coordinates
# Import necessary modules
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point# Filepath
fp = '../data/addresses.txt'
# Read the data
data = pd.read_csv(fp,sep=";")len(data)5data.head()| id | addr | |
|---|---|---|
| 0 | 1000 | Boat Club Road, 411001, Pune, Maharastra |
| 1 | 1001 | Koregaon, 415501, Pune, Maharastra |
| 2 | 1002 | Kothrud, 411038, Pune, Maharastra |
| 3 | 1003 | Balewadi, 411045, Pune, Maharastra |
| 4 | 1004 | Baner, 411047, Pune, Maharastra |
# Import the geocoding tool
from geopandas.tools import geocode
# Geocode addresses using Nominatim. Remember to provide a custom "application name" in the user_agent parameter!
geo = geocode(data["addr"], provider="nominatim", user_agent="autogis_xx", timeout=4)geo.head()| geometry | address | |
|---|---|---|
| 0 | POINT (73.87826 18.53937) | Boat Club Road, Pune City, Pune, Maharashtra, ... |
| 1 | POINT (73.89299 18.53772) | Koregaon Park, Suyojan Society, Ghorpuri, Pune... |
| 2 | POINT (73.80767 18.50389) | Kothrud, Pune City, Maharashtra, 411038, India |
| 3 | POINT (73.76912 18.57767) | Prakashgad Society, Balewadi, Perfect 10 Inter... |
| 4 | POINT (73.77686 18.56424) | Baner, Pune City, Maharashtra, 511045, India |
# Joining the original dataframe with geoencoded dataframe
join = geo.join(data)
join = join.drop(columns=['address'])join.head()| geometry | id | addr | |
|---|---|---|---|
| 0 | POINT (73.87826 18.53937) | 1000 | Boat Club Road, 411001, Pune, Maharastra |
| 1 | POINT (73.89299 18.53772) | 1001 | Koregaon, 415501, Pune, Maharastra |
| 2 | POINT (73.80767 18.50389) | 1002 | Kothrud, 411038, Pune, Maharastra |
| 3 | POINT (73.76912 18.57767) | 1003 | Balewadi, 411045, Pune, Maharastra |
| 4 | POINT (73.77686 18.56424) | 1004 | Baner, 411047, Pune, Maharastra |
Explore the area obtained from Geocoding on the map
from shapely.geometry import box
minx = 73.76
miny = 18.537
maxx = 73.89
maxy = 18.56
geom = box(minx, miny, maxx, maxy)
clipping_gdf = gpd.GeoDataFrame({"geometry": [geom]}, index=[0], crs="epsg:4326")
# Explore the extent on a map
clipping_gdf.explore()Make this Notebook Trusted to load map: File -> Trust Notebook
# Write the data to shape file
outfp = r"../data/addresses.shp"
join.to_file(outfp)