[GRASS-SVN] r59166 - grass/trunk/raster/r.li/r.li.edgedensity
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Mar 2 12:43:31 PST 2014
Author: mmetz
Date: 2014-03-02 12:43:30 -0800 (Sun, 02 Mar 2014)
New Revision: 59166
Modified:
grass/trunk/raster/r.li/r.li.edgedensity/edgedensity.c
grass/trunk/raster/r.li/r.li.edgedensity/r.li.edgedensity.html
Log:
r.li.edgedensity: rewrite, optimize, sync to fragstats
Modified: grass/trunk/raster/r.li/r.li.edgedensity/edgedensity.c
===================================================================
--- grass/trunk/raster/r.li/r.li.edgedensity/edgedensity.c 2014-03-02 16:09:18 UTC (rev 59165)
+++ grass/trunk/raster/r.li/r.li.edgedensity/edgedensity.c 2014-03-02 20:43:30 UTC (rev 59166)
@@ -9,6 +9,22 @@
*
*/
+/****************************************************************************
+ *
+ * MODULE: r.li.edgedensity
+ * AUTHOR(S): Serena Pallecchi student of Computer Science University of Pisa (Italy)
+ * Commission from Faunalia Pontedera (PI) www.faunalia.it
+ * Rewrite: Markus Metz
+ *
+ * PURPOSE: calculates edge density index
+ * COPYRIGHT: (C) 2006-2014 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
#include <stdlib.h>
#include <fcntl.h> /* for O_RDONLY usage */
#include <math.h>
@@ -22,9 +38,9 @@
#include "../r.li.daemon/daemon.h"
rli_func edgedensity;
-int calculate(int fd, struct area_entry *ad, char **valore, double *result);
-int calculateD(int fd, struct area_entry *ad, char **valore, double *result);
-int calculateF(int fd, struct area_entry *ad, char **valore, double *result);
+int calculate(int fd, struct area_entry *ad, char **par, double *result);
+int calculateD(int fd, struct area_entry *ad, char **par, double *result);
+int calculateF(int fd, struct area_entry *ad, char **par, double *result);
int main(int argc, char *argv[])
{
@@ -74,28 +90,25 @@
output->answer);
}
-int edgedensity(int fd, char **valore, struct area_entry *ad, double *result)
+int edgedensity(int fd, char **par, struct area_entry *ad, double *result)
{
- struct Cell_head hd;
int ris = -1;
double indice = 0;
- Rast_get_cellhd(ad->raster, "", &hd);
-
switch (ad->data_type) {
case CELL_TYPE:
{
- ris = calculate(fd, ad, valore, &indice);
+ ris = calculate(fd, ad, par, &indice);
break;
}
case DCELL_TYPE:
{
- ris = calculateD(fd, ad, valore, &indice);
+ ris = calculateD(fd, ad, par, &indice);
break;
}
case FCELL_TYPE:
{
- ris = calculateF(fd, ad, valore, &indice);
+ ris = calculateF(fd, ad, par, &indice);
break;
}
default:
@@ -114,65 +127,38 @@
}
-int calculate(int fd, struct area_entry *ad, char **valore, double *result)
+int calculate(int fd, struct area_entry *ad, char **par, double *result)
{
- double e = 0;
- double somma = 0;
- double area = 0;
- CELL *buf_corr, *buf_sup, *buf_inf, *buf_null;
- CELL prevCell, corrCell, supCell, infCell, nextCell;
- AVL_table *array;
- long tot = 0;
- long zero = 0;
- long m = 0;
- long bordoCorr = 0;
- avl_tree albero = NULL;
+ CELL *buf, *buf_sup, *buf_null;
+ CELL corrCell, precCell, supCell;
+ CELL ptype;
+ long nedges, area;
int i, j;
- int mask_fd = -1, *mask_corr, *mask_sup, *mask_inf, *mask_tmp;
- int masked = FALSE;
- int ris;
- generic_cell c1;
+ int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
+ struct Cell_head hd;
- buf_sup = NULL;
- c1.t = CELL_TYPE;
+ Rast_get_window(&hd);
/* open mask if needed */
- mask_corr = mask_sup = mask_inf = NULL;
+ mask_fd = -1;
+ mask_buf = mask_sup = NULL;
+ masked = FALSE;
if (ad->mask == 1) {
-
- if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0) {
- G_fatal_error("Cannot open mask file <%s>", ad->mask_name);
- return RLI_ERRORE; /* FIXME: can not return from a fatal error */
- }
-
- mask_corr = G_malloc(ad->cl * sizeof(int));
- if (mask_corr == NULL) {
- G_fatal_error("malloc mask_corr failed");
+ if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
return RLI_ERRORE;
- }
-
- mask_inf = G_malloc(ad->cl * sizeof(int));
- if (mask_inf == NULL) {
- G_fatal_error("malloc mask_inf failed");
+ mask_buf = G_malloc(ad->cl * sizeof(int));
+ if (mask_buf == NULL) {
+ G_fatal_error("malloc mask_buf failed");
return RLI_ERRORE;
}
-
mask_sup = G_malloc(ad->cl * sizeof(int));
if (mask_sup == NULL) {
- G_fatal_error("malloc mask_sup failed");
+ G_fatal_error("malloc mask_buf failed");
return RLI_ERRORE;
}
+ for (j = 0; j < ad->cl; j++)
+ mask_buf[j] = 0;
- for (i = 0; i < ad->cl; i++) {
- mask_corr[i] = 0;
- mask_sup[i] = 0;
- }
-
- if (read(mask_fd, mask_inf, (ad->cl * sizeof(int))) < 0) {
- G_fatal_error("reading mask_inf at first row");
- return RLI_ERRORE;
- }
-
masked = TRUE;
}
@@ -185,252 +171,172 @@
/* the first time buf_sup is all null */
Rast_set_c_null_value(buf_null, Rast_window_cols());
buf_sup = buf_null;
- buf_inf = buf_null;
+ if (par != NULL) { /* only 1 class */
+ char *sval;
+
+ sval = par[0];
+ ptype = atoi(sval);
+ }
+ else
+ Rast_set_c_null_value(&ptype, 1);
+
+ nedges = 0;
+ area = 0;
+
/* for each raster row */
- for (j = 0; j < ad->rl; j++) {
+ for (i = 0; i < ad->rl; i++) {
/* read row of raster */
- buf_corr = RLI_get_cell_raster_row(fd, j + ad->y, ad);
+ buf = RLI_get_cell_raster_row(fd, i + ad->y, ad);
- if (j > 0) /* not first row */
- buf_sup = RLI_get_cell_raster_row(fd, j - 1 + ad->y, ad);
+ if (i > 0) /* not first row */
+ buf_sup = RLI_get_cell_raster_row(fd, i - 1 + ad->y, ad);
- if ((j + 1) < ad->rl) { /* not last row */
- buf_inf = RLI_get_cell_raster_row(fd, 1 + j + ad->y, ad);
- }
- else {
- buf_inf = buf_null;
- }
-
/* read mask if needed */
if (masked) {
mask_tmp = mask_sup;
- mask_sup = mask_corr;
- mask_corr = mask_inf;
- mask_inf = mask_tmp;
-
- if ((j + 1) < ad->rl) { /* not last row */
- if (read(mask_fd, mask_inf, (ad->cl * sizeof(int))) < 0) {
- G_fatal_error("reading mask_inf at row j");
- return RLI_ERRORE;
- }
- }
- else {
- int z = 0;
-
- for (z = 0; z < ad->cl; z++)
- mask_inf[z] = 0;
- }
+ mask_sup = mask_buf;
+ mask_buf = mask_tmp;
+ if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
+ return RLI_ERRORE;
}
- Rast_set_c_null_value(&nextCell, 1);
- Rast_set_c_null_value(&prevCell, 1);
- Rast_set_c_null_value(&corrCell, 1);
+ Rast_set_c_null_value(&precCell, 1);
- for (i = 0; i < ad->cl; i++) {
- corrCell = buf_corr[i + ad->x];
+ for (j = 0; j < ad->cl; j++) {
+ corrCell = buf[j + ad->x];
- if (masked && mask_corr[i] == 0) {
+ if (masked && mask_buf[j] == 0) {
Rast_set_c_null_value(&corrCell, 1);
}
-
- if (!(Rast_is_null_value(&corrCell, CELL_TYPE))) {
+ else {
+ /* total sample area */
area++;
- if ((i + 1) == ad->cl) /*last cell of the row */
- Rast_set_c_null_value(&nextCell, 1);
- else if (masked && mask_corr[i + 1] == 0)
- Rast_set_c_null_value(&nextCell, 1);
- else
- nextCell = buf_corr[i + 1 + ad->x];
+ }
- if (masked && mask_sup[i] == 0)
- Rast_set_c_null_value(&supCell, 1);
- else
- supCell = buf_sup[i + ad->x];
+ supCell = buf_sup[j + ad->x];
+ if (masked && (mask_sup[j] == 0)) {
+ Rast_set_c_null_value(&supCell, 1);
+ }
- if (masked && mask_inf[i] == 0)
- Rast_set_c_null_value(&infCell, 1);
- else
- infCell = buf_inf[i + ad->x];
-
- /* calculate how many edges the cell has */
- if ((Rast_is_null_value(&prevCell, CELL_TYPE)) ||
- (corrCell != prevCell)) {
- bordoCorr++;
+ if (!Rast_is_c_null_value(&ptype)) {
+ /* only one patch type */
+ if (!Rast_is_c_null_value(&corrCell) && corrCell == ptype) {
+ if (corrCell != precCell)
+ nedges++;
+ if (corrCell != supCell)
+ nedges++;
+ /* right and bottom */
+ if (i == ad->rl - 1)
+ nedges++;
+ if (j == ad->cl - 1)
+ nedges++;
}
-
- if ((Rast_is_null_value(&supCell, CELL_TYPE)) ||
- (corrCell != supCell)) {
- bordoCorr++;
+ if (!Rast_is_c_null_value(&precCell) && precCell == ptype) {
+ if (corrCell != precCell)
+ nedges++;
}
-
- if ((Rast_is_null_value(&infCell, CELL_TYPE)) ||
- (corrCell != infCell)) {
- bordoCorr++;
+ if (!Rast_is_c_null_value(&supCell) && supCell == ptype) {
+ if (corrCell != supCell)
+ nedges++;
}
-
- if ((Rast_is_null_value(&nextCell, CELL_TYPE)) ||
- (corrCell != nextCell)) {
- bordoCorr++;
- }
-
- /* store the result in the tree */
- if (albero == NULL) {
- c1.val.c = corrCell;
- albero = avl_make(c1, bordoCorr);
- if (albero == NULL) {
- G_fatal_error("avl_make error");
- return RLI_ERRORE;
+ }
+ else {
+ /* all patch types */
+ if (!Rast_is_c_null_value(&corrCell)) {
+ if (corrCell != precCell) {
+ nedges++;
}
- m++;
- }
- else {
- c1.val.c = corrCell;
- ris = avl_add(&albero, c1, bordoCorr);
-
- switch (ris) {
- case AVL_ERR:
- {
- G_fatal_error("avl_add error");
- return RLI_ERRORE;
- }
- case AVL_ADD:
- {
- m++;
- break;
- }
- case AVL_PRES:
- {
- break;
- }
- default:
- {
- G_fatal_error("avl_add unknown error");
- return RLI_ERRORE;
- }
+ if (corrCell != supCell) {
+ nedges++;
}
+ /* right and bottom */
+ if (i == ad->rl - 1)
+ nedges++;
+ if (j == ad->cl - 1)
+ nedges++;
}
- bordoCorr = 0;
+ if (!Rast_is_c_null_value(&precCell) && corrCell != precCell) {
+ nedges++;
+ }
+ if (!Rast_is_c_null_value(&supCell) && corrCell != supCell) {
+ nedges++;
+ }
}
- prevCell = corrCell;
+ precCell = corrCell;
}
}
/* calculate index */
if (area > 0) {
- if (valore != NULL) { /* only 1 class */
- char *sval;
- int val;
- CELL cella;
+ double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
+ double elength, cell_size;
- sval = valore[0];
- val = atoi(sval);
- cella = val;
- c1.t = CELL_TYPE;
- c1.val.c = cella;
- e = (double)howManyCell(albero, c1);
- somma = e;
- }
- else { /* all classes */
- array = G_malloc(m * sizeof(AVL_tableRow));
- if (array == NULL) {
- G_fatal_error("malloc array failed");
- return RLI_ERRORE;
- }
- tot = avl_to_array(albero, zero, array);
- if (tot != m) {
- G_warning
- ("avl_to_array unexpected value. the result could be wrong");
- }
- for (i = 0; i < m; i++) {
- e = (double)array[i]->tot;
- somma = somma + e;
- }
- G_free(array);
- }
+ /* calculate distance */
+ G_begin_distance_calculations();
+ /* EW Dist at North edge */
+ EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
+ /* EW Dist at South Edge */
+ EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
+ /* NS Dist at East edge */
+ NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
+ /* NS Dist at West edge */
+ NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
- *result = somma * 10000 / area;
+ elength = ((EW_DIST1 + EW_DIST2) / (2 * hd.cols) +
+ (NS_DIST1 + NS_DIST2) / (2 * hd.rows)) / 2;
+
+ cell_size = ((EW_DIST1 + EW_DIST2) / (2 * hd.cols)) *
+ ((NS_DIST1 + NS_DIST2) / (2 * hd.rows));
+
+ *result = (double) nedges * elength * 10000 / (area * cell_size);
}
else
Rast_set_d_null_value(result, 1);
if (masked) {
close(mask_fd);
- G_free(mask_inf);
- G_free(mask_corr);
+ G_free(mask_buf);
G_free(mask_sup);
}
-
G_free(buf_null);
- avl_destroy(albero);
return RLI_OK;
}
-int calculateD(int fd, struct area_entry *ad, char **valore, double *result)
+int calculateD(int fd, struct area_entry *ad, char **par, double *result)
{
- double e = 0;
- double somma = 0;
- double area = 0;
-
- DCELL *buf_corr, *buf_sup, *buf_inf, *buf_null;
- DCELL prevCell, corrCell, supCell, infCell, nextCell;
-
- AVL_table *array;
-
- long tot = 0;
- long zero = 0;
- long m = 0;
- long bordoCorr = 0;
-
- avl_tree albero = NULL;
-
+ DCELL *buf, *buf_sup, *buf_null;
+ DCELL corrCell, precCell, supCell;
+ DCELL ptype;
+ long nedges, area;
int i, j;
- int mask_fd = -1, *mask_corr, *mask_sup, *mask_inf, *mask_tmp;
- int masked = FALSE;
- int ris;
+ int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
+ struct Cell_head hd;
- generic_cell c1;
+ Rast_get_window(&hd);
- c1.t = DCELL_TYPE;
-
/* open mask if needed */
- mask_corr = mask_sup = mask_inf = NULL;
+ mask_fd = -1;
+ mask_buf = mask_sup = NULL;
+ masked = FALSE;
if (ad->mask == 1) {
- if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0) {
- G_fatal_error("Cannot open mask file");
- return RLI_ERRORE; /* FIXME: can not return from a fatal error */
- }
-
- mask_corr = G_malloc(ad->cl * sizeof(int));
- if (mask_corr == NULL) {
- G_fatal_error("malloc mask_corr failed");
+ if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
return RLI_ERRORE;
- }
-
- mask_inf = G_malloc(ad->cl * sizeof(int));
- if (mask_inf == NULL) {
- G_fatal_error("malloc mask_inf failed");
+ mask_buf = G_malloc(ad->cl * sizeof(int));
+ if (mask_buf == NULL) {
+ G_fatal_error("malloc mask_buf failed");
return RLI_ERRORE;
}
-
mask_sup = G_malloc(ad->cl * sizeof(int));
if (mask_sup == NULL) {
- G_fatal_error("malloc mask_sup failed");
+ G_fatal_error("malloc mask_buf failed");
return RLI_ERRORE;
}
+ for (j = 0; j < ad->cl; j++)
+ mask_buf[j] = 0;
- for (i = 0; i < ad->cl; i++) {
- mask_corr[i] = 0;
- mask_sup[i] = 0;
- }
-
- if (read(mask_fd, mask_inf, (ad->cl * sizeof(int))) <= 0) {
- G_fatal_error("reading mask_corr");
- return RLI_ERRORE;
- }
-
masked = TRUE;
}
@@ -443,252 +349,172 @@
/* the first time buf_sup is all null */
Rast_set_d_null_value(buf_null, Rast_window_cols());
buf_sup = buf_null;
- buf_inf = buf_null;
+ if (par != NULL) { /* only 1 class */
+ char *sval;
+
+ sval = par[0];
+ ptype = atof(sval);
+ }
+ else
+ Rast_set_d_null_value(&ptype, 1);
+
+ nedges = 0;
+ area = 0;
+
/* for each raster row */
- for (j = 0; j < ad->rl; j++) {
+ for (i = 0; i < ad->rl; i++) {
+
/* read row of raster */
- buf_corr = RLI_get_dcell_raster_row(fd, j + ad->y, ad);
+ buf = RLI_get_dcell_raster_row(fd, i + ad->y, ad);
- /* not first row */
- if (j > 0)
- buf_sup = RLI_get_dcell_raster_row(fd, j - 1 + ad->y, ad);
+ if (i > 0) /* not first row */
+ buf_sup = RLI_get_dcell_raster_row(fd, i - 1 + ad->y, ad);
- if ((j + 1) < ad->rl) { /* not last row */
- buf_inf = RLI_get_dcell_raster_row(fd, 1 + j + ad->y, ad);
- }
- else {
- buf_inf = buf_null;
- }
-
/* read mask if needed */
if (masked) {
-
mask_tmp = mask_sup;
- mask_sup = mask_corr;
- mask_corr = mask_inf;
- mask_inf = mask_tmp;
-
- if ((j + 1) < ad->rl) { /* not last row */
- if (read(mask_fd, mask_inf, (ad->cl * sizeof(int))) < 0) {
- G_fatal_error("reading mask_inf");
- return RLI_ERRORE;
- }
- }
- else {
- int z = 0;
-
- for (z = 0; z < ad->cl; z++)
- mask_inf[z] = 0;
- }
+ mask_sup = mask_buf;
+ mask_buf = mask_tmp;
+ if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
+ return RLI_ERRORE;
}
- Rast_set_d_null_value(&nextCell, 1);
- Rast_set_d_null_value(&prevCell, 1);
- Rast_set_d_null_value(&corrCell, 1);
+ Rast_set_d_null_value(&precCell, 1);
- for (i = 0; i < ad->cl; i++) { /* for each cell in the row */
- corrCell = buf_corr[i + ad->x];
+ for (j = 0; j < ad->cl; j++) {
+ corrCell = buf[j + ad->x];
- if (masked && mask_corr[i] == 0) {
+ if (masked && mask_buf[j] == 0) {
Rast_set_d_null_value(&corrCell, 1);
}
-
- if (!(Rast_is_null_value(&corrCell, DCELL_TYPE))) {
+ else {
+ /* total sample area */
area++;
- if ((i + 1) == ad->cl) /*last cell of the row */
- Rast_set_d_null_value(&nextCell, 1);
- else if (masked && mask_corr[i + 1] == 0)
- Rast_set_d_null_value(&nextCell, 1);
- else
- nextCell = buf_corr[i + 1 + ad->x];
+ }
- if (masked && mask_sup[i] == 0)
- Rast_set_d_null_value(&supCell, 1);
- else
- supCell = buf_sup[i + ad->x];
+ supCell = buf_sup[j + ad->x];
+ if (masked && (mask_sup[j] == 0)) {
+ Rast_set_d_null_value(&supCell, 1);
+ }
- if (masked && mask_inf[i] == 0)
- Rast_set_d_null_value(&infCell, 1);
- else
- infCell = buf_inf[i + ad->x];
-
- /* calculate how many edges the cell has */
-
- if ((Rast_is_null_value(&prevCell, DCELL_TYPE)) ||
- (corrCell != prevCell)) {
- bordoCorr++;
+ if (!Rast_is_d_null_value(&ptype)) {
+ /* only one patch type */
+ if (!Rast_is_d_null_value(&corrCell) && corrCell == ptype) {
+ if (corrCell != precCell)
+ nedges++;
+ if (corrCell != supCell)
+ nedges++;
+ /* right and bottom */
+ if (i == ad->rl - 1)
+ nedges++;
+ if (j == ad->cl - 1)
+ nedges++;
}
-
- if ((Rast_is_null_value(&supCell, DCELL_TYPE)) ||
- (corrCell != supCell)) {
- bordoCorr++;
+ if (!Rast_is_d_null_value(&precCell) && precCell == ptype) {
+ if (corrCell != precCell)
+ nedges++;
}
-
- if ((Rast_is_null_value(&infCell, DCELL_TYPE)) ||
- (corrCell != infCell)) {
- bordoCorr++;
+ if (!Rast_is_d_null_value(&supCell) && supCell == ptype) {
+ if (corrCell != supCell)
+ nedges++;
}
-
- if ((Rast_is_null_value(&nextCell, DCELL_TYPE)) ||
- (corrCell != nextCell)) {
- bordoCorr++;
- }
-
- /* store the result in the tree */
- if (albero == NULL) {
- c1.val.dc = corrCell;
- albero = avl_make(c1, bordoCorr);
- if (albero == NULL) {
- G_fatal_error("avl_make error");
- return RLI_ERRORE;
+ }
+ else {
+ /* all patch types */
+ if (!Rast_is_d_null_value(&corrCell)) {
+ if (corrCell != precCell) {
+ nedges++;
}
- m++;
- }
- else {
- c1.val.dc = corrCell;
- ris = avl_add(&albero, c1, bordoCorr);
-
- switch (ris) {
- case AVL_ERR:
- {
- G_fatal_error("avl_add error");
- return RLI_ERRORE;
- }
- case AVL_ADD:
- {
- m++;
- break;
- }
- case AVL_PRES:
- {
- break;
- }
- default:
- {
- G_fatal_error("avl_add unknown error");
- return RLI_ERRORE;
- }
+ if (corrCell != supCell) {
+ nedges++;
}
+ /* right and bottom */
+ if (i == ad->rl - 1)
+ nedges++;
+ if (j == ad->cl - 1)
+ nedges++;
}
- bordoCorr = 0;
+ if (!Rast_is_d_null_value(&precCell) && corrCell != precCell) {
+ nedges++;
+ }
+ if (!Rast_is_d_null_value(&supCell) && corrCell != supCell) {
+ nedges++;
+ }
}
- prevCell = corrCell;
+ precCell = corrCell;
}
}
/* calculate index */
if (area > 0) {
- if (valore != NULL) { /* only 1 class */
- char *sval;
- double val;
- DCELL cella;
+ double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
+ double elength, cell_size;
- sval = valore[0];
- val = (double)atof(sval);
- cella = val;
- c1.val.dc = cella;
- c1.t = DCELL_TYPE;
- e = (double)howManyCell(albero, c1);
- somma = e;
- }
- else { /* all classes */
- array = G_malloc(m * sizeof(AVL_tableRow));
- if (array == NULL) {
- G_fatal_error("malloc array failed");
- return RLI_ERRORE;
- }
- tot = avl_to_array(albero, zero, array);
- if (tot != m) {
- G_warning
- ("avl_to_array unexpected value. the result could be wrong");
- }
- for (i = 0; i < m; i++) {
- e = (double)array[i]->tot;
- somma = somma + e;
- }
- G_free(array);
- }
- *result = somma * 10000 / area;
+ /* calculate distance */
+ G_begin_distance_calculations();
+ /* EW Dist at North edge */
+ EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
+ /* EW Dist at South Edge */
+ EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
+ /* NS Dist at East edge */
+ NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
+ /* NS Dist at West edge */
+ NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
+
+ elength = ((EW_DIST1 + EW_DIST2) / (2 * hd.cols) +
+ (NS_DIST1 + NS_DIST2) / (2 * hd.rows)) / 2;
+
+ cell_size = ((EW_DIST1 + EW_DIST2) / (2 * hd.cols)) *
+ ((NS_DIST1 + NS_DIST2) / (2 * hd.rows));
+
+ *result = (double) nedges * elength * 10000 / (area * cell_size);
}
else
Rast_set_d_null_value(result, 1);
if (masked) {
close(mask_fd);
- G_free(mask_inf);
- G_free(mask_corr);
+ G_free(mask_buf);
G_free(mask_sup);
}
-
G_free(buf_null);
- avl_destroy(albero);
return RLI_OK;
}
-int calculateF(int fd, struct area_entry *ad, char **valore, double *result)
+int calculateF(int fd, struct area_entry *ad, char **par, double *result)
{
- double e = 0;
- double somma = 0;
- double area = 0;
-
- FCELL *buf_corr, *buf_sup, *buf_inf, *buf_null;
- FCELL prevCell, corrCell, supCell, infCell, nextCell;
-
- AVL_table *array;
-
- long tot = 0;
- long zero = 0;
- long m = 0;
- long bordoCorr = 0;
-
- avl_tree albero = NULL;
-
+ FCELL *buf, *buf_sup, *buf_null;
+ FCELL corrCell, precCell, supCell;
+ FCELL ptype;
+ long nedges, area;
int i, j;
- int mask_fd = -1, *mask_corr, *mask_sup, *mask_inf, *mask_tmp;
- int masked = FALSE;
- int ris;
+ int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
+ struct Cell_head hd;
- generic_cell c1;
+ Rast_get_window(&hd);
- c1.t = FCELL_TYPE;
/* open mask if needed */
- mask_corr = mask_sup = mask_inf = NULL;
+ mask_fd = -1;
+ mask_buf = mask_sup = NULL;
+ masked = FALSE;
if (ad->mask == 1) {
- if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0) {
- G_fatal_error("Cannot open mask file");
- return RLI_ERRORE; /* FIXME: can not return from a fatal error */
- }
-
- mask_corr = G_malloc(ad->cl * sizeof(int));
- if (mask_corr == NULL) {
- G_fatal_error("malloc mask_corr failed");
+ if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
return RLI_ERRORE;
- }
-
- mask_inf = G_malloc(ad->cl * sizeof(int));
- if (mask_inf == NULL) {
- G_fatal_error("malloc mask_inf failed");
+ mask_buf = G_malloc(ad->cl * sizeof(int));
+ if (mask_buf == NULL) {
+ G_fatal_error("malloc mask_buf failed");
return RLI_ERRORE;
}
-
mask_sup = G_malloc(ad->cl * sizeof(int));
if (mask_sup == NULL) {
- G_fatal_error("malloc mask_sup failed");
+ G_fatal_error("malloc mask_buf failed");
return RLI_ERRORE;
}
+ for (j = 0; j < ad->cl; j++)
+ mask_buf[j] = 0;
- for (i = 0; i < ad->cl; i++) {
- mask_corr[i] = 0;
- mask_sup[i] = 0;
- }
-
- if (read(mask_fd, mask_inf, (ad->cl * sizeof(int))) <= 0) {
- G_fatal_error("reading mask_corr");
- return RLI_ERRORE;
- }
-
masked = TRUE;
}
@@ -701,184 +527,136 @@
/* the first time buf_sup is all null */
Rast_set_f_null_value(buf_null, Rast_window_cols());
buf_sup = buf_null;
- buf_inf = buf_null;
+ if (par != NULL) { /* only 1 class */
+ char *sval;
+
+ sval = par[0];
+ ptype = atof(sval);
+ }
+ else
+ Rast_set_f_null_value(&ptype, 1);
+
+ nedges = 0;
+ area = 0;
+
/* for each raster row */
- for (j = 0; j < ad->rl; j++) {
+ for (i = 0; i < ad->rl; i++) {
+
/* read row of raster */
- buf_corr = RLI_get_fcell_raster_row(fd, j + ad->y, ad);
+ buf = RLI_get_fcell_raster_row(fd, i + ad->y, ad);
- if (j > 0) /* not first row */
- buf_sup = RLI_get_fcell_raster_row(fd, j - 1 + ad->y, ad);
+ if (i > 0) /* not first row */
+ buf_sup = RLI_get_fcell_raster_row(fd, i - 1 + ad->y, ad);
- if ((j + 1) < ad->rl) { /* not last row */
- buf_inf = RLI_get_fcell_raster_row(fd, 1 + j + ad->y, ad);
- }
- else {
- buf_inf = buf_null;
- }
-
/* read mask if needed */
if (masked) {
mask_tmp = mask_sup;
- mask_sup = mask_corr;
- mask_corr = mask_inf;
- mask_inf = mask_tmp;
-
- if ((j + 1) < ad->rl) { /* not last row */
- if (read(mask_fd, mask_inf, (ad->cl * sizeof(int))) < 0) {
- G_fatal_error("reading mask_inf");
- return RLI_ERRORE;
- }
- }
- else {
- int z = 0;
-
- for (z = 0; z < ad->cl; z++)
- mask_inf[z] = 0;
- }
+ mask_sup = mask_buf;
+ mask_buf = mask_tmp;
+ if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
+ return RLI_ERRORE;
}
- Rast_set_f_null_value(&nextCell, 1);
- Rast_set_f_null_value(&prevCell, 1);
- Rast_set_f_null_value(&corrCell, 1);
+ Rast_set_f_null_value(&precCell, 1);
- for (i = 0; i < ad->cl; i++) { /* for each cell in the row */
- corrCell = buf_corr[i + ad->x];
+ for (j = 0; j < ad->cl; j++) {
+ corrCell = buf[j + ad->x];
- if (masked && mask_corr[i] == 0) {
+ if (masked && mask_buf[j] == 0) {
Rast_set_f_null_value(&corrCell, 1);
}
-
- if (!(Rast_is_null_value(&corrCell, FCELL_TYPE))) {
+ else {
+ /* total sample area */
area++;
- if ((i + 1) == ad->cl) /*last cell of the row */
- Rast_set_f_null_value(&nextCell, 1);
- else if (masked && mask_corr[i + 1] == 0)
- Rast_set_f_null_value(&nextCell, 1);
- else
- nextCell = buf_corr[i + 1 + ad->x];
+ }
- if (masked && mask_sup[i] == 0)
- Rast_set_f_null_value(&supCell, 1);
- else
- supCell = buf_sup[i + ad->x];
+ supCell = buf_sup[j + ad->x];
+ if (masked && (mask_sup[j] == 0)) {
+ Rast_set_f_null_value(&supCell, 1);
+ }
- if (masked && mask_inf[i] == 0)
- Rast_set_f_null_value(&infCell, 1);
- else
- infCell = buf_inf[i + ad->x];
-
- /* calculate how many edges the cell has */
- if ((Rast_is_null_value(&prevCell, FCELL_TYPE)) ||
- (corrCell != prevCell)) {
- bordoCorr++;
+ if (!Rast_is_f_null_value(&ptype)) {
+ /* only one patch type */
+ if (!Rast_is_f_null_value(&corrCell) && corrCell == ptype) {
+ if (corrCell != precCell)
+ nedges++;
+ if (corrCell != supCell)
+ nedges++;
+ /* right and bottom */
+ if (i == ad->rl - 1)
+ nedges++;
+ if (j == ad->cl - 1)
+ nedges++;
}
-
- if ((Rast_is_null_value(&supCell, FCELL_TYPE)) ||
- (corrCell != supCell)) {
- bordoCorr++;
+ if (!Rast_is_f_null_value(&precCell) && precCell == ptype) {
+ if (corrCell != precCell)
+ nedges++;
}
-
- if ((Rast_is_null_value(&infCell, FCELL_TYPE)) ||
- (corrCell != infCell)) {
- bordoCorr++;
+ if (!Rast_is_f_null_value(&supCell) && supCell == ptype) {
+ if (corrCell != supCell)
+ nedges++;
}
-
- if ((Rast_is_null_value(&nextCell, FCELL_TYPE)) ||
- (corrCell != nextCell)) {
- bordoCorr++;
- }
-
- /* store the result in the tree */
- if (albero == NULL) {
- c1.val.fc = corrCell;
- albero = avl_make(c1, bordoCorr);
- if (albero == NULL) {
- G_fatal_error("avl_make error");
- return RLI_ERRORE;
+ }
+ else {
+ /* all patch types */
+ if (!Rast_is_f_null_value(&corrCell)) {
+ if (corrCell != precCell) {
+ nedges++;
}
- m++;
- }
- else {
- c1.val.fc = corrCell;
- ris = avl_add(&albero, c1, bordoCorr);
-
- switch (ris) {
- case AVL_ERR:
- {
- G_fatal_error("avl_add error");
- return RLI_ERRORE;
- }
- case AVL_ADD:
- {
- m++;
- break;
- }
- case AVL_PRES:
- {
- break;
- }
- default:
- {
- G_fatal_error("avl_add unknown error");
- return RLI_ERRORE;
- }
+ if (corrCell != supCell) {
+ nedges++;
}
+ /* right and bottom */
+ if (i == ad->rl - 1)
+ nedges++;
+ if (j == ad->cl - 1)
+ nedges++;
}
- bordoCorr = 0;
+ if (!Rast_is_f_null_value(&precCell) && corrCell != precCell) {
+ nedges++;
+ }
+ if (!Rast_is_f_null_value(&supCell) && corrCell != supCell) {
+ nedges++;
+ }
}
- prevCell = corrCell;
+ precCell = corrCell;
}
}
/* calculate index */
if (area > 0) {
- if (valore != NULL) { /* only 1 class */
- char *sval;
- float val;
- FCELL cella;
+ double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
+ double elength, cell_size;
- sval = valore[0];
- val = (float)atof(sval);
- cella = val;
- c1.t = FCELL_TYPE;
- c1.val.fc = cella;
- e = (double)howManyCell(albero, c1);
- somma = e;
+ /* calculate distance */
+ G_begin_distance_calculations();
+ /* EW Dist at North edge */
+ EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
+ /* EW Dist at South Edge */
+ EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
+ /* NS Dist at East edge */
+ NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
+ /* NS Dist at West edge */
+ NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
- }
- else { /* all classes */
- array = G_malloc(m * sizeof(AVL_tableRow));
- if (array == NULL) {
- G_fatal_error("malloc array failed");
- return RLI_ERRORE;
- }
- tot = avl_to_array(albero, zero, array);
- if (tot != m) {
- G_warning
- ("avl_to_array unexpected value. the result could be wrong");
- }
- for (i = 0; i < m; i++) {
- e = (double)array[i]->tot;
- somma = somma + e;
- }
- G_free(array);
- }
- *result = somma * 10000 / area;
+ elength = ((EW_DIST1 + EW_DIST2) / (2 * hd.cols) +
+ (NS_DIST1 + NS_DIST2) / (2 * hd.rows)) / 2;
+
+ cell_size = ((EW_DIST1 + EW_DIST2) / (2 * hd.cols)) *
+ ((NS_DIST1 + NS_DIST2) / (2 * hd.rows));
+
+ *result = (double) nedges * elength * 10000 / (area * cell_size);
}
else
Rast_set_d_null_value(result, 1);
if (masked) {
close(mask_fd);
- G_free(mask_inf);
- G_free(mask_corr);
+ G_free(mask_buf);
G_free(mask_sup);
}
-
G_free(buf_null);
- avl_destroy(albero);
return RLI_OK;
}
Modified: grass/trunk/raster/r.li/r.li.edgedensity/r.li.edgedensity.html
===================================================================
--- grass/trunk/raster/r.li/r.li.edgedensity/r.li.edgedensity.html 2014-03-02 16:09:18 UTC (rev 59165)
+++ grass/trunk/raster/r.li/r.li.edgedensity/r.li.edgedensity.html 2014-03-02 20:43:30 UTC (rev 59166)
@@ -19,6 +19,8 @@
the landscape involving patch type k</li>
<li> <b>Area</b>: total landscape area</li>
</ul>
+<p>
+The unit is meters per hectare.
<h2>NOTES</h2>
@@ -30,13 +32,8 @@
<!-- TODO: verify next: -->
If the input raster map contains only NULL values then <em>r.li.edgedensity</em>
consider to have 0 patches.<br>
-If area is 0 <em>r.li.edgedensity</em> returns -1; this is only possible if input
+If area is 0 <em>r.li.edgedensity</em> returns NULL; this is only possible if input
raster is masked.
-If you want to change these -1 values to NULL, run subsequently on the resulting map:
-<div class="code"><pre>
-r.null setnull=-1 input=my_map
-</pre></div>
-after index calculation.
<h2>EXAMPLES</h2>
To calculate the edge density index on map <em>my_map</em>, using
More information about the grass-commit
mailing list