<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/how-to-quickly-enrich-a-map-with-natural-and-anthropic-details/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.francescobailo.net/2018/08/how-to-quickly-enrich-a-map-with-natural-and-anthropic-details/"> R – : Francesco Bailo :</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><a rel="nofollow" href="http://www.francescobailo.net/wordpress/wp-content/uploads/2018/08/Enrich-your-map-with-natural-and-anthropic-details.png" target="_blank"><br>
<img sizes="(max-width: 853px) 100vw, 853px" src="https://i2.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/Enrich-your-map-with-natural-and-anthropic-details-853x1024.png?w=450" alt="" class="wp-image-945" data-recalc-dims="1"></a></p>
<p>In this post I show how to enrich a ggplot map with data obtained from the Open Street Map (OSM) API. After adding elevation details to the map, I add water bodies and elements identifying human activity. To highlight the areas more densely inhabitated, I propose to use a density-based clustering algorithm of OSM features.</p>
<p><span></span></p>
<div>
<h2>Data</h2>
<p>Let’s first start with some data: a few points somewhere in Italy.</p>
<pre><code>my_points.df &lt;-
  data.frame(lon = c(9.442366, 9.30222, 9.304297, 9.322077, 
                     9.304432, 9.302683, 9.321543, 9.301541, 9.402329, 9.285115, 9.319564, 
                     9.33959, 9.33367, 9.261872, 9.308274, 9.3011, 9.45355, 9.32527, 
                     9.25892, 9.303647), 
             lat = c(45.952427, 46.041, 45.958035, 45.922367, 
                     45.953438, 45.966248, 45.919487, 45.959721, 46.05933, 46.009094, 
                     45.917286, 46.081585, 45.997231, 45.98764, 46.076529, 45.959641, 
                     45.972649, 45.9151, 46.055069, 46.07671))</code></pre>
<p>I will use the great <a rel="nofollow" href="https://cran.r-project.org/web/packages/sf/index.html" target="_blank"><code>sf</code></a> package so let’s create a <code>sf</code> spatial object with</p>
<pre><code>library(sf)
my_points.sf &lt;- st_as_sf(my_points.df, coords = c(&quot;lon&quot;, &quot;lat&quot;), crs = 4326)</code></pre>
<p><code>crs = 4326</code> sets the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Spatial_reference_system" target="_blank">Coordinate Reference System</a> using an integer <a rel="nofollow" href="https://www.epsg-registry.org/" target="_blank">EPSG</a> code. <a rel="nofollow" href="https://epsg.io/4326" target="_blank">EPSG:4326</a> uses the World Geodetic System 1984 (WGS 84) that is used in GPS and its unit is the <em>degree</em> (not the meter!).</p>
<p>In this post I show how to enrich a ggplot map with data obtained from the Open Street Map (OSM) API. After adding elevation details to the map, I add water bodies and elements identifying human activity. To highlight the areas more densely inhabitated, I propose to use a density-based clustering algorithm of OSM features.</p>
<p>As a first step let’s define a bounding box that includes my points.</p>
<pre><code>my_bbox &lt;- c(xmin = min(my_points.df$lon),
             xmax = max(my_points.df$lon),
             ymin = min(my_points.df$lat),
             ymax = max(my_points.df$lat))</code></pre>
<p>Once I have set the points, I can create a polygon with the function <code>sf::st_polygon()</code>, which takes a two-column matrix for each polygon and create a simple feature (sf) object, and the function <code>sf::st_geometry()</code>, which takes a simple future object and get its geometry. Points need to be five because we always need to close our polygons! So the last (or 5th) row will have the same pair of coordinates of the first row. <code>st_crs</code> will set the CRS for the geometry <code>my_bbox.sf</code>.</p>
<pre><code>my_bbox.m &lt;- 
  matrix(c(my_bbox['xmin'], my_bbox['xmin'], my_bbox['xmax'], my_bbox['xmax'], my_bbox['xmin'], 
           my_bbox['ymax'], my_bbox['ymin'], my_bbox['ymin'], my_bbox['ymax'], my_bbox['ymax']),
         ncol = 2)
