maps2.Rmd

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

```{r maps2-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>

# A second example #

This second example illustrates the creating of a base map for North America that conforms to the projection used for the `na10km_v2` data, a data set of present-day climate that accompanies a packrat-midden database: Strickland, L.E., Thompson, R.S., Shafer, S.L., Bartlein, P.J., Pelltier, R.T., Anderson, K.H., Schumann, R.R., McFadden, A.K., 2024. Plant macrofossil data for 48-0 ka in the USGS North American Packrat Midden Database, version 5.0. Scientific Data 11, 68.
[[https://www.nature.com/articles/s41597-023-02616-y]](https://www.nature.com/articles/s41597-023-02616-y).  As before, *Natural Earth* shapefiles are read and projected, this time using a Lambert Azimuthal Equal-Area projection, and trimmed to the appropriate region.

## Read the Natural Earth shapefiles ##

Load the appropriate packages.

```{r maps2-2, messages=FALSE, results="hide"}
library(sf)
library(stars)
library(ggplot2)
library(RColorBrewer)
library(classInt)
```

Set the shapefile names, including those for global coastlines, adminstrative units (borders).  (Note that these are smaller-spatial/larger-cartrographic resolution shapefiles, more appropriate for continental-scale as opposed to global-scale mapping.)  Set the filenames:

```{r maps2-3 }
# read Natural Earth shapefiles
shape_path <- "/Users/bartlein/Projects/RESS/data/RMaps/source/"
coast_shapefile <- paste(shape_path, "ne_10m_coastline/ne_10m_coastline.shp", sep="")
admin0_shapefile <- paste(shape_path, "ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp", sep="")
admin1_shapefile <- paste(shape_path, 
        "ne_10m_admin_1_states_provinces_lakes/ne_10m_admin_1_states_provinces_lakes.shp", sep="")
lakes_shapefile <- paste(shape_path, "ne_10m_lakes/ne_10m_lakes.shp", sep="")
grat15_shapefile <- paste(shape_path, "ne_10m_graticules_all/ne_10m_graticules_15.shp", sep="")
```

As before, query geometry types to compose object names:

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

Read and plot the shapefiles. At a higher spatial resolution, these files will take more time to load. (Note:  summary output is suppressed)

```{r maps2-5, results="hide", warnings=FALSE, cache=TRUE}
# read shapefiles
coast_lines_sf <- st_read(coast_shapefile) # note geometry type MULTILINESTRING
plot(st_geometry(coast_lines_sf), col="gray50")
admin0_poly_sf <- st_read(admin0_shapefile) # note:  geometry type MULTIPOLYGON
plot(st_geometry(admin0_poly_sf), col="gray70", border="red", add = TRUE)
admin1_poly_sf <- st_read(admin1_shapefile) # note:  geometry type MULTIPOLYGON
plot(st_geometry(admin1_poly_sf), col="gray70", border="pink", add = TRUE)
lakes_poly_sf <- st_read(lakes_shapefile) # note:  geometry type POLYGON
plot(st_geometry(lakes_poly_sf), col="blue", add = TRUE)
grat15_lines_sf <- st_read(grat15_shapefile) # note:  geometry type LINESTRING
plot(st_geometry(grat15_lines_sf), col="gray50", add = TRUE)

