Visualizing Geography and Area

Motivation

Create Maps in R
Associate Data to Maps

What are maps and how can WE create them

Abstract representation of an area, can be 2-D, 2D+ or 3-D.
Software exists to represent geographical areas with 2-D or 2-D that is made to look like 3-D.
R has multiple packages that allow you to visualize geographical areas.

First Known World Map: Imago Mundi

From Babylonia (close to Baghdad), 6th to 7th Century BC.
First known map of the known world surrounded by “bitter” water.
No scale, obvious mistakes, with cuneiform text.

First Known Map: Imago Mundi

Pre Modern Maps

Maps are not always represented in shapes that are similar to reality. Abstraction takes may take geometric forms or physical form in order to communicate additional information to spatial.
Pre Colombian Mapmaking Andean region:
When creating maps, motivation is not always navigation or precision.

Pre Modern Maps

Ptolemy 150 AD representing the Mediterranean(only 13th century edition survives)

[ By http://en.wikipedia.org/wiki/User%3APortolanero/Wikimedia Commons, Public Domain ( http://en.wikipedia.org/wiki/User%3APortolanero/Wikimedia Commons, Public Domain, https://commons.wikimedia.org/w/index.php?curid=16043748)

Modern Maps

Navigation and exploration.
- Scale of geographical objects.
- Precision in relation to other geographical objects.
GIS systems use modern mapping systems to add information.

Leading to Modern Maps

Problem with 2-D maps: How to represent a 3 D structure on a 2 D surface. “The most well-known map projection is the Mercator projection. This map projection has the property of being conformal. However, it has been criticized throughout the 20th century for enlarging regions further from the equator. To contrast, equal-area projections such as the Sinusoidal projection and the Gall–Peters projection show the correct sizes of countries relative to each other, but distort angles. The National Geographic Society and most atlases favor map projections that compromise between area and angular distortion, such as the Robinson projection and the Winkel tripel projection.” Wikipedia

Leading to Modern Maps

Latitude, Longitude:
- Latitude was easier to measure since earth does not move along what we refer to north-south axis, position relative to stars becomes practical to measure on the sea/land.
Longitude on the sea was tough to measure. No constant in the sky since the Earth is turning around its axis and no fixed geographical point such as a mountain to be used as a fixed point.

Solution to Longitude Problem

The solution was invented by John Harrison
The longitude could be measured by the time difference between the 12 a clock local time and the time at the port of origin since by 19th century, the speed that the earth traveled around its own axis was known for at least two thousand years. A difference of one hour in local solar time corresponds to 15° of longitude. The problem was to know what the time was at the port of origin while at sea.
Harrison invented special clocks that allow you to measure precise time on a ship.
Awarded more than 20,000 pounds in 19th century.

Digital Maps

We will learn how to map out USA geographical features, ofcourse the principles we discuss apply in general.
The files we will use is stored by CENSUS

TIGER Shape Files

CENSUS The web interface contains a multitude of border options (states, counties, ZCTA, etc…)
I choose to download counties. It is a zip file so extract it all.

TIGER Shape Files

The history of shape files held by census is an interesting one to read if you ever work with Geographical Information Systems.
tl_2025_us_county.shp
Each file in the extracted folder has a purpose.
TIGER shape files contain geographic boundaries as geometric objects defined by coordinates.

TIGER Shape Files

There are 4 important files in the folder.
- shp
- shx
- dbf
- prj
It is possible to to read and draw shapes containing geometries using shp and shx files. To create a map with a coordinate system information you need the prj file and dbf adds information.

TIGER Shape Files

shp file contains coordinates that can form lines, polygons and other shapes for the geographical unit (for instance counties).
shx file contains an index. It tells R, where each geographic feature’s geometry begins in the .shp file and how long that geometry record is.
dbf file contains attributes associated with each geographical unit (how many square miles is land and water etc…)
prj file contains the coordinate system information (for the file NAD83).

Reading in shp

# you need to have installed sf package
library(sf)
### Working Directory in my system specifying where the file is
### note that the name tl_2025_us_county is part of the file path and the ### name of the file itself.

wd<-'C:/Users/rm84/OneDrive - Texas State University/Desktop/teaching/3339/Maps_files/tl_2025_us_county'
setwd(wd)
### Read in the shp file but note that the st_read function uses shx dbf ### and the prj file in addition to the shp file.

counties <- st_read("tl_2025_us_county.shp")

Reading layer `tl_2025_us_county' from data source 
  `C:\Users\rm84\OneDrive - Texas State University\Desktop\teaching\3339\Maps_files\tl_2025_us_county\tl_2025_us_county.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 3235 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.43979
Geodetic CRS:  NAD83

Reading in the shape file

head(counties)

Simple feature collection with 6 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -104.9808 ymin: 26.29927 xmax: -79.13832 ymax: 44.19608
Geodetic CRS:  NAD83
  STATEFP COUNTYFP COUNTYNS GEOID        GEOIDFQ     NAME        NAMELSAD LSAD
1      40      075 01101825 40075 0500000US40075    Kiowa    Kiowa County   06
2      46      079 01265776 46079 0500000US46079     Lake     Lake County   06
3      37      033 01008542 37033 0500000US37033  Caswell  Caswell County   06
4      48      377 01383974 48377 0500000US48377 Presidio Presidio County   06
5      39      057 01074041 39057 0500000US39057   Greene   Greene County   06
6      48      489 01384030 48489 0500000US48489  Willacy  Willacy County   06
  CLASSFP MTFCC CSAFP CBSAFP METDIVFP FUNCSTAT      ALAND    AWATER    INTPTLAT
1      H1 G4020  <NA>   <NA>     <NA>        A 2629039892  40296743 +34.9214893
2      H1 G4020  <NA>   <NA>     <NA>        A 1457916151  31746795 +44.0284497
3      H1 G4020  <NA>   <NA>     <NA>        A 1102042927   8293623 +36.3943252
4      H1 G4020  <NA>   <NA>     <NA>        A 9985057447   1773188 +30.0058912
5      H1 G4020   212  19430     <NA>        A 1071302625   6798109 +39.6874785
6      H1 G4020   154  39700     <NA>        A 1529680719 502856285 +26.4818607
      INTPTLON                       geometry
1 -098.9816168 MULTIPOLYGON (((-98.95506 3...
2 -097.1232229 MULTIPOLYGON (((-96.88886 4...
3 -079.3396193 MULTIPOLYGON (((-79.14343 3...
4 -104.2616192 MULTIPOLYGON (((-104.9808 3...
5 -083.8948943 MULTIPOLYGON (((-84.10668 3...
6 -097.5947313 MULTIPOLYGON (((-97.86187 2...

Reading in the shape file

counties$geometry[1]

Geometry set for 1 feature 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -99.40929 ymin: 34.59437 xmax: -98.61982 ymax: 35.12418
Geodetic CRS:  NAD83

Code Basic Plot

map_1=plot(st_geometry(counties))

Code Basic Plot

NULL

Code Basic GGPlot

##ggplot should already been installed
library(ggplot2)

map_2=ggplot(counties) +
  geom_sf(fill = "lightgray", color = "black")

Map

Coordinates for GGPlot Lower 48

##ggplot should already been installed

map_3=ggplot(counties) +
  geom_sf() +
  coord_sf(
    xlim = c(-125, -66),
    ylim = c(24, 50)
  )

Map for GGPlot Lower 48

Filtering the data for GGPlot

# dplyr should already be installed
library(dplyr)
# |> is referred to as the piping operator. It allows you not to have to create an explicitly named intermediate object. 
# filter chooses data based on the logical arguments within the parenthesis
# ! turns the logical evaluation from TRUE to FALSE and from FALSE to TRUE %in% evaluates STATEFP values in the counties shape file as either TRUE or FALSE.
lower48 <- counties |>
  filter(!STATEFP %in% c(
    "02", # Alaska
    "15", # Hawaii
    "60", # American Samoa
    "66", # Guam
    "69", # Northern Mariana Islands
    "72", # Puerto Rico
    "78"  # U.S. Virgin Islands
  ))

map_4=ggplot(lower48) +
  geom_sf()

Filtering the data for GGPlot

Unique Id is GEOID:

# Data dictionary would give you which column has the unique ID
g_id=lower48$GEOID

# Arbitrarily generate random numbers to simulate values of interest
randoms=rnorm(n=length(g_id),mean=10,sd=3)
#binding two columns to create a data frame, both columns are char.
g_id=as.data.frame(cbind(g_id,randoms))
# convert the second column back to numeric.
g_id[,2]=as.numeric(g_id[,2])

names(g_id)=c("GEOID","rv")

Unique Id is GEOID:

# There are different type of joins, methods to use if id column have different labels but those are things you can learn oyo. left join does exactly what it says it does. 

lower48 <- lower48 |>
  left_join(g_id, by = "GEOID")
#Choropleth GGPlot
map_5 <- ggplot(lower48) +
  geom_sf(aes(fill = rv))