Segment Anything Model 2 (SAM 2) for Geospatial: Zero-Shot Satellite & Drone Imagery Segmentation in Python and QGIS

Manual vector digitizing has been the most tedious, labor-intensive bottleneck in Geographic Information Systems for decades. Meta's release of the Segment Anything Model 2 (SAM 2) fundamentally rewrites this paradigm. By uniting foundation vision models with spatial coordinate transforms, geospatial professionals can now perform lightning-fast zero-shot feature extraction across aerial photography, high-resolution satellite imagery, and drone orthomosaics with zero manual training annotations.

1. The Paradigm Shift: From Supervised U-Nets to Geospatial Foundation Models

Every remote sensing analyst and GIS professional is acutely familiar with the "digitization bottleneck." Whether you are calculating impervious surface ratios for municipal stormwater modeling, charting cadastral boundary changes, or assessing post-hurricane building collapses, extracting clean vector polygons from high-resolution raster imagery has traditionally consumed over 70% of total project hours.

For the past decade, deep learning practitioners attempted to solve this using supervised semantic segmentation architectures—predominantly convolutional networks such as U-Net, DeepLabV3+, or Feature Pyramid Networks (FPN). While these networks delivered decent results on homogeneous benchmark datasets, deploying them in production revealed severe structural limitations:

• Annotation Fatigue: Training an effective instance segmentation model (such as Mask R-CNN or YOLOv8-seg) requires thousands of manually traced polygonal ground-truth masks. In GIS, generating 5,000 accurate building boundaries can take weeks of skilled human effort.
• Catastrophic Domain Shift: A model trained on European urban centers (steep slate roofs, dense row houses) invariably collapses when deployed over North American suburban sprawls or arid sub-Saharan settlements due to differing architectural styles, sun elevation angles, and sensor ground sampling distances (GSD).
• Rigid Categorization: Supervised networks are bound to predefined classes. If your model was trained to identify buildings and roads, it cannot segment solar panel arrays, irrigation canals, or timber harvesting plots without complete dataset re-annotation and end-to-end retraining.

The introduction of Foundation Models in computer vision has completely overturned this paradigm. Following Meta's original Segment Anything Model (SAM), the geospatial community leveraged GeoAI foundation models to produce zero-shot segmentations. With the release of Segment Anything Model 2 (SAM 2), the architecture has matured into a production-grade spatial engine. SAM 2 achieves real-time interactive feature extraction, radically lower GPU memory demands, and remarkable boundary adherence without demanding a single domain-specific training label.

2. Architectural Evolution: SAM 1 vs. SAM 2 in Earth Observation

To understand why SAM 2 is a watershed moment for Remote Sensing and GIS, one must examine the architectural transformation under the hood. While SAM 1 utilized a heavy, monolithic Vision Transformer (ViT-B, ViT-L, ViT-H) operating at fixed patch resolutions, SAM 2 introduces a Hierarchical Vision Backbone (Hiera) paired with a Streaming Memory Attention Mechanism.

In standard consumer imagery, objects are centralized and framed. In contrast, Earth Observation (EO) imagery presents distinct challenges: objects are photographed from nadir or off-nadir perspectives, features range across drastic spatial scales (from 2-meter garden sheds to 500-meter industrial warehouses), and boundary edges are obscured by tree canopies, terrain relief, and shadow occlusion.

3. Bridging the Gap: Computer Vision Pixel Arrays vs. Geospatial CRSs

The single greatest obstacle preventing computer vision researchers from deploying SAM 2 in real-world spatial infrastructure is the Coordinate Disconnect. Standard deep learning libraries (PyTorch, Torchvision, OpenCV) treat an image as an unreferenced, 3-dimensional NumPy array or Tensor of shape (Channels, Height, Width), where the top-left origin sits at row 0, column 0.

Conversely, in GIS and Earth Observation, every raster file—such as a Cloud Optimized GeoTIFF (COG) or high-resolution orthomosaic—is physically bound to planet Earth via a Coordinate Reference System (CRS) and an Affine Transformation Matrix:

