Sentinel-2 NDVI Zonal Statistics
Domain: agriculture_env — Precision Agriculture

Summary: Compute a Sentinel-2 NDVI GeoTIFF and polygon-level NDVI zonal statistics for agricultural units with a CRS mismatch between vector and raster inputs.

Task prompt:
You are given a polygon shapefile of agricultural units and a 6-band Sentinel-2 GeoTIFF.

Task directory:
`base`

Inputs:
- `base\input\rotation_count.shp`
  - associated `.dbf`, `.shx`, `.prj`, and `.cpg` files are present in the same directory
  - each feature has `id_lcp`, `rot_count`, and polygon geometry
- `base\input\s2_l2a_30m_clip.tif`

Band numbering is 1-based:
- band 1 = B2
- band 2 = B3
- band 3 = B4 (Red)
- band 4 = B8 (NIR)
- band 5 = B11
- band 6 = B12

The shapefile and raster use different coordinate reference systems.

Optional Python runtime:
- Use `base\software\python_geo.bat` if you want the benchmark-provisioned geospatial Python stack.
- The package manifest is visible at `base\input\runtime_env\pyproject.toml`.

Required outputs under `base\output`:
- `ndvi.tif`
- `polygon_ndvi_stats.csv`

Required logic:
1. Load the input shapefile. Use `id_lcp` as the polygon identifier and `rot_count` exactly as provided. Do not recompute or modify crop rotation counts.
2. Load the input Sentinel-2 raster.
3. Compute NDVI using `(B8 - B4) / (B8 + B4)`, where B4 is raster band 3 and B8 is raster band 4.
4. Use the raster values exactly as provided. Do not apply reflectance scaling, offsets, normalization, clipping, cloud masking, or any other masking rule.
5. A pixel is invalid only if band 3 is non-finite, band 4 is non-finite, or `(B8 + B4) == 0`.
6. Write `ndvi.tif` as a single-band `float32` GeoTIFF with the same width, height, CRS, transform, extent, and pixel alignment as the input raster.
7. Store invalid pixels in `ndvi.tif` as `NaN` values in the raster array.
8. Do not reproject or resample the raster.
9. Reproject polygon geometries to the raster CRS before computing zonal statistics. Do not dissolve, explode, simplify, buffer, repair, or otherwise modify the input polygon geometries.
10. Use the pixel-center rule (`all_touched=False`) for polygon inclusion.
11. For every polygon, compute `valid_px`, `mean_ndvi`, and `median_ndvi` from finite NDVI pixels.
12. Include every polygon exactly once. If a polygon has no finite NDVI pixels, write `valid_px` as `0`, `mean_ndvi` as `NA`, and `median_ndvi` as `NA`.
13. Write `polygon_ndvi_stats.csv` with exactly these columns, in order: `id_lcp`, `rot_count`, `valid_px`, `mean_ndvi`, `median_ndvi`.
14. Sort the CSV by `id_lcp` in ascending lexicographic order. Treat `id_lcp` as a string identifier.
15. Compute statistics from the final `float32` NDVI raster values, then round `mean_ndvi` and `median_ndvi` to 6 decimal places only when writing the CSV.

Rules:
- Use only the provided local files.
- Do not modify input files.
- Do not change required output filenames.
- Do not ask for confirmation; execute directly.


Required actions:
Use the staged shapefile `input/rotation_count.shp` and Sentinel-2 raster `input/s2_l2a_30m_clip.tif`.
Compute NDVI from raster band 4 and band 3 without extra scaling or masking.
Write `output/ndvi.tif` as single-band float32 with the same georeferencing as the input raster.
Reproject polygons to the raster CRS for zonal statistics, using pixel-center inclusion.
Write `output/polygon_ndvi_stats.csv` with the required five columns and one row per polygon.

Software: Python
rasterio
geopandas

Input files:
name: rotation_count.shp
format: ESRI Shapefile
description: Agricultural unit polygons with `id_lcp` and precomputed `rot_count` fields
name: s2_l2a_30m_clip.tif
format: GeoTIFF
description: Six-band Sentinel-2 raster; band 3 is B4/Red and band 4 is B8/NIR
name: runtime_env/pyproject.toml
format: TOML
description: Optional dependency manifest for a Python geospatial solution