1.1 Read the Natural Earth shapefiles

Load the appropriate packages.

library(maptools)

## Loading required package: sp

## Checking rgeos availability: TRUE

library(sp)
library(rgdal)

## rgdal: version: 1.4-3, (SVN revision 828)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1

library(raster)
library(rgeos)

## rgeos version: 0.4-3, (SVN revision 595)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE

library(lattice)

Set the shapefile names, including those for global coastlines, adminstrative units (borders). Set the filenames:

# Natural Earth shape files -- global (Robinson) projections
# get shapefiles from http://www.naturalearthdata.com
shape_path <- "/Users/bartlein/Projects/ESSD/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="")

Read and plot the shapefiles (note: summary output is suppressed)

layer <- ogrListLayers(coast_shapefile)
ogrInfo(coast_shapefile, layer=layer)
coast_lines <- readOGR(coast_shapefile, layer=layer)
summary(coast_lines)  
plot(coast_lines)

layer <- ogrListLayers(admin0_shapefile)
ogrInfo(admin0_shapefile, layer=layer)
admin0_poly <- readOGR(admin0_shapefile, layer=layer)
summary(admin0_poly)  
plot(admin0_poly, bor="gray", add=TRUE)

layer <- ogrListLayers(admin1_shapefile)
ogrInfo(admin1_shapefile, layer=layer)
admin1_poly <- readOGR(admin1_shapefile, layer=layer)
summary(admin1_poly)  
plot(admin1_poly, bor="lightgreen", add=TRUE)

layer <- ogrListLayers(lakes_shapefile)
ogrInfo(lakes_shapefile, layer=layer)
lakes_poly <- readOGR(lakes_shapefile, layer=layer)
summary(lakes_poly)  
plot(lakes_poly, bor="lightblue", add=TRUE)

lrglakes_poly <- lakes_poly[as.numeric(lakes_poly$scalerank) <= 2 ,]
plot(lrglakes_poly, bor="purple", add=TRUE)

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

head(admin0_poly)

