Freshwater

1. Global Lakes and Wetlands Database (GLWD - WWF, Lehner & Dooll)

Drawing upon a variety of existing maps, data and information, WWF and the Centre for Environmental Systems Research, University of Kassel, Germany created the Global Lakes and Wetlands Database (GLWD). The combination of best available sources for lakes and wetlands on a global scale (1:1 to 1:3 million resolution), and the application of GIS functionality enabled the generation of a database which focuses in three coordinated levels on (1) large lakes and reservoirs, (2) smaller water bodies, and (3) wetlands.

Level 1 (GLWD-1) comprises the 3067 largest lakes (area ≥ 50 km2) and 654 largest reservoirs (storage capacity ≥ 0.5 km3) worldwide, and includes extensive attribute data. Access Level 1 data

Level 2 (GLWD-2) comprises permanent open water bodies with a surface area ≥ 0.1 km2 excluding the water bodies contained in GLWD-1. Access Level 2 data

The base layers:

# World map
world_map <- rnaturalearth::ne_countries(scale = 'small', returnclass = c("sf"))

####
# union by continents
world_map <- world_map %>% 
  dplyr::group_by(continent) %>% 
  summarise()
#plot(world_map)  

# Base map
world_map.s <- ggplot() +
  geom_sf(data = world_map, size = .2, fill = "gray95", col = "gray90") +
  theme(panel.grid.major = element_line(color = gray(0.95), linetype = "dashed", size = 0.5))+
  xlab("Longitude")+ 
  ylab("Latitude")+
  theme_minimal()



# Alternative basemap, for details (higher resolution):
world_map <- map_data('world')
world_map.b <- ggplot(data = world_map) +
  geom_polygon(data = world_map, aes(x=long, y = lat, group=group), fill = "gray80", colour = "gray70", alpha = .4) +
  theme(panel.grid.major = element_line(color = gray(0.95), linetype = "dashed", size = 0.5))+
  xlab("Longitude")+ 
  ylab("Latitude")+
  theme_minimal()

object.size(world_map.b)
## 4919152 bytes
object.size(world_map.s)
## 66576 bytes
# European Estuaries to fit in the details
EuTW <- sf::read_sf("D:/shapes/wise_transitional_waters_shp-20191011T113722Z-001/wise_transitional_waters_shp/Transitional_waters.shp") %>% st_transform(crs=4326)
GLWD <- sf::read_sf("D:/shapes/GLWD-level1/glwd_1.shp") 
GLWD <- GLWD %>%  
  st_set_crs(4326)

# global plot:
world_map.s +
    geom_sf(data=GLWD, aes(col=TYPE), alpha=.5)+
    ggtitle("WWF Global lakes and wetlands  (GLWD-level 1)")

# detail plot
world_map.b +
    geom_sf(data=EuTW, col="orange", alpha=.5)+
    geom_sf(data=GLWD, aes(col=TYPE, fill=TYPE), alpha=.4)+
    #theme(legend.position="none")+
    coord_sf(xlim = c(-10, 0), ylim = c(36,46))+
    ggtitle("Detail WWF Global lakes and wetlands  (GLWD-level 1)+EuTW")

##########################
# Level 2 DATA
GLWD2 <- sf::read_sf("D:/shapes/GLWD-level2/glwd_2.shp") %>%  st_set_crs(4326)
# sp::plot(GLWD2["TYPE"], col = 'blue', main="GLWD2")

# world_map.s +
#      geom_sf(data=GLWD2, aes(col=TYPE), alpha=.5)+
#      ggtitle("WWF Global lakes and wetlands  (GLWD-level 2)")

world_map.b +
    geom_sf(data=EuTW, col="orange", alpha=.5)+
    geom_sf(data=GLWD, aes(col=TYPE, fill=TYPE), alpha=.4)+
    geom_sf(data=GLWD2, aes(col=TYPE, fill=TYPE), alpha=.4)+
    #theme(legend.position="none")+
    coord_sf(xlim = c(-10, 0), ylim = c(36,46))+
    ggtitle("Detail WWF GLWD (level 1 & 2) + EuTW")

