Geospatial Vector Data in Python

14. Geospatial Vector Data in Python#

Attribution: Parts of this tutorial are developed based on the content from the following great sources: Vector data in Python; and Introduction to GeoPandas.

In this lecture, you will learn how to interact with geospatial data in Python. Our focus is on smaller datasets in this lecture, and in the next one we will learn how to handle large datasets for scalable analysis.

We will use the GeoPandas package to open, manipulate and write vector datasets.

14.1. Intro to GeoPandas#

GeoPandas extends the popular pandas library for data analysis to geospatial applications. The main pandas objects (the Series and the DataFrame) are expanded to GeoPandas objects (GeoSeries and GeoDataFrame). This extension is implemented by including geometric types, represented in Python using the shapely library, and by providing dedicated methods for spatial operations (union, intersection, etc.). The relationship between Series, DataFrame, GeoSeries and GeoDataFrame can be briefly explained as follow:

A Series is a one-dimensional array with axis, holding any data type (integers, strings, floating-point numbers, Python objects, etc.)
A DataFrame is a two-dimensional labeled data structure with columns of potentially different types.
A GeoSeries is a Series object designed to store shapely geometry objects.
A GeoDataFrame is an extened pandas.DataFrame, which has a column with geometry objects, and this column is a GeoSeries.

Each GeoSeries can contain any geometry type (you can even mix them within a single array) and has a GeoSeries.crs attribute, which stores information about the projection. Therefore, each GeoSeries in a GeoDataFrame can be in a different projection, allowing you to have, for example, multiple versions (different projections) of the same geometry.

Note: Only one GeoSeries in a GeoDataFrame is considered the active geometry, which means that all geometric operations applied to a GeoDataFrame operate on this active column.

import matplotlib.pyplot as plt
import geopandas as gpd
from cartopy import crs as ccrs
import geodatasets
import pandas as pd

14.2. Reading Files#

nybb_gdf = gpd.read_file(geodatasets.get_path("ny.bb"))

nybb_gdf.info

<bound method DataFrame.info of    BoroCode       BoroName     Shape_Leng    Shape_Area  \
       5  Staten Island  330470.010332  1.623820e+09   
       4         Queens  896344.047763  3.045213e+09   
       3       Brooklyn  741080.523166  1.937479e+09   
       1      Manhattan  359299.096471  6.364715e+08   
       2          Bronx  464392.991824  1.186925e+09   

                                            geometry  
