Hard to Measure Well: Can Feasible Policies Reduce Methane Emissions?

import sys
import itertools
from pathlib import Path
from concurrent.futures import ProcessPoolExecutor
import pandas as pd
import pyarrow
import pyarrow.parquet as pq
import logging
from functools import partial

if sys.version_info < (3, 7):
    raise AssertionError("Need python version 3.7 or higher")

# Create a logger to output any messages we might have...
logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)
handler = logging.StreamHandler()
formatter = logging.Formatter("%(message)s")
handler.setFormatter(formatter)
logger.addHandler(handler)


def dtypes_to_schema(dtype_dict):
    schema = []
    type_conversion = {
        "Int64": pyarrow.int64(),
        "Int32": pyarrow.int32(),
        "int64": pyarrow.int64(),
        "int32": pyarrow.int32(),
        "int": pyarrow.int32(),
        "float64": pyarrow.float64(),
        "float32": pyarrow.float32(),
        "float": pyarrow.float32(),
        "str": pyarrow.string(),
        "O": pyarrow.string(),
        # Note that the arguments to pyarrow.dictionary changed between v0.13.0
        # and v0.14.1
        "category": pyarrow.dictionary(
            pyarrow.int32(), pyarrow.string(), ordered=False
        ),
        "date32": pyarrow.date32(),  # Note: date32 isn't a python type
    }
    for field_name, dtype in dtype_dict.items():
        type = type_conversion[dtype]
        schema.append(pyarrow.field(field_name, type, nullable=True))
    return pyarrow.schema(schema)


def extension_int_to_float(df, exclude=[]):
    """
    You can delete this function once this issue is resolved:
    https://issues.apache.org/jira/browse/ARROW-5379
    """
    extension_int_types = {
        pd.Int8Dtype(),
        pd.Int16Dtype(),
        pd.Int32Dtype(),
        pd.Int64Dtype(),
    }
    new_types = {
        col: "float64"
        for col in df.columns
        if df[col].dtypes in extension_int_types and col not in exclude
    }
    return df.astype(new_types)


def fix_monthly_dates(df):
    df = coerce_to_date(df, ["date"], drop_bad=True)
    # Note: these are sometimes NA
    df["month"] = df["date"].dt.month.astype("Int32")
    df["year"] = df["date"].dt.year.astype("Int32")
    del df["date"]
    return df


def read_csv(filepath_or_buffer, *args, **kwargs):
    """Thin wrapper around pd.read_csv"""
    logger.info(f"  Reading {filepath_or_buffer}")
    return pd.read_csv(filepath_or_buffer, *args, **kwargs)


def convert_monthly(csv_file, out_dir):
    csv_file = Path(csv_file)
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    # Note: we're using python 3.7+, so dicts keep their order
    pd_dtypes = {
        "entity_id": "Int32",
        "api": "str",
        "api_list": "str",
        "date": "str",
        "oil": "str",  # will convert to float later
        "gas": "str",  # will convert to float later
        "water": "str",  # will convert to float later
        "well_count": "str",  # will convert to Int later
        "days": "str",  # will convert to Int later
        "daily_avg_oil": "str",  # will convert to float later
        "daily_avg_gas": "str",  # will convert to float later
        "daily_avg_water": "str",  # will convert to float later
        "reservoir": "str",
        "well_name": "str",
        "well_number": "str",
        "operator_alias": "str",
        "production_type": "str",
        "production_status": "str",
        "entity_type": "str",
    }
    # Need to read these as str, then convert to float, because of errors.
    float_cols = [
        "oil",
        "gas",
        "water",
        "daily_avg_oil",
        "daily_avg_gas",
        "daily_avg_water",
    ]
    int_cols = ["well_count", "days"]
    output_dtypes = pd_dtypes.copy()
    del output_dtypes["date"]
    output_dtypes["year"] = "int32"
    output_dtypes["month"] = "int32"
    for c in float_cols:
        output_dtypes[c] = "float"
    for c in int_cols:
        output_dtypes[c] = "int"
    pq_schema = dtypes_to_schema(output_dtypes)

    csv_iter = read_csv(
        csv_file,
        dtype=pd_dtypes,
        chunksize=3_000_000,
        index_col=False,
        names=list(pd_dtypes.keys()),
        on_bad_lines="warn",
        header=0,
        low_memory=True,
        na_values={"(N/A)"},
        keep_default_na=True,
    )
    partition_cols = ["year"]
    for chunk in csv_iter:
        chunk = (
            chunk.pipe(fix_monthly_dates)
            .pipe(coerce_to_float, cols=float_cols)
            .pipe(coerce_to_integer, cols=int_cols)
            .pipe(extension_int_to_float)
            .pipe(fill_na_for_partitioning, cols=partition_cols)
        )
        table = pyarrow.Table.from_pandas(chunk, preserve_index=False, schema=pq_schema)
        pq.write_to_dataset(
            table, root_path=str(out_dir), partition_cols=partition_cols, version="2.0"
        )
    logger.info(f"  Done with monthly production {csv_file}")


def fill_na_for_partitioning(df, cols):
    """
    Replace NaN or null values with non-NA placeholders.

    These placeholders vary by column type: "NA" for string or category,
    -9 for numeric.

    The goal here is to prevent write_to_dataset from dropping groups with
    missing-valued partition columns. This function can be cut out when this
    pyarrow issue is fixed:
    https://issues.apache.org/jira/projects/ARROW/issues/ARROW-7345
    """
    orig_order = df.columns
    assert orig_order.nunique() == df.shape[1]

    to_fill = df[cols]
    fill_str = to_fill.select_dtypes(include=["object", "category"]).fillna("NA")
    fill_num = to_fill.select_dtypes(include="number").fillna(-9)
    unfilled = to_fill.columns.difference(fill_str.columns).difference(fill_num.columns)
    if len(unfilled) > 0:
        raise NotImplementedError(f"Can't fill columns: {unfilled.tolist()}")
    # Everything else:
    remainder = df.drop(cols, axis=1)
    out = pd.concat([fill_str, fill_num, remainder], axis=1)[orig_order.tolist()]
    # Reset columns to original ordering:
    return out


def coerce_cols(df, cols, target_type, drop_bad=True, **kwargs):
    """Coerce `cols` to `target_type`, optionally dropping rows that fail.

    A very small number of rows fail to parse because they have things like
    "WELL" in the oil quantity column. Just drop these.
    But because of this, we have to read all of these column as str, then convert

    args: df a dataframe
    cols: columns, currently string, that will be coerced
    target_type: what type to coerce to?
    drop_bad: should we drop the whole row if the coersion fails? (Default True)
    kwargs: arguments passed on to pd.to_datetime or pd.to_numeric.
    """

    rename_dict = {c: c + "_str" for c in cols}
    if df.columns.isin(rename_dict.values()).any():
        raise ValueError("temp rename column already exists")
    df = df.rename(columns=rename_dict)

    if target_type == "datetime":
        conversion_fn = pd.to_datetime
    else:
        conversion_fn = partial(pd.to_numeric, downcast=target_type)

    for new_col, old_col in rename_dict.items():
        try:
            # Happy case first
            df[new_col] = conversion_fn(df[old_col], errors="raise", **kwargs)
        except ValueError:
            df[new_col] = conversion_fn(df[old_col], errors="coerce", **kwargs)
            newly_missing = df[new_col].isna() & df[old_col].notna()
            newly_missing_count = newly_missing.sum()
            logger.info(
                f"  {newly_missing_count} values failed to parse as {target_type}."
                + f" Here are a few: "
                + ", ".join(df.loc[newly_missing, old_col].iloc[:6])
            )
            if drop_bad:
                # Drop rows that failed to parse (very few of these)
                df = df.loc[~newly_missing]
        finally:
            del df[old_col]
    return df


coerce_to_date = partial(coerce_cols, target_type="datetime", format="%Y-%m-%d")
coerce_to_float = partial(coerce_cols, target_type="float")
coerce_to_integer = partial(coerce_cols, target_type="integer")


