农业地块agricultural parcel / field parcel是土地管理与农情监测的基本空间单元其边界的准确勾绘是农业补贴核发、农业保险定损、耕地面积统计、精准农业变量作业与作物产量估测等应用的基础。随着无人机UAV与高分辨率遥感影像的普及基于影像自动提取田块边界、替代人工数字化作业已成为遥感与农业信息化领域的重要研究问题#!/usr/bin/env python3 Batch-segment farmland images and estimate each parcels area. from __future__ import annotations import argparse import csv import json import math import sys from dataclasses import dataclass from pathlib import Path from typing import Iterable import cv2 import numpy as np IMAGE_EXTENSIONS {.jpg, .jpeg, .png, .bmp, .tif, .tiff, .webp} dataclass(frozenTrue) class SegmentationConfig: max_dimension: int min_area_pixels: int min_saturation: int min_excess_green: float | None max_texture: float boundary_percentile: float split_width: int dataclass(frozenTrue) class AreaCalibration: pixel_area: float | None unit: str source: str def convert(self, pixels: int) - float: if self.pixel_area is None: return float(pixels) return pixels * self.pixel_area property def output_unit(self) - str: return self.unit if self.pixel_area is not None else pixel^2 def parse_args() - argparse.Namespace: parser argparse.ArgumentParser( description( Segment cultivated parcels in one image or every image in a folder, write binary/label/overlay images, and export parcel areas. ) ) parser.add_argument( input, nargs?, default., helpInput image or folder (default: current folder)., ) parser.add_argument( -o, --output-dir, defaultoutput, helpOutput folder (default: output)., ) calibration parser.add_argument_group(area calibration) calibration.add_argument( --pixel-area, typefloat, helpKnown physical area represented by one original-image pixel., ) calibration.add_argument( --reference-area, typefloat, helpKnown physical area of a reference region., ) calibration.add_argument( --reference-pixels, typefloat, helpPixel count of the same reference region., ) calibration.add_argument( --unit, defaultm2, helpPhysical area unit used by the calibration (default: m2)., ) segmentation parser.add_argument_group(segmentation) segmentation.add_argument( --min-area-pixels, typeint, default4000, helpDiscard parcels smaller than this many original pixels (default: 4000)., ) segmentation.add_argument( --max-dimension, typeint, default1800, helpMaximum processing width/height; area is measured at original size., ) segmentation.add_argument( --min-saturation, typeint, default22, helpMinimum HSV saturation for cultivated vegetation (default: 22)., ) segmentation.add_argument( --min-excess-green, typefloat, helpManual ExG threshold. By default an Otsu threshold is used., ) segmentation.add_argument( --max-texture, typefloat, default30.0, helpLocal grayscale standard-deviation limit used to suppress woodland., ) segmentation.add_argument( --boundary-percentile, typefloat, default65.0, helpColor-gradient percentile treated as parcel boundaries (default: 65)., ) segmentation.add_argument( --split-width, typeint, default2, helpBoundary dilation width at processing resolution (default: 2)., ) parser.add_argument( --save-block-masks, actionstore_true, helpAlso save one full-resolution binary PNG for every parcel., ) return parser.parse_args() def validate_args(args: argparse.Namespace) - tuple[SegmentationConfig, AreaCalibration]: if args.max_dimension 256: raise ValueError(--max-dimension must be at least 256) if args.min_area_pixels 1: raise ValueError(--min-area-pixels must be positive) if not 0 args.min_saturation 255: raise ValueError(--min-saturation must be between 0 and 255) if args.max_texture 0: raise ValueError(--max-texture must be positive) if not 50 args.boundary_percentile 100: raise ValueError(--boundary-percentile must be in [50, 100)) if args.split_width 0: raise ValueError(--split-width cannot be negative) has_reference_area args.reference_area is not None has_reference_pixels args.reference_pixels is not None if has_reference_area ! has_reference_pixels: raise ValueError( --reference-area and --reference-pixels must be provided together ) if args.pixel_area is not None and has_reference_area: raise ValueError( use either --pixel-area or --reference-area/--reference-pixels, not both ) for name in (pixel_area, reference_area, reference_pixels): value getattr(args, name) if value is not None and value 0: raise ValueError(f--{name.replace(_, -)} must be positive) if args.pixel_area is not None: calibration AreaCalibration(args.pixel_area, args.unit, pixel_area) elif has_reference_area: pixel_area args.reference_area / args.reference_pixels calibration AreaCalibration(pixel_area, args.unit, reference_region) else: calibration AreaCalibration(None, pixel^2, uncalibrated) config SegmentationConfig( max_dimensionargs.max_dimension, min_area_pixelsargs.min_area_pixels, min_saturationargs.min_saturation, min_excess_greenargs.min_excess_green, max_textureargs.max_texture, boundary_percentileargs.boundary_percentile, split_widthargs.split_width, ) return config, calibration def find_images(input_path: Path, output_dir: Path) - list[Path]: if input_path.is_file(): if input_path.suffix.lower() not in IMAGE_EXTENSIONS: raise ValueError(funsupported image extension: {input_path.suffix}) return [input_path] if not input_path.is_dir(): raise FileNotFoundError(finput does not exist: {input_path}) output_resolved output_dir.resolve() images [] for path in input_path.rglob(*): if not path.is_file() or path.suffix.lower() not in IMAGE_EXTENSIONS: continue try: path.resolve().relative_to(output_resolved) except ValueError: images.append(path) return sorted(images) def read_image(path: Path) - np.ndarray: data np.fromfile(path, dtypenp.uint8) image cv2.imdecode(data, cv2.IMREAD_COLOR) if image is None: raise ValueError(fcannot decode image: {path}) return image def write_image(path: Path, image: np.ndarray, quality: int 95) - None: path.parent.mkdir(parentsTrue, exist_okTrue) suffix path.suffix.lower() params: list[int] [] if suffix in {.jpg, .jpeg}: params [cv2.IMWRITE_JPEG_QUALITY, quality] success, encoded cv2.imencode(suffix, image, params) if not success: raise ValueError(fcannot encode output image: {path}) encoded.tofile(path) def resize_for_processing(image: np.ndarray, max_dimension: int) - tuple[np.ndarray, float]: height, width image.shape[:2] scale min(1.0, max_dimension / max(height, width)) if scale 1.0: return image.copy(), scale resized cv2.resize( image, (max(1, round(width * scale)), max(1, round(height * scale))), interpolationcv2.INTER_AREA, ) return resized, scale def odd_kernel_size(value: int) - int: return max(3, value if value % 2 1 else value 1) def build_farmland_mask( image: np.ndarray, config: SegmentationConfig ) - tuple[np.ndarray, dict[str, float]]: blue, green, red cv2.split(image.astype(np.float32)) excess_green 2.0 * green - red - blue exg_u8 np.clip(excess_green 128.0, 0, 255).astype(np.uint8) if config.min_excess_green is None: otsu_value, _ cv2.threshold( exg_u8, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU ) exg_threshold max(4.0, float(otsu_value) - 128.0) else: exg_threshold config.min_excess_green hsv cv2.cvtColor(image, cv2.COLOR_BGR2HSV) hue, saturation, value cv2.split(hsv) vegetation ( (excess_green exg_threshold) (hue 18) (hue 100) (saturation config.min_saturation) (value 24) ) gray cv2.cvtColor(image, cv2.COLOR_BGR2GRAY).astype(np.float32) texture_window odd_kernel_size(round(min(image.shape[:2]) / 85)) local_mean cv2.boxFilter(gray, cv2.CV_32F, (texture_window, texture_window)) local_square_mean cv2.boxFilter( gray * gray, cv2.CV_32F, (texture_window, texture_window) ) local_std np.sqrt(np.maximum(local_square_mean - local_mean * local_mean, 0)) mask (vegetation (local_std config.max_texture)).astype(np.uint8) * 255 base max(1, round(min(image.shape[:2]) / 600)) open_kernel cv2.getStructuringElement( cv2.MORPH_ELLIPSE, (odd_kernel_size(base), odd_kernel_size(base)) ) close_size odd_kernel_size(base * 3) close_kernel cv2.getStructuringElement( cv2.MORPH_ELLIPSE, (close_size, close_size) ) mask cv2.morphologyEx(mask, cv2.MORPH_OPEN, open_kernel) mask cv2.morphologyEx(mask, cv2.MORPH_CLOSE, close_kernel) diagnostics { excess_green_threshold: round(float(exg_threshold), 3), vegetation_fraction: round(float(np.count_nonzero(mask) / mask.size), 6), texture_window: texture_window, } return mask, diagnostics def color_gradient(image: np.ndarray) - np.ndarray: smoothed cv2.GaussianBlur(image, (0, 0), sigmaX2.0, sigmaY2.0) lab cv2.cvtColor(smoothed, cv2.COLOR_BGR2LAB) gradients [] for channel in cv2.split(lab): gx cv2.Sobel(channel, cv2.CV_32F, 1, 0, ksize3) gy cv2.Sobel(channel, cv2.CV_32F, 0, 1, ksize3) gradients.append(cv2.magnitude(gx, gy)) return np.maximum.reduce(gradients) def segment_parcels( image: np.ndarray, farmland_mask: np.ndarray, config: SegmentationConfig, processing_scale: float, ) - tuple[np.ndarray, dict[str, float]]: gradient color_gradient(image) mask_values gradient[farmland_mask 0] if mask_values.size 0: return np.zeros(farmland_mask.shape, dtypenp.int32), { boundary_threshold: 0.0, seed_count: 0, } boundary_threshold float( np.percentile(mask_values, config.boundary_percentile) ) boundaries ((gradient boundary_threshold) (farmland_mask 0)).astype( np.uint8 ) boundary_kernel cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3, 3)) boundaries cv2.morphologyEx(boundaries, cv2.MORPH_CLOSE, boundary_kernel) if config.split_width: boundaries cv2.dilate( boundaries, boundary_kernel, iterationsconfig.split_width ) core cv2.bitwise_and(farmland_mask, cv2.bitwise_not(boundaries * 255)) core cv2.morphologyEx(core, cv2.MORPH_OPEN, boundary_kernel) component_count, component_labels, stats, _ cv2.connectedComponentsWithStats( core, connectivity8 ) min_area_working max( 12, int(math.ceil(config.min_area_pixels * processing_scale * processing_scale)) ) min_seed_area max(8, min_area_working // 8) markers np.zeros(farmland_mask.shape, dtypenp.int32) markers[farmland_mask 0] 1 next_marker 2 for component_id in range(1, component_count): if stats[component_id, cv2.CC_STAT_AREA] min_seed_area: continue markers[component_labels component_id] next_marker next_marker 1 seed_count next_marker - 2 if seed_count 0: _, fallback cv2.connectedComponents(farmland_mask, connectivity8) return fallback.astype(np.int32), { boundary_threshold: round(boundary_threshold, 3), seed_count: int(fallback.max()), } watershed_input cv2.GaussianBlur(image, (0, 0), sigmaX1.2, sigmaY1.2) watershed cv2.watershed(watershed_input, markers) raw_labels np.where((watershed 2) (farmland_mask 0), watershed - 1, 0) raw_labels assign_unlabeled_pixels(raw_labels.astype(np.int32), farmland_mask) labels remove_small_regions(raw_labels.astype(np.int32), min_area_working) diagnostics { boundary_threshold: round(boundary_threshold, 3), seed_count: seed_count, } return labels, diagnostics def assign_unlabeled_pixels( labels: np.ndarray, allowed_mask: np.ndarray ) - np.ndarray: if not np.any(labels 0): return labels distance_source np.where(labels 0, 0, 255).astype(np.uint8) _, nearest cv2.distanceTransformWithLabels( distance_source, cv2.DIST_L2, 5, labelTypecv2.DIST_LABEL_PIXEL, ) nearest_to_region np.zeros(int(nearest.max()) 1, dtypenp.int32) source_pixels distance_source 0 nearest_to_region[nearest[source_pixels]] labels[source_pixels] filled nearest_to_region[nearest] filled[allowed_mask 0] 0 return split_disconnected_regions(filled) def split_disconnected_regions(labels: np.ndarray) - np.ndarray: result np.zeros_like(labels, dtypenp.int32) next_id 1 for old_id in range(1, int(labels.max()) 1): ys, xs np.where(labels old_id) if xs.size 0: continue x0, x1 int(xs.min()), int(xs.max()) 1 y0, y1 int(ys.min()), int(ys.max()) 1 region (labels[y0:y1, x0:x1] old_id).astype(np.uint8) component_count, components cv2.connectedComponents(region, connectivity8) target result[y0:y1, x0:x1] for component_id in range(1, component_count): target[components component_id] next_id next_id 1 return result def remove_small_regions(labels: np.ndarray, min_area: int) - np.ndarray: counts np.bincount(labels.ravel()) kept_ids np.flatnonzero(counts min_area) kept_ids kept_ids[kept_ids ! 0] lookup np.zeros(len(counts), dtypenp.int32) lookup[kept_ids] np.arange(1, len(kept_ids) 1, dtypenp.int32) return lookup[labels] def upscale_and_measure( working_labels: np.ndarray, original_shape: tuple[int, int], min_area_pixels: int, ) - tuple[np.ndarray, list[dict[str, int]]]: original_height, original_width original_shape labels cv2.resize( working_labels, (original_width, original_height), interpolationcv2.INTER_NEAREST, ).astype(np.int32) labels remove_small_regions(labels, min_area_pixels) rows: list[dict[str, int]] [] for block_id in range(1, int(labels.max()) 1): ys, xs np.where(labels block_id) if xs.size 0: continue rows.append( { block_id: block_id, pixel_count: int(xs.size), centroid_x: int(round(float(xs.mean()))), centroid_y: int(round(float(ys.mean()))), bbox_x: int(xs.min()), bbox_y: int(ys.min()), bbox_width: int(xs.max() - xs.min() 1), bbox_height: int(ys.max() - ys.min() 1), } ) return labels, rows def palette_color(block_id: int) - tuple[int, int, int]: hue (block_id * 47) % 180 hsv np.uint8([[[hue, 190, 235]]]) bgr cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR)[0, 0] return int(bgr[0]), int(bgr[1]), int(bgr[2]) def render_labels(labels: np.ndarray) - np.ndarray: colored np.zeros((*labels.shape, 3), dtypenp.uint8) for block_id in range(1, int(labels.max()) 1): colored[labels block_id] palette_color(block_id) return colored def format_area_label(area: float) - str: absolute abs(area) if absolute 1_000_000: return f{area / 1_000_000:.2f}M if absolute 10_000: return f{area / 1_000:.0f}K if absolute 1_000: return f{area / 1_000:.1f}K return f{area:.2f} def find_label_position(labels: np.ndarray, row: dict[str, int | float | str]) - tuple[int, int]: block_id int(row[block_id]) x int(row[bbox_x]) y int(row[bbox_y]) width int(row[bbox_width]) height int(row[bbox_height]) region (labels[y : y height, x : x width] block_id).astype(np.uint8) distance cv2.distanceTransform(region, cv2.DIST_L2, 5) _, _, _, max_location cv2.minMaxLoc(distance) return x max_location[0], y max_location[1] def render_overlay( image: np.ndarray, labels: np.ndarray, rows: list[dict[str, int | float | str]], ) - np.ndarray: colored render_labels(labels) foreground labels 0 overlay image.copy() blended cv2.addWeighted(image, 0.62, colored, 0.38, 0) overlay[foreground] blended[foreground] contour_mask (foreground.astype(np.uint8) * 255) contours, _ cv2.findContours( contour_mask, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE ) contour_width max(1, round(min(image.shape[:2]) / 1400)) cv2.drawContours( overlay, contours, -1, (255, 255, 255), contour_width, cv2.LINE_AA ) font_scale max(0.4, min(0.72, min(image.shape[:2]) / 3000)) thickness 1 for row in rows: text f{row[block_id]}:{format_area_label(float(row[area]))} (text_width, text_height), _ cv2.getTextSize( text, cv2.FONT_HERSHEY_SIMPLEX, font_scale, thickness ) if ( text_width 8 int(row[bbox_width]) or text_height 8 int(row[bbox_height]) ): text str(row[block_id]) (text_width, text_height), _ cv2.getTextSize( text, cv2.FONT_HERSHEY_SIMPLEX, font_scale, thickness ) origin find_label_position(labels, row) text_origin ( max(0, min(origin[0] - text_width // 2, image.shape[1] - text_width)), max(text_height, min(origin[1] text_height // 2, image.shape[0] - 3)), ) cv2.putText( overlay, text, text_origin, cv2.FONT_HERSHEY_SIMPLEX, font_scale, (0, 0, 0), thickness 2, cv2.LINE_AA, ) cv2.putText( overlay, text, text_origin, cv2.FONT_HERSHEY_SIMPLEX, font_scale, (255, 255, 255), thickness, cv2.LINE_AA, ) return overlay def write_csv( path: Path, source_name: str, rows: Iterable[dict[str, int | float | str]] ) - None: fieldnames [ source_image, block_id, pixel_count, area, unit, centroid_x, centroid_y, bbox_x, bbox_y, bbox_width, bbox_height, ] with path.open(w, encodingutf-8-sig, newline) as csv_file: writer csv.DictWriter(csv_file, fieldnamesfieldnames) writer.writeheader() for row in rows: writer.writerow({source_image: source_name, **row}) def process_image( image_path: Path, image_output_dir: Path, config: SegmentationConfig, calibration: AreaCalibration, save_block_masks: bool, ) - dict[str, object]: image read_image(image_path) working_image, scale resize_for_processing(image, config.max_dimension) farmland_mask, mask_diagnostics build_farmland_mask(working_image, config) working_labels, split_diagnostics segment_parcels( working_image, farmland_mask, config, scale ) labels, measured_rows upscale_and_measure( working_labels, image.shape[:2], config.min_area_pixels ) rows: list[dict[str, int | float | str]] [] for measured in measured_rows: rows.append( { **measured, area: round(calibration.convert(measured[pixel_count]), 6), unit: calibration.output_unit, } ) image_output_dir.mkdir(parentsTrue, exist_okTrue) binary (labels 0).astype(np.uint8) * 255 write_image(image_output_dir / binary.png, binary) write_image(image_output_dir / labels.png, render_labels(labels)) write_image( image_output_dir / overlay.jpg, render_overlay(image, labels, rows), ) write_csv(image_output_dir / areas.csv, image_path.name, rows) if save_block_masks: block_dir image_output_dir / blocks block_dir.mkdir(parentsTrue, exist_okTrue) for row in rows: block_id int(row[block_id]) block_mask (labels block_id).astype(np.uint8) * 255 write_image(block_dir / fblock_{block_id:04d}.png, block_mask) total_pixels int(np.count_nonzero(labels)) summary: dict[str, object] { source_image: str(image_path), image_width: int(image.shape[1]), image_height: int(image.shape[0]), processing_scale: round(scale, 6), block_count: len(rows), foreground_pixels: total_pixels, total_area: round(calibration.convert(total_pixels), 6), area_unit: calibration.output_unit, pixel_area: calibration.pixel_area, calibration_source: calibration.source, parameters: { min_area_pixels: config.min_area_pixels, min_saturation: config.min_saturation, min_excess_green: config.min_excess_green, max_texture: config.max_texture, boundary_percentile: config.boundary_percentile, split_width: config.split_width, }, diagnostics: {**mask_diagnostics, **split_diagnostics}, } with (image_output_dir / summary.json).open(w, encodingutf-8) as file: json.dump(summary, file, ensure_asciiFalse, indent2) return summary def main() - int: args parse_args() try: config, calibration validate_args(args) input_path Path(args.input).resolve() output_dir Path(args.output_dir).resolve() images find_images(input_path, output_dir) if not images: raise FileNotFoundError(fno supported images found under: {input_path}) print(fFound {len(images)} image(s).) failures 0 batch_rows [] for index, image_path in enumerate(images, start1): print(f[{index}/{len(images)}] Processing {image_path.name} ...) try: image_output_dir output_dir / image_path.stem summary process_image( image_path, image_output_dir, config, calibration, args.save_block_masks, ) batch_rows.append(summary) print( f {summary[block_count]} blocks, ftotal {summary[total_area]} {summary[area_unit]} ) except Exception as error: failures 1 print(f ERROR: {error}, filesys.stderr) output_dir.mkdir(parentsTrue, exist_okTrue) with (output_dir / batch_summary.json).open(w, encodingutf-8) as file: json.dump(batch_rows, file, ensure_asciiFalse, indent2) if failures: print(fCompleted with {failures} failed image(s)., filesys.stderr) return 1 print(fResults written to: {output_dir}) return 0 except (FileNotFoundError, ValueError) as error: print(fERROR: {error}, filesys.stderr) return 2 if __name__ __main__: raise SystemExit(main())