Data Science II

Session 1: Quarto for Spatial Reporting, Coordinates, and CRS

Bogdan G. Popescu
bogdan.popescu@johncabot.edu

John Cabot University

What You’ll Learn Today

  1. Quarto as the container for spatial work — built from scratch, click by click
  2. The four Quarto settings that matter for maps
  3. Why the Earth’s shape forces us to think about coordinate reference systems
  4. Geographic vs. projected systems, and what every projection distorts
  5. Vector layers: points, lines, polygons
  6. Reading and mapping a real shapefile of Mexico

By the end, you will have a Quarto document that renders a map of Mexico.

The Plan: One Document, Five Weeks

Start one .qmd today. Add to it every week:

Week What you add
1 A map of Mexico
2 sf objects, your own state or municipio
3 A choropleth of an indicator
4 A spatial join and distances
5 A siting analysis

By Session 5 that document is your reto deliverable.

Before We Start: Your Folder

Every line of code in these slides runs on your machine, provided your folders look like this:

ciencia_datos_2/
├── data/
│   ├── gadm41_MEX_shp/       <- Mexico boundaries (GADM)
│   ├── denue/                <- DENUE point data
│   └── mex_formats/          <- the same states as .geojson and .gpkg
└── session1/
    └── my_notes.qmd          <- you work here

Code in these slides reads ../data/... — up one level, then into data. Keep that shape and everything copies and pastes without edits.

Part A: Quarto for Spatial Reporting

Quarto: Three Layers

Every .qmd file has exactly three kinds of content:

  1. YAML header — metadata and output settings, between --- fences
  1. Markdown — your prose, headings, tables, links
  1. Code chunks — R that runs when you render

Let’s build one together, from an empty RStudio window.

Step 1: Open RStudio

Nothing open yet — just the Console on the left and an empty Environment on the right.

Step 2: File → New File → Quarto Document

Note Quarto Document and Quarto Presentation are separate entries. Today: Document.

Step 3: The New Document Dialog

What the Dialog Options Mean

Title — goes straight into the YAML header. You can change it later.

HTML / PDF / Word — HTML is the right default. PDF needs a LaTeX install; Word needs Word.

Engine: Knitr — the R engine. Leave it.

Use visual markdown editor — a word-processor-like view.

Tip

Leave the visual editor checked for now. You can toggle between Visual and Source at any time — they edit the same file.

Step 4: The Visual Editor

RStudio gives you a template document. Headings look like headings; code chunks look like boxes.

Step 5: The Source Editor

The same file, shown as raw text. ## Quarto is a heading; ```{r} opens a chunk.

Visual or Source?

Visual

  • Easier for prose
  • Tables and links via toolbar
  • Good when writing the report

Source

  • You see exactly what is in the file
  • Necessary for chunk options
  • Good when debugging

Note

Everything in these slides is written in Source. When something goes wrong, Source is where you find out why.

Step 6: Clear the Template

Delete the sample text. Keep the YAML header, one heading, and one chunk.

Step 7: Save It

Warning

Save before you render. An unsaved document has no folder, so relative paths like data/... have nothing to be relative to.

Step 8: Now It Has a Name

The tab now reads example.qmd, and the tooltip shows the Render shortcut: Shift+Cmd+K (Shift+Ctrl+K on Windows).

Step 9: Write Some Text

Plain prose between the chunks. This is the part that makes it a report and not a script.

Step 10: Comment Your Code

Inside a chunk, # starts a comment — exactly as in any R script.

Step 11: Run the Chunk

The green ▶ runs the chunk. Output appears inline, under the chunk, and in the Console.

Step 12: Insert Another Chunk

The circled button inserts a new chunk — or press Option+Cmd+I (Alt+Ctrl+I on Windows). Learn this shortcut; you will use it constantly.

Step 13: Render

Step 14: The Result

Prose, code, and output — one document, reproducible from source. That is the whole point.

The YAML Header

---
title: "Access to Health Services"
author: "Your Name"
date: "2026-09-15"
format: html
toc: true
toc-title: "Contents"
---

# Introduction

This is an example.
  • format: html — what gets produced
  • toc: true — table of contents
  • toc-depth: 2 — how many heading levels
  • toc-title: — rename it

YAML is whitespace-sensitive. Two spaces of indentation mean something; a tab breaks it.

What toc-depth Does

The same document, rendered three times. Only toc-depth changes.

toc-depth: 1

toc-depth: 2

toc-depth: 3

Each level deeper adds the next rank of headings — #, then ##, then ###. The default is 3.

Code Chunks: The Options That Matter

