Data Science II

Session 2: Coordinates, CRS, and Your First Map

Bogdan G. Popescu
bgpopescu@tec.mx

Tec de Monterrey

Where We Left Off

Last week: Quarto — chunks, options, and sending a single self-contained .html.

We drew one map of Mexico and moved on. Today we ask where it actually is.

What You’ll Learn Today

  1. Why the Earth’s shape forces us to think about coordinate reference systems
  2. Geographic vs projected systems, and what every projection distorts
  3. The CRS error you will hit, and st_transform
  4. What an sf dataframe is, and how geom_sf() reads its geometry column
  5. Vector layers: points, lines, polygons, and the seven sf geometry types
  6. Reading a real shapefile of Mexico, and what is inside one

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.

Tip

A working rule: leave your data in WGS84 — it is what you downloaded, and it is what you hand on to someone else. Reproject only when you measure — distances, areas, buffers — because those need metres, not degrees.

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.

Or Just Stick to WGS84

Both layers have to use the same CRS, but it does not have to be 6372.

#Step1: Bring the polygons back to the CRS the files arrived in
mex_wgs <- st_transform(mex_projected, 4326)

#Step2: The same join, with nothing done to the point
st_join(cdmx, mex_wgs)$NAME_1
[1] "Distrito Federal"

Same answer. For a point-in-polygon question, WGS84 is enough.

Tip

Keep everything in WGS84 as you read it. Project to 6372 only when you need to measure something: a distance, an area, a buffer.

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.

What an sf Dataframe Actually Is

Start with an ordinary dataframe. Nothing spatial about it — two of the columns just happen to hold numbers we call coordinates.

#Step1: A plain dataframe
clinics <- data.frame(
  city = c("Monterrey", "Guadalajara", "Mérida"),
  beds = c(420, 310, 180),
  lon  = c(-100.3161, -103.3496, -89.5926),
  lat  = c(  25.6866,   20.6597,  20.9674))

#Step2: Promote it to sf -- lon FIRST, then lat, and say which CRS they are in
clinics_sf <- st_as_sf(clinics, coords = c("lon", "lat"), crs = 4326)

clinics_sf
Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -103.3496 ymin: 20.6597 xmax: -89.5926 ymax: 25.6866
Geodetic CRS:  WGS 84
         city beds                  geometry
1   Monterrey  420 POINT (-100.3161 25.6866)
2 Guadalajara  310 POINT (-103.3496 20.6597)
3      Mérida  180  POINT (-89.5926 20.9674)

lon and lat are gone. A single geometry column took their place, and the CRS is now attached to the object.

Still a Dataframe, With One Extra Column

#Step1: It is BOTH an sf object and a data.frame
class(clinics_sf)
[1] "sf"         "data.frame"
#Step2: The columns you started with, plus one
names(clinics_sf)
[1] "city"     "beds"     "geometry"
#Step3: sf remembers which column holds the geometry
attr(clinics_sf, "sf_column")
[1] "geometry"
#Step4: And that column is a list of POINTs
st_geometry(clinics_sf)
Geometry set for 3 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -103.3496 ymin: 20.6597 xmax: -89.5926 ymax: 25.6866
Geodetic CRS:  WGS 84

Tip

sf tracks the geometry column by name. That is how every st_ function knows which column to operate on, and why you never pass the coordinates around yourself.

How ggplot Finds the Geometry

#Note: no aes(x = , y = ) anywhere
ggplot() +
  geom_sf(data = mex1, fill = "#dbe6f0") +
  geom_sf(data = clinics_sf,
          colour = "#e8883a", size = 4) +
  theme_bw()

geom_sf() looks up the column named by sf_column and draws that. You never tell it where x and y are — it already knows.

Without sf you would be writing geom_point(data = clinics, aes(x = lon, y = lat)), and you would have to handle the projection yourself.

The Geometry Column Decides the Mark

The same geom_sf() call draws a dot, a line, or a filled shape. What changes is what is sitting in the geometry column.

