Extracting shapes from GeoTIFF using rasterio returns excessive number of polygons

Question

I am using the shapes interface of rasterio to extract polygons from a GeoTIFF image. rasterio returns an excessive number of polygons which is not expected. It is taking a long time to extract a relatively very simple raster. I did additional processing of the resulting shapefiles, namely Dissolve and Multipart to a single part, and eventually, get the correct polygons.

Can any experts in rasterio can shed some light on what is happening? For the sample raster, I am expecting fewer polygons whereas rasterio shapes are 77562 in number. File shared here

import pandas as pd
import geopandas
import matplotlib.pyplot as plt
import numpy as np
import os

import rasterio
from rasterio.features import shapes
mask = None
results=None
with rasterio.Env():
    with rasterio.open('/content/small_area.tif') as src:
        image = src.read(1) # first band
        results = (
     
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.transform)))
        
geoms=list(results)        

import geopandas as gp

gdf = geopandas.GeoDataFrame.from_features(geoms) #Convert to Geopandas dataframe

gdf.insert(0, 'New_ID', range(0, 0 + len(gdf)))

gdf.head()

gdf.set_crs(epsg=3857, inplace=True)
gdf = gdf.to_crs(epsg=3857)
gdf.set_geometry(col='geometry', inplace=True)
gdf.to_file("small_area.shp")

Answer 1

I am getting exactly the same number of polygons (77562) with GDAL utility gdal_polygonize.

gdal_polygonize -f gpkg small_area.tif small.gpkg

If I open the resulting vector layer with QGIS and apply random color classification by the DN field that contains the pixel values from the source raster it is easy to see why there are so many polygons in the result. This is a magnified view of the polygons and the colorful boxes are presenting a diagonal line in the source map.

The straight lines which look like solid black in the map are actually rendered with antialiasing and pixels have many different color values.

The rasterio.features.shapes function is documented in https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html

Get shapes and values of connected regions in a dataset or array.

The antialiased colors are not connected and therefore the result contain large number of single pixel polygons.

By the screen capture it looks like that there were connected blues divided into different polygons but actually there are two different dark blue colors and they have values 32 and 241 in the data.

The function works as it is supposed to work. You can try to re-process the source map with another rasterio function rasterio.features.sieve

Replace small polygons in source with value of their largest neighbor.

Polygons are found for each set of neighboring pixels of the same value.