def convert_headers(csv_file, out_dir):
    csv_file = Path(csv_file)
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    pd_dtypes = {
        "api": "str",  # API/UWI
        "operator_alias_legacy": "str",  # Operator Alias (Legacy)
        "operator_company_name": "str",  # Operator Company Name
        "operator_reported": "str",  # Operator (Reported)
        "operator_ticker": "str",  # Operator Ticker
        "well_name": "str",  # Well/Lease Name
        "well_number": "str",  # Well Number
        "entity_type": "str",  # Entity Type
        "county": "str",  # County/Parish
        "di_basin": "str",  # DI Basin
        "di_play": "str",  # DI Play
        "di_subplay": "str",  # DI Subplay
        "reservoir": "str",  # Reservoir
        "production_type": "str",  # Production Type
        "producing_status": "str",  # Producing Status
        "drill_type": "str",  # Drill Type
        "first_prod_date": "str",  # First Prod Date
        "last_prod_date": "str",  # Last Prod Date
        "cum_gas": "float",  # Cum Gas
        "cum_oil": "float",  # Cum Oil
        "cum_boe": "float",  # Cum BOE
        "cum_water": "float",  # Cum Water
        "cum_mmcfge": "float",  # Cum MMCFGE
        "cum_bcfge": "float",  # Cum BCFGE
        "daily_gas": "float",  # Daily Gas
        "daily_oil": "float",  # Daily Oil
        "first_month_oil": "float",  # First Month Oil
        "first_month_gas": "float",  # First Month Gas
        "first_month_water": "float",  # First Month Water
        "first_6_oil": "float",  # First 6 Oil
        "first_6_gas": "float",  # First 6 Gas
        "first_6_boe": "float",  # First 6 BOE
        "first_6_water": "float",  # First 6 Water
        "first_12_oil": "float",  # First 12 Oil
        "first_12_gas": "float",  # First 12 Gas
        "first_12_boe": "float",  # First 12 BOE
        "first_12_mmcfge": "float",  # First 12 MMCFGE
        "first_12_water": "float",  # First 12 Water
        "first_24_oil": "float",  # First 24 Oil
        "first_24_gas": "float",  # First 24 Gas
        "first_24_boe": "float",  # First 24 BOE
        "first_24_mmcfge": "float",  # First 24 MMCFGE
        "first_24_water": "float",  # First 24 Water
        "first_60_oil": "float",  # First 60 Oil
        "first_60_gas": "float",  # First 60 Gas
        "first_60_boe": "float",  # First 60 BOE
        "first_60_water": "float",  # First 60 Water
        "first_60_mmcfge": "float",  # First 60 MMCFGE
        "prac_ip_oil_daily": "float",  # Prac IP Oil Daily
        "prac_ip_gas_daily": "float",  # Prac IP Gas Daily
        "prac_ip_cfged": "float",  # Prac IP CFGED
        "prac_ip_boe": "float",  # Prac IP BOE
        "latest_oil": "float",  # Latest Oil
        "latest_gas": "float",  # Latest Gas
        "latest_water": "float",  # Latest Water
        "prior_12_oil": "float",  # Prior 12 Oil
        "prior_12_gas": "float",  # Prior 12 Gas
        "prior_12_water": "float",  # Prior 12 Water
        "last_test_date": "str",  # Last Test Date
        "last_flow_pressure": "float",  # Last Flow Pressure
        "last_whsip": "float",  # Last WHSIP
        "2nd_month_gor": "float",  # 2nd Month GOR
        "latest_gor": "float",  # Latest GOR
        "cum_gor": "float",  # Cum GOR
        "last_12_yield": "float",  # Last 12 Yield
        "2nd_month_yield": "float",  # 2nd Month Yield
        "latest_yield": "float",  # Latest Yield
        "peak_gas": "float",  # Peak Gas
        "peak_gas_month_no.": "Int32",  # Peak Gas Month No.
        "peak_oil": "float",  # Peak Oil
        "peak_oil_month_no.": "Int32",  # Peak Oil Month No.
        "peak_boe": "float",  # Peak BOE
        "peak_boe_month_no.": "Int32",  # Peak BOE Month No.
        "peak_mmcfge": "float",  # Peak MMCFGE
        "peak_mmcfge_month_no.": "Int32",  # Peak MMCFGE Month No.
        "upper_perforation": "float",  # Upper Perforation
        "lower_perforation": "float",  # Lower Perforation
        "gas_gravity": "float",  # Gas Gravity
        "oil_gravity": "float",  # Oil Gravity
        "completion_date": "str",  # Completion Date
        "well_count": "Int32",  # Well Count
        "max_active_wells": "Int32",  # Max Active Wells
        "months_produced": "Int32",  # Months Produced
        "gas_gatherer": "str",  # Gas Gatherer
        "oil_gatherer": "str",  # Oil Gatherer
        "lease_number": "str",  # Lease Number
        "spud_date": "str",  # Spud Date
        "measured_depth_td": "float",  # Measured Depth (TD)
        "true_vertical_depth": "float",  # True Vertical Depth
        "gross_perforated_interval": "float",  # Gross Perforated Interval
        "field": "str",  # Field
        "state": "str",  # State
        "district": "str",  # District
        "aapg_geologic_province": "str",  # AAPG Geologic Province
        "country": "str",  # Country
        "section": "str",  # Section
        "township": "str",  # Township
        "range": "str",  # Range
        "abstract": "str",  # Abstract
        "block": "str",  # Block
        "survey": "str",  # Survey
        "ocs_area": "str",  # OCS Area
        "pgc_area": "str",  # PGC Area
        "surface_latitude_wgs84": "float",  # Surface Latitude (WGS84)
        "surface_longitude_wgs84": "float",  # Surface Longitude (WGS84)
        "last_12_oil": "float",  # Last 12 Oil
        "last_12_gas": "float",  # Last 12 Gas
        "last_12_water": "float",  # Last 12 Water
        "entity_id": "int32",  # Entity ID
    }
    date_cols = [
        "first_prod_date",
        "last_prod_date",
        "last_test_date",
        "completion_date",
        "spud_date",
    ]

    output_dtypes = pd_dtypes.copy()
    for d in date_cols:
        output_dtypes[d] = "date32"
    output_dtypes["first_prod_year"] = "int32"
    pq_schema = dtypes_to_schema(output_dtypes)
    csv_iter = read_csv(
        csv_file,
        dtype=pd_dtypes,
        chunksize=1_000_000,
        index_col=False,
        names=list(pd_dtypes.keys()),
        on_bad_lines="warn",
        parse_dates=date_cols,
        header=0,
        low_memory=True,
    )
    partition_cols = ["first_prod_year"]
    for chunk in csv_iter:
        chunk["first_prod_year"] = chunk["first_prod_date"].dt.year
        chunk = chunk.pipe(extension_int_to_float, exclude=partition_cols).pipe(
            fill_na_for_partitioning, cols=partition_cols
        )
        table = pyarrow.Table.from_pandas(chunk, preserve_index=False, schema=pq_schema)
        pq.write_to_dataset(
            table, root_path=str(out_dir), partition_cols=partition_cols, version="2.0"
        )
    logger.info(f"  Done with headers: {csv_file}")


def make_uniform_bins(x, num):
    import numpy as np

    # Bins are evenly spaced in the range of unsigned 64bit ints
    bins = np.linspace(0, np.iinfo("uint64").max - 1024, num=num, dtype="uint64")
    # Hash (determanistic) of the values is a uint64
    hashes = pd.util.hash_pandas_object(x, index=False)
    assert hashes.dtypes == "uint64"
    return np.digitize(hashes, bins)


def regroup_parquet_fragments(directory, recursive=True):
    """
    Read all the parquet files in a directory and write a single file.
    If there are inner directories, `regroup_parquet_fragments` will be called
    recursively on those.

    Why:
    Because I'm calling write_to_dataset repeatedly for multiple chunks of data,
    I sometimes get multiple fragments for the same partition. That's a pain,
    so this function fixes that. After running it, there will be only one file
    (or zero if the directory was initially empty), named `file.parquet`.
    """
    assert directory.is_dir()
    pq_files = []
    dirs = []
    other = []
    for p in directory.iterdir():
        if p.is_file() and p.suffix == ".parquet":
            pq_files.append(p)
        elif p.is_dir():
            dirs.append(p)
        else:
            other.append(p)
    if len(other):
        raise ValueError("Unexpected files found: " + ", ".join(other))
    if recursive:
        [regroup_parquet_fragments(d, recursive=recursive) for d in dirs]

    outfile = directory.joinpath("file.parquet")
    assert not outfile.exists()
    if len(pq_files) == 0:
        return
    elif len(pq_files) == 1:
        pq_files[0].rename(outfile)
        return
    else:
        data = pq.ParquetDataset(pq_files, memory_map=True).read(
            use_pandas_metadata=True
        )
        pq.write_table(data, outfile, version="2.0")
        [f.unlink() for f in pq_files]


def main(snakemake):
    # Fill with partial so we can easily map with one argument.
    monthly = partial(convert_monthly, out_dir=Path(snakemake.output["monthly"]))
    headers = partial(convert_headers, out_dir=Path(snakemake.output["headers"]))

    with ProcessPoolExecutor(max_workers=snakemake.threads) as ex:
        ex.map(monthly, snakemake.input["monthly"])
        ex.map(headers, snakemake.input["headers"])
    regroup_parquet_fragments(Path(snakemake.output["headers"]))
    regroup_parquet_fragments(Path(snakemake.output["monthly"]))


if __name__ == "__main__":
    main(snakemake)

source(here::here("code/shared_functions.r"))

cite_attached_packages(snakemake@output[[1]], c(
  "base", # R base
  "arrow",
  "brms",
  "broom",
  "cmdstanr",
  "curl",
  "data.table",
  "digest",
  "dplyr",
  "forcats",
  "fs",
  "furrr",
  "ggplot2",
  "future",
  "glue",
  "here",
  "igraph",
  "jsonlite",
  "lubridate",
  "magrittr",
  "matrixStats",
  "nngeo",
  "posterior",
  "processx",
  "purrr",
  "RColorBrewer",
  "readxl",
  "rlang",
  "safejoin",
  "sf",
  "stringr",
  "tibble",
  "tidyr",
  "tidyselect",
  "units",
  "unix"
))

source(here::here("code/shared_functions.r"))

INPUT_CRS <- sf::st_crs(4326) # https://epsg.io/4326
WORKING_CRS <- sf::st_crs(6350) # https://epsg.io/6350


find_well_pads <- function(df, well_pad_number_start, verbose=FALSE) {
  if (verbose) {
    message("  ", paste0(sort(unique(df$aapg_geologic_province)), collapse=", "))
  }
  df_radius <- df %>%
    # Drop different entity_id at the exact same location
    # (distinct.sf already keeps geometry column by default, but for some
    # reason gives an error)
    dplyr::distinct(aapg_geologic_province, geometry) %>%
    dplyr::mutate(
      buffer_radius_m = dplyr::case_when(
        aapg_geologic_province == "SAN JOAQUIN BASIN" ~ units::as_units(20, "m"),
        TRUE ~ units::as_units(50, "m")
      ))
  well_pad_geometry <- df_radius %>%
    sf::st_geometry() %>%
    sf::st_buffer(dist=df_radius$buffer_radius_m) %>%
    st_union_intersection()
  new_well_pad_ids <- seq.int(
    from = well_pad_number_start,
    length.out = length(well_pad_geometry)
  )
  well_pad_info <- well_pad_geometry %>%
    sf::st_centroid() %>%
    sf::st_transform(crs=INPUT_CRS) %>%
    sf::st_coordinates() %>%
    dplyr::as_tibble() %>%
    dplyr::transmute(
      well_pad_lon = .data$X,
      well_pad_lat = .data$Y,
      well_pad_id  = new_well_pad_ids,
    )
  well_pad_df <- sf::st_sf(well_pad_info, geometry = well_pad_geometry)

  out <- dplyr::select(df, entity_id) %>%
    sf::st_join(well_pad_df, left=FALSE) %>%
    sf::st_drop_geometry() %>%
    dplyr::as_tibble()
  out
}


