Week 3: Vector and Raster Representations

PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2026

Jeff Jacobs

jj1088@georgetown.edu

Wednesday, September 9, 2026

Logistical Table-Setting

  • HW1 Public Tests will be updated ASAP
  • TA Intros!
  • Coding Workshops!

TA Intros

(In alphabetical order by surname wohoo)

Jackson Howes (DSAN)
jh2787@georgetown.edu

  • BS in Government AND Mathematics (Honors), St. Lawrence, Canton, NY
  • Survived fire and flames of this class 🤯

Hermine Wilhelmsen (McCourt)
hw623@georgetown.edu

  • Recent MSc in Sustainability and Social Innovation, HEC Paris 🥳🧑‍🎓
  • Ombuds-person, UNFCCC alum

Relevant Word of the Day: 山 (Shān) = Mountain

Where We Left Off: Let’s Make Some Dang Maps!

Our First Map: Polygons!

From OpenData DC… Last week .shp, this week .geojson!

Code 
library(tidyverse)
library(rmarkdown)
library(sf)
# Load DC tracts data
dc_sf_fpath <- "data/Census_Tracts_-_2020.geojson"
dc_sf <- st_read(dc_sf_fpath);
Reading layer `Census_Tracts_-_2020' from data source 
  `/Users/jpj/gtown-local/ppol6805/w03/data/Census_Tracts_-_2020.geojson' 
  using driver `GeoJSON'
Simple feature collection with 206 features and 315 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.11976 ymin: 38.79165 xmax: -76.9094 ymax: 38.99585
Geodetic CRS:  WGS 84
Code 
cols_to_keep <- c(
  "OBJECTID", "TRACT", "GEOID", "ALAND", "AWATER",
  "STUSAB", "SUMLEV", "GEOCODE", "STATE", "NAME",
  "POP100", "HU100", "geometry"
)
dc_sf <- dc_sf |> select(cols_to_keep)

sf Objects

dc_sf is an R object of type sf (short for “simple features”), which extends data.frame by adding a special column named geometry (containing POLYGONs)

Code 
class(dc_sf)
[1] "sf"         "data.frame"
Code 
dim(dc_sf)
[1] 206  13
Code 
rmarkdown::paged_table(dc_sf, options=list(rows.print=6))

What Are These POLYGON Objects?

Code 
head(dc_sf)
Simple feature collection with 6 features and 12 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.05714 ymin: 38.90524 xmax: -77.00865 ymax: 38.96135
Geodetic CRS:  WGS 84
  OBJECTID  TRACT       GEOID  ALAND AWATER STUSAB SUMLEV     GEOCODE STATE
1        1 002002 11001002002 849376      0     DC    140 11001002002    11
2        2 002101 11001002101 600992      0     DC    140 11001002101    11
3        3 002102 11001002102 725975      0     DC    140 11001002102    11
4        4 002201 11001002201 415173      0     DC    140 11001002201    11
5        5 002202 11001002202 698895    566     DC    140 11001002202    11
6        6 000101 11001000101 199776   5261     DC    140 11001000101    11
                NAME POP100 HU100                       geometry
1 Census Tract 20.02   4072  1532 POLYGON ((-77.03642 38.9497...
2 Census Tract 21.01   5687  2335 POLYGON ((-77.02824 38.9581...
3 Census Tract 21.02   5099  2221 POLYGON ((-77.01997 38.9581...
4 Census Tract 22.01   3485  1229 POLYGON ((-77.02743 38.9518...
5 Census Tract 22.02   3339  1454 POLYGON ((-77.01968 38.9520...
6  Census Tract 1.01   1406   999 POLYGON ((-77.05714 38.9105...

Zoom… Enhance…

Let’s try printing just the first row?

Code 
rmarkdown::paged_table(dc_sf[1,])

Printing the special geometry column?

Code 
rmarkdown::paged_table(dc_sf[1, "geometry"])

Nope, we’ll have to open the “black box” using the st_coordinates() function from the sf library →

Code 
sf::st_coordinates(dc_sf[1,"geometry"])
               X        Y L1 L2
  [1,] -77.03642 38.94975  1  1
  [2,] -77.03643 38.94827  1  1
  [3,] -77.03628 38.94827  1  1
  [4,] -77.03539 38.94827  1  1
  [5,] -77.03452 38.94827  1  1
  [6,] -77.03442 38.94827  1  1
  [7,] -77.03375 38.94827  1  1
  [8,] -77.03309 38.94827  1  1
  [9,] -77.03287 38.94827  1  1
 [10,] -77.03273 38.94827  1  1
 [11,] -77.03170 38.94824  1  1
 [12,] -77.03163 38.94819  1  1
 [13,] -77.03146 38.94802  1  1
 [14,] -77.03112 38.94769  1  1
 [15,] -77.03109 38.94767  1  1
 [16,] -77.03107 38.94765  1  1
 [17,] -77.03104 38.94762  1  1
 [18,] -77.03098 38.94753  1  1
 [19,] -77.03080 38.94735  1  1
 [20,] -77.03071 38.94726  1  1
 [21,] -77.03069 38.94724  1  1
 [22,] -77.03044 38.94699  1  1
 [23,] -77.03023 38.94679  1  1
 [24,] -77.03005 38.94661  1  1
 [25,] -77.02985 38.94642  1  1
 [26,] -77.02982 38.94640  1  1
 [27,] -77.02965 38.94613  1  1
 [28,] -77.02893 38.94614  1  1
 [29,] -77.02821 38.94613  1  1
 [30,] -77.02780 38.94613  1  1
 [31,] -77.02747 38.94613  1  1
 [32,] -77.02702 38.94613  1  1
 [33,] -77.02660 38.94613  1  1
 [34,] -77.02667 38.94641  1  1
 [35,] -77.02679 38.94720  1  1
 [36,] -77.02681 38.94738  1  1
 [37,] -77.02694 38.94826  1  1
 [38,] -77.02695 38.94834  1  1
 [39,] -77.02695 38.94840  1  1
 [40,] -77.02696 38.94847  1  1
 [41,] -77.02701 38.94877  1  1
 [42,] -77.02702 38.94884  1  1
 [43,] -77.02704 38.94900  1  1
 [44,] -77.02705 38.94909  1  1
 [45,] -77.02708 38.94934  1  1
 [46,] -77.02714 38.94974  1  1
 [47,] -77.02716 38.94993  1  1
 [48,] -77.02721 38.95028  1  1
 [49,] -77.02723 38.95037  1  1
 [50,] -77.02724 38.95046  1  1
 [51,] -77.02725 38.95055  1  1
 [52,] -77.02726 38.95063  1  1
 [53,] -77.02728 38.95072  1  1
 [54,] -77.02729 38.95081  1  1
 [55,] -77.02740 38.95163  1  1
 [56,] -77.02743 38.95188  1  1
 [57,] -77.02756 38.95281  1  1
 [58,] -77.02756 38.95283  1  1
 [59,] -77.02757 38.95291  1  1
 [60,] -77.02773 38.95400  1  1
 [61,] -77.02774 38.95409  1  1
 [62,] -77.02779 38.95451  1  1
 [63,] -77.02783 38.95476  1  1
 [64,] -77.02787 38.95505  1  1
 [65,] -77.02788 38.95516  1  1
 [66,] -77.02789 38.95523  1  1
 [67,] -77.02796 38.95572  1  1
 [68,] -77.02800 38.95602  1  1
 [69,] -77.02802 38.95612  1  1
 [70,] -77.02803 38.95622  1  1
 [71,] -77.02807 38.95653  1  1
 [72,] -77.02812 38.95687  1  1
 [73,] -77.02814 38.95703  1  1
 [74,] -77.02816 38.95718  1  1
 [75,] -77.02818 38.95729  1  1
 [76,] -77.02821 38.95750  1  1
 [77,] -77.02821 38.95758  1  1
 [78,] -77.02822 38.95765  1  1
 [79,] -77.02823 38.95781  1  1
 [80,] -77.02824 38.95797  1  1
 [81,] -77.02824 38.95805  1  1
 [82,] -77.02824 38.95813  1  1
 [83,] -77.02824 38.95829  1  1
 [84,] -77.02824 38.95836  1  1
 [85,] -77.02823 38.95851  1  1
 [86,] -77.02822 38.95872  1  1
 [87,] -77.02817 38.95935  1  1
 [88,] -77.02811 38.96003  1  1
 [89,] -77.02806 38.96074  1  1
 [90,] -77.02803 38.96102  1  1
 [91,] -77.02803 38.96105  1  1
 [92,] -77.02802 38.96120  1  1
 [93,] -77.02802 38.96123  1  1
 [94,] -77.02806 38.96127  1  1
 [95,] -77.02816 38.96130  1  1
 [96,] -77.02852 38.96135  1  1
 [97,] -77.02851 38.96130  1  1
 [98,] -77.02851 38.96128  1  1
 [99,] -77.02851 38.96126  1  1
[100,] -77.02854 38.96121  1  1
[101,] -77.02855 38.96119  1  1
[102,] -77.02859 38.96115  1  1
[103,] -77.02876 38.96100  1  1
[104,] -77.02925 38.96054  1  1
[105,] -77.02948 38.96033  1  1
[106,] -77.02952 38.96030  1  1
[107,] -77.02953 38.96028  1  1
[108,] -77.02958 38.96026  1  1
[109,] -77.02962 38.96024  1  1
[110,] -77.02965 38.96023  1  1
[111,] -77.02967 38.96022  1  1
[112,] -77.02970 38.96022  1  1
[113,] -77.02999 38.96003  1  1
[114,] -77.03003 38.96000  1  1
[115,] -77.03011 38.95994  1  1
[116,] -77.03013 38.95992  1  1
[117,] -77.03016 38.95989  1  1
[118,] -77.03038 38.95966  1  1
[119,] -77.03061 38.95938  1  1
[120,] -77.03073 38.95925  1  1
[121,] -77.03081 38.95915  1  1
[122,] -77.03115 38.95874  1  1
[123,] -77.03138 38.95847  1  1
[124,] -77.03148 38.95836  1  1
[125,] -77.03158 38.95824  1  1
[126,] -77.03172 38.95807  1  1
[127,] -77.03181 38.95796  1  1
[128,] -77.03204 38.95770  1  1
[129,] -77.03224 38.95745  1  1
[130,] -77.03231 38.95737  1  1
[131,] -77.03238 38.95729  1  1
[132,] -77.03252 38.95712  1  1
[133,] -77.03257 38.95706  1  1
[134,] -77.03266 38.95695  1  1
[135,] -77.03271 38.95689  1  1
[136,] -77.03294 38.95662  1  1
[137,] -77.03298 38.95658  1  1
[138,] -77.03318 38.95637  1  1
[139,] -77.03341 38.95620  1  1
[140,] -77.03376 38.95566  1  1
[141,] -77.03379 38.95563  1  1
[142,] -77.03420 38.95514  1  1
[143,] -77.03460 38.95466  1  1
[144,] -77.03504 38.95413  1  1
[145,] -77.03520 38.95394  1  1
[146,] -77.03560 38.95347  1  1
[147,] -77.03588 38.95314  1  1
[148,] -77.03604 38.95295  1  1
[149,] -77.03612 38.95285  1  1
[150,] -77.03625 38.95271  1  1
[151,] -77.03642 38.95254  1  1
[152,] -77.03642 38.95188  1  1
[153,] -77.03642 38.95081  1  1
[154,] -77.03642 38.95053  1  1
[155,] -77.03642 38.94975  1  1

Working With sf Objects

Code 
# Select column by name
head(dc_sf$NAME)
[1] "Census Tract 20.02" "Census Tract 21.01" "Census Tract 21.02"
[4] "Census Tract 22.01" "Census Tract 22.02" "Census Tract 1.01"
Code 
# Select column by index
head(dc_sf[,4])
Simple feature collection with 6 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.05714 ymin: 38.90524 xmax: -77.00865 ymax: 38.96135
Geodetic CRS:  WGS 84
   ALAND                       geometry
1 849376 POLYGON ((-77.03642 38.9497...
2 600992 POLYGON ((-77.02824 38.9581...
3 725975 POLYGON ((-77.01997 38.9581...
4 415173 POLYGON ((-77.02743 38.9518...
5 698895 POLYGON ((-77.01968 38.9520...
6 199776 POLYGON ((-77.05714 38.9105...

And… Actually Displaying the Map!

Code 
# We can extract *just* the geometry column with the st_geometry function...
dc_geo <- sf::st_geometry(dc_sf)
# This produces an object of type "sfc"
class(dc_geo)
[1] "sfc_POLYGON" "sfc"
Code 
# Plot the geometry with base R's plot() function
plot(dc_geo)

More “Modern” Plots with ggplot2!

Code 
library(ggplot2)
dc_sf |>
  ggplot2::ggplot() +
  ggplot2::geom_sf() +
  ggplot2::theme_classic()

Code 
dc_sf |>
  ggplot2::ggplot() +
  ggplot2::geom_sf() +
  ggplot2::theme_classic()

How to Do Things with Geometries

From the sf Cheatsheet

HW1 \(\rightarrow\) HW2

  • Once you finish HW1, you’ll know how to create geometries with sf and terra
  • So now, what can you do with them?
  • For example, we’d like to be able to say things like:
    • “The new lamppost cannot be placed at \((x, y)\), since there is already a building there!”
    • “There are \(N_1\) lampposts in County 1, and \(N_2\) lampposts in County 2”
    • “The average resident in Neighborhood A lives 2 km away from their nearest bus stop

First Things First: Loading and Saving

  • Note how there were no data files in HW1 😱
  • From HW2 onwards (and in your GIS life), we’ll:
    • Download from e.g. city Open Data Portals: geo data files, but also loading on-the-fly (this week)
    • Summarize/aggregate (this week and next week)
    • Visualize findings (“Mapping Libraries” unit)

Vector Formats

Shapefiles (.shp et al.)

A shape“file” is actually (at least) three separate files bundled together:

  • Mandatory .shp: Containing feature geometries
  • Mandatory .shx: Positional indices
  • Mandatory .dbf: Data attributes
  • Optional .prj: Coordinate reference system
  • Optional .xml: Metadata

Shapefiles

Let’s see what’s inside the shapefile we first saw in Week 1, containing data on DC’s Census Tracts: Census Tracts in 2020

DC Census Tracts (with the Georgetown campus tract highlighted!) from OpenData.DC.gov

Shapefile Anatomy

From Rodrigue (2016)

GeoJSON / TopoJSON (.geojson)

  • JavaScript Object Notation: General cross-platform format
  • Useful when data is too complex for e.g. .csv
  • TopoJSON = Memory-efficient GeoJSON
  • Bonus: Inline preview on GitHub!
my_data.geojson
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [30, 20], [45, 40],
            [10, 40], [30, 20]
          ]
        ]
      },
      "properties": {
        "color": "green",
        "area": 3565747
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [15, 5], [40, 10],
            [10, 20], [5, 10], 
            [15, 5]
          ]
        ]
      },
      "properties": {
        "color": "red",
        "area": 3272386
      }
    }
  ]
}

GeoPackage (.gpkg)

Raster Formats

  • GeoTIFF (.tif or .tiff)
    • Based on TIFF format developed at NASA
  • NetCDF (.nc4)
    • Used in earth sciences, as format for data sources measured and distributed multiple times per day over large full-country or full-continent areas.

Coordinate Reference Systems (CRS)

  • EPSG (European Petroleum Survey Group) Registry: Most common way to specify a CRS

    • For example, 4326 is the EPSG code for the WGS84 coordinate system

  • PROJ: Rather than opaque numeric code like EPSG, uses plaintext “proj-strings” containing parameter info: datum, ellipsoid, projection, and units (e.g. meters). Example: PROJ4 code EPSG:4326 is represented as

    +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

  • WKT: Lengthy but human-readable descriptions

Geospatial Operations 1: Unary Operations

Getting the Geometries

Using rnaturalearth with mapview

Code 
set.seed(6805)
library(tidyverse)
library(sf)
library(rnaturalearth)
library(mapview)
france_sf <- ne_countries(country = "France", scale = 50)
(france_map <- mapview(france_sf, label = "geounit", legend = FALSE))

Centroid of France

Code 
france_cent_sf <- sf::st_centroid(france_sf)
france_map + mapview(france_cent_sf, label = "Centroid", legend = FALSE)

One We Already Saw: Union

Computing the union of all geometries in the sf via sf::st_union()

Code 
library(leaflet.extras2)
africa_sf <- ne_countries(continent = "Africa", scale = 50)
africa_union_sf <- sf::st_union(africa_sf)
africa_map <- mapview(africa_sf, label="geounit", legend=FALSE)
africa_union_map <- mapview(africa_union_sf, label="st_union(africa)", legend=FALSE)
africa_map | africa_union_map

Helpful for Rasterizing: BBox

Code 
africa_bbox_sf <- sf::st_bbox(africa_sf)
africa_bbox_map <- mapview(africa_bbox_sf, label="st_bbox(africa)", legend=FALSE)
africa_map | africa_bbox_map

Convex Hulls by Country

Code 
africa_countries_cvx <- sf::st_convex_hull(africa_sf)
africa_countries_cvx_map <- mapview(africa_countries_cvx, label="geounit", legend=FALSE)
africa_map | africa_countries_cvx_map

Convex Hull of Continent

Use st_union() first:

Code 
africa_cvx <- africa_sf |> st_union() |> st_convex_hull()
africa_cvx_map <- mapview(africa_cvx, label="geounit", legend=FALSE)
africa_map | africa_cvx_map

One We Already Saw: Centroids

Computing the centroid of all geometries in the sf via sf::st_centroid()

Code 
africa_cents_sf <- sf::st_centroid(africa_sf)
africa_cents_map <- mapview(africa_cents_sf, label="geounit", legend=FALSE)
africa_map | africa_cents_map

Geospatial Operations 2: Binary Operations

Spatial Joins

Code 
nc <- system.file("shape/nc.shp", package="sf") |>
  read_sf() |>
  st_transform('EPSG:2264')
gr <- st_sf(
         label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = ""),
         geom = st_make_grid(nc))
gr$col <- sf.colors(10, categorical = TRUE, alpha = .3)
# cut, to verify that NA's work out:
gr <- gr[-(1:30),]
suppressWarnings(nc_j <- st_join(nc, gr, largest = TRUE))
par(mfrow = c(2,1), mar = rep(0,4))
plot(st_geometry(nc_j), border = 'grey')
plot(st_geometry(gr), add = TRUE, col = gr$col)
text(st_coordinates(st_centroid(st_geometry(gr))), labels = gr$label, cex = .85)
# the joined dataset:
plot(st_geometry(nc_j), border = 'grey', col = nc_j$col)
text(st_coordinates(st_centroid(st_geometry(nc_j))), labels = nc_j$label, cex = .7)
plot(st_geometry(gr), border = '#88ff88aa', add = TRUE)

Spatial Sampling

Code 
# Sample random points
africa_points_list <- sf::st_sample(africa_union_sf, 10)
africa_points_sf <- sf::st_sf(africa_points_list)
africa_points_map <- mapview(africa_points_sf, label="Random Point", col.regions=cb_palette[1], legend=FALSE)
africa_map + africa_points_map

The “Default” Predicate: st_intersects

Code 
countries_w_points <- africa_sf[africa_points_sf,]
mapview(countries_w_points, label="geounit", legend=FALSE) + africa_points_map

Counting with lengths()

Code 
country_inter <- sf::st_intersects(africa_sf, africa_points_sf)
# Computes point counts for each polygon
(num_intersections <- lengths(country_inter))
 [1] 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0
[39] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2
Code 
africa_sf <- africa_sf |> mutate(
  num_points = num_intersections
) |> arrange(geounit)
africa_sf |> select(geounit, num_points) |> head()
Simple feature collection with 6 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -8.68335 ymin: -26.8542 xmax: 30.81143 ymax: 37.09238
Geodetic CRS:  WGS 84
       geounit num_points                       geometry
1      Algeria          2 MULTIPOLYGON (((8.576563 36...
2       Angola          2 MULTIPOLYGON (((13.07275 -4...
3        Benin          0 MULTIPOLYGON (((1.622656 6....
4     Botswana          0 MULTIPOLYGON (((25.25879 -1...
5 Burkina Faso          0 MULTIPOLYGON (((0.9004883 1...
6      Burundi          0 MULTIPOLYGON (((30.55361 -2...

Plotting with mapview

Code 
mapview(africa_sf, zcol="num_points")

Plotting with ggplot2

Since we’re starting to get into data attributes rather than geometric features, switching to ggplot2 is recommended!

Code 
africa_sf |> ggplot(aes(fill=num_points)) +
  geom_sf() +
  theme_classic()

Getting Fancier…

  • To do fancier geospatial operations, we’ll need to start overthinking the different possible relationships between two or more geometries!
  • To this end: predicates

References

Rodrigue, Jean-Paul. 2016. The Geography of Transport Systems. Taylor & Francis.