padsd.c 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832
  1. /****************************************************************************
  2. *
  3. * MODULE: r.li.padsd
  4. * AUTHOR(S): Serena Pallecchi (original contributor)
  5. * student of Computer Science University of Pisa (Italy)
  6. * Commission from Faunalia Pontedera (PI) www.faunalia.it
  7. * Fixes: Markus Neteler <neteler itc.it>
  8. * Rewrite: Markus Metz
  9. * Patch identification: Michael Shapiro - CERL
  10. *
  11. * PURPOSE: calculates calculates standard deviation of patch areas
  12. * COPYRIGHT: (C) 2007-2014 by the GRASS Development Team
  13. *
  14. * This program is free software under the GNU General Public
  15. * License (>=v2). Read the file COPYING that comes with GRASS
  16. * for details.
  17. *
  18. *****************************************************************************/
  19. #include <stdlib.h>
  20. #include <fcntl.h>
  21. #include <math.h>
  22. #include <grass/gis.h>
  23. #include <grass/raster.h>
  24. #include <grass/glocale.h>
  25. #include "../r.li.daemon/daemon.h"
  26. #include "../r.li.daemon/GenericCell.h"
  27. /* template is patchnum */
  28. /* cell count and type of each patch */
  29. struct pst {
  30. long count;
  31. generic_cell type;
  32. };
  33. rli_func patchAreaDistributionSD;
  34. int calculate(int fd, struct area_entry *ad, double *result);
  35. int calculateD(int fd, struct area_entry *ad, double *result);
  36. int calculateF(int fd, struct area_entry *ad, double *result);
  37. int main(int argc, char *argv[])
  38. {
  39. struct Option *raster, *conf, *output;
  40. struct GModule *module;
  41. G_gisinit(argv[0]);
  42. module = G_define_module();
  43. module->description =
  44. _("Calculates standard deviation of patch area a raster map");
  45. G_add_keyword(_("raster"));
  46. G_add_keyword(_("landscape structure analysis"));
  47. G_add_keyword(_("patch index"));
  48. /* define options */
  49. raster = G_define_standard_option(G_OPT_R_INPUT);
  50. conf = G_define_standard_option(G_OPT_F_INPUT);
  51. conf->key = "config";
  52. conf->description = _("Configuration file");
  53. conf->required = YES;
  54. output = G_define_standard_option(G_OPT_R_OUTPUT);
  55. if (G_parser(argc, argv))
  56. exit(EXIT_FAILURE);
  57. return calculateIndex(conf->answer, patchAreaDistributionSD, NULL,
  58. raster->answer, output->answer);
  59. }
  60. int patchAreaDistributionSD(int fd, char **par, struct area_entry *ad, double *result)
  61. {
  62. double indice = 0;
  63. int ris = RLI_OK;
  64. switch (ad->data_type) {
  65. case CELL_TYPE:
  66. {
  67. ris = calculate(fd, ad, &indice);
  68. break;
  69. }
  70. case DCELL_TYPE:
  71. {
  72. ris = calculateD(fd, ad, &indice);
  73. break;
  74. }
  75. case FCELL_TYPE:
  76. {
  77. ris = calculateF(fd, ad, &indice);
  78. break;
  79. }
  80. default:
  81. {
  82. G_fatal_error("data type unknown");
  83. return RLI_ERRORE;
  84. }
  85. }
  86. if (ris != RLI_OK) {
  87. *result = -1;
  88. return RLI_ERRORE;
  89. }
  90. *result = indice;
  91. return RLI_OK;
  92. }
  93. int calculate(int fd, struct area_entry *ad, double *result)
  94. {
  95. CELL *buf, *buf_sup, *buf_null;
  96. CELL corrCell, precCell, supCell;
  97. long npatch, area;
  98. long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
  99. struct pst *pst;
  100. long nalloc, incr;
  101. int i, j, k;
  102. int connected;
  103. int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
  104. struct Cell_head hd;
  105. Rast_get_window(&hd);
  106. buf_null = Rast_allocate_c_buf();
  107. Rast_set_c_null_value(buf_null, Rast_window_cols());
  108. buf_sup = buf_null;
  109. /* initialize patch ids */
  110. pid_corr = G_malloc(Rast_window_cols() * sizeof(long));
  111. pid_sup = G_malloc(Rast_window_cols() * sizeof(long));
  112. for (j = 0; j < Rast_window_cols(); j++) {
  113. pid_corr[j] = 0;
  114. pid_sup[j] = 0;
  115. }
  116. /* open mask if needed */
  117. mask_fd = -1;
  118. mask_buf = mask_sup = NULL;
  119. masked = FALSE;
  120. if (ad->mask == 1) {
  121. if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
  122. return RLI_ERRORE;
  123. mask_buf = G_malloc(ad->cl * sizeof(int));
  124. if (mask_buf == NULL) {
  125. G_fatal_error("malloc mask_buf failed");
  126. return RLI_ERRORE;
  127. }
  128. mask_sup = G_malloc(ad->cl * sizeof(int));
  129. if (mask_sup == NULL) {
  130. G_fatal_error("malloc mask_buf failed");
  131. return RLI_ERRORE;
  132. }
  133. for (j = 0; j < ad->cl; j++)
  134. mask_buf[j] = 0;
  135. masked = TRUE;
  136. }
  137. /* calculate number of patches */
  138. npatch = 0;
  139. area = 0;
  140. pid = 0;
  141. /* patch size and type */
  142. incr = 1024;
  143. if (incr > ad->rl)
  144. incr = ad->rl;
  145. if (incr > ad->cl)
  146. incr = ad->cl;
  147. if (incr < 2)
  148. incr = 2;
  149. nalloc = incr;
  150. pst = G_malloc(nalloc * sizeof(struct pst));
  151. for (k = 0; k < nalloc; k++) {
  152. pst[k].count = 0;
  153. }
  154. for (i = 0; i < ad->rl; i++) {
  155. buf = RLI_get_cell_raster_row(fd, i + ad->y, ad);
  156. if (i > 0) {
  157. buf_sup = RLI_get_cell_raster_row(fd, i - 1 + ad->y, ad);
  158. }
  159. if (masked) {
  160. mask_tmp = mask_sup;
  161. mask_sup = mask_buf;
  162. mask_buf = mask_tmp;
  163. if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
  164. return 0;
  165. }
  166. ltmp = pid_sup;
  167. pid_sup = pid_corr;
  168. pid_corr = ltmp;
  169. Rast_set_c_null_value(&precCell, 1);
  170. connected = 0;
  171. for (j = 0; j < ad->cl; j++) {
  172. pid_corr[j + ad->x] = 0;
  173. corrCell = buf[j + ad->x];
  174. if (masked && (mask_buf[j] == 0)) {
  175. Rast_set_c_null_value(&corrCell, 1);
  176. }
  177. if (Rast_is_c_null_value(&corrCell)) {
  178. connected = 0;
  179. precCell = corrCell;
  180. continue;
  181. }
  182. area++;
  183. supCell = buf_sup[j + ad->x];
  184. if (masked && (mask_sup[j] == 0)) {
  185. Rast_set_c_null_value(&supCell, 1);
  186. }
  187. if (!Rast_is_c_null_value(&precCell) && corrCell == precCell) {
  188. pid_corr[j + ad->x] = pid_corr[j - 1 + ad->x];
  189. connected = 1;
  190. pst[pid_corr[j + ad->x]].count++;
  191. }
  192. else {
  193. connected = 0;
  194. }
  195. if (!Rast_is_c_null_value(&supCell) && corrCell == supCell) {
  196. if (pid_corr[j + ad->x] != pid_sup[j + ad->x]) {
  197. /* connect or merge */
  198. /* after r.clump */
  199. if (connected) {
  200. npatch--;
  201. if (npatch == 0) {
  202. G_fatal_error("npatch == 0 at row %d, col %d", i, j);
  203. }
  204. }
  205. old_pid = pid_corr[j + ad->x];
  206. new_pid = pid_sup[j + ad->x];
  207. pid_corr[j + ad->x] = new_pid;
  208. if (old_pid > 0) {
  209. /* merge */
  210. /* update left side of the current row */
  211. for (k = 0; k < j; k++) {
  212. if (pid_corr[k + ad->x] == old_pid)
  213. pid_corr[k + ad->x] = new_pid;
  214. }
  215. /* update right side of the previous row */
  216. for (k = j + 1; k < ad->cl; k++) {
  217. if (pid_sup[k + ad->x] == old_pid)
  218. pid_sup[k + ad->x] = new_pid;
  219. }
  220. pst[new_pid].count += pst[old_pid].count;
  221. pst[old_pid].count = 0;
  222. if (old_pid == pid)
  223. pid--;
  224. }
  225. else {
  226. pst[new_pid].count++;
  227. }
  228. }
  229. connected = 1;
  230. }
  231. if (!connected) {
  232. /* start new patch */
  233. npatch++;
  234. pid++;
  235. pid_corr[j + ad->x] = pid;
  236. if (pid >= nalloc) {
  237. pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
  238. for (k = nalloc; k < pid + incr; k++)
  239. pst[k].count = 0;
  240. nalloc = pid + incr;
  241. }
  242. pst[pid].count = 1;
  243. pst[pid].type.t = CELL_TYPE;
  244. pst[pid].type.val.c = corrCell;
  245. }
  246. precCell = corrCell;
  247. }
  248. }
  249. if (npatch > 0) {
  250. double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
  251. double area_units;
  252. double area_p;
  253. double cell_size_m;
  254. double mps;
  255. double psd;
  256. double diff;
  257. /* calculate distance */
  258. G_begin_distance_calculations();
  259. /* EW Dist at North edge */
  260. EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
  261. /* EW Dist at South Edge */
  262. EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
  263. /* NS Dist at East edge */
  264. NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
  265. /* NS Dist at West edge */
  266. NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
  267. cell_size_m = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
  268. (((NS_DIST1 + NS_DIST2) / 2) / hd.rows);
  269. area_units = cell_size_m * area;
  270. /* mean patch size in hectares */
  271. mps = area_units / (npatch * 10000);
  272. /* calculate patch size standard deviation */
  273. psd = 0;
  274. for (old_pid = 1; old_pid <= pid; old_pid++) {
  275. if (pst[old_pid].count > 0) {
  276. area_p = cell_size_m * pst[old_pid].count / 10000;
  277. diff = area_p - mps;
  278. psd += diff * diff;
  279. }
  280. }
  281. psd = sqrt(psd / npatch);
  282. *result = psd;
  283. }
  284. else
  285. Rast_set_d_null_value(result, 1);
  286. if (masked) {
  287. close(mask_fd);
  288. G_free(mask_buf);
  289. G_free(mask_sup);
  290. }
  291. G_free(buf_null);
  292. G_free(pid_corr);
  293. G_free(pid_sup);
  294. G_free(pst);
  295. return RLI_OK;
  296. }
  297. int calculateD(int fd, struct area_entry *ad, double *result)
  298. {
  299. DCELL *buf, *buf_sup, *buf_null;
  300. DCELL corrCell, precCell, supCell;
  301. long npatch, area;
  302. long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
  303. struct pst *pst;
  304. long nalloc, incr;
  305. int i, j, k;
  306. int connected;
  307. int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
  308. struct Cell_head hd;
  309. Rast_get_window(&hd);
  310. buf_null = Rast_allocate_d_buf();
  311. Rast_set_d_null_value(buf_null, Rast_window_cols());
  312. buf_sup = buf_null;
  313. /* initialize patch ids */
  314. pid_corr = G_malloc(Rast_window_cols() * sizeof(long));
  315. pid_sup = G_malloc(Rast_window_cols() * sizeof(long));
  316. for (j = 0; j < Rast_window_cols(); j++) {
  317. pid_corr[j] = 0;
  318. pid_sup[j] = 0;
  319. }
  320. /* open mask if needed */
  321. mask_fd = -1;
  322. mask_buf = mask_sup = NULL;
  323. masked = FALSE;
  324. if (ad->mask == 1) {
  325. if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
  326. return RLI_ERRORE;
  327. mask_buf = G_malloc(ad->cl * sizeof(int));
  328. if (mask_buf == NULL) {
  329. G_fatal_error("malloc mask_buf failed");
  330. return RLI_ERRORE;
  331. }
  332. mask_sup = G_malloc(ad->cl * sizeof(int));
  333. if (mask_sup == NULL) {
  334. G_fatal_error("malloc mask_buf failed");
  335. return RLI_ERRORE;
  336. }
  337. for (j = 0; j < ad->cl; j++)
  338. mask_buf[j] = 0;
  339. masked = TRUE;
  340. }
  341. /* calculate number of patches */
  342. npatch = 0;
  343. area = 0;
  344. pid = 0;
  345. /* patch size and type */
  346. incr = 1024;
  347. if (incr > ad->rl)
  348. incr = ad->rl;
  349. if (incr > ad->cl)
  350. incr = ad->cl;
  351. if (incr < 2)
  352. incr = 2;
  353. nalloc = incr;
  354. pst = G_malloc(nalloc * sizeof(struct pst));
  355. for (k = 0; k < nalloc; k++) {
  356. pst[k].count = 0;
  357. }
  358. for (i = 0; i < ad->rl; i++) {
  359. buf = RLI_get_dcell_raster_row(fd, i + ad->y, ad);
  360. if (i > 0) {
  361. buf_sup = RLI_get_dcell_raster_row(fd, i - 1 + ad->y, ad);
  362. }
  363. if (masked) {
  364. mask_tmp = mask_sup;
  365. mask_sup = mask_buf;
  366. mask_buf = mask_tmp;
  367. if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
  368. return 0;
  369. }
  370. ltmp = pid_sup;
  371. pid_sup = pid_corr;
  372. pid_corr = ltmp;
  373. Rast_set_d_null_value(&precCell, 1);
  374. connected = 0;
  375. for (j = 0; j < ad->cl; j++) {
  376. pid_corr[j + ad->x] = 0;
  377. corrCell = buf[j + ad->x];
  378. if (masked && (mask_buf[j] == 0)) {
  379. Rast_set_d_null_value(&corrCell, 1);
  380. }
  381. if (Rast_is_d_null_value(&corrCell)) {
  382. connected = 0;
  383. precCell = corrCell;
  384. continue;
  385. }
  386. area++;
  387. supCell = buf_sup[j + ad->x];
  388. if (masked && (mask_sup[j] == 0)) {
  389. Rast_set_d_null_value(&supCell, 1);
  390. }
  391. if (!Rast_is_d_null_value(&precCell) && corrCell == precCell) {
  392. pid_corr[j + ad->x] = pid_corr[j - 1 + ad->x];
  393. connected = 1;
  394. pst[pid_corr[j + ad->x]].count++;
  395. }
  396. else {
  397. connected = 0;
  398. }
  399. if (!Rast_is_d_null_value(&supCell) && corrCell == supCell) {
  400. if (pid_corr[j + ad->x] != pid_sup[j + ad->x]) {
  401. /* connect or merge */
  402. /* after r.clump */
  403. if (connected) {
  404. npatch--;
  405. if (npatch == 0) {
  406. G_fatal_error("npatch == 0 at row %d, col %d", i, j);
  407. }
  408. }
  409. old_pid = pid_corr[j + ad->x];
  410. new_pid = pid_sup[j + ad->x];
  411. pid_corr[j + ad->x] = new_pid;
  412. if (old_pid > 0) {
  413. /* merge */
  414. /* update left side of the current row */
  415. for (k = 0; k < j; k++) {
  416. if (pid_corr[k + ad->x] == old_pid)
  417. pid_corr[k + ad->x] = new_pid;
  418. }
  419. /* update right side of the previous row */
  420. for (k = j + 1; k < ad->cl; k++) {
  421. if (pid_sup[k + ad->x] == old_pid)
  422. pid_sup[k + ad->x] = new_pid;
  423. }
  424. pst[new_pid].count += pst[old_pid].count;
  425. pst[old_pid].count = 0;
  426. if (old_pid == pid)
  427. pid--;
  428. }
  429. else {
  430. pst[new_pid].count++;
  431. }
  432. }
  433. connected = 1;
  434. }
  435. if (!connected) {
  436. /* start new patch */
  437. npatch++;
  438. pid++;
  439. pid_corr[j + ad->x] = pid;
  440. if (pid >= nalloc) {
  441. pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
  442. for (k = nalloc; k < pid + incr; k++)
  443. pst[k].count = 0;
  444. nalloc = pid + incr;
  445. }
  446. pst[pid].count = 1;
  447. pst[pid].type.t = CELL_TYPE;
  448. pst[pid].type.val.c = corrCell;
  449. }
  450. precCell = corrCell;
  451. }
  452. }
  453. if (npatch > 0) {
  454. double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
  455. double area_units;
  456. double area_p;
  457. double cell_size_m;
  458. double mps;
  459. double psd;
  460. double diff;
  461. /* calculate distance */
  462. G_begin_distance_calculations();
  463. /* EW Dist at North edge */
  464. EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
  465. /* EW Dist at South Edge */
  466. EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
  467. /* NS Dist at East edge */
  468. NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
  469. /* NS Dist at West edge */
  470. NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
  471. cell_size_m = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
  472. (((NS_DIST1 + NS_DIST2) / 2) / hd.rows);
  473. area_units = cell_size_m * area;
  474. /* mean patch size in hectares */
  475. mps = area_units / (npatch * 10000);
  476. /* calculate patch size standard deviation */
  477. psd = 0;
  478. for (old_pid = 1; old_pid <= pid; old_pid++) {
  479. if (pst[old_pid].count > 0) {
  480. area_p = cell_size_m * pst[old_pid].count / 10000;
  481. diff = area_p - mps;
  482. psd += diff * diff;
  483. }
  484. }
  485. psd = sqrt(psd / npatch);
  486. *result = psd;
  487. }
  488. else
  489. Rast_set_d_null_value(result, 1);
  490. if (masked) {
  491. close(mask_fd);
  492. G_free(mask_buf);
  493. G_free(mask_sup);
  494. }
  495. G_free(buf_null);
  496. G_free(pid_corr);
  497. G_free(pid_sup);
  498. G_free(pst);
  499. return RLI_OK;
  500. }
  501. int calculateF(int fd, struct area_entry *ad, double *result)
  502. {
  503. FCELL *buf, *buf_sup, *buf_null;
  504. FCELL corrCell, precCell, supCell;
  505. long npatch, area;
  506. long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
  507. struct pst *pst;
  508. long nalloc, incr;
  509. int i, j, k;
  510. int connected;
  511. int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
  512. struct Cell_head hd;
  513. Rast_get_window(&hd);
  514. buf_null = Rast_allocate_f_buf();
  515. Rast_set_f_null_value(buf_null, Rast_window_cols());
  516. buf_sup = buf_null;
  517. /* initialize patch ids */
  518. pid_corr = G_malloc(Rast_window_cols() * sizeof(long));
  519. pid_sup = G_malloc(Rast_window_cols() * sizeof(long));
  520. for (j = 0; j < Rast_window_cols(); j++) {
  521. pid_corr[j] = 0;
  522. pid_sup[j] = 0;
  523. }
  524. /* open mask if needed */
  525. mask_fd = -1;
  526. mask_buf = mask_sup = NULL;
  527. masked = FALSE;
  528. if (ad->mask == 1) {
  529. if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
  530. return RLI_ERRORE;
  531. mask_buf = G_malloc(ad->cl * sizeof(int));
  532. if (mask_buf == NULL) {
  533. G_fatal_error("malloc mask_buf failed");
  534. return RLI_ERRORE;
  535. }
  536. mask_sup = G_malloc(ad->cl * sizeof(int));
  537. if (mask_sup == NULL) {
  538. G_fatal_error("malloc mask_buf failed");
  539. return RLI_ERRORE;
  540. }
  541. for (j = 0; j < ad->cl; j++)
  542. mask_buf[j] = 0;
  543. masked = TRUE;
  544. }
  545. /* calculate number of patches */
  546. npatch = 0;
  547. area = 0;
  548. pid = 0;
  549. /* patch size and type */
  550. incr = 1024;
  551. if (incr > ad->rl)
  552. incr = ad->rl;
  553. if (incr > ad->cl)
  554. incr = ad->cl;
  555. if (incr < 2)
  556. incr = 2;
  557. nalloc = incr;
  558. pst = G_malloc(nalloc * sizeof(struct pst));
  559. for (k = 0; k < nalloc; k++) {
  560. pst[k].count = 0;
  561. }
  562. for (i = 0; i < ad->rl; i++) {
  563. buf = RLI_get_fcell_raster_row(fd, i + ad->y, ad);
  564. if (i > 0) {
  565. buf_sup = RLI_get_fcell_raster_row(fd, i - 1 + ad->y, ad);
  566. }
  567. if (masked) {
  568. mask_tmp = mask_sup;
  569. mask_sup = mask_buf;
  570. mask_buf = mask_tmp;
  571. if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
  572. return 0;
  573. }
  574. ltmp = pid_sup;
  575. pid_sup = pid_corr;
  576. pid_corr = ltmp;
  577. Rast_set_f_null_value(&precCell, 1);
  578. connected = 0;
  579. for (j = 0; j < ad->cl; j++) {
  580. pid_corr[j + ad->x] = 0;
  581. corrCell = buf[j + ad->x];
  582. if (masked && (mask_buf[j] == 0)) {
  583. Rast_set_f_null_value(&corrCell, 1);
  584. }
  585. if (Rast_is_f_null_value(&corrCell)) {
  586. connected = 0;
  587. precCell = corrCell;
  588. continue;
  589. }
  590. area++;
  591. supCell = buf_sup[j + ad->x];
  592. if (masked && (mask_sup[j] == 0)) {
  593. Rast_set_f_null_value(&supCell, 1);
  594. }
  595. if (!Rast_is_f_null_value(&precCell) && corrCell == precCell) {
  596. pid_corr[j + ad->x] = pid_corr[j - 1 + ad->x];
  597. connected = 1;
  598. pst[pid_corr[j + ad->x]].count++;
  599. }
  600. else {
  601. connected = 0;
  602. }
  603. if (!Rast_is_f_null_value(&supCell) && corrCell == supCell) {
  604. if (pid_corr[j + ad->x] != pid_sup[j + ad->x]) {
  605. /* connect or merge */
  606. /* after r.clump */
  607. if (connected) {
  608. npatch--;
  609. if (npatch == 0) {
  610. G_fatal_error("npatch == 0 at row %d, col %d", i, j);
  611. }
  612. }
  613. old_pid = pid_corr[j + ad->x];
  614. new_pid = pid_sup[j + ad->x];
  615. pid_corr[j + ad->x] = new_pid;
  616. if (old_pid > 0) {
  617. /* merge */
  618. /* update left side of the current row */
  619. for (k = 0; k < j; k++) {
  620. if (pid_corr[k + ad->x] == old_pid)
  621. pid_corr[k + ad->x] = new_pid;
  622. }
  623. /* update right side of the previous row */
  624. for (k = j + 1; k < ad->cl; k++) {
  625. if (pid_sup[k + ad->x] == old_pid)
  626. pid_sup[k + ad->x] = new_pid;
  627. }
  628. pst[new_pid].count += pst[old_pid].count;
  629. pst[old_pid].count = 0;
  630. if (old_pid == pid)
  631. pid--;
  632. }
  633. else {
  634. pst[new_pid].count++;
  635. }
  636. }
  637. connected = 1;
  638. }
  639. if (!connected) {
  640. /* start new patch */
  641. npatch++;
  642. pid++;
  643. pid_corr[j + ad->x] = pid;
  644. if (pid >= nalloc) {
  645. pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
  646. for (k = nalloc; k < pid + incr; k++)
  647. pst[k].count = 0;
  648. nalloc = pid + incr;
  649. }
  650. pst[pid].count = 1;
  651. pst[pid].type.t = CELL_TYPE;
  652. pst[pid].type.val.c = corrCell;
  653. }
  654. precCell = corrCell;
  655. }
  656. }
  657. if (npatch > 0) {
  658. double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
  659. double area_units;
  660. double area_p;
  661. double cell_size_m;
  662. double mps;
  663. double psd;
  664. double diff;
  665. /* calculate distance */
  666. G_begin_distance_calculations();
  667. /* EW Dist at North edge */
  668. EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
  669. /* EW Dist at South Edge */
  670. EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
  671. /* NS Dist at East edge */
  672. NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
  673. /* NS Dist at West edge */
  674. NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
  675. cell_size_m = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
  676. (((NS_DIST1 + NS_DIST2) / 2) / hd.rows);
  677. area_units = cell_size_m * area;
  678. /* mean patch size in hectares */
  679. mps = area_units / (npatch * 10000);
  680. /* calculate patch size standard deviation */
  681. psd = 0;
  682. for (old_pid = 1; old_pid <= pid; old_pid++) {
  683. if (pst[old_pid].count > 0) {
  684. area_p = cell_size_m * pst[old_pid].count / 10000;
  685. diff = area_p - mps;
  686. psd += diff * diff;
  687. }
  688. }
  689. psd = sqrt(psd / npatch);
  690. *result = psd;
  691. }
  692. else
  693. Rast_set_d_null_value(result, 1);
  694. if (masked) {
  695. close(mask_fd);
  696. G_free(mask_buf);
  697. G_free(mask_sup);
  698. }
  699. G_free(buf_null);
  700. G_free(pid_corr);
  701. G_free(pid_sup);
  702. G_free(pst);
  703. return RLI_OK;
  704. }