How to chessboard segment landsat pixel edges in python?

Question

Does anyone know how to do a chessboard segmentation (in the same way that ecognition does) that will select all the edges of a Landsat pixel in python/a Jupyter notebook? Given that I can use the notebook to retrieve the imagery as a multi-band xarray, does anyone have a good recipe using scikit-image or similar to take the imagery array and run a chessboard segmentation on it, outputting the pixel edges as a shapefile? I want every pixel edge as opposed to areas of similarity.

I can do this in ecognition, and create a workflow. I'm looking for a python equivalent to the chessboard segmentation workflow in ecognition that I can apply on a large scale for multiple times. I'm not doing it in ecognition because I need to do large amounts of timesteps and large amounts of spatial area, and it's impractical in human-time-cost to do so for thousands of images.

Answer 1

The code below will extract each Landsat raster pixel as a polygon vector in geopandas.GeoDataFrame format. input_data represents the data in the input xarray dataset; you'll need to substitute in the dataset's geotransform and CRS data for input_transform and input_crs:

import numpy as np
import geopandas as gpd
import rasterio.features
import matplotlib.pyplot as plt
from shapely.geometry import shape

# Input array to segment and vectorise
input_array = np.random.rand(5, 5)
input_transform = (30, 0, 0, 0, -30, 0)  # replace with xarray transform
input_crs = 'EPSG:4326'  # replace with xarray CRS

plt.imshow(input_array)

Example input:

# Create array with a unique value per cell
unique_pixels = np.arange(input_array.size).reshape(input_array.shape)

# Vectorise each unique feature in array
vectors = rasterio.features.shapes(source=unique_pixels.astype(np.int16), 
                                     transform=input_transform)

# Extract polygons and values from generator
vectors = list(vectors)
values = [value for polygon, value in vectors]
polygons = [shape(polygon) for polygon, value in vectors]

# Create a geopandas dataframe populated with the polygon shapes
gdf = gpd.GeoDataFrame(data={'id': values},
                       geometry=polygons,
                       crs=input_crs)

# Plot vectors by attribute
gdf.plot(column='id', edgecolor='black')

Example output: