3

I have a geotiff file.

import xarray as xr
urbanData = xr.open_rasterio('myGeotiff.tif')
plt.imshow(urbanData)

Here the link to the file.

I can convert the file as a dataframe with coordinates as points

ur  = xr.DataArray(urbanData, name='myData')
ur  = ur.to_dataframe().reset_index() 
gdfur  = gpd.GeoDataFrame(ur, geometry=gpd.points_from_xy(ur.x, ur.y))

However I would like to get a dataframe that contains the geometry of the pixels as polygons and not as points. Is it possible?

emax
  • 6,965
  • 19
  • 74
  • 141
  • 1
    What does "geometry of the pixels as polygons" even mean? A GeoTIFF is a raster file, not a vector format, so I doubt that it would store polygon data. – jjramsey May 11 '21 at 13:34
  • @jjramsey that's the point. Is it possible to convert the raster file as a vector one? – emax May 11 '21 at 13:36
  • 1
    That's dicey. For *art*, there are tools for making vector files from raster files by doing boundary tracing and storing the path of the trace as a vector. However, the raster files in that scenario would be relatively simple, usually containing depictions of curves or regions of solid color. A decent candidate for such vector tracing would be a scan of a black & white cartoon drawing, not the sort of content that would usually be in a GeoTIFF. – jjramsey May 11 '21 at 13:47
  • 1
    Sounds to me like you want to polygonize your raster. If so, [this](https://gis.stackexchange.com/questions/187877/how-to-polygonize-raster-to-shapely-polygons) could be a good start – Val May 11 '21 at 14:09

2 Answers2

2

Somewhat to my surprise, I haven't really found a package which wrap rasterio.features to take DataArrays and produce GeoDataFrames.

These might be very useful though:

https://corteva.github.io/geocube/stable/

https://corteva.github.io/rioxarray/stable/

I generally use something like this:

import affine
import geopandas as gpd
import rasterio.features
import xarray as xr
import shapely.geometry as sg


def polygonize(da: xr.DataArray) -> gpd.GeoDataFrame:
    """
    Polygonize a 2D-DataArray into a GeoDataFrame of polygons.

    Parameters
    ----------
    da : xr.DataArray

    Returns
    -------
    polygonized : geopandas.GeoDataFrame
    """
    if da.dims != ("y", "x"):
        raise ValueError('Dimensions must be ("y", "x")')

    values = da.values
    transform = da.attrs.get("transform", None)
    if transform is None:
        raise ValueError("transform is required in da.attrs")
    transform = affine.Affine(*transform)
    shapes = rasterio.features.shapes(values, transform=transform)

    geometries = []
    colvalues = []
    for (geom, colval) in shapes:
        geometries.append(sg.Polygon(geom["coordinates"][0]))
        colvalues.append(colval)

    gdf = gpd.GeoDataFrame({"value": colvalues, "geometry": geometries})
    gdf.crs = da.attrs.get("crs")
    return gdf

Note that you should squeeze the band dimensions from your xarray first to make it 2D, after reading it with xr.open_rasterio:

urbanData = xr.open_rasterio('myGeotiff.tif').squeeze('band', drop=True)
Huite Bootsma
  • 451
  • 2
  • 6
  • 1
    this method doesnot create polygons with holes. if you want to do so, use `geometries.append(shape(geom))` instead of `geometries.append(sg.Polygon(geom["coordinates"][0]))|` – Hugo Roussaffa Dec 27 '21 at 03:46
0

With geocube 0.4+:

import rioxarray
from geocube.vector import vectorize

data = rioxarray.open_rasterio("myGeotiff.tif").squeeze()
data.name = "myData"
gdf = vectorize(data)
snowman2
  • 646
  • 4
  • 11