my_bbox.sf &lt;- st_geometry(st_polygon(x = list(my_bbox.m)))
st_crs(my_bbox.sf) &lt;- 4326</code></pre>
<p>To give some margin to my points, I create an additional buffered bbox polygon around <code>my_bbox.sf</code>. I want to specify the buffer distance in meters, so I need to transform my polygon to a meter-based coordinate system such as <a rel="nofollow" href="https://epsg.io/32632" target="_blank">EPSG:32632</a> and then back to EPSG:4326. EPSG:32632 uses the same geodetic CRS of EPSG:4326 (WGS 84) but crucially its unit is the <em>meter</em> (not the degree!) so all calculations are meter-based. EPSG:32632 is the 32N zone of the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system" target="_blank">Universal Transverse Mercator</a> (UTM) coordinate system. The world grid that can help you find the best UTM zone for your project is <a rel="nofollow" href="http://www.dmap.co.uk/utmworld.htm" target="_blank">here</a>.</p>
<pre><code>my_bbox_buff_2500.sf &lt;- 
  my_bbox.sf %>%
  st_transform(crs = 32632) %>%
  st_buffer(dist = 2500) %>% # 2.5 kilometers
  st_transform(crs = 4326)

my_bbox_buff_5000.sf &lt;- 
  my_bbox.sf %>%
  st_transform(crs = 32632) %>%
  st_buffer(dist = 5000) %>% # 5 kilometers
  st_transform(crs = 4326)

my_bbox_buff_25000.sf &lt;- 
  my_bbox.sf %>%
  st_transform(crs = 32632) %>%
  st_buffer(dist = 25000) %>% # 25 kilometers
  st_transform(crs = 4326)</code></pre>
<p>This is how my point data look like in its essence.</p>
<pre><code>library(ggplot2)

my_world_map &lt;- map_data('world')
my_world_map &lt;- my_world_map[my_world_map$region %in% c(&quot;Italy&quot;,&quot;Switzerland&quot;),]

ggplot() + 
  geom_sf(data = my_points.sf) + 
  geom_sf(data = my_bbox_buff_2500.sf, fill = NA) +
  coord_sf(xlim = c(st_bbox(my_bbox_buff_25000.sf)['xmin'], st_bbox(my_bbox_buff_25000.sf)['xmax']), 
           ylim = c(st_bbox(my_bbox_buff_25000.sf)['ymin'], st_bbox(my_bbox_buff_25000.sf)['ymax'])) +
  geom_polygon(data = my_world_map, 
               aes(x=long, y = lat, group = group,
                   fill = region), colour = 'black', alpha = .4) +
  theme_bw()</code></pre>