assign_well_pads <- function(header_df, verbose=TRUE) {
  stopifnot(!anyDuplicated(header_df$entity_id))
  # My definition of a well pad, based on Alvarez et al. (2018) SI:
  # - a group of wells that are close together
  # - if the DI record is a lease (irrelevant for CA, NM, and post-1999 CO),
  #   assume it is its own well pad; don't group it with other records
  # - if more than two wells are nearby, take the union. (See the function
  #   st_union_intersection.)
  # NB: It would be nice to use operator, but there are a bunch of wells with
  # different operators reported at the same location.
  # Steps:
  # 1. Make into an sf dataframe
  # 2. Transform to a useful projection
  # 3. Buffer all points 50m, except the San Joaquin basin, where the buffer is 20m.
  # 4. Union the buffers and assign IDs
  # 5. Take the intersection of the original points with the buffers
  # 6. Do the above separately for each state or basin, then re-aggregate.
  # Note that the well-pad definitions are sensitive to the projection and
  # assumed well spacing (eg. 50m) for a small number of wells.

  # For records that are marked as "LEASE" (as opposed to "COM", "DRIP POINT",
  # "SWD", "UNIT", or "WELL"), don't do the geographic work, just label them
  # their own singleton well pad.
  singletons <- header_df %>%
    dplyr::group_by(aapg_geologic_province, operator_company_name) %>%
    dplyr::filter(dplyr::n() == 1 | entity_type == "LEASE") %>%
    dplyr::ungroup() %>%
    dplyr::select(entity_id, surface_longitude_wgs84, surface_latitude_wgs84) %>%
    dplyr::mutate(well_pad_id = dplyr::row_number()) %>%
    dplyr::rename(
      well_pad_lon = surface_longitude_wgs84,
      well_pad_lat = surface_latitude_wgs84,
    )

  # For records that aren't LEASE or singleton, look for nearby wells within
  # groups defined by basin and operator_company_name. Note that there are a
  # small number of overlaps, where a well falls within a well pad with
  # different basin or operator_company_name. (This is true even if you only use
  # basin.)
  non_singleton <- header_df %>%
    dplyr::anti_join(singletons, by="entity_id") %>%
    sf::st_as_sf(
      crs=INPUT_CRS,
      coords=c("surface_longitude_wgs84", "surface_latitude_wgs84")
    ) %>%
    sf::st_transform(crs=WORKING_CRS)

  # well_pads is the geometry of the well pads (a POLYGON buffered around the
  # points of each group of wells)
  non_singleton_lst <- non_singleton %>%
    # Split by basin so we can run different regions in parallel
    # Well pads will not span different subgroups of the group_by, so be careful
    # before adding grouping variables.
    dplyr::group_by(aapg_geologic_province, operator_company_name) %>%
    dplyr::group_split()
  # We then also need to make sure we're not assigning duplicate well_pad_id
  # values in the different parallel processes, so calculate disjoint sets of
  # possible_well_pad_ids
  n_groups <- length(non_singleton_lst)
  subgroup_row_counts <- purrr::map_int(non_singleton_lst, nrow)
  well_pad_id_bounds <- c(nrow(singletons), subgroup_row_counts) %>% cumsum()
  well_pad_id_starts <- dplyr::lag(well_pad_id_bounds)[-1] + 1
  # possible_well_pad_ids isn't used, just checked
  possible_well_pad_ids <- purrr::map2(
      well_pad_id_starts,
      well_pad_id_bounds[-1],
      ~seq.int(from=.x, to=.y, by=1)
    ) %>%
    unlist() %>%
    c(singletons$well_pad_id)

  stopifnot(
    all(subgroup_row_counts) > 0,
    anyDuplicated(possible_well_pad_ids) == 0,
    length(non_singleton_lst) == length(well_pad_id_starts)
  )

  # Actually do the well pad creation, and re-add the singletons
  entity_pad_crosswalk <- furrr::future_map2_dfr(
      non_singleton_lst, well_pad_id_starts, find_well_pads,
      verbose=verbose,
      # we're not doing any RNG, but suppress warnings
      .options=furrr::furrr_options(seed=TRUE)
    ) %>%
    dplyr::bind_rows(singletons)
  # Note: singletons got assigned well_pad_id values 1:nrow(singletons), but it
  # doesn't matter that they're bound at the end here.

  out <- dplyr::inner_join(header_df, entity_pad_crosswalk, by="entity_id")
  stopifnot(
    anyDuplicated(header_df$entity_id) == 0,
    anyDuplicated(entity_pad_crosswalk$entity_id) == 0,
    setequal(header_df$entity_id, out$entity_id)
  )
  if (verbose) {
    n_wells <- nrow(header_df)
    n_unmatched <- n_wells - nrow(out)
    well_pad_well_counts <- dplyr::count(out, well_pad_id)
    n_wells_at_multi_well_pads <- dplyr::filter(well_pad_well_counts, n > 1) %>%
      dplyr::pull("n") %>%
      sum_() %>%
      fill_na(0)
    n_wells_at_single_well_pads <- dplyr::filter(well_pad_well_counts, n == 1) %>%
      dplyr::pull("n") %>%
      sum_() %>%
      fill_na(0)
    glue_message(
      "Of {n_wells} total wells, {n_wells_at_multi_well_pads} were matched to ",
      "multi-well pads, {n_wells_at_single_well_pads} were matched to single-",
      "well pads, and {n_unmatched} were unmatched for data quality reasons."
    )
  }
  out
}


parse_year_range <- function(year_range) {
  # Use a regex header_df to pull out the years
  # look for things like "1990-2018"
  stopifnot(length(year_range) == 1)
  years <- stringr::str_match(year_range, "^(\\d{4})-(\\d{4})$")[1, 2:3] %>%
    as.integer() %>%
    unclass()
  stopifnot(length(years) == 2, !anyNA(years))
  years
}


create_well_pad_crosswalk <- function(header_dir, output_file, year_range) {
  year_range %<>% parse_year_range()
  stopifnot(dir.exists(header_dir), length(output_file) == 1)
  header_df <- arrow::open_dataset(header_dir) %>%
    dplyr::filter(
      first_prod_year >= year_range[1],
      first_prod_year <= year_range[2],
      !is.na(aapg_geologic_province),
      aapg_geologic_province != "(N/A)",
      !is.na(surface_latitude_wgs84),
      !is.na(surface_longitude_wgs84),
    ) %>%
    dplyr::select(
      entity_id, aapg_geologic_province,
      surface_latitude_wgs84, surface_longitude_wgs84,
      # NA operator_company_name is allowed (currently coded as "(N/A)")
      operator_company_name, entity_type
    ) %>%
    dplyr::collect() %>%
    dplyr::distinct()  # distinct for multiple wells at the same entity and location

  well_pad_crosswalk <- assign_well_pads(header_df, verbose=FALSE)
  arrow::write_parquet(well_pad_crosswalk, output_file)
}


test_well_pad_creation <- function() {
  header_df <- tibble::tibble(
    entity_id = 1:4,
    aapg_geologic_province = "SAN JOAQUIN BASIN",
    entity_type = "WELL",
    operator_company_name="AAAAA",
    # First two points with identical location, third point about 9m away,
    # and fourth point 4.3km away
    surface_longitude_wgs84 = c(-119.7514, -119.7514, -119.7515, -119.8),
    surface_latitude_wgs84 = 35.48126,
  )

  expected_results <- header_df %>%
    dplyr::mutate(well_pad_id = c(1, 1, 1, 2))
  actual_results <- assign_well_pads(header_df, verbose=FALSE) %>%
    dplyr::select(-well_pad_lon, -well_pad_lat) # don't test that part
  compare_results <- all.equal(actual_results, expected_results)
  if (!isTRUE(compare_results)) {
    print(tibble::tibble(
      expected_id = expected_results$well_pad_id,
      actual_id = actual_results$well_pad_id,
    ))
    print(actual_results)
    stop(compare_results)
  }
}


if (!exists("snakemake")) {
  stop("This script is meant to be run with snakemake.")
}

# Set up resources (mem limit doesn't work on MacOS)
memory_limit(snakemake@resources[['mem_mb']])
future::plan("multicore", workers = snakemake@threads)

test_well_pad_creation()

create_well_pad_crosswalk(
  header_dir = snakemake@input$headers,
  output_file = snakemake@output[[1]],
  year_range = snakemake@wildcards$year_range
)

source(here::here("code/shared_functions.r"))

options(scipen=5, SOURCE_DATE_EPOCH=0)
set.seed(6350) # CRS as seed? Sure.

# Note on variable names:
# mcfd is mcf per day (thousands of cubic feet per day), but is reported monthly
# bbld is barrels per day, but is measured monthly.
# There are other variables that report the average on days when production was
# happening, but I don't use those (maybe I should?).

LON_LAT_CRS <- sf::st_crs(4326) # https://epsg.io/4326
OUTPUT_CRS <- sf::st_crs(6350) # https://epsg.io/6350 (CONUS Albers)
MAX_ACCEPTABLE_MATCH_DIST_METERS <- 500
# These are the source types we'll keep. There are others, like landfills and
# dairies, that we don't care about. This is a list of oil and gas types.
# I wasn't sure about gas compressors, but they seem to be far from wells.
SOURCE_TYPES <- c(
  # "gas compressor",
  # "gas distribution line",
  # "gas LNG station",
  # "gas processing plant",
  # "gas storage facility",
  "oil/gas compressor",
  "oil/gas drill rig",
  "oil/gas gathering line",
  "oil/gas possible plugged well",
  "oil/gas pumpjack",
  "oil/gas stack",
  "oil/gas tank",
  "oil/gas unknown infrastucture",
  "oil/gas waste lagoon"
)

