GeoPandas

Read File

import geopandas

path_to_data = geopandas.datasets.get_path("nybb")
gdf = geopandas.read_file(path_to_data)

Write File

gdf.to_file("my_file.shp") # save Shapely
gdf.to_file("my_file.geojson", driver="GeoJSON") # save JSON

Measuring area

gdf = gdf.set_index("BoroName")
gdf["area"] = gdf.area

gdf.head()

	BoroCode	Shape_Leng	Shape_Area	geometry	area
BoroName
Staten Island	5	330470.010332	1.623820e+09	MULTIPOLYGON (((970217.022 145643.332, 970227....	1.623822e+09
Queens	4	896344.047763	3.045213e+09	MULTIPOLYGON (((1029606.077 156073.814, 102957...	3.045214e+09
Brooklyn	3	741080.523166	1.937479e+09	MULTIPOLYGON (((1021176.479 151374.797, 102100...	1.937478e+09
Manhattan	1	359299.096471	6.364715e+08	MULTIPOLYGON (((981219.056 188655.316, 980940....	6.364712e+08
Bronx	2	464392.991824	1.186925e+09	MULTIPOLYGON (((1012821.806 229228.265, 101278...	1.186926e+09

Getting polygon boundary and centroid

gdf['boundary'] = gdf.boundary # LineString

gdf['boundary']['Manhattan']

gdf['centroid'] = gdf.centroid

gdf['centroid']['Manhattan']

Measuring distance

gdf['centroid']
gdf['centroid']['Manhattan'].distance(gdf['centroid']['Queens'])

48401.272479243045

Making maps

gdf.plot()

<AxesSubplot:>

gdf.plot("area", legend=True) # plot the area measured

<AxesSubplot:>

gdf = gdf.set_geometry("boundary") # switching the active geometry to boundary
gdf.plot("area", legend=True)

<AxesSubplot:>

# add layers by using the same plot axis
ax = gdf.plot("area", legend=True)
gdf = gdf.set_geometry("centroid")
gdf.plot("area", ax=ax)

<AxesSubplot:>

gdf = gdf.set_geometry("geometry")
gdf["buffered"] = gdf.buffer(10000) # buffering the active geometry by 10 000 feet
gdf["buffered_centroid"] = gdf["centroid"].buffer(10000)

ax = gdf["buffered"].plot(alpha=.5)
gdf["buffered_centroid"].plot(ax=ax, color="red", alpha=.5)
gdf["boundary"].plot(ax=ax, color="white", linewidth=.5)

<AxesSubplot:>

# use zorder to control the order of layers
ax = gdf["buffered"].plot(alpha=.5, zorder=2)
gdf["buffered_centroid"].plot(ax=ax, color="red", alpha=.5, zorder = 1)
gdf["boundary"].plot(ax=ax, color="white", linewidth=.5, zorder=3)

<AxesSubplot:>

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
cities = geopandas.read_file(geopandas.datasets.get_path('naturalearth_cities'))

import matplotlib.pyplot as plt
world = world[(world.pop_est>0) & (world.name!="Antarctica")]
world['gdp_per_cap'] = world.gdp_md_est / world.pop_est

fig, ax = plt.subplots(1, 1)

from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
world.plot(column='gdp_per_cap', ax=ax, legend=True, cax=cax)

<AxesSubplot:>

fig, ax = plt.subplots(1, 1)

world.plot(column='pop_est',ax=ax,legend=True,legend_kwds={'label': "Population by Country",'orientation': "horizontal"})

<AxesSubplot:>

world.plot(column='gdp_per_cap', cmap='OrRd');

# need mapclassify to be installed
# ‘box_plot’, ‘equal_interval’, ‘fisher_jenks’, ‘fisher_jenks_sampled’, ‘headtail_breaks’, ‘jenks_caspall’, ‘jenks_caspall_forced’, ‘jenks_caspall_sampled’, ‘max_p_classifier’, ‘maximum_breaks’, ‘natural_breaks’, ‘quantiles’, ‘percentiles’, ‘std_mean’ or ‘user_defined’
world.plot(column='gdp_per_cap', cmap='OrRd', scheme='box_plot')

<AxesSubplot:>

Missed Data

import numpy as np
world.loc[np.random.choice(world.index, 40), 'pop_est'] = np.nan
world.plot(column='pop_est', figsize=(15, 10), missing_kwds={'color': 'lightgrey'});

world.plot(column="pop_est",legend=True,scheme="quantiles",figsize=(15, 10),missing_kwds={"color": "lightgrey","edgecolor": "red","hatch": "///","label": "Missing values",},);

Geometry relations

brooklyn = gdf.loc["Brooklyn", "geometry"]
brooklyn

gdf["geometry"].intersects(brooklyn)

BoroName
Staten Island    False
Queens            True
Brooklyn          True
Manhattan         True
Bronx            False
dtype: bool

gdf["within"] = gdf["buffered_centroid"].within(gdf)
gdf["within"]

BoroName
Staten Island     True
Queens            True
Brooklyn         False
Manhattan        False
Bronx            False
Name: within, dtype: bool

gdf = gdf.set_geometry("buffered_centroid")
ax = gdf.plot("within", legend=True, categorical=True, legend_kwds={'loc': "upper left"})  # using categorical plot and setting the position of the legend
gdf["boundary"].plot(ax=ax, color="black", linewidth=.5)

<AxesSubplot:>

Display options

geopandas.options

Options(
  display_precision: None [default: None]
      The precision (maximum number of decimals) of the coordinates in the
      WKT representation in the Series/DataFrame display. By default (None),
      it tries to infer and use 3 decimals for projected coordinates and 5
      decimals for geographic coordinates.
  use_pygeos: True [default: True]
      Whether to use PyGEOS to speed up spatial operations. The default is
      True if PyGEOS is installed, and follows the USE_PYGEOS environment
      variable if set.
  )

geopandas.options.display_precision = 2
gdf['geometry']

BoroName
Staten Island    MULTIPOLYGON (((970217.02 145643.33, 970227.22...
Queens           MULTIPOLYGON (((1029606.08 156073.81, 1029578....
Brooklyn         MULTIPOLYGON (((1021176.48 151374.80, 1021002....
Manhattan        MULTIPOLYGON (((981219.06 188655.32, 980940.52...
Bronx            MULTIPOLYGON (((1012821.81 229228.26, 1012785....
Name: geometry, dtype: geometry

Indexing

• cx
• loc
• iloc

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
southern_world = world.cx[:, 0:]
southern_world.plot(figsize=(10, 3))

<AxesSubplot:>

Geometric Manipulations

world.rotate(angle=90, origin='center', use_radians=False).plot()

<AxesSubplot:>

world.translate(xoff=20).plot()

<AxesSubplot:>

Aggregation with dissolve

• ‘first’, ‘last’, ‘min’, ‘max’, ‘sum’, ‘mean’, ‘median’

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
world = world[['continent', 'geometry', 'pop_est']]
continents = world.dissolve(by='continent', aggfunc='sum')
continents.plot(column='pop_est', scheme='quantiles', cmap='YlOrRd')

<AxesSubplot:>

Merging Data

world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
cities = geopandas.read_file(geopandas.datasets.get_path('naturalearth_cities'))

country_shapes = world[['geometry', 'iso_a3']]
country_names = world[['name', 'iso_a3']]

countries = world[['geometry', 'name']]
countries = countries.rename(columns={'name':'country'})
countries.head()

	geometry	country
0	MULTIPOLYGON (((180.00 -16.07, 180.00 -16.56, ...	Fiji
1	POLYGON ((33.90 -0.95, 34.07 -1.06, 37.70 -3.1...	Tanzania
2	POLYGON ((-8.67 27.66, -8.67 27.59, -8.68 27.4...	W. Sahara
3	MULTIPOLYGON (((-122.84 49.00, -122.97 49.00, ...	Canada
4	MULTIPOLYGON (((-122.84 49.00, -120.00 49.00, ...	United States of America

# append
joined = world.geometry.append(cities.geometry)

# merge
country_shapes = country_shapes.merge(country_names, on='iso_a3')
country_shapes.head()

	geometry	iso_a3	name
0	MULTIPOLYGON (((180.00 -16.07, 180.00 -16.56, ...	FJI	Fiji
1	POLYGON ((33.90 -0.95, 34.07 -1.06, 37.70 -3.1...	TZA	Tanzania
2	POLYGON ((-8.67 27.66, -8.67 27.59, -8.68 27.4...	ESH	W. Sahara
3	MULTIPOLYGON (((-122.84 49.00, -122.97 49.00, ...	CAN	Canada
4	MULTIPOLYGON (((-122.84 49.00, -120.00 49.00, ...	USA	United States of America

# sjoin merged based on the spatial relationship
# how, left, right, inner
# op,intersects, within, contains
cities_with_country = geopandas.sjoin(cities, countries, how="inner", op='intersects')
cities_with_country.head()

	name	geometry	index_right	country
0	Vatican City	POINT (12.45 41.90)	141	Italy
1	San Marino	POINT (12.44 43.94)	141	Italy
192	Rome	POINT (12.48 41.90)	141	Italy
2	Vaduz	POINT (9.52 47.13)	114	Austria
184	Vienna	POINT (16.36 48.20)	114	Austria

Geocoding

• converting place names to location on Earth

boros = geopandas.read_file(geopandas.datasets.get_path("nybb"))
boros.BoroName

	BoroCode	BoroName	Shape_Leng	Shape_Area	geometry
0	5	Staten Island	330470.010332	1.623820e+09	MULTIPOLYGON (((970217.02 145643.33, 970227.22...
1	4	Queens	896344.047763	3.045213e+09	MULTIPOLYGON (((1029606.08 156073.81, 1029578....
2	3	Brooklyn	741080.523166	1.937479e+09	MULTIPOLYGON (((1021176.48 151374.80, 1021002....
3	1	Manhattan	359299.096471	6.364715e+08	MULTIPOLYGON (((981219.06 188655.32, 980940.52...
4	2	Bronx	464392.991824	1.186925e+09	MULTIPOLYGON (((1012821.81 229228.26, 1012785....

boro_locations = geopandas.tools.geocode(boros.BoroName)
boro_locations

	geometry	address
0	POINT (-74.08 40.64)	Staten Island, NY, United States
1	POINT (-73.82 40.76)	Queens County, United States
2	POINT (-73.99 40.69)	Brooklyn, NY, United States
3	POINT (-73.97 40.78)	Manhattan, New York, United States
4	POINT (-73.92 40.83)	Bronx County, United States

fig, ax = plt.subplots()
boros.to_crs("EPSG:4326").plot(ax=ax, color="white", edgecolor="black");
boro_locations.plot(ax=ax, color="red");

Reference

• API Reference
• Reference
• GeoPandas
• GeoPandas Documentation