Command Line Interface

@click.command(context_settings=SETTINGS)
@cligj.features_in_arg
@click.version_option(version=version, message="%(version)s")
@click.option("--raster", "-r", required=True)
@click.option("--all-touched/--no-all-touched", default=False)
@click.option("--band", type=int, default=1)
@click.option("--categorical/--no-categorical", default=False)
@click.option("--indent", type=int, default=None)
@click.option("--info/--no-info", default=False)
@click.option("--nodata", type=int, default=None)
@click.option("--prefix", type=str, default="_")
@click.option("--stats", type=str, default=None)
@click.option("--sequence/--no-sequence", type=bool, default=False)
@cligj.use_rs_opt
def zonalstats(features, raster, all_touched, band, categorical, indent, 
               info, nodata, prefix, stats, sequence, use_rs):
    """
    Generate summary statistics of geospatial raster datasets based on vector features.
    
    Reads GeoJSON Features and outputs GeoJSON with additional statistical properties.
    The raster is specified by the required -r/--raster argument.
    
    Parameters:
    - features: Input GeoJSON features (from stdin or file)
    - raster: Path to raster file (-r/--raster, required)
    - all_touched: Include all touched pixels vs center-point only
    - band: Raster band number (default: 1)
    - categorical: Enable categorical analysis mode
    - indent: JSON indentation for output formatting
    - info: Enable info-level logging
    - nodata: Override raster nodata value
    - prefix: Prefix for statistic property names (default: "_")
    - stats: Space-delimited statistics to calculate or "all"
    - sequence: Output as GeoJSON sequence vs FeatureCollection
    - use_rs: Use record separator for sequence output
    """

@click.command(context_settings=SETTINGS)
@cligj.features_in_arg
@click.version_option(version=version, message="%(version)s")
@click.option("--raster", "-r", required=True)
@click.option("--band", type=int, default=1)
@click.option("--nodata", type=int, default=None)
@click.option("--indent", type=int, default=None)
@click.option("--interpolate", type=str, default="bilinear")
@click.option("--property-name", type=str, default="value")
@click.option("--sequence/--no-sequence", type=bool, default=False)
@cligj.use_rs_opt
def pointquery(features, raster, band, indent, nodata, interpolate,
               property_name, sequence, use_rs):
    """
    Query raster values at points of input GeoJSON Features.
    
    Adds raster values to feature properties and outputs as GeoJSON.
    For Points, uses point geometry. For other types, queries all vertices.
    
    Parameters:
    - features: Input GeoJSON features (from stdin or file)
    - raster: Path to raster file (-r/--raster, required)
    - band: Raster band number (default: 1)
    - indent: JSON indentation for output formatting
    - nodata: Override raster nodata value
    - interpolate: Interpolation method - "bilinear" or "nearest" (default: "bilinear")
    - property_name: Property name for query results (default: "value")
    - sequence: Output as GeoJSON sequence vs FeatureCollection
    - use_rs: Use record separator for sequence output
    """

SETTINGS = {"help_option_names": ["-h", "--help"]}  # Click command settings
version: str  # Version string from rasterstats.__version__

# Basic zonal stats from shapefile to GeoJSON
fio cat watersheds.shp | rio zonalstats -r elevation.tif > watershed_elevation.geojson

# Calculate specific statistics
fio cat polygons.shp | rio zonalstats -r temperature.tif --stats "mean max min std" > temp_stats.geojson

# All available statistics
fio cat zones.shp | rio zonalstats -r data.tif --stats "all" > comprehensive_stats.geojson

# Categorical analysis with prefix
fio cat counties.shp | rio zonalstats -r landcover.tif --categorical --prefix "lc_" > landcover_stats.geojson

# All-touched rasterization
fio cat thin_polygons.shp | rio zonalstats -r raster.tif --all-touched > touched_stats.geojson

# Specific band and nodata override
fio cat areas.shp | rio zonalstats -r multiband.tif --band 3 --nodata -999 > band3_stats.geojson

# Basic point queries
fio cat points.shp | rio pointquery -r elevation.tif > points_with_elevation.geojson

