GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
GRASS Raster Library

Table of Contents

by GRASS Development Team (https://grass.osgeo.org)

TODO: Needs to be cleaned up. The upper GRASS 4.x and the lower GRASS 5.x/6.x parts need to me merged and updated to GRASS 7

GRASS Raster Map Processing

This library provides a suite of routines which process raster map file. The processing of raster map files consists of determining which raster map file or files are to be processed (specified on the module command line), locating the raster map file in the database, opening the raster files, dynamically allocating i/o buffers, reading or writing the raster files, closing the raster files, and creating support files for newly created raster maps.

Raster map data can be of type CELL, FCELL or DCELL, they are defined in "gis.h". CELL is a 32-bit signed integer, FCELL is an IEEE single-precision floating-point, and DCELL is an IEEE double-precision floating-point. 3D rasters (grid3d) is treated as DCELL (see related library).

Finding Raster Files in the Database

GRASS allows the user to specify raster map names (or any other GIS database file) either as a simple unqualified name, such as "soils", or in case of input raster maps also as a fully qualified name, such as "soils@mapset", where mapset is the mapset where the raster map is to be found. Often only the unqualified raster map name is provided on the command line and searched in all mapsets indicated in the mapset search path (SEARCH_PATH file in the actual mapset; managed with g.mapsets command).

The following routines search the database for raster map files:

Looks for the raster file in the database. If found, the mapset where the raster file lives is returned. If not found, the NULL pointer is returned. If the user specifies a fully qualified raster map name which exists, then G_find_raster() modifies name by removing the "@<I>mapset</I>".

For example, to find a raster map in all mapsets listed the mapset search path:

char name[GNAME_MAX];
char *mapset;
if ((mapset = G_find_raster(name,"")) == NULL)
/* not found */

To check that the raster map exists in the current mapset:

char name[GNAME_MAX];
if (G_find_raster(name, G_mapset()) == NULL)
/* not found */

Opening an Existing Raster File

The following routine opens the raster map file for reading.

This routine opens the raster map in given mapset for reading. A nonnegative file descriptor is returned if the open is successful. Otherwise a diagnostic message is printed and a negative value is returned. This routine does quite a bit of work. Since GRASS users expect that all raster maps will be resampled into the current region (nearest neighbor method), the resampling index for the raster map is prepared by this routine after the map is opened. The resampling is based on the active module region. Preparation required for reading the various raster map formats (CELL, FCELL, DCELL) is also done.

Creating and Opening New Raster Files

The following routines create a new raster map in the current mapset and open it for writing. G_legal_filename() should be called first to make sure that raster map name is a valid name.

Note: It is not an error for raster map to already exist. New raster maps are actually created as temporary files and moved into the cell or fcell directory when closed. This allows an existing raster map to be read at the same time that it is being rewritten. G_find_raster() could be used to see if raster map already exists.

Warning: However, there is a subtle trap. The temporary file, which is created using G_tempfile(), is named using the current process id. If the new raster file is opened by a parent process which exits after creating a child process using fork(), the raster file may never get created since the temporary file would be associated with the parent process, not the child. GRASS management automatically removes temporary files associated with processes that are no longer running. If fork() must be used, the safest course of action is to create the child first, then open the raster file (see the discussion under G_tempfile() for more details).

Creates and opens the raster map for writing by Rast_put_row() which writes the file row by row in sequential order. The raster file data will be compressed as it is written. A nonnegative file descriptor is returned if the open is successful. Otherwise a diagnostic message is printed and a negative value is returned.

Creates and opens the raster file for writing by Rast_put_row() which writes the file row by row in sequential order. The raster file will be in uncompressed format when closed. A nonnegative file descriptor is returned if the open is successful. Otherwise a warning message is printed on stderr and a negative value is returned.

General use of this routine is not recommended. This routine is provided so the r.compress module can create uncompressed raster files.

Allocating Raster I/O Buffers

Since there is no predefined limit for the number of columns in the region, buffers which are used for reading and writing raster data must be dynamically allocated.

This routine allocates a buffer of type CELL/FCELL/DCELL just large enough to hold one row of raster data (based on the number of columns in the active region).

int input_fd;
char inmap;
RASTER_MAP_TYPE data_type;
input_fd = Rast_open_old(inmap, "");
data_type = Rast_get_map_type(input_fd);
cell = Rast_allocate_buf(data_type);

FIXME 7/2009: next still true?

If larger buffers are required, the routine G_malloc() can be used.

If sufficient memory is not available, an error message is printed and exit() is called.

  • G_zero_buf()

This routines assigns each member of the raster buffer array to zero. It assumes that the buffer has been allocated using Rast_allocate_c_buf().

Reading Raster Files

Raster data can be thought of as a two-dimensional matrix. The routines described below read one full row of the matrix. It should be understood, however, that the number of rows and columns in the matrix is determined by the region, not the raster file itself. Raster data is always read resampled (nearest neighbor) into the region. This allows the user to specify the coverage of the database during analyses. It also allows databases to consist of raster files which do not cover exactly the same area, or do not have the same grid cell resolution. When raster files are resampled into the region, they all "look" the same.

Note: The rows and columns are specified "C style", i.e., starting with 0.

This routine reads the specified row from the raster file open on file descriptor (as returned by Rast_open_old()) into the buffer. The buffer must be dynamically allocated large enough to hold one full row of raster data. It can be allocated using Rast_allocate_buf(). This routine prints a diagnostic message and returns -1 if there is an error reading the raster file. Otherwise a nonnegative value is returned.

This routine reads the specified row from the raster file open on file descriptor into the buffer like Rast_get_row() does. The difference is that masking is suppressed. If the user has a mask set, Rast_get_row() will apply the mask but Rast_get_row_nomask() will ignore it. This routine prints a diagnostic message and returns -1 if there is an error reading the raster file. Otherwise a nonnegative value is returned.

Note: Ignoring the mask is not generally acceptable. Users expect the mask to be applied. However, in some cases ignoring the mask is justified. For example, the GRASS modules r.describe, which reads the raster file directly to report all data values in a raster file, and r.slope.aspect, which produces slope and aspect from elevation, ignore both the mask and the region. However, the number of GRASS modules which do this should be minimal. See Mask for more information about the mask.

Writing Raster Files

This routine writes one row of raster data from buffer to the raster file open on file descriptor. The raster file must have been opened with Rast_open_new(). The buffer must have been allocated large enough for the region, perhaps using Rast_allocate_buf(). If there is an error writing the raster file, a warning message is printed and -1 is returned. Otherwise 1 is returned.

Note: The rows are written in sequential order. The first call writes row 0, the second writes row 1, etc. The following example assumes that the raster file is to be created:

int fd, row, nrows, ncols;
CELL *buf;
RASTER_MAP_TYPE data_type;
fd = Rast_open_old(inmap, ""); /*...data for the row */
data_type = Rast_get_map_type(fd);
Rast_put_row(fd, buf, data_type);

Closing Raster Files

All raster files are closed by one of the following routines, whether opened for reading or for writing.

The raster file opened on file descriptor is closed. Memory allocated for raster processing is freed. If open for writing, skeletal support files for the new raster file are created as well.

Note: If a module wants to explicitly write support files (e.g., a specific color table) for a raster file it creates, it must do so after the raster file is closed. Otherwise the close will overwrite the support files. See Raster Map Layer Support Routines for routines which write raster support files.

The raster file opened on file descriptor is closed. Memory allocated for raster processing is freed. If open for writing, the raster file is not created and the temporary file created when the raster file was opened is removed (see Creating and Opening New Raster Files).

This routine is useful when errors are detected and it is desired to not create the new raster file. While it is true that the raster file will not be created if the module exits without closing the file, the temporary file will not be removed at module exit. GRASS database management will eventually remove the temporary file, but the file can be quite large and will take up disk space until GRASS does remove it. Use this routine as a courtesy to the user.

Raster Map Layer Support Routines

GRASS map layers have a number of support files associated with them. These files are discussed in detail in Raster_Maps. The support files are the raster header, the category file, the color table, the history file, and the range file. Each support file has its own data structure and associated routines.

Raster Header File

The raster header file (cellhd) contains information describing the geographic extent of the raster map, the grid cell resolution, the format used to store the data in the raster file, see Cell_head structure for details. The routines described below use the Cell_head structure which is shown in detail in GIS Library Data Structures.

The raster header for the raster file in the specified mapset is read into the Cell_head structure.

Note: If the raster file is a reclassified, the raster header for the referenced raster file is read instead. See Rast_is_reclass() for distinguishing reclass files from regular raster files.

Note: It is not necessary to get the raster header for a map layer in order to read the raster file data. The routines which read raster file data automatically retrieve the raster header information and use it for resampling the raster file data into the active region. If it is necessary to read the raster file directly without resampling into the active region, then the raster header can be used to set the active region using G_set_window().

This function fills in missing parts of the input cell header (or region). It also makes projection-specific adjustments.

This routine writes the information from the Cell_head structure to the raster header file for the map layer in the current mapset.

Note: Programmers should have no reason to use this routine. It is used by Rast_close() to give new raster files correct header files, and by the r.support module to give users a means of creating or modifying raster headers.

This function determines if the raster file is a reclass file. Returns 1 if raster file is a reclass file, 0 if it is not, and -1 if there was a problem reading the raster header.

This function generates a child reclass maps list from the cell_misc/reclassed_to file which stores this list. The cell_misc/reclassed_to file is written by Rast_put_reclass(). Rast_is_reclassed_to() is used by g.rename, g.remove and r.reclass to prevent accidentally deleting the parent map of a reclassed raster map.

Raster split windows

The input window determines the number of rows and columns for Rast_get_row() etc, i.e. the valid range for the "row" parameter, and the number of elements in the "buf" array.

The output window determines the number of rows and columns for Rast_put_row() etc, i.e. the number of times it should be called ("put" operations don't have a "row" parameter) and the number of elements in the "buf" array.

For most modules, the input and output windows should be the same, but e.g. r.resamp.* will typically use split windows.

The rationale was that the old GRASS 6 mechanism was inefficient. It would change the window twice for each row of output (set "read" window, read rows, set "write" window, write row). Each change caused the column mapping to be recalculated for each input map. With split windows, there are separate windows for input maps and output maps, eliminating the need to change back and forth between input and output windows. The column mapping is calculated when a map is opened and never changes. Attempting to change the input window while any input maps are open generates an error, as does attempting to change the output window while any output maps are open.

Also, certain functions assume the existence of a single window, e.g. Rast_get_window(), Rast_window_{rows,cols}(). These functions still exist, and work so long as the windows aren't split (i.e. neither Rast_set_input_window() nor Rast_set_output_window() have been called, or the windows have since been joined by calling Rast_set_window()), but will generate an error if the windows are split. If the windows are split, the code must use e.g. Rast_get_input_window(), Rast_input_window_rows() etc to read the appropriate window.

The "split window" design eliminates the most common reason for changing the window while maps are open, although there may be cases it doesn't cover (e.g. reading multiple input maps at different resolutions won't work). If we need to support that, there are a number of possible solutions.

One is to store the current window in the fileinfo structure whenever the column mapping is created, check that it matches the current window whenever Rast_read_row() etc is called, and either re-generate it or generate an error if it differs. This avoids having to needlessly re- calculate the column mapping for all input maps on every set-window operation, but requires a window comparison (which probably needs some tolerance for comparing floating-point fields) for each get-row operation.

Another is to allow each map to have a separate window, set from the current window when the map is opened. The main drawback is that Rast_window_rows() etc become meaningless. We would need to add yet another "level" of window splitting, so you'd have: single window, read/write windows, per-map windows.

As usual, the main problem isn't implementation, but design and ensuring that existing code doesn't silently break (the current implementation should ensure that any breakage results in a fatal error).

Raster Category File

GRASS map layers have category labels associated with them. The category file is structured so that each category in the raster file can have a one-line description. The format of this file is described in Raster_Category_File_Format.

The routines described below manage the category file. Some of them use the Categories structure which is described in GIS Library Data Structures.

Reading and Writing the Raster Category File

The following routines read or write the category file itself:

The category file for raster file is read into the cats structure. If there is an error reading the category file, a diagnostic message is printed and -1 is returned. Otherwise, 0 is returned.

Writes the category file for the raster file in the current mapset from the cats structure. Returns 1 if successful. Otherwise, -1 is returned (no diagnostic is printed).

  • Rast_get_title()

If only the map layer title is needed, it is not necessary to read the entire category file into memory. This routine gets the title for raster file directly from the category file, and returns a pointer to the title. A legal pointer is always returned. If the map layer does not have a title, then a pointer to the empty string "" is returned.

  • Rast_put_title()

If it is only desired to change the title for a map layer, it is not necessary to read the entire category file into memory, change the title, and rewrite the category file. This routine changes the title for the raster file in the current mapset directly in the category file. It returns a pointer to the title.

Querying and Changing the Categories Structure

The following routines query or modify the information contained in the category structure:

This routine looks up category in the cats structure and returns a pointer to a string which is the label for the category. A legal pointer is always returned. If the category does not exist in cats then a pointer to the empty string "" is returned.

Warning: The pointer that is returned points to a hidden static buffer. Successive calls to Rast_get_cat() overwrite this buffer.

Map layers store a one-line title in the category structure as well. This routine returns a pointer to the title contained in the cats structure. A legal pointer is always returned. If the map layer does not have a title, then a pointer to the empty string "" is returned.

To construct a new category file, the structure must first be initialized. This routine initializes the cats structure, and copies the title into the structure.

For example:

struct Categories cats;
Rast_init_cats ((CELL) 0, "", &cats);

Copies the label into the cats structure for category.

Copies the title is copied into the cats structure.

Frees memory allocated by Rast_read_cats(), Rast_init_cats() and Rast_set_cat().

Raster Color Table

GRASS map layers have colors associated with them. The color tables are structured so that each category in the raster file has its own color. The format of this file is described in Raster_Color_Table_Format.

The routines that manipulate the raster color file use the Colors structure which is described in detail in GIS Library Data Structures.

Reading and Writing the Raster Color File

The following routines read, create, modify, and write color tables.

The color table for the raster file in the specified mapset is read into the colors structure. If the data layer has no color table, a default color table is generated and 0 is returned. If there is an error reading the color table, a diagnostic message is printed and -1 is returned. If the color table is read ok, 1 is returned.

Write map layer color table. The color table is written for the raster file in the specified mapset from the colors structure. If there is an error, -1 is returned. No diagnostic is printed. Otherwise, 1 is returned.

The colors structure must be created properly, i.e., Rast_init_colors() to initialize the structure and Rast_add_color_rule() to set the category colors.

Note: The calling sequence for this function deserves special attention. The mapset parameter seems to imply that it is possible to overwrite the color table for a raster file which is in another mapset. However, this is not what actually happens. It is very useful for users to create their own color tables for raster files in other mapsets, but without overwriting other users' color tables for the same raster file. If mapset is the current mapset, then the color file for name will be overwritten by the new color table. But if mapset is not the current mapset, then the color table is actually written in the current mapset under the colr2 element as: colr2/mapset/name.

Lookup Up Raster Colors

These routines translates raster values to their respective colors.

Extracts colors for an array of raster values. The colors from the raster array are stored in the red, green, and blue arrays.

Note: The red, green, and blue intensities will be in the range 0-255.

Get a category color. The red, green, and blue intensities for the color associated with category are extracted from the colors structure. The intensities will be in the range 0-255.

Creating and/or Modifying the Color Table

These routines allow the creation of customized color tables as well as the modification of existing tables.

Initialize colors structure for subsequent calls Rast_add_color_rule() and Rast_set_color().

This is the heart and soul of the color logic. It adds a color rule to the colors structure. The colors defined by the red, green, and blue values are assigned to the categories respectively. Colors for data values between two categories are not stored in the structure but are interpolated when queried by Rast_lookup_colors() and Rast_get_color(). The color components must be in the range 0-255.

For example, to create a linear grey scale for the range 200-1000:

struct Colors colr;
Rast_add_color_rule ((CELL) 200, 0,0,0,(CELL) 1000, 255,255,255) ;

The programmer is encouraged to review Raster_Color_Table_Format.

Note: The colors structure must have been initialized by Rast_init_colors(). See Predefined Color Tables for routines to build some predefined color tables.

  • Rast_set_color()

Set a category color. The red, green, and blue intensities for the color associated with category The intensities must be in the range 0-255. Values below zero are set as zero, values above 255 are set as 255.

Warning: Use of this routine is discouraged because it defeats the new color logic. It is provided only for backward compatibility. Overuse can create large color tables. Rast_add_color_rule() should be used whenever possible.

Note: The colors structure must have been initialized by Rast_init_color().

Get color range. Gets the minimum and maximum raster values that have colors associated with them.

Free dynamically allocated memory associated with the colors structure.

Note: This routine may be used after Rast_read_colors() as well as after Rast_init_colors().

Predefined Color Tables

The following routines generate entire color tables. The tables are loaded into a colors structure based on a range of category values from minimum to maximum value. The range of values for a raster map can be obtained, for example, using Rast_read_range().

Note: The color tables are generated without information about any particular raster file.

These color tables may be created for a raster map, but they may also be generated for loading graphics colors. These routines return -1 if minimum value is greater than maximum value, 1 otherwise.

Generates a color table for aspect data.

Generates a color table with 3 sections: red only, green only, and blue only, each increasing from none to full intensity. This table is good for continuous data, such as elevation.

Generates a color table with 3 sections: red only, green only, and blue only, each increasing from none to full intensity and back down to none. This table is good for continuous data like elevation.

Generates a grey scale color table. Each color is a level of grey, increasing from black to white.

Generates a "shifted" rainbow color table - yellow to green to cyan to blue to magenta to red. The color table is based on rainbow colors. (Normal rainbow colors are red, orange, yellow, green, blue, indigo, and violet.) This table is good for continuous data, such as elevation.

Generates random colors. Good as a first pass at a color table for nominal data.

Generates a color table that goes from red to yellow to green.

Generates a color table that goes from green to yellow to red.

Generates a histogram contrast-stretched grey scale color table that goes from the, histogram information.

Raster History File

The history file contains documentary information about the raster file: who created it, when it was created, what was the original data source, what information is contained in the raster file, etc.

The following routines manage this file. They use the History structure which is described in GIS Library Data Structures.

Note: This structure has existed relatively unmodified since the inception of GRASS. It is in need of overhaul. Programmers should be aware that future versions of GRASS may no longer support either the routines or the data structure which support the history file.

Reads raster history file. This routine reads the history file for the raster map into the History structure.

Writes raster history file. This routine writes the history file for the raster map in the current mapset from the History structure.

Note: The History structure should first be initialized using Rast_short_history().

This routine initializes History structure, recording the date, user, module name and the raster map.

Note: This routine only initializes the data structure. It does not write the history file.

Raster Range File

The following routines manage the raster range file. This file contains the minimum and maximum values found in the raster file. The format of this file is described in Raster_Range_File_Format.

The routines below use the Range data structure which is described in GIS Library Data Structures.

Reads raster range. This routine reads the range information for the raster map into the range structure. A diagnostic message is printed and -1 is returned if there is an error reading the range file. Otherwise, 0 is returned.

Write raster range file. This routine writes the range information for the raster map in the current mapset from the range structure. A diagnostic message is printed and -1 is returned if there is an error writing the range file. Otherwise, 0 is returned.

The range structure must be initialized and updated using the following routines:

Initializes the range structure for updates by Rast_update_range() and Rast_row_update_range().

Compares the category value with the minimum and maximum values in the range structure, modifying the range if the value extends the range.

This routine updates the range data just like Rast_update_range().

The range structure is queried using the following routine:

Get range minimum and maximum value.

Raster Histograms

The following routines provide a relatively efficient mechanism for computing and querying a histogram of raster data. They use the Cell_stats structure to hold the histogram information. The histogram is a count associated with each unique raster value representing the number of times each value was inserted into the structure.

These next two routines are used to manage the Cell_stats structure:

Initialize cell stats, This routine, which must be called first, initializes the Cell_stats structure.

Free cell stats. The memory associated with structure is freed. This routine may be called any time after calling Rast_init_cell_stats().

This next routine stores values in the histogram:

Adds data to cell stats. Once all values are stored, the structure may be queried either randomly (i.e., search for a specific raster value) or sequentially (retrieve all raster values, in ascending order, and their related count):

Random query of cell stats. This routine allows a random query of the Cell_stats structure. The routine returns 1 if cat was found in the structure, 0 otherwise.

Sequential retrieval is accomplished using these next 2 routines:

Reset/rewind cell stats. The structure is rewound (i.e., positioned at the first raster category) so that sorted sequential retrieval can begin.

Retrieve sorted cell stats. Retrieves the next cat, count combination from the structure. Returns 0 if there are no more items, non-zero if there are more.

For example:

struct Cell_stats s;
CELL cat;
long count;
/* updating s occurs here */
while(Rast_next_cell_stat(&cat, &count, &s))
fprintf(stdout, "%ld %ld\n", (long) cat, count);

GRASS 5 raster API

Needs to be merged into above sections.

The CELL Data Type

GRASS integer raster map data is defined to be of type CELL. This data type is defined in the "gis.h" header file. Programmers must declare all variables and buffers which will hold raster data or category codes as type CELL. Under GRASS the CELL data type is declared to be "int", but the programmer should not assume this. What should be assumed is that CELL is a signed integer type. It may be changed sometime to short or long. This implies that use of CELL data with routines which do not know about this data type (e.g., fprintf(stdout, ...), sscanf(), etc.) must use an intermediate variable of type long. To print a CELL value, it must be cast to long. For example:

CELL c; /* raster value to be printed */
/* some code to get a value for c */
fprintf(stdout, "%ld\n", (long) c); /* cast c to long to print */

To read a CELL value, for example from user typed input, it is necessary to read into a long variable, and then assign it to the CELL variable. For example (this example does not check for valid inputs, EOF, etc., which good code must do):

char userbuf[128];
CELL c;
long x;
fprintf (stdout, "Which category? "); /* prompt user */
gets(userbuf); /* get user response * /
sscanf (userbuf,"%ld", &x); /* scan category into long variable */
c = (CELL) x; /* assign long value to CELL value */

Of course, with GRASS library routines that are designed to handle the CELL type, this problem does not arise. It is only when CELL data must be used in routines which do not know about the CELL type, that the values must be cast to or from long.

Changes to gis.h

The "gis.h" contains 5 new items:

typedef float FCELL;
typedef double DCELL;
typedef int RASTER_MAP_TYPE;
#define CELL_TYPE 0
#define FCELL_TYPE 1
#define DCELL_TYPE 2

Also "gis.h" contains the definitions for new structures:

struct FPReclass;
struct FPRange;
struct Quant;

Some of the old structures such as

struct Categories
struct Cell_stats;
struct Range;
struct _Color_Rule_;
struct _Color_Info_;
struct Colors;

were modified, so it is very important to use functional interface to access and set elements of these structures instead of accessing elements of the structures directly. Because some former elements such as for example (struct Range range.pmin) do not exist anymore. It was made sure non of the former elements have different meaning, so that the programs which do access the old elements directly either do not compile or work exactly the same way as prior to change.

NULL-value functions

Set NULL value. For raster type CELL_TYPE use Rast_set_c_null_value(), for FCELL_TYPE Rast_set_f_null_value(), and for DCELL_TYPE Rast_set_d_null_value().

Insert NULL value. If raster type is CELL_TYPE calls Rast_insert_c_null_values(), FCELL_TYPE calls Rast_insert_f_null_values(), and DCELL_TYPE calls Rast_insert_d_null_values().

Check NULL value. If raster type is CELL_TYPE, calls Rast_is_c_null_value(), FCELL_TYPE calls Rast_is_f_null_value(), and DCELL_TYPE calls Rast_is_d_null_value().

It isn't good enough to test for a particular NaN bit pattern since the machine code may change this bit pattern to a different NaN. The test will be

if(fcell == 0.0) return 0;
if(fcell > 0.0) return 0;
if(fcell < 0.0) return 0;
return 1;

or (as suggested by Mark Line)

return(FCELL != fcell) ;

Allocates an array of char based on the number of columns in the current region.

Reads a row from NULL value bitmap file for the raster map open for read. If there is no bitmap file, then this routine simulates the read as follows: non-zero values in the raster map correspond to non-NULL; zero values correspond to NULL. When MASK exists, masked cells are set to null.

Floating-point and type-independent functions

Test for current maskreturns file descriptor number if MASK is in use and -1 if no MASK is in use.

Returns true(1) if raster map is a floating-point dataset; false(0) otherwise.

Returns the storage type for raster map:

  • CELL_TYPE for integer raster map
  • FCELL_TYPE for floating-point raster map
  • DCELL_TYPE for floating-point raster map (double precision)
  • Rast_open_new()

Creates new raster map in the current mapset. Calls G_open_map_new() for raster type CELL_TYPE, G_open_fp_map_new() for floating-point raster map. The use of this routine by applications is discouraged since its use would override user preferences (what precision to use).

This controls the storage type for floating-point raster maps. It affects subsequent calls to Rast_open_fp_new(). The type must be one of FCELL_TYPE (float) or DCELL_TYPE (double). The use of this routine by applications is discouraged since its use would override user preferences.

Creates a new floating-point raster map (in .tmp) and returns a file descriptor. The storage type (float or double) is determined by the last call to Rast_set_fp_type() or the default (float - unless the GRASS environmental variable GRASS_FP_DOUBLE is set).

Allocate an arrays of CELL, FCELL, or DCELL (depending on raster type) based on the number of columns in the current region.

  • Rast_incr_void_ptr()

Advances void pointer by n bytes. Returns new pointer value. Useful in raster row processing loops, substitutes

CELL *cell;
cell += n;

Now

rast = Rast_incr_void_ptr(rast, Rast_raster_size(data_type));

Where rast is void* and data_type is RASTER_MAP_TYPE can be used instead of rast++.) Very useful to generalize the row processing

  • loop (i.e. void * buf_ptr += Rast_raster_size(data_type))
    • Rast_raster_size()

