Add Safeguard Mechanism facilities to oil & gas sector

notebooks/oil_gas_locations.py

@@ -0,0 +1,78 @@
+import datetime
+import geopandas as gpd
+import os
+import pandas as pd
+import pathlib
+from openmethane_prior.lib import PriorConfig
+from openmethane_prior.lib.data_manager.manager import DataManager
+from openmethane_prior.sectors.oil_gas.data.nopta_titles import nopta_titles_data_source
+from openmethane_prior.sectors.oil_gas.data.nopta_wells import nopta_wells_data_source
+from openmethane_prior.sectors.oil_gas.data.nsw_drillholes import nsw_drillholes_data_source
+from openmethane_prior.sectors.oil_gas.data.nsw_titles import nsw_titles_data_source
+from openmethane_prior.sectors.oil_gas.data.qld_boreholes import qld_boreholes_data_source
+from openmethane_prior.sectors.oil_gas.data.qld_leases import qld_leases_data_source
+from openmethane_prior.sectors.oil_gas.data.wa_titles import wa_titles_data_source
+from openmethane_prior.sectors.oil_gas.data.wa_wells import wa_wells_data_source
+from openmethane_prior.sectors.oil_gas.emission_sources.nsw_sources import nsw_emission_sources
+from openmethane_prior.sectors.oil_gas.emission_sources.offshore_sources import offshore_emission_sources
+from openmethane_prior.sectors.oil_gas.emission_sources.qld_sources import qld_emission_sources
+from openmethane_prior.sectors.oil_gas.emission_sources.wa_sources import wa_emission_sources
+os.environ["START_DATE"] = "2023-07-01"
+os.environ["END_DATE"] = "2023-07-01"
+prior_config = PriorConfig.from_env()
+data_manager = DataManager(
+    data_path=pathlib.Path("../data/inputs"),
+    prior_config=prior_config,
+)
+
+nsw_drillholes_da = data_manager.get_asset(nsw_drillholes_data_source)
+nsw_titles_da = data_manager.get_asset(nsw_titles_data_source)
+nsw_df = nsw_emission_sources(
+    start_date=datetime.date(2021, 1, 1),
+    end_date=datetime.date(2023, 1, 2),
+    nsw_drillholes_da=nsw_drillholes_da,
+    nsw_titles_da=nsw_titles_da,
+)
+qld_boreholes_da = data_manager.get_asset(qld_boreholes_data_source)
+qld_leases_da = data_manager.get_asset(qld_leases_data_source)
+qld_df = qld_emission_sources(
+    start_date=datetime.date(2021, 1, 1),
+    end_date=datetime.date(2023, 1, 2),
+    qld_boreholes_da=qld_boreholes_da,
+    qld_leases_da=qld_leases_da,
+)
+wa_wells_da = data_manager.get_asset(wa_wells_data_source)
+wa_titles_da = data_manager.get_asset(wa_titles_data_source)
+wa_df = wa_emission_sources(
+    start_date=datetime.date(2021, 1, 1),
+    end_date=datetime.date(2023, 1, 2),
+    wa_wells_da=wa_wells_da,
+    wa_titles_da=wa_titles_da,
+)
+offshore_wells_da = data_manager.get_asset(nopta_wells_data_source)
+offshore_titles_da = data_manager.get_asset(nopta_titles_data_source)
+offshore_df = offshore_emission_sources(
+    start_date=datetime.date(2021, 1, 1),
+    end_date=datetime.date(2023, 1, 2),
+    offshore_wells_da=offshore_wells_da,
+    offshore_titles_da=offshore_titles_da,
+)
+all_df = pd.concat([
+    nsw_df,
+    qld_df,
+    wa_df,
+    offshore_df,
+])
+
+au_shp = gpd.GeoDataFrame.from_file("~/Downloads/STE_2021_AUST_SHP_GDA2020/STE_2021_AUST_GDA2020.shx")
+au_map = au_shp.plot(color='white', edgecolor='#CCCCCC', figsize=(32,24))
+au_map.set_xlim(110, 130)
+au_map.set_ylim(-35, -10)
+nsw_titles_da.data.plot(ax=au_map, alpha=0.5)
+qld_leases_da.data.plot(ax=au_map, alpha=0.5)
+# wa_wells_da.data.plot(ax=au_map, markersize=3, column="class", legend=True)
+wa_titles_da.data.plot(ax=au_map, alpha=0.5)
+all_df.plot(ax=au_map, markersize=3, column="data_source", legend=True)
+au_map
+
+