# Nearest neighbor interpolation
fio cat sample_points.shp | rio pointquery -r categorical_raster.tif --interpolate nearest > point_categories.geojson

# Custom property name
fio cat weather_stations.shp | rio pointquery -r temperature.tif --property-name "temp_c" > stations_temp.geojson

# Query specific band
fio cat locations.shp | rio pointquery -r landsat.tif --band 4 --property-name "nir" > nir_values.geojson

# Line profile (queries all vertices)
fio cat transect_line.shp | rio pointquery -r elevation.tif --property-name "elevation" > elevation_profile.geojson

# Formatted JSON output
fio cat polygons.shp | rio zonalstats -r data.tif --indent 2 > formatted_output.geojson

# Sequence output (one feature per line)
fio cat large_dataset.shp | rio zonalstats -r raster.tif --sequence > streaming_output.geojson

# Record separator for streaming
fio cat huge_dataset.shp | rio zonalstats -r raster.tif --sequence --use-rs | head -n 100

# Chain with other rasterio/fiona commands
fio cat input.shp | \
  rio zonalstats -r dem.tif --stats "mean" --prefix "elev_" | \
  rio pointquery -r temperature.tif --property-name "temp" | \
  fio collect > combined_analysis.geojson

# Filter results with jq
fio cat watersheds.shp | \
  rio zonalstats -r precipitation.tif --stats "sum" --prefix "precip_" | \
  jq '.features[] | select(.properties.precip_sum > 1000)' > high_precip_watersheds.geojson

# Convert to other formats
fio cat regions.shp | \
  rio zonalstats -r population.tif --stats "sum count" --prefix "pop_" | \
  fio load | \
  fio dump --format CSV > population_by_region.csv

# Enable info logging
fio cat polygons.shp | rio zonalstats -r raster.tif --info > stats_with_logging.geojson

# Check version
rio zonalstats --version
rio pointquery --version

# Help information
rio zonalstats --help
rio pointquery --help

# Direct GeoJSON input
echo '{"type":"Point","coordinates":[100,200]}' | \
  rio pointquery -r elevation.tif

# Process from URL
curl -s "https://example.com/data.geojson" | \
  rio zonalstats -r remote_raster.tif

# Pipe through multiple processing steps
fio cat raw_data.shp | \
  rio zonalstats -r raster1.tif --prefix "r1_" | \
  rio zonalstats -r raster2.tif --prefix "r2_" | \
  fio dump --format Shapefile output_with_stats.shp

# Stream processing of large datasets
fio cat large_polygons.shp | \
  rio zonalstats -r high_res_raster.tif --sequence | \
  while IFS= read -r feature; do
    # Process each feature individually
    echo "$feature" | jq '.properties.mean' >> means.txt
  done

# Parallel processing (GNU parallel)
fio cat huge_dataset.shp | \
  parallel --pipe --block 10M \
    "rio zonalstats -r raster.tif --sequence" > \
  parallel_results.geojson

# Install rasterstats
pip install rasterstats

# Verify rasterio plugins are registered
rio --help
# Should show 'zonalstats' and 'pointquery' in the commands list

# Check plugin versions
rio zonalstats --version
rio pointquery --version

# From pyproject.toml
entry_points = {
    "rasterio.rio_plugins": {
        "zonalstats": "rasterstats.cli:zonalstats",
        "pointquery": "rasterstats.cli:pointquery"
    }
}

# Missing raster file
fio cat points.shp | rio pointquery -r nonexistent.tif
# Error: Unable to open raster file

# Invalid statistics
fio cat polygons.shp | rio zonalstats -r raster.tif --stats "invalid_stat"
# Error: Stat 'invalid_stat' not valid

# Coordinate system mismatch (warning, but processed)
fio cat utm_polygons.shp | rio zonalstats -r geographic_raster.tif
# Warning: Coordinate systems may not match

# Invalid interpolation method
fio cat points.shp | rio pointquery -r raster.tif --interpolate "invalid"
# Error: interpolate must be nearest or bilinear