<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/furthest-water/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://owi.usgs.gov/blog/furthest-water/"> The USGS OWI blog </a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Finding the Location Furthest from Water in the Conterminous United States</h2>
<p>The idea for this post came a few months back when I received an email<br>
that started, “I am a writer and teacher and am reaching out to you with<br>
a question related to a piece I would like to write about the place in<br>
the United States that is furthest from a natural body of surface water.<br>
My question is, in short, might you know where this theoretical place is<br>
located?”</p>
<p>As someone who works on the<br>
<a rel="nofollow" href="https://water.usgs.gov/watercensus/" target="_blank">National Water Census</a><br>
and<br>
<a rel="nofollow" href="https://www.usgs.gov/core-science-systems/ngp/national-hydrography" target="_blank">National hydrography data</a>,<br>
this question piqued my interest. I also had an idea of one way to<br>
answer the question that would use some tools and techniques I’ve<br>
recently been developing in R. What follows is just that, <strong>one</strong> way<br>
that this question might be answered. The code I used to find the<br>
solution is presented in detail after the solution to the question.</p>
<p>The solution graphics shown below were generated with the following<br>
code: The <code>furthest_water.R</code> script can be<br>
<a rel="nofollow" href="https://github.com/USGS-OWI/owi-blog/tree/master/content/static/furthest-water" target="_blank">found here</a>.</p>
<pre><code>source(&quot;furthest_water.R&quot;)
furthest_water(scenario = &quot;waterbodies&quot;)

furthest_water(scenario = c(&quot;waterbodies&quot;, &quot;filter_monthly_flow&quot;))

