import argparseimport mathimport osimport numpy as npimport rasteriofrom rasterio.enums import Resamplingfrom rasterio.warp import reprojectdef parse_args(): parser = argparse.ArgumentParser( description="Compute TVDI / CWSI from multispectral and thermal raster inputs." ) parser.add_argument( "--multispectral", default=r"D:\data\MultiSpectral\0707\result1.tif", help="Path to multispectral input raster with Red and NIR bands.", ) parser.add_argument( "--thermal", default=r"D:\data\MultiSpectral\0708s\TIR_LST20260708.tif", help="Path to thermal LST input raster.", ) parser.add_argument( "--output-dir", default=r"D:\data\MultiSpectral\output", help="Directory where TVDI and drought class outputs will be written.", ) parser.add_argument( "--nir-band", type=int, default=3, help="Band index for NIR in the multispectral file (1-based).", ) parser.add_argument( "--red-band", type=int, default=2, help="Band index for Red in the multispectral file (1-based).", ) parser.add_argument( "--min-ndvi", type=float, default=0.15, help="Minimum NDVI threshold for edge estimation and TVDI calculation.", ) parser.add_argument( "--max-ndvi", type=float, default=0.95, help="Maximum NDVI threshold for edge estimation.", ) parser.add_argument( "--bin-count", type=int, default=35, help="Number of NDVI bins used for dry/wet edge extraction.", ) parser.add_argument( "--dry-threshold", type=float, default=0.95, help="Percentile used to select the dry-edge temperature within each NDVI bin.", ) parser.add_argument( "--wet-threshold", type=float, default=0.05, help="Percentile used to select the wet-edge temperature within each NDVI bin.", ) return parser.parse_args()def reproject_to_match(source_path, target_profile): with rasterio.open(source_path) as src: dest = np.full((target_profile["height"], target_profile["width"]), np.nan, dtype=np.float32) reproject( source=rasterio.band(src, 1), destination=dest, src_transform=src.transform, src_crs=src.crs, dst_transform=target_profile["transform"], dst_crs=target_profile["crs"], resampling=Resampling.bilinear, src_nodata=src.nodata, dst_nodata=np.nan, ) return destdef calculate_ndvi(red, nir): red = red.astype(np.float32) nir = nir.astype(np.float32) denom = nir + red with np.errstate(divide="ignore", invalid="ignore"): ndvi = (nir - red) / denom ndvi[denom == 0] = np.nan return ndvidef fit_edge(ndvi_values, lst_values, bin_count, quantile): valid = ~np.isnan(ndvi_values) & ~np.isnan(lst_values) ndvi_values = ndvi_values[valid] lst_values = lst_values[valid] if ndvi_values.size == 0: raise ValueError("No valid NDVI/LST pixels for edge fitting.") bins = np.linspace(np.nanmin(ndvi_values), np.nanmax(ndvi_values), bin_count + 1) bin_centers = [] edge_lst = [] for left, right in zip(bins[:-1], bins[1:]): mask = (ndvi_values >= left) & (ndvi_values < right) if mask.sum() < 50: continue values = lst_values[mask] edge_value = np.nanpercentile(values, quantile) if np.isfinite(edge_value): bin_centers.append((left + right) / 2.0) edge_lst.append(edge_value) if len(bin_centers) < 2: raise ValueError("Not enough bins with valid values to fit a line.") slope, intercept = np.polyfit(bin_centers, edge_lst, 1) return slope, intercept, np.array(bin_centers), np.array(edge_lst)def calculate_tvdi(lst, ndvi, dry_line, wet_line): dry_slope, dry_intercept = dry_line wet_slope, wet_intercept = wet_line lst_dry = dry_slope * ndvi + dry_intercept lst_wet = wet_slope * ndvi + wet_intercept denom = lst_dry - lst_wet with np.errstate(divide="ignore", invalid="ignore"): tvdi = (lst - lst_wet) / denom tvdi[denom == 0] = np.nan tvdi = np.clip(tvdi, 0.0, 1.0) return tvdidef classify_tvdi(tvdi): classes = np.full(tvdi.shape, 0, dtype=np.uint8) valid = ~np.isnan(tvdi) classes[(tvdi >= 0.0) & (tvdi < 0.2)] = 1 classes[(tvdi >= 0.2) & (tvdi < 0.4)] = 2 classes[(tvdi >= 0.4) & (tvdi < 0.6)] = 3 classes[(tvdi >= 0.6) & (tvdi < 0.8)] = 4 classes[(tvdi >= 0.8) & (tvdi <= 1.0)] = 5 classes[~valid] = 0 return classesdef save_raster(path, array, profile, dtype, nodata=None, compress="lzw"): profile = profile.copy() profile.update( dtype=dtype, count=1, compress=compress, nodata=nodata, ) with rasterio.open(path, "w", **profile) as dst: dst.write(array.astype(dtype), 1)def main(): args = parse_args() os.makedirs(args.output_dir, exist_ok=True) with rasterio.open(args.multispectral) as ms: red = ms.read(args.red_band).astype(np.float32) nir = ms.read(args.nir_band).astype(np.float32) ms_profile = ms.profile ndvi = calculate_ndvi(red, nir) print(f"Computed NDVI from bands {args.red_band} and {args.nir_band}.") lst = reproject_to_match(args.thermal, ms_profile) print(f"Reprojected thermal raster to multispectral grid ({lst.shape}).") valid_mask = ( ~np.isnan(ndvi) & ~np.isnan(lst) & (ndvi >= args.min_ndvi) & (ndvi <= args.max_ndvi) & (lst > 0) & (lst < 500) ) if valid_mask.sum() < 1000: raise RuntimeError("Too few valid pixels for edge estimation after masking.") ndvi_for_edge = ndvi[valid_mask] lst_for_edge = lst[valid_mask] dry_line = fit_edge(ndvi_for_edge, lst_for_edge, args.bin_count, 100 * args.dry_threshold) wet_line = fit_edge(ndvi_for_edge, lst_for_edge, args.bin_count, 100 * args.wet_threshold) dry_slope, dry_intercept, dry_centers, dry_values = dry_line wet_slope, wet_intercept, wet_centers, wet_values = wet_line print("Dry edge: LST = {:.3f} * NDVI + {:.3f}".format(dry_slope, dry_intercept)) print("Wet edge: LST = {:.3f} * NDVI + {:.3f}".format(wet_slope, wet_intercept)) tvdi = calculate_tvdi(lst, ndvi, (dry_slope, dry_intercept), (wet_slope, wet_intercept)) classification = classify_tvdi(tvdi) tvdi_path = os.path.join(args.output_dir, "TVDI_20260708.tif") class_path = os.path.join(args.output_dir, "Drought_Class_20260708.tif") save_raster(tvdi_path, tvdi.astype(np.float32), ms_profile, dtype=rasterio.float32, nodata=np.nan) save_raster(class_path, classification, ms_profile, dtype=rasterio.uint8, nodata=0) print(f"Saved TVDI raster to {tvdi_path}") print(f"Saved drought class raster to {class_path}") print("Drought classes: 0=invalid, 1=no stress, 2=light, 3=moderate, 4=severe, 5=extreme")if __name__ == "__main__": main()
10个月宝宝每天需要喝多少奶粉？