If data_type is CELL_TYPE, returns sizeof(CELL), for FCELL_TYPE returns sizeof(FCELL), and for data_type is DCELL_TYPE, returns sizeof(DCELL).

Compares raster values.

Copies raster values.

  • Rast_set_value_c()

If Rast_is_c_null_value() is true, sets value to null value. Converts CELL value to raster data type value and stores the result. Used for assigning CELL values to raster cells of any type.

  • Rast_set_value_f()

If Rast_is_f_null_value() is true, sets value to null value. Converts FCELL val to raster data type and stores the result. Used for assigning FCELL values to raster cells of any type.

  • Rast_set_value_d()

If Rast_is_d_null_value() is true, sets value to null value. Converts DCELL val to raster data type and stores the result. Used for assigning DCELL values to raster cells of any type.

  • Rast_get_value_c()

Retrieves the value of raster type, converts it to CELL type and returns the result. If null value is stored, returns CELL null value. Used for retreiving CELL values from raster cells of any type.

Note: when data_type != CELL_TYPE, no quantization is used, only type conversion.

  • Rast_get_value_f()

Retrieves the value of raster type, converts it to FCELL type and returns the result. If null value is stored, returns FCELL null value. Used for retreiving FCELL values from raster cells of any type.

  • Rast_get_value_d()

