Common solutions for creating maps usually involves geographic information system (GIS) software, such as ArcGIS, QGIS, eSpatial, etc., which allow to visually prepare a map in a very similar approach as one would prepare a poster or a document layout. On the other hand, R has also advanced spatial capabilities and can be used to draw maps. As we mentioned earlier for the graphics, R allows you to create map highly customization which makes it a solution increasingly attractive. In addition, the leaflet package now offers interesting solution to add some interactivity to your maps.
Several packages will be used today. It is generally a good practice to group the package you will use at the beginning of your script:
-library(maptools)
-library(rgdal)
-library (sp)
-library(raster)
-library(ggplot2)
-library(sf)
-library(rnaturalearth)
-library(rnaturalearthdata)
-library(ggspatial)
-library(rgbif)
-library(mapr) #devtools::install_github("cran/mapr")
-library (GISTools) #devtools::install_github("cran/GISTools")
-library(marmap)
-library(leaflet)
-
-#Sys.setlocale("LC_TIME", "English")
-#Sys.setlocale(category = "LC_ALL", "Chinese (Traditional)_Taiwan.950")
-
-# IMPORTANT NOTE AS 2023/10: MANY PACKAGES ARE MIGRATING TO TERRA PACKAGE. YOU CAN CONSULT DOCUMENTATION ONLINE TO EDIT FUNCTIONS, BUT ALSO ASK CHATGPT COULD BE RELEVANTBasic map
-Getting a simple map is pretty straightforward in R using the packagemaptools
-
-
- Simplified world country polygons -
data(wrld_simpl) -
-
- “Plot the world” -
plot(wrld_simpl)
We have now a couple of options to modify a plot. Therefore you can customize this simple map with a some arguments mentioned earlier:
- -
Try to zoom on Taiwan. Quickly, you will see obvious limit in this simple map. The resolution is too bad at at fine scale. See output of str(wrld_simpl) and ?wrld_simpl.
Shapes
-gpx track
--
-
- reading -
If you like running, hiking, walking, or just being lazy in a park - many devices can now export gpx format. You can read the file directly as lines (i.e, tracks), points (i.e, track_points), and a few other formats. Look at the help on readOGR function form the rgdal package: Bindings for the ‘Geospatial’ Data Abstraction Library. Note that the same could be accomplished using the sf package and using the st_readfunction.
OGR data source with driver: GPX
-Source: "D:\Dropbox\NTU\IONTU\Class_lectures\OCEAN_5098_R\OCEAN5098\df4810487c8a45838b57b12f305a2529def7225f\Data\run.gpx", layer: "tracks"
-with 1 features
-It has 12 fieldsplot(run1, main='Line') # my running activity
-run2 <- readOGR(dsn="Data/run.gpx",layer="track_points")OGR data source with driver: GPX
-Source: "D:\Dropbox\NTU\IONTU\Class_lectures\OCEAN_5098_R\OCEAN5098\df4810487c8a45838b57b12f305a2529def7225f\Data\run.gpx", layer: "track_points"
-with 931 features
-It has 26 fieldsplot(run2, main='Points')
-dev.off()null device
- 1 Note the use of par(mfrow=c(1,2)) to put the two plots next to each others.
Those are spatial shapes. The line connects points recorded at given intervals. Resolution of the track will depend on the resolution of the records and your speed (being lazy actually increases the accuracy).
--
-
- writing -
The reverse can obviously be done. The function writeOGR exports objects as SpatialPointsDataFrame, SpatialLinesDataFrame, or SpatialPolygonsDataFrame (as defined in the sp package). This spatial information is critical because you can imagine that many projections exist for those spatial objects. writeOGR will write a shapefile (ArcGIS compatible), but many other formats are available by specifying the correct driver. Let’s look at how to export spatial data from the simple map created earlier.
writeOGR(wrld_simpl, dsn="Data", layer = "world_test", driver = "ESRI Shapefile", overwrite_layer = TRUE)
-# ESRI Shapefile is defaultshp files
-Now we could open world_test.shp in ArcGIS (or others), but we can also import shapefile (.shp) back into R. Let’s read the .shp file previously created.
world_shp <- readOGR(dsn = "Data",layer = "world_test")OGR data source with driver: ESRI Shapefile
-Source: "D:\Dropbox\NTU\IONTU\Class_lectures\OCEAN_5098_R\OCEAN5098\df4810487c8a45838b57b12f305a2529def7225f\Data", layer: "world_test"
-with 246 features
-It has 11 fieldsplot(world_shp)
Spatial data
-Creating spatial data from scratch in R seems a little bit confusing at first, but it is important to understand what is a spatial objects and how it works.
-As we saw before when making a plot, you can add many layers (vector-based objects) on the plot created: points points, lines lines, or also polygons) to your map. A transformation of numerical x and y coordinates must however be applied in the case of a map in order to recognize those values as spatial data.
-
-
SpatialPointsDataFramefor plotting points
-
plot(wrld_simpl,xlim=c(115,128) ,ylim=c(19.5,27.5),col='#D2B48C',bg='lightblue') # TW map
-coords <- matrix(c(121.537290,120.265541, 25.021335, 22.626524),ncol=2) # NTU and SYS univ.
-coords <- coordinates(coords) # assign values as spatial coordinates
-spoints <- SpatialPoints(coords) # create SpatialPoints
-df <- data.frame(location=c("NTU","SYS")) # create a dataframe
-spointsdf <- SpatialPointsDataFrame(spoints,df) # create a SpatialPointsDataFrame
-plot(spointsdf,add=T,col=c('black','black'),pch=19,cex=2.2) # plot it on our map
-text(121,24, 'TAIWAN', cex=1)
--Note the use of the argument
-add=Tin the second call of the functionplot. It allows adding an element to a plot already opened.
-
-
SpatialLinesDataFramefor plotting lines: let’s move to Canada where the provinces where drw with a ruler: an example with the province of Saskatchewan and following the same logic as ealier.
-
plot(wrld_simpl,xlim=c(-130,-60),ylim=c(45,80),col='#D2B48C',bg='lightblue')
-coords <- matrix(c(-110,-102,-102,-110,-110,60,60,49,49,60),ncol=2)
-l <- Line(coords)
-ls <- Lines(list(l),ID="1")
-sls <- SpatialLines(list(ls))
-df <- data.frame(province="Saskatchewan")
-sldf <- SpatialLinesDataFrame(sls,df)
-plot(sldf,add=T,col='#3d2402', cex=2)
-text(-114, 55, 'Saskatchewan', srt=90, cex=0.5)
-text(-114, 63, 'CANADA', cex=1)
-
-
SpatialPolygonesDataFramefor plotting polygons instead of lines. It will help in filling the spatial shpae with a selected color or pattern.
-
plot(wrld_simpl,xlim=c(-130,-60),ylim=c(45,80),col='#D2B48C',bg='lightblue')
-coords <- matrix(c(-110,-102,-102,-110,-110,60,60,49,49,60),ncol=2)
-p <- Polygon(coords)
-ps <- Polygons(list(p),ID="1")
-sps <- SpatialPolygons(list(ps))
-df <- data.frame(province="Saskatchewan")
-spdf <- SpatialPolygonsDataFrame(sps,df)
-plot(spdf,add=T,col='#45220d')
-text(-114, 55, 'Saskatchewan', srt=90, cex=0.7)
-text(-114, 63, 'CANADA', cex=1)
-text(-103, 46, 'UNITED STATES', cex=1)
-text(-40, 78, 'GREENLAND', cex=1)
-text(-35, 55, 'Atlantic Ocean', cex=1, col='#071238')
raster
-The package raster provides useful function to play with spatial: download, unzip, and import spatial shape (don’t forget to set up your working directory). Those spatial data are high resolution spatial data at the country level. Have a look at What is raster data? from the ArcGIS website.
GADM
-The function getData (rasterpackage) can download directly polygons (vectorized) shape from GADM (Global Administrative Areas, or other sources). You will simply need a country code such as TWN for Taiwan, JPN for Japan, etc. In addition, the argument level indicates the level you are trying to get. For provinces (first level subdivision) level=1, for the overall country level level=0. Check ?getData.
-Be reasonable at first, you may experience some delays because of slow connection -
-TWN1 <- getData('GADM', country="TWN", level=0) # data Taiwan
-JPN <- getData('GADM', country="JPN", level=0) # data Japan
-class(TWN1) # those datasets are SpatialPolygonsDataFrame[1] "SpatialPolygonsDataFrame"
-attr(,"package")
-[1] "sp"par(mfrow = c(1, 2))
-plot(TWN1,axes=T,bg=colors()[431],col='grey')
-plot(JPN,axes=T,bg=colors()[431],col='grey')
Don’t forget to close your graphical window:
-dev.off()null device
- 1 -
-
- zoom on a point a map -
Simply set up xlim and ylim
-
-
- Level -
TWN2 <- getData('GADM', country="TWN", level=1)
-TWN2$NAME_1[1] "Fujian" "Kaohsiung" "New Taipei" "Taichung" "Tainan"
-[6] "Taipei" "Taiwan" You can see the list of counties is not complete. It is often the case when using GADM that some ‘details’ are missing. Let’s visualize one county for which shape data are available.
-plot(TWN1,col="grey",xlim=c(119,122.5), ylim=c(21.5,25.5), bg=colors()[431], axes=T)
-KAO <- TWN2[TWN2$NAME_1=="Kaohsiung",]
-plot(KAO,col="grey 33",add=TRUE)
Note again the use of add=TRUE in the latest plot function.
-
-
- A map can further be customized by combining all seen earlier (e.g. adding population) -
# base map
-plot(TWN1,col="grey",xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431], axes=T)
-# adding spatial polygones
-TAI <- TWN2[TWN2$NAME_1=="Taipei" | TWN2$NAME_1=="New Taipei" ,]
-plot(TAI,col="black",add=TRUE)
-# adding spatial points
-coords <- matrix(cbind(lon=c(121.2,121.55,121.8),lat=c(25,25.19,24.5)),ncol=2)
-coords <- coordinates(coords)
-spoints <- SpatialPoints(coords)
-df <- data.frame(location=c("City 1","City 2","City 3"),pop=c(138644,390095,34562))
-spointsdf <- SpatialPointsDataFrame(spoints,df)
-scalefactor <- sqrt(spointsdf$pop)/sqrt(max(spointsdf$pop))
-plot(spointsdf,add=TRUE,col='white',pch=1,cex=scalefactor*3,lwd=2)
-# adding a location of NTU (not spatial point here)
-points(121.537290,25.021335, type="p", pch=18, col='white', cex=1.5)
-# adding text
-text(121.53,24.921335,"NTU", col='white', font=2)
-# adding scale
-maps::map.scale(x=120, y=25.4)
-# adding north arrow
-GISTools::north.arrow(xb=120.3,yb=24.7, len=0.06, lab='N')
Note the use of the synthax maps::map.scale. It means you wanna use the function map.scale from the package maps. It avoid passing by library(maps), thus avoiding possible conflict among functions.
Country data
-Each country will usually have open data platform with country shape data. Those data are usually more accurate and should be preferred. In Taiwan, shapefile of counties are available for download here. You could get even get a finer resolution at the township level if you wish (large file).
-First, prepare your system to receive traditional Chinese characters if like mine your system is not set to accept them. Actually, those steps have no relationship with making a map however to deal with special characters is a common issue you may experience when “mapping around”.
-Sys.setlocale(category = "LC_ALL", "Chinese (Traditional)_Taiwan.950")[1] "LC_COLLATE=Chinese (Traditional)_Taiwan.950;LC_CTYPE=Chinese (Traditional)_Taiwan.950;LC_MONETARY=Chinese (Traditional)_Taiwan.950;LC_NUMERIC=C;LC_TIME=Chinese (Traditional)_Taiwan.950"sessionInfo()R version 4.2.2 (2022-10-31 ucrt)
-Platform: x86_64-w64-mingw32/x64 (64-bit)
-Running under: Windows 10 x64 (build 19045)
-
-Matrix products: default
-
-locale:
-[1] LC_COLLATE=Chinese (Traditional)_Taiwan.950
-[2] LC_CTYPE=Chinese (Traditional)_Taiwan.950
-[3] LC_MONETARY=Chinese (Traditional)_Taiwan.950
-[4] LC_NUMERIC=C
-[5] LC_TIME=Chinese (Traditional)_Taiwan.950
-system code page: 65001
-
-attached base packages:
-[1] stats graphics grDevices utils datasets methods
-[7] base
-
-other attached packages:
- [1] leaflet_2.1.2 marmap_1.0.10
- [3] GISTools_0.7-4 rgeos_0.6-3
- [5] MASS_7.3-58.1 RColorBrewer_1.1-3
- [7] mapr_0.5.2 rgbif_3.7.7
- [9] ggspatial_1.1.8 rnaturalearthdata_0.1.0
-[11] rnaturalearth_0.3.3 sf_1.0-13
-[13] ggplot2_3.4.2 raster_3.6-20
-[15] rgdal_1.6-6 maptools_1.1-7
-[17] sp_1.6-0
-
-loaded via a namespace (and not attached):
- [1] oai_0.4.0 bit64_4.0.5 httr_1.4.6
- [4] tools_4.2.2 bslib_0.4.2 utf8_1.2.3
- [7] R6_2.5.1 KernSmooth_2.23-20 DBI_1.1.3
-[10] lazyeval_0.2.2 colorspace_2.1-0 withr_2.5.0
-[13] tidyselect_1.2.0 downlit_0.4.2 bit_4.0.5
-[16] compiler_4.2.2 cli_3.6.1 xml2_1.3.4
-[19] sass_0.4.6 scales_1.2.1 classInt_0.4-9
-[22] proxy_0.4-27 stringr_1.5.0 digest_0.6.31
-[25] foreign_0.8-83 rmarkdown_2.21 pkgconfig_2.0.3
-[28] htmltools_0.5.5 fastmap_1.1.1 highr_0.10
-[31] maps_3.4.1 htmlwidgets_1.6.2 rlang_1.1.1
-[34] rstudioapi_0.14 RSQLite_2.3.1 shape_1.4.6
-[37] jquerylib_0.1.4 generics_0.1.3 adehabitatMA_0.3.16
-[40] jsonlite_1.8.4 crosstalk_1.2.0 dplyr_1.1.2
-[43] distill_1.5 magrittr_2.0.3 gdistance_1.6.4
-[46] Matrix_1.6-0 Rcpp_1.0.10 munsell_0.5.0
-[49] fansi_1.0.4 lifecycle_1.0.3 terra_1.7-29
-[52] stringi_1.7.12 whisker_0.4.1 yaml_2.3.7
-[55] plyr_1.8.8 grid_4.2.2 blob_1.2.4
-[58] lattice_0.20-45 knitr_1.42 pillar_1.9.0
-[61] igraph_1.4.3 reshape2_1.4.4 codetools_0.2-18
-[64] glue_1.6.2 evaluate_0.21 data.table_1.14.8
-[67] vctrs_0.6.2 gtable_0.3.3 cachem_1.0.8
-[70] xfun_0.39 mime_0.12 e1071_1.7-13
-[73] class_7.3-20 ncdf4_1.21 tibble_3.2.1
-[76] memoise_2.0.1 units_0.8-2 You can go then to the website and download data locally. Alternatively, using the code source of the page your can target directly the url of the data set you want to download (I failed many time before identifyinng the right link). This is a zip file that you want to temporary store using tempfile and immediately extract it in a designated directory. Again here, nothing about mapping but you have a look on the details of the code below.
#bug file if file existes
-url <- 'https://data.moi.gov.tw/MoiOD/System/DownloadFile.aspx?DATA=72874C55-884D-4CEA-B7D6-F60B0BE85AB0'
-path1 <- tempfile(fileext = ".zip")
-if (file.exists(path1)) 'file alredy exists' else download.file(url, path1, mode="wb")
-zip::unzip(zipfile = path1,exdir = 'Data')The function readOGR reads the .shp file. Other arguments are for the file encoding and to deal with the Chinese characters because my system is in English.
taiwan <- readOGR('Data/COUNTY_MOI_1090820.shp', use_iconv=TRUE, encoding='UTF-8')OGR data source with driver: ESRI Shapefile
-Source: "D:\Dropbox\NTU\IONTU\Class_lectures\OCEAN_5098_R\OCEAN5098\df4810487c8a45838b57b12f305a2529def7225f\Data\COUNTY_MOI_1090820.shp", layer: "COUNTY_MOI_1090820"
-with 22 features
-It has 4 fieldsMake the plot:
- -
This is a map of TAIWAN using official data from the government. It includes all remote territories such as Taiping and Diaoyutai islands. County names are here accessible using:
-taiwan$COUNTYNAME [1] "連江縣" "宜蘭縣" "彰化縣" "南投縣" "雲林縣" "基隆市" "臺北市"
- [8] "新北市" "臺中市" "臺南市" "桃園市" "苗栗縣" "嘉義市" "嘉義縣"
-[15] "金門縣" "高雄市" "臺東縣" "花蓮縣" "澎湖縣" "新竹市" "新竹縣"
-[22] "屏東縣"taiwan$COUNTYENG [1] "Lienchiang County" "Yilan County" "Changhua County"
- [4] "Nantou County" "Yunlin County" "Keelung City"
- [7] "Taipei City" "New Taipei City" "Taichung City"
-[10] "Tainan City" "Taoyuan City" "Miaoli County"
-[13] "Chiayi City" "Chiayi County" "Kinmen County"
-[16] "Kaohsiung City" "Taitung County" "Hualien County"
-[19] "Penghu County" "Hsinchu City" "Hsinchu County"
-[22] "Pingtung County" ggplot2 & sf
-As we mentioned earlier, the package ggplot2 implements the grammar of graphics in R. While ggplot2 is becoming the standard for R graphs, it does not handle spatial data specifically. The current state-of-the-art of spatial objects in R relies on Spatial classes defined in the package sp, but the new package sf has implemented the “simple feature” standard, and is steadily taking over sp. Relatively recently, the package ggplot2 has allowed the use of simple features from the package sf as layers in a graph (since the version 3.0.0 of ggplot2). The combination of ggplot2 and sf therefore enables to programmatically create maps, using the grammar of graphics, just as informative or visually appealing as traditional GIS software.
Note that the conversion from sp to sf is achievable using the function st_as_sf from the package sf.
# not run
-# wrld_simpl <- st_as_sf(wrld_simpl)Theme and datasets
-After loading the basic packages necessary for all maps, i.e. ggplot2 and sf. We can change the theme of our plot. The classic dark-on-light theme for ggplot2 (theme_bw) is appropriate for maps (default: theme_grey()):
The package rnaturalearth provides a map of countries of the entire world. Use ne_countries to pull country data and choose the scale (rnaturalearthhires is necessary for scale = "large"). The function can return sp classes (default) or directly sf classes, as defined in the argument returnclass:
world <- ne_countries(scale = "medium", returnclass = "sf")
-class(world)[1] "sf" "data.frame"In addition to country shape polygones, this dataset contains information on the population in every country.
-A ggplot map
-First, let us start with creating a base map of the world using ggplot2. This base map will then be extended with different map elements, as well as zoomed in to an area of interest. We can check that the world map was properly retrieved and converted into an sf object, and plot it with ggplot2. This time we will indicate that the geometry of our plot is constrained by a sf object:
In your map of the world, the plot panel is expanded beyond the size of the earth (you can see that the graticule lines end before the edge of the plot panel), and hence no axis ticks are drawn. One way to solve the issue is to turn off the expansion while defining the coordinates.
- -
As before, the layers are added one at a time in a ggplot call, so the order of each layer is very important. All data will have to be in an sf format to be used by ggplot2; data in other formats (e.g. classes from sp) will be manually converted to sf classes if necessary.
-
-
- Title, subtitle, and axis labels using
ggtitle,xlab,ylab
-
A title and a subtitle can be added to the map using the function ggtitle, passing any valid character string (e.g. with quotation marks) as arguments. Axis names are absent by default on a map, but can be changed to something more suitable (e.g. “Longitude” and “Latitude”), depending on the map:
ggplot(data = world) +
- geom_sf() +
- coord_sf(expand = FALSE) +
- xlab("Longitude") + ylab("Latitude") +
- ggtitle("World map", subtitle = paste0("(", length(unique(world$name)), " countries)"))
-
-
- Map color in
geom_sf
-
In many ways, sf geometries are no different than regular geometries, and can be displayed with the same level of control on their attributes. Here is an example with the polygons of the countries filled with a green color (argument fill), using black for the outline of the countries (argument color):
ggplot(data = world) +
- geom_sf(color = "black", fill = "lightgreen") +
- coord_sf(expand = FALSE) 
The package ggplot2 allows the use of more complex color schemes, such as a gradient on one variable of the data. Here is another example that shows the population of each country. In this example, we use the “viridis” colorblind-friendly palette for the color gradient (with option = "magma" for the magma variant), using the square root of the population (which is stored in the variable POP_EST of the world object):
ggplot(data = world) +
- geom_sf(aes(fill = pop_est)) +
- coord_sf(expand = FALSE) +
- scale_fill_viridis_c(option = "plasma", trans = "sqrt") # sqrt transform
-
-
- Projection and extend (
coord_sf)
-
The function coord_sf allows to deal with the coordinate system, which includes both projection and extent of the map. By default, the map will use the coordinate system of the first layer that defines one (i.e. scanned in the order provided), or if none, fall back on WGS84 (latitude/longitude, the reference system used in GPS). Using the argument crs, it is possible to override this setting, and project on the fly to any projection. This can be achieved using any valid PROJ4 string (here, the European-centric ETRS89 Lambert Azimuthal Equal-Area projection):
ggplot(data = world) +
- geom_sf() +
- scale_fill_viridis_c(option = "plasma", trans = "sqrt") +
- coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")
Spatial Reference System Identifier (SRID) or an European Petroleum Survey Group (EPSG) code are available for the projection of interest, they can be used directly instead of the full PROJ4 string. The two following calls are equivalent for the ETRS89 Lambert Azimuthal Equal-Area projection, which is EPSG code 3035:
-The extent of the map can also be set in coord_sf, in practice allowing to “zoom” in the area of interest, provided by limits on the x-axis (xlim), and on the y-axis (ylim). Note that the limits are automatically expanded by a fraction to ensure that data and axes don’t overlap; it can also be turned off to exactly match the limits provided with expand = FALSE:
ggplot(data = world) +
- geom_sf(aes(fill = pop_est)) +
- coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE) +
- scale_fill_viridis_c(option = "plasma") # linear scale
ggspatial elements
--
-
- Scale bar and North arrow -
Several packages are available to create a scale bar on a map (e.g. prettymapr, vcd, ggsn, or legendMap). Here the package ggspatial provides easy-to-use functions. scale_bar that allows to add simultaneously the north symbol and a scale bar into the ggplot map. Five arguments need to be set manually: lon, lat, distance_lon, distance_lat, and distance_legend. The location of the scale bar has to be specified in longitude/latitude in the lon and lat arguments. The shaded distance inside the scale bar is controlled by the distance_lon argument while its width is determined by distance_lat. Additionally, it is possible to change the font size for the legend of the scale bar (argument legend_size, which defaults to 3). The North arrow behind the “N” north symbol can also be adjusted for its length (arrow_length), its distance to the scale (arrow_distance), or the size the N north symbol itself (arrow_north_size, which defaults to 6). Note that all distances (distance_lon, distance_lat, distance_legend, arrow_length, arrow_distance) are set to "km" by default in distance_unit; they can also be set to nautical miles with “nm”, or miles with “mi”.
ggplot(data = world) +
- geom_sf(aes(fill = pop_est)) +
- coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE) +
- scale_fill_viridis_c(option = "plasma") +
- annotation_scale(location = "br", width_hint = 0.3) +
- annotation_north_arrow(location = "br", which_north = "true",
- pad_x = unit(0.5, "cm"), pad_y = unit(1, "cm"),
- style = north_arrow_fancy_orienteering) 
Note: if you plot a larger area, you may get a warning on the inaccuracy of the scale bar.
--
-
- Country names and other names (
geom_textandannotate)
-
The world data set already contains country names and the coordinates of the centroid of each country (among more information). We can use this information to plot country names, using world as a regular data.frame in ggplot2. The function geom_text can be used to add a layer of text to a map using geographic coordinates. The function requires the data needed to enter the country names, which is the same data as the world map. Again, we have a very flexible control to adjust the text at will on many aspects:
-
-
The size (argument
size);
-The alignment, which is centered by default on the coordinates provided. The text can be adjusted horizontally or vertically using the arguments
hjustandvjust, which can either be a number between 0 (right/bottom) and 1 (top/left) or a character (“left”, “middle”, “right”, “bottom”, “center”, “top”). The text can also be offset horizontally or vertically with the argumentnudge_xandnudge_y;
-The font of the text, for instance its color (argument
color) or the type of font (fontface);
-The overlap of labels, using the argument
check_overlap, which removes overlapping text. Alternatively, when there is a lot of overlapping labels, the packageggrepelprovides ageom_text_repelfunction that moves label around so that they do not overlap.
-For the text labels, we are defining the centroid of the countries with
st_centroid, from the packagesf. Then we combined the coordinates with the centroid, in thegeometryof the spatial data frame. The packagesfis necessary for the commandst_centroid.
-
Additionally, the annotate function can be used to add a single character string at a specific location, as demonstrated here to add the “Gulf of Mexico”Pacific Ocean” and Ryukyu Archipelago
sf::sf_use_s2(FALSE) # FOR ERROR turn off the s2 processing
-world_points<- st_centroid(world)
-world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))
-
-ggplot(data = world) +
-geom_sf(aes(fill = pop_est)) +
-geom_text(data= world_points,aes(x=X, y=Y, label=name),
- color = "grey", fontface = "bold", check_overlap = FALSE) +
- scale_fill_viridis_c(option = "plasma") +
- annotate(geom = "text", x = 124, y = 21, label = "Pacific Ocean", fontface = "italic", color = "#0b3c8a", size = 5) +
- annotate(geom = "text", x = 124.2, y = 24, label = "Ryukyu archipelago", fontface = "italic", color = "#d41919", size = 3) +
- coord_sf(xlim = c(118, 128), ylim = c(17, 27), expand = FALSE)
To goo further: Finalization
-Now to make the final touches, the theme of the map can be edited to make it more appealing. We suggested the use of theme_bw for a standard theme, but there are many other themes that can be selected from (see for instance ?ggtheme in ggplot2, or the package ggthemes which provide several useful themes). Moreover, specific theme elements can be tweaked to get to the final outcome:
-
-
Position of the legend: Although not used in this example, the argument
legend.positionallows to automatically place the legend at a specific location (e.g. “topright”, “bottomleft”, etc.);
-Grid lines (graticules) on the map: by using
panel.grid.majorandpanel.grid.minor, grid lines can be adjusted. Here we set them to a gray color and dashed line type to clearly distinguish them from country borders lines;
-Map background: the argument
panel.backgroundcan be used to color the background, which is the ocean essentially, with a light blue; -Many more elements of a theme can be adjusted, which would be too long to cover here. We refer the reader to the documentation for the functiontheme.
-
ggplot(data = world) +
- geom_sf(fill= 'antiquewhite') +
- geom_text(data= world_points,aes(x=X, y=Y, label=name), color = 'darkblue', fontface = 'bold', check_overlap = FALSE) +
- annotate(geom = 'text', x = -90, y = 26, label = 'Gulf of Mexico', fontface = 'italic', color = 'grey22', size = 6) +
- annotation_scale(location = 'bl', width_hint = 0.5) +
- annotation_north_arrow(location = 'bl', which_north = 'true', pad_x = unit(0.75, 'in'), pad_y = unit(0.5, 'in'), style = north_arrow_fancy_orienteering) +
- coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97), expand = FALSE) +
- xlab('Longitude') + ylab('Latitude') +
- ggtitle('Map of the Gulf of Mexico and the Caribbean Sea') +
- theme(panel.grid.major = element_line(color = gray(.5), linetype = 'dashed', size = 0.5), panel.background = element_rect(fill = 'aliceblue'))
Exporting map
-The final map is now ready, it is very easy to save it using ggsave. This function allows a graphic (typically the last plot displayed) to be saved in a variety of formats, including the most common PNG (raster bitmap) and PDF (vector graphics), with control over the size and resolution of the outcome. For instance here, we save a PDF version of the map, which keeps the best quality, and a PNG version of it for web purposes:
-Special
-Species distribution rgbif
-Import Global Biodiversity Information Facility (GBIF) data using rgbif package and crete distribution map using the mapr package.
gbif.res <- occ_search(scientificName = "Urocissa caerulea", limit = 1200)
-
-map_ggplot(gbif.res) +
- coord_sf(xlim = c(120, 123), ylim = c(21, 26), expand = FALSE)
Note the slight difference in overlap of the distribution data on this shape file. To date, I could not figure out the problem here.
-Have a look here, to combine climate and species data.
-Bathymetric map with marmap
-marmap can query and plot NOAA’s bathymetry database
# query
-TW.bathy <- getNOAA.bathy(lon1=118,lon2=124, lat1=21,lat2=26,resolution=1) # don't put too wide / resolution: 1
-# define palette
-blues <- colorRampPalette(c("darkblue", "cyan"))
-greys <- colorRampPalette(c(grey(0.4),grey(0.99)))
-# make the plot
-plot.bathy(TW.bathy,
- image=T,
- deepest.isobath=c(-6000,-120,-30,0),
- shallowest.isobath=c(-1000,-60,0,0),
- step=c(2000,60,30,0),
- lwd=c(0.3,1,1,2),
- lty=c(1,1,1,1),
- col=c("grey","black","black","black"),
- drawlabels=c(T,T,T,F),
- bpal = list(c(0,max(TW.bathy),greys(100)),c(min(TW.bathy),0,blues(100))),
- land=T, xaxs="i"
- )
Profiles can be extract using get.transect:
tw.profile <-get.transect(TW.bathy,x1=119.5,y1=23.75, x2=122,y2=23.75, dis=TRUE)
-plotProfile(tw.profile) 
-#### Not Run: extract a profile Manually
-#### manual.profile<-get.transect (TW.bathy, loc=T,dist=T)
-#### plotProfile(manual.profile)Interactive maps
-Leaflet is one of the most popular open-source JavaScript libraries for interactive maps. You can loose hours and its functionalities are amazing. The R package leaflet makes it easy to integrate and control Leaflet maps in R.
FRE <- paste(sep = "<br/>",
- "<b><a href='https://www.dipintothereef.com/'>FRELAb TAIWAN</a></b>",
- "Functional Reef Ecology Lab",
- "Institute of Oceanography, NTU")
-
-
-leaflet(taiwan) %>%
- addPolygons(weight=0.5) %>%
- addTiles(group="Kort") %>%
- addPopups(121.53725, 25.021252, FRE, options = popupOptions(closeButton = FALSE))Combo maps
-Interactive X GBIF ?
-See: https://poldham.github.io/abs/mapgbif.html & https://data-blog.gbif.org/post/gbif-maps-api-using-r-and-leaflet/ among many many other onlin resources
-
-⚠ Practice 5.1 Create and customize an interactive map of your choice after exploring the functionality of the leaflet package. Make it obvious it is your map and not something you copy online by highlighting points of interest. Push both your .Rmd and .html files into a public repository available from your Github account.You will share with me be email [vianneydenis@g.ntu.edu.tw] the address (URL) of this repository (such as https://github.com/vianneydenis/OCEAN5098.git) before next Monday in order for me to check your work. The title of your email should be `Practice 5.1 (your name: your student no.). BE CREATIVE ;)
-