Unhandled type of field in `exactextract`

[JakubCha]

In the case of vector input with some field types set to 0 length as Decimal (double) field in the screenshot below exactextract library raises RuntimeError.

The field unknown_field_length is the one raising error. Integer fields with 0 lengths are not raising errors.
The same error is raised if String field has 0 length.

Full exception traceback:

2025-03-17T20:02:25     WARNING    An error has occurred while executing Python code:
             
             
             RuntimeError: Unhandled type for field unknown_length_field in PyFeatureBase::field_type
             
             
             Traceback (most recent call last):
              File "C:\Users//AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\zonal_exact\task_classes.py", line 76, in run
              result_stats = exact_extract(
              File "E:\OSGeo4W\apps\Python39\lib\site-packages\exactextract\exact_extract.py", line 476, in exact_extract
              processor.process()
             RuntimeError: Unhandled type for field unknown_length_field in PyFeatureBase::field_type
             
             
             
             
             Python version: 3.9.18 (heads/master:5eba59e, Feb 1 2024, 20:02:10) [MSC v.1929 64 bit (AMD64)]
             
             
             QGIS version: 3.36.0-Maidenhead Maidenhead, 09951dc0

My first intuition is that the error comes from the fact that QGIS use Qt field types (https://doc.qt.io/qt-6/qttypes.html) and at some point in exactextract it's incorrectly converted to the proper field type (0 field length).

[JakubCha]

It might be the error that some people got before but I wasn't able to find it. #31

[JakubCha]

The proper solution is to find why this error is raised with such field types in the first place but it'd probably require debugging down to the C++ level of exactextract: https://github.com/isciences/exactextract/blob/a5575814952f3be9b8f708199821c1502a19f48b/python/src/pybindings/feature_bindings.cpp#L128

Another possible solution (or rather a workaround) is to extract only ID and geometry columns from the input, pass it to exactextract and merge other columns back at the QGIS level.

[JakubCha]

exactextract seems to be failing only with QgsVectorLayer created with ogr (OGR Data Provider), QgsVectorLayer created with memory data provider with the same fields is working fine.

Minimum example:

from qgis.core import (
    QgsCoordinateReferenceSystem,
    QgsVectorLayer,
    QgsRasterLayer,
    QgsGeometry,
    QgsFeature,
    QgsField,
    QgsFields,
    QgsProject,
)
from PyQt5.QtCore import QVariant

from zonal_exact.task_classes import CalculateStatsTask

vector_layer = QgsVectorLayer("Polygon?crs=epsg:3035", "polygon_layer", "ogr")
provider = vector_layer.dataProvider()

# Define the fields
fields = QgsFields()
fields.append(QgsField("id", QVariant.Int))
# add decimal field with undefined length
fields.append(QgsField("decimal", QVariant.Double))

# Add the fields to the layer
provider.addAttributes(fields)
vector_layer.updateFields()

# Create 5 polygons in the left side of the raster
for i in range(0, 6):
    feature = QgsFeature()
    feature.setGeometry(
        QgsGeometry.fromWkt(f"POLYGON((0 {i-1}, 1 {i-1}, 1 {i}, 0 {i}, 0 {i-1}))")
    )
    feature.setAttributes([i, i+0.1])
    provider.addFeature(feature)
# Create 5 polygons in the right side of the raster
for i in range(5, 0, -1):
    feature = QgsFeature()
    feature.setGeometry(
        QgsGeometry.fromWkt(f"POLYGON((5 {i}, 4 {i}, 4 {i+1}, 5 {i+1}, 5 {i}))")
    )
    feature.setAttributes([i + 10, i+10+0.1])
    provider.addFeature(feature)
# Create huge polygon spanning over most part of the raster
feature = QgsFeature()
feature.setGeometry(QgsGeometry.fromWkt("POLYGON((1 -1, 5 -1, 1 5, 1 -1))"))
feature.setAttributes([20, 20+0.1])
provider.addFeature(feature)

# Create a raster file
raster_file = str(tmp_path / "pytest_raster.tif")
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(raster_file, 6, 6, 2, gdal.GDT_Float32)
dataset.SetGeoTransform(
    [0, 1, 0, 5, 0, -1]
)  # Define the affine transform (how raster is scaled, rotated, skewed, and/or translated)

nodata_value = -9999.0
for i in range(1, 3):
    band = dataset.GetRasterBand(i)
    band.SetNoDataValue(nodata_value)
    data_array = np.array([[(j + i) + (k + i) for j in range(6)] for k in range(6)])
    data_array[5, :] = nodata_value  # Set last row to NoData value
    data_array[:, 5] = nodata_value  # Set last column to NoData value
    band.WriteArray(data_array)

dataset.FlushCache()

# Create a raster layer from the raster file
raster_layer = QgsRasterLayer(raster_file, "raster_layer")
raster_layer.setCrs(QgsCoordinateReferenceSystem("EPSG:3035"))

raster_layer_path = raster_layer.dataProvider().dataSourceUri()

# Create a console writer
stats_to_calculate = ["mean", "max", "min"]
# Initialize the task and run it
task = CalculateStatsTask(
        "Test Task",
        QgsTask.CanCancel,
        [],
        vector_layer,
        [raster_layer_path],
        None,  # No weights
        stats_to_calculate,
        ["fid", "decimal"],
        geospatial_output=False,
        strategy="feature-sequential",
    )
task.run()

Full error output:

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
zonal_exact/task_classes.py:90: in run
    result_stats = exact_extract(
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

rast = [<exactextract.raster.GDALRasterSource object at 0x7f0ee6841670>, <exactextract.raster.GDALRasterSource object at 0x7f0ee6841790>]
vec = <exactextract.feature.QGISFeatureSource object at 0x7f0ee68417f0>
ops = [<exactextract._exactextract.Operation object at 0x7f0ee66c87f0>, <exactextract._exactextract.Operation object at 0x7f...act._exactextract.Operation object at 0x7f0ee66c86f0>, <exactextract._exactextract.Operation object at 0x7f0ee66c8730>]

    def exact_extract(
        rast,
        vec,
        ops,
        *,
        weights=None,
        include_cols=None,
        include_geom=False,
        strategy: str = "feature-sequential",
        max_cells_in_memory: int = 30000000,
        grid_compat_tol: float = 1e-3,
        output: str = "geojson",
        output_options: Optional[Mapping] = None,
        progress=False,
    ):
        """Calculate zonal statistics.
    
        Args:
           rast: A :py:class:`RasterSource` or filename that can be opened
                 by GDAL/rasterio/xarray.
           vec: A :py:class:`FeatureSource` or filename that can be opened
                 by GDAL/fiona
           ops: A list of :py:class:`Operation` objects, or strings that
                can be used to construct them (e.g., ``"mean"``, ``"quantile(q=0.33)"``).
                Check out :doc:`Available operations </operations>` for more information.
           weights: An optional :py:class:`RasterSource` or filename for
                weights to be used in weighted operations.
           include_cols: An optional list of columns from the
                input features to be included into the output.
           include_geom: Flag indicating whether the geometry
                should be copied from the input features into
                the output.
           strategy: Specifies the strategy to use when processing features.
              Detailed performance notes are available in
              :ref:`Performance - Processing strategies <performance-processing-strategies>`.
              Must be set to one of:
    
                     - ``"feature-sequential"`` (the default):
                       iterate over the features in ``vec``,
                       read the corresponding pixels from ``rast``/``weights``,
                       and compute the summary operations. This offers predictable
                       memory consumption but may be inefficient if the order of
                       features in ``vec`` causes the same, relatively large raster
                       blocks to be read and decompressed many times.
    
                     - ``"raster-sequential"``:
                       iterate over chunks of pixels in ``rast``, identify the intersecting
                       features from ``vec``, and compute statistics. This performs better
                       than ``strategy="feature-sequential"`` in some cases, but comes at
                       a cost of higher memory usage.
           max_cells_in_memory: Indicates the maximum number of raster cells that should be
                                loaded into memory at a given time.
           grid_compat_tol: require value and weight grids to align within ``grid_compat_tol`` times the smaller of the two grid resolutions
           output: An :py:class:`OutputWriter` or one of the following strings:
    
                     - "geojson" (the default): return a list of GeoJSON-like features
                     - "pandas": return a ``pandas.DataFrame`` or ``geopandas.GeoDataFrame``,
                                 depending on the value of ``include_geom``
                     - "gdal": write results to disk using GDAL/OGR as they are written. This
                               option (with ``strategy="feature-sequential"``) avoids the need
                               to maintain results for all features in memory at a single time,
                               which may be significant for operations with large result sizes
                               such as ``cell_id``, ``values``, etc.
           output_options: an optional dictionary of options passed to the :py:class:`writer.JSONWriter`, :py:class:`writer.PandasWriter`, or :py:class:`writer.GDALWriter`.
           progress: if `True`, a progress bar will be displayed. Alternatively, a
                     function may be provided that will be called with the completion fraction
                     and a status message.
        """
        rast = prep_raster(rast)
        weights = prep_raster(weights, name_root="weight")
        vec = prep_vec(vec)
    
        if output_options is None:
            output_options = {}
    
        ops = prep_ops(
            ops, rast, weights, add_unique=output_options.get("frac_as_map", False)
        )
    
        warn_on_crs_mismatch(vec, ops)
    
        if "frac_as_map" in output_options:
            output_options = copy.copy(output_options)
            del output_options["frac_as_map"]
            output_options["map_fields"] = {
                "frac": ("unique", "frac"),
                "weighted_frac": ("unique", "weighted_frac"),
            }
    
        if type(include_cols) is str:
            include_cols = [include_cols]
    
        Processor = prep_processor(strategy)
    
        writer = prep_writer(output, vec.srs_wkt(), output_options)
    
        processor = Processor(vec, writer, ops, include_cols)
        if include_geom:
            processor.add_geom()
        processor.set_max_cells_in_memory(max_cells_in_memory)
        processor.set_grid_compat_tol(grid_compat_tol)
    
        if progress:
            processor.show_progress(True)
    
            if callable(progress):
                processor.set_progress_fn(progress)
            elif progress is True:
                try:
                    import tqdm
    
                    bar = tqdm.tqdm(total=100)
    
                    def status(frac, message):
                        pct = frac * 100
                        bar.update(pct - bar.n)
                        bar.set_description(message)
                        if pct == 100:
                            bar.close()
    
                except ImportError:
    
                    def status(frac, message):
                        print(f"[{frac*100:0.1f}%] {message}")
    
                processor.set_progress_fn(status)
            else:
                raise ValueError("progress should be True or a function")
    
>       processor.process()
E       RuntimeError: Unable to cast Python instance of type <class 'int'> to C++ type '?' (#define PYBIND11_DETAILED_ERROR_MESSAGES or compile in debug mode for details)