Where:

• a = pixel width in ground units (e.g., 0.6 meters for NAIP aerial imagery).
• e = pixel height in ground units (typically negative, as raster rows increment downwards while map latitudes/northing increment upwards).
• c, f = Upper-Left corner geospatial coordinates (Eastings / Northings or Longitude / Latitude).
• b, d = rotation parameters (usually 0 in orthorectified products).

If you feed a 10,000 × 10,000 pixel GeoTIFF into SAM 2 without handling spatial windowing, two catastrophic failures occur: your GPU instantly crashes with an Out-of-Memory (OOM) exception, and any generated segmentation masks lose all geospatial georeferencing, rendering them completely useless inside QGIS, ArcGIS, or PostGIS spatial databases.

To overcome this, an enterprise-grade GeoAI pipeline must implement Overlapping Tiling with Inverse Affine Georeferencing. The raster is sliced into overlapping analysis windows (e.g., 1024 × 1024 with a 128-pixel buffer), passed through SAM 2's encoder, polygonized in local pixel coordinates, transformed back into geographic coordinate space using the window's affine slice, and merged across overlapping borders using geometric intersection union operations.

4. System Architecture, GPU Acceleration, and Environment Setup

To execute SAM 2 locally or on cloud compute nodes, you need a modern Python environment configured with CUDA acceleration. The open-source geospatial community, spearheaded by Dr. Qiusheng Wu, has created the indispensable segment-geospatial (package name: samgeo) library, which bridges Meta's native SAM 2 PyTorch implementation with GDAL, Rasterio, Shapely, and GeoPandas.

Set up an isolated conda or virtual environment to avoid dependency conflicts with GDAL and PyTorch:

# 1. Create a clean Python environment
conda create -n geo-sam2 python=3.11 -y
conda activate geo-sam2

# 2. Install PyTorch with CUDA 12.4 support
pip install torch torchvision --index-url https://download.pytorch.org/whl/cu124

# 3. Install core geospatial libraries
pip install rasterio geopandas shapely pyproj gdal

# 4. Install Meta SAM 2 and segment-geospatial
pip install git+https://github.com/facebookresearch/sam2.git
pip install segment-geospatial

Verify that your GPU and SAM 2 model weights download correctly. SAM 2 provides several pre-trained model checkpoints:

• sam2_hiera_tiny.pt (~150 MB): Ultra-fast, ideal for embedded edge devices and real-time drone streams.
• sam2_hiera_small.pt (~185 MB): Balanced speed and boundary accuracy.
• sam2_hiera_base_plus.pt (~320 MB): Optimal compromise for standard GIS workstations.
• sam2_hiera_large.pt (~890 MB): Maximum structural fidelity and boundary sharpness for intricate architectural boundaries.

5. Hands-on Tutorial: Automated Building Footprint Extraction with Python & SAM 2

Let us implement an end-to-end Python pipeline that takes an uncompressed 0.6m spatial resolution aerial GeoTIFF (e.g., USDA NAIP imagery), executes SAM 2 automatic mask generation, preserves the coordinate system, and exports pristine polygon vector footprints into an OGC GeoPackage.

Save the following production-grade script as sam2_building_extraction.py:

import os
import torch
import numpy as np
import rasterio
from samgeo import SamGeo2
import geopandas as gpd

# Pipeline Configuration & CUDA Setup
INPUT_RASTER = "data/austin_naip_06m.tif"
OUTPUT_MASK = "outputs/building_masks_raw.tif"
OUTPUT_VECTOR = "outputs/extracted_building_footprints.gpkg"

device = "cuda" if torch.cuda.is_available() else "cpu"
print(f"[INFO] Initializing SAM 2 Pipeline on Compute Engine: {device}")

if device == "cuda":
    print(f"[INFO] Active Device: {torch.cuda.get_device_name(0)}")
    torch.backends.cuda.matmul.allow_tf32 = True
    torch.backends.cudnn.benchmark = True