2. Lakes and rivers (naturalearth)

Natural earth has a database, that is accessible directly by r, for global lakes and rivers.

# Lakes
LAKES <- ne_download(scale = "large", type = 'lakes', category = 'physical', returnclass = "sf")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\lioca\AppData\Local\Temp\RtmpO21Oyt", layer: "ne_10m_lakes"
## with 1354 features
## It has 37 fields
## Integer64 fields read as strings:  scalerank ne_id
sp::plot(LAKES["name"], col = 'blue4', main="lakes")

# # Rivers
RIVERS <- ne_download(scale =  "large", type = 'rivers_lake_centerlines', category = 'physical', returnclass = "sf")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\lioca\AppData\Local\Temp\RtmpO21Oyt", layer: "ne_10m_rivers_lake_centerlines"
## with 1455 features
## It has 34 fields
## Integer64 fields read as strings:  rivernum ne_id
sp::plot(RIVERS["name"], col = 'blue', main="rivers_lake_centerlines")

# Lakes + Reservoirs (manually downloaded from https://www.naturalearthdata.com/downloads/10m-physical-vectors/10m-lakes/)
LARE <- sf::read_sf("D:/shapes/ne_50m_rivers_lake_centerlines/ne_50m_rivers_lake_centerlines.shp")
sp::plot(LARE["featurecla"], col = 'pink', main="ne_50m_rivers_lake_centerlines")

world_map.b +
    geom_sf(data=EuTW, col="orange3", fill="yellow", alpha=.8)+
    geom_sf(data=RIVERS, col="blue4", alpha=.4)+
    geom_sf(data=LAKES, col="orange", fill="orange", alpha=.5)+
    geom_sf(data=LARE, col="pink", alpha=.4)+
    #theme(legend.position="none")+
    coord_sf(xlim = c(-10, 0), ylim = c(36,46))+
    ggtitle("Rivers & Lakes (rnaturalearth)")

# further detailed map:
world_map.b +
    geom_sf(data=RIVERS, col="blue4", alpha=.4)+
    geom_sf(data=LAKES, col="orange", fill="orange", alpha=.5)+
    geom_sf(data=LARE, col="pink", alpha=.4)+
    geom_sf(data=EuTW, col="orange3", fill="yellow", alpha=.8)+
    #theme(legend.position="none")+
    coord_sf(xlim = c(-9, -6.5), ylim = c(36.5,37.5))+
    ggtitle("Rivers & Lakes (rnaturalearth)")

See also ‘osmdata’ library for adding rivers and lakes (example).

3. Water (OpenStreetMap)

OpenStreetMap represents physical features on the ground (e.g., roads or buildings) using tags attached to its basic data structures (its nodes, ways, and relations). Each tag describes a geographic attribute of the feature being shown by that specific node, way or relation. available features

The ones that can be relevant here would be:
-key=‘natural’; value= c(‘water’, ‘wetland’, ‘bay’, ‘lake’)

Note that puting key=natural and value = water, it gave all water features apperantly. Note also that osm retrieves polygons, lines, etc. each separatly.
“The natural key is used to describe wide variety of physical geography, geological and landcover features, including ones that have been modified by humans.”ref

-key=‘waterway’; value = c(‘river’,‘riverbank’,‘stream’,‘tidal_chanel’,‘canal’,‘drain’ ,‘ditch’, ‘pressurised’,‘fairway’)

“By definition, a waterway is assumed to have a direction of flow. The direction of the way should be downstream (from the waterway’s source to its mouth).
For small waterways, a way indicating a direction of flow and its path is sufficient. When a waterway becomes wide enough that mapping the area it flows through is desired.
If a waterway way is known to not have water flow year round, it should always be tagged with intermittent=yes. (Non-permanent flow can be further described with the seasonal=* tag to describe when a flow of water exists, but intermittent=yes should still be present in these cases.)
If a waterway is named from its source to its destination, it’s strongly suggested that all of its ways be placed in a waterway relation. Doing this allows Nominatim to group the ways together and return exactly one named result per named waterway that exists in OpenStreetMap.”ref.

