Point Queries

def point_query(vectors, raster, band=1, layer=0, nodata=None, affine=None,
                interpolate="bilinear", property_name="value",
                geojson_out=False, boundless=True):
    """
    Query raster values at point locations.
    
    Parameters:
    - vectors: Path to vector source or geo-like python objects
    - raster: Path to raster file or numpy array (requires affine if ndarray)
    - band: Raster band number, counting from 1 (default: 1)
    - layer: Vector layer to use, by name or number (default: 0)
    - nodata: Override raster nodata value (default: None)
    - affine: Affine transform, required for numpy arrays (default: None)
    - interpolate: Interpolation method - "bilinear" or "nearest" (default: "bilinear")
    - property_name: Property name for GeoJSON output (default: "value")
    - geojson_out: Return GeoJSON features with values as properties (default: False)
    - boundless: Allow points beyond raster extent (default: True)
    
    Returns:
    List of values (or single value for Point geometries)
    """

def gen_point_query(vectors, raster, band=1, layer=0, nodata=None, affine=None,
                    interpolate="bilinear", property_name="value",
                    geojson_out=False, boundless=True):
    """
    Generator version of point queries.
    
    Parameters are identical to point_query().
    
    Returns:
    Generator yielding values for each feature
    """

def bilinear(arr, x, y):
    """
    Bilinear interpolation on 2x2 array using unit square coordinates.
    
    Parameters:
    - arr: 2x2 numpy array of values
    - x: X coordinate on unit square (0.0 to 1.0)
    - y: Y coordinate on unit square (0.0 to 1.0)
    
    Returns:
    Interpolated value or None if masked data present
    
    Raises:
    AssertionError if array not 2x2 or coordinates outside unit square
    """

def point_window_unitxy(x, y, affine):
    """
    Calculate 2x2 window and unit square coordinates for point interpolation.
    
    Parameters:
    - x: Point X coordinate in CRS units
    - y: Point Y coordinate in CRS units  
    - affine: Affine transformation from pixel to CRS coordinates
    
    Returns:
    Tuple of (window, unit_coordinates) where:
    - window: ((row_start, row_end), (col_start, col_end))
    - unit_coordinates: (unit_x, unit_y) on 0-1 square
    """

def geom_xys(geom):
    """
    Generate flattened series of 2D points from Shapely geometry.
    Handles Point, LineString, Polygon, and Multi* geometries.
    Automatically converts 3D geometries to 2D.
    
    Parameters:
    - geom: Shapely geometry object
    
    Yields:
    (x, y) coordinate tuples for all vertices
    """

from rasterstats import point_query

# Query single point
point = {'type': 'Point', 'coordinates': (245309.0, 1000064.0)}
value = point_query(point, "elevation.tif")
print(value)  # [74.09817594635244]

# Query multiple points
points = [
    {'type': 'Point', 'coordinates': (100, 200)},
    {'type': 'Point', 'coordinates': (150, 250)},
    {'type': 'Point', 'coordinates': (200, 300)}
]
values = point_query(points, "temperature.tif")
print(values)  # [15.2, 18.7, 22.1]

# Bilinear interpolation (default) - smooth values between pixels
values_bilinear = point_query(points, "elevation.tif", interpolate="bilinear")

# Nearest neighbor - use exact pixel values
values_nearest = point_query(points, "elevation.tif", interpolate="nearest")

print(f"Bilinear: {values_bilinear[0]}")  # 74.098
print(f"Nearest: {values_nearest[0]}")   # 74.000

# Query all vertices of a linestring (useful for elevation profiles)
line = {
    'type': 'LineString',
    'coordinates': [(100, 100), (150, 150), (200, 200), (250, 250)]
}
profile_values = point_query(line, "elevation.tif")
print(profile_values)  # [100.5, 125.3, 150.8, 175.2]

# Query polygon boundary vertices
polygon = {
    'type': 'Polygon',
    'coordinates': [[(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)]]
}
boundary_values = point_query(polygon, "temperature.tif")
print(boundary_values)  # [20.1, 22.5, 25.0, 23.2, 20.1]

import numpy as np
from affine import Affine

# Create sample raster data
raster_data = np.random.rand(100, 100) * 50  # 0-50 range
transform = Affine.translation(0, 100) * Affine.scale(1, -1)

# Query points from array
points = [{'type': 'Point', 'coordinates': (25.5, 75.3)}]
values = point_query(points, raster_data, affine=transform,
                    interpolate="bilinear")
print(values)  # [23.456]

# Return results as GeoJSON features with query values as properties
points_with_elevation = point_query(
    "sample_points.shp", 
    "dem.tif",
    geojson_out=True,
    property_name="elevation"
)

# Each feature retains original geometry and properties plus query result
print(points_with_elevation[0])
# {
#   'type': 'Feature',
#   'geometry': {'type': 'Point', 'coordinates': [100, 200]},
#   'properties': {'id': 1, 'name': 'Site A', 'elevation': 1250.5}
# }

# Query different bands
red_values = point_query(points, "landsat.tif", band=1)    # Red band
green_values = point_query(points, "landsat.tif", band=2)  # Green band  
blue_values = point_query(points, "landsat.tif", band=3)   # Blue band

# Combine results
rgb_values = list(zip(red_values, green_values, blue_values))
print(rgb_values[0])  # (145, 167, 98)

# Points outside raster extent (with boundless=True)
distant_point = {'type': 'Point', 'coordinates': (999999, 999999)}
value = point_query(distant_point, "small_raster.tif", boundless=True)
print(value)  # [None] - outside extent

# Points on nodata areas
nodata_point = {'type': 'Point', 'coordinates': (100, 100)}
value = point_query(nodata_point, "raster_with_holes.tif")
print(value)  # [None] - nodata pixel

# Override nodata value
value = point_query(nodata_point, "raster.tif", nodata=-9999)

from rasterstats import gen_point_query

# Process large point datasets efficiently
large_points = "million_points.shp"
elevation_raster = "high_resolution_dem.tif"

# Generator avoids loading all results into memory
for i, value in enumerate(gen_point_query(large_points, elevation_raster)):
    if i % 10000 == 0:
        print(f"Processed {i} points, current value: {value}")
    
    # Process value immediately
    if value and value[0] > 1000:  # Points above 1000m elevation
        # Do something with high elevation points
        pass

# Query points in different coordinate systems
# (raster and vector CRS should match or be properly transformed beforehand)

# UTM coordinates
utm_points = [{'type': 'Point', 'coordinates': (500000, 4500000)}]
values = point_query(utm_points, "utm_raster.tif")

# Geographic coordinates (lat/lon)
latlon_points = [{'type': 'Point', 'coordinates': (-122.4194, 37.7749)}]
values = point_query(latlon_points, "geographic_raster.tif")

from typing import Union, List, Optional, Tuple, Generator, Any
from numpy import ndarray
from affine import Affine
from shapely.geometry.base import BaseGeometry

InterpolationMethod = Union["bilinear", "nearest"]
PointValue = Union[float, int, None]
PointResult = Union[PointValue, List[PointValue]]
Window = Tuple[Tuple[int, int], Tuple[int, int]]
UnitCoordinates = Tuple[float, float]