##   scalerank      featurecla LABELRANK     SOVEREIGNT SOV_A3 ADM0_DIF LEVEL              TYPE       ADMIN
## 0         3 Admin-0 country         5    Netherlands    NL1        1     2           Country       Aruba
## 1         0 Admin-0 country         3    Afghanistan    AFG        0     2 Sovereign country Afghanistan
## 2         0 Admin-0 country         3         Angola    AGO        0     2 Sovereign country      Angola
## 3         3 Admin-0 country         6 United Kingdom    GB1        1     2        Dependency    Anguilla
## 4         0 Admin-0 country         6        Albania    ALB        0     2 Sovereign country     Albania
## 5         3 Admin-0 country         6        Finland    FI1        1     2           Country       Aland
##   ADM0_A3 GEOU_DIF     GEOUNIT GU_A3 SU_DIF     SUBUNIT SU_A3 BRK_DIFF        NAME     NAME_LONG BRK_A3
## 0     ABW        0       Aruba   ABW      0       Aruba   ABW        0       Aruba         Aruba    ABW
## 1     AFG        0 Afghanistan   AFG      0 Afghanistan   AFG        0 Afghanistan   Afghanistan    AFG
## 2     AGO        0      Angola   AGO      0      Angola   AGO        0      Angola        Angola    AGO
## 3     AIA        0    Anguilla   AIA      0    Anguilla   AIA        0    Anguilla      Anguilla    AIA
## 4     ALB        0     Albania   ALB      0     Albania   ALB        0     Albania       Albania    ALB
## 5     ALD        0       Aland   ALD      0       Aland   ALD        0       Aland Aland Islands    ALD
##      BRK_NAME BRK_GROUP ABBREV POSTAL                    FORMAL_EN FORMAL_FR NOTE_ADM0 NOTE_BRK
## 0       Aruba      <NA>  Aruba     AW                        Aruba      <NA>     Neth.     <NA>
## 1 Afghanistan      <NA>   Afg.     AF Islamic State of Afghanistan      <NA>      <NA>     <NA>
## 2      Angola      <NA>   Ang.     AO  People's Republic of Angola      <NA>      <NA>     <NA>
## 3    Anguilla      <NA>   Ang.     AI                         <NA>      <NA>      U.K.     <NA>
## 4     Albania      <NA>   Alb.     AL          Republic of Albania      <NA>      <NA>     <NA>
## 5       Aland      <NA>  Aland     AI                Åland Islands      <NA>      Fin.     <NA>
##     NAME_SORT NAME_ALT MAPCOLOR7 MAPCOLOR8 MAPCOLOR9 MAPCOLOR13  POP_EST GDP_MD_EST POP_YEAR LASTCENSUS
## 0       Aruba     <NA>         4         2         2          9   103065     2258.0      -99       2010
## 1 Afghanistan     <NA>         5         6         8          7 28400000    22270.0      -99       1979
## 2      Angola     <NA>         3         2         6          1 12799293   110300.0      -99       1970
## 3    Anguilla     <NA>         6         6         6          3    14436      108.9      -99        -99
## 4     Albania     <NA>         1         4         1          6  3639453    21810.0      -99       2001
## 5       Aland     <NA>         4         1         4          6    27153     1563.0      -99        -99
##   GDP_YEAR                    ECONOMY              INCOME_GRP WIKIPEDIA FIPS_10_ ISO_A2 ISO_A3 ISO_N3
## 0      -99       6. Developing region 2. High income: nonOECD       -99       AA     AW    ABW    533
## 1      -99  7. Least developed region           5. Low income       -99       AF     AF    AFG    004
## 2      -99  7. Least developed region  3. Upper middle income       -99       AO     AO    AGO    024
## 3      -99       6. Developing region  3. Upper middle income       -99       AV     AI    AIA    660
## 4      -99       6. Developing region  4. Lower middle income       -99       AL     AL    ALB    008
## 5      -99 2. Developed region: nonG7    1. High income: OECD       -99      -99     AX    ALA    248
##   UN_A3 WB_A2 WB_A3   WOE_ID WOE_ID_EH                   WOE_NOTE ADM0_A3_IS ADM0_A3_US ADM0_A3_UN
## 0   533    AW   ABW 23424736  23424736 Exact WOE match as country        ABW        ABW        -99
## 1   004    AF   AFG 23424739  23424739 Exact WOE match as country        AFG        AFG        -99
## 2   024    AO   AGO 23424745  23424745 Exact WOE match as country        AGO        AGO        -99
## 3   660   -99   -99 23424751  23424751 Exact WOE match as country        AIA        AIA        -99
## 4   008    AL   ALB 23424742  23424742 Exact WOE match as country        ALB        ALB        -99
## 5   248   -99   -99 12577865  12577865 Exact WOE match as country        ALA        ALD        -99
##   ADM0_A3_WB     CONTINENT REGION_UN       SUBREGION                 REGION_WB NAME_LEN LONG_LEN
## 0        -99 North America  Americas       Caribbean Latin America & Caribbean        5        5
## 1        -99          Asia      Asia   Southern Asia                South Asia       11       11
## 2        -99        Africa    Africa   Middle Africa        Sub-Saharan Africa        6        6
## 3        -99 North America  Americas       Caribbean Latin America & Caribbean        8        8
## 4        -99        Europe    Europe Southern Europe     Europe & Central Asia        7        7
## 5        -99        Europe    Europe Northern Europe     Europe & Central Asia        5       13
##   ABBREV_LEN TINY HOMEPART
## 0          5    4      -99
## 1          4  -99        1
## 2          4  -99        1
## 3          4  -99      -99
## 4          4  -99        1
## 5          5    5      -99

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.

can_poly <- admin1_poly[admin1_poly$sov_a3 == "CAN" ,]
us_poly <- admin1_poly[admin1_poly$sov_a3 == "US1",]
plot(coast_lines)
plot(can_poly, bor="red", add=TRUE)
plot(us_poly, bor="blue", add=TRUE)

Convert the U.S. and Canada polygons to SpatialLines:

can_lines <- as(can_poly, "SpatialLines")
us_lines <- as(us_poly, "SpatialLines")

1.2 Project the shape files

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

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 = CRS(na_proj4string)

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

coast_lines_proj <-spTransform(coast_lines, na_crs)
admin0_poly_proj <-spTransform(admin0_poly, na_crs)
admin1_poly_proj <-spTransform(admin1_poly, na_crs)
lakes_poly_proj <-spTransform(lakes_poly, na_crs)
lrglakes_poly_proj <-spTransform(lrglakes_poly, na_crs)
can_poly_proj <-spTransform(can_poly, na_crs)
us_poly_proj <-spTransform(us_poly, na_crs)
can_lines_proj <-spTransform(can_lines, na_crs)
us_lines_proj <-spTransform(us_lines, na_crs)