<p><a rel="nofollow" href="https://i0.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-7-1.png" target="_blank"><img sizes="(max-width: 1344px) 100vw, 1344px" src="https://i0.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-7-1.png?w=450" alt="" class="wp-image-939" data-recalc-dims="1"></a></p>
<p>To give some geographic context to the data I plotted also the country polygons of Italy and Switzerland from the <code>ggplot2::map_data('world')</code>. <code>coord_sf</code> allows to set the limits of the plot which are extracted from the bounding box of <code>my_bbox_buff_25000.sf</code> with the function <code>sf::st_bbox()</code>.</p>
</div>
<p><span></span></p>
<p></p>
<p></p>
<p></p>
<div>
<h2>Add elevation details</h2>
<p>Depending on the specific geography of the place you need to map, it might be helpful to provide an idea of the elevation. Elevation differences (that is, hills, mountains and valleys) are crucial constraints in the development of human and animal activities. To map the elevation around our points, I need to download data of a Digital Elevation Model (DEM). There are different options, I will use here data obtained by the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Shuttle_Radar_Topography_Mission" target="_blank">Shuttle Radar Topography Mission</a> (SRTM). The package <a rel="nofollow" href="https://cran.r-project.org/web/packages/raster/index.html" target="_blank"><code>raster</code></a> ships with the function (<code>raster::getData()</code>) to download the data directly from the R console. The arguments <code>lat</code> and <code>lon</code> are used to get the data file that is needed for the project.</p>
<pre><code>library(raster)
dem.raster &lt;- getData(&quot;SRTM&quot;, lat = 46.0146, lon = 9.344197, download = TRUE)</code></pre>
<p>The SRTM data file is actually pretty large so it probably a good idea to crop it to reduce plotting time. The function <code>raster::crop()</code> takes an sp object (form the <a rel="nofollow" href="https://cran.r-project.org/web/packages/sp/index.html" target="_blank"><code>sp</code></a> package) to be used as mask to determine the extent of the resulting (cropped) raster. I need to convert my bbox first with <code>as(my_bbox_buff_25000.sf, 'Spatial')</code>.</p>
<pre><code>dem.raster &lt;- crop(dem.raster, as(my_bbox_buff_25000.sf, 'Spatial'), snap='out')</code></pre>
<p>Because I need to add my raster as a ggplot layer, I need to transform it first into a matrix of points with <code>raster::rasterToPoints()</code> and finally into a <code>data.frame</code> object. The third column of the <code>dem.df</code> contains the altitude value (<code>alt</code>).</p>
<pre><code>dem.m  &lt;-  rasterToPoints(dem.raster)
dem.df &lt;-  data.frame(dem.m)
colnames(dem.df) = c(&quot;lon&quot;, &quot;lat&quot;, &quot;alt&quot;)</code></pre>
<p>I can already plot <code>dem.df</code> with</p>
<pre><code>ggplot() +
  geom_raster(data = dem.df, aes(lon, lat, fill = alt), alpha = .45) +
  scale_fill_gradientn(colours = terrain.colors(100)) +
  geom_sf(data = my_bbox_buff_2500.sf, fill = NA) +
  coord_sf(xlim = c(st_bbox(my_bbox_buff_25000.sf)['xmin'], st_bbox(my_bbox_buff_25000.sf)['xmax']), 
           ylim = c(st_bbox(my_bbox_buff_25000.sf)['ymin'], st_bbox(my_bbox_buff_25000.sf)['ymax'])) +
  geom_polygon(data = my_world_map, 
               aes(x=long, y = lat, group = group), fill = NA, colour = 'black') +
  theme_bw()</code></pre>
<p><a rel="nofollow" href="https://i2.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-11-1.png" target="_blank"><img sizes="(max-width: 1344px) 100vw, 1344px" src="https://i2.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-11-1.png?w=450" alt="" class="wp-image-940" data-recalc-dims="1"></a></p>
<p>Now, instead of colouring each cell of the raster based on its altitude, I want to convey the same information with the hill shade. First, I extract the terrain characteristics (the <code>slope</code> and the <code>aspect</code> in this case) with <code>raster::terrain()</code>. Then, I compute the hill shade with <code>raster::hillShade()</code>, setting the elevation angle (of the sun) to <code>40</code> and the direction angle of the light to <code>270</code>. I then transform the resulting <code>hill.raster</code> into a <code>data.frame</code> as before. (Thanks to <a rel="nofollow" href="https://stackoverflow.com/a/21487844" target="_blank">tim riffe</a> for this)</p>
<pre><code>slope.raster &lt;- terrain(dem.raster, opt='slope')
aspect.raster &lt;- terrain(dem.raster, opt='aspect')
hill.raster &lt;- hillShade(slope.raster, aspect.raster, 40, 270)

hill.m &lt;- rasterToPoints(hill.raster)
hill.df &lt;-  data.frame(hill.m)
colnames(hill.df) &lt;- c(&quot;lon&quot;, &quot;lat&quot;, &quot;hill&quot;)</code></pre>
<p>And here we go</p>
<pre><code>ggplot() +
  geom_raster(data = hill.df, aes(lon, lat, fill = hill), alpha = .45) +
  scale_fill_gradientn(colours = grey.colors(100)) +
  geom_sf(data = my_bbox_buff_2500.sf, fill = NA) +
  coord_sf(xlim = c(st_bbox(my_bbox_buff_25000.sf)['xmin'], st_bbox(my_bbox_buff_25000.sf)['xmax']), 
           ylim = c(st_bbox(my_bbox_buff_25000.sf)['ymin'], st_bbox(my_bbox_buff_25000.sf)['ymax'])) +
  geom_polygon(data = my_world_map, 
               aes(x=long, y = lat, group = group), fill = NA, colour = 'black') +
  theme_bw()</code></pre>