furthest_water(scenario = c(&quot;waterbodies&quot;, &quot;remove_intermittent&quot;, &quot;filter_monthly_flow&quot;))
</code></pre>
<h2>Method</h2>
<p>To answer the question, “where is the furthest place from a natural body<br>
of water?”, I first had to decide where all the places I might want to<br>
check are. To do this, I first built a search box around the<br>
southwestern US to limit the data volume I had to process. I then<br>
created a grid of points 5km apart within the search box and the US<br>
coastline. This set of points, pictured below, will stand in for all the<br>
places we might want to look. Note that there are actually about 150k<br>
points on this image, but they draw as one area because their symbols<br>
overlap.</p>
<p><img alt="All Search
Points" title="All locations searched" src="https://i2.wp.com/owi.usgs.gov/blog/static/furthest-water/all_search_points.png?w=456&ssl=1" data-recalc-dims="1"><br>
With all possible places defined, I needed all the water bodies. The<br>
National Hydrography Dataset contains data used to map rivers and lakes<br>
across the country. As a bonus, there is a dataset known as the “NHDPlus<br>
v2.1 National Seamless” that is available as a single geospatial<br>
database for the whole conterminous US. See<br>
<a rel="nofollow" href="https://www.epa.gov/waterdata/get-data" target="_blank">this EPA page</a><br>
for lots of helpful documentation and data downloads. The NHDPlus<br>
contains a few million river segments and almost a half million water<br>
bodies mapped at 1:100,000-scale. The NHDPlus also has some useful<br>
attributes that categorize waterbodies as ephemeral (sometimes dry) or<br>
perennial (always wet) and modeled estimates of mean flow by month of<br>
the year for rivers. Geospatial data representing rivers and bodies of<br>
water like lakes along with the attributes about how likely they are to<br>
be wet is just what I needed. The maps below show how dense the rivers<br>
and water bodies in the NHDPlus are.</p>
<p><img alt="All flow lines" title="All flow lines" src="https://i1.wp.com/owi.usgs.gov/blog/static/furthest-water/flowlines.png?w=456&ssl=1" data-recalc-dims="1"><br>
<img alt="All water
bodies" title="All water bodies" src="https://i1.wp.com/owi.usgs.gov/blog/static/furthest-water/waterbodies.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The algorithm I used to figure out which of the search locations is<br>
furthest from water uses a nearest neighbor search implemented in an<br>
<a rel="nofollow" href="https://cran.r-project.org/web/packages/RANN/" target="_blank">R package called RANN</a><br>
which allowed me to find which points are within a given distance of a<br>
flowline or water body. By starting with a small distance and<br>
incrementing to larger and larger distances, removing locations as I<br>
went, I was able to narrow down to one point that was further than any<br>
other from all the water bodies. The marvel of this technique is that<br>
it’s not just 150k points and a few million rivers and lakes. To make it<br>
work, I converted the lines and polygons that represent rivers and water<br>
bodies to more than 30 million individual points that make up all the<br>
geometry!</p>
<p>In each of the three scenarios below, there are two graphics. The first,<br>
is a simple visual showing where the location is relative to rivers and<br>
water bodies that are in the vicinity. The second is an animated graphic<br>
showing each step of the nearest neighbor search as the search radius<br>
increases incrementally.</p>
<h2>Furthest from where there might be water</h2>
<p>In the first scenario, I included everything in the NHDPlus that might<br>
be a water body. This includes lakes and rivers that are categorized as<br>
ephemeral and rivers with attributes that show they are usually dry one<br>
or more months a year. As you can see below, the location is in the<br>
Bonneville salt flats. The irony is that, while not categorized as a<br>
water body in the NHD, this picture taken from I-80 shows that the salt<br>
flats do get filled with water periodically!</p>
<p><p>
<br>By<br>
<a rel="nofollow" href="https://commons.wikimedia.org/w/index.php?title=User:Ricraider&action=edit&redlink=1" title="User:Ricraider (page does not exist)" target="_blank">Ricraider</a><br>
– <span lang="en">Own work</span>,<br>
<a rel="nofollow" href="https://creativecommons.org/licenses/by-sa/3.0" title="Creative Commons Attribution-Share Alike 3.0" target="_blank">CC<br>
BY-SA 3.0</a>,<br>
<a rel="nofollow" href="https://commons.wikimedia.org/w/index.php?curid=30988413" target="_blank">Link</a>
</p>
<p><a rel="nofollow" href="https://www.google.com/maps/place/32%C2%B016%2752.9%22N+113%C2%B059&amp;%2339;03.0%22W/@32.2813555,-113.9863557,705m/data=!3m2!1e3!4b1!4m6!3m5!1s0x0:0x0!7e2!8m2!3d32.2813509!4d-113.9841667" target="_blank"><em>Explore this place in google maps.</em></a></p>
<p><img alt="Furthest from where there might be
water" title="Furthest from where there might be water" src="https://i0.wp.com/owi.usgs.gov/blog/static/furthest-water/furthest_water_scenario_1.png?w=456&ssl=1" data-recalc-dims="1"><br>
<img alt="Furthest from where there might be water
animation" title="Furthest from where there might be water animation" src="https://i2.wp.com/owi.usgs.gov/blog/static/furthest-water/animation_scenario_1.gif?w=456&ssl=1" data-recalc-dims="1"></p>
<h2>Furthest from modeled water</h2>
<p>The place we say is the furthest from water shouldn’t be a dry lake bed,<br>
should it? Let’s look at another scenario. In this one, I used the<br>
monthly average flow estimates from the NHDPlus. These estimates are<br>
calculated using the “Enhanced Runoff Method” or EROM. You can find<br>
documentation<br>
<a rel="nofollow" href="https://www.epa.gov/waterdata/learn-more" target="_blank">about the method here.</a><br>
The EROM estimates are available for each month of the year, providing a<br>
normal monthly-mean flow estimate for the period of the analysis (1971<br>
to 2000). For this scenario, I removed any flowline that had a<br>
monthly-mean flow estimate of zero for any month. As can be seen below,<br>
I came up with a location almost on the US-Mexico border in the desert<br>
in far south west Arizona.</p>
<p>This is a place I would believe is actually the furthest from a natural<br>
body of water but I noticed that there were some small waterbodies in<br>
the desert that almost certainly dry out in the summer.</p>
<p><a rel="nofollow" href="https://www.google.com/maps/place/32%C2%B016%2752.9%22N+113%C2%B059&amp;%2339;03.0%22W/@32.2813555,-113.9863557,705m/data=!3m2!1e3!4b1!4m6!3m5!1s0x0:0x0!7e2!8m2!3d32.2813509!4d-113.9841667" target="_blank"><em>Explore this place in google maps</em></a></p>
<p><img alt="Furthest from modeled
water" title="Furthest from modeled water" src="https://i2.wp.com/owi.usgs.gov/blog/static/furthest-water/furthest_water_scenario_2.png?w=456&ssl=1" data-recalc-dims="1"><br>
<img alt="Furthest from modeled water
animation" title="Furthest from modeled water animation" src="https://i0.wp.com/owi.usgs.gov/blog/static/furthest-water/animation_scenario_2.gif?w=456&ssl=1" data-recalc-dims="1"></p>
<h2>Furthest from modeled and non-ephemeral water</h2>
<p>For the third scenario, I used the same mean-monthly flow method to<br>
remove rivers that probably dry out for a month or more a year and I<br>
removed flowlines and water bodies that were categorized as ephemeral.<br>
This removed many water body polygons in the desert that are only wet<br>
part of the year. As can be seen below, now we are in the Sonoran desert<br>
west of Phoenix, AZ. By the looks of things, this really may be the<br>
furthest from a natural body of surface water in the US (over 72km/45mi)<br>
– especially after a long period with no rain.</p>
<p><a rel="nofollow" href="https://www.google.com/maps/place/33%C2%B032%2740.5%22N+113%C2%B041&amp;%2339;36.6%22W/@33.5445945,-113.6956887,695m/data=!3m2!1e3!4b1!4m5!3m4!1s0x0:0x0!8m2!3d33.54459!4d-113.6935" target="_blank"><em>Explore this place in google maps</em></a></p>
<p><img alt="Furthest from modeled and non-ephemeral
water" title="Furthest from modeled and non-ephemeral water" src="https://i2.wp.com/owi.usgs.gov/blog/static/furthest-water/furthest_water_scenario_3.png?w=456&ssl=1" data-recalc-dims="1"><br>
<img alt="Furthest from modeled and non-ephemeral water
animation" title="Furthest from modeled and non-ephemeral water animation" src="https://i2.wp.com/owi.usgs.gov/blog/static/furthest-water/animation_scenario_3.gif?w=456&ssl=1" data-recalc-dims="1"></p>
<p>As the variety of results given different inputs shows, this question<br>
doesn’t necessarily have one answer. It depends what we assume a body of<br>
water is and how often that body of water needs to be wet. The factors<br>
at play that determine if a river or water body is wet are diverse and<br>
highly interrelated. Ranging from the obvious, like how much<br>
<a rel="nofollow" href="https://water.weather.gov/precip/" target="_blank">it’s rained recently</a>,<br>
to less obvious, like how much<br>
<a rel="nofollow" href="https://water.usgs.gov/edu/rivers-contain-groundwater.html" target="_blank">groundwater is flowing into or out of a river</a>.<br>
Ecosystems have a role to play too. In desert environments,<br>
<a rel="nofollow" href="https://pubs.usgs.gov/wsp/1423/report.pdf" target="_blank">some plants depend on groundwater</a><br>
which can actually draw down groundwater around rivers, reducing flows<br>
or even causing surface flow to cease even though groundwater may<br>
continue to be available just under the surface. This is just a brief<br>
discussion of the complexities of the natural water cycle and the<br>
natural features that arise from it. For more, the<br>
<a rel="nofollow" href="https://water.usgs.gov/edu/watercycle.html" target="_blank">USGS Water Science School</a><br>
has a wealth of information and all<br>
<a rel="nofollow" href="https://pubs.er.usgs.gov/" target="_blank">USGS publications</a><br>
and<br>
<a rel="nofollow" href="https://data.usgs.gov/datacatalog/#fq=dataType%3A(collection%20OR%20non-collection)&amp;fq=keywords%3A%22NWIS%22&amp;q=*%3A*" target="_blank">data</a><br>
are free and available for everyone.</p>
<h2>Code Explanation</h2>
<p>The three scenarios above were run using a function that executes the<br>
code described below. There is a mix of comments and summary text<br>
written outside the code. To execute this, you will need to install the<br>
required packages and download the<br>
<a rel="nofollow" href="https://www.epa.gov/waterdata/nhdplus-national-data" target="_blank">NHDPlus National Seamless Database</a><br>
and a<br>
<a rel="nofollow" href="https://www.arcgis.com/home/item.html?id=b07a9393ecbd430795a6f6218443dccc" target="_blank">state boundaries layer available here.</a><br>
These two datasets are referenced in the code.</p>
<h3>Load Packages and Data</h3>
<p>In this first step, I load all the packages and data I need.</p>
<p>Packages include:</p>
<ul>
<li><a rel="nofollow" href="https://r-spatial.github.io/sf/" target="_blank"><code>sf</code></a></li>
<li><a rel="nofollow" href="https://dplyr.tidyverse.org/" target="_blank"><code>dplyr</code></a></li>
<li><a rel="nofollow" href="https://cran.r-project.org/web/packages/RANN/index.html" target="_blank"><code>RANN</code></a></li>
<li><a rel="nofollow" href="https://dblodgett-usgs.github.io/nhdplusTools/" target="_blank"><code>nhdplusTools</code></a></li>
<li><a rel="nofollow" href="https://r-quantities.github.io/units/" target="_blank"><code>units</code></a></li>
<li><a rel="nofollow" href="https://stringr.tidyverse.org/" target="_blank"><code>stringr</code></a></li>
<li><a rel="nofollow" href="https://cran.r-project.org/web/packages/gifski/index.html" target="_blank"><code>gifski</code></a></li>
<li><a rel="nofollow" href="https://cran.r-project.org/web/packages/ggmap/" target="_blank"><code>ggmap</code></a></li>
</ul>
<p>Note that <code>sf</code> and <code>dplyr</code> do the vast majority of the work and they are<br>
loaded with a <code>library()</code> command. All other packages are accessed using<br>
the <code>package::function()</code> syntax</p>
<p>Variables created here are:</p>
<ul>
<li><code>flowlines</code>: NHDPlus flowlines with lots of attributes. Loaded with<br>
<a rel="nofollow" href="https://dblodgett-usgs.github.io/nhdplusTools/" target="_blank"><code>nhdplusTools</code></a><br>
as a helper.</li>
<li><code>water_bodies</code>: NHDPlus waterbodies loaded directly from the<br>
National Seamless using<br>
<a rel="nofollow" href="https://r-spatial.github.io/sf/" target="_blank"><code>sf</code></a><br>
<code>read_sf()</code>.</li>
<li><code>min_monthlies</code>: Minimum of the twelve monthly flow estimates for<br>
each flowline. Calculated by applying the <code>min</code> function to all the<br>
monthly flow collumns from the <code>flowlines</code> attributes.</li>
<li><code>remove_fcodes</code>: NHDPlus feature codes that should be removed.</li>
<li><code>crs</code>: A coordinate reference system to perform the analysis in.</li>
<li><code>states</code>: A set of US state boundaries loaded from shapefile with<br>
<code>sf</code>.</li>
<li><code>bbox</code>: A bounding box to limit the analysis to less than the whole<br>
country. Created in the <code>sf</code> bbox format then converted in an <code>sfc</code><br>
geometry.</li>
</ul>
<pre><code>library(sf)
library(dplyr)
# to install nhdplusTools
# devtools::install_github(&quot;dblodgett-usgs/nhdplusTools&quot;)
# First we will use a nhdplusTools to load up the national seamless geodatabase.
nhdplusTools::nhdplus_path(&quot;nhdplus_data/NHDPlusV21_National_Seamless.gdb&quot;)
staged_data &lt;- nhdplusTools::stage_national_data(include = &quot;flowline&quot;, 
                                                 output_path = &quot;nhdplus_data/&quot;)