Retrieves the value of raster type, converts it to DCELL type and returns the result. If null value is stored, returns DCELL null value. Used for retreiving DCELL values from raster cells of any type.

  • Rast_get_row ()

For CELL_TYPE raster type calls Rast_get_c_row(), FCELL_TYPE calls Rast_get_f_row(), and DCELL_TYPE Rast_get_d_row().

Same as Rast_get_f_row() except no masking occurs.

Read a row from the raster map performing type conversions as necessary based on the actual storage type of the map. Masking, resampling into the current region. NULL-values are always embedded (never converted to a value).

Same as Rast_get_f_row() except no masking occurs.

Same as Rast_get_f_row() except that the array double.

Same as Rast_get_d_row() except no masking occurs.

Reads a row of raster data and leaves the NULL values intact. (As opposed to the deprecated function Rast_get_row() which converts NULL values to zero.)

Note: When the raster map is old and null file doesn't exist, it is assumed that all 0-cells are no-data. When map is floating point, uses quant rules set explicitly by Rast_set_quant_rules or stored in map's quant file to convert floats to integers.

Same as Rast_get_c_row() except no masking occurs.

If raster type is CELL_TYPE, calls Rast_put_c_row(), if FCELL_TYPE, then calls Rast_put_f_row(), and for DCELL_TYPE, calls Rast_put_d_row().

