Examples: Census + Maps

This section combines 2018 Census data with geographic visualizations. The general pattern is: aggregate microdata first, then merge with geometry.

Tip

Key rule: Always aggregate first (with pandas in Python, with dplyr in R), then merge geometry onto the summary rows. Never use geometry= directly on large microdata — it attaches a polygon to every row and is very slow.

---

Internet access by municipality

What percentage of households has internet access in each municipality of Sacatepéquez?

Python

import geoquetzal as gq

# Census data
df = gq.hogares(departamento="Sacatepequez")

# PCH9_I: Does this household have internet? (1=Yes, 2=No)
acceso_internet = (
    df.groupby("MUNICIPIO")["PCH9_I"]
    .apply(lambda x: (x == 1).mean() * 100)
    .round(1)
    .reset_index(name="pct_internet")
)

# Geometry
munis = gq.municipios("Sacatepequez")

# Join using INE municipality code
resultado = munis.merge(acceso_internet, left_on="codigo_muni", right_on="MUNICIPIO")

# Interactive map
resultado.explore(
    column="pct_internet",
    cmap="YlGnBu",
    tooltip=["municipio", "pct_internet"],
    tooltip_kwds={"aliases": ["Municipality", "% Internet"]},
    tiles="CartoDB positron",
    style_kwds={"weight": 1, "color": "white", "fillOpacity": 0.8},
)

Make this Notebook Trusted to load map: File -> Trust Notebook

R

library(geoquetzal)
library(dplyr)
library(mapview)

# Census data
df <- hogares(departamento = "Sacatepequez")

# PCH9_I: Does this household have internet? (1=Yes, 2=No)
acceso_internet <- df |>
  group_by(MUNICIPIO) |>
  summarise(pct_internet = round(mean(PCH9_I == 1, na.rm = TRUE) * 100, 1))

# Geometry + join using INE municipality code
resultado <- merge(
  municipios("Sacatepequez"), acceso_internet,
  by.x = "codigo_muni", by.y = "MUNICIPIO"
)

# Interactive map
mapview(resultado,
        zcol        = "pct_internet",
        col.regions = RColorBrewer::brewer.pal(9, "YlGnBu"),
        layer.name  = "% Internet",
        map.types   = "CartoDB.Positron",
        label       = "municipio")

---

Dominant ethnic group by municipality

Which ethnic group is the majority in each municipality of Huehuetenango?

Python

import geoquetzal as gq

df = gq.personas(departamento="Huehuetenango")

# PCP12: ethnic self-identification (mode = majority group)
etnia_dominante = (
    df.groupby("MUNICIPIO")["PCP12"]
    .agg(lambda x: x.mode()[0])
    .reset_index(name="etnia_cod")
)

# Replace codes with labels
valores = gq.describe_personas("PCP12")["valores"]
etnia_dominante["etnia"] = etnia_dominante["etnia_cod"].map(valores)

# Geometry
munis = gq.municipios("Huehuetenango")
resultado = munis.merge(etnia_dominante, left_on="codigo_muni", right_on="MUNICIPIO")

# Interactive map
resultado.explore(
    column="etnia",
    tooltip=["municipio", "etnia"],
    tooltip_kwds={"aliases": ["Municipality", "Majority Ethnic Group"]},
    tiles="CartoDB positron",
    style_kwds={"weight": 1, "color": "white"},
)

Make this Notebook Trusted to load map: File -> Trust Notebook

R

library(geoquetzal)
library(dplyr)
library(mapview)

df <- personas(departamento = "Huehuetenango")

# PCP12: ethnic self-identification (mode = majority group)
etnia_dominante <- df |>
  group_by(MUNICIPIO) |>
  summarise(etnia_cod = as.integer(names(which.max(table(PCP12))))) |>
  mutate(etnia = recode(as.character(etnia_cod),
    "1" = "Maya", "2" = "Garífuna", "3" = "Xinka",
    "4" = "Afrodescendiente", "5" = "Ladino", "6" = "Extranjero"))

# Geometry + join
resultado <- merge(
  municipios("Huehuetenango"), etnia_dominante,
  by.x = "codigo_muni", by.y = "MUNICIPIO"
)

# Interactive map
mapview(resultado,
        zcol       = "etnia",
        map.types  = "CartoDB.Positron",
        label      = "municipio",
        layer.name = "Majority Ethnic Group")

---

Cumulative emigration by department

