# Zonal Statistics Generation

The `gdptools.zonal_gen` module provides powerful tools for computing zonal statistics on raster data. This module supports both standard and area-weighted zonal statistics calculations, with multiple computational engines for scalability and performance optimization.

## Features

- **Multiple Computational Engines**: Choose between serial, parallel, Dask, and exactextract engines based on your data size and computational resources
- **High-Performance Option**: The `exactextract` engine provides optimized C++ performance with fractional pixel coverage
- **Area-Weighted Statistics**: Calculate more accurate statistics by accounting for partial pixel coverage within zones
- **Flexible Output Formats**: Export results to CSV with customizable precision and file naming
- **Categorical and Continuous Data Support**: Handle both categorical (e.g., land cover) and continuous (e.g., temperature) raster data
- **Performance Monitoring**: Built-in execution time tracking for performance optimization

## Quick Start

### Standard Zonal Statistics

```python
from gdptools.zonal_gen import ZonalGen
from gdptools.data.user_data import UserData

# Initialize user data (preprocessed raster and vector data)
user_data = UserData(...)

# Create zonal statistics generator
zonal_gen = ZonalGen(
    user_data=user_data,
    zonal_engine="parallel",
    zonal_writer="csv",
    out_path="/path/to/output",
    file_prefix="zonal_stats",
    jobs=4
)

# Calculate zonal statistics
stats = zonal_gen.calculate_zonal()
print(stats.head())
```

### Area-Weighted Zonal Statistics

```{tip}
For `WeightedZonalGen`, use an equal-area CRS for accurate area weighting. The project’s recommended default is `EPSG:6931` (EASE-Grid 2.0 Global Equal Area), which performs well globally. For regional work, consider a local equal-area CRS (e.g., `EPSG:5070` for CONUS, `EPSG:3035` for Europe).
```

```python
from gdptools.zonal_gen import WeightedZonalGen
from pyproj import CRS

# Use an appropriate equal-area projection
albers_crs = CRS.from_proj4("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96")

# Create weighted zonal statistics generator
weighted_gen = WeightedZonalGen(
    user_data=user_data,
    weight_gen_crs=albers_crs,
    zonal_engine="dask",
    zonal_writer="csv",
    out_path="/path/to/output",
    file_prefix="weighted_stats",
    precision=3,
    jobs=8
)

# Calculate area-weighted zonal statistics
weighted_stats = weighted_gen.calculate_weighted_zonal()
```

## Computational Engines

### Serial Engine

Best for small datasets or when memory is limited:

```python
zonal_gen = ZonalGen(
    user_data=user_data,
    zonal_engine="serial",
    zonal_writer="csv"
)
```

### Parallel Engine

Optimal for medium to large datasets on multi-core systems:

```python
zonal_gen = ZonalGen(
    user_data=user_data,
    zonal_engine="parallel",
    zonal_writer="csv",
    jobs=4  # Use 4 CPU cores
)
```

### Dask Engine

Best for very large datasets or distributed computing:

```python
zonal_gen = ZonalGen(
    user_data=user_data,
    zonal_engine="dask",
    zonal_writer="csv",
    jobs=8  # Use 8 workers
)
```

### Exactextract Engine

Best for high-performance zonal statistics with fractional pixel coverage:

```python
zonal_gen = ZonalGen(
    user_data=user_data,
    zonal_engine="exactextract",
    zonal_writer="csv",
    out_path="/path/to/output",
    file_prefix="zonal_stats"
)
```

