Getting the data into Databricks

Stage zero: optional, and only if you want to reproduce the example yourself

You almost certainly do not need this page. The rest of this site assumes your data is already in Databricks, because that is the situation it was written for. This page exists only so the worked example can be reproduced from scratch, and it writes to a catalog rather than reading from one.

Skip to Getting connected and finding the tables unless you specifically want to load this data yourself.

Everything on this page is shown but not run. It downloads roughly 32.5 million rows over 25 to 40 minutes, stages about 6 GB of transient files, and needs write access to a catalog. Running it as part of a render would be wrong, so the chunks here are set not to execute.

You need CREATE TABLE and CREATE VOLUME on a schema you own. If you do not have that, this page is not available to you, and nothing else on the site requires it.

The shape of the load, and why it is not dbWriteTable()

The obvious route is to build a data frame in R and write it row by row over DBI. That works up to about 50,000 rows and warns well before it, which is three orders of magnitude short of what is needed here.

The route that scales goes through a file. Write Parquet locally, upload it to a Unity Catalog volume, then have Databricks read the volume into a Delta table. The data crosses once, as a compressed file, and the table is created by the server reading its own storage.

library(DBI)
library(dplyr)
library(tibble)
library(purrr)
library(stringr)
library(glue)
library(cli)
library(readr)
library(arrow)
library(sf)
library(brickster)

catalog <- Sys.getenv("DATABRICKS_CATALOG")
schema  <- Sys.getenv("DATABRICKS_SCHEMA")

con <- dbConnect(
  odbc::databricks(),
  httpPath = Sys.getenv("DATABRICKS_PATH"),
  DefaultStringColumnLength = 65535
)

Two licences, and they are different. Getting either wrong is a licensing error rather than a tidiness one, so each is written next to the table it belongs to.

licence_ogl <- paste(
  "Contains public sector information licensed under the",
  "Open Government Licence v3.0.",
  "Source: Environment Agency hydrology data."
)

licence_ccby <- paste(
  "Storm overflow outfall locations (c) the water",
  "companies, licensed under CC BY 4.0 and published",
  "via Stream / Water UK's National Storm Overflow Hub."
)

The volume is the staging area, and creating it is idempotent so this is safe to re-run.

volume <- "raw"
volume_root <- glue("/Volumes/{catalog}/{schema}/{volume}")

have <- db_uc_volumes_list(catalog, schema)$volumes |>
  map_chr(\(v) v$name %||% NA_character_)

if (!volume %in% have) {
  db_uc_volumes_create(
    catalog, schema, volume,
    volume_type = "MANAGED",
    comment = "Parquet staging for the worked example"
  )
}

One function that does the whole round trip

Every table below goes the same way, so the mechanics are worth writing once. Note the DROP TABLE and the directory clear: read_files() reads every file in the directory it is pointed at, so a re-run with a fresh temporary name would silently double the row count rather than replacing it.