read_jpl_plumes <- function(input_files) {
  df <- data.table::fread(input_files[["duren_2019_plumes"]], data.table=FALSE) %>%
    dplyr::as_tibble() %>%
    dplyr::rename_all(make_better_names) %>%
    dplyr::filter(source_type_best_estimate %in% !!SOURCE_TYPES) %>%
    dplyr::mutate(
      datetime_of_detection = lubridate::mdy_hm(paste(date_of_detection, time_of_detection_utc), tz="UTC"),
      date_of_detection = lubridate::mdy(date_of_detection),
      flight_name = stringr::str_match(candidate_identifier, "^(ang\\d{8}t\\d{6})")[, 2L, drop=TRUE]
    ) %>%
    dplyr::rename(
      emiss_kg_hr = qplume_kg_hr_plume_emissions,
      emiss_se_kg_hr = sigma_qplume_kg_hr_uncertainty_for_plume_emissions,
    ) %>%
    dplyr::select(-sectors_ipcc, -time_of_detection_utc) %>%
    ensure_id_vars(source_identifier, datetime_of_detection) %>%
    sf::st_as_sf(
      crs=LON_LAT_CRS,
      coords=c("plume_longitude_deg", "plume_latitude_deg"),
      agr=c(
        source_identifier = "identity", candidate_identifier = "identity",
        date_of_detection = "constant", source_type_best_estimate = "constant",
        emiss_kg_hr = "constant", emiss_se_kg_hr = "constant",
        datetime_of_detection = "constant", flight_name = "constant"
      )
    ) %>%
    sf::st_transform(OUTPUT_CRS)
  stopifnot(!anyNA(df$flight_name))
  df
}

read_jpl_sites <- function(input_files) {
  # Sites are different than plumes because they revisited
  # This excel file is translated from a PDF in the Durent et al. (2019)
  # supplementary materials.
  df <- readxl::read_excel(input_files[["duren_2019_sites"]]) %>%
    dplyr::rename_all(make_better_names) %>%
    dplyr::filter(source_type %in% !!SOURCE_TYPES) %>%
    dplyr::rename(
      emiss_kg_hr = qsource_kg_hr,
      persistence_frac = source_persistence_f,
    ) %>%
    # Run this rename twice because the sigma character gets assigned a different name on windows.
    rename_cols(c("emiss_se_kg_hr" = "x_q_kg_hr"), strict=FALSE) %>%
    rename_cols(c("emiss_se_kg_hr" = "x_u_f073_q_kg_hr"), strict=FALSE) %>%
    dplyr::select(-ipcc_sector, -confidence_in_persistence) %>%
    ensure_id_vars(source_identifier) %>%
    sf::st_as_sf(
      crs=LON_LAT_CRS,
      coords=c("source_longitude_deg", "source_latitude_deg"),
      agr=c(
        source_identifier = "identity", source_type = "constant",
        n_overflights = "constant", persistence_frac = "constant",
        emiss_kg_hr = "constant", emiss_se_kg_hr = "constant"
      )
    ) %>%
    sf::st_transform(OUTPUT_CRS)
  df
}

read_headers <- function(years, states) {
  header_dir <- here::here("data/generated/production/well_headers/")
  stopifnot(dir.exists(header_dir), !anyNA(years), !anyNA(states))
  date_min <- lubridate::make_date(min(years), 1, 1)
  date_max <- lubridate::make_date(max(years) + 1, 1, 1)

  well_headers <- arrow::open_dataset(header_dir) %>%
    dplyr::filter(
      state %in% !!states,
      # !(is.na(first_60_oil) & is.na(first_60_gas)),
      # !is.na(completion_date)
    ) %>%
    dplyr::select(
      county, state, production_type, drill_type, aapg_geologic_province,
      surface_latitude_wgs84, surface_longitude_wgs84, first_prod_date, last_prod_date,
      first_60_oil, first_60_gas, completion_date, spud_date, months_produced, entity_id,
    ) %>%
    dplyr::collect() %>%
    dplyr::filter(
      !is.na(surface_latitude_wgs84),
      !is.na(surface_longitude_wgs84),
      !is.na(first_prod_date),
      !is.na(last_prod_date),
      # Drop wells that started producing after the end of observation or ended
      # before the beginning (just for efficiency of not doing spatial operations
      # on a bunch of irrelevant wells)
      first_prod_date < !!date_max,
      last_prod_date > !!date_min,
      !production_type %in% c("WATER", "STEAMFLOOD"),
    ) %>%
    dplyr::rename(basin = aapg_geologic_province) %>%
    dplyr::mutate(production_type = dplyr::case_when(
      production_type %in% c("GAS", "OIL") ~ production_type,
      production_type == "GAS STORE" ~ "GAS",
      production_type == "OIL (CYCLIC STEAM)" ~ "OIL",
      TRUE ~ NA_character_
    )) %>%
    ensure_id_vars(entity_id) %>%
    sf::st_as_sf(
      crs=LON_LAT_CRS,
      coords=c("surface_longitude_wgs84", "surface_latitude_wgs84"),
      agr=c(
        entity_id = "identity", county = "constant", production_type = "constant",
        drill_type = "constant", basin = "constant", months_produced = "constant",
        surface_latitude_wgs84 = "constant", surface_longitude_wgs84 = "constant",
        first_prod_date = "constant", last_prod_date = "constant",
        first_60_oil = "constant", first_60_gas = "constant",
        completion_date = "constant", spud_date = "constant"
      )
    ) %>%
    sf::st_transform(OUTPUT_CRS)
  well_headers
}

match_wells_to_flight_day <- function(well_headers, flight_paths) {
  wells_flown_over <- well_headers %>%
    ensure_id_vars(entity_id) %>%
    # Only consider wells that are in the flight paths *and* are producing when
    # the plane flys over. This join generates duplicate rows wells because
    # the there are multiple flights over the same areas. Keep rows that had
    # any overpass while active, then drop duplicates.
    # (doing it this way because distinct doesn't work quite the same on sf objects)
    sf::st_join(dplyr::select(flight_paths, flight_date), join=sf::st_intersects, left=FALSE) %>%
    dplyr::filter(flight_date >= first_prod_date, flight_date <= last_prod_date) %>%
    sf::st_drop_geometry() %>%
    dplyr::distinct(entity_id, flight_date)

  dplyr::inner_join(well_headers, wells_flown_over, by="entity_id")
}

st_nn <- function(...){
  # Suppress nngeo's narration
  suppressMessages(nngeo::st_nn(...))
}

load_flight_paths <- function(plume_measurements) {
  # Flight paths:
  flight_df <- data.table::fread(
      here::here("data/studies/duren_etal_2019/AVIRIS-NG Flight Lines - AVIRIS-NG Flight Lines.csv"),
      data.table=FALSE
    ) %>%
    dplyr::as_tibble() %>%
    dplyr::rename(flight_name = Name) %>%
    dplyr::rename_all(make_better_names) %>%
    dplyr::mutate(
      flight_date = lubridate::make_date(.data$year, .data$month, .data$day),
      utc_hour = dplyr::if_else(dplyr::between(utc_hour, 0, 23), utc_hour, NA_integer_),
      utc_minute = dplyr::if_else(dplyr::between(utc_minute, 0, 23), utc_minute, NA_integer_),
      flight_datetime_utc = lubridate::make_datetime(
        .data$year, .data$month, .data$day, .data$utc_hour, .data$utc_minute
      ),
    ) %>%
    dplyr::semi_join(plume_measurements, by="flight_name") %>% # keep matching flights
    dplyr::select(flight_name, flight_date, flight_datetime_utc, number_of_lines, investigator,
      dplyr::starts_with("lat"), dplyr::starts_with("lon")
    ) %>%
    dplyr::distinct() %>%
    ensure_id_vars(flight_name)
  stopifnot(setequal(flight_df$flight_name, plume_measurements$flight_name))

  flight_polygons <- flight_df %>%
    dplyr::select(flight_name, dplyr::starts_with("lon"), dplyr::starts_with("lat")) %>%
    # This is terrible.
    tidyr::pivot_longer(dplyr::starts_with("lon"), names_to="idx1", names_prefix="lon", values_to="lon") %>%
    tidyr::pivot_longer(dplyr::starts_with("lat"), names_to="idx2", names_prefix="lat", values_to="lat") %>%
    dplyr::filter(idx1 == idx2) %>%
    dplyr::filter(!is.na(lat), !is.na(lon)) %>%
    sf::st_as_sf(coords = c("lon", "lat"), crs = LON_LAT_CRS) %>%
    dplyr::group_by(flight_name) %>%
    dplyr::summarize(
      geometry = sf::st_cast(sf::st_combine(geometry), "POLYGON"),
      .groups="drop"
    ) %>%
    sf::st_transform(OUTPUT_CRS)

  # Join back the variables (flight_polygons has to be first to get the s3 dispatch right)
  out <- dplyr::inner_join(flight_polygons,
    dplyr::select(flight_df, flight_name, flight_date, flight_datetime_utc, number_of_lines, investigator),
    by="flight_name"
  )
  if ("date_of_detection" %in% colnames(plume_measurements)) {
    # flight dates are detection dates:
    stopifnot(setequal(out$flight_date, plume_measurements$date_of_detection))
  }
  stopifnot(
    # detected plumes are inside flight paths:
    all(sf::st_intersects(plume_measurements, sf::st_union(out), sparse=FALSE))
  )
  out
}