pyproject.toml

@@ -18,14 +18,16 @@ dependencies = [
     "rioxarray>=0.15.5,<0.16",
     "pyproj>=3.6.1,<4",
     "pandas>=2.2.2,<3",
-    "geopandas>=0.14.4,<0.15",
+    "geopandas>=1.1.2",
     "python-dotenv>=1.0.1,<2",
     "colorama>=0.4.6,<0.5",
     "cdsapi>=0.7.3,<0.8",
     "shapely>=2.0.4,<3",
     "environs>=11.0.0,<12",
     "prettyprinter>=0.18.0,<0.19",
     "pytest-mock>=3.15.1,<4",
+    "owslib>=0.35.0",
+    "bmi-arcgis-restapi>=2.4.15",
 ]
 
 [dependency-groups]

src/openmethane_prior/data_sources/safeguard/data.py

@@ -58,6 +58,7 @@
     "safeguard_facility_name", # exact facility_name from Safeguard Mechanism
     "data_source_name", # DataSource name with facility details
     "data_source_id", # identifier in reference DataSource
+    "notes", # notes about how this location is related to the facility
 ]
 
 
@@ -75,6 +76,7 @@ def parse_location_csv(data_source: ConfiguredDataSource):
         filepath_or_buffer=data_source.asset_path,
         header=0,
         names=safeguard_locations_csv_columns,
+        usecols=["safeguard_facility_name", "data_source_name", "data_source_id"],
     )
 
     # filter out any rows with incomplete information
@@ -102,7 +104,7 @@ def parse_safeguard_csv(data_source: ConfiguredDataSource):
 
 safeguard_locations_data_source = DataSource(
     name="safeguard-locations",
-    url="https://openmethane.s3.amazonaws.com/prior/inputs/facility-locations-v1.csv",
+    url="https://openmethane.s3.amazonaws.com/prior/inputs/facility-locations-v1.1.csv",
     parse=parse_location_csv,
 )
 

src/openmethane_prior/lib/__init__.py

@@ -23,13 +23,15 @@
     DataManager,
     DataSource,
 )
+from .grid.grid import Grid
+from .grid.geometry import polygon_cell_intersection
 from .grid.regrid import regrid_data
 from .outputs import add_sector, convert_to_timescale
 from .create_prior import create_prior
 from .raster import remap_raster
 from .sector.config import PriorSectorConfig
 from .sector.sector import PriorSector
