[GRASS-SVN] r59022 - grass/trunk/raster/r.li/r.li.patchdensity
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Feb 13 03:38:00 PST 2014
Author: mmetz
Date: 2014-02-13 03:38:00 -0800 (Thu, 13 Feb 2014)
New Revision: 59022
Modified:
grass/trunk/raster/r.li/r.li.patchdensity/main.c
grass/trunk/raster/r.li/r.li.patchdensity/r.li.patchdensity.html
Log:
r.li.patchdensity: fix number of patches
Modified: grass/trunk/raster/r.li/r.li.patchdensity/main.c
===================================================================
--- grass/trunk/raster/r.li/r.li.patchdensity/main.c 2014-02-13 11:31:59 UTC (rev 59021)
+++ grass/trunk/raster/r.li/r.li.patchdensity/main.c 2014-02-13 11:38:00 UTC (rev 59022)
@@ -59,30 +59,36 @@
int patch_density(int fd, char **par, struct area_entry *ad, double *result)
{
- CELL *buf, *sup;
- int count = 0, i, j, connected = 0, complete_line = 1, other_above = 0;
+ CELL *buf, *sup, *cnull;
+ CELL pid, old_pid, *pid_curr, *pid_sup, *ctmp;
+ int count, i, j, k, connected, other_above;
+ struct Cell_head hd;
+ int mask_fd, *mask_buf, null_count;
double area;
- struct Cell_head hd;
- CELL complete_value;
double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
- int mask_fd = -1, *mask_buf, *mask_sup, null_count = 0;
- G_debug(1, "begin patch_density() index");
-
- Rast_set_c_null_value(&complete_value, 1);
Rast_get_cellhd(ad->raster, "", &hd);
+ cnull = Rast_allocate_c_buf();
+ Rast_set_c_null_value(cnull, Rast_window_cols());
+ sup = cnull;
+
+ /* initialize patch ids */
+ pid_curr = Rast_allocate_c_buf();
+ Rast_set_c_null_value(pid_curr, Rast_window_cols());
+ pid_sup = Rast_allocate_c_buf();
+ Rast_set_c_null_value(pid_sup, Rast_window_cols());
+
/* open mask if needed */
+ mask_fd = -1;
+ mask_buf = NULL;
+ null_count = 0;
if (ad->mask == 1) {
if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
return 0;
-
- mask_buf = malloc(ad->cl * sizeof(int));
- mask_sup = malloc(ad->cl * sizeof(int));
+ mask_buf = G_malloc(ad->cl * sizeof(int));
}
- sup = Rast_allocate_c_buf();
-
/* calculate distance */
G_begin_distance_calculations();
/* EW Dist at North edge */
@@ -94,103 +100,96 @@
/* NS Dist at West edge */
NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
+ /* calculate number of patches */
+ count = 0;
+ connected = 0;
+ other_above = 0;
+ pid = 0;
- /* calculate number of patch */
-
for (i = 0; i < ad->rl; i++) {
buf = RLI_get_cell_raster_row(fd, i + ad->y, ad);
-
if (i > 0) {
sup = RLI_get_cell_raster_row(fd, i - 1 + ad->y, ad);
}
-
/* mask values */
if (ad->mask == 1) {
- int k;
-
- if (i > 0) {
- int *tmp;
-
- tmp = mask_sup;
- mask_buf = mask_sup;
- }
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) {
- Rast_set_c_null_value(mask_buf + k, 1);
+ Rast_set_c_null_value(&buf[k + ad->x], 1);
null_count++;
}
}
-
}
+
+ ctmp = pid_sup;
+ pid_sup = pid_curr;
+ pid_curr = ctmp;
+ Rast_set_c_null_value(pid_curr, Rast_window_cols());
+ connected = 0;
+ other_above = 1;
+ for (j = 0; j < ad->cl; j++) {
+
+ if (Rast_is_null_value(&(buf[j + ad->x]), CELL_TYPE)) {
+ connected = 0;
+ other_above = 1;
+ continue;
+ }
- if (complete_line) {
- if (!Rast_is_null_value(&(buf[ad->x]), CELL_TYPE) &&
- buf[ad->x] != complete_value)
- count++;
-
- for (j = 0; j < ad->cl - 1; j++) {
-
- if (buf[j + ad->x] != buf[j + 1 + ad->x]) {
- complete_line = 0;
- if (!Rast_is_null_value(&(buf[j + 1 + ad->x]), CELL_TYPE)
- && buf[j + 1 + ad->x] != complete_value)
- count++;
+ if (sup[j + ad->x] == buf[j + ad->x]) {
+
+ if (!connected) {
+ pid_curr[j + ad->x] = pid_sup[j + ad->x];
}
-
- }
- if (complete_line) {
- complete_value = buf[ad->x];
- }
- }
- else {
- complete_line = 1;
- connected = 0;
- other_above = 0;
- for (j = 0; j < ad->cl; j++) {
- if (sup[j + ad->x] == buf[j + ad->x]) {
- connected = 1;
- if (other_above) {
- other_above = 0;
+
+ if (pid_curr[j + ad->x] != pid_sup[j + ad->x]) {
+ if (connected) {
count--;
}
- }
- else {
- if (connected &&
- !Rast_is_null_value(&(buf[j + ad->x]), CELL_TYPE))
- other_above = 1;
- }
- if (j < ad->cl - 1 && buf[j + ad->x] != buf[j + 1 + ad->x]) {
- complete_line = 0;
- if (!connected &&
- !Rast_is_null_value(&(buf[j + ad->x]), CELL_TYPE)) {
-
- count++;
- connected = 0;
- other_above = 0;
+ if (other_above) {
+ pid_curr[j + ad->x] = pid_sup[j + ad->x];
+ for (k = j + ad->x - 1; k >= ad->x; k--) {
+ if (buf[k] != buf[j + ad->x])
+ break;
+ pid_curr[k] = pid_sup[j + ad->x];
+ }
}
else {
- connected = 0;
- other_above = 0;
+ old_pid = pid_sup[j + ad->x];
+ pid_sup[j + ad->x] = pid_curr[j + ad->x];
+
+ for (k = j + 1; k < ad->cl; k++) {
+ if (pid_sup[k + ad->x] == old_pid) {
+ pid_sup[k + ad->x] = pid_curr[j + ad->x];
+ }
+ }
}
}
+
+ other_above = 0;
+ connected = 1;
}
- if (!connected &&
- sup[ad->cl - 1 + ad->x] != buf[ad->cl - 1 + ad->x]) {
- if (!Rast_is_null_value
- (&(buf[ad->cl - 1 + ad->x]), CELL_TYPE)) {
- count++;
- complete_line = 0;
+
+ if (!connected) {
+ count++;
+ pid++;
+ pid_curr[j + ad->x] = pid;
+ }
+
+ if (j < ad->cl - 1) {
+ if (buf[j + ad->x] == buf[j + 1 + ad->x]) {
+
+ connected = 1;
+ pid_curr[j + 1 + ad->x] = pid_curr[j + ad->x];
}
+ else {
+ other_above = 1;
+ connected = 0;
+ }
}
-
- if (complete_line)
- complete_value = buf[ad->x];
-
}
-
}
area = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
@@ -200,9 +199,15 @@
if (area != 0)
*result = (count / area) * 1000000;
else
- *result = -1;
+ Rast_set_d_null_value(result, 1);
- G_free(sup);
+ if (ad->mask == 1) {
+ G_free(mask_buf);
+ }
+ G_free(cnull);
+ G_free(pid_curr);
+ G_free(pid_sup);
+
return RLI_OK;
}
Modified: grass/trunk/raster/r.li/r.li.patchdensity/r.li.patchdensity.html
===================================================================
--- grass/trunk/raster/r.li/r.li.patchdensity/r.li.patchdensity.html 2014-02-13 11:31:59 UTC (rev 59021)
+++ grass/trunk/raster/r.li/r.li.patchdensity/r.li.patchdensity.html 2014-02-13 11:38:00 UTC (rev 59022)
@@ -5,8 +5,10 @@
f(sample_area) = (Patch_Number/Area) * 1000000
</pre></div>
-that is 1000000 by number of patch for area unit.
-This index is calculated using a 4 neighbour algorithm.
+that is the number of patches per 1,000,000 area units.
+<p>
+This index is calculated using a 4 neighbour algorithm, diagonal cells
+are ignored when tracing a patch.
<h2>NOTES</h2>
@@ -17,11 +19,12 @@
<p>
<!-- TODO: verify next: -->
A map of NULL values is considered to have zero patches. <br>
-If raster area is 0, <em>r.li.patchdensity</em> returns -1. This is only
-possible if the raster is masked.<br>
-If you want to change these -1 values to NULL, run subsequently on the resulting map:
+If the size of the sample area is 0, <em>r.li.patchdensity</em> returns
+NULL. This is only possible if the raster is masked.<br>
+If you want to change these NULL values to some number, run subsequently
+on the resulting map:
<div class="code"><pre>
-r.null setnull=-1 input=my_out
+r.null null=-1 input=my_out
</pre></div>
after index calculation.
More information about the grass-commit
mailing list