Segment Anything on satellite and aerial imagery in Python: a samgeo tutorial

Short answer: to segment satellite or aerial imagery with Segment Anything in Python, install the CPU build of PyTorch and segment-geospatial. Save an RGB GeoTIFF of about 1024 × 1024 pixels, and create SamGeo(model_type="vit_b", automatic=True). Then call sam.generate("image.tif", output="masks.tif") and sam.tiff_to_gpkg(...) to get one polygon per object. SAM gives you shapes, not classes, so filter the polygons by area, shape and a spectral index to keep what you want. For targeted results, use automatic=False, call set_image() once, and then prompt with points or boxes in map coordinates.

Segment Anything (SAM) can outline almost any object in a photograph, including roofs, fields, car parks and trees, without being trained on your data. The segment-geospatial package, samgeo, brings it to GIS. It reads GeoTIFFs, accepts prompts in longitude and latitude, and writes georeferenced rasters and vectors. In this tutorial, we run the smallest SAM model on a CPU over 2023 aerial imagery of a neighbourhood in Ames, Iowa. We test what it finds against OpenStreetMap’s building footprints.

You will learn how to:

1. stream a patch of NAIP aerial imagery from the Planetary Computer and save it as a GeoTIFF that SAM can read;
2. run automatic mask generation with SAM ViT-B on a CPU and measure how long it takes;
3. vectorise the masks to a GeoPackage and compute area, shape and NDVI for every segment;
4. filter the segments to find buildings, and score them against OpenStreetMap;
5. use point prompts and box prompts to segment specific objects.

All outputs below are real results from sam_aerial_segmentation.py, run on 17 October 2026 with Python 3.12, segment-geospatial 1.4.2, segment-anything 1.0 and PyTorch 2.14.1 (CPU). The machine was a desktop AMD Ryzen 9 7950X, limited to 8 threads, with no GPU. With the model and imagery cached, a full run takes about a minute; the first run also downloads the 375 MB model.

Install samgeo for CPU

PyTorch’s default wheels on Linux include CUDA and weigh several gigabytes. If you have no NVIDIA GPU, install the CPU build first from PyTorch’s own index, then install everything else:

python -m venv .venv-samgeo
.venv-samgeo\Scripts\activate          # on macOS or Linux: source .venv-samgeo/bin/activate
pip install torch==2.14.1 torchvision==0.29.1 --index-url https://download.pytorch.org/whl/cpu
pip install -r requirements.txt        # segment-geospatial, rasterio, geopandas, osmnx, ...

A separate virtual environment is a good idea: samgeo pulls in many dependencies, including opencv-python-headless and huggingface_hub. When you import it, segment-geospatial 1.4.2 prints a notice about installing SamGeo 3, which you can ignore.

Step 1: a patch of NAIP aerial imagery

SAM resizes every input so its longer side is 1024 pixels. A larger image is shrunk, and small objects are lost. We therefore cut a patch of exactly 1024 × 1024 pixels. At 0.5 m per pixel, that is a 512 × 512 m patch, which covers a primary school, about 120 houses and a group of apartment blocks.

The 2023 NAIP survey of Iowa was flown at 0.3 m. Resampling to 0.5 m by averaging fits more area into one SAM tile and slightly smooths the image. NAIP is distributed as quarter-quadrangle tiles, and this patch straddles two of them. We read each tile through a WarpedVRT onto our own grid and fill the empty pixels from the next tile:

CENTRE = (-93.684, 42.016)        # lon/lat: a school, houses and apartments in west Ames, Iowa
SIZE, RES = 1024, 0.5             # 1024 x 1024 pixels of 0.5 m = 512 x 512 m, exactly one SAM input tile

catalog = pystac_client.Client.open(STAC_URL, modifier=planetary_computer.sign_inplace)
nearby = list(catalog.search(collections=["naip"], intersects=box(*CENTRE, *CENTRE).buffer(0.004)).items())
year = max(i.datetime.year for i in nearby)               # the most recent NAIP survey
items = [i for i in nearby if i.datetime.year == year]    # NAIP tiles are quarter quads: the patch may span several
...
rgbn = np.zeros((4, SIZE, SIZE), dtype="uint8")
for it in items:                                      # mosaic: fill the empty pixels from each tile in turn
    with rasterio.open(it.assets["image"].href) as src, \
            WarpedVRT(src, crs=CRS, transform=TRANSFORM, width=SIZE, height=SIZE, src_nodata=0, nodata=0,
                      resampling=Resampling.average) as vrt:
        part = vrt.read()                             # 4 bands: red, green, blue, near-infrared
    empty = rgbn.max(axis=0) == 0
    rgbn[:, empty] = part[:, empty]