<p><a rel="nofollow" href="https://i1.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-13-1.png" target="_blank"><img sizes="(max-width: 1344px) 100vw, 1344px" src="https://i1.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-13-1.png?w=450" alt="" class="wp-image-941" data-recalc-dims="1"></a></p>
</div>
<div>
<h2>Add bodies of water</h2>
<p>As you might have guessed from the presence of a few flat surfaces, there are lakes in the region. We can use the excellent <a rel="nofollow" href="https://cran.r-project.org/web/packages/osmdata/index.html" target="_blank"><code>osmdata</code></a> package to query the <a rel="nofollow" href="https://www.openstreetmap.org/" target="_blank">Open Street Map (OSM)</a> <a rel="nofollow" href="https://wiki.openstreetmap.org/wiki/Overpass_API" target="_blank">Overpass API</a> to download information about lakes and rivers. With the function <code>osmdata::opq()</code> we define the bbox for the query and with <code>osmdata::add_osm_feature()</code> we extract specific OSM features using the <a rel="nofollow" href="https://wiki.openstreetmap.org/wiki/Tags#Keys_and_values" target="_blank">tags</a> (<code>key</code> and <code>value</code>).</p>
<pre><code>library(osmdata)
osm_lakes.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_25000.sf)) %>%
  add_osm_feature(key = 'water', value = 'lake') %>%
  osmdata_sf()
osm_lakes.sf &lt;- osm_lakes.sf$osm_multipolygons

osm_rivers.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_25000.sf)) %>%
  add_osm_feature(key = 'waterway', value = 'river') %>%
  osmdata_sf()
osm_rivers.sf &lt;- osm_rivers.sf$osm_lines

ggplot() +
  geom_raster(data = hill.df, aes(lon, lat, fill = hill), alpha = .45) +
  scale_fill_gradientn(colours = grey.colors(100)) +
  geom_sf(data = osm_lakes.sf, fill = '#9ecae1', colour = NA) +
  geom_sf(data = osm_rivers.sf, colour = '#9ecae1', size = 0.05) +
  geom_sf(data = my_bbox_buff_2500.sf, fill = NA) +
  coord_sf(xlim = c(st_bbox(my_bbox_buff_25000.sf)['xmin'], st_bbox(my_bbox_buff_25000.sf)['xmax']), 
           ylim = c(st_bbox(my_bbox_buff_25000.sf)['ymin'], st_bbox(my_bbox_buff_25000.sf)['ymax'])) +
  geom_polygon(data = my_world_map, 
               aes(x=long, y = lat, group = group), fill = NA, colour = 'black') +
  theme_bw()</code></pre>
<p><a rel="nofollow" href="https://i1.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-14-1.png" target="_blank"><img sizes="(max-width: 1344px) 100vw, 1344px" src="https://i1.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-14-1.png?w=450" alt="" class="wp-image-942" data-recalc-dims="1"></a></p>
</div>
<div>
<h2>Add roads</h2>
<p>Similarly we can download details about roads. The road network on OSM might be very dense, so it is better to do some filtering using the OSM <a rel="nofollow" href="https://wiki.openstreetmap.org/wiki/Key:highway" target="_blank">tags</a>.</p>
<pre><code>osm_roads_primary.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_5000.sf)) %>%
  add_osm_feature(key = 'highway', value = 'trunk') %>%
  osmdata_sf()
osm_roads_primary.sf &lt;- osm_roads_primary.sf$osm_lines

osm_roads_secondary.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_5000.sf)) %>%
  add_osm_feature(key = 'highway', value = 'secondary') %>%
  osmdata_sf()
osm_roads_secondary.sf &lt;- osm_roads_secondary.sf$osm_lines

osm_roads_tertiary.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_5000.sf)) %>%
  add_osm_feature(key = 'highway', value = 'tertiary') %>%
  osmdata_sf()
osm_roads_tertiary.sf &lt;- osm_roads_tertiary.sf$osm_lines