load_monthly_prod <- function(wells_by_flight_day) {
  if (inherits(wells_by_flight_day, "sf")) {
    wells_by_flight_day %<>% sf::st_drop_geometry()
  }
  prod_dir <- here::here("data/generated/production/monthly_production/")
  stopifnot(
    dir.exists(prod_dir),
    is_id(wells_by_flight_day, entity_id, flight_date)
  )
  well_to_match <- wells_by_flight_day %>%
    dplyr::transmute(
      flight_date = flight_date,
      year = lubridate::year(flight_date),
      month = lubridate::month(flight_date),
      entity_id = entity_id
  )
  # These are just for speed, so we don't read anything we're positive we don't want
  desired_entity_id <- unique(well_to_match$entity_id)
  desired_year <- unique(well_to_match$year)
  desired_month <- unique(well_to_match$month)

  # Only diff from the stored schema is year, which currently registers as string
  # To select another column, add it to this schema.
  schema <- arrow::schema(
    year = arrow::int32(), month = arrow::int32(), entity_id = arrow::int32(),
    daily_avg_gas = arrow::float(), daily_avg_oil = arrow::float(), well_count = arrow::int32()
  )
  prod <- arrow::open_dataset(prod_dir, schema=schema) %>%
    dplyr::filter(
      daily_avg_gas > 0,
      year %in% !!desired_year,
      month %in% !!desired_month,
      entity_id %in% !!desired_entity_id,
    ) %>%
    dplyr::collect() %>%
    # 1:m join (U), no column conflicts (C)
    safejoin::safe_inner_join(well_to_match, by=c("entity_id", "year", "month"), check="U C")
  prod
}

load_well_records <- function(flight_paths, states, nat_gas_price) {
  years <- flight_paths$flight_date %>% lubridate::year() %>% unique()
  # First, find wells that were flown over and had production start before the
  # flight date and end after the flight date.
  # wells_by_flight_day has one row for each well each day it was flown over.
  wells_by_flight_day <- read_headers(years, states) %>%
    match_wells_to_flight_day(flight_paths)

  # Note: If the well was flown over in different months, this is the average
  # across months.
  monthly_prod_when_flown_over <- load_monthly_prod(wells_by_flight_day) %>%
    dplyr::group_by(entity_id) %>%
    dplyr::summarize(
      oil_avg_bbld = mean_(daily_avg_oil),
      gas_avg_mcfd = mean_(daily_avg_gas),
      .groups="drop"
    )
  # Recall wells_by_flight_day has one row for each well each day it was flown over.
  # Drop down to one row per well.
  # Do both well age and price here because it's easier than handling prices later
  well_age_and_price <- wells_by_flight_day %>%
    sf::st_drop_geometry() %>%
    dplyr::mutate(age_yr = as.numeric(flight_date - first_prod_date) / 365.25) %>%
    match_full_state_names() %>%
    dplyr::rename(state = state_full) %>%
    harmonize_basin_name() %>%
    match_commodity_prices(nat_gas_price) %>%
    dplyr::group_by(entity_id) %>%
    dplyr::summarize(
      age_yr = mean(age_yr),
      gas_price_per_mcf = mean(gas_price_per_mcf),
      gas_frac_methane = mean(gas_frac_methane),
      .groups="drop"
    )
  stopifnot(noNAs(well_age_and_price))

  wells <- wells_by_flight_day %>%
    dplyr::group_by(entity_id) %>%
    dplyr::select(-flight_date) %>%
    dplyr::distinct(.keep_all=TRUE) %>% # could speed this up here by not doing a spatial distinct
    ensure_id_vars(entity_id) %>%
    dplyr::inner_join(well_age_and_price, by="entity_id") %>% # 1:1 join
    harmonize_basin_name(group_small_CA_basins=TRUE)
  dplyr::inner_join(wells, monthly_prod_when_flown_over, by="entity_id") # 1:1 join
}

match_with_wells <- function(observed_sites, wells) {
  sites_matched <- sf::st_join(
      observed_sites, wells,
      join=st_nn, k=1, maxdist=MAX_ACCEPTABLE_MATCH_DIST_METERS,
      left=FALSE,
      progress=FALSE
    ) %>%
    # drop doubly-matched sites.
    dplyr::group_by(entity_id) %>%
    dplyr::filter(dplyr::n() == 1) %>%
    dplyr::ungroup() %>%
    ensure_id_vars(entity_id)
  if (nrow(sites_matched) != nrow(observed_sites)) {
    warning(
      nrow(observed_sites) - nrow(sites_matched),
      " leaks didn't have a matching well"
    )
  }
  sites_matched
}

make_pairs_plots <- function(jpl_sites_matched) {
  if (!rlang::is_installed("GGally")) {
    stop("GGally package required for this function")
  }
  pseudo_log <- scales::pseudo_log_trans() # basically asinh

  plt_matched <- jpl_sites_matched %>%
    dplyr::mutate_at(dplyr::vars(oil_avg_bbld, gas_avg_mcfd, emiss_kg_hr, emiss_se_kg_hr), list(log=pseudo_log$transform)) %>%
    sf::st_drop_geometry() %>%
    GGally::ggpairs(
      ggplot2::aes(color=production_type, alpha=0.7),
      columns = c("emiss_kg_hr_log", "oil_avg_bbld_log", "gas_avg_mcfd_log", "persistence_frac"),
      progress=FALSE
    ) +
    ggplot2::theme_bw()
  save_plot(plt_matched, here::here("graphics/pairs_plot_matched_wells_jpl_di.pdf"), scale=3)

}

make_plots <- function(wells_all, jpl_sites_matched, ground_studies) {
  # GAS_SCALE <- ggplot2::scale_x_continuous(trans="pseudo_log", breaks=c(0, 10, 100, 1000, 10000, 100000))
  # EMISS_SCALE <- ggplot2::scale_y_continuous(trans="pseudo_log", breaks=c(0, 10, 30, 100, 1000))

  make_plot_qq <- function(df, add_line=TRUE, ...) {
    df <- dplyr::filter(df, emiss_kg_hr > 0)

    params <- as.list(MASS::fitdistr(df$emiss_kg_hr, "log-normal")$estimate)
    # The `...` here is just so I can pass in a color aesthetic in the combined case
    plt <- ggplot2::ggplot(df, ggplot2::aes(sample=emiss_kg_hr, ...)) +
      ggplot2::geom_qq(distribution=stats::qlnorm, dparams=params) +
      ggplot2::scale_x_continuous(
        trans="pseudo_log",
        limits=c(-0.2, 3000),
        breaks=c(0, 10, 30, 100, 1000)
      ) +
      ggplot2::scale_y_continuous(
        trans="pseudo_log",
        limits=c(-0.2, 3000),
        breaks=c(0, 10, 30, 100, 1000)
      ) +
      ggplot2::theme_bw() +
      ggplot2::labs(
        x="Theoretical distribution (kg/hr)",
        y="Observed distribution (kg/hr)"
      )
    if (add_line) {
      plt <- plt + ggplot2::geom_qq_line(distribution=stats::qlnorm, dparams=params)
    }
    plt
  }
  make_plot_emiss_prod_point_density <- function(df, adjust=0.1) {
    min_emiss <- min_(dplyr::filter(df, emiss_kg_hr > 0)$emiss_kg_hr)
    plt <- df %>%
      dplyr::mutate(emiss_kg_hr_filled =
        dplyr::if_else(is.na(emiss_kg_hr), min_emiss - (0.1 * min_emiss), emiss_kg_hr)) %>%
      dplyr::filter(gas_avg_mcfd > 0) %>%
      ggplot2::ggplot(ggplot2::aes(x=gas_avg_mcfd, y=emiss_kg_hr_filled)) +
      ggpointdensity::geom_pointdensity(adjust=adjust, alpha=0.8) +
      # options are "magma", "plasma", "viridis", or "cividis"
      ggplot2::scale_color_viridis_c(option="inferno")+
      ggplot2::geom_hline(yintercept=min_emiss) +
      ggplot2::scale_x_continuous(
        trans="pseudo_log",
        breaks=c(0, 10^(1:6)),
        limits=c(0, 1.5e6)
      ) +
      ggplot2::scale_y_continuous(trans="pseudo_log", breaks=c(0, 10, 30, 100, 1000)) +
      ggplot2::theme_bw() +
      ggplot2::theme(legend.position="none") +
      ggplot2::labs(
        x="Average gas production (mcf/mo)",
        y="Measured emissions (kg/hr)"#,
        # color="Well\ncount"
      )
    plt
  }
  #
  # plt_ecdf_jpl <- ggplot2::ggplot(jpl_sites_matched, ggplot2::aes(x=emiss_kg_hr)) +
  #   ggplot2::stat_ecdf(geom="step") +
  #   ggplot2::scale_x_log10() +
  #   ggplot2::theme_bw()

  plt_obs_count_jpl <- (make_plot_emiss_prod_point_density(wells_all) +
    ggplot2::labs(
      title="Most wells have no JPL-observed emissisons"
    )) %>%
    save_plot(here::here("graphics/observation_count_jpl_flights.pdf"))

  plt_qq_jpl <- (make_plot_qq(jpl_sites_matched) +
    ggplot2::labs(
      title="Observed JPL measures are distributed log-normal"
    ))
    save_plot(plt_qq_jpl, here::here("graphics/jpl_flights_qq_plot.pdf"))

  ground_studies <- dplyr::rename(ground_studies,
    gas_avg_mcfd = gas_production_mcfd
  )
  plt_obs_count_ground <- (make_plot_emiss_prod_point_density(ground_studies, adjust=0.05) +
    ggplot2::labs(
      title="Ground measures have many fewer zeros"
    )) %>%
    save_plot(here::here("graphics/observation_count_ground_studies.pdf"))

  plt_qq_ground <- (make_plot_qq(ground_studies) +
    ggplot2::labs(
      title="Ground-based measurements"
    )) %>%
    save_plot(here::here("graphics/ground_studies_qq_plot.pdf"))

  to_plot_combined_qq <- dplyr::bind_rows(
      dplyr::mutate(ground_studies, src = "Ground studies"),
      sf::st_drop_geometry(dplyr::mutate(jpl_sites_matched, src = "JPL flights")),
    )
  plot_qq_combined <- make_plot_qq(to_plot_combined_qq, add_line = TRUE, color=src) +
    # ggplot2::facet_grid(rows="src") +
    ggplot2::labs(
      title="Comparison QQ plots"
    )
  save_plot(plot_qq_combined, here::here("graphics/combined_ground_jpl_qq_plot.pdf"))

}

