Add Python bindings

Co-authored-by: Jon Dalrymple <[email protected]>

Author: dbaston

.gitlab-ci.yml

@@ -17,7 +17,7 @@ variables:
     - geos-config --version
     - mkdir build-coverage
     - cd build-coverage
-    - cmake -DCMAKE_BUILD_TYPE=Coverage -DBUILD_CLI=NO ..
+    - cmake -DCMAKE_BUILD_TYPE=Coverage -DBUILD_CLI=NO -DBUILD_PYTHON=NO ..
     - make -j2 catch_tests
     - valgrind --leak-check=full --error-exitcode=1 ./catch_tests
   after_script:
@@ -72,14 +72,30 @@ test:cli:
     - geos-config --version
     - mkdir build-coverage
     - cd build-coverage
-    - cmake -DCMAKE_BUILD_TYPE=Coverage -DBUILD_CLI=YES ..
+    - cmake -DCMAKE_BUILD_TYPE=Coverage -DBUILD_CLI=YES -DBUILD_PYTHON=NO ..
     - make -j2 exactextract_bin subdivide
     - pytest ../test
   after_script:
     - cd build-coverage
     - lcov --capture --directory CMakeFiles --output-file coverage.info
     - bash <(curl -s https://codecov.io/bash)
 
+test:python:
+  stage: test
+  image: isciences/exactextract-test-env:geos312
+  script:
+    - apt-get install -y python3-dev pybind11-dev
+    - geos-config --version
+    - mkdir build-coverage
+    - cd build-coverage
+    - cmake -DCMAKE_BUILD_TYPE=Coverage -DBUILD_CLI=NO -DBUILD_PYTHON=YES ..
+    - cmake --build . -j2
+    - ctest -R pybindings --output-on-failure
+  after_script:
+    - cd build-coverage
+    - lcov --capture --directory CMakeFiles --output-file coverage.info
+    - bash <(curl -s https://codecov.io/bash)
+
 build:
   stage: build
   script:

CMakeLists.txt

@@ -9,6 +9,8 @@ set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
+
+include(CTest)
 include(GNUInstallDirs)
 
 set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake)
@@ -19,6 +21,7 @@ find_package(GEOS REQUIRED)
 #Configure some options the various components this module can build
 option(BUILD_CLI "Build the exactextract cli binary" ON) #requires gdal, cli11
 option(BUILD_TEST "Build the exactextract tests" ON) #requires catch
+option(BUILD_PYTHON "Build the exactextract Python bindings" ON) # requires pybind11
 option(BUILD_DOC "Build documentation" ON) #requires doxygen
 
 if(BUILD_CLI)
@@ -246,6 +249,7 @@ set(PROJECT_SOURCES
 
 add_library(${LIB_NAME} ${PROJECT_SOURCES})
 
+
 # Check matrix bounds for debug builds
 set_target_properties(${LIB_NAME}
         PROPERTIES COMPILE_DEFINITIONS $<$<CONFIG:Debug>:MATRIX_CHECK_BOUNDS>)
@@ -277,6 +281,10 @@ target_link_libraries(
 
 set_target_properties(${LIB_NAME} PROPERTIES OUTPUT_NAME ${LIB_NAME})
 
+if (BUILD_PYTHON)
+    add_subdirectory(python)
+endif() # BUILD_PYTHON
+
 if(BUILD_DOC)
     # Doxygen configuration from https://vicrucann.github.io/tutorials/quick-cmake-doxygen/
     # check if Doxygen is installed

python/CMakeLists.txt

@@ -0,0 +1,38 @@
+# Build the Python bindings, which require pybind11
+find_package(pybind11 REQUIRED)
+if (pybind11_FOUND)
+    # TODO: which pybind11 version
+endif() #pybind11_FOUND
+if (NOT pybind11_FOUND)
+    message(FATAL_ERROR
+    "pybind11 was not found. It is still possible to build and test libexactextract, but the "
+    "exactextract Python bindings cannot be built or installed.")
+endif() #NOT pybind11_FOUND
+
+
+set(PYBIND_SOURCES
+src/pybindings/bindings.cpp
+src/pybindings/feature_bindings.cpp
+src/pybindings/feature_bindings.h
+src/pybindings/feature_source_bindings.cpp
+src/pybindings/feature_source_bindings.h
+src/pybindings/operation_bindings.cpp
+src/pybindings/operation_bindings.h
+src/pybindings/processor_bindings.cpp
+src/pybindings/processor_bindings.h
+src/pybindings/raster_source_bindings.cpp
+src/pybindings/raster_source_bindings.h
+src/pybindings/writer_bindings.cpp
+src/pybindings/writer_bindings.h)
+
+pybind11_add_module(_exactextract MODULE ${PYBIND_SOURCES})
+target_include_directories(_exactextract PRIVATE
+    ${CMAKE_SOURCE_DIR}/src
+    ${GEOS_INCLUDE_DIR})
+    set_property(TARGET ${LIB_NAME} PROPERTY POSITION_INDEPENDENT_CODE ON)
+    target_link_libraries(_exactextract PRIVATE ${LIB_NAME} ${GEOS_LIBRARY})
+
+add_test(NAME "pybindings" 
+         COMMAND ${CMAKE_COMMAND} -E env
+         PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}:${CMAKE_CURRENT_LIST_DIR}/src
+         python3 -m pytest ${CMAKE_CURRENT_LIST_DIR}/tests)

python/src/exactextract/__init__.py

@@ -0,0 +1,10 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+""" Python bindings for exactextract """
+
+from .feature import Feature, GDALFeature, JSONFeature
+from .feature_source import FeatureSource, GDALFeatureSource, JSONFeatureSource
+from .operation import Operation
+from .processor import FeatureSequentialProcessor, RasterSequentialProcessor
+from .raster_source import RasterSource, GDALRasterSource, NumPyRasterSource
+from .writer import Writer, JSONWriter, GDALWriter

python/src/exactextract/feature.py

@@ -0,0 +1,53 @@
+from _exactextract import Feature
+
+
+class GDALFeature(Feature):
+    def __init__(self, f):
+        Feature.__init__(self)
+        self.feature = f
+
+    def set(self, name, value):
+        self.feature.SetField(name, value)
+
+    def get(self, name):
+        return self.feature.GetField(name)
+
+    def geometry(self):
+        return bytes(self.feature.GetGeometryRef().ExportToWkb())
+
+    def fields(self):
+        defn = self.feature.GetDefnRef()
+        return [defn.GetFieldDefn(i).GetName() for i in range(defn.GetFieldCount())]
+
+
+class JSONFeature(Feature):
+    def __init__(self, f=None):
+        Feature.__init__(self)
+        self.feature = f or {}
+
+    def set(self, name, value):
+        if name == "id":
+            self.feature["id"] = value
+        else:
+            if "properties" not in self.feature:
+                self.feature["properties"] = {}
+            self.feature["properties"][name] = value
+
+    def get(self, name):
+        if name == "id":
+            return self.feature["id"]
+        else:
+            return self.feature["properties"][name]
+
+    def geometry(self):
+        import json
+
+        return json.dumps(self.feature["geometry"])
+
+    def fields(self):
+        fields = []
+        if "id" in self.feature:
+            fields.append("id")
+        if "properties" in self.feature:
+            fields += self.feature["properties"].keys()
+        return fields

python/src/exactextract/feature_source.py

@@ -0,0 +1,26 @@
+from _exactextract import FeatureSource
+
+from .feature import GDALFeature, JSONFeature
+
+
+class GDALFeatureSource(FeatureSource):
+    def __init__(self, src, id_field):
+        super().__init__(id_field)
+        self.src = src
+
+    def __iter__(self):
+        for f in self.src:
+            yield GDALFeature(f)
+
+
+class JSONFeatureSource(FeatureSource):
+    def __init__(self, src, id_field):
+        super().__init__(id_field)
+        if type(src) is dict:
+            self.src = [src]
+        else:
+            self.src = src
+
+    def __iter__(self):
+        for f in self.src:
+            yield JSONFeature(f)

python/src/exactextract/operation.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+from __future__ import annotations
+
+from typing import Optional
+
+from _exactextract import Operation as _Operation
+
+from .raster_source import RasterSource
+
+
+class Operation(_Operation):
+    """Binding class around exactextract Operation"""
+
+    def __init__(
+        self,
+        stat_name: str,
+        field_name: str,
+        raster: RasterSource,
+        weights: Optional[RasterSource] = None,
+    ):
+        """
+        Create Operation object from stat name, field name, raster, and weighted raster
+
+        Args:
+            stat_name (str): Name of the stat. Refer to docs for options.
+            field_name (str): Field name to use. Output of operation will have title \'{field_name}_{stat_name}\'
+            raster (RasterSource): Raster to compute over.
+            weights (Optional[RasterSource], optional): Weight raster to use. Defaults to None.
+        """
+        super().__init__(stat_name, field_name, raster, weights)

python/src/exactextract/processor.py

@@ -0,0 +1,47 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from typing import List, Union
+
+from _exactextract import (
+    FeatureSequentialProcessor as _FeatureSequentialProcessor,
+    RasterSequentialProcessor as _RasterSequentialProcessor,
+)
+
+from .feature_source import FeatureSource
+from .operation import Operation
+from .writer import Writer
+
+
+class FeatureSequentialProcessor(_FeatureSequentialProcessor):
+    """Binding class around exactextract FeatureSequentialProcessor"""
+
+    def __init__(self, ds: FeatureSource, writer: Writer, op_list: List[Operation]):
+        """
+        Create FeatureSequentialProcessor object
+
+        Args:
+            ds (FeatureSource): Dataset to use
+            writer (Writer): Writer to use
+            op_list (List[Operation]): List of operations
+        """
+        super().__init__(ds, writer)
+        for op in op_list:
+            self.add_operation(op)
+
+
+class RasterSequentialProcessor(_RasterSequentialProcessor):
+    """Binding class around exactextract RasterSequentialProcessor"""
+
+    def __init__(self, ds: FeatureSource, writer: Writer, op_list: List[Operation]):
+        """
+        Create RasterSequentialProcessor object
+
+        Args:
+            ds (FeatureSource): Dataset to use
+            writer (Writer): Writer to use
+            op_list (List[Operation]): List of operations
+        """
+        super().__init__(ds, writer)
+        for op in op_list:
+            self.add_operation(op)

python/src/exactextract/raster_source.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from __future__ import annotations
+
+import pathlib
+from typing import Optional, Tuple, Union
+from osgeo import gdal
+
+from _exactextract import RasterSource
+
+
+class GDALRasterSource(RasterSource):
+    def __init__(self, ds, band_idx: int = 1):
+        super().__init__()
+        self.ds = ds
+
+        # Sanity check inputs
+        if band_idx is not None and band_idx <= 0:
+            raise ValueError("Raster band index starts from 1!")
+
+        self.band = self.ds.GetRasterBand(band_idx)
+
+    def res(self):
+        gt = self.ds.GetGeoTransform()
+        return gt[1], abs(gt[5])
+
+    def extent(self):
+        gt = self.ds.GetGeoTransform()
+
+        dx, dy = self.res()
+
+        left = gt[0]
+        right = left + dx * self.ds.RasterXSize
+        top = gt[3]
+        bottom = gt[3] - dy * self.ds.RasterYSize
+
+        return (left, bottom, right, top)
+
+    def read_window(self, x0, y0, nx, ny):
+        return self.band.ReadAsArray(xoff=x0, yoff=y0, win_xsize=nx, win_ysize=ny)
+
+
+class NumPyRasterSource(RasterSource):
+    def __init__(self, mat, xmin, ymin, xmax, ymax):
+        super().__init__()
+        self.mat = mat
+        self.ext = (xmin, ymin, xmax, ymax)
+
+    def res(self):
+        ny, nx = self.mat.shape
+        dy = (self.ext[3] - self.ext[1]) / ny
+        dx = (self.ext[2] - self.ext[0]) / nx
+
+        return (dx, dy)
+
+    def extent(self):
+        return self.ext
+
+    def read_window(self, x0, y0, nx, ny):
+        return self.mat[y0 : y0 + ny, x0 : x0 + ny]