```{r}
#| echo: true        # show the code
#| eval: true        # run the code
#| warning: false    # hide warnings
#| message: false    # hide package chatter
summary(mtcars)
```

echo and eval are independent. echo: true, eval: false shows code without running it — useful for an expensive download step you do not want repeated on every render.

Chunk Options for Maps

This is the part that is not in a normal Quarto tutorial.

```{r}
#| fig-width: 8       # plot canvas, in inches
#| fig-height: 5      # taller than wide is usually wrong for Mexico
#| out-width: "100%"  # how big it appears in the document
#| fig-dpi: 150       # resolution; 96 for slides, 150+ for print
```

Tip

Mexico is wider than it is tall — about 32 degrees of longitude by 18 of latitude. Start at fig-width: 8, fig-height: 5 and adjust from there.

The Noise Problem

Every package announces itself. Here is dplyr in a rendered report:

Nobody wants to read “The following objects are masked from ‘package:stats’” in a policy report.

Installing What You Need

Run this once, in the Console — not in your document:

install.packages(c("sf", "ggplot2", "dplyr", "tidyverse",
                   "spData", "cowplot", "ggspatial",
                   "classInt", "viridis", "leaflet", "patchwork"))

Warning

sf is the one that can fail. It needs GDAL, GEOS and PROJ underneath. On Mac and Windows the binary from CRAN includes them; on Linux install libgdal-dev, libgeos-dev and libproj-dev first.

Check it worked:

library(sf)
sf_extSoftVersion()[1:3]

Getting the Data

Everything the code in these slides reads — the Mexico shapefiles, DENUE, and the .geojson and .gpkg copies — is in one Dropbox folder:

Download the data/ folder

Unzip it so that data/ sits next to your session1/ folder, exactly as on the folder slide. Then ../data/gadm41_MEX_shp/... resolves and every chunk runs unedited.

Silencing It

sf is worse than dplyr — it reports the geometry type, the bounding box, and the full CRS on every single read.

```{r}
#| warning: false
#| message: false
mex <- st_read(
  "../data/gadm41_MEX_shp/gadm41_MEX_1.shp",
  quiet = TRUE)
```

Important

You need both. quiet = TRUE silences sf’s own reporting; the chunk options silence R’s warnings and messages. They are different channels.

Caching Slow Reads

Municipio-level shapefiles are large. Mexico’s level-2 file is 53 MB.

Without caching, you re-read it on every render, every time you fix a typo.

```{r}
#| cache: true
mex2 <- st_read("../data/gadm41_MEX_shp/gadm41_MEX_2.shp", quiet = TRUE)
```

Warning

cache: true notices when the code changes — not when the data file changes. If you replace the shapefile, delete the _cache folder by hand.

embed-resources: The Trap

You build an interactive map. It works on your laptop. You email your reader the .html file — just that one file, the way you would send any attachment.

They see a blank white box.

Why It Happens

An HTML file by default links to its JavaScript, CSS, and images, which sit in a _files folder next to it. Send the HTML alone and every link breaks.

---
format:
  html:
    embed-resources: true
---

This inlines everything into a single, larger, self-sufficient file.

Important

One rule, two halves: embed the resources so everything lives inside the file, then send that single .html file. Nothing else travels with it — no _files folder, no images. For anything you send to anyone, embed-resources: true is not optional. We build a leaflet map in Session 3 — that is exactly when this bites.

Tables and Callouts

Table

| State   | Municipios |
|:--------|-----------:|
| Oaxaca  |        570 |
| Puebla  |        217 |
| Yucatán |        106 |
State Municipios
Oaxaca 570
Puebla 217
Yucatán 106

Callout

:::{.callout-note}
Source: INEGI, Marco
Geoestadístico 2020.
:::

Note

Source: INEGI, Marco Geoestadístico 2020.

Types: note, tip, warning, caution, important

Your Starting Template

In the course folder, one level up: template_reto.qmd — every setting from this part, already in place.

---
title: "My Project"
format:
  html:
    toc: true
    embed-resources: true
execute:
  warning: false
  message: false
---

execute: at the top level applies to every chunk, so you stop repeating yourself.

Open it now. This is the document you carry to Session 5.

Part B: Coordinates and Coordinate Reference Systems

From Documents to Maps

We have a container. Now: how does a place on Earth become a number we can plot?

This is the part people skip. Skipping it is why maps come out empty, sideways, or refusing to combine.

The Earth Is Not a Sphere

  • The Earth is not a perfect sphere
  • It bulges at the equator — an oblate spheroid
  • Any flat map is therefore an approximation of an approximation

Meridians