The script saves two files: a 4-band GeoTIFF (with near-infrared, for NDVI later) and a 3-band RGB copy for SAM, which expects ordinary colour images.

2 NAIP tiles from 2023 overlap the patch, e.g. ia_m_4209359_sw_15_030_20230902_20240104
flown 02 Sep 2023 | 0.3 m | EPSG:26915 | bands: red, green, blue, near-infrared
mosaic of 2 tiles; empty pixels: 0
patch: 1024 x 1024 pixels at 0.5 m = 512 x 512 m, resampled from 0.3 m; vegetated (NDVI > 0.2): 18%

Step 2: reference buildings from OpenStreetMap

To judge SAM’s results, we need building outlines that someone has drawn by hand. OpenStreetMap has them for this neighbourhood. Overpass, the usual OSM query server, was unavailable while we wrote this. For a box this small, the main OSM editing API serves the raw data directly, and OSMnx can read the file:

r = requests.get("https://api.openstreetmap.org/api/0.6/map", params={"bbox": ",".join(map(str, patch_ll))},
                 headers={"User-Agent": "technicalgis-tutorial"}, timeout=120)
osm_file.write_bytes(r.content)
osm = ox.features_from_xml(osm_file, tags={"building": True})
osm = osm[osm.geometry.type.isin(["Polygon", "MultiPolygon"])].to_crs(CRS).reset_index()
osm = osm[osm.within(patch.buffer(-1))].reset_index(drop=True)    # whole buildings only

The editing API is meant for editors, not bulk downloads. It refuses boxes larger than 0.25 square degrees, and its usage policy asks for light use. For larger areas, use ox.features_from_bbox through Overpass, or a Geofabrik extract.

146 building footprints entirely inside the patch; area median 212 m², range 7-6,449 m²
building
house                 121
semidetached_house     14
apartments              9
school                  1
yes                     1

OpenStreetMap is not a perfect reference. Several apartment buildings along the bottom edge of Figure 1 are missing from it, so SAM gets no credit for finding them. OSM outlines are also traced from different imagery, and NAIP is not a true orthophoto: roofs lean away from the camera by a metre or two. Even a perfect segmentation would not reach an IoU of 1.

Step 3: automatic mask generation with SAM ViT-B

SAM comes in three sizes: ViT-H (636 M parameters, 2.6 GB), ViT-L (308 M, 1.2 GB) and ViT-B (94 M, 375 MB). The smallest is the only practical choice on a CPU. In automatic mode, SAM places a grid of point prompts over the image, predicts masks for each, keeps the masks it rates as good and stable, and removes duplicates:

sam = SamGeo(model_type="vit_b", automatic=True, device="cpu", checkpoint_dir=str(CACHE),
             sam_kwargs={"points_per_side": 32,            # a 32 x 32 grid of point prompts = 1,024 prompts
                         "pred_iou_thresh": 0.86,          # keep masks the model rates as good...
                         "stability_score_thresh": 0.90,   # ...and that do not change much with the threshold
                         "min_mask_region_area": 40})      # remove holes and islands under 40 pixels (10 m²)
sam.generate(str(rgb_tif), output=str(masks_tif), foreground=True, unique=True)
sam.tiff_to_gpkg(str(masks_tif), str(masks_gpkg), simplify_tolerance=None)
segs = gpd.read_file(masks_gpkg).dissolve(by="value", as_index=False)   # one (multi)polygon per mask

unique=True gives every mask its own integer value in the output GeoTIFF. tiff_to_gpkg turns each run of equal pixels into a polygon, and dissolving by value gives one row per mask. Where masks overlap, samgeo paints the larger ones first, so smaller objects stay on top.

model loaded in 0.5 s: sam_vit_b_01ec64.pth, 375 MB, 94 M parameters
generate(): 306 masks in 47 s
vectorised to 306 polygons in sam_masks.gpkg; 34% of the patch is inside a mask

