Skills › Scientific Agent Skills › geomaster › raster_workflows.py
raster_workflows.py
153 lines · skills/geomaster/scripts/raster_workflows.py at 5f2bc03
1"""Small, in-memory geospatial recipes; caller supplies scientific provenance.
3Import from the skill's scripts directory. No downloads or model fitting occur on import.
4"""
5from pathlib import Path
6import numpy as np
7import rasterio
8from rasterio.features import rasterize
11def as_float(values):
12 """Preserve masked cells and nonfinite values as NaN before arithmetic."""
13 return np.ma.asarray(values, dtype=np.float64).filled(np.nan)
16def ratio(numerator, denominator):
17 numerator, denominator = np.broadcast_arrays(as_float(numerator), as_float(denominator))
18 result = np.full(numerator.shape, np.nan)
19 valid = np.isfinite(numerator) & np.isfinite(denominator) & (denominator != 0)
20 np.divide(numerator, denominator, out=result, where=valid)
21 return result
24def normalized_difference(a, b):
25 a, b = as_float(a), as_float(b)
26 if a.shape != b.shape:
27 raise ValueError('Bands must have identical shapes on a verified common grid')
28 return ratio(a - b, a + b)
31def spectral_indices(blue, green, red, nir, swir1, swir2):
32 """Input is aligned, masked physical reflectance, NOT unscaled digital numbers."""
33 bands = [as_float(x) for x in (blue, green, red, nir, swir1, swir2)]
34 if len({x.shape for x in bands}) != 1:
35 raise ValueError('Bands must have identical shapes')
36 blue, green, red, nir, swir1, swir2 = bands
37 return {
38 'NDVI': normalized_difference(nir, red),
39 'EVI': 2.5 * ratio(nir - red, nir + 6 * red - 7.5 * blue + 1),
40 'SAVI': 1.5 * ratio(nir - red, nir + red + 0.5),
41 'NDWI': normalized_difference(green, nir),
42 'MNDWI': normalized_difference(green, swir1),
43 'NBR': normalized_difference(nir, swir2),
44 'NDBI': normalized_difference(swir1, nir),
45 }
48def terrain_metrics(dem, transform, crs, *, azimuth=315, altitude=45):
49 """Finite-difference slope/aspect/hillshade for north-up metric DEMs.
51 Elevations MUST be metres. Aspect is downslope, clockwise from north; flat
52 pixels have undefined (NaN) aspect. This is not a hydrological flow algorithm.
53 """
54 from pyproj import CRS
55 crs = CRS.from_user_input(crs)
56 if not crs.is_projected or any(abs(a.unit_conversion_factor - 1) > 1e-12 for a in crs.axis_info[:2]):
57 raise ValueError('A projected CRS with metre horizontal units is required')
58 if not np.isfinite(tuple(transform)).all() or transform.b != 0 or transform.d != 0 or transform.a <= 0 or transform.e >= 0:
59 raise ValueError('A finite north-up, unrotated transform is required')
60 if not np.isfinite([azimuth, altitude]).all() or not 0 <= altitude <= 90:
61 raise ValueError('Finite azimuth and altitude in [0, 90] required')
62 z = as_float(dem)
63 if z.ndim != 2 or min(z.shape) < 2:
64 raise ValueError('DEM must be a two-dimensional array with at least 2 cells per axis')
65 dz_dy, dz_dx = np.gradient(z, transform.e, transform.a)
66 norm = np.hypot(dz_dx, dz_dy)
67 slope = np.arctan(norm)
68 aspect = np.mod(np.arctan2(-dz_dx, -dz_dy), 2 * np.pi)
69 az, alt = np.deg2rad([azimuth, altitude])
70 shade = np.sin(alt) * np.cos(slope) + np.cos(alt) * np.sin(slope) * np.cos(az - aspect)
71 valid = np.isfinite(z) & np.isfinite(norm)
72 return (np.where(valid, np.rad2deg(slope), np.nan),
73 np.where(valid & (norm > 0), np.rad2deg(aspect), np.nan),
74 np.where(valid, np.clip(shade, 0, 1), np.nan))
77def write_ndvi(input_path, output_path, *, red_band, nir_band, scale, offset):
78 """Verified stack band indices and provider-specific common scale/offset required.
80 Input mask must already incorporate cloud/shadow QA. Rasterio does not infer
81 Sentinel band positions. A supplied scale/offset is applied exactly once.
82 """
83 if Path(output_path).exists() or Path(output_path).is_symlink():
84 raise ValueError('Output already exists; select a new path')
85 if red_band == nir_band:
86 raise ValueError('Red and NIR must be distinct bands')
87 if Path(input_path).resolve() == Path(output_path).resolve():
88 raise ValueError('Input and output paths must differ')
89 if not np.isfinite([scale, offset]).all() or scale <= 0:
90 raise ValueError('Scale must be positive and scale/offset must be finite')
91 with rasterio.open(input_path) as src:
92 if src.crs is None:
93 raise ValueError('Input CRS is missing')
94 for band in (red_band, nir_band):
95 if isinstance(band, bool) or not isinstance(band, int) or not 1 <= band <= src.count:
96 raise ValueError('Band positions are one-based integers within the stack')
97 red = src.read(red_band, masked=True).astype('float64') * scale + offset
98 nir = src.read(nir_band, masked=True).astype('float64') * scale + offset
99 result = normalized_difference(nir, red).astype('float32')
100 profile = src.profile.copy()
101 profile.update(driver='GTiff', count=1, dtype='float32', nodata=np.nan)
102 with rasterio.open(output_path, 'w', **profile) as dst:
103 dst.write(result, 1)
104 dst.set_band_description(1, 'NDVI')
105 return result
108def classify_imagery(raster_path, training_gdf, output_path, *, random_state=42):
109 """Small in-memory RF fit/predict recipe. Accuracy must be assessed separately.
111 Reject overlapping training pixel centres and require labels 1..65534; 0 is output
112 nodata. Align labels to raster CRS, preserve masks, and return the fitted model.
113 """
114 from sklearn.ensemble import RandomForestClassifier
115 if Path(output_path).exists() or Path(output_path).is_symlink():
116 raise ValueError('Output already exists; select a new path')
117 if Path(raster_path).resolve() == Path(output_path).resolve():
118 raise ValueError('Input and output paths must differ')
119 with rasterio.open(raster_path) as src:
120 if src.crs is None or training_gdf.crs is None:
121 raise ValueError('Both raster and training data need known CRS')
122 if training_gdf.empty or training_gdf.geometry.isna().any() or training_gdf.geometry.is_empty.any() or not training_gdf.is_valid.all():
123 raise ValueError('Training geometries must be nonempty and valid')
124 if not training_gdf.geom_type.isin(['Polygon', 'MultiPolygon']).all():
125 raise ValueError('Training geometries must be polygons')
126 labels = np.asarray(training_gdf['class_id'], dtype=float)
127 if not np.isfinite(labels).all() or np.any(labels != np.floor(labels)) or np.any((labels < 1) | (labels > 65534)):
128 raise ValueError('class_id must be an integer from 1 to 65534')
129 training = training_gdf.to_crs(src.crs)
130 image = as_float(src.read(masked=True))
131 valid = np.isfinite(image).all(axis=0)
132 targets = np.zeros(image.shape[1:], dtype='uint16')
133 occupied = np.zeros_like(targets, dtype=bool)
134 for geom, label in zip(training.geometry, labels):
135 cells = rasterize([(geom, 1)], out_shape=targets.shape, transform=src.transform, dtype='uint8').astype(bool)
136 if not np.any(cells & valid):
137 raise ValueError('A training polygon contains no valid pixel centres')
138 if np.any(occupied & cells):
139 raise ValueError('Training polygons overlap at raster pixel centres')
140 occupied |= cells
141 targets[cells] = int(label)
142 train_mask = occupied & valid
143 if np.unique(targets[train_mask]).size < 2:
144 raise ValueError('At least two observed classes are required')
145 model = RandomForestClassifier(n_estimators=100, random_state=random_state, n_jobs=1)
146 model.fit(image[:, train_mask].T, targets[train_mask])
147 prediction = np.zeros(targets.shape, dtype='uint16')
148 prediction[valid] = model.predict(image[:, valid].T)
149 profile = src.profile.copy()
150 profile.update(driver='GTiff', dtype='uint16', nodata=0, count=1)
151 with rasterio.open(output_path, 'w', **profile) as dst:
152 dst.write(prediction, 1)
153 return model