<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/groundwater-levels-across-us-affected-by-mexico-earthquake/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://owi.usgs.gov/blog/earthquake/"> The USGS OWI blog </a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h1>What happened?</h1>
<p>A magnitude 8.1 earthquake was measured in the Pacific region offshore Chiapas, Mexico at 04:49:21 UTC on September 8th 2017.</p>
<p>Working in the Water Mission Area for the USGS, our <a rel="nofollow" href="https://owi.usgs.gov/datascience/" target="_blank">Data Science Team</a> was already hard at work polishing up water visualizations for the devastating <a rel="nofollow" href="https://owi.usgs.gov/vizlab/hurricane-harvey/" target="_blank">Hurricane Harvey</a>, as well as creating a new visualization for <a rel="nofollow" href="https://owi.usgs.gov/vizlab/hurricane-irma/" target="_blank">Hurricane Irma</a>. That is to say, the resources available for a timely visualization were already spent.</p>
<p>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 <a rel="nofollow" href="https://waterdata.usgs.gov/nwis" target="_blank">NWIS</a> pages.</p>
<p>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:</p>
<p><a rel="nofollow" href="https://waterdata.usgs.gov/nwis/uv?cb_72019=on&format=gif_default&sites&site_no=402411077374801&site_no=364818094185302&site_no=405215084335400&site_no=370812080261901&site_no=393851077343001&site_no=444302070252401&site_no=324307117063502&site_no=401804074432601&site_no=292618099165901&site_no=421157075535401&site_no=373904118570701&site_no=343457096404501&period=&begin_date=2017-09-07&end_date=2017-09-08" target="_blank">Affected Groundwater Levels</a></p>
<p>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 blog 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.</p>
<h2>Aggregate USGS Data</h2>
<p>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.</p>
<p>Using the <a rel="nofollow" href="https://github.com/USGS-R/dataRetrieval" target="_blank">dataRetrieval</a> R package, it then becomes trivial to access USGS data and import it into the R environment.</p>
<p>The following code retrieves the data:</p>
<pre><code>library(dataRetrieval)

sites &lt;- c(&quot;402411077374801&quot;,
           &quot;364818094185302&quot;,
           &quot;405215084335400&quot;,
           &quot;370812080261901&quot;,
           &quot;393851077343001&quot;,
           &quot;444302070252401&quot;,
           &quot;324307117063502&quot;,
           &quot;421157075535401&quot;,
           &quot;373904118570701&quot;,
           &quot;343457096404501&quot;,
           &quot;401804074432601&quot;,
           &quot;292618099165901&quot;)

gw_data &lt;- readNWISuv(sites, 
                      parameterCd = &quot;72019&quot;, 
                      startDate = &quot;2017-09-07&quot;,
                      endDate = &quot;2017-09-08&quot;)
</code></pre>
<p>By default, <em>all</em> 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.</p>
<p><code>dataRetrieval</code> attaches an attribute <code>siteInfo</code> 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 <code>dplyr</code>.</p>
<pre><code>library(dplyr)

unique_sites &lt;- attr(gw_data, &quot;siteInfo&quot;)

south_to_north &lt;- unique_sites %>%
  arrange(desc(dec_lat_va))
</code></pre>
<p>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 <code>south_to_north</code> data frame to create ordered levels for the sites:</p>
<pre><code>gw_data_sub &lt;- gw_data %>%
            filter(dateTime > as.POSIXct(&quot;2017-09-07 00:00:00&quot;, tz=&quot;UTC&quot;),
                   dateTime &lt; as.POSIXct(&quot;2017-09-09 00:00:00&quot;, tz=&quot;UTC&quot;)) %>%
  mutate(site_no = factor(site_no, 
                          levels = south_to_north$site_no))
</code></pre>
<h2>Plot Explorations</h2>
<p>Since the y-axis is quite different between sites, it seemed to make the most sense to use <code>ggplot2</code>’s faceting. The first iteration was as follows:</p>
<pre><code>library(ggplot2)

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

depth_plot
</code></pre>
<p><img alt="Initial water level plot" title="Initial plot of water levels" src="https://i1.wp.com/owi.usgs.gov/blog/static/earthquake/first_try-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The <code>site_no</code> 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.</p>
<p>So, first, let’s map some colors to the unique <code>site_no</code>:</p>
<pre><code>cbValues &lt;- c(&quot;#DCDA4B&quot;,&quot;#999999&quot;,&quot;#00FFFF&quot;,&quot;#CEA226&quot;,&quot;#CC79A7&quot;,&quot;#4E26CE&quot;,
              &quot;#0000ff&quot;,&quot;#78C15A&quot;,&quot;#79AEAE&quot;,&quot;#FF0000&quot;,&quot;#00FF00&quot;,&quot;#B1611D&quot;,
              &quot;#FFA500&quot;,&quot;#F4426e&quot;, &quot;#800000&quot;, &quot;#808000&quot;)
cbValues &lt;- cbValues[1:length(sites)]
names(cbValues) &lt;-   sites
</code></pre>
<h2>Map Legend</h2>
<p>The following code was used to map the sites to a basic lower-48 map of the United States:</p>
<pre><code>library(maptools)
library(maps)
library(sp)

proj.string &lt;- &quot;+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs&quot;

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

conus &lt;- to_sp('state')

wgs84 &lt;- &quot;+init=epsg:4326&quot;

sites_names &lt;- select(unique_sites, dec_lon_va, dec_lat_va, site_no)

coordinates(sites_names) &lt;- c(1,2)
proj4string(sites_names) &lt;- CRS(wgs84)

sites_names = sites_names %>% 
  spTransform(CRS(proj4string(conus)))
sites_names.df &lt;- as.data.frame(sites_names)

gsMap &lt;- ggplot() +
  geom_polygon(aes(x = long, y = lat, group = group),
               data = conus, fill = &quot;grey90&quot;,
               color = &quot;white&quot;)+
  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(),
        legend.position=&quot;none&quot;)

