6

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?

3
  • It is not an error but a warning
    – gene
    Nov 17, 2020 at 10:29
  • data_file is a binary raster, I want to create patches of value 1. I know that it's a warning but I get a gdf without 'geometry' column so it didn't realize that it was polygons (I eddited the question to match your comments) Nov 17, 2020 at 11:12
  • @user2856 can you stop editing my question. SO meta says that errors of more that 50 of width should be written in quote block to improve the quality of reading on small devices. (one ref out of many) Nov 22, 2020 at 11:16

2 Answers 2

9

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"
)
0
6

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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