ggplot() +
  geom_raster(data = hill.df, aes(lon, lat, fill = hill), alpha = .45) +
  scale_fill_gradientn(colours = grey.colors(100)) +
  geom_sf(data = osm_lakes.sf, fill = '#9ecae1', colour = NA) +
  geom_sf(data = osm_rivers.sf, colour = '#9ecae1', size = 0.05) +
  geom_sf(data = osm_roads_primary.sf, colour = '#636363', size = 0.1) +
  geom_sf(data = osm_roads_secondary.sf, colour = '#636363', size = 0.05) +
  geom_sf(data = osm_roads_tertiary.sf, colour = '#636363', size = 0.02) +
  geom_sf(data = my_bbox_buff_2500.sf, fill = NA) +
  coord_sf(xlim = c(st_bbox(my_bbox_buff_5000.sf)['xmin'], st_bbox(my_bbox_buff_5000.sf)['xmax']), 
           ylim = c(st_bbox(my_bbox_buff_5000.sf)['ymin'], st_bbox(my_bbox_buff_5000.sf)['ymax'])) +
  geom_polygon(data = my_world_map, 
               aes(x=long, y = lat, group = group), fill = NA, colour = 'black') +
  theme_bw()</code></pre>
<p><a rel="nofollow" href="https://i0.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-15-1.png" target="_blank"><img sizes="(max-width: 1344px) 100vw, 1344px" src="https://i0.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-15-1.png?w=450" alt="" class="wp-image-943" data-recalc-dims="1"></a></p>
</div>
<div>
<h2>Identify densely populated areas</h2>
<p>Densely populated areas are certainly very important in every map that wants to describe human activities. Again, I use OSM data to identify these areas. First, I query for features that are usually associated with human residency: residential roads and buildings. On OSM roads and buildings are <a rel="nofollow" href="https://stefanodesabbata.com/2014/05/28/openstreetmap/" target="_blank">not mapped with the same completeness everywhere</a> but the density of their presence might provide a good enough approximation of where people actually live. The number of such features (polygons and lines) is usually very large. But actually I don’t need <em>individual features</em>, I just need <em>the geography of densely populated areas</em>.</p>
<p>First, I query the OSM API for all the buildings (<code>osm_polygons</code>), residential roads and unclassified roads (<code>osm_lines</code>). Second, I transform polygons and lines into points with <code>sf::st_centroid()</code> (which requires a CRS transformation). Finally, I store all the point simple feature objects into <code>osm_residential_areas_pnt.df</code>, with two columns respectively for longitude and latitude, and create from it the spatial object <code>osm_residential_areas_pnt.sf</code>.</p>
<pre><code>osm_buildings.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'building') %>%
  osmdata_sf()
osm_buildings.sf &lt;- osm_buildings.sf$osm_polygons
osm_buildings_pnt.sf &lt;- 
  osm_buildings.sf %>%
  st_transform(crs = 32632) %>%
  st_centroid() %>%
  st_transform(crs = 4326)

osm_roads_residential.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_5000.sf)) %>%
  add_osm_feature(key = 'highway', value = 'residential') %>%
  osmdata_sf()
osm_roads_residential.sf &lt;- osm_roads_residential.sf$osm_lines
osm_roads_residential_pnt.sf &lt;- 
  osm_roads_residential.sf %>%
  st_transform(crs = 32632) %>%
  st_centroid() %>%
  st_transform(crs = 4326)

osm_roads_unclassified.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_5000.sf)) %>%
  add_osm_feature(key = 'highway', value = 'unclassified') %>%
  osmdata_sf()
osm_roads_unclassified.sf &lt;- osm_roads_unclassified.sf$osm_lines
osm_roads_unclassified_pnt.sf &lt;- 
  osm_roads_unclassified.sf %>%
  st_transform(crs = 32632) %>%
  st_centroid() %>%
  st_transform(crs = 4326)

osm_residential_areas_pnt.df &lt;- 
  rbind(
    do.call(rbind, 
          st_geometry(osm_buildings_pnt.sf)),
        do.call(rbind, 
          st_geometry(osm_roads_residential_pnt.sf)),
        do.call(rbind, 
          st_geometry(osm_roads_unclassified_pnt.sf))) %>%
  as.data.frame() 
colnames(osm_residential_areas_pnt.df) &lt;- c('lon', 'lat')

osm_residential_areas_pnt.sf &lt;- 
  st_as_sf(osm_residential_areas_pnt.df, coords = c(&quot;lon&quot;, &quot;lat&quot;), crs = 4326)
osm_residential_areas_pnt.sf &lt;- 
  st_transform(osm_residential_areas_pnt.sf, 32632)</code></pre>