Write the next row of the raster map performing type conversion to the actual storage type of the resultant map. Keep track of the range of floating-point values. Also writes the NULL-value bitmap from the NULL-values embedded in the data array.

Same as Rast_put_f_row() except that the array is double.

Writes a row of raster data and a row of the null-value bitmap, only treating NULL as NULL. (As opposed to the deprecated function Rast_put_row() which treats zero values also as NULL.)

  • Rast_zero_row()

Depending on raster type zeroes out Rast_window_cols() CELLs, FCELLs, or DCELLs stored in cell buffer.

Extracts a cell value from raster map at given position with nearest neighbor interpolation, bilinear interpolation or cubic interpolation.

Upgrades to Raster Functions (comparing to GRASS 4.x)

These routines will be modified (internally) to work with floating-point and NULL-values.

If the map is a new floating point, move the .tmp file into the fcell element, create an empty file in the cell directory; write the floating-point range file; write a default quantization file quantization file is set here to round fp numbers (this is a default for now). Create an empty category file, with max cat = max value (for backwards compatibility). Move the .tmp NULL-value bitmap file to the cell_misc directory.

Arrange for the NULL-value bitmap to be read as well as the raster map. If no NULL-value bitmap exists, arrange for the production of NULL-values based on zeros in the raster map.