gsMap
</code></pre>
<p><img alt="Initial map work" title="Site locations on US map" src="https://i1.wp.com/owi.usgs.gov/blog/static/earthquake/map-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The initial line graphs needed to be updated with those colors and dropping the facet labels:</p>
<pre><code>depth_plot &lt;- 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 = &quot;reverse&quot;) +
  facet_grid(site_no ~ ., scales = &quot;free&quot;) +
  ylab(label = attr(gw_data, &quot;variableInfo&quot;)$variableName) +
  xlab(label = &quot;UTC time&quot;) +
  scale_color_manual(values = cbValues) +
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        legend.position=&quot;none&quot;,
        panel.spacing.y=unit(0.04, &quot;lines&quot;))
</code></pre>
<p>Finally, the graphs were combined using the <code>viewport</code> function from the <code>grid</code> package:</p>
<pre><code>library(grid)
vp &lt;- viewport(width = 0.35, height = 0.22, 
               x = 0.26, y = 0.87)

png(&quot;static/earthquake/earthquake.png&quot;,width = 8, height = 8, units = &quot;in&quot;, res = 200)
print(depth_plot)
print(gsMap, vp = vp)
dev.off()
</code></pre>
<pre><code>## png 
##   2
</code></pre>
<p><img alt="Water levels in US affected by Mexico earthquake" title="Water levels in US affected by Mexico earthquake" src="https://i1.wp.com/owi.usgs.gov/blog/static/earthquake/earthquake.png?w=456&ssl=1" data-recalc-dims="1"></p>
<h2>What does the plot mean?</h2>
<p>While this blog 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 <a rel="nofollow" href="https://www.facebook.com/USGSNaturalHazards/posts/136010120350602" target="_blank">USGSNaturalHazards</a>:</p>
<p>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: <a rel="nofollow" href="https://groundwaterwatch.usgs.gov/" target="_blank">https://groundwaterwatch.usgs.gov/</a> 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.</p>
<p>Read about how the 2011 earthquake in Mineral, VA affected ground water levels here: <a rel="nofollow" href="https://water.usgs.gov/ogw/eq/VAquake2011.html" target="_blank">https://water.usgs.gov/ogw/eq/VAquake2011.html</a></p>
<p>Find more information about groundwater response to earthquakes. <a rel="nofollow" href="https://earthquake.usgs.gov/learn/topics/groundwater.php" target="_blank">https://earthquake.usgs.gov/learn/topics/groundwater.php</a></p>
<p>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.</p>
<h2>Have Questions on Groundwater?</h2>
<ul>
<li>Rodney A. Sheets <a rel="nofollow" href="mailto:rasheets@usgs.gov" target="_blank"><i></i></a> <a rel="nofollow" href="https://www.usgs.gov/staff-profiles/rodney-a-sheets" target="_blank"><i></i></a></li>
<li>Charles Schalk <a rel="nofollow" href="mailto:cwschalk@usgs.gov" target="_blank"><i></i></a> <a rel="nofollow" href="https://www.usgs.gov/staff-profiles/charles-schalk" target="_blank"><i></i></a></li>
<li>Leonard Orzol <a rel="nofollow" href="mailto:llorzol@usgs.gov" target="_blank"><i></i></a> <a rel="nofollow" href="https://www.usgs.gov/staff-profiles/leonard-orzol" target="_blank"><i></i></a></li>
</ul>
<h2>Have Questions on NWIS?</h2>
<ul>
<li>Brad Garner <a rel="nofollow" href="mailto:bdgarner@usgs.gov" target="_blank"><i></i></a> <a rel="nofollow" href="https://www.usgs.gov/staff-profiles/bradley-garner" target="_blank"><i></i></a></li>
</ul>
<h2>Have Questions on the R package: dataRetrieval?</h2>
<table>
<colgroup>
<col width="19%">
<col width="73%">
</colgroup>
<tbody>
<tr>
<td><a rel="nofollow" href="https://github.com/USGS-R/dataRetrieval" target="_blank"><img alt="USGS-R image icon" src="https://i1.wp.com/owi.usgs.gov/blog/images/USGS_R.png?w=456&ssl=1" data-recalc-dims="1"></a></td>
<td><a rel="nofollow" href="https://github.com/USGS-R/dataRetrieval/issues" target="_blank">dataRetrieval</a>: This R package is designed to obtain USGS or EPA water quality sample data, streamflow data, and metadata directly from web services.</td>
</tr>
</tbody>
</table>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://owi.usgs.gov/blog/earthquake/"> The USGS OWI blog </a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/UdY1ke_05kI" height="1"></a></div>