That “0.5 s” is with the checkpoint already cached. On the first run, samgeo downloaded the 375 MB file from Meta’s servers, and loading took 16 s on our connection. generate() took 45–47 s across our runs. Almost all of that time is the 1,024 point prompts plus the image encoder. Halving points_per_side to 16 cuts the prompts by four, at the cost of missing small objects. On a typical laptop CPU, expect these timings to be noticeably longer. A GPU brings them down to a few seconds.

Two things stand out in Figure 2. First, SAM segments objects of every kind: roofs, but also trees, the turning circles of the cul-de-sacs, car parks and the playing field. Second, it leaves most of the image (66%) unsegmented. Continuous backgrounds such as lawns and roads have no clear edges, and the default thresholds drop masks the model is unsure of.

Step 4: filter the masks to find buildings

Because SAM does not label anything, we add the labels ourselves with a few attributes per polygon: area, two shape measures and the mean NDVI from the near-infrared band.

segs["area_m2"] = segs.area
segs["compactness"] = 4 * np.pi * segs.area / segs.length ** 2         # 1 for a circle, ~0.785 for a square
mrr = segs.minimum_rotated_rectangle()
segs["rectangularity"] = segs.area / mrr.area                          # 1 for a perfect rectangle
segs["elongation"] = [e[-1] / max(e[0], 0.1) for e in edges]           # long side / short side of mrr
segs["ndvi"] = ndimage.mean(ndvi, mask_ids, segs["value"].to_numpy())

RULES = {"40 m² <= area <= 12,000 m²": segs["area_m2"].between(40, 12000),
         "NDVI < 0.15 (not vegetation)": segs["ndvi"] < 0.15,
         "rectangularity >= 0.75": segs["rectangularity"] >= 0.75,
         "elongation <= 4 (not a road)": segs["elongation"] <= 4}

Roofs are small to medium, not green, and fairly rectangular. Trees are green and round; roads and paths are long and thin. Each rule removes one kind of false positive:

  + 40 m² <= area <= 12,000 m²     ->  284 segments left
  + NDVI < 0.15 (not vegetation)   ->  211 segments left
  + rectangularity >= 0.75         ->  142 segments left
  + elongation <= 4 (not a road)   ->  139 segments left
precision: 81 of 139 kept segments match an OSM building with IoU >= 0.5 (58%)
recall before filtering: 107 of 146 OSM buildings have a SAM mask with IoU >= 0.5 (73%)
recall after filtering:  81 of 146 (55%); median IoU of matched buildings 0.73

The scores look modest, but the reasons for the errors in Figure 4 are instructive:

• Missing references. Some of the 58 orange “false positives” are real buildings that OpenStreetMap lacks, such as the apartment buildings along the bottom of the patch.
• One building, many masks. SAM split the school’s flat roof into several sections at its parapets and level changes. Each piece fails the IoU test against the single OSM outline. You can fix this by merging touching segments with similar colour, or by using box prompts (Step 6).
• Building-like non-buildings. A playing field, a car park and the school’s play areas are flat, not very green and rectangular, so they passed every rule. Add rules for them, such as a maximum area or a brightness test, or accept them and clean up by hand.
• Houses under trees. Many of the misses are houses partly hidden by trees in this mature suburb. The mask included the canopy, so the IoU fell below 0.5, or the NDVI rule dropped the mask.

Filtering cost 26 of the 107 buildings SAM had found (73% → 55% recall). Rules are a blunt tool. With a few hundred labelled examples, a small classifier, such as a random forest on these same attributes, usually separates the classes better.

Step 5: point prompts

Automatic mode answers “what objects are here?” Prompts answer “what is the object here?” In prompt mode, SAM runs its heavy image encoder once in set_image(). Each prompt then runs only the small mask decoder. With point_crs, samgeo converts longitude/latitude to pixel coordinates for you:

sp = SamGeo(model_type="vit_b", automatic=False, device="cpu", checkpoint_dir=str(CACHE))
sp.set_image(str(rgb_tif))                       # runs the image encoder once; every prompt after this is fast
POINTS = {"school roof": (-93.685575, 42.015863), "playing field": (-93.685261, 42.014694),
          "playground": (-93.684486, 42.015689), "apartment block": (-93.682507, 42.015250)}
for name, lonlat in POINTS.items():
    m, scores, _ = sp.predict(point_coords=[list(lonlat)], point_labels=1, point_crs="EPSG:4326",
                              multimask_output=True, return_results=True)