If the map is floating-point, arrange for quantization to integer for Rast_get_c_row(), et. al., by reading the quantization rules for the map using Rast_read_quant().

If the programmer wants to read the floating point map using uing quant rules other than the ones stored in map's quant file, he/she should call Rast_set_quant_rules() after the call to Rast_open_old().

If the map is floating-point, quantize the floating-point values to integer using the quantization rules established for the map when the map was opened for reading (this quantization is read from cell_misc/name/f_quant file, but can be reset after opening raster map by Rast_set_quant_rules()). NULL values are converted to zeros. This routine is deprecated!!

Zero values are converted to NULLs. Write a row of the NULL value bit map. This routine is deprecated!!

NULL (no data) handling

The null file is stored in cell_misc/name/null file.

-2^31 (= 0x80000000 = -2147483648) is the null value for the CELL type, so you'll never see that value in a map.

The FP nulls are the all-ones bit patterns. These corresponds to NaN according to the IEEE-754 formats, although it isn't the "default" NaN pattern generated by most architectures (which is usually 7fc00000 or ffc00000 for float and 7ff8000000000000 or fff8000000000000 for double, i.e. an all-ones exponent, the top-bit of the mantissa set, and either sign).

So far as arithmetic is concerned, any value with an all-ones exponent and a non-zero mantissa is treated as NaN. Rast_set_d_null_value() and Rast_set_f_null_value() use the all-ones bit pattern. This is one of the many NaN values (anything with an all-ones exponent and a non-zero mantissa is NaN). As the topmost bit (i.e. the sign bit) is set, it is possible that some code would consider that to be "-NaN". E.g. code which writes a leading "-" based upon the sign bit before considering the other components would do so.

