Segmenting remote sensing imagery with point prompts¶

This notebook shows how to generate object masks from point prompts with the Segment Anything Model 2 (SAM 2).

Make sure you use GPU runtime for this notebook. For Google Colab, go to Runtime -> Change runtime type and select GPU as the hardware accelerator.

Install dependencies¶

Uncomment and run the following cell to install the required dependencies.

1	`# %pip install -U segment-geospatial`

Import libraries¶

import leafmap
from samgeo import SamGeo2
from samgeo.common import regularize

Create an interactive map¶

m = leafmap.Map(center=[47.653287, -117.588070], zoom=16, height="800px")
m.add_basemap("Satellite")
m

Download a sample image¶

Pan and zoom the map to select the area of interest. Use the draw tools to draw a polygon or rectangle on the map. If no geometry is drawn, the default bounding box will be used.

if m.user_roi is not None:
    bbox = m.user_roi_bounds()
else:
    bbox = [-117.6029, 47.65, -117.5936, 47.6563]

image = "satellite.tif"
leafmap.map_tiles_to_geotiff(
    output=image, bbox=bbox, zoom=18, source="Satellite", overwrite=True
)

You can also use your own image. Uncomment and run the following cell to use your own image.

1	`# image = '/path/to/your/own/image.tif'`

Display the downloaded image on the map.

m.layers[-1].visible = False
m.add_raster(image, layer_name="Image")
m

Initialize SAM class¶

Set automatic=False to enable the SAM2ImagePredictor.

sam = SamGeo2(
    model_id="sam2-hiera-large",
    automatic=False,
)

Specify the image to segment.

sam.set_image(image)

Segment the image¶

Use the predict_by_points() method to segment the image with specified point coordinates. You can use the draw tools to add place markers on the map. If no point is added, the default sample points will be used.

if m.user_rois is not None:
    point_coords_batch = m.user_rois
else:
    point_coords_batch = [
        [-117.599896, 47.655345],
        [-117.59992, 47.655167],
        [-117.599928, 47.654974],
        [-117.599518, 47.655337],
    ]

Segment the objects using the point prompts and save the output masks.

sam.predict_by_points(
    point_coords_batch=point_coords_batch,
    point_crs="EPSG:4326",
    output="mask.tif",
    dtype="uint8",
)

Display the result¶

Add the segmented image to the map.

m.add_raster("mask.tif", cmap="viridis", nodata=0, opacity=0.7, layer_name="Mask")
m

Use an existing vector dataset as points prompts¶

Alternatively, you can specify a file path or HTTP URL to a vector dataset containing point geometries.

geojson = "https://github.com/opengeos/datasets/releases/download/places/wa_building_centroids.geojson"

Display the vector data on the map.

m = leafmap.Map()
m.add_raster(image, layer_name="Image")
m.add_circle_markers_from_xy(
    geojson, radius=3, color="red", fill_color="yellow", fill_opacity=0.8
)
m

Segment image with a vector dataset¶

Segment the image using the specified file path to the vector dataset.

output_masks = "building_masks.tif"

sam.predict_by_points(
    point_coords_batch=geojson,
    point_crs="EPSG:4326",
    output=output_masks,
    dtype="uint8",
    multimask_output=False,
)

Display the segmented masks on the map.

m.add_raster(
    output_masks, cmap="jet", nodata=0, opacity=0.7, layer_name="Building masks"
)
m

Clean up the result¶

Remove small objects from the segmented masks, fill holes, and compute geometric properties.

out_vector = "building_vector.geojson"
out_image = "buildings.tif"

array, gdf = sam.region_groups(
    output_masks, min_size=200, out_vector=out_vector, out_image=out_image
)

gdf.head()

Regularize building footprints¶

Regularize the building footprints using the regularize() method.

output_regularized = "building_regularized.geojson"
regularize(out_vector, output_regularized)

Display the regularized building footprints on the map.

m = leafmap.Map()
m.add_raster(image, layer_name="Image")
style = {
    "color": "#ffff00",
    "weight": 2,
    "fillColor": "#7c4185",
    "fillOpacity": 0,
}
m.add_raster(out_image, cmap="tab20", opacity=0.7, nodata=0, layer_name="Buildings")
m.add_vector(
    output_regularized, style=style, layer_name="Building regularized", info_mode=None
)
m

Interactive segmentation¶

sam.show_map()