-key=‘water’; value=‘lake’

Note also that there is a limit for data requests in osm.
Note also OSM has an excellent coastline.

# make a boundary box:
kk2 <- st_bbox(c(xmin = c(-10), xmax=c(-4), ymin = 36, ymax=40))

# save to use after - I did not managed to make this work with rmarkdown active. Thus, after some tries, decided to downlowd the data, save it and open it within rmarkdown (but the code used is bellow).

# load the saved data:
load("D:/shapes/osm_water_data.RData")

# extract data:
# osm_lakes.sf <-
#   opq(bbox = kk2) %>%
#   add_osm_feature(key = 'natural',value= c('water')) %>%
#   osmdata_sf()
# table(osm_lakes.sf$osm_multipolygons$water)
# sp::plot(osm_lakes.sf$osm_multipolygons["water"], col = 'blue2', main="Open Street Map - natural")

# plot_
world_map.b +
    geom_sf(data=osm_lakes.sf$osm_multipolygons, aes(col=name, fill=water), alpha=.4)+
    geom_sf(data=osm_lakes.sf$osm_polygons, aes(col=name, fill=water), alpha=.4)+
    geom_sf(data=osm_lakes.sf$osm_multilines, alpha=.4)+
    theme(legend.position="none")+
    coord_sf(xlim = c(-10, -4), ylim = c(36,40))+
    ggtitle("Natural: water (OpenStreetMap)")

## Wetland
# osm_wet.sf <-
#   opq(bbox = kk2) %>%
#   add_osm_feature(key = 'natural',value= c('wetland')) %>%
#   osmdata_sf()

world_map.b +
    geom_sf(data=osm_wet.sf$osm_multipolygons, aes(col=name, fill=wetland), alpha=.4)+
    geom_sf(data=osm_wet.sf$osm_polygons, aes(col=name, fill=wetland), alpha=.4)+
    geom_sf(data=osm_wet.sf$osm_lines, alpha=.4)+
    theme(legend.position="none")+
    coord_sf(xlim = c(-10, -4), ylim = c(36,40))+
    ggtitle("Natural: wetand (OpenStreetMap)")

## Tidal chanel (estuaries): no features.
## Bay 
# osm_bay.sf <-
#   opq(bbox = kk2) %>%
#   add_osm_feature(key = 'natural',value= c('bay')) %>%
#   osmdata_sf()
# osm_bay.sf

world_map.b +
    geom_sf(data=osm_bay.sf$osm_multipolygons, aes(col=name, fill=place), alpha=.4)+
    geom_sf(data=osm_bay.sf$osm_polygons, aes(col=name, fill=place), alpha=.4)+
    geom_sf(data=osm_bay.sf$osm_lines, alpha=.4)+
    theme(legend.position="none")+
    coord_sf(xlim = c(-10, -4), ylim = c(36,40))+
    ggtitle("Natural: bay (OpenStreetMap)")

# # data extraction code:
# osm_rivers.sf <-
#   opq(bbox = kk2) %>%
#   add_osm_feature(key = 'waterway', value=c('river','riverbank')) %>%
#   osmdata_sf()

# sp::plot(osm_rivers.sf$osm_lines["name"], col = 'blue', main="Open Street Map - waterway")
# sp::plot(osm_rivers.sf$osm_multipolygons["name"], col = 'orange', main="Open Street Map - waterway", add=TRUE)

# plot_
world_map.b +
    geom_sf(data=osm_rivers.sf$osm_multipolygons, aes(col=water, fill=water), alpha=.4)+
    geom_sf(data=osm_rivers.sf$osm_polygons, aes(col=name, fill=water), alpha=.4)+
    geom_sf(data=osm_rivers.sf$osm_multilines, aes(col=name, fill=water), alpha=.4)+
    theme(legend.position="none")+
    coord_sf(xlim = c(-10, -4), ylim = c(36,40))+
    ggtitle("Waterway (OpenStreetMap)")

