Building an Automated Python ETL Pipeline for Japan's 1,741 Municipalities Open Data (SSDSE-A)

Japan's Statistics Bureau (Ministry of Internal Affairs and Communications) and the National Statistics Center (NSTAC) publish an extraordinary open dataset called SSDSE (System of Social and Demographic Statistics).

Specifically, SSDSE-A covers all 1,741 municipalities across Japan, aggregating 128 comprehensive public indicators including demographics, households, economic activities, labor forces, municipal finance, healthcare facilities, and education.

However, if you simply pass this dataset into pd.read_csv(), you immediately encounter classic government data traps:

• A multi-level 3-row header that breaks schema inference.
• Legacy Shift-JIS (CP932) encoding that throws UnicodeDecodeError.
• Prefixed municipality codes (R01100) and dropping leading zeros, destroying relational keys.
• Missing derived analytical KPIs like habitable population density and aging ratios.

In this article, I'll walk through building a clean, automated Python pipeline to ingest, normalize, enrich, and export Japan's census open data into flat, modern formats (Parquet, BOM-UTF8 CSV, SQLite).

---

1. The Anatomy of the Raw CSV (And Its 4 Traps)

Here is a glimpse of the raw SSDSE-A-2026.csv published by NSTAC:

地域コード,都道府県,市区町村,A1101,A110101,A110102,...
,,,2020,2020,2020,...
,,,総人口,総人口（男）,総人口（女）,...
R01100,北海道,札幌市,1973395,918682,1054713,...
R01202,北海道,函館市,251084,113965,137119,...

The 4 Pitfalls:

1. 3-Row Composite Header:
   • Row 0: Statistical indicator codes (A1101)
   • Row 1: Survey survey years (2020)
   • Row 2: Human-readable Japanese names (総人口 = Total Population)
   • Running pd.read_csv() treats rows 1 and 2 as data, coercing every integer column into an object dtype.
2. The "R" Prefix and Zero-Dropping:
   • Standard Japanese municipality codes are 5 digits (JIS X 0402). For instance, Sapporo is 01100.
   • The raw data prefixes an "R" (R01100). If parsed as numeric, it drops the leading zero and becomes 1100, breaking relational joins with map shapefiles and corporate registries.
3. Encoding Quirks:
   • Distributed in CP932 (Windows-31J / Shift-JIS).
4. Raw Counts Without Relative Metrics:
   • The raw data has total population and habitable area, but lacks pre-calculated density or demographic ratios essential for GIS and spatial modeling.

---

2. Step 1: Automated Fetch Pipeline

We start by fetching the raw dataset directly from NSTAC's server. We enforce explicit timeouts and an authentic User-Agent to avoid connection blocks.

import os
from pathlib import Path
import requests

RAW_DIR = Path("data/raw")
RAW_DIR.mkdir(parents=True, exist_ok=True)
RAW_FILE = RAW_DIR / "SSDSE-A-2026.csv"
DATA_URL = "https://www.nstac.go.jp/files/SSDSE-A-2026.csv"

HEADERS = {
    "User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/120.0.0.0 Safari/537.36"
}

def fetch_raw_data() -> Path:
    if RAW_FILE.exists():
        print(f"File already exists: {RAW_FILE}")
        return RAW_FILE

    print(f"Downloading from {DATA_URL}...")
    res = requests.get(DATA_URL, headers=HEADERS, timeout=60)
    res.raise_for_status()
    with open(RAW_FILE, "wb") as f:
        f.write(res.content)
    print(f"Downloaded: {RAW_FILE} ({RAW_FILE.stat().st_size:,} bytes)")
    return RAW_FILE

---

3. Step 2: Unpacking Multi-Headers & Normalizing Codes

To handle the 3-row header, we parse the first 3 rows separately to build clean column names, and then read the data body starting from row 4.

import pandas as pd
import numpy as np