-from .units import kg_to_period_cell_flux
+from .units import kg_to_period_cell_flux, days_in_period
 from .utils import (
     area_of_rectangle_m2,
     datetime64_to_datetime,

src/openmethane_prior/lib/data_manager/parsers.py

@@ -15,10 +15,22 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 #
+import geopandas as gpd
 import pandas as pd
 
 from openmethane_prior.lib.data_manager.source import ConfiguredDataSource
 
 
 def parse_csv(data_source: ConfiguredDataSource) -> pd.DataFrame:
+    """Read and parse a ConfiguredDataSource CSV asset as a pandas DataFrame."""
     return pd.read_csv(data_source.asset_path)
+
+
+def parse_geo(data_source: ConfiguredDataSource):
+    """Read and parse a file containing a collection of geometry vector data
+    into a geopandas GeoDataFrame. Asset file type can be anything supported by
+    pyogrio, which includes GeoJSON, GeoPackage, Shapefiles, etc."""
+    geo_df = gpd.read_file(data_source.asset_path)
+    if geo_df.crs is not None and geo_df.crs != "EPSG:4326":
+        geo_df = geo_df.to_crs(epsg=4326)
+    return geo_df

src/openmethane_prior/lib/grid/geometry.py

@@ -0,0 +1,62 @@
+#
+# Copyright 2025 The Superpower Institute Ltd.
+#
+# This file is part of Open Methane.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+import numpy as np
+from shapely import Polygon, MultiPolygon
+import typing
+
+from .grid import Grid
+
+# a cell intersection is a tuple containing a coordinate tuple, and a fraction
+# of intersection, representing how much of a shape falls within each cell.
+CellIntersection = typing.TypeVar("CellIntersection", bound=tuple[tuple[float, float], float])
+
+def polygon_cell_intersection(
+    geometry: Polygon | MultiPolygon,
+    grid: Grid,
+) -> list[CellIntersection]:
+    """Determine which grid cells a polygon intersects with, returning a list
+    of grid cell indexes and their ratio of intersection with the polygon.
+
+    Polygon shapes should have their coordinates specified in lat/lon in
+    the same CRS used by the Grid."""
+    # find the smallest area of the grid to examine based on the bounding box
+    bounds_llc_lon, bounds_llc_lat, bounds_urc_lon, bounds_urc_lat = geometry.bounds
+    bbox_ix, bbox_iy, _ = grid.lonlat_to_cell_index(
+        np.array([bounds_llc_lon, bounds_urc_lon, bounds_urc_lon, bounds_llc_lon]),
+        np.array([bounds_llc_lat, bounds_llc_lat, bounds_urc_lat, bounds_urc_lat]),
+    )
+    min_ix, max_ix = np.min(bbox_ix), np.max(bbox_ix)
+    min_iy, max_iy = np.min(bbox_iy), np.max(bbox_iy)
+
+    cell_bounds_lon, cell_bounds_lat = grid.cell_bounds_lonlat()
+
+    intersections: list[CellIntersection] = []
+    for iy in range(min_iy, max_iy + 1):
+        for ix in range(min_ix, max_ix + 1):
+            # if part of the shape falls outside the grid, skip those areas
+            if not grid.valid_cell_coords(ix, iy):
+                continue
+
+            # find the intersection between the geometry and the grid cell
+            cell_polygon = Polygon(zip(cell_bounds_lon[iy, ix], cell_bounds_lat[iy, ix]))
+            intersection = geometry.intersection(cell_polygon)
+
+            if intersection is not None and intersection.area > 0:
+                intersections.append(((ix, iy), intersection.area / geometry.area))
+
+    return intersections

src/openmethane_prior/lib/utils.py

@@ -29,6 +29,8 @@
 import typing
 
 import numpy as np
+import pandas as pd
+import geopandas as gpd
 from numpy.typing import ArrayLike
 from urllib.parse import urlparse
 import xarray as xr
@@ -196,3 +198,23 @@ def is_url(maybe_url: str) -> bool:
     """
     parsed = urlparse(maybe_url)
     return parsed.scheme != "" and parsed.netloc != ""
+
+
+DF = typing.TypeVar("DF", bound=pd.DataFrame | gpd.GeoDataFrame)
+def rows_in_period(
+    df: DF,
+    start_date: datetime.date,
+    end_date: datetime.date,
+    start_field: str = "start_date",
+    end_field: str = "end_date",
+) -> DF:
+    """Return the rows of the DataFrame df where the period between df[start_field]
+    and df[end_field] overlaps with the period between start_date and end_date."""
+    midnight = datetime.time(0, 0)
+    period_start = datetime.datetime.combine(date=start_date, time=midnight)
+    period_end = datetime.datetime.combine(date=end_date + datetime.timedelta(days=1), time=midnight)
+
+    return df[
+        (df[start_field] <= np.datetime64(period_end))
+        & ((np.isnat(df[end_field])) | (df[end_field] >= np.datetime64(period_start)))
+    ]

src/openmethane_prior/sectors/oil_gas/README.md

@@ -0,0 +1,94 @@
+
+# Oil and Gas Sector
+
+The oil and gas sector is the most complicated sector implementation in the
+Australia prior. There's no pre-existing high quality dataset that can be
+used to both locate and scale emission sources from the sector, and no clear
+spatial proxy.
+
+Oil and gas extraction occurs in "fields" of wells or bores, somtimes located
+offshore. Extracted resources are transported through pipelines to various
+pieces of infrastructure to be compressed or combined, and eventually to
+refinery facilities. Methane can be lost to the atmosphere for various reasons
+at any of these locations, resulting in physically dispersed infrastructure
+which is not well documented by the private companies which operate it.
+
+## Methodology
+
+Emissions are likely to occur in three possible parts of the oil and gas
+infrastructure:
+- wells or boreholes where oil or gas are extracted
+- processing facilities
+- pipelines connecting wells and other infrastructure
+
+### Wells and boreholes
+
+#### Locations
+
+Oil and gas extraction are covered by petroleum mining laws in Australia, and
+most Australian states provide public datasets of petroleum mining titles
+(aka leases/licenses/tenements). Some states also provide public datasets
+of drill/bore/well locations.
+
+Drillhole datasets typically include enough detail to determine which bores
+or wells are used for petroleum production, giving us precise locations for
+emission sources.
+
+#### Activity period
+
+Drillhole datasets typically include the approximate date a bore was drilled
+(sometimes just the year), but don't include the date a bore was capped or
+depleted. Although some bore datasets do include the bore status, ie CAPPED,
+this is the status **at the time the dataset was last updated**, and not
+necessarily the status in the period of interest for the prior.
+
+We use a very rough method to determine whether a single bore may have been
+producing emissions during a target period (`start_date` to `end_date`).
+- if the bore drill date is later than `end_date`, no emissions are possible
+- correlate the bore coordinates with all petroleum production titles it
+  is contained by
+- if the range from `start_date` to `end_date` overlaps with the range from
+  title grant date to title expiry date, emissions are possible
+
+This method has obvious problems:
+
+1. emissions are unlikely to start immediately after drilling, however without
+   a more accurate start date, this is at least a lower bound
+2. most bores/wells will be depleted or capped well before the petroleum title
+   expires, however without an accurate capping date, this is at least an upper
+   bound.
+3. there is potentially a long period (years) where an entire field / title
+   has stopped producing, but the title is still active.
+
+On a coarse grid such as our 10x10km grid, 1. and 2. probably aren't serious,
+because emissions will be dispersed amongst a number of points within each grid
+cell.
+
+The biggest issue is 3., as it will place emissions in some fields, possibly
+for years after production has ceased. Assuming that capped wells don't
+continue to emit methane, this will misallocate potentially significant
+emissions. Unfortunately we currently don't have a better solution.
+
+### Processing facilities
+
+Not yet implemented.
+
+### Pipelines
+
+Public datasets do exist which record the shape and status of oil and gas
+pipelines in Australia, however there are several issues using these shape
+files to estimate emissions.
+
+1. Lack of temporal data
+
+Pipeline datasets list the current status of pipelines, but don't necessarily
+include the dates pipelines actually began operating.
+
+2. Lack of leakage estimates
+
+We're not currently aware of any research or data on how much methane leaks
+from pipelines, or where leakage occurs. This makes it hard to attribute a
+specific volume of emission to specific locations along any given pipeline.
+
+For now, pipelines are not included in our bottom up estimate. If better data
+or research becomes available in the future, this would be a welcome addition.

src/openmethane_prior/sectors/oil_gas/data/esri_types.py

@@ -0,0 +1,39 @@
+#
+# Copyright 2026 The Superpower Institute Ltd.
+#
+# This file is part of OpenMethane.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+import datetime
+import numpy as np
+
+
+def map_esri_date_to_date(esri_date_milliseconds) -> datetime.date | None:
+    """Convert an esriFieldTypeDate value to UTC datetime.date."""
+    if esri_date_milliseconds is None \
+            or np.isnan(esri_date_milliseconds) \
+            or esri_date_milliseconds <= 0:
+        return None
+
+    date_value = datetime.datetime.fromtimestamp(
+        esri_date_milliseconds / 1000,
+        tz=datetime.timezone.utc,
+    )
+    return date_value.date()
+
+
+def map_esri_date_to_str(esri_date_milliseconds) -> str | None:
+    """Convert an esriFieldTypeDate value to RFC3339 date string."""
+    date = map_esri_date_to_date(esri_date_milliseconds)
+    return date.isoformat() if date is not None else None