The `exactextract` engine uses the [exactextract](https://github.com/isciences/exactextract) library,
which provides high-performance zonal statistics through optimized C++ code. Key benefits include:

- **Fractional pixel coverage**: Accounts for partial pixel overlap with polygon boundaries
- **Memory efficiency**: Processes data without loading entire rasters into memory
- **Speed**: Often significantly faster than pure Python implementations for large datasets

```{note}
The `exactextract` engine requires the `exactextract` Python package (`pip install exactextract`).
It produces output in the same format as the serial/parallel/dask engines, but values may differ
slightly due to the fractional coverage algorithm vs. whole-pixel counting.
```

## Data Types and Statistics

### Continuous Data

For continuous raster data (temperature, precipitation, elevation):

```python
stats = zonal_gen.calculate_zonal(categorical=False)
# Returns: count, mean, std, min, 25%, 50%, 75%, max, sum
```

### Categorical Data

For categorical raster data (land cover, soil types):

```python
stats = zonal_gen.calculate_zonal(categorical=True)
# Returns: fraction columns for each category + count
```

## Best Practices

### Coordinate Reference Systems

- **For Area-Weighted Statistics**: Use equal-area projections appropriate for your study region
- **Geographic Coordinates**: Can introduce area distortions; consider reprojecting to equal-area
- **Local Projections**: Often provide the most accurate area calculations

```python
# Good: Equal-area projection for continental US
albers_crs = CRS.from_proj4("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96")

# Acceptable: UTM zone for smaller regions
utm_crs = CRS.from_epsg(32618)  # UTM Zone 18N

# Caution: Geographic coordinates may introduce area distortions
geographic_crs = CRS.from_epsg(4326)
```

### Performance Optimization

1. **Choose the Right Engine**: Start with parallel for most use cases
2. **Tune Job Count**: Typically use number of CPU cores minus 1. If you request
   more workers than available CPUs, the zonal engines automatically cap the
   value and emit a warning so you can adjust expectations.
3. **Memory Management**: Monitor memory usage with large datasets
4. **Data Preprocessing**: Ensure data is properly aligned and preprocessed

### Output Management

```python
# Organized output with timestamps
zonal_gen = ZonalGen(
    user_data=user_data,
    zonal_engine="parallel",
    zonal_writer="csv",
    out_path="/results/zonal_stats",
    file_prefix="temperature_zones",
    append_date=True,  # Creates: 2024_01_15_14_30_00_temperature_zones.csv
    precision=2  # Round to 2 decimal places
)
```

## Error Handling

### Common Issues and Solutions

**FileNotFoundError**: Output path doesn't exist

```python
from pathlib import Path
output_path = Path("/path/to/output")
output_path.mkdir(parents=True, exist_ok=True)
```

**Memory Issues**: Dataset too large for available memory

```python
# Switch to Dask engine for out-of-core processing
zonal_gen = ZonalGen(
    user_data=user_data,
    zonal_engine="dask",
    zonal_writer="csv",
    jobs=4
)
```

**CRS Mismatch**: Ensure data alignment

```python
# Verify CRS compatibility before processing
print(f"Raster CRS: {user_data.raster_crs}")
print(f"Vector CRS: {user_data.vector_crs}")
print(f"Weight CRS: {weight_gen_crs}")
```

## Advanced Examples

### Multi-Temporal Analysis

```python
import pandas as pd
from pathlib import Path

# Process multiple time periods
results = []
for year in range(2010, 2021):
    user_data = UserData(raster_path=f"/data/temp_{year}.tif", ...)

    zonal_gen = ZonalGen(
        user_data=user_data,
        zonal_engine="parallel",
        zonal_writer="csv",
        out_path=f"/results/{year}",
        file_prefix=f"temp_stats_{year}",
        jobs=4
    )

    stats = zonal_gen.calculate_zonal()
    stats['year'] = year
    results.append(stats)

# Combine all years
all_stats = pd.concat(results, ignore_index=True)
```

### Comparison of Standard vs. Weighted Statistics

```python
# Calculate both standard and weighted statistics
standard_stats = zonal_gen.calculate_zonal()
weighted_stats = weighted_gen.calculate_weighted_zonal()

# Compare results
comparison = pd.merge(
    standard_stats[['zone_id', 'mean']],
    weighted_stats[['zone_id', 'weighted_mean']],
    on='zone_id',
    suffixes=('_standard', '_weighted')
)

comparison['difference'] = (
    comparison['weighted_mean'] - comparison['mean']
)
```

## Type Definitions

```{eval-rst}
.. autodata:: gdptools.zonal_gen.ZONAL_ENGINES
   :annotation: = Literal["serial", "parallel", "dask"]

   Computational engine options for zonal statistics calculation.

.. autodata:: gdptools.zonal_gen.ZONAL_WRITERS
   :annotation: = Literal["csv"]

   Output format options for zonal statistics results.
```

## API Reference

### Classes

```{eval-rst}
.. autoclass:: gdptools.zonal_gen.ZonalGen
   :members:
   :show-inheritance:

   Standard zonal statistics calculator for raster data.

.. autoclass:: gdptools.zonal_gen.WeightedZonalGen
   :members:
   :show-inheritance:

   Area-weighted zonal statistics calculator for raster data.
```

## See Also

- {doc}`weight_gen_classes`: Weight generation for raster-vector operations
- {doc}`agg_gen_classes`: Aggregation methods for processed weights
- {doc}`user_input_data_classes`: Data preparation and validation tools

---

_For questions or contributions, please refer to the project's GitHub repository or contact the development team._