plForestFragment

Initial commit
netzel

netzel commited on 2024-02-21 22:51:55
Showing 3 changed files, with 385 additions and 0 deletions.

Makefile
View file @ 1ed00b7
... ... 
@@ -0,0 +1,28 @@
1 
+
2 
+PROG = plforestfragment
3 
+
4 
+CC = gcc
5 
+
6 
+CFLAGS=-Wall -fopenmp -O3 -msse2 -mfpmath=sse -Wl,-rpath,'$$ORIGIN/' -I/usr/include/gdal -lgomp -lgdal -lezgdal -largtable3 -lm
7 
+
8 
+COMMON_C = $(shell ls *.c)
9 
+COMMON_C_O = $(COMMON_C:%.c=%.o)
10 
+
11 
+ifndef PREFIX
12 
+  PREFIX = /usr/local
13 
+endif
14 
+
15 
+
16 
+default: $(PROG)
17 
+
18 
+$(PROG): $(COMMON_C_O) $(HEADERS)
19 
+	$(CC) $(CFLAGS) -o $(PROG) $(COMMON_O) $(COMMON_C_O)
20 
+
21 
+$(COMMON_C_O): $(COMMON_C) $(HEADERS_C)
22 
+	$(CC) $(CFLAGS) -c  $(COMMON_C)
23 
+
24 
+install:
25 
+	cp -f $(PROG) $(PREFIX)/bin
26 
+
27 
+clean:
28 
+	rm -f *.o $(PROG)
Makefile.win
View file @ 1ed00b7
... ... 
@@ -0,0 +1,20 @@
1 
+PROG = plforestfragment.exe
2 
+
3 
+CC = x86_64-w64-mingw32-gcc
4 
+
5 
+CFLAGS= -I../../Win64_cross/gdal/include -I../../Win64_cross/ezgdal -I../../Win64_cross/argtable3 -L../../Win64_cross/argtable3 -L../../Win64_cross/gdal/bin -L../../Win64_cross/ezgdal -Wall -fopenmp -O3 -msse2 -mfpmath=sse -lgomp -lgdal -llibezgdal -llibargtable3 -lm
6 
+
7 
+COMMON_C = $(shell ls *.c)
8 
+COMMON_C_O = $(COMMON_C:%.c=%.o)
9 
+HEADERS = $(shell ls *.h)
10 
+
11 
+default: $(PROG)
12 
+
13 
+$(PROG): $(COMMON_C_O) $(HEADERS)
14 
+	$(CC) $(CFLAGS) -o $(PROG) $(COMMON_C_O)
15 
+
16 
+$(COMMON_C_O): $(COMMON_C) $(HEADERS)
17 
+	$(CC) $(CFLAGS) -c  $(COMMON_C)
18 
+
19 
+clean:
20 
+	rm -f *.o $(PROG)
main.c
View file @ 1ed00b7
... ... 
@@ -0,0 +1,337 @@
1 
+/****************************************************************************
2 
+ *
3 
+ * PROGRAM:      forestfragment
4 
+ *
5 
+ * AUTHOR(S):    Pawel Netzel - URK
6 
+ *
7 
+ * PURPOSE:      Program calculates multiscale fores fragmentation measure
8 
+ *               with using Jensen-Shannon Similarity measuring similarity
9 
+ *               to fully forested areas.
10 
+ *
11 
+ * COPYRIGHT:    (C) 2023 by Pawel Netzel
12 
+ *
13 
+ *               This program is a free software under the GNU General Public
14 
+ *               License (>=v2).
15 
+ *
16 
+ ***************************************************************************/
17 
+
18 
+#include <stdlib.h>
19 
+#include <string.h>
20 
+#include <omp.h>
21 
+#include <gdal.h>
22 
+#include <ezgdal.h>
23 
+#include <argtable3.h>
24 
+
25 
+/****
26 
+ *  The procedure displays program usage
27 
+ *  and returns to the system
28 
+ */
29 
+void usage(char *progname, void *argtable) {
30 
+
31 
+	printf("\nplForestFragment calculates forest fragmentation index based on \n");
32 
+	printf("similarity approach.\n");
33 
+	printf("It can calculate the index for one scale. It can also produce \n");
34 
+	printf("the multiscale fragmentation index. \n");
35 
+	printf("plForestFragment takes a categorical raster layer as an input.\n");
36 
+	printf("A user can point a category that represents the forest class.\n");
37 
+
38 
+	printf("\nUsage:\n\t%s", progname);
39 
+	arg_print_syntax(stdout,argtable,"\n");
40 
+	printf("\n");
41 
+	arg_print_glossary_gnu(stdout,argtable);
42 
+	printf("\n");
43 
+	exit(0);
44 
+}
45 
+
46 
+/****
47 
+ *  Comparator to sort radiuses
48 
+ */
49 
+int compare(const void* numA, const void* numB)
50 
+{
51 
+    const int* num1 = (const int*)numA;
52 
+    const int* num2 = (const int*)numB;
53 
+
54 
+    if (*num1 > *num2) {
55 
+        return 1;
56 
+    }
57 
+    else {
58 
+        if (*num1 == *num2)
59 
+            return 0;
60 
+        else
61 
+            return -1;
62 
+    }
63 
+}
64 
+
65 
+/****
66 
+ *  Calculating index in the coocurence matrix
67 
+ *  based on values of categories
68 
+ */
69 
+int calc_cooc(double v1, double v2, double forest) {
70 
+    if(v1==forest && v2==forest) return 0;
71 
+    else if(v1==forest) return 1;
72 
+    else if(v2==forest) return 2;
73 
+    else return 3;
74 
+}
75 
+
76 
+/****
77 
+ *  Creating multiscale matrix that will be used as a mask
78 
+ *  to get circular neighbourhood
79 
+ */
80 
+int** create_mask(int* rad, int count) {
81 
+    int radius, radius_max=rad[count-1];
82 
+    int size = 2*radius_max+1;
83 
+    int i, r, c;
84 
+    int** mask = (int**)malloc(size*sizeof(int*));
85 
+    double x, y;
86 
+
87 
+    for(i=0; i<size; i++)
88 
+        mask[i] = (int*)calloc(size, sizeof(int));
89 
+    for(r=0; r<size; r++) {
90 
+        y = pow(r - radius_max, 2);
91 
+        for(c=0; c<size; c++) {
92 
+            x = pow(c - radius_max, 2);
93 
+            if(sqrt(x+y)<=radius_max)
94 
+                mask[r][c] = count-1;
95 
+            else
96 
+                mask[r][c] = count+1;
97 
+        }
98 
+    }
99 
+    for(i=count-2; i>-1; i--) {
100 
+        radius = rad[i];
101 
+        for(r=0; r<size; r++) {
102 
+            y = pow(r - radius_max, 2);
103 
+            for(c=0; c<size; c++) {
104 
+                x = pow(c - radius_max, 2);
105 
+                if(sqrt(x+y)<=radius)
106 
+                    mask[r][c] = i;
107 
+            }
108 
+        }
109 
+    }
110 
+    return mask;
111 
+}
112 
+
113 
+void free_mask(int** mask, int size) {
114 
+    int i;
115 
+    for(i=0; i<size; i++) free(mask[i]);
116 
+    free(mask);
117 
+}
118 
+
119 
+/****
120 
+ *  The function calculates Jensen-Shannon Similarity - JSS
121 
+ *  between compund coocurence matrix and reference matrix
122 
+ */
123 
+double calc_similarity(int** cooc, int count, double weight, double* reference) {
124 
+    int i, j, pos=0;
125 
+    double tsum=0.0, sum, M=0.0, HM=0.0, HA=0.0, HB=0.0;
126 
+    double* norm =(double*)calloc(count*4,sizeof(double));
127 
+    // normalize
128 
+    for(i=0; i<count; i++) {
129 
+        sum = 0.0;
130 
+        for(j=0; j<4; j++)
131 
+            sum += cooc[i][j];
132 
+        tsum += sum;
133 
+        sum /= weight;
134 
+        for(j=0; j<4; j++)
135 
+            if(sum>0)
136 
+                norm[pos++] = (double)cooc[i][j]/sum;
137 
+    }
138 
+    if(tsum==0.0) return 0.0;
139 
+    for(i=0; i<count*4; i++) {
140 
+        M = 0.5*(norm[i] + reference[i]);
141 
+        HM -= (M>0.0)?M*log2(M):0.0;
142 
+        HA -= (norm[i]>0)?norm[i]*log2(norm[i]):0.0;
143 
+        HB -= (reference[i]>0.0)?reference[i]*log2(reference[i]):0.0;
144 
+    }
145 
+    free(norm);
146 
+    return 1.0-sqrt(HM-0.5*(HA+HB));
147 
+}
148 
+
149 
+/****
150 
+ *   Main calculation procedure
151 
+ *   It fills output layer with calculated similarities.
152 
+ *   The procedure uses one or more scales defines by
153 
+ *   variables: count and rad
154 
+ */
155 
+void calc_forest_fragmentation(EZGDAL_LAYER* in_fd, double forest,
156 
+                               EZGDAL_LAYER* out_fd,
157 
+                               int* rad, int count) {
158 
+    EZGDAL_STRIPE* stripe = in_fd->stripe;
159 
+    int layer_cols = in_fd->cols, layer_rows = in_fd->rows;
160 
+    double* buffer = out_fd->buffer;
161 
+    double in_no_data = in_fd->no_data;
162 
+    double out_no_data = out_fd->no_data;
163 
+    double weight = 1.0/(double)count;
164 
+    int i, row, col;
165 
+    int radius_max = rad[count-1];
166 
+    double* reference = (double*)calloc(count*4,sizeof(double));
167 
+    int** mask = create_mask(rad, count);
168 
+
169 
+    for(i=0; i<count; i++)
170 
+        reference[i*4] = weight;
171 
+
172 
+    ezgdal_show_message_nolf(stdout,"Calculations...   0%");
173 
+    /****
174 
+     *  For each row of input raster
175 
+     */
176 
+    for(row=0; row<layer_rows; row++) {
177 
+        ezgdal_load_stripe_data(stripe, row-radius_max);
178 
+        ezgdal_show_progress(stderr,row,layer_rows);
179 
+        /****
180 
+         *  For each cell int the row - parallel processing
181 
+         */
182 
+#pragma omp parallel for default(shared) schedule(dynamic, 100)
183 
+        for(col=0; col<layer_cols; col++) {
184 
+            int i, c, r;
185 
+            EZGDAL_FRAME* frame = &(stripe->frame[col]);
186 
+            double** frame_buffer = frame->buffer;
187 
+            if(frame_buffer[radius_max][radius_max]!=in_no_data) {
188 
+                int** cooc = (int **)malloc(count*sizeof(int*));
189 
+                for(i=0; i<count; i++)
190 
+                    cooc[i] = (int*)calloc(4, sizeof(int));
191 
+                /****
192 
+                 *  For each scale
193 
+                 */
194 
+                for(i=0; i<count; i++) {
195 
+                    int radius = rad[i];
196 
+                    int* cooc_scale = cooc[i];
197 
+                    /****
198 
+                     *  Calculation of the coocurence matrix
199 
+                     */
200 
+                    for(r=radius_max-radius; r<radius_max+radius; r++) {
201 
+                        for(c=radius_max-radius; c<radius_max+radius; c++) {
202 
+                            if(mask[r][c]<=i) {
203 
+                                double v = frame_buffer[r][c];
204 
+                                if(v!=in_no_data) {
205 
+                                    double v_right = frame_buffer[r][c+1];
206 
+                                    double v_down = frame_buffer[r+1][c];
207 
+                                    if(v_right!=in_no_data)
208 
+                                        cooc_scale[calc_cooc(v, v_right, forest)]++;
209 
+                                    if(v_down!=in_no_data)
210 
+                                        cooc_scale[calc_cooc(v, v_down, forest)]++;
211 
+                                }
212 
+                            }
213 
+                        }
214 
+                    }
215 
+                }
216 
+                /****
217 
+                 *  Calculate of similarity based on counpound coocurence matrix
218 
+                 */
219 
+                buffer[col] = calc_similarity(cooc, count, weight, reference);
220 
+                for(i=0; i<count; i++)
221 
+                    free(cooc[i]);
222 
+                free(cooc);
223 
+            } else
224 
+                buffer[col] = out_no_data;
225 
+        }
226 
+        ezgdal_write_buffer(out_fd, row);
227 
+    }
228 
+    ezgdal_show_progress(stderr,row,layer_rows);
229 
+
230 
+    free(reference);
231 
+    free_mask(mask, 2*radius_max+1);
232 
+}
233 
+
234 
+
235 
+int main(int argc, char *argv[])
236 
+{
237 
+    EZGDAL_LAYER *in_fd, *out_fd;
238 
+    char *OUTPUT;
239 
+    char *INPUT;
240 
+    double forest = 1.0;
241 
+
242 
+    struct arg_str  *inp   = arg_str1("i","input","<file_name>","name of input files (GeoTIFF)");
243 
+    struct arg_dbl  *frst  = arg_dbl0("f","forest-value","<dbl>","value for forest (default: 1.0)");
244 
+    struct arg_str  *out   = arg_str1("o","output","<file_name>","name of output file with forest fragmentation index (JSS) (GeoTIFF)");
245 
+    struct arg_int  *rad   = arg_intn("r","radius","<n>",1,8,"radius(es) in cells defining scale(s)");
246 
+    struct arg_int  *th    = arg_int0("t",NULL,"<n>","number of threads (default: 1)");
247 
+    struct arg_lit  *help  = arg_lit0("h","help","print help and exit");
248 
+    struct arg_end  *end   = arg_end(20);
249 
+
250 
+    void* argtable[] = {inp,frst,out,rad,th,help,end};
251 
+
252 
+    int nerrors = arg_parse(argc,argv,argtable);
253 
+
254 
+    if (help->count > 0)
255 
+      usage(argv[0],argtable);
256 
+
257 
+    if (nerrors > 0) {
258 
+      arg_print_errors(stdout,end,argv[0]);
259 
+      usage(argv[0],argtable);
260 
+    }
261 
+
262 
+    /****
263 
+     *  Setting the number of threads for OpenMP - default 1
264 
+     */
265 
+    if(th->count > 0)
266 
+      omp_set_num_threads(th->ival[0]);
267 
+    else
268 
+      omp_set_num_threads(1);
269 
+
270 
+    /****
271 
+     *  Setting the value for forested cells - default 1
272 
+     */
273 
+    if(frst->count > 0)
274 
+      forest = frst->dval[0];
275 
+
276 
+    INPUT = (char *)inp->sval[0];
277 
+    OUTPUT = (char *)out->sval[0];
278 
+
279 
+    /****
280 
+     *  Sorting radiuses from the smallest to the largest
281 
+     */
282 
+    qsort(rad->ival, rad->count, sizeof(int), compare);
283 
+
284 
+    /****
285 
+     *  Opening input raster layer
286 
+     */
287 
+    if(!ezgdal_file_exists(INPUT)) {
288 
+      printf("\nThe file '%s' does not exist!\n\n",INPUT);
289 
+      exit(0);
290 
+    }
291 
+
292 
+    ezgdal_show_message_nolf(stdout,"Opening input file...");
293 
+    in_fd = ezgdal_open_layer(INPUT);
294 
+
295 
+    if(in_fd==NULL) {
296 
+      printf("\nThe file '%s' can't be opened!\n\n",INPUT);
297 
+      exit(0);
298 
+    } else
299 
+      ezgdal_show_message(stdout," OK");
300 
+
301 
+    ezgdal_create_stripe(in_fd, -rad->ival[rad->count-1], 2*rad->ival[rad->count-1]+1);
302 
+    ezgdal_create_all_frames(in_fd->stripe, -rad->ival[rad->count-1], 1);
303 
+
304 
+    /****
305 
+     *  Creating output raster layer
306 
+     */
307 
+    int nodata = -9999;
308 
+    ezgdal_show_message_nolf(stdout,"Creating output file...");
309 
+    out_fd = ezgdal_create_layer(OUTPUT,
310 
+                                 ezgdal_layer_get_wkt(in_fd),
311 
+                                 "Float32",
312 
+                                 ezgdal_layer_get_at(in_fd),
313 
+                                 in_fd->rows,
314 
+                                 in_fd->cols,
315 
+                                 &nodata, EZGDAL_COMPRESS_DEFLATE
316 
+                                );
317 
+    if(out_fd==NULL) {
318 
+      printf("\nThe file '%s' can't be created!\n\n",OUTPUT);
319 
+      exit(0);
320 
+    } else
321 
+      ezgdal_show_message(stdout," OK");
322 
+
323 
+    /****
324 
+     *  Calculation of the forest fragmentation measure
325 
+     */
326 
+    calc_forest_fragmentation(in_fd, forest, out_fd, rad->ival, rad->count);
327 
+
328 
+
329 
+    /****
330 
+     *  Cleaning up and closing files
331 
+     */
332 
+    ezgdal_close_layer(out_fd);
333 
+    ezgdal_close_layer(in_fd);
334 
+    arg_free(argtable);
335 
+
336 
+    exit(EXIT_SUCCESS);
337 
+}
0 338