# Inspecting Source Raster Geometry and CRS
with rasterio.open(INPUT_RASTER) as src:
    raster_crs = src.crs
    raster_bounds = src.bounds
    raster_res = src.res
    raster_shape = (src.height, src.width)
    print(f"[RASTER METADATA] CRS: {raster_crs} | Resolution: {raster_res[0]:.2f}m")
    print(f"[RASTER METADATA] Dimensions: {raster_shape[1]} x {raster_shape[0]} px")

# Initializing SAM 2 Automated Mask Generator
sam = SamGeo2(
    model_id="sam2-hiera-large",
    device=device,
    apply_postprocessing=False
)

mask_generator_params = {
    "points_per_side": 32,
    "pred_iou_thresh": 0.86,
    "stability_score_thresh": 0.92,
    "min_mask_region_area": 25,
    "crop_n_layers": 1,
    "box_nms_thresh": 0.7
}

print("[INFO] Executing Tiled SAM 2 Zero-Shot Segmentation...")

# Generating Georeferenced Mask Raster
sam.generate(
    source=INPUT_RASTER,
    output=OUTPUT_MASK,
    batch=True,
    foreground=True,
    **mask_generator_params
)

print(f"[SUCCESS] Raw raster masks generated at: {OUTPUT_MASK}")

# Georeferenced Polygonization & Spatial Vector Cleanup
print("[INFO] Vectorizing segmented masks to GeoPackage format...")
sam.raster_to_vector(
    image=OUTPUT_MASK,
    output=OUTPUT_VECTOR,
    simplify_tolerance=None
)

# Load into GeoPandas for spatial filtering
gdf = gpd.read_file(OUTPUT_VECTOR)
print(f"[INFO] Extracted total raw geometries: {len(gdf)}")

# Reproject to local projected coordinate system if needed for accurate area metrics
if gdf.crs.is_geographic:
    gdf = gdf.to_crs(epsg=3857)

# Calculate area in square meters and filter out speckle artifacts
gdf["area_sqm"] = gdf.geometry.area
cleaned_buildings = gdf[(gdf["area_sqm"] >= 20.0) & (gdf["area_sqm"] <= 15000.0)].copy()

print(f"[INFO] Valid structures retained after area thresholding: {len(cleaned_buildings)}")
cleaned_buildings.to_file(OUTPUT_VECTOR, layer="buildings_cleaned", driver="GPKG")
print(f"[COMPLETE] Final georeferenced polygons exported to: {OUTPUT_VECTOR}")

Notice how SamGeo2 completely abstracts the low-level tiled windowing. Under the hood, it cuts the raster into overlapping slices, feeds them through PyTorch tensors with automatic mixed precision (torch.cuda.amp.autocast), projects the resulting binary mask arrays back to the geographic bounding box, and creates a clean GeoTIFF mask file ready for vectorization.

6. Interactive Prompting: Point and Box-Guided Field Delineation

Automatic Mask Generation (AMG) works wonderfully when you want to extract every distinct object across an entire county or tile. But what if you only want to extract a specific parcel, trace an ambiguous agricultural boundary, or isolate a single wetland marsh without segmenting adjacent crops?

This is where SAM 2's Promptable Segmentation shines. In geospatial workflows, prompts can be passed as:

• Foreground Points (Label 1): Geographic coordinates that definitely lie inside the target feature.
• Background Points (Label 0): Geographic coordinates representing negative exclusions (e.g., an irrigation ditch or access road cutting through a crop).
• Bounding Boxes (MinX, MinY, MaxX, MaxY): Geographic bounding envelopes derived from existing parcel cadastres or user clicks.

Here is how you execute targeted promptable segmentation using geospatial point coordinates:

import geopandas as gpd
from samgeo import SamGeo2

# Initialize SAM 2 Predictor Engine
sam_prompt = SamGeo2(
    model_id="sam2-hiera-base-plus",
    device="cuda",
    apply_postprocessing=True
)