<p>Once I have the coordinate <code>data.frame</code>, with one row for each OSM feature, I pass it to the function <code>dbscan::dbscan()</code> of the <a rel="nofollow" href="https://cran.r-project.org/web/packages/sf/index.html" target="_blank"><code>dbscan</code></a> package, which implements a <em>Density-based spatial clustering</em> algorithm to identify areas with the highest density of points . <code>eps</code> and <code>minPts</code>, respectively the size of the neighbourhood and the number of minimum points for each region, should be adjusted accordingly to the local geography. Once the number the of clusters are identified, I can create a polygon with the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Convex_hull" target="_blank"><code>convex hull</code></a> of the clustered points with the function <code>grDevices::chull()</code>. Each created polygon is stored the spatial object <code>pop_dense_areas.sf</code>.</p>
<pre><code>library(dbscan)
res &lt;-  dbscan(osm_residential_areas_pnt.df, eps = 0.0005, minPts = 10)

pop_dense_areas.sf &lt;- 
  st_sf(id = 1:max(res$cluster), 
        geometry = st_sfc(lapply(1:max(res$cluster), function(x) st_geometrycollection())))
st_crs(pop_dense_areas.sf) &lt;- 4326
pop_dense_areas.sf &lt;- st_transform(pop_dense_areas.sf, 32632)

for (cl in 1:max(res$cluster)) {
  these_points &lt;- osm_residential_areas_pnt.df[which(res$cluster == cl),]
  this_chull &lt;- chull(these_points)
  this_mat &lt;- as.matrix(these_points[this_chull,])
  this_mat &lt;- rbind(this_mat, this_mat[1,])
  this_polygon &lt;- st_geometry(st_polygon(x = list(this_mat)))
  st_crs(this_polygon) &lt;- 4326
  this_polygon &lt;- st_transform(this_polygon, 32632)
  pop_dense_areas.sf$geometry[cl] &lt;- this_polygon
}

pop_dense_areas.sf &lt;- st_transform(pop_dense_areas.sf, 4326)</code></pre>
<p>To give an idea of the actual result, you see below a densely populated area inferred from the cluster of points derived from building and residential roads. The algorithm successfully identified relatively more anthorpised regions. Tweaking of the two arguments of the density-based clustering function might be necessary depending on the size and human geography of the final map.</p>
<p><a rel="nofollow" href="https://i1.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-18-1.png" target="_blank"><img sizes="(max-width: 1344px) 100vw, 1344px" src="https://i1.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/unnamed-chunk-18-1.png?w=450" alt="" class="wp-image-944" data-recalc-dims="1"></a></p>
<p>The final enrichment of my map is the addition of a few names of populated places (<em>towns</em> and <em>villages</em> in the OSM tagging system). Again I use the OSM API. I filter the results for the more important places based on their relative population. To avoid overcrowding the map with labels, I only select the 10 larger places. <code>geom_sf</code> doesn’t allow to label points (yet). It is then necessary to convert the spatial object into a simple <code>data.frame</code> with two coordinate columns. Names might be too long, so I add a <code>\n</code> (newline) after 10 characters with <code>gsub('(.{1,10})(\\s|$)', '\\1\n', str)</code>. (Thanks to <a rel="nofollow" href="https://stackoverflow.com/a/2352006" target="_blank">xiechao</a> for this)</p>
<pre><code>library(dplyr)
osm_villages.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'place', value = 'village') %>%
  osmdata_sf()
osm_villages.sf &lt;- osm_villages.sf$osm_points

osm_towns.sf &lt;- 
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'place', value = 'town') %>%
  osmdata_sf()
osm_towns.sf &lt;- osm_towns.sf$osm_points

osm_larger_places.df &lt;- 
  rbind(cbind(as.data.frame(osm_villages.sf)[,c('name','population')],
              st_coordinates(osm_villages.sf)),
        cbind(as.data.frame(osm_towns.sf)[,c('name','population')], 
              st_coordinates(osm_towns.sf))) %>%
  mutate(population = as.numeric(as.character(population))) %>%
  top_n(10, population)

