Creating GeoPandas GeoDataFrame from rasterio features

Question

I try to create a GeoDataFrame of paches of the same values in a binary raster using rasterio and GeoPandas. data_file is a binary TIFF raster with 0 where there is no deforestation and 1 if there is a deforestation event.

I tried the following :

# read the data and create the shapes 
with rasterio.open(data_file) as f:
    
    data = data.astype('int16')
    shapes = rasterio.features.shapes(data)

# read the shapes as separate lists
deforestation = []
geometry = []
for shape, value in shapes:
    deforestation.append(value)
    geometry.append(shape)

# build the gdf object over the two lists
gdf = gpd.GeoDataFrame(
    {'deforestation': deforestation, 'geometry': geometry },
    crs="EPSG:4326"
)

I get the following warning :

/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:9:
FutureWarning: Assigning CRS to a GeoDataFrame without a geometry column is now deprecated and will not be supported in the future.
if name == 'main':

From what I understand, the shapes GeoJSON-like dicts are not interpreted as shapely shape by GeoPandas, so I end up with a Gdf without 'geometry' column which cannot be used to create a shapefile.

What is the more efficient way to transform them into shapely polygon?

Answer 1

Try building a shapely Polygon from the geojson-like dicts returned by rasterio.features.shapes using the shapely.geometry.shape function.

Completely untested example:

import geopandas as gpd
import rasterio
from shapely.geometry import shape

# read the data and create the shapes 
with rasterio.open(data_file) as f:
    data = f.read(1)
    data = data.astype('int16')
    shapes = rasterio.features.shapes(data, transform=f.transform)

# read the shapes as separate lists
deforestation = []
geometry = []
for shapedict, value in shapes:
    deforestation.append(value)
    geometry.append(shape(shapedict))

# build the gdf object over the two lists
gdf = gpd.GeoDataFrame(
    {'deforestation': deforestation, 'geometry': geometry },
    crs="EPSG:4326"
)

Answer 2

With geocube 0.4+:

import rioxarray
from geocube.vector import vectorize

data = rioxarray.open_rasterio(data_file).squeeze().astype('int16')
data.name = "deforestation"
gdf = vectorize(data)