12 Streamline your code

In this lesson we will learn how to extract geospatial data from a CSV to create a geopandas.GeoDataFrame, introduce more customizations for maps and matplotlib figures, and go over strategies to streamline our code.

By the end of this lesson, students will be able to:

  • create a geopandas.GeoDataFrame from a CSV with longitude and latitude columns
  • format printed output using f-strings
  • use for loops, alone and with zip(), to avoid repeating code
  • check that their code does what they expect using assert
  • create figures with multiple axes

About the data

The U.S. energy landscape relies on a mix of fossil fuels and renewables, each with unique environmental and economic impacts. As the nation works toward sustainability and energy security, understanding this energy mix is essential for informed policy and progress toward cleaner energy.

In this lesson, we will use data from the U.S. Energy Information Administration (EIA) about operable electric generating plants in the United States by energy source, as of February 2025 [1]. The dataset includes information on plant types and energy sources, offering insights into the diversity of power sources—from fossil fuels to renewables—that supply electricity nationwide. The dataset’s metadata can be accessed here. The EIA data on electric plants has been downloaded as a CSV and reprojected into the EPSG:4269 CRS for this lesson. It can be accessed here.

Each row in the data is one power plant. These are the columns we will use in this lesson, with descriptions adapted from the EIA metadata:

Column Description
state State where the power plant is located. Also includes the District of Columbia and Puerto Rico
primsource Primary energy source of the power plant, for example solar, natural gas or wind. It is the source with the largest net summer capacity in the plant
install_mw Installed (nameplate) capacity, in megawatts (MW): the maximum output of the plant’s generators as rated by the manufacturer
longitude, latitude Location of the power plant, in the NAD83 CRS (EPSG:4269)

Additionally, we will use a TIGER shapefile of the US states from the United States Census Bureau [2]. TIGER stands for Topologically Integrated Geographic Encoding and Referencing. This used to be the data format in which the US Census distributed geospatial data, but since 2008 TIGER files are converted to shapefiles. You can view the metadata for all the TIGER shapefiles here.

Follow these steps to download the shapefile with the United States’ states:

  1. At the 2026 TIGER page, under Download, click on “Web Interface”
  2. For year, select 2026, and for layer type select “States (and equivalent)”. Click submit.
  3. Click on “Download national file”.

The column descriptions for the US states shapefile, from the TIGER/Line Shapefiles technical documentation, are:

Column Type Description
REGION String Region code
DIVISION String Division code
STATEFP String State FIPS code
STATENS String State NS code
GEOID String State identifier; state FIPS code
GEOIDFQ String Fully qualified GEOID; used to join data tables from data.census.gov
STUSPS String Two-character United States Postal Service (USPS) state postal abbreviation
NAME String State name
LSAD String LSAD code for the state
MTFCC String MAF/TIGER Feature Class Code (G4000)
FUNCSTAT String Functional status
ALAND Number Land area
AWATER Number Water area
INTPTLAT String Latitude of the internal point
INTPTLON String Longitude of the internal point

CSV to geopandas.GeoDataFrame

Let’s start by importing packages and updating viewing options:

import os

import pandas as pd
from pandas.api.types import is_string_dtype, is_numeric_dtype
import geopandas as gpd
import matplotlib.pyplot as plt


# Display all columns when looking at dataframes
pd.set_option("display.max_columns", None)

Next, we import the power plants dataset. In this lesson, we read it directly from its URL.

# Import power plants data
URL = 'https://raw.githubusercontent.com/carmengg/eds-220-book/refs/heads/main/data/power_plants_epsg4269.csv'
power_plants = pd.read_csv(URL)

