Dot gov

Official websites use .gov
A .gov website belongs to an official government organization in the United States.


Secure .gov websites use HTTPS
A lock ( ) or https:// means you’ve safely connected to the .gov website. Share sensitive information only on official, secure websites.

Groundwater levels across US affected by Mexico earthquake

Using ggplot2 to inspect water levels affected by earthquake.

Date Posted September 12, 2017 Last Updated July 20, 2023
Author Laura DeCicco
Reading Time 6 minutes Share

What happened?

A magnitude 8.1 earthquake was measured in the Pacific region offshore Chiapas, Mexico at 04:49:21 UTC on September 8th 2017.

Working in the Water Mission Area for the USGS, our Data Science Team was already hard at work polishing up water visualizations for the devastating Hurricane Harvey, as well as creating a new visualization for Hurricane Irma. That is to say, the resources available for a timely visualization were already spent.

By 10 am on Friday however, a team of USGS groundwater scientist had realized that the affects of the earthquake were being observed in groundwater water level measurements across the United States, and were collecting a list of sites that clearly showed the phenomena. It was interesting to view the data directly from the National Water Information System NWIS pages.

While the NWIS site allowed a very easy way to observe the data, it was not quite as dramatic and interesting since it was only plotted in local time zones (so the earthquake peaks were not lining up), and the y-axis scales were so different, they could not be plotted on a single graph:

Affected Groundwater Levels

The challenge was on: could we very quickly get a graph out that aggregates this data, unifies the time zones, and visualize the information in an intuitive way for the general population? This post documents the process of rapidly creating the plot in order to get the information out to social media in a timely manner. The final plot could certainly be improved, but with the limited resources and time constraints, it still managed to be a fun and interesting graphic.

Aggregate USGS Data

Without question, the hardest part is to find the right data. Luckily for me, this was already done by 10am by the groundwater scientists Rodney A. Sheets, Charles Schalk, and Leonard Orzol (contact information is provided at the end of this post). They provided me with a list of unique site identification numbers from groundwater wells across the nation, and a “parameter code” that represents water levels measurements.

Using the dataRetrieval R package, it then becomes trivial to access USGS data and import it into the R environment.

The following code retrieves the data:


sites <- c("402411077374801",

gw_data <- readNWISuv(sites,
                      parameterCd = "72019",
                      startDate = "2017-09-07",
                      endDate = "2017-09-08")

By default, all the data in this data frame comes back in the “UTC” time zone. So, we’ve already crossed one of the hurdles - unifying time zones.

dataRetrieval attaches an attribute siteInfo to the returned data. The information in that attribute includes station names and location. We can extract that information and sort the sites from south to north using a bit of dplyr.


unique_sites <- attr(gw_data, "siteInfo")

south_to_north <- unique_sites %>%

Now, we can do a bit of filtering so that we will only plot the 24 hours of on September 8th (UTC). Also, we can use the south_to_north data frame to create ordered levels for the sites:

gw_data_sub <- gw_data %>%
            filter(dateTime > as.POSIXct("2017-09-07 00:00:00", tz="UTC"),
                   dateTime < as.POSIXct("2017-09-09 00:00:00", tz="UTC")) %>%
  mutate(site_no = factor(site_no,
                          levels = south_to_north$site_no))

Plot Explorations

Since the y-axis is quite different between sites, it seemed to make the most sense to use ggplot2’s faceting. The first iteration was as follows:


depth_plot <- ggplot(data = gw_data_sub) +
  geom_line(aes(x = dateTime, y = X_72019_00000)) +
  theme_bw() +
  scale_y_continuous(trans = "reverse") +
  facet_grid(site_no ~ ., scales = "free") +
  ylab(label = attr(gw_data, "variableInfo")$variableName) +
  xlab(label = "UTC time") +
  theme(strip.background = element_blank(),
        strip.text.y = element_text(angle = 0),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

Initial water level plot

Initial plot of water levels

The site_no column (the unique site identification numbers) are not terribly useful for a general audience, so we decided to eliminate those facet labels on the right and add a map as a legend. Also, we wanted to use consistent colors across the map and lines.

So, first, let’s map some colors to the unique site_no:

cbValues <- c("#DCDA4B","#999999","#00FFFF","#CEA226","#CC79A7","#4E26CE",
              "#FFA500","#F4426e", "#800000", "#808000")
cbValues <- cbValues[1:length(sites)]
names(cbValues) <-   sites

Map Legend

The following code was used to map the sites to a basic lower-48 map of the United States:


proj.string <- "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

to_sp <- function(...){
  map <- maps::map(..., fill=TRUE, plot = FALSE)
  IDs <- sapply(strsplit(map$names, ":"), function(x) x[1])
  map.sp <- map2SpatialPolygons(map, IDs=IDs, proj4string=CRS("+proj=longlat +datum=WGS84"))
  map.sp.t <- spTransform(map.sp, CRS(proj.string))

conus <- to_sp('state')

wgs84 <- "+init=epsg:4326"

sites_names <- select(unique_sites, dec_lon_va, dec_lat_va, site_no)

coordinates(sites_names) <- c(1,2)
proj4string(sites_names) <- CRS(wgs84)

sites_names = sites_names %>%
sites_names.df <-

gsMap <- ggplot() +
  geom_polygon(aes(x = long, y = lat, group = group),
               data = conus, fill = "grey90",
               color = "white")+
  geom_point(data = sites_names.df,
             aes(x = dec_lon_va, y=dec_lat_va, color = site_no),
             size = 2) +
  scale_color_manual(values = cbValues) +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        panel.background = element_rect(fill = 'white',
                                        colour = 'black'),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.border = element_blank(),

Initial map work

Site locations on US map

The initial line graphs needed to be updated with those colors and dropping the facet labels:

depth_plot <- ggplot(data = gw_data_sub) +
  geom_line(aes(x = dateTime, y = X_72019_00000,
                color = site_no), size = 1.5) +
  theme_bw() +
  scale_y_continuous(trans = "reverse") +
  facet_grid(site_no ~ ., scales = "free") +
  ylab(label = attr(gw_data, "variableInfo")$variableName) +
  xlab(label = "UTC time") +
  scale_color_manual(values = cbValues) +
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.spacing.y=unit(0.04, "lines"))

Finally, the graphs were combined using the viewport function from the grid package:

vp <- viewport(width = 0.35, height = 0.22,
               x = 0.26, y = 0.87)

png("static/earthquake/earthquake.png",width = 8, height = 8, units = "in", res = 200)
print(gsMap, vp = vp)
## png
##   2
Water levels in US affected by Mexico earthquake

Water levels in US affected by Mexico earthquake

What does the plot mean?

While this post is primarily focused on the process used to create the plot, far more effort was expended by others in understanding and describing the data. The following information was released with the graph, initially as a Facebook post USGSNaturalHazards:

Did you know? We often see a response to large (and sometimes not so large) earthquakes in groundwater levels in wells. The USGS maintains a network of wells for monitoring various things like natural variability in water levels and response to pumping and climate change across the U.S. The well network can be seen here: The M8.1 earthquake that took place 9-8-2017 in the Pacific region offshore Chiapas, Mexico was observed across the U.S. in confined aquifer wells within the network. The wells in this plot are available from the National Water Information System (NWIS) web are from: TX, CA (southern), OK, MO, VA, CA (central), MD, NJ, PA, OH, NY and ME (from south to north). This is just a small sampling of the wells. USGS hydrologists are poring over the data from other states to see what the full extent of response might have been.

Read about how the 2011 earthquake in Mineral, VA affected ground water levels here:

Find more information about groundwater response to earthquakes.

The frequency of data in these plots varies – some are collected every minute, some every 15 minutes, and some hourly – depending upon the purpose of the monitoring, giving different responses to the earthquake. The magnitude of the response depends on the local geology and hydrology, when the data point was collected, and when the seismic wave hit the aquifer that the well penetrates.

Have Questions on Groundwater?

  • Rodney A. Sheets
  • Charles Schalk
  • Leonard Orzol

Have Questions on NWIS?

  • Brad Garner

Have Questions on the R package: dataRetrieval?

USGS-R image icondataRetrieval: This R package is designed to obtain USGS or EPA water quality sample data, streamflow data, and metadata directly from web services.


Related Posts

  • Calculating Moving Averages and Historical Flow Quantiles

    October 25, 2016

    This post will show simple way to calculate moving averages, calculate historical-flow quantiles, and plot that information. The goal is to reproduce the graph at this link: PA Graph.

  • Hydrologic Analysis Package Available to Users

    July 26, 2022

    A new R computational package was created to aggregate and plot USGS groundwater data, providing users with much of the functionality provided in Groundwater Watch and the Florida Salinity Mapper.

  • Large sample pulls using dataRetrieval

    July 26, 2022

    dataRetrieval is an R package that provides user-friendly access to water data from either the Water Quality Portal (WQP) or the National Water Information Service (NWIS).

  • dataRetrieval Tutorial - Using R to Discover Data

    January 8, 2021

    R is an open-source programming language. It is known for extensive statistical capabilities, and also has powerful graphical capabilities. Another benefit of R is the large and generally helpful user-community.

  • The Hydro Network-Linked Data Index

    November 2, 2020

    Introduction updated 11-2-2020 after updates described here. The Hydro Network-Linked Data Index (NLDI) is a system that can index data to NHDPlus V2 catchments and offers a search service to discover indexed information.