flowlines &lt;- readRDS(staged_data$flowline)

# Now lets read in the waterbodies directly from the national seamless database.
if(&quot;waterbodies&quot; %in% scenario) {
  water_bodies &lt;- read_sf(nhdplus_path(), &quot;NHDWaterbody&quot;)
}

if(&quot;filter_monthly_flow&quot; %in% scenario) {
  monthlies &lt;- which(grepl(&quot;QA_[0-1][0-9]&quot;, names(flowlines)))
  min_monthlies &lt;- apply(st_set_geometry(flowlines, NULL)[monthlies], 1, min)
}

if(&quot;remove_intermittent&quot; %in% scenario) {
  fcodes &lt;- read_sf(nhdplus_path(), &quot;NHDFCode&quot;)
  remove_fcodes &lt;- fcodes[which(fcodes$Hydrograph == &quot;Intermittent&quot;), ]
}

# http://spatialreference.org/ref/sr-org/epsg-5070/
crs &lt;- st_crs(5070)

# See: https://www.arcgis.com/home/item.html?id=b07a9393ecbd430795a6f6218443dccc for this file
states &lt;- read_sf(&quot;states_21basic/states.shp&quot;)

# http://bboxfinder.com/#25.284438,-127.265625,43.707594,-94.438477
bbox &lt;- c(-127.265625,25.284438,-94.438477,43.7075949)
names(bbox) &lt;- c(&quot;xmin&quot;, &quot;ymin&quot;, &quot;xmax&quot;, &quot;ymax&quot;)
class(bbox) &lt;- &quot;bbox&quot;
bbox &lt;- st_as_sfc(bbox)
st_crs(bbox) &lt;- 4326
</code></pre>
<h3>Transform and Filter</h3>
<p>In this next step, I transform all the data into a consistent coordinate<br>
reference system and get everything filtered. This is all the data I<br>
want to use for the actual nearest neighbor search. Processing includes:<br>
– <code>states</code>: Transform to analysis coordinate system and remove small<br>
islands. – <code>flowlines</code> and <code>water_bodies</code>: Filter out flowlines with min<br>
monthly flow of 0 and feature codes indicating intermittent. Transform<br>
to analysis coordinate system and intersect with analysis bounds and<br>
simplify geometry to remove un-needed data and precision.</p>
<p>At this point, I save the input flowlines and waterbodies to disk for<br>
use later because the process is destructive and I need to be careful<br>
how much system memory I’m using.</p>
<pre><code>bbox &lt;- st_transform(bbox, crs)