aggregate_to_well_pads <- function(wells_all, well_pad_mapping_file) {
  stopifnot(
    anyDuplicated(wells_all$entity_id) == 0,
    length(well_pad_mapping_file) == 1
  )
  well_pad_mapping <- arrow::read_parquet(
    well_pad_mapping_file,
    col_select=c("entity_id", "well_pad_id", "well_pad_lon", "well_pad_lat")
  )
  # We drop study-specific columns here, but might want to bring them back.
  # n_overflights, persistence_frac, ...
  # NOTE: we're doing an inner join with well_pad_id here. That means wells that
  # aren't part of the well pad mapping will not be part of the output.
  well_pad_df <- wells_all %>%
    # Check that the join is unique by right side (V), types match (T),
    # columns don't collide (C), and all rows in left side are matched (m)
    # Note: not currently true that all LHS are matched -- m instead of M warns.
    safejoin::safe_inner_join(well_pad_mapping, by="entity_id", check="V m T C") %>%
    dplyr::group_by(well_pad_id) %>%
    dplyr::summarize(
      # These aggregating functions are defined in shared_functions.r
      emiss_kg_hr = mean_(emiss_kg_hr),
      emiss_se_kg_hr = mean_(emiss_se_kg_hr),
      production_type = Mode(production_type),
      drill_type = Mode(drill_type),
      first_60_oil = sum_(first_60_oil),
      first_60_gas = sum_(first_60_gas),
      oil_avg_bbld = sum_(oil_avg_bbld),
      gas_avg_mcfd = sum_(gas_avg_mcfd),
      months_produced = mean_(months_produced),
      county = Mode(county),
      # Like Lyon et al. 2016, we'll define pad age by age of the most recently
      # drilled well.
      age_yr = min(age_yr),
      gas_price_per_mcf = mean(gas_price_per_mcf), # probably unique by well pad
      # already unique by well_pad_id; the particular aggregation doesn't matter
      well_pad_lon = min(well_pad_lon),
      well_pad_lat = min(well_pad_lat),
      basin = min(basin),
      gas_frac_methane = min(gas_frac_methane),
      .groups="drop"
    )

  stopifnot(!anyNA(well_pad_df$well_pad_id))

  # NOTES on well pads:
  # - There's still a vast number of well pads with no measurements
  # (currently 112 with measurements and 17006 without. This could be improved
  # slightly by thinking about deduplication, but it's minimal.)
  well_pad_df
}

load_ground_studies <- function(input_files) {
  # Drop studies that can't be used for our application. See the notes in
  # notes/measurement_paper_notes.csv
  studies_that_dont_work_for_us <- c("Rella et al.")

  df_alvarez <- readxl::read_excel(input_files[["alvarez_2018"]], sheet="inputs_meas_sites") %>%
    dplyr::rename_all(make_better_names) %>%
    dplyr::select(basin, study, methane_emissions_kgh, gas_production_mcfd) %>%
    dplyr::rename(emiss_kg_hr = methane_emissions_kgh) %>%
    dplyr::mutate(
      study = dplyr::case_when(
        study == "Omara et al." ~ "Omara et al. (2016)",
        study == "Robertson et al." ~ "Robertson et al. (2017)",
        TRUE ~ study
      ),
      basin = dplyr::case_when(
        basin == "SWPA" ~ "SW Pennsylvania",
        basin == "Weld County" ~ "Denver Julesburg",
        TRUE ~ basin
      )
    )

  df_omara_2018 <- data.table::fread(input_files[["omara_2018"]], data.table=FALSE) %>%
    dplyr::as_tibble() %>%
    dplyr::transmute(
      basin = dplyr::case_when(
        site == "Uinta Basin (Uintah County, UT)" ~ "Uinta", # match with Alvarez et al.
        site == "Denver Julesburg Basin (Weld County, CO)" ~ "Denver Julesburg",
        site == "NE PA (Bradford, Susquehanna, Wyoming, Sullivan Counties)" ~ "NE Pennsylvania",
        TRUE ~ site
      ),
      study = "Omara et al. (2018)",
      gas_production_mcfd = tot_gas_prod_mcfd,
      emiss_kg_hr = avg_emissions_kg_per_h,
    )

  df_zavala_araiza_2018 <- readxl::read_excel(input_files[["zavala_araiza_2018"]]) %>%
    dplyr::rename_all(make_better_names) %>%
    dplyr::filter(method == "tracer flux") %>% # different sampling strategies
    dplyr::transmute(
      emiss_kg_hr = ch4_emission_rate_kgh,
      gas_production_mcfd = gas_production_mcfd,
      study = "Zavala-Araiza et al. (2018)",
      basin = "Alberta",
      # other variables that could be interesting:
      # ch4_lb, ch4_ub, oil_production_bbld, age_yr, wells_per_site,
      # reported_emissions_kgh, gas_composition_c1_percent
    )
  out <- dplyr::bind_rows(df_alvarez, df_omara_2018, df_zavala_araiza_2018)
  stopifnot(all(studies_that_dont_work_for_us %in% out$study)) # guard against typos
  out <- dplyr::filter(out, !study %in% !!studies_that_dont_work_for_us)
  out
}

compare_ground_studies_with_jpl <- function(wells_all, ground_studies) {
  # Table:
  # - N non-zero
  # - N zero
  # - mean, median leak rate
  # - correlation with size, among non-zero leaks
  # - age?
  wells_num_zero <- dplyr::filter(wells_all, is.na(emiss_kg_hr) | emiss_kg_hr < 5) %>% nrow()
  wells_all_stats <- dplyr::filter(wells_all, emiss_kg_hr >= 5, !is.na(gas_avg_mcfd)) %>%
    dplyr::summarize(
      source = "JPL flights",
      n_positive = dplyr::n(),
      mean_if_positive = mean_(emiss_kg_hr),
      corr_with_size = stats::cor(emiss_kg_hr, gas_avg_mcfd),
      .groups="drop"
    ) %>%
    dplyr::mutate(n_zero = !!wells_num_zero)
  ground_studies_num_zero <- dplyr::filter(ground_studies, is.na(emiss_kg_hr) | emiss_kg_hr == 0) %>% nrow()
  ground_studies_stats <- dplyr::filter(ground_studies, emiss_kg_hr > 0) %>%
    dplyr::summarize(
      source = "Ground studies",
      n_positive = dplyr::n(),
      mean_if_positive = mean_(emiss_kg_hr),
      # corr_with_size = stats::cor(emiss_kg_hr, gas_avg_mcfd)
      corr_with_size = 0,
      .groups="drop"
    ) %>%
    dplyr::mutate(n_zero = !!ground_studies_num_zero)
  tab <- cbind(t(wells_all_stats), t(ground_studies_stats))
  return(tab)
}

read_lyon_etal_2016 <- function(input_files) {
  locations <- readxl::read_excel(input_files[["lyon_etal_2016_locations"]]) %>%
    dplyr::as_tibble() %>%
    dplyr::rename_all(make_better_names) %>%
    ensure_id_vars(pad_id) %>%
    dplyr::select(-basin) %>%
    dplyr::mutate(longitude = as.double(stringr::str_trim(longitude))) # excel issues

  df <- readxl::read_excel(input_files[["lyon_etal_2016_measures"]], sheet=1) %>%
    dplyr::as_tibble() %>%
    dplyr::rename_all(make_better_names) %>%
    dplyr::select(-video_id) %>%
    ensure_id_vars(pad_id) %>%
    dplyr::left_join(locations, by="pad_id") %>%
    # Rename to harmonize with other datasets
    dplyr::rename(
      gas_avg_mcfd = gas_production_mcf_pad_day,
      oil_avg_bbld = oil_production_bbl_pad_day,
      detect_emiss = emissions_detected_0_no_1_yes,
    ) %>%
    dplyr::mutate(
      # This is approx the same as our definition
      age_yr = well_age_months_since_initial_production_of_newest_well / 12,
      flight_date = as.Date("2014-08-01"),
      state = basin_to_state(basin),
    )

  df
}

basin_to_state <- function(basin) {
  stopifnot(is.character(basin))
  # These aren't perfect, but for basins that cross states, we pick the majority
  # state. (But if you're worried about this, you should worry about using
  # state citygate prices instead.)
  conversions = c(
    "Bakken" = "North Dakota",
    "Barnett" = "Texas",
    "EagleFord" = "Texas",
    "Fayetteville" = "Arkansas",
    "Marcellus" = "Pennsylvania",
    "PowderRiver" = "Wyoming",
    "Uintah" = "Utah",
    "San Juan" = "New Mexico",
    "San Joaquin" = "California"
  )
  missing_basins <- setdiff(unique(basin), names(conversions))
  if (length(missing_basins) > 0) {
    stop("Missing states for these basins: ", paste(missing_basins, collapse=", "))
  }
  conversions[basin]
}