# save(osm_rivers.sf, osm_wet.sf, osm_lakes.sf,osm_bay.sf, file="D:/shapes/osm_water_data.RData")

# All together:
world_map.b +
    geom_sf(data=osm_lakes.sf$osm_multipolygons, aes(col=water, fill=water), alpha=.4)+
    geom_sf(data=osm_lakes.sf$osm_polygons, aes(col=water, fill=water), alpha=.4)+
    geom_sf(data=osm_lakes.sf$osm_multilines, alpha=.4)+
    geom_sf(data=osm_wet.sf$osm_multipolygons, aes(col=wetland, fill=wetland), alpha=.4)+
    geom_sf(data=osm_wet.sf$osm_polygons, aes(col=wetland, fill=wetland), alpha=.4)+
    geom_sf(data=osm_wet.sf$osm_lines, alpha=.4)+
    geom_sf(data=EuTW, col="black", fill="yellow", alpha=.8)+
    geom_sf(data=osm_rivers.sf$osm_multipolygons, aes(col=water, fill=water), alpha=.4)+
    geom_sf(data=osm_rivers.sf$osm_polygons, aes(col=water, fill=water), alpha=.4)+
    theme(legend.position="none")+
    coord_sf(xlim = c(-10, -4), ylim = c(36,40))+
    ggtitle("Natural & Waterway (OpenStreetMap)+European transition waters (yellow)")

# ggplotly(p5)

rm(osm_bay.sf,osm_wet.sf,osm_lakes.sf)

4. Tidal Flat Ecosystems (Murray et al. 2019)