Meridians run north–south.

They measure longitude — east or west of Greenwich, from -180 to 180.

The prime meridian is in red. Mexico, in orange, sits west of it.

Parallels

Parallels run east–west.

They measure latitude — north or south of the equator, from -90 to 90.

The equator is in red. Mexico, in orange, sits north of it.

The Geographic Graticule

Together they form a grid. Mexico sits roughly between -118 and -87 longitude, 14 and 33 latitude.

Warning

sf expects longitude first, then latitude — x before y. Reversing them is the single most common way to put Mexico in the Indian Ocean.

What Is a Coordinate Reference System?

To place an object on a map you need two things:

  • the coordinates of the object
  • a system of reference saying how those coordinates relate to a physical location on Earth

A CRS has three components:

  • Sphere and ellipsoid description
  • Geoid
  • Datum

CRS Component 1: Sphere and Ellipsoid

  • We can assume the Earth is a perfect sphere — this simplifies the mathematics
  • For more accurate measurement we use an ellipsoid
  • Distance calculations differ between the two

CRS Component 2: The Geoid

  • The Earth’s surface is not smooth
  • These undulations are invisible to the eye but matter for local measurement
  • Below, the undulations are exaggerated 4,000 times

CRS Component 3: The Datum

To reconcile a simple mathematical model with the undulating shape of the Earth, we:

  • align the geoid with the ellipsoid
  • map the Earth’s surface features onto that ellipsoid

A datum is how you choose to align them.

Local and Geocentric Datums

Local datums fit one region well:

  • NAD27 — continental US
  • ED50 — Western Europe

Geocentric datums fit the whole world:

  • WGS84
  • NAD83

Two Families

There are two kinds of CRS, and the difference decides what you can do:

Geographic Projected
Units degrees metres
Shape curved surface flat plane
Example WGS84 (EPSG:4326) Mexico LCC (EPSG:6372)
Good for storing, sharing measuring, mapping

Going from geographic to projected requires a mathematical transformation — and every such transformation distorts something.

WGS84, or EPSG:4326

The most common geographic system is the World Geodetic System of 1984.

It was established because navigation, aviation, and geography needed global maps.

Almost everything you download — GADM, INEGI, OpenStreetMap — arrives in WGS84.

EPSG:4326 is its code. EPSG codes let you name a CRS with one short integer, and you will use them constantly.

Every projected system has one too — and choosing between them is the next slide.

Finding an EPSG Code

Search any country at epsg.io — the screenshot searches Italy; for Mexico the answer is EPSG:6372.

Projected Coordinate Systems

  • A projected coordinate system is a reference system for identifying locations and measuring features on a flat (map) surface

  • Going from a geographic system to a projected one requires mathematical transformations

  • To perform spatial analyses correctly, all your data has to be in the same coordinate system

There is no single correct projection. Each one preserves something and sacrifices something else.

Types of Projected Coordinate Systems

#Step1: Load the libraries
library(sf); library(tidyverse); library(cowplot); library(ggspatial)
sf_use_s2(FALSE)

#Step2: Get world polygons from the spData package
dunia <- spData::world %>% st_as_sf()

#Step3: Draw them in one projection
ggplot() +
  geom_sf(data = dunia) +
  coord_sf(crs = "EPSG: 3857") +
  ggtitle("World Mercator") +
  theme_minimal_grid()

Types of Projected Coordinate Systems

Types of Projected Coordinate Systems

#Only the CRS code and the title change
ggplot() +
  geom_sf(data = dunia) +
  coord_sf(crs = "ESRI:54030") +
  ggtitle("World Robinson") +
  theme_minimal_grid()

Types of Projected Coordinate Systems

Types of Projected Coordinate Systems

#Only the CRS code and the title change
ggplot() +
  ggspatial::layer_spatial(data = dunia) +
  coord_sf(crs = "ESRI:54032") +
  ggtitle("World Azimuthal Equidistant") +
  theme_minimal_grid()

Types of Projected Coordinate Systems

Types of Projected Coordinate Systems

#Only the CRS code and the title change
ggplot() +
  ggspatial::layer_spatial(data = dunia) +
  coord_sf(crs = "ESRI:54002") +
  ggtitle("World Equidistant Cylindrical") +
  theme_minimal_grid()

Types of Projected Coordinate Systems

Types of Projected Coordinate Systems

#Only the CRS code and the title change
ggplot() +
  ggspatial::layer_spatial(data = dunia) +
  coord_sf(crs = "ESRI:54010") +
  ggtitle("World Eckert VI") +
  theme_minimal_grid()