states &lt;- states %>%
  st_transform(crs) %>%
  st_cast(&quot;POLYGON&quot;) %>%
  mutate(AREA = st_area(.)) %>%
  filter(AREA > units::set_units(2500000000, &quot;m^2&quot;)) %>%
  select(-AREA)

if(&quot;filter_monthly_flow&quot; %in% scenario) {
  flowlines &lt;- flowlines[which(min_monthlies !=0 ), ]
}

if(&quot;remove_intermittent&quot; %in% scenario) {
  flowlines &lt;- flowlines[which(!flowlines$FCODE %in% remove_fcodes$FCode), ] 
  if(&quot;waterbodies&quot; %in% scenario) {
    water_bodies &lt;- water_bodies[which(!water_bodies$FCODE %in% remove_fcodes$FCode), ]
  }
}

flowlines &lt;- flowlines %>%
  st_transform(crs) %>%
  st_simplify(1000) %>%
  st_intersection(bbox)

if(&quot;waterbodies&quot; %in% scenario) {
  water_bodies &lt;-
    st_transform(water_bodies, crs) %>%
    st_buffer(0) %>%
    st_simplify(500) %>%
    st_intersection(bbox)
}

# Save some intermediate artifacts that we'll read back in later.
saveRDS(flowlines, &quot;temp_flowlines.rds&quot;)
if(&quot;waterbodies&quot; %in% scenario) saveRDS(water_bodies, &quot;temp_water_bodies.rds&quot;)
</code></pre>
<h3>Convert to coordinates</h3>
<p>Now that I have all my data in the right coordinate system and have<br>
removed all the data that I don’t want, I can set up the data for the<br>
actual nearest neighbor search. This code converts the data to<br>
coordinate pairs with a set of identifiers for each feature and feature<br>
part. Some of this code is not used, but is included to demonstrate how<br>
identifiers work in what comes out of <code>sf::st_coordinates</code>. The<br>
“MULTILINESTRING” coordinates have columns X, Y, L1, and L2 where L1 is<br>
the part number and L2 is the overall feature. The “MULTIPOLYGON”<br>
coordinates have columns X, Y, L1, L2, and L3 where L1 is for the main<br>
ring or holes, L2 for each polygon, and L3 for the overall feature. In<br>
this case, I just extract the identifier for the overall feature and<br>
could have used that information to track back to the actual water body<br>
or flowline in particular that I end up with at the end.</p>
<pre><code>if(&quot;waterbodies&quot; %in% scenario) wb_COMID &lt;- water_bodies$COMID
fl_COMID &lt;- flowlines$COMID