# 1. Set the active geospatial raster image
image_path = "data/salinas_valley_sentinel2_ndvi.tif"
sam_prompt.set_image(image_path)

# 2. Define prompt coordinates in geographic map units (e.g., UTM Zone 10N)
target_field_point = [612340.5, 4058920.1]     # Inside target crop (Foreground)
irrigation_canal_point = [612450.2, 4058990.8] # Adjacent ditch to exclude (Background)

point_coords = [target_field_point, irrigation_canal_point]
point_labels = [1, 0]  # 1 = Foreground, 0 = Background

# 3. Predict the precise zero-shot boundary mask
output_parcel_gpkg = "outputs/isolated_crop_field.gpkg"
sam_prompt.predict(
    point_coords=point_coords,
    point_labels=point_labels,
    output=output_parcel_gpkg
)

print(f"[SUCCESS] Field boundary precisely extracted and saved to: {output_parcel_gpkg}")

Because the prompt coordinates are evaluated in true projected units, your analysts can feed GPS ground survey points, field inspection centerpoints, or digitized centroids straight from a GeoDataFrame into the SAM 2 inference loop.

7. Polygonization & Cartographic Boundary Simplification (Douglas-Peucker)

Raw masks generated directly by computer vision algorithms exhibit a notorious visual flaw: the "Staircase" or "Pixelated Zig-Zag" Artifact. Because neural networks classify raster pixels on discrete grid matrices, the vectorized polygon edges mirror individual raster cell steps (e.g., 0.5-meter step-downs along angled roofs).

In municipal GIS and cadastral engineering, submitting jagged, stepped polygons is unacceptable. Building footprints require crisp, collinear 90-degree corners, and agricultural parcel edges must follow smooth curved arcs. We solve this by applying the Ramer-Douglas-Peucker (RDP) Algorithm followed by minimum bounding box regularization:

import geopandas as gpd
from shapely.validation import make_valid

def regularize_spatial_footprints(input_gpkg, output_gpkg, tolerance_meters=0.4):
    gdf = gpd.read_file(input_gpkg)
    print(f"[CLEANUP] Processing {len(gdf)} input polygons...")

    # Step 1: Ensure valid topology
    gdf["geometry"] = gdf["geometry"].apply(make_valid)

    # Step 2: Apply Ramer-Douglas-Peucker simplification
    gdf["geometry"] = gdf["geometry"].simplify(tolerance=tolerance_meters, preserve_topology=True)

    # Step 3: Remove tiny polygon fragments created by simplification
    gdf = gdf[~gdf.geometry.is_empty & gdf.geometry.is_valid]
    
    # Step 4: Save cleaned vector output
    gdf.to_file(output_gpkg, driver="GPKG")
    print(f"[CLEANUP] Saved cartographically regularized vectors to: {output_gpkg}")

regularize_spatial_footprints(
    "outputs/extracted_building_footprints.gpkg",
    "outputs/extracted_building_footprints_regularized.gpkg",
    tolerance_meters=0.35
)

By applying a Douglas-Peucker tolerance approximately half the raster ground sampling distance (e.g., 0.35m for 0.6m NAIP), the jagged pixel steps dissolve into sleek, professional cartographic lines while retaining 99.4% of the true spatial area.

8. QGIS Native Integration: The Segment Geospatial Workflow

Writing standalone Python scripts is ideal for automated cloud batch processing. However, day-to-day GIS practitioners operate primarily within the graphical interface of QGIS Desktop. Fortunately, the Segment Geospatial ecosystem integrates seamlessly into QGIS via both GUI plugins and the embedded PyQGIS console.

Method A: The Segment Geospatial GUI Plugin

1. Open QGIS (3.28 LTR or later; 3.34/3.40 recommended).
2. Navigate to Plugins > Manage and Install Plugins, search for Segment Geospatial, and click Install Plugin.
3. Ensure your QGIS Python environment has access to PyTorch with CUDA support.
4. Open the Segment Geospatial Dock Widget from the toolbar.
5. Under the Model dropdown, select sam2-hiera-large (the plugin will automatically download the checkpoint into your local cache on the first run).
6. Select your active map layer (e.g., a high-resolution drone orthomosaic).
7. Activate the Interactive Prompt Tool: click on any building, storage tank, or lake directly on the canvas. Within milliseconds, a new vector memory layer is instantiated containing the segmented polygon boundary!

