1 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. As before, Natural Earth shapefiles are read and projected, this time using a Lambert Azimuthal Equal-Area projection, and trimmed to the appropriate region.

1.1 Read the Natural Earth shapefiles

Load the appropriate packages.

## Loading required package: sp
## Checking rgeos availability: TRUE
## 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
## 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

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

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

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

##   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.

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

1.2 Project the shape files

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

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

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]

Now do the admin shape files:

Plot the projected shapefiles.

1.3 Write out the projected and trimmed shape files

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

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

## 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

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 Map of North American shaded relief

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.

##          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

## 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
## Warning in points2grid(points, tolerance, round): grid has empty column/rows in dimension 1
## 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.

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