# filter out the small lakes
lrglakes_poly_sf <- lakes_poly_sf[as.numeric(lakes_poly_sf$scalerank) <= 2 ,]
plot(lrglakes_poly_sf, bor="purple", add=TRUE)
```


Take a look at the `admin1_poly` dataframe, to figure out the codes for Canadian and U.S. provincial and state borders.

```{r maps2-6 }
head(admin0_poly_sf)
```

The approprate code to extract the U.S. and Canada data is `admin1_poly$sr_sov_a3 == "CAN"` and `admin1_poly$sr_sov_a3 == "US1"`.  Extract the borders, and plot the resulting shapefiles.

```{r maps2-7 }
# extract US and Canada boundaries
can_poly_sf <- admin1_poly_sf[admin1_poly_sf$sov_a3 == "CAN" ,]
us_poly_sf <- admin1_poly_sf[admin1_poly_sf$sov_a3 == "US1",]
plot(st_geometry(coast_lines_sf))
plot(st_geometry(can_poly_sf), border = "red", add=TRUE)
plot(st_geometry(us_poly_sf), border = "blue", add=TRUE)
```

Convert the U.S. and Canada polygons to `MULTILINESTRINGS`:

```{r maps2-8 }
# convert geometries from polygons to lines
admin0_lines_sf <- st_cast(admin0_poly_sf, "MULTILINESTRING")
can_lines_sf <- st_cast(can_poly_sf, "MULTILINESTRING")
us_lines_sf <- st_cast(us_poly_sf, "MULTILINESTRING")
```

## Project the shape files ##

Set the `proj4string` value and the coordinate reference system for the na10km_v2 grid:

```{r maps2-9 }
# Lambert Azimuthal Equal Area
na_proj4string <- "+proj=laea +lon_0=-100 +lat_0=50 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
na_crs = st_crs(na_proj4string)
```

Project the various shapefiles (and plot the coastlines as an example):

```{r maps2-10 }
coast_lines_proj <- st_transform(coast_lines_sf, crs = na_crs)
admin0_poly_proj <- st_transform(admin0_poly_sf, crs = na_crs)
admin0_lines_proj <- st_transform(admin0_lines_sf, crs = na_crs)
admin1_poly_proj <- st_transform(admin1_poly_sf, crs = na_crs)
lakes_poly_proj <- st_transform(lakes_poly_sf, crs = na_crs)
lrglakes_poly_proj <- st_transform(lrglakes_poly_sf, crs = na_crs)
can_poly_proj <- st_transform(can_poly_sf, crs = na_crs)
us_poly_proj <- st_transform(us_poly_sf, crs = na_crs)
can_lines_proj <- st_transform(can_lines_sf, crs = na_crs)
us_lines_proj <- st_transform(us_lines_sf, crs = na_crs)
grat15_lines_proj <- st_transform(grat15_lines_sf, crs = na_crs)
```

```{r maps2-11, cache=TRUE}
plot(st_geometry(admin0_lines_proj), col = "gray70")
plot(st_geometry(coast_lines_proj), col = "black", add=TRUE)
plot(st_geometry(grat15_lines_proj), col = "gray70", add=TRUE)
```

Define a bounding box for trimming the polygon and line shape files to the area covered by the na10km_v2 grid. The extent of the area is known from the definition of the grid, but could also be determined by reading an na10km_v2 netCDF file.  The projected `admin` shape files are quite complicated, and create "topology exception errors".  These can be fixed using an approach discussed on StackExchange [[link]](http://gis.stackexchange.com/questions/163445/r-solution-for-topologyexception-input-geom-1-is-invalid-self-intersection-er)

```{r maps2-12 }
# North America Bounding box (coords in metres)
na10km_bb <- st_bbox(c(xmin = -5770000, ymin = -4510000, xmax = 5000000, ymax = 4480000), crs = na_proj4string)
na10km_bb <- st_as_sfc(na10km_bb)
```

Clip or trim the coastlines to the `na10km_bb` boundary-box object, using the `st_intersection()` function:

```{r maps2-13 }
# clip the projected objects to the na10km_bb
na10km_coast_lines_proj <- st_intersection(coast_lines_proj, na10km_bb)
plot(st_geometry(na10km_coast_lines_proj))
plot(st_geometry(na10km_bb), add = TRUE)
```


```{r maps2-14, warning=FALSE}
# trim the other projected objects
na10km_lakes_poly_proj <- st_intersection(lakes_poly_proj, na10km_bb)
na10km_lrglakes_poly_proj <- st_intersection(lrglakes_poly_proj, na10km_bb)
na10km_can_poly_proj <- st_intersection(can_poly_proj, na10km_bb)
na10km_us_poly_proj <- st_intersection(us_poly_proj, na10km_bb)
na10km_can_lines_proj <- st_intersection(can_lines_proj, na10km_bb)
na10km_us_lines_proj <- st_intersection(us_lines_proj, na10km_bb)
na10km_grat15_lines_proj <- st_intersection(grat15_lines_proj, na10km_bb)
na10km_admin0_poly_proj <- st_intersection(admin0_poly_proj, na10km_bb)
na10km_admin1_poly_proj <- st_intersection(admin1_poly_proj, na10km_bb)
```

Plot the projected shapefiles.

```{r maps2-15 }
# plot the projected objects
plot(st_geometry(na10km_coast_lines_proj))
plot(st_geometry(na10km_can_lines_proj), col = "pink", add = TRUE)
plot(st_geometry(na10km_us_lines_proj), col = "lightblue", add = TRUE)
plot(st_geometry(na10km_lrglakes_poly_proj), col = "lightblue", bor = "blue", add=TRUE)
plot(st_geometry(na10km_admin0_poly_proj), bor = "gray70", add = TRUE)
plot(st_geometry(na10km_grat15_lines_proj), col = "gray70", add = TRUE)
plot(st_geometry(na10km_bb), border = "purple", add = TRUE)
```

## Write out the projected and trimmed shape files ##

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

```{r maps2-16, warning=FALSE}
# write out the various shapefiles 
outpath <- "/Users/bartlein/Projects/RESS/data/RMaps/derived/na10km_10m/"
```

```{r maps2-17 }
# coast lines 
outshape <- na10km_coast_lines_proj
outfile <- "na10km_10m_coast_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```

