0

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.

  • Can you provide a bit more information on what 'chessboard segmentation' is? – Alex Leith Mar 17 at 2:21
  • Have you seen the scikit implementation: scikit-image.org/docs/dev/api/… – Aaron Mar 17 at 3:28
  • @Aaron I had not, thanks :) – somewheresouth Mar 17 at 5:32
  • @AlexLeith it's a segmentation that picks out pixel edges, essentially just chopping your image up by areas of a certain size. Ideally in this case grabbing pixel edges. – somewheresouth Mar 17 at 5:34
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:

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:

Example output

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.