How many people emigrated from each department between 2002 and 2018?

Python

import geoquetzal as gq

df = gq.emigracion()

emigracion = (
    df.groupby("DEPARTAMENTO")
    .size()
    .reset_index(name="emigrants")
)

deptos = gq.departamentos()
resultado = deptos.merge(emigracion, left_on="codigo_depto", right_on="DEPARTAMENTO")

resultado.explore(
    column="emigrants",
    cmap="OrRd",
    tooltip=["departamento", "emigrants"],
    tooltip_kwds={"aliases": ["Department", "Emigrants"]},
    tiles="CartoDB positron",
    style_kwds={"weight": 1, "color": "white", "fillOpacity": 0.8},
)

Make this Notebook Trusted to load map: File -> Trust Notebook

R

library(geoquetzal)
library(dplyr)
library(mapview)

df <- emigracion()

emigracion_agg <- df |>
  group_by(DEPARTAMENTO) |>
  summarise(emigrants = n())

resultado <- merge(
  departamentos(), emigracion_agg,
  by.x = "codigo_depto", by.y = "DEPARTAMENTO"
)

mapview(resultado,
        zcol        = "emigrants",
        col.regions = RColorBrewer::brewer.pal(9, "OrRd"),
        layer.name  = "Emigrants",
        map.types   = "CartoDB.Positron",
        label       = "departamento")

---

Sub-municipal choropleth with Voronoi

Internet access at lugar poblado level in Antigua Guatemala — the first sub-municipal census visualization available for Guatemala.

Python

import geoquetzal as gq

# Voronoi polygons (approximation of sub-municipal boundaries)
vor = gq.voronoi_lugares_poblados(municipio="Antigua Guatemala")

   ℹ 2 lugares poblados excluidos por coordenadas nulas (códigos terminados en 999 — asentamientos sin nombre oficial).
   ✓ 57 polígonos Voronoi generados

# Pre-aggregated census data by lugar poblado
lp = gq.lugares_poblados(municipio="Antigua Guatemala")
lp = lp.drop(columns=["nombre", "lat", "longitud", "area"])

# Join
gdf = vor.merge(lp, on=["departamento", "municipio", "lugar_poblado"])

# % of households with internet
gdf["pct_internet"] = (
    gdf["pch9_i_si"] / (gdf["pch9_i_si"] + gdf["pch9_i_no"]) * 100
).round(1)

# Interactive map with hover tooltips
gdf.explore(
    column="pct_internet",
    cmap="YlGnBu",
    tooltip=["nombre", "pct_internet"],
    tooltip_kwds={"aliases": ["Lugar Poblado", "% with Internet"]},
    popup=["nombre", "pct_internet", "pch9_i_si", "pch9_i_no"],
    popup_kwds={"aliases": [
        "Lugar Poblado", "% with Internet",
        "HH with Internet", "HH without Internet"
    ]},
    legend=True,
    tiles="CartoDB positron",
    style_kwds={"weight": 0.5, "color": "white"},
)

Make this Notebook Trusted to load map: File -> Trust Notebook

R

library(geoquetzal)
library(mapview)

# Voronoi polygons (approximation of sub-municipal boundaries)
vor <- voronoi_lugares_poblados(municipio = "Antigua Guatemala")

# Pre-aggregated census data by lugar poblado
lp <- lugares_poblados(municipio = "Antigua Guatemala")
lp <- lp[, !names(lp) %in% c("nombre", "lat", "longitud", "area")]

# Join
gdf <- merge(vor, lp, by = c("departamento", "municipio", "lugar_poblado"))

# % of households with internet
gdf$pct_internet <- round(
  gdf$pch9_i_si / (gdf$pch9_i_si + gdf$pch9_i_no) * 100, 1
)

# Interactive map
mapview(gdf,
        zcol        = "pct_internet",
        col.regions = RColorBrewer::brewer.pal(9, "YlGnBu"),
        layer.name  = "% Internet",
        map.types   = "CartoDB.Positron",
        lwd         = 0.5,
        color       = "white",
        label       = "nombre")

Note

Voronoi polygons are approximations for visualization — INE does not publish official lugar poblado boundaries. Each polygon represents the zone of influence of the lugar poblado centroid within its municipio.

---

Maya population percentage by lugar poblado

Python

import geoquetzal as gq

vor = gq.voronoi_lugares_poblados(departamento="Sacatepequez")

   ℹ 22 lugares poblados excluidos por coordenadas nulas (códigos terminados en 999 — asentamientos sin nombre oficial).
   ✓ 237 polígonos Voronoi generados

