Definitions: what is an estuary?

The definition of estuaries is complex, in view of its strong dynamic in time and space. However, it generally passes by a water mass with periodic salinity gradient and characterised by a particular topography. Within this broad definition, different water bodies can be included, like transition zones, wetlands, bays, fjords, coastal hydrosystems, etc.

I found several papers reviewing the subject, in spite of the lack of a global consensus.

Furthermore, I could only find a global dataset (presented latter) that had insuficient resolution and cocrdinate problems. Thus, I decided that the subject deserved further research, and though it would beuseful to make a document summarizing the major steps to reach it.

Whitfield, A., & Elliott, M. (2011). 1.07-Ecosystem and biotic classifications of estuaries and coasts. Treatise on estuarine and coastal science, 99-124.

In 2011 it was published in ‘Estuaries and Coasts’, Worldwide Typology of Nearshore Coastal Systems: Defining the Estuarine Filter of River Inputs to the Oceans, by Durr et al.
Fundamental to the definition of an estuary is the necessity to distinguish between an ‘estuary’ and ‘estuarine waters’, with the former requiring a physiographic form and the latter (e.g., an estuarine plume on an open coast) not constrained by this requirement.
The seaward boundary of an estuary can be defined geographically by a line between the land masses on each side of the entrance to the estuary from the coastal waters and where there is a notable break in the line of the coast. However, according to Ketchum (1983), this can be ambiguous where deltaic structures are formed at the mouth. The above-mentioned problem can be overcome by including geomorphological discontinuities (e.g., an ebb tidal delta) along the coast as part of the estuary.. It has also been suggested that the seaward limits of an estuary be defined by salinity distribution, but this view is widely contested.
A recent definition of an estuary using the salinity theme is provided by Elliott and McLusky (2002), who suggested that the estuaries are represented by “an area receiving freshwater inputs where the waters on a depth-averaged basis have a salinity of less than 95% of the adjacent local offshore seawater for 95% of the time.” This chapter favors the geomorphological criterion for the seaward boundary of an estuary, with plumes of estuarine water entering the sea being regarded as part of the TWs within the coastal marine environment rather than the estuary sensu stricto.
McHugh (1967) has defined these offshore estuarine zones of the world as those waters bounded by a surface salinity of 33.5 (Figure 6) and Frankignoulle et al. (1998) defined this zone as being ‘outer estuaries’ (as opposed to ‘inner estuaries’ that are semi-enclosed by land).

European Union (EU) Water Framework Directive (WFD; European Union, 2000)

Perhaps the value of a comprehensive system of estuarine definition and classification is epitomized by the European Union (EU) Water Framework Directive (WFD; European Union, 2000), which requires all surface waters to be classified, monitored, and compared against reference conditions as components of catchment management plans. The Directive includes estuaries in the term ‘transitional waters’ (TWs); hence, European estuarine scientists will have to use TWs as a means towards understanding and managing estuaries and their related habitats (Elliott and McLusky, 2002).

The EU definition considers that an estuary includes reedbeds, saltmarsh, and other intertidal habitats, and forms an ecological unit with the surrounding terrestrial coastal habitat types. Although this started as a convenient administrative term to encompass all waters that are neither the open coast nor true freshwaters, TWs have now been accepted as a well-defined set of ecosystems, each with its own set of characteristics different from the end members; and where all estuaries are TWs but not all TWs are estuaries.

Classification: CMECS Us standards - US

There are US standards to classify the areas, ‘coastal and marine classification standard’. These are very detailed gathering all information, and have recently been applied in the Gulf of Mexico.
In this document the areas are classified first into 3 aquatic types: Lacustrine, Estuarine and Marine, based on salinity and geomorphological characteristics of the setting. Then, within each aquatic system, there are subsystems, based on depth and position relative to the shoreline.
-Lacustrine Systems include the Subsystems Littoral and Limnetic, conforming, respectively to shoreline and deepwater habitats.
-The Estuarine System has four Subsystems: Coastal, Open Water, Tidal Riverine Coastal and Tidal Riverine Open Water.
-The Marine System is comprised of three Subsystems: Nearshore, Offshore, and Oceanic, distinguished by total water depth.
-Bay, Fjord, etc. are part of the geomorphological classification (tectonic settings, physiographic settings describe landscape level geomorphological features from the coast to mid-ocean spreading centers), and do not necessarily correspond to the common names (e.g. Blue Hill Bay is actually a Fjord).