# Now turn both flowlines and water_bodies into coordinate pairs only
flowlines &lt;-
  st_cast(flowlines, &quot;MULTILINESTRING&quot;) %>%
  st_coordinates()
if(&quot;waterbodies&quot; %in% scenario) {
  water_bodies &lt;-
    st_cast(water_bodies, &quot;MULTIPOLYGON&quot;) %>%
    st_coordinates()
}

# Extract the identifier of features from the coordinates
flowlines &lt;- flowlines[,c(1,2,4)] %>%
  data.frame() %>%
  rename(ID = L2)
if(&quot;waterbodies&quot; %in% scenario) {
  water_bodies &lt;- water_bodies[,c(1,2,5)] %>%
    data.frame() %>%
    rename(ID = L3)
}

# Switch the new ID column to the &quot;COMID&quot; values from the source data.
if(&quot;waterbodies&quot; %in% scenario) {
  water_bodies[[&quot;COMID&quot;]] &lt;- wb_COMID[water_bodies[[&quot;ID&quot;]]]
}
flowlines[[&quot;COMID&quot;]] &lt;- fl_COMID[flowlines[[&quot;ID&quot;]]]

# Bind together into one huge set of coordinates.
if(&quot;waterbodies&quot; %in% scenario) {
  coords &lt;- rbind(water_bodies, flowlines) %>%
    select(-ID)
  rm(flowlines, water_bodies)
} else {
  coords &lt;- flowlines %>%
    select(-ID)
  rm(flowlines)
}
</code></pre>
<h3>Search Points</h3>
<p>In this short step, I create my set of all points using the function<br>
<a rel="nofollow" href="https://www.rdocumentation.org/packages/base/versions/3.5.1/topics/expand.grid?tap_a=5644-dce66f&tap_s=10907-287229" target="_blank"><code>expand.grid()</code></a><br>
and the extents of the bounding box I created above. Using the complete<br>
grid of locations, I convert it to as <code>sf</code> <code>data.frame</code> and remove all<br>
the points outside the analysis box and outside the states. I then<br>
convert the set of points inside the analysis area to a <code>data.frame()</code><br>
of coordinates.</p>
<pre><code># Create a set of search locations
search_bbox &lt;- st_bbox(bbox)
x &lt;- seq(search_bbox$xmin, search_bbox$xmax, 5000)
y &lt;- seq(search_bbox$ymin, search_bbox$ymax, 5000)
search &lt;- expand.grid(x,y)

