[GRASS-SVN] r59169 - grass/trunk/raster/r.li/r.li.shape
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Mar 2 14:34:11 PST 2014
Author: mmetz
Date: 2014-03-02 14:34:11 -0800 (Sun, 02 Mar 2014)
New Revision: 59169
Modified:
grass/trunk/raster/r.li/r.li.shape/main.c
grass/trunk/raster/r.li/r.li.shape/r.li.shape.html
Log:
r.li.shape: calculate shape index
Modified: grass/trunk/raster/r.li/r.li.shape/main.c
===================================================================
--- grass/trunk/raster/r.li/r.li.shape/main.c 2014-03-02 21:03:37 UTC (rev 59168)
+++ grass/trunk/raster/r.li/r.li.shape/main.c 2014-03-02 22:34:11 UTC (rev 59169)
@@ -18,12 +18,24 @@
#include <stdlib.h>
#include <fcntl.h>
+#include <math.h>
#include <grass/gis.h>
#include <grass/raster.h>
#include <grass/glocale.h>
#include "../r.li.daemon/daemon.h"
+#include "../r.li.daemon/GenericCell.h"
+/* type, cell count and edge count of each patch */
+struct pst {
+ generic_cell type;
+ long cells;
+ long edges;
+};
+
rli_func shape_index;
+int calculate(int fd, struct area_entry *ad, double *result);
+int calculateD(int fd, struct area_entry *ad, double *result);
+int calculateF(int fd, struct area_entry *ad, double *result);
int main(int argc, char *argv[])
{
@@ -59,55 +71,754 @@
int shape_index(int fd, char **par, struct area_entry *ad, double *result)
{
- double area;
- struct Cell_head hd;
- CELL complete_value;
- double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
- int mask_fd = -1, null_count = 0;
- int i = 0, k = 0;
- int *mask_buf = NULL;
+ int ris = RLI_OK;
+ double indice = 0;
- Rast_set_c_null_value(&complete_value, 1);
+ switch (ad->data_type) {
+ case CELL_TYPE:
+ {
+ ris = calculate(fd, ad, &indice);
+ break;
+ }
+ case DCELL_TYPE:
+ {
+ ris = calculateD(fd, ad, &indice);
+ break;
+ }
+ case FCELL_TYPE:
+ {
+ ris = calculateF(fd, ad, &indice);
+ break;
+ }
+ default:
+ {
+ G_fatal_error("data type unknown");
+ return RLI_ERRORE;
+ }
+ }
- Rast_get_window(&hd);
+ if (ris != RLI_OK) {
+ return RLI_ERRORE;
+ }
+ *result = indice;
+
+ return RLI_OK;
+}
+
+
+int calculate(int fd, struct area_entry *ad, double *result)
+{
+ CELL *buf, *buf_sup, *buf_null;
+ CELL corrCell, precCell, supCell;
+ long npatch, area;
+ long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
+ struct pst *pst;
+ long nalloc, incr;
+ int i, j, k;
+ int connected;
+ int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
+
+ buf_null = Rast_allocate_c_buf();
+ Rast_set_c_null_value(buf_null, Rast_window_cols());
+ buf_sup = buf_null;
+
+ /* initialize patch ids */
+ pid_corr = G_malloc(ad->cl * sizeof(long));
+ pid_sup = G_malloc(ad->cl * sizeof(long));
+
+ for (j = 0; j < ad->cl; j++) {
+ pid_corr[j] = 0;
+ pid_sup[j] = 0;
+ }
+
/* open mask if needed */
+ 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)
- return 0;
+ return RLI_ERRORE;
+ 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_buf failed");
+ return RLI_ERRORE;
+ }
+ for (j = 0; j < ad->cl; j++)
+ mask_buf[j] = 0;
+ masked = TRUE;
+ }
+
+ /* calculate number of patches */
+ npatch = 0;
+ area = 0;
+ pid = 0;
+
+ /* patch size and type */
+ incr = 1024;
+ if (incr > ad->rl)
+ incr = ad->rl;
+ if (incr > ad->cl)
+ incr = ad->cl;
+ if (incr < 2)
+ incr = 2;
+ nalloc = incr;
+ pst = G_malloc(nalloc * sizeof(struct pst));
+ for (k = 0; k < nalloc; k++) {
+ pst[k].cells = 0;
+ pst[k].edges = 0;
+ }
+
+ for (i = 0; i < ad->rl; i++) {
+ buf = RLI_get_cell_raster_row(fd, i + ad->y, ad);
+ if (i > 0) {
+ buf_sup = RLI_get_cell_raster_row(fd, i - 1 + ad->y, ad);
+ }
+
+ if (masked) {
+ mask_tmp = mask_sup;
+ mask_sup = mask_buf;
+ mask_buf = mask_tmp;
+ if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
+ return 0;
+ }
+
+ ltmp = pid_sup;
+ pid_sup = pid_corr;
+ pid_corr = ltmp;
+
+ Rast_set_c_null_value(&precCell, 1);
+
+ connected = 0;
+ for (j = 0; j < ad->cl; j++) {
+ pid_corr[j] = 0;
+
+ corrCell = buf[j + ad->x];
+ if (masked && (mask_buf[j] == 0)) {
+ Rast_set_c_null_value(&corrCell, 1);
+ }
+ else {
+ /* total sample area */
+ area++;
+ }
+
+ supCell = buf_sup[j + ad->x];
+ if (masked && (mask_sup[j] == 0)) {
+ Rast_set_c_null_value(&supCell, 1);
+ }
+
+ if (Rast_is_c_null_value(&corrCell)) {
+ if (!Rast_is_c_null_value(&precCell))
+ pst[pid_corr[j - 1]].edges++;
+ if (!Rast_is_c_null_value(&supCell))
+ pst[pid_sup[j]].edges++;
+ connected = 0;
+ precCell = corrCell;
+ continue;
+ }
+
+ if (!Rast_is_c_null_value(&precCell) && corrCell == precCell) {
+ pid_corr[j] = pid_corr[j - 1];
+ connected = 1;
+ pst[pid_corr[j]].cells++;
+ }
+ else {
+ connected = 0;
+ }
+
+ if (!Rast_is_c_null_value(&supCell) && corrCell == supCell) {
+
+ if (pid_corr[j] != pid_sup[j]) {
+ /* connect or merge */
+ /* after r.clump */
+ if (connected) {
+ npatch--;
+
+ if (npatch == 0) {
+ G_fatal_error("npatch == 0 at row %d, col %d", i, j);
+ }
+ }
+
+ old_pid = pid_corr[j];
+ new_pid = pid_sup[j];
+ pid_corr[j] = new_pid;
+ if (old_pid > 0) {
+ /* merge */
+ /* update left side of the current row */
+ for (k = 0; k < j; k++) {
+ if (pid_corr[k] == old_pid)
+ pid_corr[k] = new_pid;
+ }
+ /* update right side of the previous row */
+ for (k = j + 1; k < ad->cl; k++) {
+ if (pid_sup[k] == old_pid)
+ pid_sup[k] = new_pid;
+ }
+ pst[new_pid].cells += pst[old_pid].cells;
+ pst[old_pid].cells = 0;
+ pst[new_pid].edges += pst[old_pid].edges;
+ pst[old_pid].edges = 0;
+
+ if (old_pid == pid)
+ pid--;
+ }
+ else {
+ pst[new_pid].cells++;
+ }
+ }
+ connected = 1;
+ }
+
+ if (!connected) {
+ /* start new patch */
+ npatch++;
+ pid++;
+ pid_corr[j] = pid;
+
+ if (pid >= nalloc) {
+ pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
+
+ for (k = nalloc; k < pid + incr; k++) {
+ pst[k].cells = 0;
+ pst[k].edges = 0;
+ }
+
+ nalloc = pid + incr;
+ }
+
+ pst[pid].cells = 1;
+ pst[pid].type.t = CELL_TYPE;
+ pst[pid].type.val.c = corrCell;
+ }
+ /* update edge count for corr */
+ if (Rast_is_c_null_value(&precCell) || precCell != corrCell)
+ pst[pid_corr[j]].edges++;
+ if (Rast_is_c_null_value(&supCell) || supCell != corrCell)
+ pst[pid_corr[j]].edges++;
+ if (i == ad->rl - 1)
+ pst[pid_corr[j]].edges++;
+ if (j == ad->cl - 1)
+ pst[pid_corr[j]].edges++;
+ /* update edge count for prec */
+ if (!Rast_is_c_null_value(&precCell) && precCell != corrCell)
+ pst[pid_corr[j - 1]].edges++;
+ /* update edge count for sup */
+ if (!Rast_is_c_null_value(&supCell) && supCell != corrCell)
+ pst[pid_sup[j]].edges++;
+
+ precCell = corrCell;
+ }
+ }
+
+ if (npatch > 0) {
+ double edges, cells;
+
+ edges = cells = 0;
+
+ for (i = 1; i <= pid; i++) {
+ cells += pst[i].cells;
+ edges += pst[i].edges;
+ }
+
+ *result = 0.25 * edges / sqrt(cells);
+ }
+ else {
+ Rast_set_d_null_value(result, 1);
+ }
+
+ if (masked) {
+ close(mask_fd);
+ G_free(mask_buf);
+ G_free(mask_sup);
+ }
+ G_free(buf_null);
+ G_free(pid_corr);
+ G_free(pid_sup);
+ G_free(pst);
+
+ return RLI_OK;
+}
+
+
+int calculateD(int fd, struct area_entry *ad, double *result)
+{
+ DCELL *buf, *buf_sup, *buf_null;
+ DCELL corrCell, precCell, supCell;
+ long npatch, area;
+ long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
+ struct pst *pst;
+ long nalloc, incr;
+ int i, j, k;
+ int connected;
+ int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
+
+ buf_null = Rast_allocate_d_buf();
+ Rast_set_d_null_value(buf_null, Rast_window_cols());
+ buf_sup = buf_null;
+
+ /* initialize patch ids */
+ pid_corr = G_malloc(ad->cl * sizeof(long));
+ pid_sup = G_malloc(ad->cl * sizeof(long));
+
+ for (j = 0; j < ad->cl; j++) {
+ pid_corr[j] = 0;
+ pid_sup[j] = 0;
+ }
+
+ /* open mask if needed */
+ 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)
+ return RLI_ERRORE;
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_buf failed");
+ return RLI_ERRORE;
+ }
+ for (j = 0; j < ad->cl; j++)
+ mask_buf[j] = 0;
- for (i = 0; i < ad->rl; i++) {
+ masked = TRUE;
+ }
+
+ /* calculate number of patches */
+ npatch = 0;
+ area = 0;
+ pid = 0;
+
+ /* patch size and type */
+ incr = 1024;
+ if (incr > ad->rl)
+ incr = ad->rl;
+ if (incr > ad->cl)
+ incr = ad->cl;
+ if (incr < 2)
+ incr = 2;
+ nalloc = incr;
+ pst = G_malloc(nalloc * sizeof(struct pst));
+ for (k = 0; k < nalloc; k++) {
+ pst[k].cells = 0;
+ pst[k].edges = 0;
+ }
+
+ for (i = 0; i < ad->rl; i++) {
+ buf = RLI_get_dcell_raster_row(fd, i + ad->y, ad);
+ if (i > 0) {
+ buf_sup = RLI_get_dcell_raster_row(fd, i - 1 + ad->y, ad);
+ }
+
+ if (masked) {
+ mask_tmp = mask_sup;
+ mask_sup = mask_buf;
+ mask_buf = mask_tmp;
if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
return 0;
+ }
- for (k = 0; k < ad->cl; k++) {
- if (mask_buf[k] == 0) {
- null_count++;
+ ltmp = pid_sup;
+ pid_sup = pid_corr;
+ pid_corr = ltmp;
+
+ Rast_set_d_null_value(&precCell, 1);
+
+ connected = 0;
+ for (j = 0; j < ad->cl; j++) {
+ pid_corr[j] = 0;
+
+ corrCell = buf[j + ad->x];
+ if (masked && (mask_buf[j] == 0)) {
+ Rast_set_d_null_value(&corrCell, 1);
+ }
+ else {
+ /* total sample area */
+ area++;
+ }
+
+ supCell = buf_sup[j + ad->x];
+ if (masked && (mask_sup[j] == 0)) {
+ Rast_set_d_null_value(&supCell, 1);
+ }
+
+ if (Rast_is_d_null_value(&corrCell)) {
+ if (!Rast_is_d_null_value(&precCell))
+ pst[pid_corr[j - 1]].edges++;
+ if (!Rast_is_d_null_value(&supCell))
+ pst[pid_sup[j]].edges++;
+ connected = 0;
+ precCell = corrCell;
+ continue;
+ }
+
+ if (!Rast_is_d_null_value(&precCell) && corrCell == precCell) {
+ pid_corr[j] = pid_corr[j - 1];
+ connected = 1;
+ pst[pid_corr[j]].cells++;
+ }
+ else {
+ connected = 0;
+ }
+
+ if (!Rast_is_d_null_value(&supCell) && corrCell == supCell) {
+
+ if (pid_corr[j] != pid_sup[j]) {
+ /* connect or merge */
+ /* after r.clump */
+ if (connected) {
+ npatch--;
+
+ if (npatch == 0) {
+ G_fatal_error("npatch == 0 at row %d, col %d", i, j);
+ }
+ }
+
+ old_pid = pid_corr[j];
+ new_pid = pid_sup[j];
+ pid_corr[j] = new_pid;
+ if (old_pid > 0) {
+ /* merge */
+ /* update left side of the current row */
+ for (k = 0; k < j; k++) {
+ if (pid_corr[k] == old_pid)
+ pid_corr[k] = new_pid;
+ }
+ /* update right side of the previous row */
+ for (k = j + 1; k < ad->cl; k++) {
+ if (pid_sup[k] == old_pid)
+ pid_sup[k] = new_pid;
+ }
+ pst[new_pid].cells += pst[old_pid].cells;
+ pst[old_pid].cells = 0;
+ pst[new_pid].edges += pst[old_pid].edges;
+ pst[old_pid].edges = 0;
+
+ if (old_pid == pid)
+ pid--;
+ }
+ else {
+ pst[new_pid].cells++;
+ }
}
+ connected = 1;
}
+
+ if (!connected) {
+ /* start new patch */
+ npatch++;
+ pid++;
+ pid_corr[j] = pid;
+
+ if (pid >= nalloc) {
+ pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
+
+ for (k = nalloc; k < pid + incr; k++) {
+ pst[k].cells = 0;
+ pst[k].edges = 0;
+ }
+
+ nalloc = pid + incr;
+ }
+
+ pst[pid].cells = 1;
+ pst[pid].type.t = CELL_TYPE;
+ pst[pid].type.val.c = corrCell;
+ }
+ /* update edge count for corr */
+ if (Rast_is_d_null_value(&precCell) || precCell != corrCell)
+ pst[pid_corr[j]].edges++;
+ if (Rast_is_d_null_value(&supCell) || supCell != corrCell)
+ pst[pid_corr[j]].edges++;
+ if (i == ad->rl - 1)
+ pst[pid_corr[j]].edges++;
+ if (j == ad->cl - 1)
+ pst[pid_corr[j]].edges++;
+ /* update edge count for prec */
+ if (!Rast_is_d_null_value(&precCell) && precCell != corrCell)
+ pst[pid_corr[j - 1]].edges++;
+ /* update edge count for sup */
+ if (!Rast_is_d_null_value(&supCell) && supCell != corrCell)
+ pst[pid_sup[j]].edges++;
+
+ precCell = corrCell;
}
+ }
+
+ if (npatch > 0) {
+ double edges, cells;
+
+ edges = cells = 0;
+
+ for (i = 1; i <= pid; i++) {
+ cells += pst[i].cells;
+ edges += pst[i].edges;
+ }
+
+ *result = 0.25 * edges / sqrt(cells);
+ }
+ else {
+ Rast_set_d_null_value(result, 1);
+ }
+
+ if (masked) {
+ close(mask_fd);
G_free(mask_buf);
+ G_free(mask_sup);
}
+ G_free(buf_null);
+ G_free(pid_corr);
+ G_free(pid_sup);
+ G_free(pst);
- /* 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);
+ return RLI_OK;
+}
- area = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
- (((NS_DIST1 + NS_DIST2) / 2) / hd.rows) *
- (ad->rl * ad->cl - null_count);
+int calculateF(int fd, struct area_entry *ad, double *result)
+{
+ FCELL *buf, *buf_sup, *buf_null;
+ FCELL corrCell, precCell, supCell;
+ long npatch, area;
+ long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
+ struct pst *pst;
+ long nalloc, incr;
+ int i, j, k;
+ int connected;
+ int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
- *result = area;
+ buf_null = Rast_allocate_f_buf();
+ Rast_set_f_null_value(buf_null, Rast_window_cols());
+ buf_sup = buf_null;
- return 1;
+ /* initialize patch ids */
+ pid_corr = G_malloc(ad->cl * sizeof(long));
+ pid_sup = G_malloc(ad->cl * sizeof(long));
+
+ for (j = 0; j < ad->cl; j++) {
+ pid_corr[j] = 0;
+ pid_sup[j] = 0;
+ }
+
+ /* open mask if needed */
+ 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)
+ return RLI_ERRORE;
+ 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_buf failed");
+ return RLI_ERRORE;
+ }
+ for (j = 0; j < ad->cl; j++)
+ mask_buf[j] = 0;
+
+ masked = TRUE;
+ }
+
+ /* calculate number of patches */
+ npatch = 0;
+ area = 0;
+ pid = 0;
+
+ /* patch size and type */
+ incr = 1024;
+ if (incr > ad->rl)
+ incr = ad->rl;
+ if (incr > ad->cl)
+ incr = ad->cl;
+ if (incr < 2)
+ incr = 2;
+ nalloc = incr;
+ pst = G_malloc(nalloc * sizeof(struct pst));
+ for (k = 0; k < nalloc; k++) {
+ pst[k].cells = 0;
+ pst[k].edges = 0;
+ }
+
+ for (i = 0; i < ad->rl; i++) {
+ buf = RLI_get_fcell_raster_row(fd, i + ad->y, ad);
+ if (i > 0) {
+ buf_sup = RLI_get_fcell_raster_row(fd, i - 1 + ad->y, ad);
+ }
+
+ if (masked) {
+ mask_tmp = mask_sup;
+ mask_sup = mask_buf;
+ mask_buf = mask_tmp;
+ if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
+ return 0;
+ }
+
+ ltmp = pid_sup;
+ pid_sup = pid_corr;
+ pid_corr = ltmp;
+
+ Rast_set_f_null_value(&precCell, 1);
+
+ connected = 0;
+ for (j = 0; j < ad->cl; j++) {
+ pid_corr[j] = 0;
+
+ corrCell = buf[j + ad->x];
+ if (masked && (mask_buf[j] == 0)) {
+ Rast_set_f_null_value(&corrCell, 1);
+ }
+ else {
+ /* total sample area */
+ area++;
+ }
+
+ supCell = buf_sup[j + ad->x];
+ if (masked && (mask_sup[j] == 0)) {
+ Rast_set_f_null_value(&supCell, 1);
+ }
+
+ if (Rast_is_f_null_value(&corrCell)) {
+ if (!Rast_is_f_null_value(&precCell))
+ pst[pid_corr[j - 1]].edges++;
+ if (!Rast_is_f_null_value(&supCell))
+ pst[pid_sup[j]].edges++;
+ connected = 0;
+ precCell = corrCell;
+ continue;
+ }
+
+ if (!Rast_is_f_null_value(&precCell) && corrCell == precCell) {
+ pid_corr[j] = pid_corr[j - 1];
+ connected = 1;
+ pst[pid_corr[j]].cells++;
+ }
+ else {
+ connected = 0;
+ }
+
+ if (!Rast_is_f_null_value(&supCell) && corrCell == supCell) {
+
+ if (pid_corr[j] != pid_sup[j]) {
+ /* connect or merge */
+ /* after r.clump */
+ if (connected) {
+ npatch--;
+
+ if (npatch == 0) {
+ G_fatal_error("npatch == 0 at row %d, col %d", i, j);
+ }
+ }
+
+ old_pid = pid_corr[j];
+ new_pid = pid_sup[j];
+ pid_corr[j] = new_pid;
+ if (old_pid > 0) {
+ /* merge */
+ /* update left side of the current row */
+ for (k = 0; k < j; k++) {
+ if (pid_corr[k] == old_pid)
+ pid_corr[k] = new_pid;
+ }
+ /* update right side of the previous row */
+ for (k = j + 1; k < ad->cl; k++) {
+ if (pid_sup[k] == old_pid)
+ pid_sup[k] = new_pid;
+ }
+ pst[new_pid].cells += pst[old_pid].cells;
+ pst[old_pid].cells = 0;
+ pst[new_pid].edges += pst[old_pid].edges;
+ pst[old_pid].edges = 0;
+
+ if (old_pid == pid)
+ pid--;
+ }
+ else {
+ pst[new_pid].cells++;
+ }
+ }
+ connected = 1;
+ }
+
+ if (!connected) {
+ /* start new patch */
+ npatch++;
+ pid++;
+ pid_corr[j] = pid;
+
+ if (pid >= nalloc) {
+ pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
+
+ for (k = nalloc; k < pid + incr; k++) {
+ pst[k].cells = 0;
+ pst[k].edges = 0;
+ }
+
+ nalloc = pid + incr;
+ }
+
+ pst[pid].cells = 1;
+ pst[pid].type.t = CELL_TYPE;
+ pst[pid].type.val.c = corrCell;
+ }
+ /* update edge count for corr */
+ if (Rast_is_f_null_value(&precCell) || precCell != corrCell)
+ pst[pid_corr[j]].edges++;
+ if (Rast_is_f_null_value(&supCell) || supCell != corrCell)
+ pst[pid_corr[j]].edges++;
+ if (i == ad->rl - 1)
+ pst[pid_corr[j]].edges++;
+ if (j == ad->cl - 1)
+ pst[pid_corr[j]].edges++;
+ /* update edge count for prec */
+ if (!Rast_is_f_null_value(&precCell) && precCell != corrCell)
+ pst[pid_corr[j - 1]].edges++;
+ /* update edge count for sup */
+ if (!Rast_is_f_null_value(&supCell) && supCell != corrCell)
+ pst[pid_sup[j]].edges++;
+
+ precCell = corrCell;
+ }
+ }
+
+ if (npatch > 0) {
+ double edges, cells;
+
+ edges = cells = 0;
+
+ for (i = 1; i <= pid; i++) {
+ cells += pst[i].cells;
+ edges += pst[i].edges;
+ }
+
+ *result = 0.25 * edges / sqrt(cells);
+ }
+ else {
+ Rast_set_d_null_value(result, 1);
+ }
+
+ if (masked) {
+ close(mask_fd);
+ G_free(mask_buf);
+ G_free(mask_sup);
+ }
+ G_free(buf_null);
+ G_free(pid_corr);
+ G_free(pid_sup);
+ G_free(pst);
+
+ return RLI_OK;
}
Modified: grass/trunk/raster/r.li/r.li.shape/r.li.shape.html
===================================================================
--- grass/trunk/raster/r.li/r.li.shape/r.li.shape.html 2014-03-02 21:03:37 UTC (rev 59168)
+++ grass/trunk/raster/r.li/r.li.shape/r.li.shape.html 2014-03-02 22:34:11 UTC (rev 59169)
@@ -1,8 +1,16 @@
<h2>DESCRIPTION</h2>
-<em>r.li.shape</em> calculates the "shape index" as:<br>
-<I> f(sample_area)= Area </I><br>
+<em>r.li.shape</em> calculates the landscape shape index as:
+<div class="code"><pre>
+LSI = 0.25 * E / sqrt(A)
+</pre></div>
+with:
+<ul>
+ <li><b>E</b>: sum of all edges</li>
+ <li><b>A</b>: sum of all patch areas</li>
+</ul>
+
<h2>NOTES</h2>
Do not use absolute path names for the <b>output</b> map/file.
@@ -12,13 +20,7 @@
<p>
<!-- TODO: verify next: -->
If the input raster map contains only NULL values then <em>r.li.shape</em>
-returns -1.<br>
-<!-- does that mean the program exit code or raster values?? -->
-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.
+returns NULL.
<h2>EXAMPLES</h2>
To calculate the shape index on map <em>my_map</em>, using
More information about the grass-commit
mailing list