The Coastal and Marine Ecological Classification Standard (CMECS) was adopted by the Federal Geographic Data Committee (FGDC) in 2012 to fill the gap in the U.S. classification system [3]. Other similar national classification systems exist throughout the world. For example, the European Nature Information System (EUNIS) assigns hierarchical classifications to the aquatic, terrestrial, and marine habitats within and adjacent to Europe [4]. The United Kingdom, Australia, and New Zealand also have national classification systems [5] (Kingon, 2018).

1. Global Estuaries Database (GED: UNEP + Sea Around Us, 2003)

The UNEP-WCMC provides the shape file of a global estuaries database, within the Ocean Data Viewer, that shows the spatial distribution of over 1300 estuaries, including some lagoon systems and fjords. This dataset was developed as part of a “Sea Around Us” project. download.
It does not include all estuaries (not exhaustive), but certainly the many of them. Different areas however, have different resolutions and more importantly, some coordinates appear to be deviated (perhaps they were done ignoring the coordinate system).

Alder (2003). Putting the coast in the “Sea Around Us”. The Sea Around Us Newsletter 15: 1-2.

This is Simply the code for the base layer:

# 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 = "gray80", colour = "gray70",) +
  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
GED <- sf::read_sf("D:/shapes/UBC003_SAU_Estuaries2003_v2/DownloadPack-14_001_UBC003_SAU_Estuaries2003_v2/01_Data/14_001_UBC003_SAU_Estuaries2003_v2.shp")
names(GED)
##  [1] "ID_NUMBER"  "LABEL"      "CODE"       "COUNTRY"    "CNUMBER"   
##  [6] "CONTINENT"  "DISCHARGE"  "RIVER_SYS"  "INPUT_LAT"  "INPUT_LON" 
## [11] "REF1"       "REF2"       "YEAR_START" "YEAR_END"   "YEAR_PUB"  
## [16] "SUB"        "Area_km2"   "geometry"
# estuary[1,] %>% data.frame()
# plot(estuary[2], main="UN estuaries data-base")
 
# see map
world_map.s +  
  geom_sf(data = GED, aes(fill = LABEL, col = LABEL))+
  ggtitle("UN-S global database")+
  theme(legend.position="none")+
  coord_sf(expand = FALSE)

# # interactive map of details
# GED %>%
#   # st_buffer(dist=.01) %>%
#   # sf::st_simplify() %>% 
#   st_crop(xmin=-10, xmax=10, ymin=36, ymax=53) %>% 
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  zcol = "LABEL",
#                  col.regions = brewer.pal(11, "Spectral"),
#                  alpha.regions = 0.3,
#                  burst = TRUE,
#                  legend = FALSE)

rm(GED)

