Building Raster Attribute Table (RAT) from GeoPandas dataframe with rasterio or GDAL

Question

I am trying to figure out how to generate a RAT for a 1-band raster with rasterio or GDAL given a GeoPandas dataframe.

My GeoPandas dataframe has 4 columns (X, Y, PERC_VALUE, LANDUSE_LABEL) which I would need to use to build my raster as so:

• X, Y defines the coordinates of the pixel centroid
• PERC_VALUE is the value which has to be stored in the pixels
• LANDUSE_LABEL has information associated to the pixels which I need to store in the RAT.

The Pandas dataframe that I am using to build the GeoPandas dataframe looks like this:

landuse value   Y   X
FOREST  0.0005  1386050     4559950
CITY    0.0001  1386250     4559650
FOREST  0.0009  1386250     4559750
FOREST  0.0025  1386250     4559850
CITY    0.0059  1386250     4561950
...

I am able to create a raster to store the coverage percentage values at each corresponding X,Y pixel (code is shown at the bottom of the question). What I don't know is how can I associate each pixel to the corresponding LANDUSE_LABEL information in a RAT?

The idea is that when I open the final result with (e.g.) QGIS (if it can read RATs), I would be able to symbolize the raster per landuse label, meanwhile being able to access the percentage values of the pixels.

Is this possible with rasterio or GDAL and if so how? Do I have to use a specific raster extension for this aim (I am currently using GeoTIFF)?

Code I am using to write the PERC_VALUE with rasterio as GeoTIFF:

# convert Pandas df to GeoPandas DF (Points)
gdf = GeoDataFrame(
    dataframe.drop(['X', 'Y'], axis=1),
    crs={'init': 'epsg:3035'},
    geometry=[Point(xy) for xy in zip(dataframe.X.astype(int), dataframe.Y.astype(int))])

# write the raster
with rasterio.open(r'output.tiff',
               'w+',
               driver= 'GTiff',
               height= arr.shape[0],
               width= arr.shape[1],
               count= 1,
               dtype= numpy.dtype('float32'),
               crs= SR,
               transform= transform) as out:

    out_arr = out.read(1)

    # this is where we create a generator of geom, value pairs to use in rasterizing
    shapes = ((geom,value) for geom, value in zip(gdf['geometry'], gdf['value']))
    burned = features.rasterize(shapes=shapes, fill=0, out=out_arr, transform=transform)
    out.write(burned, 1)

    out.close()

Answer 1

I took inspiration from the code at this GitHub discussion to solve my problem and write the RAT.

The result is an .aux.xml file with the same filename as the tif produced by the code in my question.

Unfortunately, QGIS is not capable of reading the RATs, but at least I can use the xml to find out the correspondence between pixel values and their labels and symbolize the map manually. There is an open Feature Request in QGIS project to add RAT reading in QGIS but it is still a proposal: https://github.com/qgis/QGIS/issues/22427.

Here is the code (using GDAL):

# https://gdal.org/python/osgeo.gdal.RasterAttributeTable-class.html
# https://gdal.org/python/osgeo.gdalconst-module.html
ds = gdal.Open(OUTPUT_RASTER)
rb = ds.GetRasterBand(1)

# Get unique values in the band
u = np.unique(rb.ReadAsArray())
u = u.tolist()

# https://chrisalbon.com/python/data_wrangling/pandas_list_unique_values_in_column/
# extract list of unique 'landuse' values
landuse = dataframe.landuse.unique()

# Create and populate the RAT
rat = gdal.RasterAttributeTable()

rat.CreateColumn('VALUE', gdal.GFT_Integer, gdal.GFU_Generic)
rat.CreateColumn('LANDUSE', gdal.GFT_String, gdal.GFU_Generic)

for i in range(len(u)):
    rat.SetValueAsInt(i, 0, u[i])
    rat.SetValueAsString(i, 1, landuse[i])

# Associate with the band
rb.SetDefaultRAT(rat)

# Close the dataset and persist the RAT
ds = None