I wrote a little module that does this sort of calculation a while ago.
If I remember right, it worked reasonably well and also does
normalization either by neighborhood or map maximum. I have attached
the code and Makefile.
Have fun with it,
Ben
chris carleton wrote:
I've encountered another challenge with GRASS that I need some expert
advice with. I want to run a neighborhood analysis on a few DEMs, but
the options in r.neighbors don't quite seem to suit my needs. I want an
analysis that calculates the average difference between the central cell
and its neighbors (using that average as the new cell cat value). The
idea is that I'll end up with an approximation of terrain 'ruggedness'
by looking at average differences in elevation within a given
neighborhood. Any suggestions would be welcome.
Chris
------------------------------------------------------------------------
------------------------------------------------------------------------
_______________________________________________
grass-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-dev
--
Benjamin Ducke
Senior Applications Support and Development Officer
Oxford Archaeological Unit Limited
Janus House
Osney Mead
OX2 0ES
Oxford, U.K.
Tel: +44 (0)1865 263 800 (switchboard)
Tel: +44 (0)1865 980 758 (direct)
Fax :+44 (0)1865 793 496
[EMAIL PROTECTED]
------
Files attached to this email may be in ISO 26300 format (OASIS Open Document
Format). If you have difficulty opening them, please visit http://iso26300.info
for more information.
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <math.h>
#include <grass/gis.h>
#define PROGVERSION 1.0
int
main (int argc, char *argv[])
{
struct GModule *module;
struct
{
struct Option *input;
struct Option *output;
struct Option *radius;
} parm;
struct
{
struct Flag *localnorm;
struct Flag *globalnorm;
struct Flag *rect;
struct Flag *absolute;
struct Flag *quiet;
} flag;
char *sysstr;
char *mapset;
int i, j, radius, n, x, y;
int a, b, dist;
int firstrun;
int nullvalue;
double centerval;
double sum, prominence;
double min, max;
double from, to;
int nrows, ncols;
int fd, fd_out;
struct Cell_head region;
struct Colors colors;
DCELL *diskrow = NULL;
DCELL *outrow = NULL;
sysstr = G_calloc (512, sizeof (char));
module = G_define_module ();
module->description = "Calculates Llobera's prominence index";
/* setup some basic GIS stuff */
G_gisinit (argv[0]);
/* do not pause after a warning message was displayed */
G_sleep_on_error (0);
parm.input = G_define_standard_option(G_OPT_R_INPUT);
parm.input->key = "input";
parm.input->type = TYPE_STRING;
parm.input->required = YES;
parm.input->gisprompt = "old,fcell";
parm.input->description = "FP raster map for which to calculate index";
parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
parm.output->key = "output";
parm.output->type = TYPE_STRING;
parm.output->required = YES;
parm.output->gisprompt = "fcell";
parm.output->description = "FP raster to store output";
parm.radius = G_define_option ();
parm.radius->key = "radius";
parm.radius->type = TYPE_INTEGER;
parm.radius->required = YES;
parm.radius->description = "Radius of neighbourhood in map cells";
flag.absolute = G_define_flag ();
flag.absolute->key = 'a';
flag.absolute->description = "Calculate absolute differences";
flag.localnorm = G_define_flag ();
flag.localnorm->key = 'l';
flag.localnorm->description = "Local data normalisation";
flag.globalnorm = G_define_flag ();
flag.globalnorm->key = 'g';
flag.globalnorm->description = "Global data normalisation";
flag.rect = G_define_flag ();
flag.rect->key = 'r';
flag.rect->description = "Use quadratic neighbourhood";
flag.quiet = G_define_flag ();
flag.quiet->key = 'q';
flag.quiet->description = "Disable on-screen progress display";
/* parse command line */
if (G_parser (argc, argv))
{
exit (-1);
}
G_get_window (®ion);
nrows = G_window_rows ();
ncols = G_window_cols ();
/* check parameters for validity */
radius = atoi (parm.radius->answer);
if (( radius < 1 ) || ( radius > nrows/2 ) || ( radius > ncols/2 )) {
G_fatal_error ("Radius must be > 1 and smaller than half the map's extent in both directions.");
}
if ((flag.localnorm->answer) && (flag.globalnorm->answer)) {
G_fatal_error ("Please choose either local OR global normalisation.");
}
mapset = G_calloc (512, sizeof (char));
mapset = G_find_cell (parm.input->answer, "");
if ( mapset == NULL) {
G_fatal_error ("Input map does not exist in the current location.");
}
if (!G_raster_map_is_fp (parm.input->answer, mapset)) {
G_fatal_error ("Input map is not a floating point map.\nOnly integer maps are allowed as input maps.");
}
fd = G_open_cell_old (parm.input->answer, mapset);
if (fd < 0) {
G_fatal_error ("Could not open input map for reading!\n");
}
if (!G_legal_filename (parm.output->answer)) {
G_fatal_error ("%s is not a legal filename for an output map.\n",
parm.output->answer);
}
fd_out = G_open_raster_new (parm.output->answer, DCELL_TYPE);
if ( fd_out < 0 ) {
G_fatal_error ("Could not open output map for writing!\n");
}
diskrow = G_allocate_d_raster_buf ();
outrow = G_allocate_d_raster_buf ();
/* initialize statistics */
n = 0;
firstrun = 1;
min = 0;
max = 0;
sum = 0.0;
if (!flag.quiet->answer) {
fprintf (stdout,"Calculating prominence index:\n");
fflush (stdout);
}
/* main loop: read raster row from disk, add to values in neighbourhood */
i = 0;
j = 0;
nullvalue = 0;
for (y = 0; y < nrows; y ++) {
for (x = 0; x < ncols; x++) {
G_get_d_raster_row (fd, diskrow, y);
centerval = diskrow [x];
if (G_is_d_null_value ( &diskrow[x] )) {
nullvalue = 1;
} else {
nullvalue = 0;
for (j = y-radius; j <= y+radius; j++) {
if ((j>0) && (j<nrows)) {
G_get_d_raster_row (fd, diskrow, j);
nullvalue = 0;
for (i = x-radius; i <= x+radius; i++) {
if ( ( (i!=x) || (j!=y) ) && (i>0) && (i<ncols)) {
if (!G_is_d_null_value (&diskrow[i])) {
if (flag.rect->answer) {
/* use rectangular neighbourhood */
n = n + 1;
if (flag.absolute->answer) {
sum = sum + fabs ( centerval - diskrow [i]);
} else {
sum = sum + ( centerval - diskrow [i]);
}
} else {
/* use circular neighbourhood */
a = abs ( x - i );
b = abs ( y - j );
dist = (int) sqrt ( (double) (a*a) + (b*b) );
if (dist <= radius) {
n = n + 1;
if (flag.absolute->answer) {
sum = sum + fabs ( centerval - diskrow [i]);
} else {
sum = sum + ( centerval - diskrow [i]);
}
}
}
}
}
}
}
}
}
/* write one cell to output map */
prominence = sum/n;
if (!nullvalue) {
outrow [x] = prominence;
} else {
G_set_d_null_value ( &outrow[x],1 );
}
if (firstrun) {
min = prominence;
max = prominence;
firstrun = 0;
} else {
if ( prominence < min ) min = prominence;
if ( prominence > max ) max = prominence;
}
n = 0;
sum = 0.0;
}
/* write one row to output map on disk */
G_put_raster_row (fd_out, outrow, DCELL_TYPE);
if (!flag.quiet->answer) {
G_percent (y, nrows-1, 1);
fflush (stdout);
}
}
G_close_cell (fd);
G_close_cell (fd_out);
if (flag.localnorm->answer) {
if (!flag.quiet->answer) {
fprintf (stdout,"\nNormalising data by neighbourhood maximum:\n");
fflush (stdout);
}
/* re-open output map for reading */
fd = G_open_cell_old (parm.output->answer, G_mapset());
if (fd < 0) {
G_fatal_error ("Could not open output map for normalisation pass!\n");
}
G_strcpy (sysstr, parm.output->answer);
G_strcat (sysstr, ".localnorm");
fd_out = G_open_raster_new (sysstr, DCELL_TYPE);
if ( fd_out < 0 ) {
G_fatal_error ("Could not open temporary map for normalisation pass (write access)!\n");
}
i = 0;
j = 0;
for (y = 0; y < nrows; y ++) {
for (x = 0; x < ncols; x++) {
G_get_d_raster_row (fd, diskrow, y);
centerval = diskrow [x];
min = diskrow [x];
max = diskrow [x];
if (G_is_d_null_value ( &diskrow[x] )) {
nullvalue = 1;
} else {
nullvalue = 0;
for (j = y-radius; j <= y+radius; j++) {
if ((j>0) && (j<nrows)) {
G_get_d_raster_row (fd, diskrow, j);
nullvalue = 0;
for (i = x-radius; i <= x+radius; i++) {
if ( ( (i!=x) || (j!=y) ) && (i>0) && (i<ncols)) {
if (!G_is_d_null_value (&diskrow[i])) {
if (flag.rect->answer) {
/* use rectangular neighbourhood */
if ( diskrow [i] > max ) {
max = diskrow [i];
}
if ( diskrow [i] < min ) {
min = diskrow [i];
}
} else {
/* use circular neighbourhood */
a = abs ( x - i );
b = abs ( y - j );
dist = (int) sqrt ( (double) (a*a) + (b*b) );
if (dist <= radius) {
if ( diskrow[i] > max ) {
max = diskrow [i];
}
if ( diskrow[i] < min ) {
min = diskrow [i];
}
}
}
}
}
}
}
}
}
/* write one cell to output map */
if (!nullvalue) {
if (flag.absolute->answer) {
outrow [x] = (centerval - min) / ( max - min );
} else {
if ( min < 0 ) {
if ( centerval < 0 ) {
outrow [x] = centerval / fabs (min);
}
if ( centerval > 0 ) {
outrow [x] = centerval / max;
}
if ( centerval == 0 ) {
outrow [x] = 0;
}
} else {
outrow [x] = (centerval - min) / ( max - min );
}
}
} else {
G_set_d_null_value ( &outrow[x],1 );
}
}
/* write one row to output map on disk */
G_put_raster_row (fd_out, outrow, DCELL_TYPE);
if (!flag.quiet->answer) {
G_percent (y, nrows-1, 1);
fflush (stdout);
}
}
G_close_cell (fd);
G_close_cell (fd_out);
}
if (flag.globalnorm->answer) {
/* normalise evidence over total map */
if (!flag.quiet->answer) {
fprintf (stdout,"\nNormalising data by map maximum:\n");
fflush (stdout);
}
/* re-open output map for reading */
fd = G_open_cell_old (parm.output->answer, G_mapset());
if (fd < 0) {
G_fatal_error ("Could not open output map for normalisation pass!\n");
}
G_strcpy (sysstr, parm.output->answer);
G_strcat (sysstr, ".globalnorm");
fd_out = G_open_raster_new (sysstr, DCELL_TYPE);
if ( fd_out < 0 ) {
G_fatal_error ("Could not open temporary map for normalisation pass (write access)!\n");
}
/* normalise evidence within map */
for (y = 0; y < nrows; y ++) {
G_get_d_raster_row (fd, diskrow, y);
for (x = 0; x < ncols; x++) {
/* write one cell to output map */
if (!G_is_d_null_value (&diskrow[x])) {
if (flag.absolute->answer) {
outrow [x] = (diskrow [x] - min) / ( max - min );
} else {
if ( min < 0 ) {
if ( diskrow [x] < 0 ) {
outrow [x] = diskrow [x] / fabs (min);
}
if ( diskrow [x] > 0 ) {
outrow [x] = diskrow [x] / max;
}
if ( diskrow [x] == 0 ) {
outrow [x] = 0;
}
} else {
outrow [x] = (diskrow [x] - min) / ( max - min );
}
}
} else {
G_set_d_null_value ( &outrow[x],1 );
}
}
/* write one row to output map on disk */
G_put_raster_row (fd_out, outrow, DCELL_TYPE);
if (!flag.quiet->answer) {
G_percent (y, nrows-1, 1);
fflush (stdout);
}
}
G_close_cell (fd);
G_close_cell (fd_out);
}
/* if normalisation was done, we delete the original map, and rename the normalised one */
if ( ( flag.localnorm->answer ) || ( flag.globalnorm->answer ) ) {
G_remove ("cats", parm.output->answer);
G_remove ("cell", parm.output->answer);
G_remove ("cell_hd", parm.output->answer);
G_remove ("cell_misc", parm.output->answer);
G_remove ("colr", parm.output->answer);
G_remove ("colr2", parm.output->answer);
G_remove ("fcell", parm.output->answer);
G_remove ("hist", parm.output->answer);
G_rename ("cats", sysstr, parm.output->answer);
G_rename ("cell", sysstr, parm.output->answer);
G_rename ("cell_hd", sysstr, parm.output->answer);
G_rename ("cell_misc", sysstr, parm.output->answer);
G_rename ("colr", sysstr, parm.output->answer);
G_rename ("colr2", sysstr, parm.output->answer);
G_rename ("fcell", sysstr, parm.output->answer);
G_rename ("hist", sysstr, parm.output->answer);
}
/* if the data was normalised, we need to rescale to get a good grey range */
if ((flag.localnorm->answer) || (flag.globalnorm->answer)) {
if ( flag.absolute->answer) {
min = 0;
max = 1.0;
} else {
if ( min < 0 ) {
min = -1.0;
} else {
min = 0.0;
}
max = 1.0;
}
}
/* write a colormap in grey shades: lowest prominence gets black, highest white */
G_init_colors (&colors);
if ( (flag.absolute->answer) ) {
from = min;
to = max/2;
G_add_d_raster_color_rule (&from, 0, 0, 0, &to, 127, 127, 127, &colors);
from = max/2;
to = max;
G_add_d_raster_color_rule (&from, 128, 128, 128, &to, 255, 255, 255, &colors);
} else {
if ( min < 0 ) {
from = min;
to = min/2;
G_add_d_raster_color_rule (&from, 255, 0, 0, &to, 127, 0, 0, &colors);
from = min/2;
to = 0;
G_add_d_raster_color_rule (&from, 127, 0, 0, &to, 0, 0, 0, &colors);
from = 0;
}
else {
from = min;
}
to = max/2;
G_add_d_raster_color_rule (&from, 0, 0, 0, &to, 127, 127, 127, &colors);
from = max/2;
to = max;
G_add_d_raster_color_rule (&from, 128, 128, 128, &to, 255, 255, 255, &colors);
}
G_write_colors (parm.output->answer, G_mapset(), &colors);
G_free_colors (&colors);
G_free (diskrow);
G_free (outrow);
free (mapset);
if (!flag.quiet->answer) {
fprintf (stdout,"\nDone.\n");
fflush (stdout);
}
return (0);
}
MODULE_TOPDIR = ../..
PGM = r.prominence
GTLIB = -lgt
LIBES = $(GISLIB) $(DATETIMELIB)
EXTRA_CFLAGS = -I../dst/include \
include $(MODULE_TOPDIR)/include/Make/Module.make
default: cmd
_______________________________________________
grass-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-dev