One click is ambiguous: does a point on a roof mean the roof section, the whole roof or the whole school? So SAM returns three candidate masks at different scales, each with a predicted quality score:

set_image(): image embedding computed in 3.9 s (model load included)
school roof        96 ms | 3 candidate masks: 5,698 m² (score 0.74), 1,320 m² (score 0.87), 637 m² (score 0.95)
playing field      34 ms | 3 candidate masks: 55,036 m² (score 0.63), 11,163 m² (score 0.72), 5,292 m² (score 0.95)
playground         33 ms | 3 candidate masks: 3,664 m² (score 0.70), 1,328 m² (score 0.90), 627 m² (score 0.62)
apartment block    31 ms | 3 candidate masks: 2,237 m² (score 0.69), 1,394 m² (score 0.69), 275 m² (score 0.81)
best mask for the school point written to school_point_mask.tif

The highest-scoring mask is SAM’s guess, not necessarily yours. For the school, the 0.95 mask is one roof section, but the 0.74 mask is nearly the whole building. For the apartment block, the best mask (0.81) is one 275 m² roof section, but the whole block is the 0.69 candidate. Add more points to resolve the ambiguity: positive points (label 1) on every part you want and negative points (label 0) on what you do not. The prediction after set_image took 31–96 ms per prompt, fast enough for interactive use.

Step 6: box prompts from OpenStreetMap

A box is a much stronger prompt than a point. If you already know roughly where the objects are, for example from an older map, from OSM or from an object detector, SAM can trace their exact outlines in the image. Here we draw each OSM footprint’s bounding box, padded by 2 m, and let SAM trace the building inside it:

PAD = 2.0                                       # metres of margin around each footprint's bounding box
for k, g in enumerate(osm.geometry):
    minx, miny, maxx, maxy = g.bounds
    c0, r0 = to_px * (minx - PAD, maxy + PAD)
    c1, r1 = to_px * (maxx + PAD, miny - PAD)
    m, score, _ = sp.predictor.predict(box=np.array([c0, r0, c1, r1]), multimask_output=False)

We call the underlying SamPredictor directly, with the box in pixel coordinates. In segment-geospatial 1.4.2, the multi-box path of SamGeo.predict() passes the boxes on incorrectly. The script then vectorises all the masks into sam_box_buildings.gpkg.

146 box prompts in 4.0 s (28 ms per box)
IoU with the OSM footprint: median 0.71, mean 0.71; >= 0.5: 99%, >= 0.75: 28%
          buildings  median_iou
area_m2
<50               1        0.19
100-300         131        0.70
300-1000          8        0.80
>1000             6        0.84

With a box, 145 of 146 buildings scored IoU ≥ 0.5, against 73% for automatic masks before filtering, and every one took 28 ms. The single failure is a 7 m² shed that is too small at 0.5 m. Large buildings score highest (0.84), because a metre of misalignment matters less on a big outline.

Why not higher than about 0.7 for houses? In the left panel of Figure 6, SAM’s outlines are slightly larger than OSM’s and shifted. SAM traces the roof as seen from the aeroplane, including eaves and the roof’s lean. OSM mappers trace the footprint, often from different imagery. A 212 m² house with an outline 1 m too wide on each side scores an IoU of about 0.77, so much of the gap is the reference rather than SAM.

This score is also optimistic about real-world use. The boxes came from the same OSM outlines we compare against, so they were already well placed. Boxes from an object detector would be looser.

Limitations and how to improve them

• No semantic labels. SAM segments; it does not classify. Every use needs a second step, such as rules, a classifier or prompts, to decide what each mask is.
• One tile at a time. SAM’s input is 1024 × 1024. For larger areas, samgeo’s tms_to_geotiff and split_raster help, but objects cut by tile edges need merging afterwards. Overlapping tiles and merging by IoU is the usual fix.
• Resolution. At 0.5 m a house is about 30 × 30 pixels, which SAM handles well. At Sentinel-2’s 10 m, a house is a pixel, and only fields, lakes and large buildings are segmentable. The 0.3 m original would help with sheds and garages, but four tiles of 1024 × 1024 would cover the same area.
• Off-nadir and shadows. NAIP is not a true orthophoto: tall buildings lean, and shadows can merge with roofs. Expect lower accuracy in city centres than in this suburb.
• CPU speed. ViT-B automatic mode took 47 s for 0.26 km². A town of 25 km² would take over an hour on 8 CPU threads. A GPU or fewer grid points speed it up. SAM 2 and fine-tuned remote sensing models are faster or more accurate options, but samgeo needs extra packages for them.
• Licences. SAM’s code and weights are Apache-2.0, so you can use them commercially. Check the licence of your imagery separately.