MULTIPOLYGON (((970217.022 145643.332, 970227....  
MULTIPOLYGON (((1029606.077 156073.814, 102957...  
MULTIPOLYGON (((1021176.479 151374.797, 102100...  
MULTIPOLYGON (((981219.056 188655.316, 980940....  
MULTIPOLYGON (((1012821.806 229228.265, 101278...  >

You can also select to read only parts of the data using the bbox, rows or mask argument:

nybb_partial_gdf = gpd.read_file(geodatasets.get_path("ny.bb"), rows=2)

nybb_partial_gdf

	BoroCode	BoroName	Shape_Leng	Shape_Area	geometry
0	5	Staten Island	330470.010332	1.623820e+09	MULTIPOLYGON (((970217.022 145643.332, 970227....
1	4	Queens	896344.047763	3.045213e+09	MULTIPOLYGON (((1029606.077 156073.814, 102957...

14.3. Writing Files#

You can write any GeoDataFrame to local disk and set the format using the driver argument.

nybb_gdf.to_file("nybb.geojson", driver="GeoJSON")

14.4. Constructing a GeoDataFrame manually#

from shapely.geometry import Point

points_gdf = gpd.GeoDataFrame({
    'geometry': [Point(1, 1), Point(2, 2)],
    'attribute1': [1, 2],
    'attribute2': [0.1, 0.2]})

points_gdf

	geometry	attribute1	attribute2
0	POINT (1.00000 1.00000)	1	0.1
1	POINT (2.00000 2.00000)	2	0.2

14.5. Creating a GeoDataFrame from an existing dataframe#

cities_df = pd.DataFrame(
    {'City': ['Buenos Aires', 'Brasilia', 'Santiago', 'Bogota', 'Caracas'],
     'Country': ['Argentina', 'Brazil', 'Chile', 'Colombia', 'Venezuela'],
     'Latitude': [-34.58, -15.78, -33.45, 4.60, 10.48],
     'Longitude': [-58.66, -47.91, -70.66, -74.08, -66.86]})

cities_gdf = gpd.GeoDataFrame(
    cities_df, geometry=gpd.points_from_xy(cities_df.Longitude, cities_df.Latitude))

cities_gdf

	City	Country	Latitude	Longitude	geometry
0	Buenos Aires	Argentina	-34.58	-58.66	POINT (-58.66000 -34.58000)
1	Brasilia	Brazil	-15.78	-47.91	POINT (-47.91000 -15.78000)
2	Santiago	Chile	-33.45	-70.66	POINT (-70.66000 -33.45000)
3	Bogota	Colombia	4.60	-74.08	POINT (-74.08000 4.60000)
4	Caracas	Venezuela	10.48	-66.86	POINT (-66.86000 10.48000)

14.6. Working with Attributes#

All the attributes that are defined in shapely for objects, are available in GeoSeries. When you retrieve these from a GeoDataFrame however, you should note that they will be calculated based on the active geometry column.

area: shape area (units of projection – see projections)

bounds: tuple of max and min coordinates on each axis for each shape

total_bounds: tuple of max and min coordinates on each axis for entire GeoSeries

boundary:a lower dimensional object representing the object’s set-theoretic boundary. (The boundary of a polygon is a line, the boundary of a line is a collection of points. The boundary of a point is an empty collection.)

geom_type: type of geometry.

is_valid: tests if coordinates make a shape that is reasonable geometric shape according to the Simple Feature Access standard.

nybb_gdf.columns

Index(['BoroCode', 'BoroName', 'Shape_Leng', 'Shape_Area', 'geometry'], dtype='object')

nybb_gdf["boundary"] = nybb_gdf.boundary
nybb_gdf["boundary"]

  MULTILINESTRING ((970217.022 145643.332, 97022...
  MULTILINESTRING ((1029606.077 156073.814, 1029...
  MULTILINESTRING ((1021176.479 151374.797, 1021...
  MULTILINESTRING ((981219.056 188655.316, 98094...
  MULTILINESTRING ((1012821.806 229228.265, 1012...
Name: boundary, dtype: geometry

nybb_gdf["centroid"] = nybb_gdf.centroid
nybb_gdf["centroid"]

   POINT (941639.450 150931.991)
  POINT (1034578.078 197116.604)
   POINT (998769.115 174169.761)
   POINT (993336.965 222451.437)
  POINT (1021174.790 249937.980)
Name: centroid, dtype: geometry

With columnes boundary and centroid saved to the nybb_gdf, we now have three geometry comulmns in the same GeoDataFrame

14.7. Applying Basic Methods#

distance(): returns Series with minimum distance from each entry to other

centroid: returns GeoSeries of centroids

representative_point(): returns GeoSeries of points that are guaranteed to be within each geometry. It does NOT return centroids.

to_crs(): change coordinate reference system. See projections

plot(): plot GeoSeries. See mapping.

We can measure the distance of each centroid from the first centroid location:

first_point = nybb_gdf["centroid"].iloc[0]
nybb_gdf["distance"] = nybb_gdf["centroid"].distance(first_point)
nybb_gdf["distance"]

       0.000000
  103781.535276
   61674.893421
   88247.742789
  126996.283623
Name: distance, dtype: float64

14.8. Plotting GeoPandas#

Similar to pandas, if you call the plot() method of GeoDataFrame it will plot the active geometry column of it.

nybb_gdf.plot()

<Axes: >

../_images/cee2fdf0b19b35abc936c3a7a81173dcfca211ff43b51eb4ced6f4f052e70d64.png

You can customize this and plot specific attributes of the GeoDataFrame. In the following, we will calculate the area of each object in the GeoDataFrame and plot it.

nybb_gdf["area"] = nybb_gdf.area

nybb_gdf.plot("area", legend=True)

<Axes: >

../_images/a37880c35ae28f1e87f4502b97e2ed7a5795377e48c7cd4bda9e9abb3743a183.png

There is also an interactive way to plot the data using the explore function. This uses folium/leaflet.js to plot the data.

nybb_gdf.explore("area", legend=False)

Make this Notebook Trusted to load map: File -> Trust Notebook

Now, let’s calculate the centroid of each borough, and set that as the active geometry of the GeoDataFrame

nybb_gdf["centroid"] = nybb_gdf.centroid

nybb_gdf.set_geometry("centroid", inplace = True)

nybb_gdf.plot()

<Axes: >

../_images/eee59755bb691d8f2d651483b79c347c2047ddcc768ac64d7db1e27a1d3784bb.png

As you can see, the plot() method in this case plots the centroids which is the active geometry of the GeoDataFrame.

nybb_gdf.columns

Index(['BoroCode', 'BoroName', 'Shape_Leng', 'Shape_Area', 'geometry',
       'boundary', 'centroid', 'distance', 'area'],
      dtype='object')

nybb_gdf.geometry

   POINT (941639.450 150931.991)
  POINT (1034578.078 197116.604)
   POINT (998769.115 174169.761)
   POINT (993336.965 222451.437)
  POINT (1021174.790 249937.980)
Name: centroid, dtype: geometry

nybb_gdf["geometry"]

  MULTIPOLYGON (((970217.022 145643.332, 970227....
  MULTIPOLYGON (((1029606.077 156073.814, 102957...
  MULTIPOLYGON (((1021176.479 151374.797, 102100...
  MULTIPOLYGON (((981219.056 188655.316, 980940....
  MULTIPOLYGON (((1012821.806 229228.265, 101278...
Name: geometry, dtype: geometry

Note: By default when you use the read_file() command, the column containing spatial objects from the file is named “geometry” by default, and will be set as the active geometry column. However, despite using the same term for the name of the column and the name of the special attribute that keeps track of the active column, they are distinct. You can easily shift the active geometry column to a different GeoSeries with the set_geometry() command. Further, gdf.geometry will always return the active geometry column, not the column named geometry. If you wish to call a column named “geometry”, and a different column is the active geometry column, use gdf['geometry'], not gdf.geometry.

# Let's set the geometry back to its original column
nybb_gdf.set_geometry("geometry", inplace = True)

# Plot centroids on the same map as polygons
ax = nybb_gdf.plot()
nybb_gdf["centroid"].plot(ax = ax, color = "red")

<Axes: >

../_images/3b8b0f780732f8025bcbbafe73582442eb7a17c2035db1d8e0a2cca775a2f1d6.png

nybb_gdf.crs

<Projected CRS: EPSG:2263>
Name: NAD83 / New York Long Island (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - New York - counties of Bronx; Kings; Nassau; New York; Queens; Richmond; Suffolk.
- bounds: (-74.26, 40.47, -71.8, 41.3)
Coordinate Operation:
- name: SPCS83 New York Long Island zone (US Survey feet)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

14.9. Plotting with CartoPy and GeoPandas#

Let’s load another dataset, and plot it with both GeoPandas and CartoPy. Our goal is to use the features of CartoPy to create a better visualization of the data in a GeoDataFrame.

We are going to also transform the projection of a GeoDataFrame to a different crs.

land_gdf = gpd.read_file(geodatasets.get_url("naturalearth.land"))

land_gdf.plot()

<Axes: >

../_images/7e940f2203d61f6d4674ec05a203c18b1e4d86c2dcdbe3a69a255153691eb2fb.png

We can retrieve the CRS of the GeoDataFrame:

land_gdf.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

We are going to a define a new CRS using CartoPy:

crs = ccrs.AzimuthalEquidistant()

And convert it to a proj4 string/dict compatible with GeoPandas

crs_proj4 = crs.to_proj4()

/Users/hamed/miniconda3/envs/vector-env/lib/python3.11/site-packages/pyproj/crs/crs.py:1293: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  proj = self._crs.to_proj4(version=version)

Now, we project the land_gdf GeoDataFrame to the new crs

land_gdf_ae = land_gdf.to_crs(crs_proj4)

land_gdf_ae.plot()

<Axes: >

../_images/47d01df091e7bff2f24da1a9931b7a41b76e25ede8306c0b14aba4a0289169f4.png

We could have done this with CartoPy as well now that our data is in a crs understandable to CartoPy

fig, ax = plt.subplots(subplot_kw={"projection": crs})
ax.add_geometries(land_gdf_ae["geometry"], crs=crs)

<cartopy.mpl.feature_artist.FeatureArtist at 0x157bf4350>

../_images/fd60486e9e43612578d505d01517ef2d0f6af8c77d3e82ac0503ed217b53ba07.png

14.10. Merging Data#

We can merge different GeoPandas datasets using attribute joins or spatial joins.

In attribute join, a GeoSeries or GeoDataFrame is merged with a regular pandas.Series or pandas.DataFrame based on common variables (this is similar to merging in pandas).

In a spatial join, observations from two GeoSeries or GeoDataFrame are combined based on their spatial relationships to one another.

Let’s import two example datasets to work with:

chicago = gpd.read_file(geodatasets.get_path("geoda.chicago_commpop"))

groceries = gpd.read_file(geodatasets.get_path("geoda.groceries"))

chicago_shapes = chicago[['geometry', 'NID']]

chicago_names = chicago[['community', 'NID']]

chicago = chicago[['geometry', 'community']].to_crs(groceries.crs)

chicago.head()

	geometry	community
0	MULTIPOLYGON (((1181573.250 1886828.039, 11815...	DOUGLAS
1	MULTIPOLYGON (((1186289.356 1876750.733, 11862...	OAKLAND
2	MULTIPOLYGON (((1176344.998 1871187.546, 11763...	FULLER PARK
3	MULTIPOLYGON (((1182322.043 1876674.730, 11823...	GRAND BOULEVARD
4	MULTIPOLYGON (((1186289.356 1876750.733, 11862...	KENWOOD

groceries.head()

	OBJECTID	Ycoord	Xcoord	Status	Address	Chain	Category	geometry
0	16	41.973266	-87.657073	OPEN	1051 W ARGYLE ST, CHICAGO, IL. 60640	VIET HOA PLAZA	NaN	MULTIPOINT (1168268.672 1933554.350)
1	18	41.696367	-87.681315	OPEN	10800 S WESTERN AVE, CHICAGO, IL. 60643-3226	COUNTY FAIR FOODS	NaN	MULTIPOINT (1162302.618 1832900.224)
2	22	41.868634	-87.638638	OPEN	1101 S CANAL ST, CHICAGO, IL. 60607-4932	WHOLE FOODS MARKET	NaN	MULTIPOINT (1173317.042 1895425.426)
3	23	41.877590	-87.654953	OPEN	1101 W JACKSON BLVD, CHICAGO, IL. 60607-2905	TARGET/SUPER	new	MULTIPOINT (1168996.475 1898801.406)
4	27	41.737696	-87.625795	OPEN	112 W 87TH ST, CHICAGO, IL. 60620-1318	FOOD 4 LESS	NaN	MULTIPOINT (1176991.989 1847262.423)

14.10.1. Appending#

Appending GeoDataFrame and GeoSeries uses pandas append() methods. Keep in mind, that appended geometry columns needs to have the same CRS.

joined_geometries = pd.concat([chicago.geometry, groceries.geometry])

joined_geometries

    MULTIPOLYGON (((1181573.250 1886828.039, 11815...
    MULTIPOLYGON (((1186289.356 1876750.733, 11862...
    MULTIPOLYGON (((1176344.998 1871187.546, 11763...
    MULTIPOLYGON (((1182322.043 1876674.730, 11823...
    MULTIPOLYGON (((1186289.356 1876750.733, 11862...
                             ...                        
               MULTIPOINT (1171065.063 1899839.376)
               MULTIPOINT (1165217.798 1914159.975)
               MULTIPOINT (1166186.713 1883581.309)
               MULTIPOINT (1175778.816 1892214.445)
               MULTIPOINT (1185013.734 1832012.356)
Name: geometry, Length: 225, dtype: geometry

joined_geometries.iloc[9]

../_images/824736708fb8dce943837b66b91d8494f696fe6c14b4f68b40ae22510eb2f38e.svg

douglas = chicago[chicago.community == 'DOUGLAS']

oakland = chicago[chicago.community == 'OAKLAND']

douglas_oakland = pd.concat([douglas, oakland])

douglas_oakland

	geometry	community
0	MULTIPOLYGON (((1181573.250 1886828.039, 11815...	DOUGLAS
1	MULTIPOLYGON (((1186289.356 1876750.733, 11862...	OAKLAND

14.10.2. Attribute Joins#

Attribute joins are accomplished using the merge() method. We will use chicago_shapes and chicago_names which both have the NID attribute:

chicago_shapes.head()

	geometry	NID
0	MULTIPOLYGON (((-87.60914 41.84469, -87.60915 ...	35
1	MULTIPOLYGON (((-87.59215 41.81693, -87.59231 ...	36
2	MULTIPOLYGON (((-87.62880 41.80189, -87.62879 ...	37
3	MULTIPOLYGON (((-87.60671 41.81681, -87.60670 ...	38
4	MULTIPOLYGON (((-87.59215 41.81693, -87.59215 ...	39

chicago_names.head()

	community	NID
0	DOUGLAS	35
1	OAKLAND	36
2	FULLER PARK	37
3	GRAND BOULEVARD	38
4	KENWOOD	39

chicago_shapes = chicago_shapes.merge(chicago_names, on='NID')

chicago_shapes.head()

	geometry	NID	community
0	MULTIPOLYGON (((-87.60914 41.84469, -87.60915 ...	35	DOUGLAS
1	MULTIPOLYGON (((-87.59215 41.81693, -87.59231 ...	36	OAKLAND
2	MULTIPOLYGON (((-87.62880 41.80189, -87.62879 ...	37	FULLER PARK
3	MULTIPOLYGON (((-87.60671 41.81681, -87.60670 ...	38	GRAND BOULEVARD
4	MULTIPOLYGON (((-87.59215 41.81693, -87.59215 ...	39	KENWOOD

type(chicago_shapes)

geopandas.geodataframe.GeoDataFrame

chicago_shapes.plot("NID")

<Axes: >

../_images/2a76f2a64f3b7fa6efea13d213f1bf5dce4359ff7a1061c1dcca79125a36f491.png

14.10.3. Spatial Joins#

In a spatial join, two geometry objects are merged based on their spatial relationship to one another.

GeoPandas provides two spatial-join functions:

GeoDataFrame.sjoin(): joins based on binary predicates (intersects, contains, etc.)
GeoDataFrame.sjoin_nearest(): joins based on proximity, with the ability to set a maximum search radius.

14.10.3.1. Binary predicate joins#

Binary predicate joins are available via GeoDataFrame.sjoin() which has two core arguments: how and predicate.

predicate

The predicate argument specifies how GeoPandas decides whether or not to join the attributes of one object to another, based on their geometric relationship. The default spatial index in GeoPandas currently supports the following values for predicate which are defined in the Shapely documentation:

intersects
contains
within
touches
crosses
overlaps

The following figures from the amazing pygis.io webiste provide an intuitive illustration of different predicate methods:

how

The how argument specifies the type of join that will occur and which geometry is retained in the resultant GeoDataFrame. It accepts the following options:

left: All features from the first or left GeoDataFrame are kept, regardless if the features meet the specified spatial relationship criteria for a join. As all attribute fields are combined, rows that do not have a match with the right dataset may have null values in the fields that originated from the right GeoDataFrame.
right: All features from the second or right GeoDataFrame are kept, regardless if the features meet the specified spatial relationship criteria for a join. As all attribute fields are combined, rows that do not have a match with the left dataset may have null values in the fields that originated from the left GeoDataFrame.
inner: Only features from both datasets that meet the spatial relationship for the joined are kept. The geometries from the first or left GeoDataFrame are used for the join.

The following figure from pygis.io shows how these three join options operate (note that the “Outer” join is not implemented in GeoPandas:

Let’s try this using the Chicago and Groceries GeoDataFrames:

chicago.head()

	geometry	community
0	MULTIPOLYGON (((1181573.250 1886828.039, 11815...	DOUGLAS
1	MULTIPOLYGON (((1186289.356 1876750.733, 11862...	OAKLAND
2	MULTIPOLYGON (((1176344.998 1871187.546, 11763...	FULLER PARK
3	MULTIPOLYGON (((1182322.043 1876674.730, 11823...	GRAND BOULEVARD
4	MULTIPOLYGON (((1186289.356 1876750.733, 11862...	KENWOOD

groceries.head()

	OBJECTID	Ycoord	Xcoord	Status	Address	Chain	Category	geometry
0	16	41.973266	-87.657073	OPEN	1051 W ARGYLE ST, CHICAGO, IL. 60640	VIET HOA PLAZA	NaN	MULTIPOINT (1168268.672 1933554.350)
1	18	41.696367	-87.681315	OPEN	10800 S WESTERN AVE, CHICAGO, IL. 60643-3226	COUNTY FAIR FOODS	NaN	MULTIPOINT (1162302.618 1832900.224)
2	22	41.868634	-87.638638	OPEN	1101 S CANAL ST, CHICAGO, IL. 60607-4932	WHOLE FOODS MARKET	NaN	MULTIPOINT (1173317.042 1895425.426)
3	23	41.877590	-87.654953	OPEN	1101 W JACKSON BLVD, CHICAGO, IL. 60607-2905	TARGET/SUPER	new	MULTIPOINT (1168996.475 1898801.406)
4	27	41.737696	-87.625795	OPEN	112 W 87TH ST, CHICAGO, IL. 60620-1318	FOOD 4 LESS	NaN	MULTIPOINT (1176991.989 1847262.423)

Let’s first plot the data:

ax = chicago.plot()
groceries.plot(ax = ax, color="red")

<Axes: >

../_images/c97ca2b39020b723a40dc9a5eb1b21f145a88b6b3c4e67aa734c65b5d9a25e3a.png

Now, we will use the groceries GeoDataFrame to find communities from the chicago GeoDataFrame that intersect with the geometries of each grocery store:

groceries_with_community = groceries.sjoin(chicago, how="inner", predicate='intersects')
groceries_with_community.head()

	OBJECTID	Ycoord	Xcoord	Status	Address	Chain	Category	geometry	index_right	community
0	16	41.973266	-87.657073	OPEN	1051 W ARGYLE ST, CHICAGO, IL. 60640	VIET HOA PLAZA	NaN	MULTIPOINT (1168268.672 1933554.350)	30	UPTOWN
87	365	41.961707	-87.654058	OPEN	4355 N SHERIDAN RD, CHICAGO, IL. 60613-1497	JEWEL OSCO	NaN	MULTIPOINT (1168837.980 1929246.962)	30	UPTOWN
90	373	41.963131	-87.656352	OPEN	4466 N BROADWAY ST, CHICAGO, IL. 60640-5660	TARGET	NaN	MULTIPOINT (1168471.227 1929825.061)	30	UPTOWN
140	582	41.969131	-87.674882	Chicago-Ravenswood	1800 W Lawrence Ave, Chicago, IL 60640	Mariano's	NaN	MULTIPOINT (1163502.978 1932264.462)	30	UPTOWN
1	18	41.696367	-87.681315	OPEN	10800 S WESTERN AVE, CHICAGO, IL. 60643-3226	COUNTY FAIR FOODS	NaN	MULTIPOINT (1162302.618 1832900.224)	73	MORGAN PARK

Question: If we use left as the value for how in the sjoin() method, how will the result be different?

14.10.3.2. Nearest joins#

Proximity-based joins can be done via GeoDataFrame.sjoin_nearest().

GeoDataFrame.sjoin_nearest() shares the how argument with GeoDataFrame.sjoin(), and includes two additional arguments: max_distance and distance_col.

max_distance

The max_distance argument specifies a maximum search radius for matching geometries. This can have a considerable performance impact in some cases. If you can, it is highly recommended that you use this parameter.

distance_col

If set, the resultant GeoDataFrame will include a column with this name containing the computed distances between an input geometry and the nearest geometry.