Upload raster labels to a point layer attribute table

Reading raster labels

This is an update of a post I wrote way back in 2016. In GRASS GIS you can sample and upload raster values at point locations defined by a vector point layer using the v.what.rast function. If you are also interested in the raster category labels, you can use r.what, which lets you query a raster map on its category values and labels.

However, the results of r.what are written to a text file. Of course, you can use v.in.ascii to import the text file as a point vector layer in GRASS GIS. But back in 2016 I had some time at hand, so I wrote a small addon that queries the category values and labels of a raster layer for locations defined by a point layer. The category values and labels are directly written to the attribute table of a new point layer.

When I recently wanted to use it again, I had to solve a few Python 3 compatibility issues. Which made it a good time to submit it to the grass addon repository. For those running GRASS GIS 8.0 RC1 or newer, the addon is available as v.what.rast.label.

Installing the addon

You can install the add-on, as usual, using g.extension. To do so, in the menu, go to Settings > Addons extensions > install extension from addon. Alternatively, you can open the g.extension window by typing g.extension on the command line or in the console.

Installing the v.what.rast.label addon using the g.extension gui

Of course, you can install the addon directly by typing the following on the command line or in the console.

g.extension extension=v.what.rast.label operation=add

After installation, it is available in the toolbox under Addons (see the Tools tab). Or perhaps easier, just type in v.what.rast.label on the command line or in the console.

What it does

As input, you need one or more categorical raster layers with category labels, and a vector layer with points that represent the locations for which you want to sample the raster values and labels.

The function writes the raster categories and labels to the attribute table of a new point layer. The name of the column to which the values are written starts with ID_, followed by the name of the sampled raster layer. The name of the column with the raster labels is the same as the name of the sample raster layer.

Let’s use the North Carolina basic data set to illustrate the use of the addon. You can download the sample data from here (or you can install it directly from the menu if you use grass gis 8.0 or newer).

The data set includes a landuse raster layer with labels for each of the land-use categories. It also includes a vector layer points_of_interest (POI). The POIs cover a much larger area than the landuse map, so let’s first create a new point layer that only includes the POI within the bounds of the landuse map.

g.region raster=landuse
v.in.region output=regionbounds
v.select ainput=points_of_interest binput=regionbounds \
  output=POI_select operator=within

Example 1

Now, I would like to know the land-use categories for the POI locations. The screenshot below shows the use of the gui.

First, the function makes a copy of the input point layer. By default, it does not include the original attribute table. Next, the land-use category values and labels at the POI locations are uploaded to the attribute table of this new point layer POIex1.

Example of using v.what.rast.label, uploading the raster values and labels of the land-use map to the point layer’s attribute table. Note that the corresponding code is also provided (see the red arrow). A great way to familiarize oneself with the code from the comfort of a gui.

The attribute table of the resulting point layer contains three columns; the cat column with the features id, the column landuse_ID with the categories, and the column landuse with the category labels of the landuse raster layer.

The attribute table of the point layer POIex1.

Example 2

You can upload the categories and category labels of more than one raster at once. In the example below I use the command line to get the raster categories and labels of the landuse and geology maps for the POI locations.

v.what.rast.label vector=POI_select raster=landuse,geology output=POIex2

We can use the db.columns function to list the columns in the attribute table of the newly created POIex2. It shows that besides the cat column, we now have four columns holding the category values and labels of the landuse and geology map.

db.columns table=POIex2
cat
landuse_ID
landuse
geology_ID
geology

Example 3

It is also possible to upload the values of a raster layer without category labels, using the raster2 parameter. In the example below, the landuse, geology, and elevation maps are queried. The latter has no labels. This time I use the Modules module from the pyGRASS scripting library to run the command in Python. As you can see, it is the code you would use on the command line, but wrapped in a Module call.

from grass.pygrass.modules import Module
Module("v.what.rast.label", vector="POI_select", raster=["landuse", "geology"],
       raster2="elevation", output="POIex3")
Module("db.columns", table="POIex3")
cat
landuse_ID
landuse
geology_ID
geology
elevation

Example 4

There are two additional options. The first is to copy the columns of the attribute table of the input point vector layer to the attribute table of the new point layer. Incidentally, being able to do this easily was the reason I wrote this addon. You can enable this option using the -o flag.

The other option is to write the coordinates of the POI’s to the attribute table of the new point layer. I can’t even remember why I included that option. I guess to make it easier to export it as a csv file and import it in another program. Anyway, if you want to include the coordinates, use the -c flag.

from grass.pygrass.modules import Module
Module("v.what.rast.label", flags= ["o", "c"], vector="POI_select",
       raster="landuse", output="POIex4")

Printing the columns names results in a long list. Instead, we can use PIPE functionality (the python equivalent of bash pipe) to get the output, replace the line endings with a comma, and print the names on one line to the console.

from subprocess import PIPE
cols = Module("db.columns", table="POIex4", stdout_=PIPE).outputs.stdout
colsl = cols.split("\r\n")
print(*colsl, sep=", ")
cat, feature_id, featurenam, class, st_alpha, st_num, county, county_num, primlat_dm, primlon_dm, primlatdec, primlondec, srclat_dms, srclon_dms, 
srclatdec, srclondec, elev_m, map_name, x, y, landuse_ID, landuse, 

And there you have it, all the columns of the input vector layer, the two columns with the coordinates, and the last two columns with the land-use categories and category labels.



Paulo van Breugel
Paulo van Breugel
Lecturer & researcher

My interests range from biodiversity and ecology to spatial data analysis. I am also what one could describe as a lifelong learner; I enjoy to learn new things and widen my horizon, both professionally and personally.

comments powered by Disqus
Next
Previous