# Convert to sf and intersect with the state boundaries.
search &lt;- st_as_sf(data.frame(search), coords = c(&quot;Var1&quot;, &quot;Var2&quot;), crs = crs) %>%
  st_intersection(states) %>%
  st_intersection(bbox) %>%
  st_coordinates() %>%
  data.frame()
</code></pre>
<h3>Plot function</h3>
<p>In order to create the animated gif, I needed to be able to create a<br>
similar plot over and over. So I created a function to do the job. This<br>
function creates an output file name, sets up a png, gets the data ready<br>
to plot, then plots a few layers using base R graphics and<br>
<a rel="nofollow" href="https://r-spatial.github.io/sf/reference/plot.html" target="_blank"><code>plot.sf</code></a>.</p>
<pre><code># Function to plot results in a .png
plot_result &lt;- function(search, radius, crs, bbox, states) {
  fname &lt;- stringr::str_pad(paste0(as.character(radius), &quot;.png&quot;), 10, &quot;left&quot;, &quot;0&quot;)
  png(fname,
      width = 1000, height = 800)
  result &lt;- st_as_sf(search, coords = c(&quot;X&quot;, &quot;Y&quot;), crs = crs)
  
  plot(bbox, main = paste(&quot;Distance to nearest water:&quot;, radius, &quot;meters&quot;))
  plot(states$geometry, add = TRUE)
  plot(result$geometry, pch = 19, add = TRUE)
  dev.off()
}
</code></pre>
<h3>Run the nearest neighbor search</h3>
<p>Finally, we are ready to run the analysis. To me, while loops are the<br>
forbidden fruit of software development, but sometimes you just don’t<br>
know how many times you need to run a loop! In this case, we run the<br>
<a rel="nofollow" href="https://www.rdocumentation.org/packages/RANN/versions/2.6/topics/nn2?tap_a=5644-dce66f&tap_s=10907-287229" target="_blank"><code>RANN::nn2()</code></a><br>
function which returns the index of the matching nearest neighbor<br>
(<code>nn.idx</code>), if one is found. If no nearest neighbor is found in the<br>
search radius, it returns <code>nn.idx</code> of 0. So in every loop, I only keep<br>
matches where <code>nn.idx</code> == 0 returns true. I plot the filtered set to a<br>
png and increment the radius up for the next loop. The if/else block at<br>
the bottom of this step just slows down the rate of increase of the<br>
radius as we approach only one location left to avoid overshooting and<br>
missing the last one!</p>
<p>Once the while loop finishes, I pass the list of pngs created to<br>
<a rel="nofollow" href="https://www.rdocumentation.org/packages/gifski/versions/0.8.3/topics/gifski?tap_a=5644-dce66f&tap_s=10907-287229" target="_blank"><code>gifski</code></a><br>
a great little package for creating animated gifs. Finally I printed out<br>
the lat/lon of the location found so I could go look at the location in<br>
google.</p>
<pre><code>radius &lt;- 0
num_left &lt;- nrow(search)