load_table <- function(df, table, con, comment = NULL) {
  stopifnot(is.data.frame(df), nrow(df) > 0)

  local_pq <- tempfile(fileext = ".parquet")
  on.exit(unlink(local_pq), add = TRUE)
  write_parquet(df, local_pq, compression = "zstd")

  remote_dir <- glue("{volume_root}/{table}")
  tryCatch(
    db_volume_list(remote_dir)$contents,
    error = \(e) NULL
  ) |>
    keep(\(f) isFALSE(f$is_directory)) |>
    walk(\(f) db_volume_delete(f$path))

  db_volume_write(
    glue("{remote_dir}/data.parquet"),
    file = local_pq,
    overwrite = TRUE,
    progress = FALSE
  )

  fq <- glue("`{catalog}`.`{schema}`.`{table}`")
  dbExecute(con, glue("DROP TABLE IF EXISTS {fq}"))
  dbExecute(con, glue("
    CREATE TABLE {fq}
    AS SELECT * EXCEPT (_rescued_data)
    FROM read_files(
      '{remote_dir}/',
      format => 'parquet'
    )
  "))

  # The comment is where the licence lives, so it
  # travels with the table rather than only with the
  # code that loaded it.
  if (!is.null(comment)) {
    quoted <- dbQuoteString(con, comment)
    dbExecute(
      con,
      glue("COMMENT ON TABLE {fq} IS {quoted}")
    )
  }

  n <- glue("SELECT count(*) AS n FROM {fq}") |>
    dbGetQuery(conn = con) |>
    pull(n) |>
    as.numeric()

  stopifnot(n == nrow(df))
  cli_alert_success(
    "{table}: {format(n, big.mark = ',')} rows loaded"
  )
  invisible(n)
}

The fact table, in windows

The hydrology API caps a single response at two million rows, so ten years of daily readings has to be paged. Six-month windows come to roughly 1.5 million rows each, comfortably under the cap.

Each window is streamed straight to Parquet rather than accumulated in memory, which is what keeps peak R memory small while handling 32.5 million rows.

stage_dir <- Sys.getenv(
  "WATER_STAGE_DIR",
  "~/water-demo-staging"
) |>
  path.expand()

dir.create(
  file.path(stage_dir, "parquet"),
  recursive = TRUE,
  showWarnings = FALSE
)

windows <- tibble(year = 2015:2024) |>
  reframe(
    tag  = glue("{year}{c('H1', 'H2')}"),
    from = glue("{year}-{c('01-01', '07-01')}"),
    to   = glue("{year}-{c('06-30', '12-31')}"),
    .by  = year
  )

windows |>
  pwalk(\(tag, from, to, ...) {
    csv <- file.path(stage_dir, glue("{tag}.csv"))
    if (!file.exists(csv)) {
      url <- glue(
        "https://environment.data.gov.uk/hydrology/",
        "data/readings.csv?mineq-date={from}",
        "&max-date={to}&period=86400&_limit=2000000"
      )
      curl::curl_download(url, csv)
    }

    csv |>
      read_csv(show_col_types = FALSE, progress = FALSE) |>
      write_parquet(
        file.path(
          stage_dir,
          "parquet",
          glue("{tag}.parquet")
        )
      )
  })

The twenty Parquet files then load as one table, again by pointing read_files() at a directory rather than sending rows.

file.path(stage_dir, "parquet") |>
  list.files(full.names = TRUE) |>
  walk(\(f) db_volume_write(
    glue("{volume_root}/hydrology_readings/{basename(f)}"),
    file = f, overwrite = TRUE, progress = FALSE))

fq <- glue("`{catalog}`.`{schema}`.hydrology_readings")
dbExecute(con, glue("DROP TABLE IF EXISTS {fq}"))
readings_dir <- glue("{volume_root}/hydrology_readings/")

dbExecute(con, glue("
  CREATE TABLE {fq}
  CLUSTER BY (measure)
  AS SELECT * EXCEPT (_rescued_data)
  FROM read_files('{readings_dir}', format => 'parquet')
"))

quoted <- dbQuoteString(con, licence_ogl)
dbExecute(con, glue("COMMENT ON TABLE {fq} IS {quoted}"))

CLUSTER BY (measure) is worth the one extra clause. Every query in the example joins readings to their measure, so clustering on that column is what makes the joins in Exploring and reducing the big table cheap rather than a full scan.

The dimension tables

Stations and measures are small and arrive as ordinary CSV. The station table is British National Grid: easting and northing, with no longitude or latitude anywhere in it, which is the first thing to catch anyone who assumes otherwise.

fetch_csv <- function(url, ...) {
  tf <- tempfile(fileext = ".csv")
  curl::curl_download(
    url,
    tf,
    handle = curl::new_handle(accept_encoding = "gzip")
  )
  on.exit(unlink(tf), add = TRUE)
  read_csv(
    tf,
    show_col_types = FALSE,
    progress = FALSE,
    ...
  )
}

stations_url <- paste0(
  "https://environment.data.gov.uk/hydrology/id/",
  "stations.csv?_limit=20000"
)

stations <- fetch_csv(stations_url) |>
  transmute(
    station_id = notation, station_name = label,
    easting = as.numeric(easting),
    northing = as.numeric(northing),
    river_name = riverName,
    catchment_name = catchmentName,
    town,
    catchment_area_km2 = as.numeric(catchmentArea),
    date_opened = dateOpened, status = `status.label`
  ) |>
  filter(!is.na(easting), !is.na(northing))

load_table(
  stations,
  "hydrology_stations",
  con,
  glue(
    "EA hydrology stations, EPSG:27700 ",
    "easting/northing. {licence_ogl}"
  )
)

station_prefix <- paste0(
  "http://environment.data.gov.uk/",
  "hydrology/id/stations/"
)
measures_url <- paste0(
  "https://environment.data.gov.uk/hydrology/id/",
  "measures.csv?_limit=100000"
)

measures <- fetch_csv(measures_url) |>
  transmute(
    notation,
    station_id = str_remove(station, fixed(station_prefix)),
    parameter, parameter_name = parameterName,
    period = as.integer(period), period_name = periodName,
    value_statistic = valueStatistic, unit_name = unitName
  )

load_table(
  measures,
  "hydrology_measures",
  con,
  glue(
    "EA measure catalogue. `notation` joins to ",
    "hydrology_readings.measure. {licence_ogl}"
  )
)

notation joining to measure is the join the whole example depends on, and it is worth noticing now: station_id is on the measures table, not on the readings. There is no shortcut from a reading to a station.

The polygons, stored as WKB

The catchments arrive as a GeoPackage and have to become a column in a table. Geometry goes in as WKB in a BINARY column rather than as text.

That choice is worth a sentence because it has consequences later. Encoding 4,080 multipolygons takes a fraction of a second as WKB against roughly two minutes as WKT, and the result is smaller. The cost is that reading it back needs a connection that carries BINARY intact, which is the subject of the next page.

gpkg <- Sys.glob(file.path(stage_dir, "spatial", "*.gpkg"))
catch_sf <- st_read(gpkg[1], quiet = TRUE)
stopifnot(st_crs(catch_sf)$epsg == 27700)

catchments <- catch_sf |>
  st_drop_geometry() |>
  transmute(
    water_body_id, water_body_name, water_body_type,
    operational_catchment,
    opcat_id,
    management_catchment,
    mancat_id,
    river_basin_district, ea_area_name,
    classification_year = as.integer(classification_year),
    overall_water_body_class,
    ecological_class,
    chemical_class,
    water_body_area_m2 = as.numeric(water_body_area_m2)
  )

# Assigned rather than mutate()d: a WKB list column is
# not something dplyr::mutate() accepts as a vector.
catchments$geom_wkb <- catch_sf |>
  st_geometry() |>
  st_as_binary() |>
  unclass()

load_table(
  catchments,
  "wfd_catchments",
  con,
  glue(
    "WFD River Water Body Catchments Cycle 3. geom_wkb ",
    "is WKB MULTIPOLYGON in EPSG:27700. {licence_ogl}"
  )
)

water_body_area_m2 is worth keeping even though it is derivable from the geometry. It becomes the independent check that geometry survived the trip back, which When the simulation carries the geometry relies on.

The overflows, under a different licence

Nine English water companies publish outfall locations through a shared hub, each as its own service, and the field names differ in case between them. These are longitude and latitude in EPSG:4326, unlike everything else here, and that mismatch is deliberate: reconciling it is part of what The spatial join has to do.

Scottish Water is excluded, because it reserves all rights.

# Nine FeatureServers, differing only in host, org id and
# service name, so the varying parts are the table and the
# shared shape is written once.
hub <- tribble(
  ~company, ~host, ~org, ~service,
  "Severn Trent Water",
    "1", "NO7lTIlnxRMMG9Gw",
    "Severn_Trent_Water_Storm_Overflow_Activity",
  "United Utilities",
    "5", "5eoLvR0f8HKb7HWP",
    "United_Utilities_Storm_Overflow_Activity",
  "Yorkshire Water",
    "-eu1", "1WqkK5cDKUbF0CkH",
    "Yorkshire_Water_Storm_Overflow_Activity",
  "Northumbrian Water",
    "-eu1", "MSNNjkZ51iVh8yBj",
    "Northumbrian_Water_Storm_Overflow_Activity_2_view",
  "Anglian Water",
    "3", "VCOY1atHWVcDlvlJ",
    "stream_service_outfall_locations_view",
  "Wessex Water",
    "", "3SZ6e0uCvPROr4mS",
    "Wessex_Water_Storm_Overflow_Activity",
  "South West Water",
    "-eu1", "OMdMOtfhATJPcHe3",
    "NEH_outlets_PROD",
  "Southern Water",
    "-eu1", "6qJmARkS2dt2IjVA",
    "SouthernWater_StormOverflowActivity_PROD_view",
  "Thames Water",
    "2", "g6o32ZDQ33GpCIu3",
    "Thames_Water_Storm_Overflow_Activity_(Production)_view"
) |>
  mutate(
    url = glue(
      "https://services{host}.arcgis.com/{org}/",
      "arcgis/rest/services/{service}/FeatureServer"
    )
  )

# These public services time out intermittently, so
# every page is retried rather than allowed to fail
# the whole load.
read_json_once <- function(url) {
  tf <- tempfile()
  on.exit(unlink(tf), add = TRUE)
  curl::curl_download(
    url,
    tf,
    handle = curl::new_handle(timeout = 180)
  )
  jsonlite::fromJSON(tf, simplifyVector = FALSE)
}

# insistently() wraps the call in a backoff rather than
# open-coding the retry, so the function above stays
# about fetching and says nothing about failure.
get_json <- insistently(
  read_json_once,
  rate = rate_backoff(pause_base = 3, max_times = 5),
  quiet = FALSE
)

as_overflow_rows <- function(feats) {
  feats |>
    map(\(f) {
      set_names(f$attributes, tolower(names(f$attributes)))
    }) |>
    map(\(a) tibble(
      overflow_id = as.character(a[["id"]] %||% NA),
      longitude   = as.numeric(a[["longitude"]] %||% NA),
      latitude    = as.numeric(a[["latitude"]] %||% NA),
      status      = as.integer(a[["status"]] %||% NA),
      receiving_water_course =
        as.character(a[["receivingwatercourse"]] %||% NA)
    )) |>
    list_rbind()
}

# The API pages, and does not report a total, so the
# only way to know you have reached the end is a short
# page. accumulate() is the wrong shape for that: the
# stopping condition depends on what came back rather
# than on a counter.
fetch_company <- function(company, base, page = 1000L) {
  out <- list()
  offset <- 0L
  repeat {
    js <- get_json(glue(
      "{base}/0/query?where=1%3D1&outFields=*",
      "&outSR=4326&returnGeometry=false&f=json",
      "&resultRecordCount={page}&resultOffset={offset}"
    ))
    feats <- js$features
    if (!length(feats)) break
    out <- append(out, list(as_overflow_rows(feats)))
    offset <- offset + length(feats)
    if (length(feats) < page) break
  }
  list_rbind(out) |> mutate(company = company)
}

overflows <- hub |>
  pmap(\(company, url, ...) fetch_company(company, url)) |>
  list_rbind() |>
  filter(!is.na(longitude), !is.na(latitude)) |>
  select(
    overflow_id,
    company,
    longitude,
    latitude,
    status,
    receiving_water_course
  )

load_table(
  overflows,
  "storm_overflows",
  con,
  glue(
    "Storm overflow outfall locations, nine English ",
    "water companies. EPSG:4326 lon/lat: needs ",
    "transforming to 27700 to join to the catchments. ",
    "{licence_ccby}"
  )
)

The licence on this table is CC BY 4.0, not OGL, and it is operational rather than regulatory data. That difference has to survive into anything you publish from it, which is why it goes into the table comment here rather than living only in a script.

Check it landed

Two checks, and the second is the one that matters. Row counts confirm the load; decoding a sample of geometry and comparing its area against the column the table already carries confirms the geometry itself is intact rather than merely present.

glue("SHOW TABLES IN `{catalog}`.`{schema}`") |>
  dbGetQuery(conn = con) |>
  as_tibble()

probe <- glue("
  SELECT
    base64(geom_wkb) AS g,
    water_body_area_m2 AS declared
  FROM `{catalog}`.`{schema}`.wfd_catchments
  WHERE geom_wkb IS NOT NULL
    AND water_body_area_m2 IS NOT NULL
  LIMIT 200
") |>
  dbGetQuery(conn = con) |>
  as_tibble()

geom <- probe$g |>
  map(jsonlite::base64_dec) |>
  structure(class = "WKB") |>
  st_as_sfc(EWKB = FALSE, crs = 27700)

declared <- sum(probe$declared)
measured <- sum(as.numeric(st_area(geom)))
err <- abs(measured - declared) / declared

cli_alert_info(
  "decoded area differs from declared by \\
   {round(100 * err, 4)}%"
)

A difference of a few thousandths of a percent means the geometry round-tripped. A large difference, or a decode error, means truncation rather than corruption, and the next page explains why.

dbDisconnect(con)

Next

Getting connected and finding the tables.

---

This page rests on: dbWriteTable() through the DBI backend works to about 50,000 rows and warns from about 20,000, so bulk loading goes through a volume and read_files(); read_files() reads every file in the directory it is given, so a re-run with a fresh filename doubles the row count; encoding catchment geometry as WKB is far faster than WKT and produces a smaller value; hydrology_stations carries easting/northing in EPSG:27700 and no lon/lat, while storm_overflows carries lon/lat in EPSG:4326; station_id lives on the measures table rather than the readings table; the storm overflow locations are CC BY 4.0 rather than OGL.