Unit 13 - PyGRASS Raster Access

PyGRASS allows to access GRASS raster and vector maps in the sense of Python objects. This unit shows how to deal with GRASS raster data by PyGRASS API.

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')
ndvi.open()

Raster map is open by open() method. The sample code also prints basic information as number of columns and rows, min and max values, range.

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

Fig. 81 Running PyGRASS code from Python tab.

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 python3

import numpy

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 numpy.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)
)

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

Optimized version based on NumPy library:

#!/usr/bin/env python3

import numpy

from grass.pygrass.raster import RasterRow

ndvi = RasterRow('ndvi')
ndvi.open()
array = numpy.array(ndvi)
ndvi.close()

print("min={0:.6f} max={1:.6f} count={2} (no-data: {3})".format(
    numpy.nanmin(array), numpy.nanmax(array), array.size,
    numpy.count_nonzero(numpy.isnan(array)))
)

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

Úkol

Compare computation speed of two sample scripts with C-based r.univar module.

Writing raster data

PyGRASS allows also writing raster data. In the example below a NDVI map will be computed from Sentinel-2 red and near-infrated channels similarly as in Unit 05 - Raster processing.

#!/usr/bin/env python3

import numpy

from grass.pygrass.raster import RasterRow

b04 = RasterRow('L2A_T32UPB_20170706T102021_B04_10m')
b04.open()
b08 = RasterRow('L2A_T32UPB_20170706T102021_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()

Sample script to download: ndvi-pygrass.py

Úkol

Create a raster map from scratch. Raster row is represented by Buffer object in PyGRASS.

#!/usr/bin/env python3

import numpy

from grass.pygrass.raster import RasterRow

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

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

# import buffer 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[:] = numpy.random.randint(0, 1001, size=newrow.size)
    newscratch.put_row(newrow)
          
newscratch.close()

Fig. 82 Example of created raster map.

Sample script to download: pygrass-write-raster.py