def clean_ssdse_dataset(csv_path: Path) -> pd.DataFrame:
    # 1. Parse header metadata
    df_meta = pd.read_csv(csv_path, encoding="cp932", nrows=3, header=None)
    codes = df_meta.iloc[0].values
    names = df_meta.iloc[2].values

    columns = ["municipality_code", "prefecture", "municipality_name"]
    for i in range(3, len(codes)):
        columns.append(str(names[i]).strip())

    # 2. Ingest data body (explicitly enforce string types on key columns to preserve leading zeros)
    dtype_dict = {0: str, 1: str, 2: str}
    df = pd.read_csv(
        csv_path,
        encoding="cp932",
        skiprows=3,
        header=None,
        dtype=dtype_dict
    )
    df.columns = columns

    # 3. Clean municipality codes: Strip 'R' and zero-pad to exactly 5 digits
    df["municipality_code"] = (
        df["municipality_code"]
        .astype(str)
        .str.strip()
        .str.lstrip("R")
        .str.zfill(5)
    )

    # 4. Safe numerical coercion (convert '-' or missing symbols to 0)
    for col in df.columns[3:]:
        df[col] = pd.to_numeric(df[col], errors="coerce").fillna(0)

    print(f"Cleaned {len(df):,} municipalities with {len(df.columns)} columns.")
    return df

---

4. Step 3: Engineering 10 Essential Demographic & Spatial KPIs

Roughly 70% of Japan is mountainous terrain. Comparing municipal density using gross surface area heavily skews commercial and spatial evaluations. Calculating Habitable Population Density (可住地人口密度) provides an accurate measure of living and market congestion.

We calculate 10 key metrics while guarding against division-by-zero errors:

def enrich_kpis(df: pd.DataFrame) -> pd.DataFrame:
    # Temporary zero-to-NaN conversion to prevent ZeroDivisionError
    pop = df["総人口"].replace(0, np.nan)
    area_tot = df["総面積"].replace(0, np.nan)
    area_hab = df["可住地面積"].replace(0, np.nan)
    hholds = df["世帯数"].replace(0, np.nan)

    kpis = pd.DataFrame(index=df.index)

    # Spatial densities (Persons per km2)
    kpis["kpi_pop_density_total"] = (df["総人口"] / area_tot).round(1).fillna(0)
    kpis["kpi_pop_density_habitable"] = (df["総人口"] / area_hab).round(1).fillna(0)

    # Age demographic ratios (%)
    kpis["kpi_aging_rate_pct"] = ((df["65歳以上人口"] / pop) * 100).round(2).fillna(0)
    kpis["kpi_late_elderly_rate_pct"] = ((df["75歳以上人口"] / pop) * 100).round(2).fillna(0)
    kpis["kpi_youth_rate_pct"] = ((df["15歳未満人口"] / pop) * 100).round(2).fillna(0)
    kpis["kpi_working_age_rate_pct"] = ((df["15～64歳人口"] / pop) * 100).round(2).fillna(0)

    # Single-person household ratio (%)
    kpis["kpi_single_household_ratio_pct"] = ((df["単独世帯数"] / hholds) * 100).round(2).fillna(0)

    # Net migration & natural balance (persons)
    kpis["kpi_net_migration_persons"] = (df["転入者数"] - df["転出者数"]).astype(int)
    kpis["kpi_natural_change_persons"] = (df["出生数"] - df["死亡数"]).astype(int)

    # Healthcare capacity (clinics + hospitals per 1,000 residents)
    kpis["kpi_medical_facilities_per_1k"] = (
        ((df["一般診療所数"] + df["一般病院数"]) / pop) * 1000
    ).round(2).fillna(0)

    # Place KPIs directly next to identifiers
    return pd.concat([df.iloc[:, :3], kpis, df.iloc[:, 3:]], axis=1)

---

5. Step 4: Multi-Format Export (Excel, Parquet & SQLite)

For analytical environments (DuckDB, PySpark, Polars), columnar formats like Parquet offer 10x-50x faster read speeds and built-in compression. Meanwhile, business users frequently open files in Microsoft Excel.

