maps1.Rmd

---
title: "Maps (1)"
output: 
  html_document:
    fig_caption: no
    number_sections: yes
    toc: yes
    toc_float: false
    collapsed: no
---

```{r maps1-1, echo=FALSE}
options(width = 105)
knitr::opts_chunk$set(dev='png', dpi=300, cache=TRUE, out.width = "80%", out.height = "80%", verbose=TRUE)
pdf.options(useDingbats = TRUE)
klippy::klippy(position = c('top', 'right'))
```
<p><span style="color: #00cc00;">NOTE:  This page has been revised for Winter 2024, but may undergo further edits.</span></p>

# Introduction #

R has the ability through the `{maps}` package and the base graphics to generate maps, but such "out-of-the-box" maps, like other base graphics-generated illustrations, these may not be suitable for immediate publication.  Other options exit, for example, using packages like `{ggmaps}`, `{mapview},` or `{tmap}` However, the ability of the `{sf}` package to handle the diverse kinds of spatial data, including shapefiles, provides a facility for creating publishable maps.  The examples here use a set of shapefiles downloaded from the *Natural Earth* web page [[http://www.naturalearthdata.com]](http://www.naturalearthdata.com).  

There are two R packages that can be used to directly download Natural Earth data within a script: `{rnaturalearth}`, and `{rnaturalearthdata}`, but it's usually the case that one visits the Natural Earth web page anyway, and it's easy to download (and save for reuse) files from the web page. Shapefiles and rasters are available at three scales 1:10m (small spatial scale, large cartographic scale), 1:50m, and 1:110m (large spatial scale, small cartographic scale).

This first example recreates Figure 1 from  Marlon, J.R., et al., 2016,  Reconstructions of biomass burning from sediment-charcoal records to improve data–model comparisons. *Biogeosciences* 13:3225-3244. [[http://www.biogeosciences.net/13/3225/2016/]](http://www.biogeosciences.net/13/3225/2016/).  There are two general steps in the production of this map, the first involves 1) projecting the *Natural Earth* shapefiles into the Robinson projection (saving the shapefiles, which can be reused for other applications), and the second the projection of the data and assembly into the finished map.

# Set up outlines and polygons for global (Robinson projection) maps #

## Read Natural Earth shapefiles ##

The first step is to read a set of shapefiles downloaded from the Natural Earth web page [[http://www.naturalearthdata.com]](http://www.naturalearthdata.com), and project these into the Robinson projection.  Two CRAN packages, `{rnaturalearth}` and `{rnaturalearthdata}` can be used to manage and download the data within R, but the *Natural Earth* web pages are worth looking at.

Begin by loading the appropriate packages:

```{r maps1-2, message=FALSE}
library(sf)
library(stars)
library(ggplot2)
library(RColorBrewer)
library(classInt)
```

Now read the shape files, including those for global coastlines, adminstrative units (borders), large lakes, and a graticule.  Set the filenames:

```{r maps1-3 }
# Natural Earth shape files -- global (Robinson) projections
# get shapefiles from http://www.naturalearthdata.com
shape_path <- "/Users/bartlein/Projects/RESS/data/RMaps/source/"
coast_shapefile <- paste(shape_path, "ne_50m_coastline/ne_50m_coastline.shp", sep="")
ocean_shapefile <- paste(shape_path, "ne_50m_ocean/ne_50m_ocean.shp", sep="")
land_shapefile <- paste(shape_path, "ne_50m_land/ne_50m_land.shp", sep="")
admin0_shapefile <- paste(shape_path, "ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp", sep="")
admin1_shapefile <- paste(shape_path, 
      "ne_50m_admin_1_states_provinces_lakes/ne_50m_admin_1_states_provinces_lakes.shp", sep="")
lakes_shapefile <- paste(shape_path, "ne_50m_lakes/ne_50m_lakes.shp", sep="")
bb_shapefile <- paste(shape_path, "ne_50m_graticules_all/ne_50m_wgs84_bounding_box.shp", sep="")
grat30_shapefile <- paste(shape_path, "ne_50m_graticules_all/ne_50m_graticules_30.shp", sep="")
```

Before reading each file, determine what kind it is (e.g. lines vs. polygons). This information will be used to (manually, but it could be coded) to name the `sf` objects. read and plot it. As an example:

```{r maps1-4 }
# query geometry type
coast_geometry <- as.character(st_geometry_type(st_read(coast_shapefile), by_geometry = FALSE))
coast_geometry
```

So the `coast_shapefile` has a "MULTILINESTRING" geometry type, so we'll call the input `sf` object "`coast_lines_sf`"

In the 'plot()' functions below, use the `col` argument for lines, `border` for polygon border colors. Read and plot the shapefiles:

```{r maps1-5}
# read coastline
coast_lines_sf <- st_read(coast_shapefile) # note geometry type MULTILINESTRING
plot(st_geometry(coast_lines_sf), col="gray50")
```

Read and plot the other shapefiles. (Note that the summary output is suppressed.)

```{r maps1-6, results="hide", message=FALSE, cache=TRUE}
# read and plot other shape files, noting geometry types
ocean_poly_sf <- st_read(ocean_shapefile) # note:  geometry type MULTIPOLYGON
plot(st_geometry(ocean_poly_sf), col="gray80")
land_poly_sf <- st_read(land_shapefile) # note:  geometry type MULTIPOLYGON
plot(st_geometry(land_poly_sf), col="gray80")
admin0_poly_sf <- st_read(admin0_shapefile) # note:  geometry type MULTIPOLYGON
plot(st_geometry(admin0_poly_sf), col="gray70", border="red")
admin1_poly_sf <- st_read(admin1_shapefile) # note:  geometry type MULTIPOLYGON
plot(st_geometry(admin1_poly_sf), col="gray70", border="pink")
lakes_poly_sf <- st_read(lakes_shapefile) # note:  geometry type POLYGON
plot(st_geometry(lakes_poly_sf), col="blue")
bb_poly_sf <- st_read(bb_shapefile) # note:  geometry type POLYGON
plot(st_geometry(bb_poly_sf), col="gray70")
grat30_lines_sf <- st_read(grat30_shapefile) # note:  geometry type LINESTRING
plot(st_geometry(grat30_lines_sf), col="black")
```

Filter the `lakes_poly_sf` object to the the largest lakes only:

```{r maps1-7 }
# get large lakes only
lrglakes_poly_sf <- lakes_poly_sf[as.numeric(lakes_poly_sf$scalerank) <= 2 ,]
plot(st_geometry(lrglakes_poly_sf), col="lightblue")
```


Plot everything.


```{r maps1-8, cache=TRUE}
# plot everything
plot(st_geometry(coast_lines_sf), col="black")
plot(st_geometry(ocean_poly_sf), col="skyblue1", add=TRUE)
plot(st_geometry(land_poly_sf), col="palegreen", add=TRUE)
plot(st_geometry(admin0_poly_sf), border="red", add=TRUE)
plot(st_geometry(admin1_poly_sf), border="pink", add=TRUE)
plot(st_geometry(lakes_poly_sf), border="lightblue", add=TRUE)
plot(st_geometry(lrglakes_poly_sf), border="blue", add=TRUE)
plot(st_geometry(grat30_lines_sf), col="gray50", add=TRUE)
plot(st_geometry(bb_poly_sf), border="black", add=TRUE)
plot(st_geometry(coast_lines_sf), col="purple", add=TRUE)
```

## Project the shapefiles ##

Next, transform (project) the individual shapefiles into the Robinson projection using `st_transform()`.  Get the current coordinate reference system (i.e. the `PROJ4` string).

```{r maps1-9 }
# get CRS and projstring
unproj_projstring <- st_crs(coast_lines_sf)
unproj_projstring
```
Set the `PROJ` string for the projected (Robinson) `sf` objects:

```{r maps1-10}
# set new projstring
robinson_projstring <- "+proj=robin +lon_0=0w"
robinson_projstring
```

Now do the projections.

```{r maps1-11 }
# do projections
coast_lines_proj <- st_transform(coast_lines_sf, crs = st_crs(robinson_projstring))
ocean_poly_proj <- st_transform(ocean_poly_sf, crs = st_crs(robinson_projstring))
land_poly_proj <- st_transform(land_poly_sf, crs = st_crs(robinson_projstring))
admin0_poly_proj <- st_transform(admin0_poly_sf, crs = st_crs(robinson_projstring))
admin1_poly_proj <- st_transform(admin1_poly_sf, crs = st_crs(robinson_projstring))
lakes_poly_proj <- st_transform(lakes_poly_sf, crs = st_crs(robinson_projstring))
lrglakes_poly_proj <- st_transform(lrglakes_poly_sf, crs = st_crs(robinson_projstring))
grat30_lines_proj <- st_transform(grat30_lines_sf, crs = st_crs(robinson_projstring))
bb_poly_proj <- st_transform(bb_poly_sf, crs = st_crs(robinson_projstring))
```

Plot the projected shapefiles.

```{r maps1-12 }
# plot projected files
plot(st_geometry(bb_poly_proj), col="gray90", border="black", lwd=2)
plot(st_geometry(coast_lines_proj), col="black", add=TRUE)
plot(st_geometry(ocean_poly_proj), col="skyblue1", add=TRUE)
plot(st_geometry(land_poly_proj), col="palegreen", add=TRUE)
plot(st_geometry(admin0_poly_proj), border="red", add=TRUE)
plot(st_geometry(admin1_poly_proj), border="pink", add=TRUE)
plot(st_geometry(lakes_poly_proj), col="lightblue", add=TRUE)
plot(st_geometry(lrglakes_poly_proj), col="blue", add=TRUE)
plot(st_geometry(grat30_lines_proj), col="gray50", add=TRUE)
plot(st_geometry(coast_lines_proj), col="purple", add=TRUE)
plot(st_geometry(bb_poly_proj), border="black", lwd=2, add=TRUE)
```

Because the map we're creating will use the borders and lakes only as outlines (as opposed to polygons that might be filled with a specific color), those polygons can be converted to `MULTILINESTRING` objects.

```{r maps1-13 }
# convert MULTIPOLYGON to MULTILINESTRING
admin0_lines_proj <- st_cast(admin0_poly_proj, "MULTILINESTRING")
admin1_lines_proj <- st_cast(admin1_poly_proj, "MULTILINESTRING")
lakes_lines_proj <- st_cast(lakes_poly_proj, "MULTILINESTRING")
lrglakes_lines_proj <- st_cast(lrglakes_poly_proj, "MULTILINESTRING")
bb_lines_proj <- st_cast(bb_poly_proj, "MULTILINESTRING")
```

Plot the new `MULTILINESTRING` objects

```{r maps1-14 }
# test the MULTILINESTRING objects
plot(st_geometry(bb_lines_proj), col="black")
plot(st_geometry(coast_lines_proj), col="green", add=TRUE)
plot(st_geometry(admin0_lines_proj), col="lightblue", add=TRUE)
plot(st_geometry(admin1_lines_proj), col="lightblue", add=TRUE)
plot(st_geometry(lakes_lines_proj), col="blue", add=TRUE)
plot(st_geometry(lrglakes_lines_proj), col="purple", add=TRUE)
plot(st_geometry(grat30_lines_proj), col="gray", add=TRUE)
plot(st_geometry(coast_lines_proj), col="black", add=TRUE)
plot(st_geometry(bb_lines_proj), border="black", lwd = 2, add=TRUE)
```

## Write out the projected shapefiles ##

Next, write out the projected shapefiles, first setting the output path.

```{r maps1-15}
# write out the various projected shapefiles 
outpath <- "/Users/bartlein/Projects/RESS/data/RMaps/derived/glRob_50m/"
```

Write out the projected coastlines.

```{r maps1-16}
# coast lines 
outshape <- coast_lines_proj
outfile <- "glRob_50m_coast_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(coast_lines_proj, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```
It's always good practice to test whether the shapefile has indeed been written out correctly.  Read it back in and plot it.

```{r maps1-17 }
# test
test_sf <- st_read(outshapefile)
test_outline <- st_geometry(test_sf)
plot(test_outline, col="gray")
```

Write out the other shapefiles.

```{r maps1-18, eval=TRUE, results = "hide", message=FALSE, warning=FALSE}
# write out the other objects as shapefiles
outshape <- bb_poly_proj
outfile <- "glRob_50m_bb_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- bb_lines_proj
outfile <- "glRob_50m_bb_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- ocean_poly_proj
outfile <- "glRob_50m_ocean_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- land_poly_proj
outfile <- "glRob_50m_land_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- admin0_poly_proj
outfile <- "glRob_50m_admin0_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- admin0_lines_proj
outfile <- "glRob_50m_admin0_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- admin1_poly_proj
outfile <- "glRob_50m_admin1_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- admin1_lines_proj
outfile <- "glRob_50m_admin1_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- lakes_poly_proj
outfile <- "glRob_50m_lakes_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- lakes_lines_proj
outfile <- "glRob_50m_lakes_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- lrglakes_poly_proj
outfile <- "glRob_50m_lrglakes_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- lrglakes_lines_proj
outfile <- "glRob_50m_lrglakes_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

outshape <- grat30_lines_proj
outfile <- "glRob_50m_grat30_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```

## Clip out a polygon for the Caspian ## 

Another setup step, and again one that creates some reusable files, is to make a polygon shape file for the Caspian Sea, which can be plotted and filled with color. Some terrestrial databases, including the treecover data that will be plotted below include non-missing data for gridpoints that cover the Caspian Sea (or sometimes the Black Sea). It's useful to be able to knock out those points by plotting a filled polygon with the same color that is used for the ocean. The steps involved include trimming to coastlines file to the area around the Caspian, and then turning those lines into a polygon (which can be projected)>  

The first step is to creat a "bounding box" that surrouds the Caspian.

```{r maps1-19 }
# Caspian
caspian_bb <- st_bbox(c(xmin = 45, ymin = 35, xmax = 56, ymax = 50), crs = unproj_projstring)
caspian_bb <- st_as_sfc(caspian_bb)
```

```{r maps1-20 }
# get the points that define the ouline of the Caspian
# get the points that define the outline of the Caspian
caspian_lines <- st_intersection(coast_lines_sf, caspian_bb )
caspian_lines
```

"Recast" into a polygon:

```{r maps1-21}
caspian_poly <- st_cast(caspian_lines, "POLYGON")
caspian_poly
```

Plot the Caspian polygon

```{r maps1-22}
plot(st_geometry(caspian_poly), col = "skyblue")
```

Project and plot the Caspian polygon:

```{r maps1-23 }
# project the Caspian outline
caspian_poly_proj <- st_transform(caspian_poly, crs = st_crs(robinson_projstring))
plot(st_geometry(caspian_poly_proj), col = "skyblue")
```

Write out the Caspian polygon.

```{r maps1-24 }
# save the Caspian outlines
outshape <- caspian_poly
outfile <- "glRob_50m_caspian_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(caspian_poly, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# projected lines
outshape <- caspian_poly_proj
outfile <- "glRob_50m_caspian_poly_proj/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(caspian_poly_proj, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```

[[Back to top]](Rmaps.html)

# Set up tree-cover data #

One version of Fig. 1 in Marlon et al. (2016) includes a gray-shade "layer" of the "UMD" tree-cover data (Defries et al., 2000, A new global 1-km dataset of percentage tree cover derived from remote sensing, *Global Change Biology* 6:247-254)  [[http://glcf.umd.edu/data/treecover/data.shtml]](http://glcf.umd.edu/data/treecover/data.shtml), used to provide a background context for the location of the GCDv3 sites.  The data were converted to a netCDF file, and this is read, and in turn converted to an `sf` POINTS and POLYGON

## Read the tree-cover data ##

```{r maps1-25, eval=TRUE, echo=TRUE}
# read a single-variable netCDF dataset using stars to read_ncdf
tree_path <- "/Users/bartlein/Projects/RESS/data/nc_files/"
tree_name <- "treecov.nc"
tree_file <- paste(tree_path, tree_name, sep="")
tree <- read_ncdf(tree_file, var = "treecov")
tree
```

Convert to points:

```{r maps1-26}
# convert stars object to sf POINTS
treecov_pts <- st_as_sf(tree, as_points = TRUE)
plot(treecov_pts, pch = 16, cex = 0.5)

# also plot a zoomed-in view of the western U.S.
plot(treecov_pts, xlim = c(-125, -100), ylim = c(30, 50), pch = 16, main = "treecov_pts")
```

Convert to polygons, and plot a zoomed-in look at the western U.S:
```{r maps1-27, cache=TRUE}
# convert stars object to sf POLYGONS
treecov_poly <- st_as_sf(tree, as_points = FALSE)
treecov_poly

# check western U.S.
plot(treecov_poly, xlim = c(-125, -100), ylim = c(30, 50), main = "treecov_poly")
```

Write out the shapefiles.

```{r maps1-28 }
# write out the unprojected points
outpath <- "/Users/bartlein/Projects/RESS/data/RMaps/derived/treecov/"
outfile <- "treecov_pts/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(treecov_pts, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```

```{r maps1-29}
# write out the unprojected polygons
outshape <- treecov_poly
outfile <- "treecov_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(treecov_poly, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```

## Project the treecover points and polygons ##
 
Project the tree-cover data:

```{r maps1-30 }
# project the treecover data
treecov_pts_proj <- st_transform(treecov_pts, crs = st_crs(robinson_projstring))
treecov_poly_proj <- st_transform(treecov_poly, crs = st_crs(robinson_projstring))
```

Write out the projected points and polygons:

```{r maps1-31}
# write out the projected points
outshape <- treecov_pts_proj
outfile <- "treecov_pts_proj/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(treecov_pts_proj, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```
```{r maps1-32}
# write out the projected polygons
outshape <- treecov_poly_proj
outfile <- "treecov_poly_proj/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(treecov_poly_proj, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```

[[Back to top]](rmaps1.html)

# Map the GCDv3 charcoal records #

Figure 1 of Marlon et al. (2016) shows the distribution of charcoal records, and additionally shows by means of symbol color the number of samples in each record.  Two versions of the figure will be produced:  1) as published, and 2) with the tree-cover data as a gray-shaded background.

## Read the projected shape files ##

We'll assume that the necessary projected map and tree cover files are available in the current Environment. Otherwise, they could be read in as usual here.

Set cutpoints and gray-scale color numbers for plotting the tree-cover data.  These were calculated to follow a "Munsell sequence" of equally perceived class intervals.

## Read the GCDv3 data ##

Read the GCDv3 charcoal site-location data.  

```{r maps1-33 }
# read the data
csvpath <- "/Users/bartlein/Projects/RESS/data/csv_files/"
csvname <- "GCDv3_MapData_Fig1.csv"
gcdv3_csv <- read.csv(paste(csvpath, csvname, sep=""))
gcdv3_csv <- data.frame(cbind(gcdv3_csv$Long, gcdv3_csv$Lat, gcdv3_csv$samples22k))
names(gcdv3_csv) <- c("lon", "lat", "samples22k")
head(gcdv3_csv)
```

Convert to POINTS:

```{r maps1-34}
# convert to sf POINTS
gcdv3_pts <- st_as_sf(gcdv3_csv, coords = c("lon", "lat"))
gcdv3_pts <- st_set_crs(gcdv3_pts, unproj_projstring)
gcdv3_pts
```

```{r maps1-35}
# project the gcdv3 data
gcdv3_pts_proj <- st_transform(gcdv3_pts, crs = st_crs(robinson_projstring))
```

Set colors for tree cover data:

```{r maps1-36}
# colors for tree cover
treecov_clr_upper <- c(   1,    15,    20,    25,    30,    35,    40,    45,    50,    55,    60,    65,    70,    75,   999)
treecov_clr_gray <- c(1.000, 0.979, 0.954, 0.926, 0.894, 0.859, 0.820, 0.778, 0.733, 0.684, 0.632, 0.576, 0.517, 0.455, 0.389)
colnum <- findInterval(treecov_poly_proj$treecov, treecov_clr_upper)+1
clr <- gray(treecov_clr_gray[colnum], alpha=NULL)
```

Set colors and symbol sizes for characoal data:

```{r maps1-37}
# set symbol size and colors gcdv3
nsamp_cutpts <- c(10,100,1000,10000)
nsamp_colors <- c("green3","deepskyblue2","slateblue4","purple")
nsamp_cex <- c(0.5,0.5,0.5,0.5)
nsamp_num <- findInterval(gcdv3_pts$samples22k, nsamp_cutpts)+1
```

## Plot the maps ##

First, plot the sites with no background.  The code below creates a .pdf file in the working directory.  

```{r maps1-38, echo=TRUE, eval=FALSE}
# version with no background
pdffile <- "gcdv3_nsamp.pdf"
pdf(paste(pdffile,sep=""), paper="letter", width=8, height=8)

# basemap -- version with no background
plot(st_geometry(bb_poly_proj), col="gray90", border="black", lwd=0.1)
plot(st_geometry(grat30_lines_proj), col="gray50", lwd = 0.4, add=TRUE)
plot(st_geometry(admin0_poly_proj), border="gray50", lwd = 0.4, col = "white", add=TRUE)
plot(st_geometry(lrglakes_poly_proj), border="black", lwd = 0.2, col = "gray90", add=TRUE)
# plot(st_geometry(admin1_poly_proj), border="pink", add=TRUE)
plot(st_geometry(coast_lines_proj), col = "black", lwd = 0.4, add = TRUE)
plot(st_geometry(bb_lines_proj), border="black", lwd = 1.0, add = TRUE)

plot(gcdv3_pts_proj, pch=2, col=nsamp_colors[nsamp_num], cex=nsamp_cex[nsamp_num], lwd=0.6, add=TRUE)

text(-17000000, 9100000, pos=4, cex=0.8, "GCDv3 -- Number of Samples (Since 22 ka)")
legend(-17000000, -5000000, legend=c("< 10","10 - 100","100 - 1000","> 1000"), bg="white",
       title="Number of Samples", pch=2, pt.lwd=0.6, col=nsamp_colors, cex=0.6)

dev.off()
```

The resulting plot will look like this:  
![](images/gcdv3_nsamp.png)

Now plot the version with the version with the tree cover background data:

```{r maps1-39, echo=TRUE, eval=FALSE}
# version with treecover background
pdffile <- "gcdv3_nsamp_treecov.pdf"
pdf(pdffile, paper="letter", width=8, height=8)

# basemap -- version with treecover background
plot(st_geometry(bb_poly_proj), col="aliceblue", border="black", lwd=0.1)
plot(st_geometry(grat30_lines_proj), col="gray50", lwd = 0.4, add=TRUE)
plot(treecov_poly_proj, col=clr, bor=clr, lwd=0.01, ljoin="bevel", add=TRUE)
plot(st_geometry(admin0_lines_proj), col="gray50", lwd = 0.4, add=TRUE)
plot(st_geometry(lrglakes_lines_proj), col="gray50", bor = "black", lwd = 0.2, add=TRUE)
plot(st_geometry(caspian_poly_proj), col="aliceblue", bor = "black", lwd=0.2, add=TRUE)
# plot(st_geometry(admin1_poly_proj), border="pink", add=TRUE)
plot(st_geometry(coast_lines_proj), col = "black", lwd = 0.4, add = TRUE)
plot(st_geometry(bb_lines_proj), border="black", lwd = 0.5, add = TRUE)

plot(gcdv3_pts_proj, pch=2, col=nsamp_colors[nsamp_num], cex=nsamp_cex[nsamp_num], lwd=0.6, add=TRUE)

text(-17000000, 9100000, pos=4, cex=0.8, "GCDv3 -- Number of Samples (Since 22 ka)")
legend(-17000000, -5000000, legend=c("< 10","10 - 100","100 - 1000","> 1000"), bg="white",
       title="Number of Samples", pch=2, pt.lwd=0.6, col=nsamp_colors, cex=0.6)
dev.off()
```

The resulting plot will look like this:

![](images/gcdv3_nsamp_treecov.png)


Note the ocean-colored Caspian knockout.

[[Back to top]](Rmaps.html)