power_plants.head(3)
objectid plant_code plant_name utility_id utility_name sector_name street_address city county state zip primsource source_desc tech_desc install_mw total_mw bat_mw bio_mw coal_mw geo_mw hydro_mw hydrops_mw ng_mw nuclear_mw crude_mw solar_mw wind_mw other_mw source period longitude latitude
0 62897 1873 Tower 20910 Wolverine Power Supply Coop Electric Utility 9663 Cooperative Road Tower Cheboygan Michigan 49792 petroleum Petroleum = 19.3 MW Petroleum Liquids 21.3 19.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 19.3 0.0 0.0 0.0 EIA-860, EIA-860M and EIA-923 202502 -84.294970 45.362818
1 62898 1874 Caro 18895 Thumb Electric Coop of Mich Electric Utility 2437 East Dayton Road Caro Tuscola Michigan 48723 petroleum Petroleum = 8.5 MW Petroleum Liquids 9.4 8.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 8.5 0.0 0.0 0.0 EIA-860, EIA-860M and EIA-923 202502 -83.327806 43.465899
2 62899 1875 Ubly 18895 Thumb Electric Coop of Mich Electric Utility 4159 Ubly Road Ubly Huron Michigan 48475 petroleum Petroleum = 12.3 MW Petroleum Liquids 12.4 12.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 12.3 0.0 0.0 0.0 EIA-860, EIA-860M and EIA-923 202502 -82.938806 43.721699

The power plants file is a CSV. Unlike shapefiles or other geospatial file formats, a CSV doesn’t store a geometry column or a CRS, so geopandas cannot reliably build them on its own. We will need to create the geometry and set the CRS ourselves.

To do so we will use the longitude and latitude columns in the CSV. These indicate the location of the power plants in the NAD83 CRS (EPSG:4269). We can use this information to create a new geopandas.GeoDataFrame from the pandas.DataFrame using the geopandas function points_from_xy():

# Create points from latitude and longitude
points = gpd.points_from_xy(power_plants['longitude'], 
                            power_plants['latitude'])

# Create geodataframe
power_plants = gpd.GeoDataFrame(power_plants,    # Data
                                geometry=points, # Specify geometry column
                                crs='EPSG:4269'  # Specify CRS
                                )

Let’s check that we now have a geometry column:

power_plants.geometry
0         POINT (-84.29497 45.36282)
1         POINT (-83.32781 43.46590)
2         POINT (-82.93881 43.72170)
3         POINT (-85.42600 43.83930)
4         POINT (-85.85510 42.73100)
                    ...             
13440    POINT (-118.10900 34.60800)
13441     POINT (-79.28172 42.50275)
13442     POINT (-77.74272 43.18914)
13443     POINT (-75.90345 42.33782)
13444    POINT (-120.63443 37.08391)
Name: geometry, Length: 13445, dtype: geometry

With the geometry column and CRS, we can plot our dataset:

power_plants.plot()

f-strings

So far, we have printed variables by passing several arguments to the print() function. This means that we write commas between every string and variable we want to print, and then the print() function joins these into a single string, with a space between each one. For example:

print('CRS:', power_plants.crs)
CRS: EPSG:4269

Another popular way of mixing strings and variables in print statements is by creating an f-string which stands for “formatted string”. The simplest syntax for an f-string is:

f" some text {replace}"

where replace can be a variable, an expression, or a function or method call. For example:

# Explore CRS
print(f"ellipsoid: {power_plants.crs.ellipsoid}")
print(f"datum: {power_plants.crs.datum}")
ellipsoid: GRS 1980
datum: North American Datum 1983

We just created a string replacing the value inside the curly brackets {}.

One of the advantages of using f-strings is that they offer customization for formatting the output:

# Any expression can go inside the curly brackets
# :, adds thousands separators
print(f"The data has {len(power_plants):,} power plants in {power_plants['state'].nunique()} states and territories.")

# :.2f rounds to 2 decimals
print(f"Average installed capacity: {power_plants['install_mw'].mean():.2f} MW")
The data has 13,445 power plants in 52 states and territories.
Average installed capacity: 99.39 MW
Use f-strings or not?

Whether you use an f-string or simply concatenate strings with variables inside your print statements depends entirely on the application. For quickly checking a variable, a print statement might be enough, while using f-strings can be better to include custom messages during runtime. The best tool can be different depending on the task!

These are some good resources to learn more about f-string formatting:

Import shapefile

Let’s import the TIGER shapefile

# Import states data
fp = os.path.join('data','tl_2026_us_state','tl_2026_us_state.shp')
states = gpd.read_file(fp)