Types of Projected Coordinate Systems

Types of Projected Coordinate Systems

#Only the CRS code and the title change
ggplot() +
  ggspatial::layer_spatial(data = dunia) +
  coord_sf(crs = "ESRI:54031") +
  ggtitle("World Two Point Equidistant") +
  theme_minimal_grid()

Types of Projected Coordinate Systems

So Which One for Mexico?

Mercator is the one you have seen most — it is what web maps use — and it is the worst choice for a national map. It inflates everything away from the equator.

For Mexico, use EPSG:6372 (Mexico ITRF2008 / LCC):

  • a Lambert Conformal Conic, fitted to Mexico’s latitudes
  • units in metres, so distances and areas are measurable
  • the standard for official Mexican cartography
ggplot() + geom_sf(data = mex) + coord_sf(crs = 6372)

Tip

coord_sf(crs = ...) reprojects for display only. st_transform() changes the data itself. Use st_transform when you are going to measure something.

Why This Matters: The Error You Will Hit

Modern sf is good at geodesic mathematics. Ask for a distance between two points stored in degrees and it still returns a sensible answer in metres.

So CRS mistakes rarely corrupt your numbers quietly. They stop you loudly, the moment you combine two layers:

#Step1: Make a point for Mexico City, in WGS84 (degrees)
cdmx <- st_as_sf(data.frame(lon = -99.1332, lat = 19.4326),
                 coords = c("lon", "lat"), crs = 4326)

#Step2: Project the states to EPSG:6372 (metres) -- a DIFFERENT crs
mex_projected <- st_transform(mex1, 6372)

#Step3: Try to combine them
st_join(cdmx, mex_projected)
Error in `st_geos_binop()`:
! st_crs(x) == st_crs(y) is not TRUE

You will meet this exact error in Session 4. Now you know what it means.

st_transform Is the Fix

Put both layers in the same CRS. That is the whole solution.

st_crs(cdmx)$input
[1] "EPSG:4326"
st_crs(mex_projected)$input
[1] "EPSG:6372"
#Step1: Put the point into the SAME crs as the polygons, then join
result <- st_join(st_transform(cdmx, 6372), mex_projected)

#Step2: Which state did the point land in?
result$NAME_1
[1] "Distrito Federal"

Tip

Habit worth forming: transform everything to one CRS immediately after reading it, before you do anything else.

A Warning About Boundary Data

Look at that answer again: “Distrito Federal”.

GADM still uses the name retired in 2016, when the Distrito Federal became Ciudad de México.

Boundary files carry their own vintage and their own politics:

  • names change
  • municipios split and merge
  • state borders are occasionally disputed

Important

Always check when your boundary file was made, and whether its keys match the dataset you plan to join to it. This is one more reason to prefer INEGI for official Mexican work.

Vector and Raster

  • Spatial data comes as a vector or a raster
  • A river as a vector (left) and as a raster (right)

In these five sessions we work with vectors. Rasters are a course of their own.

Vector Layers

A vector layer is a set of geometries attached to a table of non-spatial attributes.

Geometries are sequences of coordinates forming points, lines, or polygons.

Points

A point is a single coordinate pair.

A hospital, a school, a DENUE establishment.

Lines

A line is an ordered sequence of points.

A road, a river, a bus route.

Polygons

A polygon is a closed line enclosing an area.

A municipio, a state, a colonia.

Part C: Your First Map of Mexico

Where to Get Shapefiles: GADM

Country and subnational boundaries for the whole world: https://gadm.org

Choosing a Country

Choosing a Country

Choosing a Country

Choosing a Country

GADM gives you nested levels: _0 country, _1 states, _2 municipios. For Mexico that is 1, 32, and 2,466 features.

What You Download

Warning

A “shapefile” is not one file. It is at least .shp (geometry), .dbf (attributes), .shx (index), and .prj (the CRS). Move them together or the layer breaks — and a missing .prj means R has no idea where on Earth your data is.

Prefer INEGI for Official Work

GADM is convenient and global. But for Mexican policy work, INEGI’s Marco Geoestadístico is better:

  • it carries CVE_ENT and CVE_MUN, the keys every other INEGI dataset joins on
  • the boundaries are the official ones
  • it goes down to localidad and AGEB

https://www.inegi.org.mx/temas/mg/

We use GADM today because it is one click. Use INEGI for your reto.

An sf Object Is a Dataframe

#Step1: Load sf
library(sf)

#Step2: Read the state boundaries
mex_states <- st_read("../data/gadm41_MEX_shp/gadm41_MEX_1.shp", quiet = TRUE)