read_frankenberg_etal_2016 <- function(input_files) {
  measures <- readxl::read_excel(input_files[["frankenberg_etal_2016_measures"]]) %>%
    dplyr::transmute(
      # Do the rename and select in one step. Drop "Thorpe ID" and
      # "Rank of size of flux for 245 sources (Frankenberg et al., 2016)."
      plume_id = `Thorpe ID`,
      latitude = Latitude,
      longitude = Longitude,
      source_type = `Source designation (Thorpe, performed after Frankenberg et al., 2016, total 273 sources with some repeat observations)`,
      emiss_kg_hr = `Estimated flux for 245 sources (Frankenberg et al., 2016). Units kg CH4/hr`
    ) %>%
    dplyr::slice(-1) %>% # omit spacer row
    dplyr::filter(
      not_na(emiss_kg_hr), # sources that had a plume, but it wasn't quantified.
      source_type %in% c("Tank", "Tanks", "Unknown", "Unknown facility",
        "Unknown infrastructure", "Well completion", "Wellpad infrastructure"
      ),
      emiss_kg_hr == 0 | emiss_kg_hr > 0.001 # one *very* small value seems like a mistake
    )
  # Source types in the data:
  # source_type                n
  # Coal mine vent shaft       1
  # Gas processing plant      10
  # Natural                    3
  # Pipeline                   3
  # Tank                      61
  # Tanks                      1
  # Unknown                   10
  # Unknown facility          17
  # Unknown infrastructure     2
  # Well completion            1
  # Wellpad infrastructure   135

  sources <- readxl::read_excel(input_files[["frankenberg_etal_2016_sources"]]) %>%
    dplyr::rename_all(make_better_names)
  multiple_divider_row <- which(sources$x == "Multiple Overpasses")
  total_rows <- nrow(sources)
  stopifnot(identical(multiple_divider_row, 179L), total_rows > 179)
  sources <- dplyr::slice(sources, -multiple_divider_row) %>%
    dplyr::select(longitude, latitude, file) %>%
    dplyr::mutate(
      flight_name =
        stringr::str_match(.data$file, "^(ang\\d{8}t\\d{6})_ch4_v1e_img$")[, 2, drop=TRUE],
      flight_date = lubridate::ymd(
        stringr::str_match(.data$file, "^ang(\\d{8})t\\d{6}_ch4_v1e_img$")[, 2, drop=TRUE]
      )
    ) %>%
    ensure_id_vars(longitude, latitude)
  stopifnot(
    !anyNA(sources$flight_name),
    nrow(dplyr::anti_join(measures, sources, by=c("longitude", "latitude"))) == 0
  )
  # For these two, longitude and latitude match exactly because they're from the
  # same underyling measurements. Some sources were detected but don't have
  # quantified measurements.
  out <- dplyr::left_join(measures, sources, by=c("longitude", "latitude")) %>%
    # Drop 4 obs that don't have associated flight paths in the online data.
    # (They're just east of the flight main zone; not sure if they were recorded the same)
    dplyr::filter(! flight_name %in% c("ang20150421t160633", "ang20150423t150648")) %>%
    dplyr::select(-file) %>%
    sf::st_as_sf(
      crs=LON_LAT_CRS,
      coords=c("longitude", "latitude"),
      agr=c(
        plume_id = "identity", source_type = "constant", emiss_kg_hr = "aggregate",
        flight_name = "constant", flight_date = "constant"
      )
    ) %>%
    sf::st_transform(OUTPUT_CRS) %>%
    # Add all-NA variables to make it easier to share some functions
    dplyr::mutate(
      emiss_se_kg_hr = NA_real_,
    )
  out
}

write_datasets <- function(data_lst, output_lst) {
  data_names <- names(data_lst)
  stopifnot(
    length(data_names) == length(data_lst),
    all(data_names %in% names(output_lst))
  )
  for (nm in data_names) {
    df <- data_lst[[nm]]
    if (inherits(df, "sf")) {
      df <- geometry_to_lonlat(df)
    }
    arrow::write_parquet(df, output_lst[[nm]])
  }
  NULL
}

standardize_columns_for_plotting <- function(df_list, censor_threshold=5) {
  # NOTE: in the output of emiss_kg_hr, zero means emissions could have been
  # detected and quantified but were not. NA means emissions could not have been
  # quantified.
  df_list$jpl_wells_all %<>%
    dplyr::transmute(
      emiss_kg_hr = dplyr::if_else(is.na(emiss_kg_hr), 0, emiss_kg_hr),
      detect_emiss = emiss_kg_hr > 0,
      src = "California",
    )
  df_list$four_corners_all_wells %<>%
    dplyr::transmute(
      emiss_kg_hr = dplyr::if_else(is.na(emiss_kg_hr), 0, emiss_kg_hr),
      detect_emiss = emiss_kg_hr > 0,
      src = "Four Corners",
    )
  df_list$lyon %<>%
    dplyr::transmute(
      emiss_kg_hr = NA_real_,
      detect_emiss = detect_emiss == 1,
      src = "Lyon et al.",
    )
  df_list$ground_studies_censored_5kgh <- df_list$ground_studies %>%
    dplyr::transmute(
      emiss_kg_hr = dplyr::if_else(.data$emiss_kg_hr > 5, .data$emiss_kg_hr, 0),
      detect_emiss = .data$emiss_kg_hr > 0,
      src = "Ground studies censored at 5 kg/hr",
    )
  df_list$ground_studies_censored_10kgh <- df_list$ground_studies %>%
    dplyr::transmute(
      emiss_kg_hr = dplyr::if_else(.data$emiss_kg_hr > 10, .data$emiss_kg_hr, 0),
      detect_emiss = .data$emiss_kg_hr > 0,
      src = "Ground studies censored at 10 kg/hr",
    )
  df_list$ground_studies %<>%
    # dplyr::filter(emiss_kg_hr > 0) %>%
    dplyr::transmute(
      emiss_kg_hr = .data$emiss_kg_hr,
      detect_emiss = .data$emiss_kg_hr > 0,
      src = "Ground studies",
    )

  df_list
}

filter_within_distance <- function(x, y, max_dist=units::as_units(1000, "m")) {
  # Keep elements of x that are within 5000 m of any element of y
  y_buffer <- sf::st_geometry(y) %>%
    sf::st_union() %>%
    sf::st_buffer(dist=max_dist)
  stopifnot(length(y_buffer) == 1)
  sf::st_intersection(x, y_buffer)
}

data_to_check_matches <- function(snakemake) {
  outfile <- snakemake@output[["data_to_check_matches"]] %||% stop("Need outfile")
  # lyon <- read_lyon_etal_2016(snakemake@input)
  well_pad_mapping <- snakemake@input[["well_pad_crosswalk"]] %>%
    arrow::read_parquet(col_select=c("entity_id", "well_pad_id")) %>%
    ensure_id_vars(entity_id)

  jpl_sites <- read_jpl_sites(snakemake@input) %>%
    dplyr::rename(source_id = source_identifier)
  wells_ca <- read_headers(2016:2017, "CA") %>%
    filter_within_distance(jpl_sites) %>%
    dplyr::inner_join(well_pad_mapping, by="entity_id")

  four_corners_df <- read_frankenberg_etal_2016(snakemake@input) %>%
    dplyr::rename(source_id = plume_id)

  wells_nm_co <- read_headers(2015, c("CO", "NM")) %>%
    filter_within_distance(four_corners_df) %>%
    dplyr::inner_join(well_pad_mapping, by="entity_id")

  df_list <- readRDS(snakemake@output[["cleaned_matched_obs"]])
  jpl_matches <- df_list$jpl_wells_all %>%
    dplyr::filter(!is.na(emiss_kg_hr)) %>%
    dplyr::select(source_identifier, entity_id) %>%
    dplyr::rename(source_id = source_identifier)
  four_corners_matches <- df_list$four_corners_all_wells %>%
    dplyr::filter(!is.na(emiss_kg_hr)) %>%
    dplyr::select(plume_id, entity_id) %>%
    dplyr::rename(source_id = plume_id)
  stopifnot(
    anyDuplicated(jpl_matches$source_id) == 0,
    anyDuplicated(jpl_matches$entity_id) == 0,
    anyDuplicated(four_corners_matches$source_id) == 0,
    anyDuplicated(four_corners_matches$entity_id) == 0
  )
  jpl_sites %<>% dplyr::left_join(jpl_matches, by="source_id")
  four_corners_df %<>% dplyr::left_join(four_corners_matches, by="source_id")
  out <- list(
    ca    = list(measures = jpl_sites,       wells = wells_ca),
    co_nm = list(measures = four_corners_df, wells = wells_nm_co)
  )
  saveRDS(out, outfile)
  invisible(out)
}

match_full_state_names <- function(df) {
  state_names <- tibble::tibble(state_abb = state.abb, state_full = state.name)
  safejoin::safe_inner_join(state_names, df,
    # 1:m merge, require a match of all rows in df
    by=c("state_abb"="state"), check="B C U L T N"
  )
}

match_commodity_prices <- function(df, price_file) {
  stopifnot("basin" %in% names(df), length(price_file) == 1)
  prices <- arrow::read_parquet(price_file, col_select=c("date", "basin", "price_real")) %>%
    dplyr::rename(gas_price_per_mcf = price_real)
  out <- dplyr::mutate(df,
      date_monthly = first_of_month(flight_date),
    ) %>%
    safejoin::safe_inner_join(
      prices,
      # Check that the join is m:1, all rows of left have a match, and there are
      # no conflicts in columns or types.
      by=c("basin"="basin", "date_monthly"="date"), check="B C V M L T"
    ) %>%
    dplyr::select(-date_monthly) %>%
    dplyr::mutate(
      gas_frac_methane = 0.95, # APPROX!
    )
  out
}

write_match_percent <- function(count_total, count_match, outfile) {
  stopifnot(
    count_total > 100, count_total < 10000, count_match > 0,
    count_match <= count_total, length(outfile) == 1
  )
  count_unmatch <- count_total - count_match
  pct_drop <- signif(100 * count_unmatch / count_total, 2)
  writeLines(paste0(pct_drop, "\\%%"), outfile)
  invisible(NULL)
}