Murray, N. J., Phinn, S. R., DeWitt, M., Ferrari, R., Johnston, R., Lyons, M. B., … & Fuller, R. A. (2019). The global distribution and trajectory of tidal flats. Nature, 565(7738), 222.
“Here we present an analysis of over 700,000 satellite images that maps the global extent of and change in tidal flats over the course of 33 years (1984–2016). We find that tidal flats, defined as sand, rock or mud flats that undergo regular tidal inundation(…)”
“Our 30-m spatial-resolution dataset of global tidal flats is in the public domain, and can be accessed interactively via Google Earth Engine or via direct download (http://intertidal.app).”
“Landsat 4, 5, 7 and 8 satellite image acquired within 1-km of the coastline”
The dataset can be downloaded from UNEP webpage also.

Intertidal Change Explorer

The data is organized as geotiff images for each XX lat/long. I was able to open these in r using raster package (and make a stack of the time series), but not sure how to work with these.

This dataset, although very interesting, will not be useful to determine the location of the estuaries (thanks to José P. Granadeiro).

Ocean data viewer has also maps of the global distribution of sea-grasses, and Global Distribution of Salt marshes, which I think also do not sort out our issues, but are potentially interesting for other works.

5. Surface Waters Dynamics (Pekel et al., 2016)

Pekel, J. F., Cottam, A., Gorelick, N., & Belward, A. S. (2016). High-resolution mapping of global surface water and its long-term changes. Nature, 540(7633), 418.

In this paper the authors identify the location and persistence of global surface water (inland and coastal) and its dynamics, using three million Landsat satellite images13, over the past 32 years at 30-metre resolution, recording the months and years when water was present, where occurrence changed and what form changes took in terms of seasonality and persistence. Between 1984 and 2015 (…).

6. WISE Large rivers and large lakes

(WISE Large rivers and large lakes extracted)[https://www.eea.europa.eu/data-and-maps/data/wise-large-rivers-and-large-lakes/zipped-shapefile-with-wise-other-large-rivers-and-tributaries-vector-line/zipped-shapefile-with-wise-other-large-rivers-and-tributaries-vector-line]

See also: https://hydrosheds.cr.usgs.gov/dataavail.php

# open downloaded shape file:

WISE <- sf::read_sf("D:/shapes/wise_other_large_rivers_and_tributaries/Other_large_rivers_and_tributaries.shp") 
WISE <- WISE %>%  
  st_transform(4326)

# global plot:
world_map.b +
    geom_sf(data=WISE, aes(col=NAME), alpha=.5)+
    theme(legend.position="none")+
    coord_sf(xlim = c(-10, 10), ylim = c(36,52))+
    ggtitle("EU WISE Large rivers and large lakes (detail)")

# global plot:
world_map.b +
    geom_sf(data=WISE, aes(col=NAME), alpha=.5)+
    theme(legend.position="none")+
    coord_sf(xlim = c(-21, 34), ylim = c(36,70))+
    ggtitle("EU WISE Large rivers and large lakes (all)")

rm(WISE)

7. Global Self-consistant Hierarchical High-resolution (Shorelines) Geography, release 2.3.0 of February 1, 2014 (http://www.soest.hawaii.edu/pwessel/gshhg/gshhg-bin-2.3.0.zip)

Note this is the one used in GEBCO

  • WDBII_river__L.* shapefile

  • The shoreline data are distributed in 6 levels:

Level 1: Continental land masses and ocean islands, except Antarctica. Level 2: Lakes

(…)

# shoreline data
gshhg.l1 <- sf::read_sf("D:/shapes/gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp") %>% st_crop(xmin=-10, xmax=10, ymin=36, ymax=44)
plot(gshhg.l1["id"])

gshhg.l2 <- sf::read_sf("D:/shapes/gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L2.shp") %>% st_crop(xmin=-10, xmax=10, ymin=36, ymax=44)
plot(gshhg.l2["id"])

world_map.b +
    geom_sf(data=gshhg.l1, col="orange", fill="orange", alpha=.5)+
    geom_sf(data=gshhg.l2, col=2, fill="orange", alpha=.5)+
    coord_sf(xlim = c(-10, 0), ylim = c(36,46))+
    ggtitle("Borders & Lakes (gshhg)")

## interactive map
# gshhg.l1 %>% 
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  zcol = "id",
#                  alpha.regions = 0.3,
#                  legend = FALSE)+
# gshhg.l2 %>%  mapview::mapview(zcol = "id",legend = FALSE)


# WDBII
gshhg.river.l1 <- sf::read_sf("D:/shapes/gshhg-shp-2.3.7/WDBII_shp/f/WDBII_river_f_L01.shp") #%>% st_crop(xmin=-10, xmax=10, ymin=36, ymax=44)
gshhg.river.l2 <- sf::read_sf("D:/shapes/gshhg-shp-2.3.7/WDBII_shp/f/WDBII_river_f_L02.shp") 

world_map.b +
    geom_sf(data=gshhg.river.l1, col=4, alpha=.5)+
    geom_sf(data=gshhg.river.l2, col=3, alpha=.5)

rm(gshhg.river.l2, gshhg.river.l1)

Small example with reflief and population incorporated (just because I feel like it)

# adding rivers and relief to IBERIA (xlim = c(-10, 0), ylim = c(36,44)) Or in detail c(-10, -6), ylim = c(36,39)) detail http://www.francescobailo.net/2018/08/how-to-quickly-enrich-a-map-with-natural-and-anthropic-details/

# note that if we just want to change the coordinates, we just touh this chunk:
my_bbox <- c(xmin = -10, xmax = -6, ymin = 36, ymax = 39)

my_bbox.m <- 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 <- st_geometry(st_polygon(x = list(my_bbox.m)))
st_crs(my_bbox.sf) <- 4326
my_bbox_buff_2500.sf <- my_bbox.sf %>%
  st_transform(crs = 32632) %>%
  st_buffer(dist = 2500) %>% # 2.5 kilometers
  st_transform(crs = 4326)


##################
# relief data
##################
require(raster)
dem.raster <- getData("SRTM", lat = 39, lon = -9, download = TRUE)
dem.raster <- crop(dem.raster, as(my_bbox_buff_2500.sf, 'Spatial'), snap='out')
slope.raster <- terrain(dem.raster, opt='slope')
aspect.raster <- terrain(dem.raster, opt='aspect')
hill.raster <- hillShade(slope.raster, aspect.raster, 40, 270)
hill.m <- rasterToPoints(hill.raster)
hill.df <-  data.frame(hill.m)
colnames(hill.df) <- c("lon", "lat", "hill")


##################
##### RIVERs DATA
##################

library(osmdata)
osm_lakes.sf <-
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'water', value = 'lake') %>%
  osmdata_sf()
osm_lakes.sf <- osm_lakes.sf$osm_multipolygons

osm_rivers.sf <-
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'waterway', value = 'river') %>%
  osmdata_sf()
osm_rivers.sf <- osm_rivers.sf$osm_lines


#########
## Roads DATA:
#########

osm_roads_primary.sf <-
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'highway', value = 'trunk') %>%
  osmdata_sf()