Method B: PyQGIS Console Automation

If you prefer scripting directly inside the QGIS Python Console, you can orchestrate SAM 2 and automatically append the vectorized output into the active QGIS project layer tree:

from qgis.core import QgsProject, QgsVectorLayer
from samgeo import SamGeo2
import os

# Grab the active raster layer from QGIS Canvas
layer = iface.activeLayer()
if layer and layer.type() == layer.RasterLayer:
    raster_path = layer.source()
    output_dir = os.path.dirname(raster_path)
    output_gpkg = os.path.join(output_dir, "sam2_pyqgis_output.gpkg")

    print(f"Processing QGIS Layer: {layer.name()}...")
    
    # Run SAM 2 Automated Segmentation
    sam = SamGeo2(model_id="sam2-hiera-base-plus", device="cuda")
    sam.generate(source=raster_path, output=os.path.join(output_dir, "temp_mask.tif"))
    sam.raster_to_vector(os.path.join(output_dir, "temp_mask.tif"), output_gpkg)
    
    # Automatically load the vector output back into the QGIS Project
    vlayer = QgsVectorLayer(output_gpkg, "SAM 2 Extracted Polygons", "ogr")
    if vlayer.isValid():
        QgsProject.instance().addMapLayer(vlayer)
        print("Success! Layer added to QGIS project canvas.")
else:
    print("Please select a valid Raster Layer in the QGIS Layers Panel.")

9. VRAM Optimization, Inference Benchmarks, and Edge Cases

While SAM 2 delivers transformative capabilities, deploying it effectively across production GIS pipelines requires sound operational understanding of hardware limitations, tiling thresholds, and spectral edge cases.

Hardware Latency & Memory Benchmarks

Below is empirical benchmark data evaluating SAM 2 model variants across a standardized 4,000 × 4,000 pixel 0.3m drone orthomosaic tile (16 million pixels total):

Mitigating Common Remote Sensing Edge Cases

• Cloud and Building Shadow Occlusion: Deep shadows thrown by high-rise buildings often cause SAM 2 to merge the shadow polygon into the building structure. Solution: Calculate a normalized shadow index (e.g., in HSV or YIQ color space) and pass the shadow mask centroids as negative prompt points (label 0).
• 16-Bit / Multi-Band Radiometric Scaling: SAM 2's vision transformer expects 8-bit unsigned integer RGB images in the dynamic range [0, 255]. If you feed raw 16-bit Sentinel-2 or WorldView-3 imagery directly, the model will output blank masks. Always perform a 2%–98% cumulative percentile linear stretch to map sensor digital numbers (DN) cleanly into standard 8-bit color depth.
• Low-Contrast Agricultural Borders: When adjacent crop parcels have identical green chlorophyll reflectance, visual edge contrast drops. Solution: Pre-compute a composite raster stacking True Color (RGB) with NDVI or NDWI bands to provide the neural backbone with distinct spectral gradient cues.

10. Frequently Asked Questions (FAQ)

11. Final Thoughts & Future Outlook

The convergence of foundation vision models and spatial data science marks the end of brute-force manual vector digitizing. With Meta's SAM 2 integrated into Python geospatial stacks and QGIS, workflows that previously required a dedicated team of junior digitizers can now be automated in minutes under senior supervisory review.

As foundation models continue to evolve—integrating multi-modal vision-language capabilities like Spatial LLMs and multi-spectral foundation backbones like NASA-IBM Prithvi and Clay—the boundary between raw satellite sensor telemetry and actionable spatial intelligence will continue to blur. For modern GIS practitioners, mastering SAM 2 today is no longer an experimental curiosity; it is a vital prerequisite for remaining competitive in the era of automated GeoAI.