2. World Estuary Delta (WED, MarineRegions

I also found this data base, which is excellent in terms of resolution, but only have the main estuaries. The ones found here, I replaced them in the UNEP file by these ones, which are much better. This file is not available in the marine regions site (due to the absence of the exact methology of each individual estuary, pers. com. L. Schepers), but got it from the Belgian geoserver.

MarineRegions:world_estuary_delta

WED <- sf::read_sf("D:/shapes/MarineRegions-world_estuary_delta.kml") %>% 
  st_transform(crs = 3857) %>% 
  st_buffer(dist=0.001) %>% #to avoid "Ring Self-intersection"
  st_transform(crs = 4326) 

world_map.s +
  geom_sf()+
  geom_sf(data=WED, col=2)+
  ggtitle("world_estuary_delta")

world_map.b +
    geom_sf(data=WED, col=2, fill=5, alpha=.5)+
    coord_sf(xlim = c(-10, 0), ylim = c(36,44))+
    ggtitle("Detail of world_estuary_delta")

# # interactive map of details
# WED %>%
#   # st_buffer(dist=.01) %>%
#   # sf::st_simplify() %>% 
#   st_crop(xmin=-10, xmax=10, ymin=36, ymax=53) %>% 
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  col.regions = brewer.pal(11, "Spectral"),
#                  alpha.regions = 0.3,
#                  burst = TRUE,
#                  legend = FALSE)

rm(WED)

3. European Transitional Waters (EuTW)

Data kindly provided by Pedro Morais (CCMAR).
A newer version should be found in this site, but I could not find it yet.

Note that the different areas have dieverse resolution and mostly, some areas are missing, for example, few areas were found arround the Mediterranean (from Marseille to Greece), for example.

“Transitional waters: Bodies of surface water in the vicinity of river mouths which are partly saline in character as a result of their proximity to coastal waters but which are substantially influenced by freshwater flows.”

EuTW <- sf::read_sf("D:/shapes/wise_transitional_waters_shp-20191011T113722Z-001/wise_transitional_waters_shp/Transitional_waters.shp") %>% st_transform(crs=4326)
names(EuTW)
##  [1] "CODE"       "ID"         "NAME"       "LGE_ID"     "LABEL"     
##  [6] "ENGLISH"    "CTY_ID"     "BSN_ID"     "UPDATED_WH" "UPDATED_BY"
## [11] "GEOBOX_LLX" "GEOBOX_LLY" "GEOBOX_URX" "GEOBOX_URY" "PFAFSTETTE"
## [16] "Shape_Leng" "Shape_Area" "geometry"
# plot(EUW)
(kk2 <- st_bbox(EuTW)) # map limits
##      xmin      ymin      xmax      ymax 
## -27.96824  35.83989  21.29961  57.89458
# Coastline
cl = ne_coastline(scale = 'large', returnclass = "sf")

# Plot
world_map.s +
    geom_sf(data=EuTW, aes(col=NAME))+
    theme(legend.position="none")+
    coord_sf(xlim = c(kk2["xmin"], kk2["xmax"]), ylim = c(kk2["ymin"], kk2["ymax"]))+
  ggtitle("European transitional waters")

# countries, large scale
pt <- ne_countries(country = c('portugal', 'spain'), scale = 'large', returnclass = "sf") %>% st_transform(crs=4326)
world_map.b +
    geom_sf(data=pt)+
    geom_sf(data=EuTW, aes(col=NAME, fill=NAME), alpha=.4)+
    theme(legend.position="none")+
    coord_sf(xlim = c(-10, -6), ylim = c(36,39))+
  ggtitle("European transitional waters (detail)")

# # Interactive map of Iberian details:
# EuTW %>%
#   st_crop(xmin=-10, xmax=10, ymin=36, ymax=53) %>% 
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  zcol = "LABEL",
#                  col.regions = brewer.pal(11, "Spectral"),
#                  alpha.regions = 0.5,
#                  burst = TRUE,
#                  legend = FALSE)

rm(EuTW)

4. New Zealand coastal hydrosystems (NZ)

Classification of New Zealand coastal hydrosystems

Dataset to download

# https://www.spatialanalytics.co.nz/post/2018/11/01/plotting-new-zealand-with-r/

NZ <- sf::read_sf("D:/shapes/NZ coastal hydrosystems/NZ_Coastal_Hydrosystems_poly.kml") %>% dplyr::select(Name)

# plot(nz_hydro)
kk2 <- st_bbox(NZ) # map limits
kk2["xmin"] = 165 #replace for easier vizualition 
  
# countries, large scale
nz <- ne_countries(country = 'new zealand', scale = 'large', returnclass = "sf") %>% st_transform(crs=4326)

world_map.b +
    geom_sf(data=nz)+
    geom_sf(data=NZ, aes(col=Name))+
    theme(legend.position="none")+
    coord_sf(xlim = c(165, kk2["xmax"]), ylim = c(kk2["ymin"], kk2["ymax"]))+
  ggtitle("New Zealand hydrosystems")

# world_map.b +
#     geom_sf(data=nz)+
#     geom_sf(data=NZ, aes(col=Name, fill=Name), alpha=.4)+
#     theme(legend.position="none")+
#     coord_sf(xlim = c(172.5, 175), ylim = c(-36, -34))+
#   ggtitle("New Zealand hydrosystems (detail)")


# # Interactive map of details
# NZ %>%
#   st_crop(kk2) %>%
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  zcol = "Name",
#                  col.regions = brewer.pal(11, "Spectral"),
#                  alpha.regions = 0.3,
#                  burst = TRUE,
#                  legend = FALSE)


rm(NZ)

5. Australia

Other links to explore in AU (thanks to OEH BIS Data Broker):

coastal adaptation to climate change.
https://www.ga.gov.au/scientific-topics/marine/coasts-estuaries from this last one I went to:
marine sediments (MARS)
marine spatial planning (AMSIS)
environmental layers
Australia’s coastal information- email sent.

Finally got the shape files from the Geomorphic habitat datasets on the OZcoast page [http://www.ozcoasts.gov.au/search_data/datasets.jsp] of the different areas.
Then I binded it together, although the files include different habitats (i.e. tidal flat, rocky shore, etc.)

-The ones I will use:
Central Basin: “Central basins are uniform, lower energy environments in the deeper and quieter parts of estuaries, and are often formed (…)”
Saltmarsh/Saltfalt (SM)-“A salt marsh is a marshy area found near estuaries and sounds. The water in salt marshes varies from completely saturated with salt to freshwater(…)”
Flood- and Ebb-tide Delta (FED)
Chanel (are tidal channels or fluvial flow river channels)
Mangroves (?)

The other ones will be excluded, as seen in the script.

# Find all shape files within the folder:
(kk1 <- list.files("D:/shapes/australia estuaries/", pattern = "*.shp", full.names=TRUE, recursive = TRUE))
## [1] "D:/shapes/australia estuaries//nsw_geohab_av/geohab_nsw_v2.shp"                                    
## [2] "D:/shapes/australia estuaries//nt_geohab_av/czm/pristine_ests/release/nt/shape/geohab_nt_v2.shp"   
## [3] "D:/shapes/australia estuaries//qld_geohab_av/czm/pristine_ests/release/qld/shape/geohab_qld_v2.shp"
## [4] "D:/shapes/australia estuaries//sa_geohab_av/shape/geohab.shp"                                      
## [5] "D:/shapes/australia estuaries//tas_geohab_av/czm/pristine_ests/release/tas/shape/geohab_tas.shp"   
## [6] "D:/shapes/australia estuaries//vic_geohab_av/czm/pristine_ests/release/vic/shape/geohab_vic.shp"   
## [7] "D:/shapes/australia estuaries//wa_geohab_av/czm/pristine_ests/release/wa/shape/geohab_wa_v2.shp"
# read the shape files and bind it together  
AU <- rbind(sf::read_sf(kk1[1]), 
        sf::read_sf(kk1[2]),
        sf::read_sf(kk1[3]),
        sf::read_sf(kk1[4]),
        sf::read_sf(kk1[5]),
        sf::read_sf(kk1[6]),
        sf::read_sf(kk1[7])
        )

# however, this file has several habitats and we only want saline gradients ones
# types of habitats:
unique(AU$GH_TYPE)
##  [1] "Channel"                   "Mangrove"                 
##  [3] "Flood- and Ebb-tide Delta" "Intertidal Flats"         
##  [5] "Barrier/back-barrier"      "Saltmarsh/Saltflat"       
##  [7] "Central Basin"             "Fluvial (bay-head) Delta" 
##  [9] "Unassigned"                "Rocky Reef"               
## [11] "Tidal Sand Banks"          "Bedrock"                  
## [13] "Coral"
# thus, we will exclude which are clearly out of scope 
# and change to a planar coordinate system 
# and simplify the shape:
AU <- AU %>% 
    dplyr::filter(!GH_TYPE %in% c("Rocky Reef", "Bedrock", "Coral", "Tidal Sand Banks","Intertidal Flats","Barrier/back-barrier","Unassigned")) %>% 
    st_transform(crs = 3857) %>%    
    st_simplify() %>% 
    st_transform(4326) %>% 
    rename(Name = "EST_NAME")

# # aggregate by area (thanks to Roger Bivand)
AU <- aggregate(AU, by = list(AU$Name, AU$GH_TYPE), head, n=1)

# Countries, large scale, map limits...
au <- ne_countries(country = 'australia', scale = 'large', returnclass = "sf")
kk2 <- st_bbox(AU) 

world_map.s +
    geom_sf(data=au)+
    geom_sf(data=AU, aes(col = GH_TYPE, fill = GH_TYPE), alpha=.5)+
    coord_sf(xlim = c(kk2["xmin"], kk2["xmax"]), ylim = c(kk2["ymin"], kk2["ymax"]))+
    theme(legend.position="none")+
    ggtitle("Australia's geomorphic habitats")

world_map.b +
    geom_sf(data=au)+
    geom_sf(data=AU, aes(col = GH_TYPE, fill = GH_TYPE), alpha=.5)+
    # theme(legend.position="none")+
    coord_sf(xlim = c(144, 148), ylim = c(-38, -40))+
  ggtitle("Australia's geomorphic habitats (detail)")

# # # interactive map of details
# AU %>%
#   st_crop(kk2) %>%
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  zcol = "GH_TYPE",
#                  col.regions = brewer.pal(11, "Spectral"),
#                  alpha.regions = 0.3,
#                  burst = TRUE,
#                  legend = FALSE)

## Export the kml to merge with the others:
# AU %>%
#   st_transform(4326) %>% # just because wgs84...
#   select(Description = Name) %>%
#   st_write("D:/shapes/australia estuaries/australia_geomorfic0.kml",
#            driver='kml',
#            fid_column_name="Name",
#            delete_dsn = TRUE)

rm(AU)

The reference for this dataset:
Heap A., Bryce S., Ryan D., Radke L., Smith C., Smith R., Harris P. and Heggie D. 2001 Australian Estuaries & Coastal Waterways: A Geoscience Perspective for Improved and Integrated Resource Management A Report to the National Land & Water Resources Audit Theme 7: Ecosystem Health AGSO Record 2001/07 Australian Geological Survey Organisation Commonwealth of Australia, Canberra (Available from the OzEstuaries website

6. United States

EDM program form EPA

Q: Where does EDM get its estuary polygons?
A: We started with the estuary and coastal system polygons created for the Environmental Monitoring & Assessment Program’s National Coastal Condition Assessment sampling. We edited that collection to identify systems that were truly estuaries based on both freshwater and saltwater connections and partial closure on the seaward boundary. For example, the National Coastal Assessment sample frame includes open systems that are not true estuaries as they are too open to the ocean. EDM’s estuarine polygon set includes all systems for which a watershed can be constructed based on an aggregation of NHDPlus EXIT catchments. EDM’s estuarine polygon set was circulated among EPA regions for review and input. In EDM, individual estuaries are identified by the ESTCODE and Estuary Name and a SUBCODE and Subestuary Name when subsystems could be identified.
I downloaded the program, and extracted each chunk of the coast separately, and binded it together.

Note also that sound systems were removed from the file, which appear to have no salinity changes (“A sound is often formed by the seas flooding a river valley. This produces a long inlet where the sloping valley hillsides descend to sea-level and continue beneath the water to form a sloping sea floor. The Marlborough Sounds in New Zealand are good examples of this type of formation.”).

# Find all shape files within the folder:
(kk1 <- dir("D:/shapes/EPA_gov US estuaries/", pattern = "*.shp", full.names=TRUE, recursive = TRUE))
## [1] "D:/shapes/EPA_gov US estuaries//atlantic/estuaries.shp"
## [2] "D:/shapes/EPA_gov US estuaries//gulf/estuaries.shp"    
## [3] "D:/shapes/EPA_gov US estuaries//missi/estuaries.shp"   
## [4] "D:/shapes/EPA_gov US estuaries//pacific/estuaries.shp"
# read the shape files and bind it together  (epsg (SRID):4269)
US <- rbind(sf::read_sf(kk1[1]) %>% st_transform(4326), 
        sf::read_sf(kk1[2]) %>% st_set_crs(4269) %>% st_transform(4326),
        sf::read_sf(kk1[3]) %>% st_transform(4326),
        sf::read_sf(kk1[4]) %>% st_transform(4326)
        )
head(US)
## Simple feature collection with 6 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -67.1755 ymin: 44.81625 xmax: -66.96859 ymax: 44.96975
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 10
##   ESTCODE STATE ESTUARY NCA_NAME SUBCODE SUBEMBAYMT PROVINCE AREA_SQKM
##   <chr>   <chr> <chr>   <chr>    <chr>   <chr>      <chr>        <dbl>
## 1 COBS    ME    Cobsco~ Western~ <NA>    <NA>       Acadian~      7.43
## 2 COBS    ME    Cobsco~ East Bay <NA>    <NA>       Acadian~      5.55
## 3 COBS    ME    Cobsco~ Friar R~ <NA>    <NA>       Acadian~      6.39
## 4 COBS    ME    Cobsco~ Cobscoo~ <NA>    <NA>       Acadian~     43.3 
## 5 COBS    ME    Cobsco~ Pennama~ <NA>    <NA>       Acadian~      4.28
## 6 COBS    ME    Cobsco~ South B~ <NA>    <NA>       Acadian~     15.9 
## # ... with 2 more variables: HECTARES <dbl>, geometry <MULTIPOLYGON [°]>
# exclude a clearly rive chunk:
US <- US[US$AREA_SQKM!=273.874,]

# what names to use :
# US %>% select(ESTUARY, NCA_NAME) %>% data.frame() %>% head(15)
US <- US[!US$NCA_NAME %in% US$NCA_NAME[grep("sound",US$NCA_NAME, ignore.case=TRUE)],]

# lets replace the 'NA' ESTUARY code, by the NCA_NAME
# the remaining NCA_NAME are actually chunks of the estuaries.
US <- US %>% 
  mutate(Name = ifelse(is.na(ESTUARY), NCA_NAME, ESTUARY)) %>%   #group_by(Name) %>% 
  #summarize(st_union(geometry)) %>% 
  dplyr::mutate(Source="US") %>% # add a source col
  dplyr::select(Name, Source, ESTUARY) 

# Make a variable with type of system:
# table(sub('^.* ([[:alnum:]]+)$', '\\1', US$Name))
US$Type <- NA
US$Type[grep("River", US$Name)] = "River"
US$Type[grep("Bays$", US$Name)] = "Bays"
US$Type[grep("Bay$", US$Name)] = "Bay"
US$Type[grep("Sound$", US$Name)] = "Sound"
US$Type[grep("Sounds$", US$Name)] = "Sound"
US$Type[grep("Creek", US$Name)] = "Creek"
US$Type[grep("Lagoon$", US$Name)] = "Lagoon"
US$Type[grep("Lake", US$Name)] = "Lake"
US$Type[grep("Lakes", US$Name)] = "Lake"
US$Type[grep("Inlet$", US$Name)] = "Inlet"
US$Type[grep("River Estuary$", US$Name)] = "River Estuary"
US$Type[grep("Cove$", US$Name)] = "Cove"
US$Type[grep("Canal$", US$Name)] = "Canal"
US$Type[grep("Harbor$", US$Name)] = "Harbor"
# sort(table(US$Name[is.na(US$Type)]))

# # aggregate by area (thanks to Roger Bivand)
# US <- aggregate(US, by = list(US$Name,US$ESTUARY, US$Type), head, n=1)

# Plot:
cl = ne_coastline(scale = 'large', returnclass = "sf")
kk2 <- st_bbox(US) # map limits

world_map.s +
    geom_sf(data=US, aes(col = Type, fill=Type))+
    #theme(legend.position="none")+
    coord_sf(xlim = c(kk2["xmin"], kk2["xmax"]), ylim = c(kk2["ymin"], kk2["ymax"]))+
  ggtitle("US's estuaries")

world_map.b +
    geom_sf(data=cl)+
    geom_sf(data=US, aes(col = Type, fill=Type), alpha=.4)+
    # theme(legend.position="none")+
    coord_sf(xlim = c(-80, -74), ylim = c(34, 40))+
  ggtitle("US's estuaries (detail)")

# # # # For an interactive map of details (not working with rmarkdown)
# US %>%
#   st_crop(kk2) %>%
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  zcol = "Name",
#                  col.regions = brewer.pal(11, "Spectral"),
#                  alpha.regions = 0.3,
#                  burst = TRUE,
#                  legend = FALSE)

## Export the kml to merge with the others:
# US %>%
#    st_transform(4326) %>% # just because wgs84...
#    select(Description = Name, Type, ESTUARY, Source) %>%
#    st_write("D:/shapes/EPA_gov US estuaries/US_estuaries.kml",
#             driver='kml',
#             fid_column_name="Name",
#             delete_dsn = TRUE)

rm(US)

7. Mexico ‘lindo y querido’

I did not found yet a good compilation, but ‘Humedales de México’ has the kml and shapefile of 142 RAMSAR sites in Mexico. We will use this resource, downloading each shape and merging into a kml file for posterior merge in the large database.

See also the (book)[http://201.116.60.25/publicaciones/AAM_2016.pdf]: “El estudio “Humedales de la República Mexicana” (2012) generó el Inventario Nacional de Humedales (INH), que incluye 6 331 humedales y complejos de humedales, cubriendo un 5% de la supericie del país (tabla 4.9). Los humedales están clasiicados en palustres (relacionados a lagunas o pantanos), lacustres (lagos), luviales (ríos), estuarinos (estuarios) y creados por la acción antropogénica.”

MX <- sf::read_sf("D:/shapes/RAMSAR142_Mexico_geoITRF92_Ene15.kml")
head(MX)
## Simple feature collection with 6 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XYZ
## bbox:           xmin: -116.6608 ymin: 21.7162 xmax: -102.2948 ymax: 32.71865
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 3
##   Name             Description                                     geometry
##   <chr>            <chr>                                 <MULTIPOLYGON [°]>
## 1 Estero El Solda~ "<html xmlns:fo=\"htt~ Z (((-110.9697 27.96954 0, -110.~
## 2 (El Jagüey) Bue~ "<html xmlns:fo=\"htt~ Z (((-102.3024 21.7174 0, -102.3~
## 3 Sistema de Hume~ "<html xmlns:fo=\"htt~ Z (((-114.7884 32.0299 0, -114.7~
## 4 Humedales del D~ "<html xmlns:fo=\"htt~ Z (((-115.1631 32.15339 0, -115.~
## 5 Agua Dulce       "<html xmlns:fo=\"htt~ Z (((-113.0211 31.92458 0, -113.~
## 6 Estero de Punta~ "<html xmlns:fo=\"htt~ Z (((-116.6021 31.766 0, -116.59~
# Plot:
cl = ne_coastline(scale = 'large', returnclass = "sf")
kk2 <- st_bbox(MX) # map limits

world_map.s +
    geom_sf(data=MX, aes(col = Name, fill=Name))+
    theme(legend.position="none")+
    coord_sf(xlim = c(kk2["xmin"], kk2["xmax"]), ylim = c(kk2["ymin"], kk2["ymax"]))+
  ggtitle("Humedales de Mexico (selected)")

# # # # Interactive map of details
# MX %>%
#   st_crop(kk2) %>%
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  zcol = "Name",
#                  col.regions = brewer.pal(11, "Spectral"),
#                  alpha.regions = 0.3,
#                  burst = TRUE,
#                  legend = FALSE)

8. Canada

Estuaries data from Canada has been hard to get. Just managed to get a chunk of the coast from Michelle Greenlaw.

CA <- sf::read_sf("D:/shapes/BaysNew/BaysNew.shp") %>% 
  st_transform(4326) %>% 
  rename(Name = "NameText")
head(CA)
## Simple feature collection with 6 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -65.61818 ymin: 43.40703 xmax: -61.31864 ymax: 45.1898
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 5
##   OBJECTID Name    Shape_Leng Shape_Area                           geometry
##      <int> <chr>        <dbl>      <dbl>                 <MULTIPOLYGON [°]>
## 1        1 Necum ~     41711.  15983467. (((-62.19349 44.95147, -62.19354 ~
## 2        2 Indian~     38302.  16843336. (((-61.84586 45.10788, -61.8456 4~
## 3        3 Aspoto~     23137.  11324948. (((-64.00654 44.50663, -64.00646 ~
## 4        4 Barrin~     99183.  86091895. (((-65.57807 43.56762, -65.57809 ~
## 5        5 Beaver~     28238.  14044547. (((-62.3799 44.90652, -62.37968 4~
## 6        6 Before~     12400.   2109668. (((-61.3395 45.18967, -61.33951 4~
# Plot:
kk2 <- st_bbox(CA) # map limits

world_map.b +
    geom_sf(data=CA, aes(col = Name, fill=Name))+
    theme(legend.position="none")+
    coord_sf(xlim = c(kk2["xmin"], kk2["xmax"]), ylim = c(kk2["ymin"], kk2["ymax"]))+
  ggtitle("Bays and Estuaries from Nova Scotia (Canada)")

# # # # Interactive map of details
# CA %>%
#   st_crop(kk2) %>%
# mapview::mapview(map.types = "Esri.OceanBasemap",
#                  zcol = "Name",
#                  col.regions = brewer.pal(11, "Spectral"),
#                  alpha.regions = 0.3,
#                  burst = TRUE,
#                  legend = FALSE)
rm(CA)

7. Individual estuaries

Musquash Estuary Marine Proctected Area

plot(sf::read_sf("D:/shapes/World-musquash.kml"))


# NOTE that here we should get the estuarine seas: Baltic and White sea, kml obtained from marineregions.org

Final map: merge all

Now, I reopen (yep, silly but easier for users) all partial files mentioned above in r and select just larger areas (areas larger than Mira estuary which is ~10 km^2.), simplify the shapes and perhaps merge it when they overlap, so that finally obtaining a polygon for all world estuaries.

## [1] 464   3
## [1] 458   3
## [1] 220
## [1] "Adur"        "Swansea Bay" "Ore"         "Ain"         "Alt"        
## [6] "Odeloura"
## [1] 90  3
## [1] 90  3
## [1] 421   3
##  [1] "Buller"           "Clutha"           "Hokianga Harbour"
##  [4] "Kaipara Harbour"  "Lake Ellesmere"   "Manawatu"        
##  [7] "Manukau Harbour"  "Cleddau"          "Hauraki Gulf"    
## [10] "Oreti"            "Bays of Islands"  "Pelorus Sound"   
## [13] "Hauraki Gulf"     "Awanui"           "Bays of Islands" 
## [16] "Waituna Wetlands" "Waiau"            "Hauraki Gulf"    
## [19] "Port Waikato"     "Hauraki Gulf"     "Waitaki"         
## [22] "Waituna Wetlands" "Wanganui"         "Hutt"
## [1] 100
## Simple feature collection with 6 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -65.61818 ymin: 43.40703 xmax: -61.31864 ymax: 45.1898
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 3
##   Name         Source                                              geometry
##   <chr>        <chr>                                     <MULTIPOLYGON [°]>
## 1 Necum Teuch~ CA     (((-62.19349 44.95147, -62.19354 44.95158, -62.1937 ~
## 2 Indian Harb~ CA     (((-61.84586 45.10788, -61.8456 45.10783, -61.84545 ~
## 3 Aspotogan B~ CA     (((-64.00654 44.50663, -64.00646 44.50667, -64.00647~
## 4 Barrington ~ CA     (((-65.57807 43.56762, -65.57809 43.56758, -65.57811~
## 5 Beaver Harb~ CA     (((-62.3799 44.90652, -62.37968 44.90648, -62.3793 4~
## 6 Before Tor ~ CA     (((-61.3395 45.18967, -61.33951 45.1895, -61.33945 4~
## [1] 149
## [1] "Purari"                  "Yukon"                  
## [3] "Canche"                  "Baie de Somme"          
## [5] "Estuaire de Seine aval"  "Estuaire de Seine moyen"

OTHERS/TO-DO

  • search Mediterranean and Canada estuaries.

  • add Musquash estuary shape (Bay of Foundy):
    descrip: Musquash Marine Protected Area, zone / outer boundary source: Oceans and Coastal Management Division, DFO Maritimes Region last_modif: 20061201 contact: Kristian Curran, Oceans and Coastal Management Division, Bedford Institute of Oceanography.

doc on global databases