0
# Point Queries
1

2
Extract raster values at specific point locations or along vector geometry vertices. Point queries enable interpolated sampling of raster data at precise coordinates using nearest neighbor or bilinear interpolation methods.
3

4
## Capabilities
5

6
### Primary Functions
7

8
Query raster values at point locations with interpolation options.
9

10
```python { .api }
11
def point_query(vectors, raster, band=1, layer=0, nodata=None, affine=None,
12
                interpolate="bilinear", property_name="value",
13
                geojson_out=False, boundless=True):
14
    """
15
    Query raster values at point locations.
16
    
17
    Parameters:
18
    - vectors: Path to vector source or geo-like python objects
19
    - raster: Path to raster file or numpy array (requires affine if ndarray)
20
    - band: Raster band number, counting from 1 (default: 1)
21
    - layer: Vector layer to use, by name or number (default: 0)
22
    - nodata: Override raster nodata value (default: None)
23
    - affine: Affine transform, required for numpy arrays (default: None)
24
    - interpolate: Interpolation method - "bilinear" or "nearest" (default: "bilinear")
25
    - property_name: Property name for GeoJSON output (default: "value")
26
    - geojson_out: Return GeoJSON features with values as properties (default: False)
27
    - boundless: Allow points beyond raster extent (default: True)
28
    
29
    Returns:
30
    List of values (or single value for Point geometries)
31
    """
32

33
def gen_point_query(vectors, raster, band=1, layer=0, nodata=None, affine=None,
34
                    interpolate="bilinear", property_name="value",
35
                    geojson_out=False, boundless=True):
36
    """
37
    Generator version of point queries.
38
    
39
    Parameters are identical to point_query().
40
    
41
    Returns:
42
    Generator yielding values for each feature
43
    """
44
```
45

46
### Interpolation Functions
47

48
Core interpolation algorithms for point value extraction.
49

50
```python { .api }
51
def bilinear(arr, x, y):
52
    """
53
    Bilinear interpolation on 2x2 array using unit square coordinates.
54
    
55
    Parameters:
56
    - arr: 2x2 numpy array of values
57
    - x: X coordinate on unit square (0.0 to 1.0)
58
    - y: Y coordinate on unit square (0.0 to 1.0)
59
    
60
    Returns:
61
    Interpolated value or None if masked data present
62
    
63
    Raises:
64
    AssertionError if array not 2x2 or coordinates outside unit square
65
    """
66

67
def point_window_unitxy(x, y, affine):
68
    """
69
    Calculate 2x2 window and unit square coordinates for point interpolation.
70
    
71
    Parameters:
72
    - x: Point X coordinate in CRS units
73
    - y: Point Y coordinate in CRS units  
74
    - affine: Affine transformation from pixel to CRS coordinates
75
    
76
    Returns:
77
    Tuple of (window, unit_coordinates) where:
78
    - window: ((row_start, row_end), (col_start, col_end))
79
    - unit_coordinates: (unit_x, unit_y) on 0-1 square
80
    """
81
```
82

83
### Geometry Processing
84

85
Functions for extracting coordinates from vector geometries.
86

87
```python { .api }
88
def geom_xys(geom):
89
    """
90
    Generate flattened series of 2D points from Shapely geometry.
91
    Handles Point, LineString, Polygon, and Multi* geometries.
92
    Automatically converts 3D geometries to 2D.
93
    
94
    Parameters:
95
    - geom: Shapely geometry object
96
    
97
    Yields:
98
    (x, y) coordinate tuples for all vertices
99
    """
100
```
101

102
## Usage Examples
103

104
### Basic Point Queries
105

106
```python
107
from rasterstats import point_query
108

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

114
# Query multiple points
115
points = [
116
    {'type': 'Point', 'coordinates': (100, 200)},
117
    {'type': 'Point', 'coordinates': (150, 250)},
118
    {'type': 'Point', 'coordinates': (200, 300)}
119
]
120
values = point_query(points, "temperature.tif")
121
print(values)  # [15.2, 18.7, 22.1]
122
```
123

124
### Interpolation Methods
125

126
```python
127
# Bilinear interpolation (default) - smooth values between pixels
128
values_bilinear = point_query(points, "elevation.tif", interpolate="bilinear")
129

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

133
print(f"Bilinear: {values_bilinear[0]}")  # 74.098
134
print(f"Nearest: {values_nearest[0]}")   # 74.000
135
```
136

137
### Querying Line and Polygon Vertices
138

139
```python
140
# Query all vertices of a linestring (useful for elevation profiles)
141
line = {
142
    'type': 'LineString',
143
    'coordinates': [(100, 100), (150, 150), (200, 200), (250, 250)]
144
}
145
profile_values = point_query(line, "elevation.tif")
146
print(profile_values)  # [100.5, 125.3, 150.8, 175.2]
147

148
# Query polygon boundary vertices
149
polygon = {
150
    'type': 'Polygon',
151
    'coordinates': [[(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)]]
152
}
153
boundary_values = point_query(polygon, "temperature.tif")
154
print(boundary_values)  # [20.1, 22.5, 25.0, 23.2, 20.1]
155
```
156

157
### Working with NumPy Arrays
158

159
```python
160
import numpy as np
161
from affine import Affine
162

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

167
# Query points from array
168
points = [{'type': 'Point', 'coordinates': (25.5, 75.3)}]
169
values = point_query(points, raster_data, affine=transform,
170
                    interpolate="bilinear")
171
print(values)  # [23.456]
172
```
173

174
### GeoJSON Output
175

176
```python
177
# Return results as GeoJSON features with query values as properties
178
points_with_elevation = point_query(
179
    "sample_points.shp", 
180
    "dem.tif",
181
    geojson_out=True,
182
    property_name="elevation"
183
)
184

185
# Each feature retains original geometry and properties plus query result
186
print(points_with_elevation[0])
187
# {
188
#   'type': 'Feature',
189
#   'geometry': {'type': 'Point', 'coordinates': [100, 200]},
190
#   'properties': {'id': 1, 'name': 'Site A', 'elevation': 1250.5}
191
# }
192
```
193

194
### Multi-band Queries
195

196
```python
197
# Query different bands
198
red_values = point_query(points, "landsat.tif", band=1)    # Red band
199
green_values = point_query(points, "landsat.tif", band=2)  # Green band  
200
blue_values = point_query(points, "landsat.tif", band=3)   # Blue band
201

202
# Combine results
203
rgb_values = list(zip(red_values, green_values, blue_values))
204
print(rgb_values[0])  # (145, 167, 98)
205
```
206

207
### Handling Edge Cases
208

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

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

220
# Override nodata value
221
value = point_query(nodata_point, "raster.tif", nodata=-9999)
222
```
223

224
### Generator Usage for Large Datasets
225

226
```python
227
from rasterstats import gen_point_query
228

229
# Process large point datasets efficiently
230
large_points = "million_points.shp"
231
elevation_raster = "high_resolution_dem.tif"
232

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

244
### Advanced Coordinate Systems
245

246
```python
247
# Query points in different coordinate systems
248
# (raster and vector CRS should match or be properly transformed beforehand)
249

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

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

259
## Types
260

261
```python { .api }
262
from typing import Union, List, Optional, Tuple, Generator, Any
263
from numpy import ndarray
264
from affine import Affine
265
from shapely.geometry.base import BaseGeometry
266

267
InterpolationMethod = Union["bilinear", "nearest"]
268
PointValue = Union[float, int, None]
269
PointResult = Union[PointValue, List[PointValue]]
270
Window = Tuple[Tuple[int, int], Tuple[int, int]]
271
UnitCoordinates = Tuple[float, float]
272
```