#Step1: Three objects, three geometry types, identical attributes
p <- st_sf(kind = "POINT",      geometry = st_sfc(st_point(c(5, 4))))
l <- st_sf(kind = "LINESTRING", geometry = st_sfc(
             st_linestring(rbind(c(1,2), c(5,6), c(9,3)))))
a <- st_sf(kind = "POLYGON",    geometry = st_sfc(
             st_polygon(list(rbind(c(2,2), c(8,2), c(8,6), c(2,2))))))

#Step2: One geom_sf() call each -- nothing else differs
patchwork::wrap_plots(
  ggplot() + geom_sf(data = p, size = 5) + ggtitle("POINT")      + theme_void(),
  ggplot() + geom_sf(data = l, linewidth = 1.2) + ggtitle("LINESTRING") + theme_void(),
  ggplot() + geom_sf(data = a, fill = "#dbe6f0") + ggtitle("POLYGON") + theme_void())

This is why the next slides matter: the geometry type is what your map will look like, before you choose a single colour.

The Geometry Column Decides the Mark

Attributes Colour It, Geometry Places It

aes() still maps the ordinary columns. The geometry column decides where; the attributes decide how it looks.

ggplot() +
  geom_sf(data = mex1, fill = "#f2f5f8") +
  geom_sf(data = clinics_sf,
          aes(size = beds, colour = city)) +
  scale_size(range = c(3, 8),
             breaks = clinics_sf$beds) +
  theme_bw()

beds and city are plain columns — exactly the aes() you already know from ggplot2. Nothing spatial about that half of the grammar.

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.

The Seven Simple Feature Types

sf builds every geometry from one of seven constructors:

one piece several pieces
st_point st_multipoint
st_linestring st_multilinestring
st_polygon st_multipolygon

st_geometrycollection holds a mixture of the other six — single or multi, in any combination.

POINT and MULTIPOINT

x <- st_sf(name = c("one clinic", "three clinics"),
           geometry = st_sfc(
             st_point(c(5, 4.2)),
             st_multipoint(rbind(c(2.6,5.6), c(5,2.6), c(7.6,5.6)))))

x$type <- st_geometry_type(x)
as.data.frame(x)[, c("name", "type", "geometry")]
           name       type                       geometry
