2

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()
  • AFAIK, QGIS does not read/write RATs. I've confirmed by loading an Arc-created integer raster that had an attribute table into QGIS. QGIS could display the raster along with its cell values, but not the associated tabular attributes. This is a shame; I look forward to the day when QGIS is RAT-aware! – Stu Smith Dec 10 '19 at 20:32
  • I found it myself after being able to succesfully write the RAT. This is very sad, but still, I am glad QGIS is there! I will answer my own question with the code for writing the RAT. – umbe1987 Dec 11 '19 at 8:47
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

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.