lp  = gq.lugares_poblados(departamento="Sacatepequez")
lp  = lp.drop(columns=["nombre", "lat", "longitud", "area"])

gdf = vor.merge(lp, on=["departamento", "municipio", "lugar_poblado"])

# Use total answered as denominator (exclude non-responses)
pcp12_cols = [c for c in gdf.columns if c.startswith("pcp12_") and c != "pcp12_nulo"]
gdf["pcp12_total_respondido"] = gdf[pcp12_cols].sum(axis=1)
gdf["pct_maya"] = (
    gdf["pcp12_maya"] / gdf["pcp12_total_respondido"] * 100
).round(1)

gdf.explore(
    column="pct_maya",
    cmap="Purples",
    tooltip=["nombre", "pct_maya"],
    tooltip_kwds={"aliases": ["Lugar Poblado", "% Maya Population"]},
    legend=True,
    tiles="CartoDB positron",
    style_kwds={"weight": 0.5, "color": "white"},
)

Make this Notebook Trusted to load map: File -> Trust Notebook

R

library(geoquetzal)
library(mapview)

vor <- voronoi_lugares_poblados(departamento = "Sacatepequez")
lp  <- lugares_poblados(departamento = "Sacatepequez")
lp  <- lp[, !names(lp) %in% c("nombre", "lat", "longitud", "area")]

gdf <- merge(vor, lp, by = c("departamento", "municipio", "lugar_poblado"))

# Use total answered as denominator (exclude non-responses)
pcp12_cols <- grep("^pcp12_", names(gdf), value = TRUE)
pcp12_cols <- pcp12_cols[pcp12_cols != "pcp12_nulo"]
gdf$pcp12_total_respondido <- rowSums(as.data.frame(gdf)[, pcp12_cols], na.rm = TRUE)
gdf$pct_maya <- round(gdf$pcp12_maya / gdf$pcp12_total_respondido * 100, 1)

mapview(gdf,
        zcol        = "pct_maya",
        col.regions = RColorBrewer::brewer.pal(9, "Purples"),
        layer.name  = "% Maya Population",
        map.types   = "CartoDB.Positron",
        lwd         = 0.5,
        color       = "white",
        label       = "nombre")

---

Dirt floor as a socioeconomic indicator

Python

import geoquetzal as gq

vor = gq.voronoi_lugares_poblados(departamento="Sacatepequez")

   ℹ 22 lugares poblados excluidos por coordenadas nulas (códigos terminados en 999 — asentamientos sin nombre oficial).
   ✓ 237 polígonos Voronoi generados

lp  = gq.lugares_poblados(departamento="Sacatepequez")
lp  = lp.drop(columns=["nombre", "lat", "longitud", "area"])

gdf = vor.merge(lp, on=["departamento", "municipio", "lugar_poblado"])

# pcv5_tierra = housing units with dirt floor
gdf["pct_dirt_floor"] = (
    gdf["pcv5_tierra"] / gdf["poblacion_total"] * 100
).round(1)

gdf.explore(
    column="pct_dirt_floor",
    cmap="YlOrBr",
    tooltip=["nombre", "pct_dirt_floor"],
    tooltip_kwds={"aliases": ["Lugar Poblado", "% Dirt Floor"]},
    legend=True,
    tiles="CartoDB positron",
    style_kwds={"weight": 0.5, "color": "white"},
)

Make this Notebook Trusted to load map: File -> Trust Notebook

R

library(geoquetzal)
library(mapview)

vor <- voronoi_lugares_poblados(departamento = "Sacatepequez")
lp  <- lugares_poblados(departamento = "Sacatepequez")
lp  <- lp[, !names(lp) %in% c("nombre", "lat", "longitud", "area")]

gdf <- merge(vor, lp, by = c("departamento", "municipio", "lugar_poblado"))

# pcv5_tierra = housing units with dirt floor
gdf$pct_dirt_floor <- round(gdf$pcv5_tierra / gdf$poblacion_total * 100, 1)

mapview(gdf,
        zcol        = "pct_dirt_floor",
        col.regions = RColorBrewer::brewer.pal(9, "YlOrBr"),
        layer.name  = "% Dirt Floor",
        map.types   = "CartoDB.Positron",
        lwd         = 0.5,
        color       = "white",
        label       = "nombre")