It's always good practice to test whether the shapefile has ideed been written out correctly.  Read it back in and plot it.

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

Write out the other shape files (output is suppressed):

```{r maps2-19, messages = FALSE, results="hide", warning=FALSE}
# write out the various shapefiles 
outshape <- na10km_bb
outfile <- "na10km_10m_bb/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# lakes poly
outshape <- na10km_lakes_poly_proj
outfile <- "na10km_10m_lakes_poly/"
outshapefile <- paste(outpath, outfile, sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# lrglakes poly
outshape <- na10km_lrglakes_poly_proj
outfile <- "na10km_10m_lrglakes_poly"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# can poly
outshape <- na10km_can_poly_proj
outfile <- "na10km_10m_can_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# us poly
outshape <- na10km_us_poly_proj
outfile <- "na10km_10m_us_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# can lines
outshape <- na10km_can_lines_proj
outfile <- "na10km_10m_can_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# us lines
outshape <- na10km_us_lines_proj
outfile <- "na10km_10m_us_lines/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# admin0 poly
outshape <- na10km_admin0_poly_proj
outfile <- "na10km_10m_admin0_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)

# admin1 poly
outshape <- na10km_admin1_poly_proj
outfile <- "na10km_10m_admin1_poly/"
outshapefile <- paste(outpath,outfile,sep="")
st_write(outshape, outshapefile, driver = "ESRI Shapefile", append = FALSE)
```
[[Back to top]](Rmaps2.html)


# Map of North American shaded relief #

## PlotØ the projected and trimmed shapefiles ##

If necessary (i.e. if starting here), read the projected shapefiles: 

```{r maps2-20, cache=TRUE}
# plot the basemap (note different order than before)
plot(st_geometry(na10km_bb), border = "black", col = "aliceblue")
plot(st_geometry(na10km_grat15_lines_proj), col = "gray70", add = TRUE)
plot(st_geometry(na10km_admin0_poly_proj), col = "white",  border = "gray", add = TRUE)
plot(st_geometry(na10km_can_lines_proj), col = "gray", add = TRUE)
plot(st_geometry(na10km_us_lines_proj), col = "gray", add = TRUE)
plot(st_geometry(na10km_lrglakes_poly_proj), col = "aliceblue", border = "black", lwd = 0.5, add = TRUE)
plot(st_geometry(na10km_coast_lines_proj), col = "black", lwd = 0.5, add = TRUE)
plot(st_geometry(na10km_bb), boder = "black", add = TRUE)
```

## Read a shaded relief file ##

Read a pre-computed shaded relief file.  This could also be crearted by reading the na10km_v2 grid-point elevations and using the `shade` function from the `terra` package.  Note that in this file, the coordinates are in km, and so they must be multiplied by 1000.