Rast_is_d_null_value() and Rast_is_f_null_value() treat any NaN as null (specifically, they test whether the value is unequal to itself).

At one time, these functions (or rather, their predecessors) checked explicitly for the all-ones pattern, but this was changed (in r33717 and r33752) to improve robustness. Apart from code explicitly setting a value to "null", NaNs can arise from calculations (0.0/0.0, sqrt(x) or log(x) for x<0, asin(x) or acos(x) for abs(x)>1, etc), and there's no guarantee as to exactly which NaN representation will result.

Presence or absence of null file:

For an integer map, any cells which were null will become zero, but any zeroes (cells which were previously either null or zero) will be treated as nulls (this is for compatibility with GRASS 4.x, which didn't have a null file, but typically used zero to indicate a null value).

For a floating-point map, any cells which were null will become zero (when writing FP data, a null has a zero written to the fcell/<map> file, and the corresponding bit is set in the null file).

Color Functions (new and upgraded)

Upgraded Colors structures

{
struct
{
int version; /* set by read_colors: -1=old,1=new */
DCELL shift;
int invert;
int is_float; /* defined on floating point raster data? */
int null_set; /* the colors for null are set? */
unsigned char null_red, null_grn, null_blu;
int undef_set; /* the colors for cells not in range are set? */
unsigned char undef_red, undef_grn, undef_blu;
struct _Color_Info_ fixed, modular;
DCELL cmin, cmax;
};

New functions to support colors for floating-point

If raster type is CELL_TYPE, calls Rast_lookup_colors(), if FCELL_TYPE, calls Rast_lookup_f_colors(), and for DCELL_TYPE Rast_lookup_d_colors().

Converts the floating-point values in the float data array to their r,g,b color components. Embedded NULL-values are handled properly as well.

Converts the floating-point values in the double data array to their r,g,b color components. Embedded NULL-values are handled properly as well.

If raster type is CELL_TYPE, calls Rast_add_c_color_rule(), if FCELL_TYPE, calls Rast_add_f_color_rule(), and for DCELL_TYPE calls Rast_add_d_color_rule().

Looks up the rgb colors for the value in the color table.

Sets a flag in the colors structure that indicates that these colors should only be looked up using floating-point raster data (not integer data).

New functions to support a color for the NULL-value

Sets the color (in colors) for the NULL-value to r,g,b.

Puts the red, green, and blue components of the color for the NULL-value into r,g,b.

New functions to support a default color

Sets the default color (in colors) to r,g,b. This is the color for values which do not have an explicit rule.

Puts the red, green, and blue components of the "default" color into r,g,b.

New functions to support treating a raster layer as a color image

Reads a row of raster data and converts it to red, green and blue components according to the colors.

This provides a convenient way to treat a raster layer as a color image without having to explicitly cater for each of CELL, FCELL and DCELL types.

Upgraded color functions

This routine reads the rules from the color file. If the input raster map is is a floating-point map (FCELL or DCELL) it calls Rast_mark_colors_as_fp().

The rules are written out using floating-point format, removing trailing zeros (possibly producing integers). The flag marking the colors as floating-point is not written.

  • Rast_get_colors_min_max()

If the color table is marked as "float", then return the minimum as -(255&#94;3 * 128) and the maximum as (255&#94;3 * 128). This is to simulate a very large range so that GRASS doesn't attempt to use colormode float to allow interactive toggling of colors.

Modified to return a color for NULL-values.

Modified to return a color for the NULL-value.

Changes to the Colors structure

Modifications to the Colors structure to support colors for floating-point data and the NULL-value consist of

  • the Color_Rule struct was changed to have DCELL value (instead of CELL cat) to have the range be floating-point values instead of integer cats.
  • a color for NULL was added
  • the special color for zero was eliminated
  • a default color for values which have no assigned color was added
  • a flag was added to the Colors structure to indicate if either the map itself is floating-point (If the map is integer and the floating point functions are used to lookup colors, the values are checked to see if they are integer, and if they are, the integer mechanism is used)
  • fp_lookup - a lookup table for floating point numbers is added. It orders the end points of fp intervals into array with a pointer to a color rule for each interval, and the binary search is then used when looking up colors instead of linearly searching through all color rules.

Changes to the \c colr file

  • The rules are written out using floating-point format, removing trailing zeros (possibly producing integers) . For example, to ramp from red to green for the range [1.3,5.0]:
                1.3:255:0:0  5:0:255:0
    
  • The NULL-value color is written as:
                nv:red:grn:blu
    
  • The default color (for values that don't have an explicit rule) is written as:
                *:red:grn:blu
    

Range functions (new and upgraded)

Modified range functions

Old range file (those with 4 numbers) should treat zeros in this file as NULL-values. New range files (those with just 2 numbers) should treat these numbers as real data (zeros are real data in this case).

An empty range file indicates that the min, max are undefined. This is a valid case, and the result should be an initialized range struct with no defined min/max.

If the range file is missing and the map is a floating-point map, this function will create a default range by calling Rast_construct_default_range().

Must set a flag in the range structure that indicates that no min/max have been defined - probably a "first" boolean flag.

NULL-values must be detected and ignored.

If the range structure has no defined min/max (first!=0) there will not be a valid range. In this case the min and max returned must be the NULL-value.

This routine only writes 2 numbers (min,max) to the range file, instead of the 4 (pmin,pmax,nmin,nmax) previously written. If there is no defined min,max, an empty file is written.

New range functions

Sets the integer range to [1,255].

If raster type is CELL_TYPE, calls Rast_read_range(), otherwise calls Rast_read_fp_range().

Read the floating point range file f_range. This file is written in binary using XDR format. If there is no defined min/max in FPRange structure, an empty f_range file is created.

An empty range file indicates that the min, max are undefined. This is a valid case, and the result should be an initialized range struct with no defined min/max.

If the range file is missing and the map is a floating-point map, this function will create a default range by calling Rast_construct_default_range().

If raster type is CELL_TYPE, calls Rast_init_range(), otherwise calls Rast_init_fp_range().

Must set a flag in the range structure that indicates that no min/max have been defined - probably a "first" boolean flag.

  • Rast_update_f_range()
  • Rast_update_d_range()

Updates the floating-point range from the values in NULL-values must be detected and ignored.

Extract the min/max from the FPRange structure. If the range structure has no defined min/max (first!=0) there will not be a valid range. In this case the min and max returned must be the NULL-value.

Write the floating point range file f_range. This file is written in binary using XDR format. If there is no defined min/max in r, an empty f_range file is created.

New and Upgraded Cell_stats functions

Modified Cell_stats functions to handle NULL-values:

Set the count for NULL-values to zero.

Look for NULLs and update the NULL-value count.

Do not return a record for the NULL-value.

Allow finding the count for the NULL-value.

Get a number of null values from stats structure. Note: when reporting values which appear in a map using Rast_next_cell_stats() , to get stats for null, call Rast_get_stats_for_null_value() first, since Rast_next_cell_stats() does not report stats for null.

New Quantization Functions

New functions to support quantization of floating-point to integer:

Writes the f_quant file for the raster map. If mapset==Rast_mapset() i.e. the map is in current mapset, then the original quant file in cell_misc/map/f_quant is written. Othewise is written into quant2/mapset/name (much like colr2 element). This results in map&#64;mapset being read using quant rules stored Rast_mapset(). See Rast_read_quant() for detailes.

Sets quant translation rules for raster map opened for reading. After calling this function, Rast_get_c_row() and Rast_get_row() will use defined rules (instead of using rules defined in map's quant file) to convert floats to ints.

Reads quantization rules for raster map and stores them in the quantization structure. If the map is in another mapset, first checks for quant2 table for this map in current mapset.

Initializes the Quant struct.

Frees any memory allocated in Quant structure and re-initializes the structure by calling Rast_quant_init().

Sets the quant rules to perform simple truncation on floats.

Sets the quant rules to perform simple rounding on floats.

  • Rast_quant_organize_fp_lookup()

Organizes fp_lookup table for faster (logarithmic) lookup time Rast_quant_organize_fp_lookup() creates a list of min and max for each quant rule, sorts this list, and stores the pointer to quant rule that should be used in between any 2 numbers in this list Also it stores extreme points for 2 infinite rules, if exist After the call to Rast_quant_organize_fp_lookup() instead of linearly searching through list of rules to find a rule to apply, quant lookup will perform a binary search to find an interval containing floating point value, and then use the rule associated with this interval. When the value doesn't fall within any interval, check for the infinite rules.

Add the rule that the floating-point range produces an integer in the range by linear interpolation. Rules that are added later have higher precedence when searching. If any of of the values is the NULL-value, this rule is not added and 0 is returned. Otherwise return 1. if the fp_lookup is organized, destroy it.

  • Rast_quant_set_positive_infinite_rule()
  • Rast_quant_get_positive_infinite_rule()

This rule has lower precedence than rules added with Rast_quant_add_rule().

  • Rast_quant_set_negative_infinite_rule()

This rule has lower precedence than rules added with Rast_quant_add_rule().

Extracts the minimum and maximum floating-point and integer values from all the rules.

  • Rast_quant_nrules()

Returns the number of rules, excluding the negative and positive "infinite" rules.

  • Rast_quant_get_rule()

The order of the rules returned by increasing n is the order in which the rules are applied when quantizing a value - the first rule applicable is used.

These next two functions are convenience functions to allow applications to easily create quantization rules other than the defaults:

Writes the f_quant file for the raster map.

Writes the f_quant file for the raster map with one rule.

Categories Labeling Functions (new and upgraded)

Upgraded Categories structure

All the new programs which are using Categories structure directly have to be modified to use API functions to update and retrieve info from Categories structure. Both new and old API function can be used, since old functions still have exact same functionality (even though internally they are implemented very differently). New function names end with raster_cats(); old function names end with _cats().

We made sure that all old fields in Categories structure are either missing in new Categories structure or have exactly the same meaning. We did it so that the modules using Categories structure directly either do not compile with new gis library or work exactly the same as bnefore. A programmer might want to read the data in a floating point map in a way that each cell value stores index of it's category label and data range. The way to do it is to call Rast_set_quant_rules() after opening the map.

This is helpful when trying to collect statistics (how many cells of each category are in the map. (although there is another new mechanism to collect such stats - see Rast_mark_cats()) . Another reason to get a category index instead of fp values is that this index will be the FID into GRASS-DBMS link. Also he can use Rast_get_ith_cat() to get the category information for each cell using this index.

Here is the new Categories structure defined in gis.h:

struct Categories
{
CELL ncats ; /* total number of categories */
CELL num ; /* the highest cell values. Only exists
for backwards compatibility = (CELL)
max_fp_values in quant rules */
char *title ; /* name of data layer */
char *fmt ; /* printf-like format to generate labels */
float m1 ; /* Multiplication coefficient 1 */
float a1 ; /* Addition coefficient 1 */
float m2 ; /* Multiplication coefficient 2 */
float a2 ; /* Addition coefficient 2 */
struct Quant q ; /* rules mapping cell values to index in
list of labels */
char **labels ; /* array of labels of size num */
int * marks ; /* was the value with this label was used? */
int nalloc;
} ;

Changes to the \c cats file

The format of explicit label entries is the same for integer maps.

    cat:description

In addition label entries of new format is supported for floating point maps.

    val:descr (where val is a floating point number) 

or

    val1:val2:descr (where val1, val2 is a floating point range) 

Internally the labels are stored for fp ranges of data. However when the cats file is written, all the decimal zeros are stripped so that integer values appear as integers in the file. Also if values are the same, only 1 value is written (i.e. first format).

This way even though the old cats files will be processed differently internally, the user or application programmer will not notice this difference as long as the proper api is used and the elements of Categories structure are not accessed directly without API calls.

Range functions (new and upgraded)

New Functions to read/write access and modify Categories structure

Is the same as existing Rast_read_cats().

Returns pointer to a string describing category.

Adds the label for range in category structure.

Returns the number of labels. DO NOT use Rast_number_of_cats() (it returns max cat number).

Returns i-th description and i-th data range from the list of category descriptions with corresponding data ranges.

Returns pointer to a string with title.

Sets marks for all categories to 0. This initializes Categories structure for subsequent calls to Rast_mark_cats () for each row of data, where non-zero mark for i-th label means that some of the cells in rast_row are labeled with i-th label and fall into i-th data range.

These marks help determine from the Categories structure which labels were used and which weren't.

Finds the next label and corresponding data range in the list of marked categories. The category (label + data range) is marked by Rast_mark_cats(). End points of the data range are converted to raster type and returned. The number of times value from i-th cat. data range appeared so far is returned in stats. See Rast_unmark_cats(), Rast_rewind_cats() and Rast_mark_cats ().

Looks up the category label for each raster value in the raster row (row of raster cell value) and updates the marks for labels found.

Note: non-zero mark for i-th label stores the number of of raster cells read so far which are labeled with i-th label and fall into i-th data range.

After call to this function Rast_get_next_marked_cat() returns the first marked cat label.

Same as existing Rast_init_cats() only ncats argument is missign. ncats has no meaning in new Categories structure and only stores (int) largets data value for backwards compatibility.

Same as existing Rast_set_cats_fmt().

Same as existing Rast_set_cats_title().

Same as existing Rast_write_cats().

Same as existing Rast_free_cats().

Library Functions that are Deprecated

These functions are deprecated, since they imply that the application that uses them has not been upgraded to handle NULL-values and should be eliminated from GRASS code.

To be replaced by Rast_get_c_row().

To be replaced by Rast_get_c_row_nomask().

To be replaced by Rast_put_c_row().

These functions are deprecated, since they can not be upgraded to support NULL-values, and should be eliminated from GRASS code.

  • Rast_open_map_new_random()
  • Rast_put_row_random()

Also, no support for random writing of floating-point rasters will be provided.

Guidelines for upgrading GRASS 4.x Modules

  • Modules that process raster maps as continuous data should read raster maps as floating-point. Modules that process raster maps as nominal data should read raster maps as integer.

Exception: Modules that process raster colors or the modules which report on raster categories labels should either always read the maps as floating-point, or read the maps as integer if the map is integer and floating-point if the map is floating-point.

  • The quantization of floating-point to integer should NOT change the color table. The color lookup should have its own separate quantization.
  • The quantization of floating-point to integer should NOT change the Categories table. The Categories structure should have its own separate quantization.
  • Modules that read or write floating-point raster maps should use double (DCELL) arrays instead of float (FCELL) arrays.
  • Modues should process NULL values in a well defined (consistent) manner. Modules that processed zero as the pseudo NULL-value should be changed to use the true NULL-value for this and process zero as normal value.
  • Modules should process non-NULL values as normal numbers and not treat any particular numbers (e.g. zero) as special.

Important hints for upgrades to raster modules

In general modules that use Rast_get_row(). Should use Rast_get_c_row() instead.

Modules that use Rast_put_row(). Should use Rast_put_c_row() instead.

Loading the Raster Library

The library is loaded by specifying

$(RASTERLIB)

in the Makefile.

List of functions

TODO: Reorder functions

Memory allocation

TODO: used in r.null and r.support, why Rast__ and not Rast_?

Raster Mask

Raster statistics

Raster category management

Raster file management

Color functions

Raster floating point management

Raster map reading

Raster map histogram

Raster map history

Raster map interpolation

Raster map quantization

yet unsorted functions