Unit 13 - PyGRASS Raster Access

PyGRASS allows directly accessing native GRASS raster and vector maps in the sense of Python objects. This unit shows how to deal with GRASS raster data by PyGRASS API, see Unit 14 - PyGRASS Vector Access for vector data.

Raster data

Raster map can be treated by RasterRow for reading raster data row by row. There is also RasterSegment which allows reading data by user-defined segments (tiles).

from grass.pygrass.raster import RasterRow
ndvi = RasterRow('ndvi')

Raster map is open by open() method. Basic information like number of columns and rows, min and max values, range printed.

ndvi.open()
print (ndvi.info.cols, ndvi.info.rows)
min, max = ndvi.info.range
print (min, max)
print (max - min)

Fig. 85 Running PyGRASS code from Python tab of Layer Manager.

Don’t forget to close the raster map at the end.

ndvi.close()

Raster statistics example

A simple PyGRASS script for computing basic univariate raster statistics below.

#!/usr/bin/env python

import numpy as np
from grass.pygrass.raster import RasterRow

ndvi = RasterRow('ndvi')
ndvi.open()

min = max = None
count = ncount = 0
for row in ndvi:
    for value in row:
        if np.isnan(value):
            ncount += 1
        else:
            if min is None:
                min = max = value
            else:
                if min > value:
                    min = value
                elif max < value:
                    max = value
        count += 1

ndvi.close()

print ("min={0:.6f} max={1:.6f} count={2} (no-data: {3})".format(
    min, max, count, ncount)
)

Tip

Compare computation speed of your simple script with C-based r.univar module…

Fig. 86 PyGRASS script and r.univar comparision.

Sample script to download: ndvi-pygrass-univar.py

Writing raster data

It is possible to write raster map using PyGRASS, it is possible to start from existing maps or from scratch. In the next example we are going to create ndvi map from the existing Sentinel bands

import numpy

b04 = RasterRow('L2A_T32VNM_20170705T105031_B04_10m')
b04.open()
b08 = RasterRow('L2A_T32VNM_20170705T105031_B08_10m')
b08.open()
ndvi = RasterRow('ndvi_pyrass')
ndvi.open('w', mtype='FCELL', overwrite=True)

for i in range(len(b04)):
    row_b04 = b04[i]
    row_b08 = b08[i]
    rowb04 = row_b04.astype(numpy.float32)
    rowb08 = row_b08.astype(numpy.float32)
    row_new = (rowb08 - rowb04) / (rowb08 + rowb04)
    ndvi.put_row(row_new)

ndvi.close()
b04.close()
b08.close()

Tip

To create a raster from scratch you need to remember to get the right number of rows and columns and create an empty row (Buffer object in PyGRASS)

newscratch = RasterRow('newscratch')
newscratch.open('w', overwrite=True)

# get computational region info
from grass.pygrass.gis.region import Region
reg = Region()

# import bufffer and create empty row
from grass.pygrass.raster.buffer import Buffer
newrow = Buffer((reg.cols,), mtype='CELL')

# we create a raster to fill all the GRASS GIS region
for r in range(reg.rows):
    newrow[:] = np.random.random_integers(0, 1000, size=newrow.size)
    newscratch.put_row(newrow)

newscratch.close()