```{r maps2-21 }
# get 10km shade values
datapath <- "/Users/bartlein/Projects/RESS/data/csv_files/"
datafile <- "na10km_shade.csv"
shade <- read.csv(paste(datapath,datafile,sep=""))
str(shade)
shade$x <- shade$x*1000
shade$y <- shade$y*1000
head(shade)
```

Convert the dataframe to a POINTS `sf` object

```{r maps2-22 }
# convert to POINTS
shade_pixels <- st_as_sf(shade, coords = c("x", "y"), crs = na_crs)
shade_pixels
```

Read some predetermined (gray-scale) colors for the shading.

```{r maps2-23 }
# read hillshade color numbers
colorfile <- "shade40_clr.csv"
shade_rgb <- read.csv(paste(datapath, colorfile, sep=""))
shade_clr <- rgb(shade_rgb)
```

Set the (gray-scale) color numbers for each pixel:

```{r maps2-24 }
shade_int <- as.integer(((shade$shade+1)/2)*40)+1
shade_colnum <- shade_clr[shade_int]
```

## Make the map ##

Plot the shaded-relief colors and the various shape files.  The location of the text string was determined by plotting an initial version of the map, and using the `locator()` function.  The `cex=0.09` argument in plotting `shade_pixels` was determined by trial and error. Redirect the output to a .pdf file.

```{r maps2-25, cache = TRUE, echo=TRUE, eval=FALSE, messages=FALSE}
pdf(file = "na_shade01b.pdf")
plot(st_geometry(na10km_bb), border = "black", col = "aliceblue")
plot(st_geometry(na10km_grat15_lines_proj), col = "gray70", add = TRUE)
plot(st_geometry(shade_pixels), pch=15, cex=0.09, col=shade_colnum, add = TRUE)
plot(st_geometry(na10km_admin0_poly_proj), lwd=0.2, bor="gray50", add=TRUE)
plot(st_geometry(na10km_can_lines_proj), lwd=0.2, col="gray50", add=TRUE)
plot(st_geometry(na10km_us_lines_proj), lwd=0.2, col="gray50", add=TRUE)
plot(st_geometry(na10km_lrglakes_poly_proj), lwd=0.2, bor="black", col="aliceblue", add=TRUE)
plot(st_geometry(na10km_coast_lines_proj), lwd=0.5, add=TRUE)
text(-5770000, 4620000, pos=c(4), offset=0.0, cex=1.0, "NA10km_v2 -- 10m Natural Earth Outlines")
plot(st_geometry(na10km_bb), add=TRUE)
dev.off()
```


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

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

# Mapping the western U.S. precipitation-ratio data #

