Spatial Data & Informative Map: Part 2
May 1, 2022
1 Clear Workspace
rm(list = ls())2 Load library
library(ggplot2) # ggplot() fortify()
library(dplyr) # %>% select() filter() bind_rows()
library(rgdal) # readOGR() spTransform()
library(raster) # intersect()
library(ggsn) # north2() scalebar()
library(rworldmap) # getMap()3 Data input
GBIF data for the Brown Trout (Ca hoi nau)
brown_trout <- read.csv("Brown_Trout_GBIF_clip.csv")4 Plot ggplot
(trout_check <- ggplot(brown_trout, mapping = aes(x = decimallongitude, y = decimallatitude)) +
geom_point(alpha = 0.5))
# add colour
(trout_check <- ggplot(brown_trout, mapping = aes(x = decimallongitude, y = decimallatitude,
colour = scientificname)) +
geom_point())
5 Preliminary map & ggplot
We can roughly see the outline of Scandinavia and maybe the Northern Mediterranean if you squint.
To plot a preliminary map, crop the world map provided by the rworldmap package using:
# Get the world map
world <- getMap(resolution = "low") # SpatialPolygonsDataFrameclipper_europe <- as(extent(-10, 32, 30, 72), "SpatialPolygons")
proj4string(clipper_europe) <- CRS(proj4string(world))
world_clip <- raster::intersect(world, clipper_europe)
world_clip_f <- fortify(world_clip)Note: The first line uses extent() to make a SpatialPolygons object which defines a bounding box inside which to crop the world map polygons. The arguments in extent() are: extent(min_longitude, max_longitude, min_latitude, max_latitude). The second line changes the coordinate reference systems of both the counding box and the world map to be equal. The third line uses intersect() to clip world by the area of the bounding box, clipper_europe. The fourth line converts the SpatialPolygonsDataFrame to a normal flat dataframe for use in ggplot()
Next, we can map tiles using geom_polygon:
(trout_map <- ggplot() +
geom_polygon(data = world_clip_f, # world map, polygons
aes(x = long, y = lat, group = group),
fill = NA, colour = "black") +
geom_point(colour = "blue", alpha = 0.5,
aes(x = decimallongitude, y = decimallatitude),
data = brown_trout) + # brown_trout data
theme_bw() +
xlab("Longitude") +
ylab("Latitude") +
coord_quickmap())
Note for next step: The country outlines work well, but to tell us more about the habitat the Brown Trout lives in we can also plot the ecoregions data on the map.
6 Shapefile
6.1 Basic infor on Shapfile
Shapefiles are a data format developed by ESRI used to hold information on spatial objects. They are pretty ubiquitous and can be used by a lot of GIS packages Shapefiles can hold polygon, line or point data.
Despite the name, a shapefile consists of a few different files:
Mandatory files: .shp = The main file containing the geometry data .shx = An index file .dbf = An attribute file holding information on each object
Common additional files: .prj = A file containing information on the Coordinate Reference system .shp.xml = a file containing object metadata, citations for data, etc.
6.2 Read Shapefile
Use readOGR(), which converts a shapefile into a SpatialPolygons object.
shpdata_FEOW <- readOGR(dsn = "/Users/admin/Desktop/Linh Data Studio/Spatial Data/Spatial Modelling/FEOW-TNC",
layer = "FEOWv1_TNC")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/admin/Desktop/Linh Data Studio/Spatial Data/Spatial Modelling/FEOW-TNC", layer: "FEOWv1_TNC"
## with 449 features
## It has 6 fields
## Integer64 fields read as strings: ECO_IDNote: dsn = “FEOW-TNC” gives the name of the folder where the shapefile can be found; layer = “FEOWv1_TNC” gives the name of the files to read in
6.3 Check Co-ordinate Reference System (CRS)
proj4string(shpdata_FEOW)## [1] NA6.4 Transform CRS
To transform the CRS, we can use spTransform and specify the correct CRS, which in this case is EPSG:WGS84 (+proj=longlat +datum=WGS84). WGS84 is normallly used to display large maps of the world.:
# shpdata_FEOW <- spTransform(shpdata_FEOW, CRS("+proj=longlat +datum=WGS84")) doesn't work, need to add if (FALSE) {}
if (FALSE) {
shpdata_FEOW <- spTransform(shpdata_FEOW, CRS("+proj=longlat +datum=WGS84"))
}6.5 Crop data to Europe
The shapefile contains ecoregions for the entire world, but we only want to plot the ecoregions where the brown trout is found. shpdata_FEOW is a SpatialPolygonsDataFrame, so we can use the same method as before to crop the object to the extent of a bounding box, using intersect():
shpdata_FEOW_clip <- raster::intersect(shpdata_FEOW, clipper_europe) # only choose a specific areas in Euro where brown trout live plot(shpdata_FEOW_clip)
6.6 ggplot & world map
ECOREGION contains all the data for the different types of ecoregions, they have names like “Aegean Drainages” and “Central Prairie”. Now we can use ECOREGION as an identifier in the fortify() command to transform the spatial object to a dataframe, where each polygon will be given an id of which ECOREGION it is from.
Transforming to a dataframe is necessary to keep the metadata for each polygon, which is necessary to colour the points in the ggplot.
shpdata_FEOW_clip_f <- fortify(shpdata_FEOW_clip, region = "ECOREGION")(map_FEOW <- ggplot() +
geom_polygon(data = shpdata_FEOW_clip_f,
aes(x = long, y = lat, group = group, fill = id),
color = "black", size = 0.5) +
geom_point(colour = "red", alpha = 0.5, size = 0.5,
aes(x = decimallongitude, y = decimallatitude),
data = brown_trout) +
theme_classic() +
theme(legend.position="bottom") +
theme(legend.title = element_text(size = 4),
legend.text = element_text(size = 4)) +
xlab("Longitude") +
ylab("Latitude") +
coord_quickmap())
6.7 Add other element 1
We want to indicate a potential area for a trout re-introduction program:
Finland and Estonia have hardly any trout, but would probably have the right conditions according to the ecoregions:
(map_FEOW_annot <- map_FEOW +
annotate("rect", # Adding a Shaded Rectangle
xmin = 20 , xmax = 35, ymin = 55, ymax = 65, # Finland and Estonia location
fill="yellow", alpha=0.5) +
annotate("text", x = 27.5, y = 61, size = 3, label = "Restock\nArea")) # size of text 
6.8 Add other element 2
To further explore which ecoregions would be suitable for a trout re-introduction program, you can also check which trout records fall within which ecoregion polygons using the rgdal package.
6.8.1 Create a SpatialPoints object
First, create a SpatialPoints object from the Brown Trout records:
brown_trout_sp <- SpatialPoints(
coords = data.frame(brown_trout$decimallongitude, brown_trout$decimallatitude),
proj4string = CRS(proj4string(shpdata_FEOW_clip)))coords = uses the coordinates from the brown trout dataset formatted as a dataframe. the CRS (proj4string) is set to be the CRS of the ecoregions spatial object.
6.8.2 Create a dataframe from Sp objective
over() from the sp package (loaded by default by many other R spatial analysis packages) then creates a dataframe with the same number of rows as brown_trout_sp, where each row contains the data of the polygon in which the data point is found.
point_match <- over(
brown_trout_sp, # sparial point
shpdata_FEOW_clip) # spatial polygon 6.8.3 Create summary table count rows
Create a summary table counting the number of rows grouped by ECOREGION.
point_match %>%
group_by(ECOREGION) %>%
tally() %>%
arrange(desc(n))## # A tibble: 15 × 2
## ECOREGION n
## <chr> <int>
## 1 Northern Baltic Drainages 43840
## 2 Norwegian Sea Drainages 28285
## 3 <NA> 14455
## 4 Eastern Iberia 11892
## 5 Barents Sea Drainages 3616
## 6 Cantabric Coast - Languedoc 3054
## 7 Gulf of Venice Drainages 1021
## 8 Upper Danube 547
## 9 Western Iberia 435
## 10 Central & Western Europe 267
## 11 Dniester - Lower Danube 162
## 12 Northern British Isles 112
## 13 Dalmatia 31
## 14 Italian Peninsula & Islands 8
## 15 Atlantic Northwest Africa 1point_match %>%
group_by(ECOREGION) %>%
tally() %>%
arrange(desc(n))## # A tibble: 15 × 2
## ECOREGION n
## <chr> <int>
## 1 Northern Baltic Drainages 43840
## 2 Norwegian Sea Drainages 28285
## 3 <NA> 14455
## 4 Eastern Iberia 11892
## 5 Barents Sea Drainages 3616
## 6 Cantabric Coast - Languedoc 3054
## 7 Gulf of Venice Drainages 1021
## 8 Upper Danube 547
## 9 Western Iberia 435
## 10 Central & Western Europe 267
## 11 Dniester - Lower Danube 162
## 12 Northern British Isles 112
## 13 Dalmatia 31
## 14 Italian Peninsula & Islands 8
## 15 Atlantic Northwest Africa 16.8.4 Map
Make a map with scale bar and a north arrow
map_FEOW_scale <- map_FEOW_annot +
scalebar(data = shpdata_FEOW_clip_f,
transform = TRUE, dist = 500, dist_unit = "km", model='WGS84',
height = 0.01, location = "bottomright", anchor = c(x = 25, y = 32))
north2(map_FEOW_scale, x = 0.3, y = 0.85, scale = 0.1, symbol = 1)