Pro Tip for Windows Excel: Saving as plain UTF-8 often triggers mojibake in Windows Excel. Always encode with utf-8-sig (BOM) and use \r\n line terminators.

import sqlite3

def export_pipeline(df: pd.DataFrame, out_dir: Path):
    out_dir.mkdir(parents=True, exist_ok=True)

    # 1. Windows Excel-Ready CSV (BOM-UTF8 with CRLF)
    csv_file = out_dir / "japan_municipal_census_ssdse_master.csv"
    df.to_csv(csv_file, index=False, encoding="utf-8-sig", lineterminator="\r\n")

    # 2. High-Performance Columnar Parquet (Snappy compression)
    pq_file = out_dir / "japan_municipal_census_ssdse_master.parquet"
    df.to_parquet(pq_file, index=False, compression="snappy")

    # 3. Indexed SQLite DB for Rapid SQL Queries
    db_file = out_dir / "municipal_master.db"
    conn = sqlite3.connect(db_file)
    df.to_sql("municipalities", conn, if_exists="replace", index=False)
    conn.execute("CREATE INDEX IF NOT EXISTS idx_muni_code ON municipalities (municipality_code);")
    conn.close()

    print(f"Export complete across all 3 formats in {out_dir}")

---

6. Advanced Spatial Engineering: Bridging Census Masters with Uber H3 Hexagonal Indexing

In modern spatial data science, joining point observations (e.g., store locations, medical clinics, real estate listings) with municipal administrative polygons is computationally expensive. Traditional GIS spatial operations like ST_Contains or gpd.sjoin require heavy geometric intersection math ($O(N \times M)$ complexity).

To enable lightning-fast spatial joins and hexagonal clustering, we can enrich each municipality with Uber's H3 Discrete Global Grid System (DGGS) based on municipal centroids (city hall coordinates):

import h3

def enrich_uber_h3_indices(df: pd.DataFrame, lat_col="latitude", lon_col="longitude") -> pd.DataFrame:
    """
    Assigns Uber H3 hierarchical hexagonal cells:
    - Resolution 7: ~5.16 km² (approx. municipal ward / town cluster)
    - Resolution 8: ~0.74 km² (neighborhood commercial scale)
    """
    df["h3_res7"] = df.apply(
        lambda r: h3.latlng_to_cell(r[lat_col], r[lon_col], 7) 
        if pd.notna(r.get(lat_col)) and pd.notna(r.get(lon_col)) else None, 
        axis=1
    )
    df["h3_res8"] = df.apply(
        lambda r: h3.latlng_to_cell(r[lat_col], r[lon_col], 8) 
        if pd.notna(r.get(lat_col)) and pd.notna(r.get(lon_col)) else None, 
        axis=1
    )
    return df

Why Uber H3 is a Game-Changer for Census Analysis:

1. $O(1)$ Hash Spatial Joins: Instead of running heavy point-in-polygon calculations across thousands of geo-coordinates, you simply match on integer/string H3 cell IDs in DuckDB or Pandas.
2. Seamless Neighbor Traversal (k-ring): Analyze surrounding demographic gravity by fetching immediate adjacent hexagons in one line: h3.grid_disk(h3_cell, k=1).
3. Web-Scale Interactive Mapping: H3 IDs render natively in high-performance WebGL frameworks like deck.gl and Kepler.gl without transferring heavy GeoJSON polygon payloads.

---

Summary & Impact

By executing this automated pipeline:

• We transform an awkward 3-header Japanese government CSV into a clean 1,741 rows × 140 columns flat dataset.
• The 5-digit JIS municipality code is standardized, allowing effortless joins with Japan's land valuation registries, healthcare establishment directories, and GIS boundary GeoJSONs.
• By optionally projecting municipal coordinates into Uber H3 hexagonal grid cells, spatial analysts unlock $O(1)$ relational joins and high-performance WebGL heatmaps.

Happy coding with Japanese open data!