match_jpl_california <- function(input_files) {
  # Note:
  # Some wells are flown over multiple times. Ideally, we would use that
  # information to estimate leak probabilities for these wells. However, for the
  # short run, we're going to count wells that weren't leaking on their first
  # overpass as not leaking, becuase we definitely want to avoid overstating the
  # leak rate (as using ever-leak status would do)
  #
  # The steps we follow are:
  # - Match each flight polygon with wells.
  # - For each matched well, only keep one randomly drawn flight
  # - Within groups of flight ID, match wells to detected plumes.
  # - Report counts of all plumes, matched plumes, and plumes that would have
  #   matched if we weren't dropping the later overpasses.
  plumes <- read_jpl_plumes(input_files) %>%
    dplyr::filter(!is.na(emiss_kg_hr), !is.na(emiss_se_kg_hr))
  flight_paths <- load_flight_paths(plumes)

  # Load the well covariates.
  # Note: this is could be improved a bit, because for repeat flyovers we want the
  # specific draw we have, but load_well_records returns the average for the
  # well across all visits.
  well_info <- load_well_records(flight_paths, "CA", input_files$nat_gas_prices)


  # This is like match_wells_to_flight_day, but different because we pick out
  # individual flights instead of flight days (some days have multiple flights)
  # All flights are in CA
  wells_flight_match <- read_headers(
      years=unique(lubridate::year(flight_paths$flight_date)),
      states="CA"
    ) %>%
    # Keep only wells that we kept in load_well_records (dropping gas == 0)
    dplyr::filter(.data$entity_id %in% well_info$entity_id) %>%
    # Only consider wells that are in the flight paths *and* are producing when
    # the plane flys over. This join generates duplicate rows wells because
    # the there are multiple flights over the same areas. Keep rows that had
    # any overpass while active, then drop duplicates.
    # (doing it this way because distinct doesn't work quite the same on sf objects)
    sf::st_join(
      dplyr::select(flight_paths, flight_date, flight_name),
      join=sf::st_intersects,
      left=FALSE
    ) %>%
    dplyr::filter(flight_date >= first_prod_date, flight_date <= last_prod_date) %>%
    dplyr::select(entity_id, flight_date, flight_name) %>%
    ensure_id_vars(entity_id, flight_date, flight_name)

  flight_names <- unique(wells_flight_match$flight_name)
  stopifnot(noNAs(wells_flight_match), length(flight_names) > 10)

  wells_rand_overpass <- wells_flight_match %>%
    dplyr::group_by(entity_id) %>%
    dplyr::slice_sample(n = 1) %>%
    dplyr::ungroup()

  .match_one_overpass <- function(flight_name, plumes, wells) {
    plumes %<>% dplyr::filter(.data$flight_name == !!flight_name) %>%
      dplyr::select(-flight_name)
    wells  %<>% dplyr::filter(.data$flight_name == !!flight_name)
    if (nrow(plumes) == 0 || nrow(wells) == 0) {
      return(NULL)
    }
    matched <- sf::st_join(plumes, wells,
        join=st_nn, k=1, progress=FALSE, left=FALSE
      ) %>%
      sf::st_drop_geometry() %>%
      # Now, if there are multiple plumes observed and matched to a well
      # _from the same flight_, then average them.
      dplyr::group_by(entity_id) %>%
      dplyr::summarize(
        emiss_kg_hr = mean_(emiss_kg_hr),
        emiss_se_kg_hr = mean_(emiss_se_kg_hr),
        .groups="drop"
      )
    all_wells_flown_over <- wells %>%
      sf::st_drop_geometry() %>%
      safejoin::safe_left_join(matched, by="entity_id", check="B C U V N L T")
    all_wells_flown_over
  }
  # Loop over flight_name and match wells and plumes for each flight.
  # Note: there's absolutely a better way to do this.
  all_wells_flown_over <- purrr::map_dfr(
    flight_names, .match_one_overpass,
    plumes=plumes, wells=wells_rand_overpass
  )
  stopifnot(
    anyDuplicated(all_wells_flown_over$entity_id) == 0,
    all(all_wells_flown_over$entity_id %in% well_info$entity_id)
  )
  all_wells_flown_over %<>% dplyr::inner_join(well_info, by="entity_id")
  matched_wells <- dplyr::filter(all_wells_flown_over, !is.na(emiss_kg_hr))

  observed_well_pads <- aggregate_to_well_pads(all_wells_flown_over, input_files$well_pad_crosswalk)

  could_have_been_matches <- sf::st_join(
    plumes, wells_flight_match,
    join=st_nn, k=1,
    maxdist=MAX_ACCEPTABLE_MATCH_DIST_METERS,
    left=FALSE,
    progress=FALSE
  )

  list(
    observed_well_pads = observed_well_pads,
    matched_wells = matched_wells,
    count_total = nrow(plumes),
    count_matched = sum(!is.na(all_wells_flown_over$emiss_kg_hr)),
    count_would_have_matched = nrow(could_have_been_matches)
  )
}


if (!exists("snakemake")) {
  snakemake <- SnakemakePlaceholder(
    input = list(
      well_pad_crosswalk = "data/generated/production/well_pad_crosswalk_1970-2018.parquet",
      headers = glue::glue("data/generated/production/well_headers/first_prod_year={year}/file.parquet", year=1990:2018),
      prod = glue::glue("data/generated/production/monthly_production/year={year}/file.parquet", year=1990:2018),
      alvarez_2018 = "data/studies/alvarez_etal_2018/aar7204_Database_S1.xlsx",
      omara_2018 = "data/studies/omara_etal_2018/Omara_etal_SI_tables.csv",
      duren_2019_plumes = "data/studies/duren_etal_2019/Plume_list_20191031.csv",
      duren_2019_sites = "data/studies/duren_etal_2019/41586_2019_1720_MOESM3_ESM.xlsx",
      lyon_etal_2016_locations = "data/studies/lyon_etal_2016/es6b00705_si_005.xlsx",
      lyon_etal_2016_measures = "data/studies/lyon_etal_2016/es6b00705_si_004.xlsx",
      frankenberg_etal_2016_sources = "data/studies/frankenberg_etal_2016/AVNG_sources_all2.xlsx",
      frankenberg_etal_2016_measures = "data/studies/frankenberg_etal_2016/FourCorners_AV_NG_detections_Werner.xlsx",
      zavala_araiza_2018 = "data/studies/zavala-araiza_etal_2018/elementa-6-284-s1.xlsx",
      nat_gas_prices = "data/generated/nat_gas_prices_by_basin.parquet"
    ),
    output = list(
      plot_obs_count_jpl = "graphics/observation_count_jpl_flights.pdf",
      plot_jpl_flights_qq = "graphics/jpl_flights_qq_plot.pdf",
      lyon_etal_2016 = "data/generated/methane_measures/lyon_etal_2016.parquet",
      data_to_check_matches = "data/generated/methane_measures/data_to_check_matches.rds",
      cleaned_matched_obs = "data/generated/methane_measures/matched_wells_all.rds",
      aviris_match_fraction_dropped = "output/tex_fragments/aviris_match_fraction_dropped.tex"
    ),
    threads = 4,
    resources = list(mem_mb = 7000),
    rule = ""
  )
}

main <- function(snakemake, make_extra_plots=FALSE) {
  ground_studies <- load_ground_studies(snakemake@input)
  lyon <- read_lyon_etal_2016(snakemake@input) %>% match_commodity_prices(snakemake@input$nat_gas_prices)
  four_corners_df <- read_frankenberg_etal_2016(snakemake@input)

  jpl_ca_list <- match_jpl_california(snakemake@input)
  jpl_wells_all <- jpl_ca_list$observed_well_pads
  jpl_sites_matched <- jpl_ca_list$matched_wells

  # Keep track of the total number of sites and how many were matched.
  # (Note that read_jpl_sites already drops some non-O&G sites)
  # Note: these numbers are used to determine geographic match quality, so I'm
  # including counts that would have matched but for the random-overpass filter
  # (see details in match_jpl_california)
  count_aviris_total <- jpl_ca_list$count_total
  count_aviris_match <- jpl_ca_list$count_would_have_matched


  wells_nm_co <- load_flight_paths(four_corners_df) %>%
    load_well_records(c("CO", "NM"), snakemake@input$nat_gas_prices)

  four_corners_matched <- match_with_wells(four_corners_df, wells_nm_co)
  # Keep track of the total number of sites and how many were matched.
  count_aviris_total <- count_aviris_total + nrow(four_corners_df)
  count_aviris_match <- count_aviris_match + nrow(four_corners_matched)

  four_corners_vars_to_keep <- c(setdiff(colnames(four_corners_df), "geometry"), "entity_id")
  four_corners_all_wells <- four_corners_matched %>%
    sf::st_drop_geometry() %>%
    dplyr::select(!!four_corners_vars_to_keep) %>%
    dplyr::right_join(wells_nm_co, by="entity_id") %>%
    aggregate_to_well_pads(snakemake@input[["well_pad_crosswalk"]])

  df_list <- list(
    jpl_wells_all = jpl_wells_all,
    four_corners_all_wells = four_corners_all_wells,
    lyon = lyon,
    ground_studies = ground_studies
  )
  saveRDS(df_list, snakemake@output[["cleaned_matched_obs"]])

  # NOTE: not using these plots.
  if (make_extra_plots) {
    # Make plots of JPL measures and ground studies.
    make_plots(jpl_wells_all, jpl_sites_matched, ground_studies)
    # Make plots of distributions across studies.
  }

  write_match_percent(count_aviris_total, count_aviris_match,
    snakemake@output$aviris_match_fraction_dropped
  )

  # Create output data:
  # To add to this list, make sure the list names match names in snakemake@output
  # list(lyon_etal_2016 = read_lyon_etal_2016(snakemake@input)) %>%
    # write_datasets(snakemake@output)
}

arrow::set_cpu_count(snakemake@threads)
memory_limit(snakemake@resources[["mem_mb"]])


main(snakemake)