while(num_left > 1) {
  num_left &lt;- nrow(search)
  
  # This is where the magic happens.
  matched &lt;- RANN::nn2(coords[,1:2],
                       search,
                       k = 1,
                       searchtype = &quot;radius&quot;,
                       radius = radius) %>%
    data.frame()
  
  # This while loop is destructive. Only keep unmatched.
  search &lt;- filter(search, matched$nn.idx == 0)
  
  plot_result(search, radius, crs, bbox, states)
  
  if(num_left > 50) {
    radius &lt;- radius + 1000
  } else if (num_left > 3) {
    radius &lt;- radius + 500
  } else {
    radius &lt;- radius + 100
  }
}

rm(coords)

gifski::gifski(list.files(pattern = &quot;*0.png&quot;))

# Where in the world is the result?
result &lt;- st_as_sf(search, coords = c(&quot;X&quot;, &quot;Y&quot;), crs = crs)
print(st_transform(result$geometry, 4326))
</code></pre>
<h3>Plot the result location</h3>
<p>Finally, I plot up the results at a local scale. Steps for this include:<br>
– I want the original data I started with, so I just load it back from<br>
disk. – I set up an area for the plot by buffering around the result<br>
location by 150km and converting the result to a simple box. –<br>
Increasingly, I like using the web mercator projection for visualization<br>
because it’s familiar and compatible with web-map tiles. – I project all<br>
the data to my plotting projection and subset it and simplify the<br>
geometry to speed things up and not overwhelm <code>plot()</code> with data for the<br>
whole country. – I grab a background map from google with ggmap::getmap<br>
to be passed to the <code>bgMap</code> input of <code>plot.sf()</code>. – I plot each layer I<br>
want into a png. – Finally I delete the png files and temporary geometry<br>
I wrote to disk.</p>
<pre><code># Load geospatial data again.
flowlines &lt;- readRDS(&quot;temp_flowlines.rds&quot;)
if(&quot;waterbodies&quot; %in% scenario) water_bodies &lt;- readRDS(&quot;temp_water_bodies.rds&quot;)