A final example describes mapping the western U.S. precipitation data via the `{ggplot2}` package, as seen in the [[Using R: Examples]](https://pjbartlein.github.io/REarthSysSci/usingR.html#ggplot2-maps). Here, the basemap will be extracted from the `{maps}` package (and later compared with the Natural Earth shapefiles basemap versions).

## Get and project basemap outlies ##

Get oulines for the western states. the `map()` function retrieves the outlines from the `maps()` package database, and `st_as_sf()`, converts those to and (unprojected) MULTIPOLYGONS `sf` object:

```{r maps2-26}
library(maps)

# get outlines of western states from {maps} package
wus_sf <- st_as_sf(map("state", regions=c("washington", "oregon", "california", "idaho", "nevada",
    "montana", "utah", "arizona", "wyoming", "colorado", "new mexico", "north dakota", "south dakota",
    "nebraska", "kansas", "oklahoma", "texas"), plot = FALSE, fill = TRUE))
head(wus_sf); tail(wus_sf)
class(wus_sf)
plot(wus_sf)
```

Here's a quick `ggplot2` map of the outlines:

```{r maps2-27}
ggplot() + geom_sf(data = wus_sf) + theme_bw()
```

Project the outlines and a graticule (extracted from the `wus_sf` object using `st_graticule()`):

```{r maps2-28}
# projection
laea = st_crs("+proj=laea +lat_0=30 +lon_0=-110") # Lambert Azimuthal Equal Area
wus_sf_proj = st_transform(wus_sf, laea)
plot(st_geometry(wus_sf_proj), axes = TRUE)
plot(st_geometry(st_graticule(wus_sf_proj)), axes = TRUE, add = TRUE)
```

Note the coordinates are now in meters.

Here's a `ggplot2` version:

```{r maps2-29}
ggplot() + geom_sf(data = wus_sf_proj) + theme_bw()
```

Get two other outlines from the `maps` database, one for all of North America, and the other for the U.S."

```{r maps2-30}
# get two other outlines
na_sf <- st_as_sf(map("world", regions=c("usa", "canada", "mexico"), plot = FALSE, fill = TRUE))
conus_sf <- st_as_sf(map("state", plot = FALSE, fill = TRUE))
```

Because `sf` objects are basically dataframes, they can be manipulated using some of the common R functions, like `rbind` (row bind). Combine the `na_sf` and `conus_sf` objects, and plot them out:

```{r maps2-31}
# na2_df merged sf's
na2_sf <- rbind(na_sf, conus_sf)
head(na2_sf)
plot(st_geometry(na2_sf))
```

Note that the map is slightly strange--the Aleutians cross the dateline (180W), so the map wraps around.

```{r maps2-32}
# plot again, with longitude limited
plot(st_geometry(na2_sf), xlim = c(-180, -40))
```

Project the additional outlines:

```{r maps2-33}
# project the other outlines
na_sf_proj = st_transform(na_sf, laea)
conus_sf_proj = st_transform(conus_sf, laea)
na2_sf_proj = st_transform(na2_sf, laea)
```

```{r maps2-34}
# plot everthing so far
plot(st_geometry(na_sf_proj), col = "red", axes = TRUE)
plot(st_geometry(wus_sf_proj), axes = TRUE, add = TRUE)
plot(st_geometry(st_graticule(wus_sf_proj)), add = TRUE)
```

## Get and project the *pratio* data ##

```{r maps2-35}
# pratios
csv_path <- "/Users/bartlein/Projects/RESS/data/csv_files/"
csv_file <- "wus_pratio.csv"
wus_pratio <- read.csv(paste(csv_path, csv_file, sep = ""))
```

Make an `sf` object:

```{r maps2-36}
# make an sf object
wus_pratio_sf <- st_as_sf(wus_pratio, coords = c("lon", "lat"))
plot(st_geometry(wus_pratio_sf), pch=16, cex = 0.5)
```

Recode the the `pjulpann` ratio to a factor:

```{r maps2-37}
# recode pjulpjan to a factor
library(classInt)
wus_pratio$pjulpjan <- wus_pratio$pjulpann/wus_pratio$pjanpann  # pann values cancel out
cutpts <- c(0.0, .100, .200, .500, .800, 1.0, 1.25, 2.0, 5.0, 10.0, 9999.0)
pjulpjan_factor <- factor(findInterval(wus_pratio$pjulpjan, cutpts))
head(cbind(wus_pratio$pjulpjan, pjulpjan_factor, cutpts[pjulpjan_factor]))
```

Here's an (unprojected) map of the precipitation ratio:

```{r maps2-38}
## ggplot2 map of pjulpjan
ggplot() +
  geom_sf(data = wus_sf, fill=NA) +
  scale_color_brewer(type = "div", palette = "PRGn", aesthetics = "color", direction = 1,
                    labels = c("0.0 - 0.1", "0.1 - 0.2", "0.2 - 0.5", "0.5 - 0.8", "0.8 - 1.0",
                                "1.0 - 1.25", "1.25 - 2.0", "2.0 - 5.0", "5.0 - 10.0", "> 10.0"),
                    name = "Jul:Jan Ppt. Ratio") +
  labs(title = "Precipitation Seasonality", x = "Longitude", y = "Latitude") +
  geom_point(aes(wus_pratio$lon, wus_pratio$lat, color = pjulpjan_factor), size = 1.5 ) +
  theme_bw()
```

Project the precipitation ratio point

```{r maps2-39}
# projection
st_crs(wus_pratio_sf) <- st_crs("+proj=longlat")
laea = st_crs("+proj=laea +lat_0=30 +lon_0=-110") # Lambert equal area
wus_pratio_sf_proj = st_transform(wus_pratio_sf, laea)
```

Plot the projected points:

```{r maps2-40}
# plot the projected points
plot(st_geometry(wus_pratio_sf_proj), pch=16, cex=0.5, axes = TRUE)
plot(st_geometry(st_graticule(wus_pratio_sf_proj)), axes = TRUE, add = TRUE)
```
Plot everything

```{r maps2-41}
# ggplot2 map of pjulpman projected
ggplot() +
  geom_sf(data = na2_sf_proj, fill=NA, color="gray") +
  geom_sf(data = wus_sf_proj, fill=NA) +
  geom_sf(data = wus_pratio_sf_proj, aes(color = pjulpjan_factor), size = 1.0 ) +
  scale_color_brewer(type = "div", palette = "PRGn", aesthetics = "color", direction = 1,
              labels = c("0.0 - 0.1", "0.1 - 0.2", "0.2 - 0.5", "0.5 - 0.8", "0.8 - 1.0",
                        "1.0 - 1.25", "1.25 - 2.0", "2.0 - 5.0", "5.0 - 10.0", "> 10.0"),
              name = "Jul:Jan Ppt. Ratio") +
              labs(x = "Longitude", y = "Latitude") +
  guides(color = guide_legend(override.aes = list(size = 3))) +
  coord_sf(crs = st_crs(wus_sf_proj), xlim = c(-1400000, 1500000), ylim = c(-400000, 2150000)) +
  labs(title = "Precipitation Seasonality", x = "Longitude", y = "Latitude") +
  theme_bw()
```

The `guides()` function expands the size (diameter) of the legend key (`size = 3`) leaving the points at their original size.

# Compare shapefiles

While we have two differnet sources of shapefiles in the environment (if all of the code chunks have been run), it would be useful to compare them. Copy the following code into a script window. Then run the script a line at a time. The first line will just set up a plotting region. Then the next three lines will plot in turn the Natural Earth outlines of the coastlines and the U.S. and Canada. Note that the outlines overlay one another *exactly* (there's no black visible). Then run the next two lines, one at a time to plot `{maps}` package outlines. Notice what happens along the coast, and in particular around Vancouver Is. and the Puget Sound. The fact that the Natural Earth outlines are still visible, as are the separate `{maps}` package outlines implies that the various shape files are not coregisterd.

```{r maps2-42, echo =TRUE, eval=FALSE}
# outlines 
# from the Natural Earth shapefile collection
plot(NULL, NULL, xlim = c(-130.0, -110), ylim = c(40.0, 50.0), 
     xlab = "Longitude", ylab = "Latitude")
plot(st_geometry(coast_lines_sf), border = "black", add = TRUE)
plot(st_geometry(can_poly_sf), border = "purple", add = TRUE)
plot(st_geometry(us_poly_sf), border = "magenta", add = TRUE)

# from the maps package
plot(st_geometry(na_sf), border = "red", add = TRUE)
plot(st_geometry(wus_sf), border = "blue", add = TRUE)
```

![](images/outlines.png)
Do the same thing for the following code, which plots the vertices along the MULTILINESTRINGS. 

```{r maps2-43, echo = TRUE, eval = FALSE}
# verticies (points) 
# from the Natural Earth shapefile collection
plot(NULL, NULL, xlim = c(-130.0, -110), ylim = c(40.0, 50.0), 
     xlab = "Longitude", ylab = "Latitude")
points(st_coordinates(coast_lines_sf), col = "black", pch = 16, cex = 0.4)
points(st_coordinates(can_poly_sf), col = "purple", pch = 16, cex = 0.5)
points(st_coordinates(us_poly_sf), col = "magenta", pch = 16, cex = 0.5)

# from the maps package
points(st_coordinates(na_sf), col = "red", pch = 16, cex = 0.5)
points(st_coordinates(wus_sf), col = "blue", pch = 16, cex = 0.5)
```

![](images/points.png)

Ugh!

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