#Step3: Check what we got back
class(mex_states)
[1] "sf"         "data.frame"
mex_states[1:3, c("NAME_1", "TYPE_1")]
Simple feature collection with 3 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -118.3665 ymin: 21.62227 xmax: -101.8353 ymax: 32.71863
Geodetic CRS:  WGS 84
               NAME_1 TYPE_1                       geometry
1      Aguascalientes Estado MULTIPOLYGON (((-102.0659 2...
2     Baja California Estado MULTIPOLYGON (((-114.1228 2...
3 Baja California Sur Estado MULTIPOLYGON (((-109.9104 2...

Why That Matters

This is the most important idea in the next four weeks:

Everything you know from dplyr still works.

filter, mutate, group_by, summarise, left_join — all of it. There is simply an extra geometry column that travels along with every row.

Tip

st_drop_geometry() gives you the plain table back when you want to inspect the attributes without pages of coordinates.

st_simplify: Making Maps Render Fast

Detailed boundaries are slow to draw, and at national scale you cannot see the detail anyway.

#Step1: Simplify the boundaries (0.02 degrees is roughly 2 km)
mex_simple <- st_simplify(mex_states, dTolerance = 0.02)

#Step2: Compare the sizes before and after
format(object.size(mex_states), units = "MB")
[1] "11.9 Mb"
format(object.size(mex_simple), units = "MB")
[1] "0.1 Mb"

11.9 MB to 0.1 MB — a hundredfold reduction, with no visible difference at this scale.

Warning

dTolerance is in the units of your CRS. In WGS84 that is degrees, so 0.02 is roughly 2 km. sf warns you about this — the warning is correct, and simplifying in a projected CRS is more defensible.

st_simplify: Before and After

The middle map is what the code above produces — invisible at this scale. Turn the dial further and the coastline collapses into straight lines.

Mapping Mexico

#Step1: Load ggplot2
library(ggplot2)

#Step2: Read the country outline
mex_country <- st_read(
  "../data/gadm41_MEX_shp/gadm41_MEX_0.shp",
  quiet = TRUE)

#Step3: Map it, projected to EPSG:6372
ggplot() +
  geom_sf(data = mex_country,
          fill = "#dbe6f0",
          colour = "#33475b",
          linewidth = 0.3) +
  coord_sf(crs = 6372) +
  theme_bw()

Adding the States

ggplot() +
  geom_sf(data = mex_simple,
          fill  = "#dbe6f0",
          colour = "#5a6b7d",
          linewidth = 0.2) +
  coord_sf(crs = 6372) +
  labs(title = "States of Mexico",
       caption = "Source: GADM 4.1") +
  theme_bw()

Two layers, same ggplot. geom_sf stacks like any other geom.

Zooming In

#Step1: Keep only the three peninsula states
peninsula <- mex_simple[
  mex_simple$NAME_1 %in%
    c("Yucatán", "Campeche",
      "Quintana Roo"), ]

#Step2: Map them, coloured by state
ggplot() +
  geom_sf(data = peninsula,
          aes(fill = NAME_1),
          colour = "white") +
  scale_fill_brewer(palette = "Blues") +
  labs(fill = "State") +
  theme_bw()

filter the sf object like a dataframe, and the map follows. No special spatial subsetting needed.

Putting It in Your Document

Everything from this session, as one chunk in template_reto.qmd:

```{r}
#| label: fig-mexico
#| fig-cap: "States of Mexico"
#| fig-width: 8
#| fig-height: 5
#| cache: true
library(sf); library(ggplot2)

mex <- st_read("../data/gadm41_MEX_shp/gadm41_MEX_1.shp", quiet = TRUE)
mex <- st_transform(mex, 6372)          # project immediately
mex <- st_simplify(mex, dTolerance = 2000)   # metres now, not degrees

ggplot() + geom_sf(data = mex, fill = "#dbe6f0") + theme_bw()
```

Note dTolerance = 2000: once you have transformed to EPSG:6372 the units are metres, so the tolerance changes meaning. Same idea, different number.

What Have We Learned?

  • Built a Quarto document from an empty RStudio window, and rendered it
  • The settings that matter for maps: fig-width, cache, execute: warning/message, and embed-resources
  • The Earth is an oblate spheroid, located by latitude and longitude
  • A CRS combines an ellipsoid, a geoid, and a datum
  • Geographic systems use degrees, projected systems use metres — every projection distorts something
  • CRS mistakes announce themselves when you combine layers; st_transform is the fix
  • An sf object is a dataframe with a geometry column

Next session: sf geometries in depth — points, lines, polygons, and building maps layer by layer.