# Simplify column names 
states.columns = states.columns.str.lower()

states.head(3)
region division statefp statens geoid geoidfq stusps name lsad mtfcc funcstat aland awater intptlat intptlon geometry
0 1 2 36 01779796 36 0400000US36 NY New York 00 G4000 A 122049230845 19256680224 +42.9133974 -075.5962723 MULTIPOLYGON (((-74.04075 40.70017, -74.04073 ...
1 4 8 32 01779793 32 0400000US32 NV Nevada 00 G4000 A 284537080522 1839846018 +39.3310928 -116.6151469 POLYGON ((-119.32418 41.99392, -119.32362 41.9...
2 4 9 02 01785533 02 0400000US02 AK Alaska 00 G4000 A 1479942808337 244276695297 +63.3469818 -152.8385657 MULTIPOLYGON (((-167.55823 60.22436, -167.5567...

and obtain some preliminary geospatial information about the states geodataframe:

print(states.crs)
states.plot()
EPSG:4269

for loops

The region, division, and statefp columns contain codes. The US Census Bureau groups the states into four regions (for example, the West), which are further split into nine divisions (for example, the Pacific division). Each state also has a FIPS (Federal Information Processing Standards) code, a numeric code that identifies it. For example, California’s FIPS code is 06. We will use these codes to select the states we need.

It can be easier to work with the codes as numbers instead of strings, so let’s update the corresponding columns in the states geodataframe. We start by checking the data type of the region, division, and statefp columns:

code_cols = ['region', 'division', 'statefp']

# Check whether codes columns are strings
for column in code_cols: 
    print(f"{column} is string dtype? {is_string_dtype(states[column])}")
region is string dtype? True
division is string dtype? True
statefp is string dtype? True

Remember for loops execute a block of code a fixed number of times, iterating over a set of objects. In this case, we iterate over the list of column names code_cols = ['region', 'division', 'statefp'].

DRY code

We could have checked whether all the region, division, and statefp columns were of string data type by using the following code:

print(f"region is string dtype? {is_string_dtype(states['region'])}")
print(f"division is string dtype? {is_string_dtype(states['division'])}")
print(f"statefp is string dtype? {is_string_dtype(states['statefp'])}")

However, this is inconvenient as it repeats the same pieces of code, only changing the column name. Instead, using the for loop allows us to succinctly print the same information:

code_cols = ['region', 'division', 'statefp']

for column in code_cols: 
    print(f"{column} is string dtype? {is_string_dtype(states[column])}")

Don’t Repeat Yourself (DRY) is a core programming principle that encourages reducing redundancy and consolidating repeated logic. Try implementing it as much as possible! If you need to repeat the “same” code more than twice, you likely need a for loop.

assert

Next, we update the data type of the code columns to be integers. This time, we check the data type of the column using the is_numeric_dtype() function inside an assert statement:

# Update code columns into integers
for column in code_cols:
    states[column] = states[column].astype('int')
    assert is_numeric_dtype(states[column])  # Check conversion

The assert keyword does nothing if the expression next to it evaluates to True. If it evaluates to False, it raises an AssertionError exception and stops your code from running any further. For example,

# Does nothing if statement is True
assert 2+2 == 4

# Raises an error if statement is False
assert 2+2 == 3
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[13], line 5
      2 assert 2+2 == 4
      4 # Raises an error if statement is False
----> 5 assert 2+2 == 3

AssertionError:

In our data type conversion code, since no AssertionError was raised, we can be confident that the data type was updated.

Data selection

For this lesson, we want to use only the contiguous states. As seen in the plot, the data covers a bigger area.

Check-in

From the TIGER shapefiles metadata we know that:

In addition to the fifty states, the Census Bureau treats the District of Columbia, Puerto Rico, and the Island Areas (American Samoa, the Commonwealth of the Northern Mariana Islands, Guam, and the U.S. Virgin Islands) as statistical equivalents of states for the purpose of data presentation.

In this US Census Bureau file we can see what each code for the region, division, and state corresponds to.

  • What are the unique values for region, division, or state codes in the data?
  • Which codes should we select to keep only states in the contiguous US?

Let’s go ahead and select the data:

# Select contiguous US states
contiguous = states[(states['region']!=9) & (~states['statefp'].isin([2,15]))]

In this code we used the syntax

~df['column'].isin([val1, val2, val3])

The ~ tilde symbol is used in Python to negate a statement. So the previous line could be read as “the values in df’s column which are not in the list [val1, val2, val3].”

Check-in

Select the data in the power_plants data frame for the contiguous US states.

Let’s check our answer:

# Select power plants in the contiguous US states
power_plants = power_plants[~power_plants['state'].isin(['Puerto Rico','Hawaii','Alaska'])]

Plotting

Before we plot our data, let’s make sure they are in the same CRS:

assert contiguous.crs == power_plants.crs
Code 
fig, ax = plt.subplots(figsize=(9, 5)) # Update figure size

# Remove the axis for a cleaner map
ax.axis('off')

# Title for the plot
ax.set_title('Operable electric generating plants in the contiguous United States', 
              fontsize=15)

# Add states
contiguous.plot(ax=ax,
               color='none',
               edgecolor='#362312')

# Add electric power plants colored by energy source
power_plants.plot(ax=ax, 
                  column='primsource',
                  legend=True,
                  markersize=4,
                  cmap='tab20',
                  alpha=0.5,
                  legend_kwds={
                      'title': 'Primary energy source',
                      'title_fontsize': 'small',
                      'fontsize': 'small',
                      'loc': 'upper left',
                      'bbox_to_anchor': (0, 0),
                      'ncols': 6  
                  })
                  
plt.show()

In the map above we specified the figure size when creating the plot. This size is given in inches, but can be updated to other units (pixels, cm, etc).

We also controlled the legend location using loc and bbox_to_anchor in the legend_kwds:

  • loc indicates the corner of the legend we want to use for placement, and
  • bbox_to_anchor is a tuple with coordinates indicating where to place the corner specified in loc relative to the axes. Values between 0 and 1 are within the axes.

matplotlib uses a variety of ways to locate elements within the graph and it is best to check the documentation to not spend too much time fiddling with locations.

for with zip

Often, we need to iterate simultaneously over two lists (or other iterables). The zip() function in Python allows you to combine two or more lists (or other iterables) so that you can iterate over their elements in pairs. When used with a for loop, it lets you process elements from each list together, like this example:

# Iterate over a single list
numbers = [1, 2, 3]
for num in numbers:
    print(num)

print('\n')  # Blank line

# Iterate over two lists in pairs using zip()
letters = ['a', 'b', 'c']
for num, letter in zip(numbers, letters):
    print(num, letter)
1
2
3


1 a
2 b
3 c

Let’s see a practical application of for loops and zip() with matplotlib subplots. A common situation when code gets repeated is when creating subplots. For example:

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(7, 3))

axes[0].set_title('This is axis 0')
axes[1].set_title('This is axis 1')
axes[2].set_title('This is axis 2')

plt.show()

In this example, notice that the axes variable returned by the plt.subplots() function is actually an array of axes we can iterate over. Remember that the figure and the axes are separate elements in a matplotlib plot.

Image source: Getting Started with Matplotlib
Check-in

Use for and zip() to create the same subplots and avoid redundancy.

Exercise

  1. Select the power plants in California in a variable named ca_power_plants.

  2. Create a list named top_sources with California’s top 3 primary sources by number of power plants. That is, the 3 values of primsource that appear in the most California power plants.

  3. Isolate the California state boundary in a variable named ca_boundary.

  4. Recreate the following plot:

References

[1]
U.S. Energy Information Administration, “Power Plants in the U.S.” Esri U.S. Federal Datasets, 2025. Available: https://hub.arcgis.com/datasets/fedmaps::power-plants-in-the-u-s--2/about. [Accessed: Oct. 24, 2026]
[2]
U.S. Census Bureau, “TIGER/Line Shapefiles: States (2026).” U.S. Census Bureau, 2026. Available: https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.2026.html. [Accessed: Oct. 24, 2026]