osm_larger_places.df$name &lt;- gsub('(.{1,10})(\\s|$)', '\\1\n', osm_larger_places.df$name)</code></pre>
<p>And now let’s plot everything along with my initial data points:</p>
<pre><code>ggplot() +
  
  geom_raster(data = hill.df, aes(lon, lat, fill = hill), alpha = .45) +
  scale_fill_gradientn(colours = grey.colors(100)) +
  
  geom_sf(data = pop_dense_areas.sf, fill = '#fee08b', colour = NA, alpha = 0.9) +
  geom_sf(data = osm_lakes.sf, fill = '#9ecae1', colour = NA) +
  geom_sf(data = osm_rivers.sf, colour = '#9ecae1', size = 0.4) +
  geom_sf(data = osm_roads_primary.sf, colour = '#636363', size = 0.1) +
  geom_sf(data = osm_roads_secondary.sf, colour = '#636363', size = 0.08) +
  geom_sf(data = osm_roads_tertiary.sf, colour = '#636363', size = 0.08) +
  geom_sf(data = my_points.sf, shape = 5, colour = 'black', size = .5) +
  
  geom_text(data = osm_larger_places.df, aes(x=X,y=Y,label=name), size = 2.5, alpha = .65) +
  
  coord_sf(xlim = c(st_bbox(my_bbox_buff_2500.sf)['xmin'], st_bbox(my_bbox_buff_2500.sf)['xmax']), 
           ylim = c(st_bbox(my_bbox_buff_2500.sf)['ymin'], st_bbox(my_bbox_buff_2500.sf)['ymax'])) +
  guides(fill=FALSE)+
  labs(x=NULL, y=NULL) +
  theme_bw()</code></pre>
<p><a rel="nofollow" href="https://i0.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/Enrich-your-map-with-natural-and-anthropic-details.png" target="_blank"><img sizes="(max-width: 2952px) 100vw, 2952px" src="https://i0.wp.com/www.francescobailo.net/wordpress/wp-content/uploads/2018/08/Enrich-your-map-with-natural-and-anthropic-details.png?w=450" alt="" class="wp-image-945" data-recalc-dims="1"></a></p>
</div>
<div>
<h1>References</h1>
<div>
<div>
<p>Hahsler, Michael, and Matthew Piekenbrock. 2018. <em>Dbscan: Density Based Clustering of Applications with Noise (Dbscan) and Related Algorithms</em>. <a rel="nofollow" href="https://cran.r-project.org/package=dbscan" target="_blank">https://CRAN.R-project.org/package=dbscan</a>.</p>
</div>
<div>
<p>Hijmans, Robert J. 2017. <em>Raster: Geographic Data Analysis and Modeling</em>. <a rel="nofollow" href="https://cran.r-project.org/package=raster" target="_blank">https://CRAN.R-project.org/package=raster</a>.</p>
</div>
<div>
<p>Padgham, Mark, Bob Rudis, Robin Lovelace, and Maëlle Salmon. 2018. <em>Osmdata: Import ’Openstreetmap’ Data as Simple Features or Spatial Objects</em>. <a rel="nofollow" href="https://cran.r-project.org/package=osmdata" target="_blank">https://CRAN.R-project.org/package=osmdata</a>.</p>
</div>
<div>
<p>Pebesma, Edzer. 2018. <em>Sf: Simple Features for R</em>. <a rel="nofollow" href="https://cran.r-project.org/package=sf" target="_blank">https://CRAN.R-project.org/package=sf</a>.</p>
</div>
<div>
<p>Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, and Kara Woo. 2018. <em>Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics</em>. <a rel="nofollow" href="https://cran.r-project.org/package=ggplot2" target="_blank">https://CRAN.R-project.org/package=ggplot2</a>.</p>
</div>
<div>
<p>Wickham, Hadley, Romain Francois, Lionel Henry, and Kirill Müller. 2017. <em>Dplyr: A Grammar of Data Manipulation</em>. <a rel="nofollow" href="https://cran.r-project.org/package=dplyr" target="_blank">https://CRAN.R-project.org/package=dplyr</a>.</p>
</div>
<div>
<p>Xie, Yihui. 2018. <em>Knitr: A General-Purpose Package for Dynamic Report Generation in R</em>. <a rel="nofollow" href="https://cran.r-project.org/package=knitr" target="_blank">https://CRAN.R-project.org/package=knitr</a>.</p>
</div>
</div>
</div>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.francescobailo.net/2018/08/how-to-quickly-enrich-a-map-with-natural-and-anthropic-details/"> R – : Francesco Bailo :</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/WqK8o9TNhfE" height="1"></a></div>