1    one clinic      POINT                  POINT (5 4.2)
2 three clinics MULTIPOINT MULTIPOINT ((2.6 5.6), (5 2...

LINESTRING and MULTILINESTRING

x <- st_sf(name = c("one road", "two roads"),
           geometry = st_sfc(
             st_linestring(rbind(c(1,2), c(4,5.8), c(9,3))),
             st_multilinestring(list(rbind(c(1,6.5), c(4,4.6), c(9,6.2)),
                                     rbind(c(1,2.2), c(4.5,1.5), c(9,2.8))))))

x$type <- st_geometry_type(x)
as.data.frame(x)[, c("name", "type", "geometry")]
       name            type                       geometry
1  one road      LINESTRING   LINESTRING (1 2, 4 5.8, 9 3)
2 two roads MULTILINESTRING MULTILINESTRING ((1 6.5, 4 ...

POLYGON and MULTIPOLYGON

x <- st_sf(name = c("one state", "state + islands"),
           geometry = st_sfc(
             st_polygon(list(rbind(c(2,1.8), c(7,1.6), c(8,4.8), c(5,6.8), c(2,1.8)))),
             st_multipolygon(list(
               list(rbind(c(1.4,1.8), c(5.4,1.6), c(6.1,4.8), c(3.4,6.2), c(1.4,1.8))),
               list(rbind(c(7,5.2), c(8.3,4.9), c(8.6,6.4), c(7.2,6.6), c(7,5.2)))))))

x$type <- st_geometry_type(x)
as.data.frame(x)[, c("name", "type", "geometry")]
             name         type                       geometry
1       one state      POLYGON POLYGON ((2 1.8, 7 1.6, 8 4...
2 state + islands MULTIPOLYGON MULTIPOLYGON (((1.4 1.8, 5....

GEOMETRYCOLLECTION

ring <- rbind(c(5.2,4.6), c(9,4.6), c(9,7.2), c(5.2,7.2), c(5.2,4.6))

x <- st_sf(name = "one feature",
           geometry = st_sfc(st_geometrycollection(list(
             st_point(c(2.2, 6.2)),
             st_linestring(rbind(c(1,1.6), c(4.2,3.6), c(9,1.8))),
             st_polygon(list(ring))))))

x$type <- st_geometry_type(x)
as.data.frame(x)[, c("name", "type", "geometry")]
         name               type                       geometry
1 one feature GEOMETRYCOLLECTION GEOMETRYCOLLECTION (POINT (...

A point, a line and a polygon in one row. st_intersection() returns one when two layers meet at both an edge and a corner.

It Holds the MULTI Types Too

The parts do not have to be single geometries — MULTIPOINT, MULTILINESTRING and MULTIPOLYGON go in just as happily.

st_geometrycollection(list(
  st_multipoint(rbind(c(1,1), c(2,2))),
  st_multilinestring(list(rbind(c(0,0), c(1,1)), rbind(c(2,2), c(3,3)))),
  st_multipolygon(list(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0)))))))

You can mix the two families in one collection — a POINT next to a MULTIPOINT is fine:

st_geometrycollection(list(st_point(c(0,0)),
                           st_multipoint(rbind(c(1,1), c(2,2)))))

Two Things It Will Not Do

1. A collection cannot contain another collection. sf refuses:

st_geometrycollection(list(st_point(c(0,0)),
                           st_geometrycollection(list(st_point(c(1,1))))))
Error in `st_geometrycollection()`:
! all(cls[2, ] != "GEOMETRYCOLLECTION") is not TRUE

2. st_collection_extract() pulls out a family, not an exact type. Asking for "POLYGON" also returns any MULTIPOLYGON inside:

gc <- st_sfc(st_geometrycollection(list(
  st_multipolygon(list(list(rbind(c(0,0), c(1,0), c(1,1), c(0,0))))))))

st_geometry_type(st_collection_extract(gc, "POLYGON"))
[1] MULTIPOLYGON
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

Note

That is a convenience, not a bug — the three arguments "POINT", "LINESTRING" and "POLYGON" each mean “everything of that dimension”. Use st_cast() afterwards if you need one exact type.

Which Type Did You Get?

#Step1: Every geometry in a layer shares one type
table(as.character(st_geometry_type(mex1)))

MULTIPOLYGON 
          32 
#Step2: One row can still be many pieces
parts <- lengths(st_geometry(mex1))
names(parts) <- mex1$NAME_1
head(sort(parts, decreasing = TRUE), 3)
         Tamaulipas             Sinaloa Baja California Sur 
                261                 199                 180

Note

15 of the 32 states are a single piece, yet the type is still MULTIPOLYGON — every geometry in one layer must share a type. st_cast() converts when a function insists on one or the other.

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.

A Shapefile Is at Least Four Files

File What it holds Without it
.shp the geometry — every coordinate nothing to draw
.shx an index into .shp the file will not open
.dbf the attribute table, one row per shape no names, no values
.prj the CRS, written out as text R cannot place it on Earth
.cpg (optional) the text encoding of .dbf accents arrive mangled

.shp, .shx and .dbf are mandatory; .prj is needed in practice. They are linked only by a shared basename, gadm41_MEX_1. Rename one and the layer stops opening.

The .shp is 11.5 MB of an 11.6 MB download. The rest are a few kilobytes, but you still need them.

The .prj Is Just a Line of Text

prj <- readLines("../data/gadm41_MEX_shp/gadm41_MEX_1.prj", warn = FALSE)
prj <- gsub("],", "],\n", prj, fixed = TRUE)
prj <- gsub(",DATUM", ",\n  DATUM", prj, fixed = TRUE)
cat(gsub(",SPHEROID", ",\n    SPHEROID", prj, fixed = TRUE))
GEOGCS["GCS_WGS_1984",
  DATUM["D_WGS_1984",
    SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]]

That text is how st_read() knows the data is in WGS84. Without the .prj the geometry still loads, but with no CRS attached to it.

The .cpg and Your Accents

#The .cpg says how the text in the .dbf is encoded
readLines("../data/gadm41_MEX_shp/gadm41_MEX_1.cpg")
[1] "UTF-8"

UTF-8 is why Yucatán and Querétaro survive the trip.

A wrong .cpg gives you Yucat�n, and then a join on state names matches nothing.

The Shapefile Is Not the Only Format

Type Format Extension
Binary ESRI Shapefile .shp, .shx, .dbf, .prj
GeoPackage (GPKG) .gpkg
Plain text GeoJSON .json or .geojson
GPS Exchange Format (GPX) .gpx
Keyhole Markup Language (KML) .kml
Spatial database PostgreSQL / PostGIS

The shapefile has real limits: four files to keep together, and column names truncated to 10 characters.

st_read() reads all of them, so the format matters less than you might expect.

Same States, Three Formats

#Step1: st_read() does not care which container it is given
gj <- st_read("../data/mex_formats/mex_states.geojson", quiet = TRUE)
gp <- st_read("../data/mex_formats/mex_states.gpkg",    quiet = TRUE)

#Step2: the same sf object comes back either way
c(geojson = nrow(gj), gpkg = nrow(gp))
geojson    gpkg 
     32      32

Same 32 states, same columns: 0.45 MB as GeoJSON against 0.27 MB as GeoPackage. Plain text takes more room than a binary database.

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.

The Same Thing, With 32 Real Rows

The clinics were three rows we typed by hand. A real shapefile is the identical structure — just wider, and with polygons where we had points.

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

#Step2: The same two classes, the same sf_column -- 32 rows instead of 3
class(mex_states)
[1] "sf"         "data.frame"
dim(mex_states)
[1] 32 12
attr(mex_states, "sf_column")
[1] "geometry"

Nothing new to learn. You already know how to read this object.

The Same Thing, With 32 Real Rows

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

MULTIPOLYGON in the geometry column, so geom_sf() will draw filled shapes — exactly as on the POINT / LINESTRING / POLYGON slide earlier.

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. Across Mexico a degree of longitude varies by about 15% between Chiapas and the northern border — invisible at national scale. Project first only if you need the tolerance to mean exact metres.

st_simplify: Before and After

At 0.02, the value used above, the change is invisible at this scale. At 0.25 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 -- the data is already WGS84, so no coord_sf needed
ggplot() +
  geom_sf(data = mex_country)

Adding the States

ggplot() +
  geom_sf(data = mex_simple) +
  labs(title = "States of Mexico")

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

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 your own document:

```{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_simplify(mex, dTolerance = 0.02)   # degrees: the CRS it arrived in

ggplot() + geom_sf(data = mex)
```

The data stays in WGS84 and so does the map. There is only one CRS to keep track of.

Where AI Gets This Wrong

Four things it cannot see, all of which you met today:

  • it learned on sp, retired in 2023, so expect readOGR()
  • it cannot see your CRS, so it guesses the units: dTolerance = 2000, or 0.02?
  • it does not know your file’s vintage — GADM still says Distrito Federal
  • longitude and latitude reversed runs perfectly well, and puts Mexico in the Indian Ocean

None of these throws an error. Every one produces a map.

Note

Use it. Then check st_crs(), plot the result, and confirm the units. The disclosure rules are in the syllabus.

What Have We Learned?

  • 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
  • Keep data in WGS84; project only when you measure
  • An sf object is a dataframe with a geometry column, and geom_sf() draws whatever is in it
  • A shapefile is at least four files, and the geometry type decides what your map looks like