/*! \page rasterlib GRASS Raster Library
by GRASS Development Team (https://grass.osgeo.org)
\tableofcontents
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
\section gisrastintro 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).
\section Finding_Raster_Files_in_the_Database 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:
- G_find_raster()
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 "@mapset".
For example, to find a raster map in all mapsets listed the mapset search
path:
\code
char name[GNAME_MAX];
char *mapset;
if ((mapset = G_find_raster(name,"")) == NULL)
/* not found */
\endcode
To check that the raster map exists in the current mapset:
\code
char name[GNAME_MAX];
if (G_find_raster(name, G_mapset()) == NULL)
/* not found */
\endcode
\section Opening_an_Existing_Raster_File Opening an Existing Raster File
The following routine opens the raster map file for reading.
- Rast_open_old()
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.
\section Creating_and_Opening_New_Raster_Files 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).
- Rast_open_new()
- Rast_open_c_new()
- Rast_open_fp_new()
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.
- Rast_open_new_uncompressed()
- Rast_open_c_new_uncompressed()
- Rast_open_fp_new_uncompressed()
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.
\section Allocating_Raster_I_O_Buffers 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.
- Rast_allocate_buf()
- Rast_allocate_c_buf()
- Rast_allocate_f_buf()
- Rast_allocate_d_buf()
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).
\code
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);
\endcode
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().
\section Reading_Raster_Files 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.
- Rast_get_row()
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.
- Rast_get_row_nomask()
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 \ref
Mask for more information about the mask.
\subsection Writing_Raster_Files Writing Raster Files
- Rast_put_row()
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:
\code
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);
\endcode
\subsection Closing_Raster_Files Closing Raster Files
All raster files are closed by one of the following routines, whether
opened for reading or for writing.
- Rast_close()
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 \ref Raster_Map_Layer_Support_Routines
for routines which write raster support files.
- Rast_unopen()
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 \ref 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.
\section Raster_Map_Layer_Support_Routines Raster Map Layer Support Routines
GRASS map layers have a number of support files associated with
them. These files are discussed in detail in \ref 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.
\section Raster_Header_File 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 \ref
GIS_Library_Data_Structures.
- Rast_get_cellhd()
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().
- G_adjust_Cell_head()
This function fills in missing parts of the input cell header (or
region). It also makes projection-specific adjustments.
- Rast_put_cellhd()
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.
- Rast_is_reclass()
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.
- Rast_is_reclassed_to()
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.
\section Raster_split_windows 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).
\section Raster_Category_File 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 \ref Raster_Category_File_Format.
The routines described below manage the category file. Some of them
use the Categories structure which is described in \ref
GIS_Library_Data_Structures.
\section Reading_and_Writing_the_Raster_Category_File Reading and Writing the Raster Category File
The following routines read or write the category file itself:
- Rast_read_cats()
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.
- Rast_write_cats()
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.
\section Querying_and_Changing_the_Categories_Structure Querying and Changing the Categories Structure
The following routines query or modify the information contained in
the category structure:
- Rast_get_cat()
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.
- Rast_get_cats_title()
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.
- Rast_init_cats()
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:
\code
struct Categories cats;
Rast_init_cats ((CELL) 0, "", &cats);
\endcode
- Rast_set_cat()
Copies the label into the cats structure for category.
- Rast_set_cats_title()
Copies the title is copied into the cats structure.
- Rast_free_cats()
Frees memory allocated by Rast_read_cats(), Rast_init_cats() and
Rast_set_cat().
\section Raster_Color_Table 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
\ref Raster_Color_Table_Format.
The routines that manipulate the raster color file use the
Colors structure which is described in detail in
\ref GIS_Library_Data_Structures.
\section Reading_and_Writing_the_Raster_Color_File Reading and Writing the Raster Color File
The following routines read, create, modify, and write color tables.
- Rast_read_colors()
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.
- Rast_write_colors()
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.
\section Lookup_Up_Raster_Colors Lookup Up Raster Colors
These routines translates raster values to their respective colors.
- Rast_lookup_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.
- Rast_get_color()
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.
\section Creating_and_or_Modifying_the_Color_Table Creating and/or Modifying the Color Table
These routines allow the creation of customized color tables as well as the
modification of existing tables.
- Rast_init_colors()
Initialize colors structure for subsequent calls
Rast_add_color_rule() and Rast_set_color().
- Rast_add_color_rule()
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:
\code
struct Colors colr;
Rast_init_colors (&colr) ;
Rast_add_color_rule ((CELL) 200, 0,0,0,(CELL) 1000, 255,255,255) ;
\endcode
The programmer is encouraged to review \ref Raster_Color_Table_Format.
Note: The colors structure must have been initialized
by Rast_init_colors(). See \ref 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().
- Rast_get_color_range()
- Rast_get_d_color_range()
Get color range. Gets the minimum and maximum raster values that have
colors associated with them.
- Rast_free_colors()
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().
\section Predefined_Color_Tables 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.
- Rast_make_aspect_colors()
Generates a color table for aspect data.
- Rast_make_ramp_colors()
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.
- Rast_make_wave_colors()
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.
- Rast_make_grey_scale_colors()
Generates a grey scale color table. Each color is a level of grey,
increasing from black to white.
- Rast_make_rainbow_colors()
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.
- Rast_make_random_colors()
Generates random colors. Good as a first pass at a color table for
nominal data.
- Rast_make_ryg_colors()
Generates a color table that goes from red to yellow to green.
- Rast_make_gyr_colors()
Generates a color table that goes from green to yellow to red.
- Rast_make_histogram_eq_colors()
Generates a histogram contrast-stretched grey scale color table that goes
from the, histogram information.
\section Raster_History_File 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 \ref 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.
- Rast_read_history()
Reads raster history file. This routine reads the history file for the
raster map into the History structure.
- Rast_write_history()
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().
- 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.
\section Raster_Range_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 \ref Raster_Range_File_Format.
The routines below use the Range data structure which is
described in \ref GIS_Library_Data_Structures.
- Rast_read_range()
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.
- Rast_write_range()
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:
- Rast_init_range()
Initializes the range structure for updates by Rast_update_range()
and Rast_row_update_range().
- Rast_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.
- Rast_row_update_range()
This routine updates the range data just like Rast_update_range().
The range structure is queried using the following routine:
- Rast_get_range_min_max()
Get range minimum and maximum value.
\section Raster_Histograms 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:
- Rast_init_cell_stats()
Initialize cell stats, This routine, which must be called first,
initializes the Cell_stats structure.
- Rast_free_cell_stats()
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:
- Rast_update_cell_stats()
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):
- Rast_find_cell_stat()
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:
- Rast_rewind_cell_stats()
Reset/rewind cell stats. The structure is rewound (i.e., positioned at
the first raster category) so that sorted sequential retrieval can
begin.
- Rast_next_cell_stat()
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:
\code
struct Cell_stats s;
CELL cat;
long count;
/* updating s occurs here */
Rast_rewind_cell_stats(&s);
while(Rast_next_cell_stat(&cat, &count, &s))
fprintf(stdout, "%ld %ld\n", (long) cat, count);
\endcode
\section GRASS_5_raster_API GRASS 5 raster API
Needs to be merged into above sections.
\subsection The_CELL_Data_Type 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:
\code
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 */
\endcode
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):
\code
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 */
\endcode
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.
\subsection GRASS_5_raster_API_gish Changes to gis.h
The "gis.h" contains 5 new items:
\code
typedef float FCELL;
typedef double DCELL;
typedef int RASTER_MAP_TYPE;
#define CELL_TYPE 0
#define FCELL_TYPE 1
#define DCELL_TYPE 2
\endcode
Also "gis.h" contains the definitions for new structures:
\code
struct FPReclass;
struct FPRange;
struct Quant;
\endcode
Some of the old structures such as
\code
struct Categories
struct Cell_stats;
struct Range;
struct _Color_Rule_;
struct _Color_Info_;
struct Colors;
\endcode
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.
\subsection NULL_value_functions NULL-value functions
- Rast_set_null_value()
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().
- Rast_insert_null_values()
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().
- Rast_is_null_value()
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
\code
if(fcell == 0.0) return 0;
if(fcell > 0.0) return 0;
if(fcell < 0.0) return 0;
return 1;
\endcode
or (as suggested by Mark Line)
\code
return(FCELL != fcell) ;
\endcode
- Rast_allocate_null_buf()
Allocates an array of char based on the number of columns in the
current region.
- Rast_get_null_value_row()
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.
\subsection Floating_point_and_type_independent_functions Floating-point and type-independent functions
- Rast_maskfd()
Test for current maskreturns file descriptor number if MASK is in use
and -1 if no MASK is in use.
- Rast_map_is_fp()
Returns true(1) if raster map is a floating-point dataset; false(0)
otherwise.
- Rast_map_type()
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).
- Rast_set_fp_type()
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.
- Rast_open_fp_new()
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).
- Rast_allocate_buf()
- Rast_allocate_c_buf()
- Rast_allocate_f_buf()
- Rast_allocate_d_buf()
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
\code
CELL *cell;
cell += n;
\endcode
Now
\code
rast = Rast_incr_void_ptr(rast, Rast_raster_size(data_type));
\endcode
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).
- Rast_raster_cmp()
Compares raster values.
- Rast_raster_cpy()
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().
- Rast_get_row_nomask().
Same as Rast_get_f_row() except no masking occurs.
- Rast_get_f_row()
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).
- Rast_get_f_row_nomask()
Same as Rast_get_f_row() except no masking occurs.
- Rast_get_d_row()
Same as Rast_get_f_row() except that the array double.
- Rast_get_d_row_nomask()
Same as Rast_get_d_row() except no masking occurs.
- Rast_get_c_row()
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.
- Rast_get_c_row_nomask()
Same as Rast_get_c_row() except no masking occurs.
- Rast_put_row()
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().
- Rast_put_f_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.
- Rast_put_d_row()
Same as Rast_put_f_row() except that the array is double.
- Rast_put_c_row()
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.
- Rast_get_sample()
Extracts a cell value from raster map at given position with nearest neighbor interpolation, bilinear interpolation or cubic interpolation.
\section Upgrades_to_Raster_Functions Upgrades to Raster Functions (comparing to GRASS 4.x)
These routines will be modified (internally) to work with
floating-point and NULL-values.
- Rast_close()
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.
- Rast_open_old()
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().
- Rast_get_row()
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!!
- Rast_put_row()
Zero values are converted to NULLs. Write a row of the NULL value bit
map. This routine is deprecated!!
\section Null_no_data 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/\