osm_roads_primary.sf <- osm_roads_primary.sf$osm_lines

osm_roads_secondary.sf <-
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'highway', value = 'secondary') %>%
  osmdata_sf()
osm_roads_secondary.sf <- osm_roads_secondary.sf$osm_lines

osm_roads_tertiary.sf <-
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'highway', value = 'tertiary') %>%
  osmdata_sf()
osm_roads_tertiary.sf <- osm_roads_tertiary.sf$osm_lines

#########
## Buildings:
#########

osm_buildings.sf <-
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'building') %>%
  osmdata_sf()
osm_buildings.sf <- osm_buildings.sf$osm_polygons
osm_buildings_pnt.sf <-
  osm_buildings.sf %>%
  st_transform(crs = 32632) %>%
  st_centroid() %>%
  st_transform(crs = 4326)

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

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

osm_residential_areas_pnt.df <-
  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) <- c('lon', 'lat')

osm_residential_areas_pnt.sf <-
  st_as_sf(osm_residential_areas_pnt.df, coords = c("lon", "lat"), crs = 4326)
osm_residential_areas_pnt.sf <-
  st_transform(osm_residential_areas_pnt.sf, 32632)

###################
## polygons around the population areas
###################

library(dbscan)
res <-  dbscan(osm_residential_areas_pnt.df, eps = 0.0005, minPts = 10)

pop_dense_areas.sf <-
  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) <- 4326
pop_dense_areas.sf <- st_transform(pop_dense_areas.sf, 32632)

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

pop_dense_areas.sf <- st_transform(pop_dense_areas.sf, 4326)

library(dplyr)
osm_villages.sf <-
  opq(bbox = st_bbox(my_bbox_buff_2500.sf)) %>%
  add_osm_feature(key = 'place', value = 'village') %>%
  osmdata_sf()
osm_villages.sf <- osm_villages.sf$osm_points

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

osm_larger_places.df <-
  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 <- gsub('(.{1,10})(\\s|$)', '\\1\n', osm_larger_places.df$name)


##################
# plot everything
##################
p1 <- 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 = 'grey85', size = 0.1) +
  geom_sf(data = osm_roads_secondary.sf, colour = 'grey90', size = 0.08) +
  geom_sf(data = osm_roads_tertiary.sf, colour = 'grey95', size = 0.08) +
  #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()

p1

# p1 +
#       geom_sf(data=EuTW, col="blue2", alpha=.4)+ # EU transitional waters
#     geom_sf(data=LARE, col="pink", size=3, alpha=.4)+ # rnaturalearth
#     geom_sf(data=RIVERS, col="blue", alpha=.4)+ # rnaturalearth rivers
#     geom_sf(data=LAKES, col="red", alpha=.4)+ # rnaturalearth lakes
#     geom_sf(data=GLWD, aes(col=TYPE, fill=TYPE), alpha=.4)+
#     geom_sf(data=GLWD2, aes(col=TYPE, fill=TYPE), alpha=.4)+
#     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']))