# Set up a plot area around our result.
plot_area &lt;- st_buffer(result$geometry, 150000) %>%
  st_bbox() %>%
  st_as_sfc()
st_crs(plot_area) &lt;- crs

# Use 3857 (web mercator) for visualization.
plot_proj &lt;- 3857
plot_area &lt;- st_transform(plot_area, plot_proj)

# Subset data to an area that satisfies plot_area.
data_area &lt;- st_bbox(plot_area) %>%
  st_as_sfc()
st_crs(data_area) &lt;- st_crs(data_area)
data_area &lt;- st_transform(data_area, crs)

if(&quot;waterbodies&quot; %in% scenario) {
  local_water &lt;- st_intersection(water_bodies, data_area) %>%
    st_simplify(1000) %>%
    st_transform(plot_proj)
}
local_flowlines &lt;- st_intersection(flowlines, data_area) %>%
  st_simplify(5000) %>%
  st_transform(plot_proj)

bgmap &lt;- plot_area %>%
  st_transform(4326) %>%
  st_bbox() %>%
  setNames(c(&quot;left&quot;, &quot;bottom&quot;, &quot;right&quot;, &quot;top&quot;)) %>%
  ggmap::get_map(zoom = 7)

# write out a png of the local area.
png(&quot;furthest_water.png&quot;, width = 1024, height = 1024)
par(omi = c(0,0,0,0), mai = c(0,0,0,0))
plot(plot_area, col = NA, border = NA, xaxs = 'i', yaxs = 'i', bgMap = bgmap)
plot(st_transform(states$geometry, plot_proj), add = TRUE)
plot(st_transform(result$geometry, plot_proj), col = &quot;red&quot;, cex = 3, lwd = 4, add = TRUE)
plot(local_flowlines$Shape, add = TRUE, col = &quot;blue&quot;)
if(&quot;waterbodies&quot; %in% scenario) plot(local_water$Shape, add = TRUE, col = &quot;azure2&quot;)
dev.off()

unlink(&quot;temp_water_bodies.rds&quot;, force = TRUE)
unlink(&quot;temp_flowlines.rds&quot;, force = TRUE)
unlink(list.files(pattern = &quot;*0.png&quot;), force = TRUE)
</code></pre>
<p>That’s it! If you made it this far, thanks for taking the time! I hope<br>
this was helpful one way or another.</p>
	<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/furthest-water/"> 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/2_VTqByMvf8" height="1"></a></div>