plot(coast_lines_proj)
plot(admin0_poly_proj, bor="gray", add=TRUE)
plot(coast_lines_proj, 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]

na10km_bb <- as(extent(-5770000,5000000,-4510000,4480000), "SpatialPolygons")
proj4string(na10km_bb) <- na_proj4string

na10km_coast_lines_proj <- gIntersection(coast_lines_proj, na10km_bb)
na10km_lakes_poly_proj <- gIntersection(lakes_poly_proj, na10km_bb)
na10km_lrglakes_poly_proj <- gIntersection(lrglakes_poly_proj, na10km_bb)
na10km_can_poly_proj <- gIntersection(can_poly_proj, na10km_bb)
na10km_us_poly_proj <- gIntersection(us_poly_proj, na10km_bb)
na10km_can_lines_proj <- gIntersection(can_lines_proj, na10km_bb)
na10km_us_lines_proj <- gIntersection(us_lines_proj, na10km_bb)

Now do the admin shape files:

na10km_bb <- gBuffer(na10km_bb, byid=TRUE, width=0)
admin0_poly_proj <- gSimplify(admin0_poly_proj, tol = 0.00001)
na10km_admin0_poly_proj <- gBuffer(admin0_poly_proj, byid=TRUE, width=0)
na10km_admin0_poly_proj <- gIntersection(admin0_poly_proj, byid=TRUE, na10km_bb)
admin1_poly_proj <- gSimplify(admin1_poly_proj, tol = 0.00001)
na10km_admin1_poly_proj <- gBuffer(admin1_poly_proj, byid=TRUE, width=0)
na10km_admin1_poly_proj <- gIntersection(admin1_poly_proj, byid=TRUE, na10km_bb)

Plot the projected shapefiles.

plot(na10km_coast_lines_proj)
plot(na10km_admin0_poly_proj, bor="gray", add=TRUE)
plot(na10km_can_lines_proj, col="pink", add=TRUE)
plot(na10km_us_lines_proj, col="lightblue", add=TRUE)
plot(na10km_lrglakes_poly_proj, bor="blue", add=TRUE)
plot(na10km_bb, bor="purple", add=TRUE)

1.3 Write out the projected and trimmed shape files

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

# write out the various shapefiles 
outpath <- "/Users/bartlein/Projects/ESSD/data/RMaps/derived/na10km_10m/"
outshape <- na10km_coast_lines_proj
outfile <- "na10km_10m_coast_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

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

test <- readOGR(outshapefile)

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/bartlein/Projects/ESSD/data/RMaps/derived/na10km_10m/na10km_10m_coast_lines", layer: "na10km_10m_coast_lines"
## with 1 features
## It has 1 fields

plot(test, col="black")

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

# write out the various shapefiles 
outshape <- na10km_bb
outfile <- "na10km_10m_bb"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

outshape <- na10km_lakes_poly_proj
outfile <- "na10km_10m_lakes_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

outshape <- na10km_lrglakes_poly_proj
outfile <- "na10km_10m_lrglakes_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

outshape <- na10km_can_poly_proj
outfile <- "na10km_10m_can_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

outshape <- na10km_us_poly_proj
outfile <- "na10km_10m_us_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

outshape <- na10km_can_lines_proj
outfile <- "na10km_10m_can_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

outshape <- na10km_us_lines_proj
outfile <- "na10km_10m_us_lines"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialLinesDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

outshape <- na10km_admin0_poly_proj
outfile <- "na10km_10m_admin0_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

outshape <- na10km_admin1_poly_proj
outfile <- "na10km_10m_admin1_poly"
outshapefile <- paste(outpath,outfile,sep="")
spdf <- data.frame(as.numeric(row.names(outshape)))
row.names(spdf) <- row.names(outshape)
outshape <- SpatialPolygonsDataFrame(outshape, spdf)
writeOGR(outshape, outshapefile, outfile, driver="ESRI Shapefile", overwrite_layer=TRUE)

[Back to top]

2.1 Read the read the projected and trimmed shapefiles

Load the appropriate packages.

library(sp)

Read the shapefiles:

shapepath <- "/Users/bartlein/Projects/ESSD/data/RMaps/derived/na10km_10m/"
na10km_bb <- readOGR(paste(shapepath, "na10km_10m_bb.shp", sep=""))
na10km_coast_lines_proj <- readOGR(paste(shapepath, "na10km_10m_coast_lines.shp", sep=""))
na10km_admin0_poly_proj <- readOGR(paste(shapepath, "na10km_10m_admin0_poly.shp", sep=""))
na10km_lakes_poly_proj <- readOGR(paste(shapepath, "na10km_10m_lakes_poly.shp", sep=""))
na10km_lrglakes_poly_proj <- readOGR(paste(shapepath, "na10km_10m_lrglakes_poly.shp", sep=""))
na10km_can_lines_proj <- readOGR(paste(shapepath, "na10km_10m_can_lines.shp", sep=""))
na10km_us_lines_proj <- readOGR(paste(shapepath, "na10km_10m_us_lines.shp", sep=""))

Plot the projected and trimmed shapefiles:

plot(na10km_admin0_poly_proj, bor="gray")
plot(na10km_can_lines_proj, col="gray", add=TRUE)
plot(na10km_us_lines_proj, col="gray", add=TRUE)
plot(na10km_lrglakes_poly_proj, col="lightblue", add=TRUE)
plot(na10km_coast_lines_proj, col="black", add=TRUE)
plot(na10km_bb, bor="black", add=TRUE)

2.2 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 hillshade function from the raster package. Note that in this file, the coordinates are in km, and so they must be multiplied by 1000.

datapath <- "/Users/bartlein/Projects/ESSD/data/csv_files/"
datafile <- "na10km_shade.csv"
shade <- read.csv(paste(datapath,datafile,sep=""))
shade$x <- shade$x*1000
shade$y <- shade$y*1000
head(shade)

##          x       y    shade      lat       lon
## 1 -5770000 -820000 0.691208 21.83661 -160.1851
## 2 -5760000 -820000 0.691208 21.90175 -160.1023
## 3 -5720000 -860000 0.310095 21.90119 -159.5651
## 4 -5720000 -850000 0.333957 21.96637 -159.6163
## 5 -5720000 -840000 0.559432 22.03148 -159.6676
## 6 -5720000 -830000 0.956591 22.09651 -159.7191

Convert the dataframe to a SpatialPixelsData Frame

coordinates(shade) <- c("x","y")
points2grid(shade)

## Warning in points2grid(shade): grid has empty column/rows in dimension 1

##                          x        y
## cellcentre.offset -5770000 -4510000
## cellsize             10000    10000
## cells.dim             1078      900

shade_pixels <- as(shade, "SpatialPixelsDataFrame")

## Warning in points2grid(points, tolerance, round): grid has empty column/rows in dimension 1

summary(shade_pixels)

## Object of class SpatialPixelsDataFrame
## Coordinates:
##        min     max
## x -5775000 5005000
## y -4515000 4485000
## Is projected: NA 
## proj4string : [NA]
## Number of points: 282017
## Grid attributes:
##   cellcentre.offset cellsize cells.dim
## x          -5770000    10000      1078
## y          -4510000    10000       900
## Data attributes:
##      shade              lat              lon          
##  Min.   :-0.5064   Min.   :-2.757   Min.   :-203.188  
##  1st Qu.: 0.6297   1st Qu.:35.292   1st Qu.:-113.434  
##  Median : 0.6958   Median :49.353   Median : -96.309  
##  Mean   : 0.6556   Mean   :46.929   Mean   : -96.031  
##  3rd Qu.: 0.7318   3rd Qu.:62.580   3rd Qu.: -74.614  
##  Max.   : 1.0000   Max.   :83.617   Max.   :   4.453

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

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:

shade_int <- as.integer(((shade$shade+1)/2)*40)+1
shade_colnum <- shade_clr[shade_int]

2.3 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 ofthe map, and using the locate() function. The cex=0.09 argument in the points() function was detrmined by trial and error.

pdf(file = "na_shade01b.pdf")
plot(na10km_bb, col="gray95")
points(shade_pixels, pch=15, cex=0.09, col=shade_colnum)
plot(na10km_admin0_poly_proj, lwd=0.2, bor="gray50", add=TRUE)
plot(na10km_can_lines_proj, lwd=0.2, col="gray50", add=TRUE)
plot(na10km_us_lines_proj, lwd=0.2, col="gray50", add=TRUE)
plot(na10km_lrglakes_poly_proj, lwd=0.2, bor="black", col="gray90", add=TRUE)
plot(na10km_coast_lines_proj, lwd=0.3, add=TRUE)
text(-5770000, 4620000, pos=c(4), offset=0.0, cex=1.0, "NA10km_v2 -- 10m Natural Earth Outlines")
plot(na10km_bb, add=TRUE)
dev.off()

The resulting plot will look like this:

[Back to top]