Frequently asked questions

Can Segment Anything be used on satellite imagery?

Yes, as long as the objects you want are at least a few dozen pixels across. SAM was trained on ordinary photographs, not satellite images, but it generalises well to sharp, high-resolution imagery such as aerial photography or 30–50 cm satellite imagery. It finds buildings, fields, water bodies, trees and car parks. On 10 m imagery, such as Sentinel-2, it only works for large objects such as fields and lakes.

What is samgeo?

samgeo is the import name of segment-geospatial, an open-source Python package by Qiusheng Wu. It wraps Meta’s Segment Anything models for geospatial data. It reads GeoTIFFs and web map tiles, accepts point and box prompts in map coordinates, and writes masks as GeoTIFF, GeoPackage or shapefile. It also offers text prompts through LangSAM and has interactive tools for Jupyter.

Can I run Segment Anything without a GPU?

Yes. Use the smallest model, ViT-B, and the CPU build of PyTorch. In this tutorial, automatic segmentation of a 1024 × 1024 tile took 47 seconds on 8 CPU threads, the image encoder took about 4 seconds, and each point or box prompt took about 30 ms. The larger ViT-H model is several times slower and is better suited to a GPU.

How do I convert SAM masks to polygons?

In samgeo, save the masks to a GeoTIFF with unique=True so that each mask has its own value. Then call sam.tiff_to_gpkg(masks_tif, "masks.gpkg") (or tiff_to_vector and tiff_to_shp). Dissolve the result by the value column with GeoPandas to get one polygon per mask. Without samgeo, rasterio.features.shapes does the same job.

Is SAM free for commercial use?

The original Segment Anything model, code and weights are released under the Apache-2.0 licence, which allows commercial use. SAM 2 is also Apache-2.0. segment-geospatial is MIT-licensed. Separate licences apply to other models samgeo can call, such as GroundingDINO for text prompts, and to your imagery.

Does SAM detect buildings automatically?

Not by itself. SAM finds the outlines of objects but does not know what they are. In this example, it outlined 73% of the buildings but also trees, car parks and a sports field. To extract buildings, filter the masks by area, shape and colour, or give SAM boxes from a detector or an existing map. For dedicated building extraction at scale, models trained on building footprints are usually more accurate.

Next steps

• Prompt SAM with your own data: draw boxes in QGIS, export them as a GeoPackage, and pass their bounds to predictor.predict(box=...).
• Fetch building footprints, roads and land use for any place with OSMnx in Street networks with OSMnx in Python.
• Learn STAC searches on the Planetary Computer in Measuring Lake Mead’s shrinking shoreline with Sentinel-2 and STAC.
• Use field boundaries to monitor crops through cloud in Sentinel-1 crop monitoring in Python, or map flood water in Sentinel-1 SAR flood mapping in Python.
• Classify every pixel into land cover classes instead of segmenting objects in Land cover classification with a Sentinel-2 random forest.
• Get building footprints and heights from lidar instead of imagery in Extract building footprints and heights from LiDAR in Python.

Data and licences

• NAIP aerial imagery: US Department of Agriculture, Farm Service Agency, National Agriculture Imagery Program, 2023. Public domain. Accessed from the Microsoft Planetary Computer on 17 October 2026.
• Building footprints: © OpenStreetMap contributors, available under the Open Database Licence. Downloaded through the OSM API.
• Segment Anything (ViT-B checkpoint sam_vit_b_01ec64.pth): Meta AI, Apache-2.0 licence. Kirillov et al. (2023), Segment Anything, arXiv:2304.02643.
• segment-geospatial: Qiusheng Wu, MIT licence. Wu, Q. and Osco, L. (2023), samgeo: A Python package for segmenting geospatial data with the Segment Anything Model (SAM), Journal of Open Source Software, 8(89), 5663.
• Code: sam_aerial_segmentation.py, tested with Python 3.12.0, PyTorch 2.14.1 (CPU), segment-geospatial 1.4.2, segment-anything 1.0, rasterio 1.5.2, geopandas 1.2.0 and osmnx 2.1.1.