1 Introduction

In this note, we illustrate how to add updated information to the existing map file and use it to plot maps with updated information.

2 Existing Base Map File and New Information

We use the built-in world map file, World , from the package sf(simple features), a standardized way to encode spatial vector data. The World data set has 17 variables including country names and the corresponding abbreviations, population densities, geometry (information of geo-polygon with longitude and latitude), etc. The data set has 177 countries. We will use the country names to define a primary key merge the new data (with updated information) with the existing data.

The new information to be included in the new data set is based on several data sets available at https://projectdat.s3.amazonaws.com/datasets.html#week12

The information we will add to the World data set is:

  1. Income from years 2000, 2005, 2010, 2015 and name them as inc00, inc05, inc10, inc15.

  2. Life expectancy from years 2000, 2005, 2010, 2015 and name them as life00, life05, life10, life15.

  3. Primary defined from the country names and their abbreviations.

#data(World)
World0 <- st_read(system.file("shapes/world.gpkg", package="spData"))
## Reading layer `world' from data source 
##   `/Library/Frameworks/R.framework/Versions/4.1/Resources/library/spData/shapes/world.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
inc <- read.csv("https://projectdat.s3.amazonaws.com/income_per_person.csv")
income <- data.frame(country = gsub(" ", "", inc$geo), inc00 = inc$X2000, inc05 = inc$X2005, inc10 = inc$X2010, inc15 = inc$X2015)
###
lifexp <- read.csv("https://projectdat.s3.amazonaws.com/life_expectancy_years.csv")
life.exp <- data.frame(country =  gsub(" ", "", lifexp$geo), lif00 = lifexp$X2000, lif05 = lifexp$X2005, lif10 = lifexp$X2010, lif15 = lifexp$X2015)
###
pop <- read.csv("https://projectdat.s3.amazonaws.com/population_total.csv")
popsize <- data.frame(country = gsub(" ", "", pop$geo), pop00 =pop$X2000, pop05 = pop$X2005, pop10 = pop$X2010, pop15 = pop$X2015)
###
region <- read.csv("https://projectdat.s3.amazonaws.com/countries_total.csv")
regions <- data.frame(country =  gsub(" ", "", region$name), iso_a3 = region$alpha.3)
###
IncLifeExp <- merge(income, life.exp, by = 'country')
IncLifeRegion <- merge(IncLifeExp, regions, by = 'country')
IncLifRegPop <- merge(IncLifeRegion, popsize, by = 'country')
IncLifRegPop$iso_a2 <- substr(IncLifRegPop$iso_a3, 1,2)
###
myWorld <- merge(World0, IncLifRegPop, by = 'iso_a2')

The next data set contains the geocode of the world capital cities. This information will be used to create pop-ups to include specific information in the data. The geocode of the capital city can be found at https://www.kaggle.com/nikitagrec/world-capitals-gps. I also placed a copy of the data set at https://raw.githubusercontent.com/pengdsci/sta553/main/map/WorldCapitalGeocode.csv

geocode <-read.csv("https://raw.githubusercontent.com/pengdsci/sta553/main/map/WorldCapitalGeocode.csv")
#geometry = paste('POINT (',CapitalLongitude ,',',CapitalLatitude,')')
#geocode$geometry = geometry
geocode$country <- gsub(" ", "", geocode$CountryName)
capital <- st_as_sf(geocode, coords = c("CapitalLongitude", "CapitalLatitude"), crs = 4326)
###
IncLifeRegionCap <- merge(capital, IncLifeRegion, by = 'country')
IncLifeRegCapPop <- merge(IncLifeRegionCap, popsize, by = 'country')

3 Thematics Map

library(tmap)
##
tmap_mode("view")  # "view" gives interactive map; 
#tmap_style("classic") ## tmap_style set to "classic"
## other available styles are: "white", "gray", "natural", 
## "cobalt", "col_blind", "albatross", "beaver", "bw", "watercolor"
tmap_options(bg.color = "skyblue", 
             legend.text.color = "white")
##
tm_shape(myWorld) +
      tm_polygons("lifeExp", 
                  legend.title = "Life Expectancy") +
      tm_layout(bg.color = "gray", 
                inner.margins = c(0, .02, .02, .02)) + 
tm_shape(IncLifeRegCapPop) +
      tm_symbols(col = "purple", 
                 size = "pop15", 
                 scale = .5,
                 alpha = 0.5,
                 popup.vars=c("CapitalName", "pop15", "inc00", "inc05", "inc10","inc15", "lif00","lif05","lif10", "lif15"))

4 Philadelphia Neighborhood Shapefiles

We can use the shapefile of any place to draw the base map of the place. For example, we can find the shapefile of the Philadelphia neighborhood at https://www.opendataphilly.org/dataset/covid-vaccinations/resource/473c9589-111b-43c9-a4a2-2dbe91f6dd7b?inner_span=True and draw the map of the Philadelphia neighborhood using the following map.

library(sf)
philly <- st_read("/Users/chengpeng/Downloads/covid_vaccines_by_census_tract/covid_vaccines_by_census_tract.shp")
## Reading layer `covid_vaccines_by_census_tract' from data source 
##   `/Users/chengpeng/Downloads/covid_vaccines_by_census_tract/covid_vaccines_by_census_tract.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 384 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.2803 ymin: 39.86746 xmax: -74.95575 ymax: 40.13793
## Geodetic CRS:  WGS 84
library(tmap)
tm_shape(philly) +
     tmap_options(check.and.fix = TRUE) +
      tm_polygons(border.col = "red",
                  border.alpha = 0.5) +
      tm_layout(bg.color = "skyblue",
                aes.color = c(fill = "skyblue", borders = "grey40",
  dots = "black", lines = "red", text = "black", na = "grey70"),
                inner.margins = c(0, .02, .02, .02))