1 // Copyright (c) 2018-2021 Robert J. Hijmans
2 //
3 // This file is part of the "spat" library.
4 //
5 // spat is free software: you can redistribute it and/or modify it
6 // under the terms of the GNU General Public License as published by
7 // the Free Software Foundation, either version 2 of the License, or
8 // (at your option) any later version.
9 //
10 // spat is distributed in the hope that it will be useful, but
11 // WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with spat. If not, see <http://www.gnu.org/licenses/>.
17
18 //#include <vector>
19 #include "spatRasterMultiple.h"
20 #include "recycle.h"
21 #include "vecmath.h"
22 //#include "vecmath.h"
23 #include <cmath>
24 #include "math_utils.h"
25 #include "file_utils.h"
26 #include "string_utils.h"
27
28
29 /*
30 std::vector<double> flat(std::vector<std::vector<double>> v) {
31 unsigned s1 = v.size();
32 unsigned s2 = v[0].size();
33
34 std::size_t s = s1 * s2;
35 std::vector<double> result(s);
36 for (size_t i=0; i<s1; i++) {
37 for (size_t j=0; j<s2; j++) {
38 result[i*s2+j] = v[i][j];
39 }
40 }
41 return result;
42 }
43 */
44
45 /*
46 SpatRaster SpatRaster::selectHighest(size_t n, bool low, SpatOptions &opt) {
47
48 SpatVector out;
49
50 if (nlyr() > 1) {
51 SpatOptions ops(opt);
52 out.addWarning("only processing the first layer");
53 std::vector<unsigned> lyr = {0};
54 *this = subset(lyr, ops);
55 }
56 if (!hasValues()) {
57 return(out);
58 }
59 if (n >= ncell()) {
60 return isnotnan(opt);
61 }
62
63 std::vector<double> sel;
64
65 if (!readStart()) {
66 return(out);
67 }
68
69 BlockSize bs = getBlockSize(opt);
70 for (size_t i = 0; i < bs.n; i++) {
71 std::vector<double> v = readBlock(bs, i);
72 for (size_t j=0; j<v.size(); j++) {
73 }
74
75 readStop();
76 return(out);
77 }
78 */
79
get_aggregate_dims(std::vector<unsigned> & fact,std::string & message)80 bool SpatRaster::get_aggregate_dims(std::vector<unsigned> &fact, std::string &message ) {
81
82 unsigned fs = fact.size();
83 if ((fs > 3) | (fs == 0)) {
84 message = "argument 'fact' should have length 1, 2, or 3";
85 return false;
86 }
87 auto min_value = *std::min_element(fact.begin(),fact.end());
88 if (min_value < 1) {
89 message = "values in argument 'fact' should be > 0";
90 return false;
91 }
92 auto max_value = *std::max_element(fact.begin(),fact.end());
93 if (max_value == 1) {
94 message = "all values in argument 'fact' are 1, nothing to do";
95 return false;
96 }
97
98 fact.resize(6);
99 if (fs == 1) {
100 fact[1] = fact[0];
101 fact[2] = 1;
102 } else if (fs == 2) {
103 fact[2] = 1;
104 }
105 // int dy = dim[0], dx = dim[1], dz = dim[2];
106 fact[0] = fact[0] < 1 ? 1 : fact[0];
107 fact[0] = fact[0] > nrow() ? nrow() : fact[0];
108 fact[1] = fact[1] < 1 ? 1 : fact[1];
109 fact[1] = fact[1] > ncol() ? ncol() : fact[1];
110 fact[2] = std::max(unsigned(1), std::min(fact[2], nlyr()));
111 // new dimensions: rows, cols, lays
112 fact[3] = std::ceil(double(nrow()) / fact[0]);
113 fact[4] = std::ceil(double(ncol()) / fact[1]);
114 fact[5] = std::ceil(double(nlyr()) / fact[2]);
115 return true;
116 }
117
118
get_aggregate_dims2(std::vector<unsigned> fact)119 std::vector<unsigned> SpatRaster::get_aggregate_dims2(std::vector<unsigned> fact) {
120 // for use with R
121 std::string message = "";
122 get_aggregate_dims(fact, message);
123 return(fact);
124 }
125
126
get_aggregates(std::vector<double> & in,size_t nr,std::vector<unsigned> dim)127 std::vector<std::vector<double>> SpatRaster::get_aggregates(std::vector<double> &in, size_t nr, std::vector<unsigned> dim) {
128
129 // dim 0, 1, 2, are the aggregations factors dy, dx, dz
130 // and 3, 4, 5 are the new nrow, ncol, nlyr
131
132 // adjust for chunk
133 //dim[3] = std::ceil(double(nr) / dim[0]);
134 //size_t bpC = dim[3];
135 size_t bpC = std::ceil(double(nr) / dim[0]);
136
137 size_t dy = dim[0], dx = dim[1], dz = dim[2];
138 size_t bpR = dim[4];
139 size_t bpL = bpR * bpC;
140
141 // new number of layers
142 size_t newNL = dim[5];
143
144 // new number of rows, adjusted for additional (expansion) rows
145 size_t adjnr = bpC * dy;
146
147 // number of aggregates
148 size_t nblocks = (bpR * bpC * newNL);
149 // cells per aggregate
150 size_t blockcells = dx * dy * dz;
151
152 // output: each row is a block
153 std::vector< std::vector<double> > a(nblocks, std::vector<double>(blockcells, std::numeric_limits<double>::quiet_NaN()));
154
155 size_t nc = ncol();
156 // size_t ncells = ncell();
157 size_t ncells = nr * nc;
158 size_t nl = nlyr();
159 size_t lstart, rstart, cstart, lmax, rmax, cmax, f, lj, cell;
160
161 for (size_t b = 0; b < nblocks; b++) {
162 lstart = dz * (b / bpL);
163 rstart = (dy * (b / bpR)) % adjnr;
164 cstart = dx * (b % bpR);
165
166 lmax = std::min(nl, (lstart + dz));
167 rmax = std::min(nr, (rstart + dy)); // nrow -> nr
168 cmax = std::min(nc, (cstart + dx));
169
170 f = 0;
171 for (size_t j = lstart; j < lmax; j++) {
172 lj = j * ncells;
173 for (size_t r = rstart; r < rmax; r++) {
174 cell = lj + r * nc;
175 for (size_t c = cstart; c < cmax; c++) {
176 a[b][f] = in[cell + c];
177 f++;
178 }
179 }
180 }
181 }
182 return(a);
183 }
184
185
compute_aggregates(std::vector<double> & in,size_t nr,size_t nc,size_t nl,std::vector<unsigned> dim,std::function<double (std::vector<double> &,bool)> fun,bool narm)186 std::vector<double> compute_aggregates(std::vector<double> &in, size_t nr, size_t nc, size_t nl, std::vector<unsigned> dim, std::function<double(std::vector<double>&, bool)> fun, bool narm) {
187
188 // dim 0, 1, 2, are the aggregations factors dy, dx, dz
189 // and 3, 4, 5 are the new nrow, ncol, nlyr
190
191 size_t dy = dim[0], dx = dim[1], dz = dim[2];
192 // size_t bpC = dim[3];
193 // adjust for chunk
194 size_t bpC = std::ceil(double(nr) / dim[0]);
195 size_t bpR = dim[4];
196 size_t bpL = bpR * bpC;
197
198 // new number of layers
199 size_t newNL = dim[5];
200
201 // new number of rows, adjusted for additional (expansion) rows
202 size_t adjnr = bpC * dy;
203
204 // number of aggregates
205 size_t nblocks = (bpR * bpC * newNL);
206 // cells per aggregate
207 size_t blockcells = dx * dy * dz;
208
209 // output: each row is a block
210 std::vector<double> out(nblocks, NAN);
211
212 // size_t nl = nlyr();
213 // size_t nc = ncol();
214 size_t ncells = nr * nc;
215 size_t lstart, rstart, cstart, lmax, rmax, cmax, f, lj, cell;
216
217 for (size_t b = 0; b < nblocks; b++) {
218 lstart = dz * (b / bpL);
219 rstart = (dy * (b / bpR)) % adjnr;
220 cstart = dx * (b % bpR);
221
222 lmax = std::min(nl, (lstart + dz));
223 rmax = std::min(nr, (rstart + dy)); // nrow -> nr
224 cmax = std::min(nc, (cstart + dx));
225
226 f = 0;
227 std::vector<double> a(blockcells, NAN);
228 for (size_t j = lstart; j < lmax; j++) {
229 lj = j * ncells;
230 for (size_t r = rstart; r < rmax; r++) {
231 cell = lj + r * nc;
232 for (size_t c = cstart; c < cmax; c++) {
233 a[f] = in[cell + c];
234 f++;
235 }
236 }
237 }
238 out[b] = fun(a, narm);
239 }
240 return(out);
241 }
242
243
244
aggregate(std::vector<unsigned> fact,std::string fun,bool narm,SpatOptions & opt)245 SpatRaster SpatRaster::aggregate(std::vector<unsigned> fact, std::string fun, bool narm, SpatOptions &opt) {
246
247 SpatRaster out;
248 std::string message = "";
249 bool success = get_aggregate_dims(fact, message);
250
251 // fact 0, 1, 2, are the aggregation factors dy, dx, dz
252 // and 3, 4, 5 are the new nrow, ncol, nlyr
253 if (!success) {
254 out.setError(message);
255 return out;
256 }
257
258 SpatExtent extent = getExtent();
259 double xmax = extent.xmin + fact[4] * fact[1] * xres();
260 double ymin = extent.ymax - fact[3] * fact[0] * yres();
261 SpatExtent e = SpatExtent(extent.xmin, xmax, ymin, extent.ymax);
262 out = SpatRaster(fact[3], fact[4], fact[5], e, "");
263 out.source[0].srs = source[0].srs;
264 // there is much more. categories, time. should use geometry and then
265 // set extent and row col
266 if (fact[5] == nlyr()) {
267 out.setNames(getNames());
268 }
269
270 if (!source[0].hasValues) {
271 return out;
272 }
273
274 if (!haveFun(fun)) {
275 out.setError("unknown function argument");
276 return out;
277 }
278
279 /*
280 size_t ifun = std::distance(f.begin(), it);
281 std::string gstring = "";
282 if (ifun > 0) {
283 std::vector<std::string> gf {"average", "min", "max", "med", "mode"};
284 gstring = gf[ifun-1];
285 }
286
287 #ifdef useGDAL
288 #if GDAL_VERSION_MAJOR >= 3
289 if (gstring != "") {
290 out = warper(out, "", gstring, opt);
291 return out;
292 }
293 #endif
294 #endif
295 */
296
297 std::function<double(std::vector<double>&, bool)> agFun = getFun(fun);
298
299 unsigned outnc = out.ncol();
300
301 //BlockSize bs = getBlockSize(4, opt.get_memfrac());
302 BlockSize bs = getBlockSize(opt);
303 //bs.n = floor(nrow() / fact[0]); # ambiguous on solaris
304 bs.n = std::floor(static_cast <double> (nrow() / fact[0]));
305
306 bs.nrows = std::vector<size_t>(bs.n, fact[0]);
307 bs.row.resize(bs.n);
308 for (size_t i =0; i<bs.n; i++) {
309 bs.row[i] = i * fact[0];
310 }
311 size_t lastrow = bs.row[bs.n - 1] + bs.nrows[bs.n - 1] + 1;
312 if (lastrow < nrow()) {
313 bs.row.push_back(lastrow);
314 bs.nrows.push_back(std::min(bs.nrows[bs.n-1], nrow()-lastrow));
315 bs.n += 1;
316 }
317 if (!readStart()) {
318 out.setError(getError());
319 return(out);
320 }
321
322 opt.steps = bs.n;
323 opt.minrows = fact[0];
324
325 if (fun == "modal") {
326 if (nlyr() == out.nlyr()) {
327 out.source[0].hasColors = hasColors();
328 out.source[0].cols = getColors();
329 out.source[0].hasCategories = hasCategories();
330 out.source[0].cats = getCategories();
331 }
332 }
333
334 if (!out.writeStart(opt)) {
335 readStop();
336 return out;
337 }
338
339 size_t nc = ncol();
340 for (size_t i = 0; i < bs.n; i++) {
341 std::vector<double> vin = readValues(bs.row[i], bs.nrows[i], 0, nc);
342 std::vector<double> v = compute_aggregates(vin, bs.nrows[i], nc, nlyr(), fact, agFun, narm);
343 if (!out.writeValues(v, i, 1, 0, outnc)) return out;
344 }
345 out.writeStop();
346 readStop();
347 return(out);
348 }
349
350
351
352
weighted_mean(SpatRaster w,bool narm,SpatOptions & opt)353 SpatRaster SpatRaster::weighted_mean(SpatRaster w, bool narm, SpatOptions &opt) {
354 SpatRaster out;
355 if (nlyr() != w.nlyr()) {
356 out.setError("nlyr of data and weights are different");
357 return out;
358 }
359
360 SpatOptions topt(opt);
361 out = arith(w, "*", topt);
362 out = out.summary("sum", narm, topt);
363 SpatRaster wsum = w.summary("sum", narm, topt);
364 return out.arith(wsum, "/", opt);
365
366 }
367
368
weighted_mean(std::vector<double> w,bool narm,SpatOptions & opt)369 SpatRaster SpatRaster::weighted_mean(std::vector<double> w, bool narm, SpatOptions &opt) {
370 SpatOptions topt(opt);
371 recycle(w, nlyr());
372 SpatRaster out = arith(w, "*", false, topt);
373 out = out.summary("sum", narm, topt);
374 double wsum = vsum(w, narm);
375 return out.arith(wsum, "/", false, opt);
376 }
377
378
separate(std::vector<double> classes,double keepvalue,double othervalue,SpatOptions & opt)379 SpatRaster SpatRaster::separate(std::vector<double> classes, double keepvalue, double othervalue, SpatOptions &opt) {
380
381 SpatRaster out;
382 if (nlyr() > 1) {
383 out.setError("input may only have one layer");
384 return out;
385 }
386 if (classes.size() == 0) {
387 SpatOptions topt(opt);
388 std::vector<std::vector<double>> rc = unique(false, topt);
389 classes = rc[0];
390 }
391
392 std::vector<int> uc(classes.size());
393 for (size_t i=0; i<classes.size(); i++) {
394 uc[i] = round(classes[i]);
395 }
396 std::sort(uc.begin(), uc.end());
397 uc.erase(std::unique(uc.begin(), uc.end()), uc.end());
398
399 size_t n = uc.size();
400 if (n == 0) {
401 out.setError("no valid classes");
402 return out;
403 }
404 out = geometry(n);
405 std::vector<std::string> snms(n);
406 for (size_t i=0; i<n; i++) {
407 snms[i] = std::to_string(uc[i]);
408 }
409 out.setNames(snms);
410
411 if (!readStart()) {
412 out.setError(getError());
413 return(out);
414 }
415 if (!out.writeStart(opt)) {
416 readStop();
417 return out;
418 }
419
420 for (size_t i = 0; i < out.bs.n; i++) {
421 std::vector<double> v = readBlock(out.bs, i);
422 size_t nn = v.size();
423 std::vector<double> vv(nn * n, NAN);
424 for (size_t j=0; j<nn; j++) {
425 if (!std::isnan(v[j])) {
426 for (size_t k=0; k<uc.size(); k++) {
427 if (v[j] == uc[k]) {
428 if (keepvalue) {
429 vv[j + k*nn] = uc[k];
430 } else {
431 vv[j + k*nn] = 1; // true
432 }
433 } else {
434 vv[j + k*nn] = othervalue;
435 }
436 }
437
438 }
439 }
440 if (!out.writeValues(vv, out.bs.row[i], out.bs.nrows[i], 0, ncol())) {
441 readStop();
442 return out;
443 }
444 }
445 readStop();
446 out.writeStop();
447 return(out);
448 }
449
450
is_in(std::vector<double> m,SpatOptions & opt)451 SpatRaster SpatRaster::is_in(std::vector<double> m, SpatOptions &opt) {
452
453 SpatRaster out = geometry();
454 if (m.size() == 0) {
455 out.setError("no matches supplied");
456 return(out);
457 }
458 if (!hasValues()) {
459 out.setError("input has no values");
460 return(out);
461 }
462
463 int hasNAN = 0;
464 for (size_t i=0; i<m.size(); i++) {
465 if (std::isnan(m[i])) {
466 hasNAN = 1;
467 m.erase(m.begin()+i);
468 break;
469 }
470 }
471 if (m.size() == 0) { // only NA
472 return isnan(opt);
473 }
474
475
476 // if m is very long, perhaps first check if the value is in range?
477
478 if (!readStart()) {
479 out.setError(getError());
480 return(out);
481 }
482
483 if (!out.writeStart(opt)) {
484 readStop();
485 return out;
486 }
487 for (size_t i = 0; i < out.bs.n; i++) {
488 std::vector<double> v = readBlock(out.bs, i);
489 std::vector<double> vv(v.size(), 0);
490 for (size_t j=0; j<v.size(); j++) {
491 if (std::isnan(v[j])) {
492 vv[j] = hasNAN;
493 } else {
494 for (size_t k=0; k<m.size(); k++) {
495 if (v[j] == m[k]) {
496 vv[j] = 1;
497 break;
498 }
499 }
500 }
501 }
502
503 if (!out.writeValues(vv, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
504 }
505 readStop();
506 out.writeStop();
507 return(out);
508 }
509
510
511
is_in_cells(std::vector<double> m,SpatOptions & opt)512 std::vector<std::vector<double>> SpatRaster::is_in_cells(std::vector<double> m, SpatOptions &opt) {
513
514 std::vector<std::vector<double>> out(nlyr());
515
516 if (m.size() == 0) {
517 return(out);
518 }
519 if (!hasValues()) {
520 return(out);
521 }
522 bool hasNAN = false;
523 for (size_t i=0; i<m.size(); i++) {
524 if (std::isnan(m[i])) {
525 hasNAN = true;
526 m.erase(m.begin()+i);
527 break;
528 }
529 }
530 // if (m.size() == 0) { // only NA
531 //nanOnly=true;
532 // }
533
534 if (!readStart()) {
535 return(out);
536 }
537
538 BlockSize bs = getBlockSize(opt);
539 size_t nc = ncol();
540 for (size_t i = 0; i < bs.n; i++) {
541 std::vector<double> v = readBlock(bs, i);
542 size_t cellperlayer = bs.nrows[i] * nc;
543 for (size_t j=0; j<v.size(); j++) {
544 size_t lyr = j / cellperlayer;
545 size_t cell = j % cellperlayer + bs.row[i] * nc;
546 if (std::isnan(v[j])) {
547 if (hasNAN) out[lyr].push_back(cell);
548 } else {
549 for (size_t k=0; k<m.size(); k++) {
550 if (v[j] == m[k]) {
551 out[lyr].push_back(cell);
552 break;
553 }
554 }
555 }
556 }
557 }
558 readStop();
559 return(out);
560 }
561
562
563
564
stretch(std::vector<double> minv,std::vector<double> maxv,std::vector<double> minq,std::vector<double> maxq,std::vector<double> smin,std::vector<double> smax,SpatOptions & opt)565 SpatRaster SpatRaster::stretch(std::vector<double> minv, std::vector<double> maxv, std::vector<double> minq, std::vector<double> maxq, std::vector<double> smin, std::vector<double> smax, SpatOptions &opt) {
566
567 SpatRaster out = geometry();
568 if (!hasValues()) return(out);
569
570 size_t nl = nlyr();
571 recycle(minv, nl);
572 recycle(maxv, nl);
573 recycle(minq, nl);
574 recycle(maxq, nl);
575 recycle(smin, nl);
576 recycle(smax, nl);
577
578 std::vector<std::vector<double>> q(nl);
579 std::vector<bool> useS(nl, false);
580 std::vector<double> mult(nl);
581
582 for (size_t i=0; i<nl; i++) {
583 if (minv[i] >= maxv[i]) {
584 out.setError("maxv must be larger than minv");
585 return out;
586 }
587 if ((!std::isnan(smin[i])) && (!std::isnan(smax[i]))) {
588 if (smin[i] >= smax[i]) {
589 out.setError("smax must be larger than smin");
590 return out;
591 }
592 useS[i] = true;
593 q[i] = {smin[i], smax[i]};
594 } else {
595 if (minq[i] >= maxq[i]) {
596 out.setError("maxq must be larger than minq");
597 return out;
598 }
599 if ((minq[i] < 0) || (maxq[i] > 1)) {
600 out.setError("minq and maxq must be between 0 and 1");
601 return out;
602 }
603 }
604 }
605
606 std::vector<bool> hR = hasRange();
607 for (size_t i=0; i<nl; i++) {
608 if (!useS[i]) {
609 if ((minq[i]==0) && (maxq[i]==1) && hR[i]) {
610 std::vector<double> rmn = range_min();
611 std::vector<double> rmx = range_max();
612 q[i] = {rmn[i], rmx[i]};
613 } else {
614 std::vector<double> probs = {minq[i], maxq[i]};
615 SpatOptions xopt(opt);
616 std::vector<double> v = getValues(i, xopt);
617 q[i] = vquantile(v, probs, true);
618 }
619 }
620 mult[i] = maxv[i] / (q[i][1]-q[i][0]);
621 }
622
623 if (!readStart()) {
624 out.setError(getError());
625 return(out);
626 }
627
628 if (!out.writeStart(opt)) {
629 readStop();
630 return out;
631 }
632 for (size_t i = 0; i < out.bs.n; i++) {
633 std::vector<double> v = readBlock(out.bs, i);
634 size_t nc = out.bs.nrows[i] * ncol();
635 for (size_t j=0; j<v.size(); j++) {
636 size_t lyr = j / nc;
637 v[j] = mult[lyr] * (v[j] - q[lyr][0]);
638 if (v[j] < minv[lyr]) v[j] = minv[lyr];
639 if (v[j] > maxv[lyr]) v[j] = maxv[lyr];
640 }
641 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
642 }
643 readStop();
644 out.writeStop();
645
646 return(out);
647 }
648
649
apply(std::vector<unsigned> ind,std::string fun,bool narm,std::vector<std::string> nms,SpatOptions & opt)650 SpatRaster SpatRaster::apply(std::vector<unsigned> ind, std::string fun, bool narm, std::vector<std::string> nms, SpatOptions &opt) {
651
652 recycle(ind, nlyr());
653 std::vector<unsigned> ui = vunique(ind);
654 unsigned nl = ui.size();
655 SpatRaster out = geometry(nl);
656 recycle(nms, nl);
657 out.setNames(nms);
658
659 if (!haveFun(fun)) {
660 out.setError("unknown function argument");
661 return out;
662 }
663
664 if (!hasValues()) return(out);
665
666 if (!readStart()) {
667 out.setError(getError());
668 return(out);
669 }
670 if (!out.writeStart(opt)) {
671 readStop();
672 return out;
673 }
674 out.bs = getBlockSize(opt);
675 // #ifdef useRcpp
676 // out.pbar = new Progress(out.bs.n+2, opt.show_progress(bs.n));
677 // out.pbar->increment();
678 // #endif
679
680 std::vector<std::vector<double>> v(nl);
681 std::vector<unsigned> ird(ind.size());
682 std::vector<unsigned> jrd(ind.size());
683 for (size_t i=0; i<nl; i++) {
684 for (size_t j=0; j<ind.size(); j++) {
685 if (ui[i] == ind[j]) {
686 v[i].push_back(0);
687 ird[j] = i;
688 jrd[j] = v[i].size()-1;
689 }
690 }
691 }
692
693 std::function<double(std::vector<double>&, bool)> theFun = getFun(fun);
694
695 for (size_t i=0; i<out.bs.n; i++) {
696 std::vector<double> a = readBlock(out.bs, i);
697 unsigned nc = out.bs.nrows[i] * ncol();
698 std::vector<double> b(nc * nl);
699 for (size_t j=0; j<nc; j++) {
700 for (size_t k=0; k<ird.size(); k++) {
701 v[ird[k]][jrd[k]] = a[j+k*nc];
702 }
703 for (size_t k=0; k<ui.size(); k++) {
704 size_t off = k * nc + j;
705 b[off] = theFun(v[k], narm);
706 }
707 }
708 if (!out.writeValues(b, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
709 }
710 readStop();
711 out.writeStop();
712 return(out);
713 }
714
715
716
mask(SpatRaster x,bool inverse,double maskvalue,double updatevalue,SpatOptions & opt)717 SpatRaster SpatRaster::mask(SpatRaster x, bool inverse, double maskvalue, double updatevalue, SpatOptions &opt) {
718
719 unsigned nl = std::max(nlyr(), x.nlyr());
720 SpatRaster out = geometry(nl, true);
721
722 if (!out.compare_geom(x, false, true, opt.get_tolerance(), true, true, true, false)) {
723 return(out);
724 }
725
726 if (!readStart()) {
727 out.setError(getError());
728 return(out);
729 }
730 if (!x.readStart()) {
731 out.setError(x.getError());
732 return(out);
733 }
734 if (!out.writeStart(opt)) {
735 readStop();
736 return out;
737 }
738 std::vector<double> v, m;
739 for (size_t i = 0; i < out.bs.n; i++) {
740 v = readValues(out.bs.row[i], out.bs.nrows[i], 0, ncol());
741 m = x.readValues(out.bs.row[i], out.bs.nrows[i], 0, ncol());
742 recycle(v, m);
743 if (inverse) {
744 if (std::isnan(maskvalue)) {
745 for (size_t i=0; i < v.size(); i++) {
746 if (!std::isnan(m[i])) {
747 v[i] = updatevalue;
748 }
749 }
750 } else {
751 for (size_t i=0; i < v.size(); i++) {
752 if (m[i] != maskvalue) {
753 v[i] = updatevalue;
754 }
755 }
756 }
757 } else {
758 if (std::isnan(maskvalue)) {
759 for (size_t i=0; i < v.size(); i++) {
760 if (std::isnan(m[i])) {
761 v[i] = updatevalue;
762 }
763 }
764 } else {
765 for (size_t i=0; i < v.size(); i++) {
766 if (m[i] == maskvalue) {
767 v[i] = updatevalue;
768 }
769 }
770 }
771 }
772 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
773
774 }
775 out.writeStop();
776 readStop();
777 x.readStop();
778 return(out);
779 }
780
781
782
mask(SpatRaster x,bool inverse,std::vector<double> maskvalues,double updatevalue,SpatOptions & opt)783 SpatRaster SpatRaster::mask(SpatRaster x, bool inverse, std::vector<double> maskvalues, double updatevalue, SpatOptions &opt) {
784
785 maskvalues = vunique(maskvalues);
786 if (maskvalues.size() == 1) {
787 return mask(x, inverse, maskvalues[0], updatevalue, opt);
788 }
789
790 unsigned nl = std::max(nlyr(), x.nlyr());
791 SpatRaster out = geometry(nl, true);
792
793 if (!out.compare_geom(x, false, true, opt.get_tolerance(), true, true, true, false)) {
794 return(out);
795 }
796
797 if (!readStart()) {
798 out.setError(getError());
799 return(out);
800 }
801 if (!x.readStart()) {
802 out.setError(x.getError());
803 return(out);
804 }
805
806 bool maskNA = false;
807 for (int i = maskvalues.size()-1; i>=0; i--) {
808 if (std::isnan(maskvalues[i])) {
809 maskNA = true;
810 maskvalues.erase(maskvalues.begin()+i);
811 }
812 }
813
814 if (!out.writeStart(opt)) {
815 readStop();
816 return out;
817 }
818 std::vector<double> v, m;
819 for (size_t i = 0; i < out.bs.n; i++) {
820 v = readValues(out.bs.row[i], out.bs.nrows[i], 0, ncol());
821 m = x.readValues(out.bs.row[i], out.bs.nrows[i], 0, ncol());
822 recycle(v, m);
823 if (inverse) {
824 for (size_t i=0; i < v.size(); i++) {
825 if (maskNA && std::isnan(m[i])) {
826 v[i] = updatevalue;
827 } else {
828 for (size_t j=0; j < maskvalues.size(); j++) {
829 if (m[i] != maskvalues[j]) {
830 v[i] = updatevalue;
831 break;
832 }
833 }
834 }
835 }
836 } else {
837 for (size_t i=0; i < v.size(); i++) {
838 if (maskNA && std::isnan(m[i])) {
839 v[i] = updatevalue;
840 } else {
841 for (size_t j=0; j < maskvalues.size(); j++) {
842 if (m[i] == maskvalues[j]) {
843 v[i] = updatevalue;
844 break;
845 }
846 }
847 }
848 }
849 }
850 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
851
852 }
853 out.writeStop();
854 readStop();
855 x.readStop();
856 return(out);
857 }
858
859
mask(SpatVector x,bool inverse,double updatevalue,bool touches,SpatOptions & opt)860 SpatRaster SpatRaster::mask(SpatVector x, bool inverse, double updatevalue, bool touches, SpatOptions &opt) {
861
862 SpatRaster out;
863 if (!hasValues()) {
864 out.setError("SpatRaster has no values");
865 return out;
866 }
867 if (inverse) {
868 out = rasterizeLyr(x, updatevalue, NAN, touches, true, opt);
869 } else {
870 SpatOptions topt(opt);
871 out = rasterizeLyr(x, 1, 0, touches, false, topt);
872 if (out.hasError()) {
873 return out;
874 }
875 if (std::isnan(updatevalue)) {
876 out = mask(out, false, 0, updatevalue, opt);
877 } else {
878 out = mask(out, false, 0, updatevalue, topt);
879 out = out.mask(*this, false, NAN, NAN, opt);
880 }
881 }
882 return(out);
883 }
884
885
886
transpose(SpatOptions & opt)887 SpatRaster SpatRaster::transpose(SpatOptions &opt) {
888
889 SpatRaster out = geometry(nlyr(), true);
890 SpatExtent eold = getExtent();
891 SpatExtent enew = getExtent();
892 enew.xmin = eold.ymin;
893 enew.xmax = eold.ymax;
894 enew.ymin = eold.xmin;
895 enew.ymax = eold.xmax;
896 out.setExtent(enew, false, "");
897 out.source[0].ncol = nrow();
898 out.source[0].nrow = ncol();
899 if (!hasValues()) return out;
900 if (!readStart()) {
901 out.setError(getError());
902 return(out);
903 }
904 if (!out.writeStart(opt)) {
905 readStop();
906 return out;
907 }
908 for (size_t i=0; i < out.bs.n; i++) {
909 unsigned nr = nrow();
910 unsigned nc = out.bs.nrows[i];
911 std::vector<double> v = readValues(0, nr, out.bs.row[i], nc);
912 std::vector<double> vv(v.size());
913 for (size_t lyr=0; lyr<nlyr(); lyr++) {
914 size_t off = lyr*ncell();
915 for (size_t r = 0; r < nr; r++) {
916 size_t rnc = off + r * nc;
917 for (size_t c = 0; c < nc; c++) {
918 vv[c*nr+r+off] = v[rnc+c];
919 }
920 }
921 }
922 if (!out.writeValues(vv, out.bs.row[i], out.bs.nrows[i], 0, out.ncol())) return out;
923 }
924 out.writeStop();
925 readStop();
926 return(out);
927 }
928
929
930
931
932
trim(double value,unsigned padding,SpatOptions & opt)933 SpatRaster SpatRaster::trim(double value, unsigned padding, SpatOptions &opt) {
934
935 long nrl = nrow() * nlyr();
936 long ncl = ncol() * nlyr();
937
938 std::vector<double> v;
939 size_t r;
940 size_t nr = nrow();
941 bool rowfound = false;
942 if (!readStart()) {
943 SpatRaster out;
944 out.setError(getError());
945 return(out);
946 }
947
948
949 size_t firstrow, lastrow, firstcol, lastcol;
950 if (std::isnan(value)) {
951 for (r=0; r<nr; r++) {
952 v = readValues(r, 1, 0, ncol());
953 if (std::count_if( v.begin(), v.end(), [](double d) { return std::isnan(d); } ) < ncl) {
954 rowfound = true;
955 break;
956 }
957 }
958
959 if (!rowfound) {
960 SpatRaster out;
961 out.setError("only cells with NA found");
962 return out;
963 }
964
965 firstrow = std::max(r - padding, size_t(0));
966
967 for (r=nrow()-1; r>firstrow; r--) {
968 v = readValues(r, 1, 0, ncol());
969 if (std::count_if(v.begin(), v.end(), [](double d) { return std::isnan(d); } ) < ncl) {
970 break;
971 }
972 }
973
974 lastrow = std::max(std::min(r+padding, nrow()), size_t(0));
975
976 if (lastrow < firstrow) {
977 std::swap(firstrow, lastrow);
978 }
979 size_t c;
980 for (c=0; c<ncol(); c++) {
981 v = readValues(0, nrow(), c, 1);
982 if (std::count_if( v.begin(), v.end(), [](double d) { return std::isnan(d); } ) < nrl) {
983 break;
984 }
985 }
986 firstcol = std::min(std::max(c-padding, size_t(0)), ncol());
987
988 for (c=ncol()-1; c>firstcol; c--) {
989 v = readValues(0, nrow(), c, 1);
990 if (std::count_if( v.begin(), v.end(), [](double d) { return std::isnan(d); } ) < nrl) {
991 break;
992 }
993 }
994 lastcol = std::max(std::min(c+padding, ncol()), size_t(0));
995 } else {
996 for (r=0; r<nr; r++) {
997 v = readValues(r, 1, 0, ncol());
998 if (std::count( v.begin(), v.end(), value) < ncl) {
999 rowfound = true;
1000 break;
1001 }
1002 }
1003
1004 if (!rowfound) {
1005 SpatRaster out;
1006 out.setError("only cells with value: " + std::to_string(value) + " found");
1007 return out;
1008 }
1009
1010 firstrow = std::max(r - padding, size_t(0));
1011
1012 for (r=nrow()-1; r>firstrow; r--) {
1013 v = readValues(r, 1, 0, ncol());
1014 if (std::count( v.begin(), v.end(), value) < ncl) {
1015 break;
1016 }
1017 }
1018
1019 lastrow = std::max(std::min(r+padding, nrow()), size_t(0));
1020
1021 if (lastrow < firstrow) {
1022 std::swap(firstrow, lastrow);
1023 }
1024 size_t c;
1025 for (c=0; c<ncol(); c++) {
1026 v = readValues(0, nrow(), c, 1);
1027 if (std::count( v.begin(), v.end(), value) < nrl) {
1028 break;
1029 }
1030 }
1031 firstcol = std::min(std::max(c-padding, size_t(0)), ncol());
1032
1033
1034 for (c=ncol()-1; c>firstcol; c--) {
1035 v = readValues(0, nrow(), c, 1);
1036 if (std::count( v.begin(), v.end(), value) < nrl) {
1037 break;
1038 }
1039 }
1040 lastcol = std::max(std::min(c+padding, ncol()), size_t(0));
1041
1042 }
1043 readStop();
1044 if (lastcol < firstcol) {
1045 std::swap(firstcol, lastcol);
1046 }
1047
1048 std::vector<double> res = resolution();
1049 double xr = res[0];
1050 double yr = res[1];
1051 SpatExtent e = SpatExtent(xFromCol(firstcol)-0.5*xr, xFromCol(lastcol)+0.5*xr, yFromRow(lastrow)-0.5*yr, yFromRow(firstrow)+0.5*yr);
1052
1053 return( crop(e, "near", opt) ) ;
1054 }
1055
1056
1057
1058
1059
clamp_vector(std::vector<double> & v,double low,double high,bool usevalue)1060 void clamp_vector(std::vector<double> &v, double low, double high, bool usevalue) {
1061 size_t n = v.size();
1062 if (usevalue) {
1063 for (size_t i=0; i<n; i++) {
1064 if ( v[i] < low ) {
1065 v[i] = low;
1066 } else if ( v[i] > high ) {
1067 v[i] = high;
1068 }
1069 }
1070 } else {
1071 for (size_t i=0; i<n; i++) {
1072 if ( (v[i] < low )| (v[i] > high)) {
1073 v[i] = NAN;
1074 }
1075 }
1076 }
1077 }
1078
1079
1080
clamp(double low,double high,bool usevalue,SpatOptions & opt)1081 SpatRaster SpatRaster::clamp(double low, double high, bool usevalue, SpatOptions &opt) {
1082
1083 SpatRaster out = geometry(nlyr(), true);
1084 if (low > high) {
1085 out.setError("lower clamp value cannot be larger than the higher clamp value");
1086 return out;
1087 }
1088 if (!hasValues()) {
1089 out.setError("cannot clamp a raster with no values");
1090 return out;
1091 }
1092
1093 if (!readStart()) {
1094 out.setError(getError());
1095 return(out);
1096 }
1097
1098 if (!out.writeStart(opt)) {
1099 readStop();
1100 return out;
1101 }
1102 for (size_t i = 0; i < out.bs.n; i++) {
1103 std::vector<double> v = readBlock(out.bs, i);
1104 clamp_vector(v, low, high, usevalue);
1105 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1106 }
1107 readStop();
1108 out.writeStop();
1109 return(out);
1110 }
1111
1112
1113
1114
selRange(SpatRaster x,int z,int recycleby,SpatOptions & opt)1115 SpatRaster SpatRaster::selRange(SpatRaster x, int z, int recycleby, SpatOptions &opt) {
1116
1117 int nl = nlyr();
1118 z = std::max(1, std::min(z, nl));
1119 size_t nrecs = 1;
1120 if (recycleby > 1 && recycleby < nl) {
1121 nrecs = nl / recycleby;
1122 } else {
1123 recycleby = 0;
1124 }
1125 SpatRaster out = geometry( z * nrecs );
1126 if (!out.compare_geom(x, false, false, opt.get_tolerance())) {
1127 return(out);
1128 }
1129 if (!hasValues()) return(out);
1130
1131 if (x.nlyr() > 1) {
1132 out.setError("index raster must have only one layer");
1133 return out;
1134 }
1135 if (!x.hasValues()) {
1136 out.setError("index raster has no values");
1137 return out;
1138 }
1139
1140 if (!readStart()) {
1141 out.setError(getError());
1142 return(out);
1143 }
1144 if (!x.readStart()) {
1145 out.setError(x.getError());
1146 return(out);
1147 }
1148
1149 if (!out.writeStart(opt)) {
1150 readStop();
1151 return out;
1152 }
1153 for (size_t i=0; i<out.bs.n; i++) {
1154 std::vector<double> v = readBlock(out.bs, i);
1155 std::vector<double> idx = x.readBlock(out.bs, i);
1156 size_t is = idx.size();
1157 std::vector<double> vv(is*z*nrecs, NAN);
1158 size_t ncell = out.bs.nrows[i] * ncol(); // same as is?
1159
1160 for (size_t j=0; j<is; j++) { //index.size (each cell)
1161 for (size_t k=0; k<nrecs; k++) {
1162 int start = idx[j] - 1 + k * recycleby; // first layer for this cell
1163 int offbase = (k*z) * ncell;
1164 if ((start >= 0) && (start < nl)) {
1165 int zz = std::min(nl-start, z); // do not surpass the last layer
1166 for (int i=0; i<zz; i++){
1167 size_t offin = (start + i) * ncell + j;
1168 size_t offout = offbase + i * ncell + j;
1169 vv[offout] = v[offin];
1170 }
1171 }
1172 }
1173 }
1174 //for (size_t j=0; j<is; j++) {
1175 // int index = idx[j] - 1;
1176 // if ((index >= 0) && (index < nl)) {
1177 // vv[j] = v[j + index * ncell];
1178 // }
1179 //}
1180 if (!out.writeValues(vv, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1181 }
1182 readStop();
1183 x.readStop();
1184 out.writeStop();
1185 return(out);
1186 }
1187
1188
rapply(SpatRaster x,double first,double last,std::string fun,bool clamp,bool narm,SpatOptions & opt)1189 SpatRaster SpatRaster::rapply(SpatRaster x, double first, double last, std::string fun, bool clamp, bool narm, SpatOptions &opt) {
1190
1191 SpatRaster out = geometry(1);
1192 if (!haveFun(fun)) {
1193 out.setError("unknown function argument");
1194 return out;
1195 }
1196
1197 bool sval = !std::isnan(first);
1198 bool eval = !std::isnan(last);
1199 if (sval && eval) {
1200 out.setError("arguments `first` or `last` must be NA. See `app` for other cases");
1201 return out;
1202 }
1203 int start = sval ? first-1 : -99;
1204 int end = eval ? last-1 : -999;
1205
1206 if (!out.compare_geom(x, false, false, opt.get_tolerance())) {
1207 return(out);
1208 }
1209 if (!x.hasValues()) {
1210 out.setError("index raster has no values");
1211 return out;
1212 }
1213 unsigned expnl = 2 - (sval + eval);
1214 if (x.nlyr() != expnl) {
1215 out.setError("index raster must have " + std::to_string(expnl) + "layer(s)");
1216 return out;
1217 }
1218 if (!hasValues()) {
1219 out.setError("no values in input");
1220 return(out);
1221 }
1222
1223 std::function<double(std::vector<double>&, bool)> theFun = getFun(fun);
1224
1225 int nl = nlyr();
1226 if (!out.writeStart(opt)) {
1227 readStop();
1228 return out;
1229 }
1230 if (!readStart()) {
1231 out.setError(getError());
1232 return(out);
1233 }
1234 if (!x.readStart()) {
1235 out.setError(x.getError());
1236 return(out);
1237 }
1238
1239 for (size_t i=0; i<out.bs.n; i++) {
1240 std::vector<double> v = readBlock(out.bs, i);
1241 std::vector<double> idx = x.readBlock(out.bs, i);
1242 size_t ncell = out.bs.nrows[i] * ncol();
1243 std::vector<double> vv(ncell, NAN);
1244 for (size_t j=0; j<ncell; j++) {
1245 if (std::isnan(idx[j])) continue;
1246 if (sval) {
1247 end = idx[j] - 1;
1248 } else if (eval) {
1249 start = idx[j] - 1;
1250 } else {
1251 start = idx[j] - 1;
1252 double dend = idx[j+ncell]-1;
1253 if (std::isnan(dend)) continue;
1254 end = dend;
1255 }
1256 if (clamp) {
1257 start = start < 0 ? 0 : start;
1258 end = end >= nl ? (nl-1) : end;
1259 }
1260 if ((start <= end) && (end < nl) && (start >= 0)) {
1261 std::vector<double> se;
1262 se.reserve(end-start+1);
1263 for (int k = start; k<=end; k++){
1264 size_t off = k * ncell + j;
1265 se.push_back(v[off]);
1266 }
1267 vv[j] = theFun(se, narm);
1268 }
1269 }
1270 if (!out.writeValues(vv, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1271 }
1272 readStop();
1273 x.readStop();
1274 out.writeStop();
1275 return(out);
1276 }
1277
1278
rappvals(SpatRaster x,double first,double last,bool clamp,bool all,double fill,size_t startrow,size_t nrows)1279 std::vector<std::vector<double>> SpatRaster::rappvals(SpatRaster x, double first, double last, bool clamp, bool all, double fill, size_t startrow, size_t nrows) {
1280
1281 std::vector<std::vector<double>> r;
1282
1283 bool sval = !std::isnan(first);
1284 bool eval = !std::isnan(last);
1285 if (sval && eval) {
1286 setError("first or last must be NA. See `app` for other cases");
1287 return r;
1288 }
1289 int start = sval ? first-1 : 0;
1290 int end = eval ? last-1 : 0;
1291
1292 if (!compare_geom(x, false, false, 0.1)) {
1293 return(r);
1294 }
1295 if (!hasValues()) {
1296 return r;
1297 }
1298 if (!x.hasValues()) {
1299 setError("index raster has no values");
1300 return r;
1301 }
1302 unsigned expnl = 2 - (sval + eval);
1303 if (x.nlyr() != expnl) {
1304 setError("index raster must have " + std::to_string(expnl) + "layer(s)");
1305 return r;
1306 }
1307
1308 int nl = nlyr();
1309 if (!readStart()) {
1310 return(r);
1311 }
1312 if (!x.readStart()) {
1313 setError(x.getError());
1314 return(r);
1315 }
1316
1317 std::vector<double> v = readValues(startrow, nrows, 0, ncol());
1318 std::vector<double> idx = x.readValues(startrow, nrows, 0, ncol());
1319 size_t ncell = nrows * ncol();
1320 r.resize(ncell);
1321
1322 for (size_t j=0; j<ncell; j++) {
1323 if (std::isnan(idx[j])) {
1324 if (all) {
1325 r[j].resize(nl, NAN);
1326 } else {
1327 r[j].push_back(NAN);
1328 }
1329 continue;
1330 }
1331 if (sval) {
1332 end = idx[j] - 1;
1333 //end = idx[j];
1334 } else if (eval) {
1335 start = idx[j] - 1;
1336 } else {
1337 start = idx[j] - 1;
1338 //double dend = idx[j+ncell];
1339 double dend = idx[j+ncell]-1;
1340 end = std::isnan(dend) ? -999 : (int) dend;
1341 }
1342 if (clamp) {
1343 start = start < 0 ? 0 : start;
1344 end = end >= nl ? (nl-1) : end;
1345 }
1346
1347 bool inrange = (start <= end) && (end < nl) && (start >= 0);
1348 if (all) {
1349 if (inrange) {
1350 r[j].resize(nl, fill);
1351 for (int k = start; k<=end; k++){
1352 size_t off = k * ncell + j;
1353 r[j][k] = v[off];
1354 }
1355 } else {
1356 r[j].resize(nl, NAN);
1357 }
1358 } else if (inrange) {
1359 r[j].reserve(end-start+1);
1360 for (int k=start; k<=end; k++){
1361 size_t off = k * ncell + j;
1362 r[j].push_back(v[off]);
1363 }
1364 } else {
1365 r[j].push_back(NAN);
1366 }
1367
1368 }
1369 readStop();
1370 x.readStop();
1371 return(r);
1372 }
1373
1374
disaggregate_dims(std::vector<unsigned> & fact,std::string & message)1375 bool disaggregate_dims(std::vector<unsigned> &fact, std::string &message ) {
1376
1377 unsigned fs = fact.size();
1378 if ((fs > 3) | (fs == 0)) {
1379 message = "argument 'fact' should have length 1, 2, or 3";
1380 return false;
1381 }
1382 auto min_value = *std::min_element(fact.begin(),fact.end());
1383 if (min_value < 1) {
1384 message = "values in argument 'fact' should be > 0";
1385 return false;
1386 }
1387 auto max_value = *std::max_element(fact.begin(),fact.end());
1388 if (max_value == 1) {
1389 message = "all values in argument 'fact' are 1, nothing to do";
1390 return false;
1391 }
1392
1393 fact.resize(3);
1394 if (fs == 1) {
1395 fact[1] = fact[0];
1396 }
1397 fact[2] = 1;
1398 return true;
1399 }
1400
1401
1402
disaggregate(std::vector<unsigned> fact,SpatOptions & opt)1403 SpatRaster SpatRaster::disaggregate(std::vector<unsigned> fact, SpatOptions &opt) {
1404
1405 SpatRaster out = geometry(nlyr(), true);
1406 std::string message = "";
1407 bool success = disaggregate_dims(fact, message);
1408 if (!success) {
1409 out.setError(message);
1410 return out;
1411 }
1412
1413 out.source[0].nrow = out.source[0].nrow * fact[0];
1414 out.source[0].ncol = out.source[0].ncol * fact[1];
1415 out.source[0].nlyr = out.source[0].nlyr * fact[2];
1416
1417 if (!hasValues()) {
1418 return out;
1419 }
1420
1421 opt.ncopies = 2*fact[0]*fact[1]*fact[2];
1422 BlockSize bs = getBlockSize(opt);
1423 //opt.set_blocksizemp();
1424 std::vector<double> v, vout;
1425 unsigned nc = ncol();
1426 unsigned nl = nlyr();
1427 std::vector<double> newrow(nc*fact[1]);
1428 if (!readStart()) {
1429 out.setError(getError());
1430 return(out);
1431 }
1432
1433 if (!out.writeStart(opt)) {
1434 readStop();
1435 return out;
1436 }
1437 for (size_t i = 0; i < bs.n; i++) {
1438 v = readValues(bs.row[i], bs.nrows[i], 0, nc);
1439 vout.resize(0);
1440 vout.reserve(v.size() * fact[0] * fact[1] * fact[2]);
1441
1442 for (size_t lyr=0; lyr<nl; lyr++) {
1443 for (size_t row=0; row<bs.nrows[i]; row++) {
1444 unsigned rowoff = row*nc + lyr*nc*bs.nrows[i];
1445 // for each new column
1446 unsigned jfact = 0;
1447 for (size_t j=0; j<nc; j++) {
1448 unsigned coloff = rowoff + j;
1449 for (size_t k=0; k<fact[1]; k++) {
1450 newrow[jfact+k] = v[coloff];
1451 }
1452 jfact += fact[1];
1453 }
1454 // for each new row
1455 for (size_t j=0; j<fact[0]; j++) {
1456 vout.insert(vout.end(), newrow.begin(), newrow.end());
1457 }
1458 }
1459 }
1460 if (!out.writeValues(vout, bs.row[i]*fact[0], bs.nrows[i]*fact[0], 0, out.ncol())) return out;
1461 }
1462 vout.resize(0);
1463 out.writeStop();
1464 readStop();
1465 return(out);
1466 }
1467
1468
1469
1470 /*
1471 SpatRaster SpatRaster::oldinit(std::string value, bool plusone, SpatOptions &opt) {
1472
1473 SpatRaster out = geometry(1);
1474
1475 std::vector<std::string> f {"row", "col", "cell", "x", "y", "chess"};
1476 bool test = std::find(f.begin(), f.end(), value) == f.end();
1477 if (test) {
1478 out.setError("not a valid init option");
1479 return out;
1480 }
1481
1482 size_t nr = nrow();
1483 size_t steps = nr; // for the pbar
1484 if (value == "chess") {
1485 steps = steps / 2;
1486 }
1487 opt.set_steps(steps);
1488 if (!out.writeStart(opt)) { return out; }
1489
1490 if (value == "row") {
1491 std::vector<double> v(ncol());
1492 for (size_t i = 0; i < nr; i++) {
1493 std::fill(v.begin(), v.end(), i+plusone);
1494 if (!out.writeValues(v, i, 1, 0, ncol())) return out;
1495 }
1496 } else if (value == "col") {
1497 std::vector<double> col(ncol());
1498 double start = plusone ? 1 : 0;
1499 std::iota(col.begin(), col.end(), start);
1500 for (unsigned i = 0; i < nr; i++) {
1501 if (!out.writeValues(col, i, 1, 0, ncol())) return out;
1502 }
1503 } else if (value == "cell") {
1504 std::vector<long> col(ncol());
1505 std::iota(col.begin(), col.end(), 0);
1506 std::vector<long> row(1);
1507 for (unsigned i = 0; i < nr; i++) {
1508 row[0] = i;
1509 std::vector<double> v = cellFromRowCol(row, col);
1510 if (plusone) for(double& d : v) d=d+1;
1511 if (!out.writeValues(v, i, 1, 0, ncol())) return out;
1512 }
1513 } else if (value == "x") {
1514 std::vector<long> col(ncol());
1515 std::iota(col.begin(), col.end(), 0);
1516 std::vector<double> x = xFromCol(col);
1517 for (unsigned i = 0; i < nr; i++) {
1518 if (!out.writeValues(x, i, 1, 0, ncol())) return out;
1519 }
1520 } else if (value == "y") {
1521 std::vector<double> v(ncol());
1522 for (unsigned i = 0; i < nr; i++) {
1523 double y = yFromRow(i);
1524 std::fill(v.begin(), v.end(), y);
1525 if (!out.writeValues(v, i, 1, 0, ncol())) return out;
1526 }
1527 } else if (value == "chess") {
1528 std::vector<double> a(ncol());
1529 std::vector<double> b(ncol());
1530 size_t nr = nrow();
1531 size_t nc = ncol();
1532 a[0] = 1;
1533 b[0] = 0;
1534 for (size_t i=1; i<nc; i++) {
1535 bool test = i%2 == 0;
1536 a[i] = test;
1537 b[i] = !test;
1538 }
1539 out.bs.n = nr/2; // for the pbar
1540 for (unsigned i=0; i<(nr-1); i=i+2) {
1541 if (!out.writeValues(a, i, 1, 0, ncol())) return out;
1542 if (!out.writeValues(b, i+1, 1, 0, ncol())) return out;
1543 }
1544 if (nr%2 == 0) {
1545 if (!out.writeValues(a, nr-2, 1, 0, ncol())) return out;
1546 if (!out.writeValues(b, nr-1, 1, 0, ncol())) return out;
1547 } else {
1548 if (!out.writeValues(a, nr-1, 1, 0, ncol())) return out;
1549 }
1550 }
1551
1552 out.writeStop();
1553 return(out);
1554 }
1555 */
1556
1557
1558
init(std::string value,bool plusone,SpatOptions & opt)1559 SpatRaster SpatRaster::init(std::string value, bool plusone, SpatOptions &opt) {
1560
1561 SpatRaster out = geometry(1);
1562
1563 std::vector<std::string> f {"row", "col", "cell", "x", "y", "chess"};
1564 bool test = std::find(f.begin(), f.end(), value) == f.end();
1565 if (test) {
1566 out.setError("not a valid init option");
1567 return out;
1568 }
1569
1570 opt.ncopies = std::max(opt.ncopies, (unsigned) 6);
1571 if (!out.writeStart(opt)) {
1572 readStop();
1573 return out;
1574 }
1575
1576 size_t nc = ncol();
1577 std::vector<double> v;
1578 if (value == "row") {
1579 for (size_t i = 0; i < out.bs.n; i++) {
1580 v.resize(nc * out.bs.nrows[i]);
1581 for (size_t j = 0; j < out.bs.nrows[i]; j++) {
1582 size_t r = out.bs.row[i] + j + plusone;
1583 for (size_t k = 0; k < nc; k++) {
1584 v[j*nc+k] = r;
1585 }
1586 }
1587 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
1588 }
1589 } else if (value == "col") {
1590 std::vector<double> cnn(nc);
1591 double start = plusone ? 1 : 0;
1592 std::iota(cnn.begin(), cnn.end(), start);
1593 size_t oldnr = 0;
1594 for (size_t i = 0; i < out.bs.n; i++) {
1595 if (oldnr != out.bs.nrows[i]) {
1596 v = cnn;
1597 recycle(v, out.bs.nrows[i] * nc);
1598 }
1599 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
1600 }
1601 } else if (value == "cell") {
1602 for (size_t i = 0; i < out.bs.n; i++) {
1603 v.resize(nc * out.bs.nrows[i]);
1604 size_t firstcell = cellFromRowCol(out.bs.row[i], 0);
1605 firstcell = plusone ? firstcell + 1 : firstcell;
1606 std::iota(v.begin(), v.end(), firstcell);
1607 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
1608 }
1609 } else if (value == "x") {
1610 std::vector<int_64> col(nc);
1611 std::iota(col.begin(), col.end(), 0);
1612 std::vector<double> xcoords = xFromCol(col);
1613 size_t oldnr = 0;
1614 for (size_t i = 0; i < out.bs.n; i++) {
1615 if (oldnr != out.bs.nrows[i]) {
1616 v = xcoords;
1617 recycle(v, out.bs.nrows[i] * nc);
1618 }
1619 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
1620 }
1621 } else if (value == "y") {
1622
1623 for (size_t i = 0; i < out.bs.n; i++) {
1624 v.resize(out.bs.nrows[i] * nc );
1625 for (size_t j = 0; j < out.bs.nrows[i]; j++) {
1626 double y = yFromRow(out.bs.row[i] + j);
1627 for (size_t k = 0; k < nc; k++) {
1628 v[j*nc+k] = y;
1629 }
1630 }
1631 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
1632 }
1633
1634 } else if (value == "chess") {
1635 std::vector<double> a(nc);
1636 std::vector<double> b(nc);
1637 for (size_t i=0; i<nc; i++) {
1638 bool even = i%2 == 0;
1639 a[i] = even;
1640 b[i] = !even;
1641 }
1642 std::vector<double> v;
1643 for (size_t i = 0; i < out.bs.n; i++) {
1644 if ((out.bs.row[i]%2) == 0) {
1645 v = a;
1646 v.insert(v.end(), b.begin(), b.end());
1647 } else {
1648 v = b;
1649 v.insert(v.end(), b.begin(), b.end());
1650 }
1651 recycle(v, out.bs.nrows[i] * nc);
1652 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
1653 }
1654 }
1655
1656 out.writeStop();
1657 return(out);
1658 }
1659
1660
init(std::vector<double> values,SpatOptions & opt)1661 SpatRaster SpatRaster::init(std::vector<double> values, SpatOptions &opt) {
1662 SpatRaster out = geometry();
1663 if (!out.writeStart(opt)) { return out; }
1664 unsigned nc = ncol();
1665 unsigned nl = nc * nlyr();
1666 if (values.size() == 1) {
1667 std::vector<double> v;
1668 for (size_t i = 0; i < out.bs.n; i++) {
1669 v.resize(out.bs.nrows[i]*nc*nl, values[0]);
1670 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
1671 }
1672 } else {
1673 int over = 0;
1674 for (size_t i = 0; i < out.bs.n; i++) {
1675 if (over > 0) {
1676 std::vector<double> newv(values.begin()+over, values.end());
1677 newv.insert(newv.end(), values.begin(), values.begin()+over);
1678 values = newv;
1679 }
1680 std::vector<double> v = values;
1681 recycle(v, out.bs.nrows[i]*nc);
1682 recycle(v, out.bs.nrows[i]*nc*nl);
1683 over = v.size() % values.size();
1684 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
1685 }
1686 }
1687 out.writeStop();
1688 return(out);
1689 }
1690
1691
isnan(SpatOptions & opt)1692 SpatRaster SpatRaster::isnan(SpatOptions &opt) {
1693 SpatRaster out = geometry();
1694 if (!hasValues()) return out;
1695 if (!readStart()) {
1696 out.setError(getError());
1697 return(out);
1698 }
1699
1700 if (!out.writeStart(opt)) {
1701 readStop();
1702 return out;
1703 }
1704 for (size_t i=0; i<out.bs.n; i++) {
1705 std::vector<double> v = readBlock(out.bs, i);
1706 for (double &d : v) d = std::isnan(d);
1707 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1708 }
1709 readStop();
1710 out.writeStop();
1711 return(out);
1712 }
1713
1714
isnotnan(SpatOptions & opt)1715 SpatRaster SpatRaster::isnotnan(SpatOptions &opt) {
1716 SpatRaster out = geometry();
1717 if (!hasValues()) return out;
1718
1719 if (!readStart()) {
1720 out.setError(getError());
1721 return(out);
1722 }
1723 if (!out.writeStart(opt)) {
1724 readStop();
1725 return out;
1726 }
1727 for (size_t i=0; i<out.bs.n; i++) {
1728 std::vector<double> v = readBlock(out.bs, i);
1729 for (double &d : v) d = ! std::isnan(d);
1730 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1731 }
1732 readStop();
1733 out.writeStop();
1734 return(out);
1735 }
1736
1737
isfinite(SpatOptions & opt)1738 SpatRaster SpatRaster::isfinite(SpatOptions &opt) {
1739 SpatRaster out = geometry();
1740 if (!hasValues()) return out;
1741
1742 if (!readStart()) {
1743 out.setError(getError());
1744 return(out);
1745 }
1746 if (!out.writeStart(opt)) {
1747 readStop();
1748 return out;
1749 }
1750 for (size_t i=0; i<out.bs.n; i++) {
1751 std::vector<double> v = readBlock(out.bs, i);
1752 for (double &d : v) d = std::isfinite(d);
1753 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1754 }
1755 readStop();
1756 out.writeStop();
1757 return(out);
1758 }
1759
1760
isinfinite(SpatOptions & opt)1761 SpatRaster SpatRaster::isinfinite(SpatOptions &opt) {
1762 SpatRaster out = geometry();
1763 if (!hasValues()) return out;
1764
1765 if (!readStart()) {
1766 out.setError(getError());
1767 return(out);
1768 }
1769 if (!out.writeStart(opt)) {
1770 readStop();
1771 return out;
1772 }
1773 for (size_t i=0; i<out.bs.n; i++) {
1774 std::vector<double> v = readBlock(out.bs, i);
1775 for (double &d : v) d = std::isinf(d);
1776 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1777 }
1778 readStop();
1779 out.writeStop();
1780 return(out);
1781 }
1782
1783
rotate(bool left,SpatOptions & opt)1784 SpatRaster SpatRaster::rotate(bool left, SpatOptions &opt) {
1785
1786 unsigned nc = ncol();
1787 unsigned nl = nlyr();
1788 unsigned hnc = (nc / 2);
1789 double addx = hnc * xres();
1790 if (left) {
1791 addx = -addx;
1792 }
1793 SpatRaster out = geometry(nlyr(), true);
1794 SpatExtent outext = out.getExtent();
1795 outext.xmin = outext.xmin + addx;
1796 outext.xmax = outext.xmax + addx;
1797 out.setExtent(outext, true);
1798
1799 if (!hasValues()) return out;
1800
1801 if (!readStart()) {
1802 out.setError(getError());
1803 return(out);
1804 }
1805 if (!out.writeStart(opt)) {
1806 readStop();
1807 return out;
1808 }
1809 std::vector<double> b;
1810 for (size_t i=0; i < out.bs.n; i++) {
1811 std::vector<double> a = readBlock(out.bs, i);
1812 for (size_t j=0; j < nl; j++) {
1813 for (size_t r=0; r < out.bs.nrows[i]; r++) {
1814 unsigned s1 = j * out.bs.nrows[i] * nc + r * nc;
1815 unsigned e1 = s1 + hnc;
1816 unsigned s2 = e1;
1817 unsigned e2 = s1 + nc;
1818 b.insert(b.end(), a.begin()+s2, a.begin()+e2);
1819 b.insert(b.end(), a.begin()+s1, a.begin()+e1);
1820 }
1821 }
1822 if (!out.writeValues(b, out.bs.row[i], nrow(), 0, ncol())) return out;
1823 b.resize(0);
1824 }
1825 out.writeStop();
1826 readStop();
1827 return(out);
1828 }
1829
1830
1831
shared_basegeom(SpatRaster & x,double tol,bool test_overlap)1832 bool SpatRaster::shared_basegeom(SpatRaster &x, double tol, bool test_overlap) {
1833 if (!compare_origin(x.origin(), tol)) return false;
1834 if (!about_equal(xres(), x.xres(), xres() * tol)) return false;
1835 if (!about_equal(yres(), x.yres(), yres() * tol)) return false;
1836 if (test_overlap) {
1837 SpatExtent e = x.getExtent();
1838 e.intersect(getExtent());
1839 if (!e.valid()) return false;
1840 }
1841 return true;
1842 }
1843
1844
1845
1846
1847
cover(SpatRaster x,std::vector<double> values,SpatOptions & opt)1848 SpatRaster SpatRaster::cover(SpatRaster x, std::vector<double> values, SpatOptions &opt) {
1849
1850 unsigned nl = std::max(nlyr(), x.nlyr());
1851 SpatRaster out = geometry(nl, true);
1852
1853 bool rmatch = false;
1854 if (out.compare_geom(x, false, false, opt.get_tolerance(), true)) {
1855 rmatch = true;
1856 } else {
1857 // if (!shared_basegeom(x, 0.1, true)) {
1858 out.setError("raster dimensions do not match");
1859 return(out);
1860 // } else {
1861 // out.msg.has_error = false;
1862 // out.msg.error = "";
1863 // SpatExtent e = getExtent();
1864 // SpatExtent xe = x.getExtent();
1865 // double prec = std::min(xres(), yres())/1000;
1866 // if (!xe.compare(e, "<=", prec)) {
1867 // SpatOptions xopt(opt);
1868 // x = x.crop(e, "near", xopt);
1869 // }
1870 // }
1871 }
1872
1873
1874 if (!x.hasValues()) {
1875 return *this;
1876 }
1877 if (!hasValues()) {
1878 if (rmatch) {
1879 return x.deepCopy();
1880 } else {
1881 SpatExtent e = getExtent();
1882 return x.extend(e, opt);
1883 }
1884 }
1885
1886 if (!readStart()) {
1887 out.setError(getError());
1888 return(out);
1889 }
1890 if (!x.readStart()) {
1891 out.setError(x.getError());
1892 return(out);
1893 }
1894
1895 if (!out.writeStart(opt)) {
1896 readStop();
1897 x.readStop();
1898 return out;
1899 }
1900 if (values.size() == 1) {
1901 double value=values[0];
1902 for (size_t i = 0; i < out.bs.n; i++) {
1903 std::vector<double> v = readValues(out.bs.row[i], out.bs.nrows[i], 0, ncol());
1904 std::vector<double> m = x.readValues(out.bs.row[i], out.bs.nrows[i], 0, ncol());
1905 recycle(v, m);
1906 if (std::isnan(value)) {
1907 for (size_t i=0; i < v.size(); i++) {
1908 if (std::isnan(v[i])) {
1909 v[i] = m[i];
1910 }
1911 }
1912 } else {
1913 for (size_t i=0; i < v.size(); i++) {
1914 if (v[i] == value) {
1915 v[i] = m[i];
1916 }
1917 }
1918 }
1919 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1920 }
1921 } else {
1922
1923 values = vunique(values);
1924 bool hasNA = false;
1925 for (int i = values.size()-1; i>=0; i--) {
1926 if (std::isnan(values[i])) {
1927 hasNA = true;
1928 values.erase(values.begin()+i);
1929 }
1930 }
1931
1932 for (size_t i = 0; i < out.bs.n; i++) {
1933 std::vector<double> v = readValues(out.bs.row[i], out.bs.nrows[i], 0, ncol());
1934 std::vector<double> m = x.readValues(out.bs.row[i], out.bs.nrows[i], 0, ncol());
1935 recycle(v, m);
1936 for (size_t i=0; i < v.size(); i++) {
1937 if (hasNA) {
1938 if (std::isnan(v[i])) {
1939 v[i] = m[i];
1940 continue;
1941 }
1942 }
1943 for (size_t i=0; i<values.size(); i++) {
1944 if (v[i] == values[i]) {
1945 v[i] = m[i];
1946 continue;
1947 }
1948 }
1949 }
1950 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
1951 }
1952 }
1953
1954 out.writeStop();
1955 readStop();
1956 x.readStop();
1957 return(out);
1958 }
1959
1960
1961
1962
extend(SpatExtent e,SpatOptions & opt)1963 SpatRaster SpatRaster::extend(SpatExtent e, SpatOptions &opt) {
1964
1965 SpatRaster out = geometry(nlyr(), true);
1966 e = out.align(e, "near");
1967 SpatExtent extent = getExtent();
1968 e.unite(extent);
1969
1970 out.setExtent(e, true);
1971 if (!hasValues() ) {
1972 if (opt.get_filename() != "") {
1973 out.addWarning("ignoring filename argument because there are no cell values");
1974 }
1975 return(out);
1976 }
1977
1978 double tol = std::min(xres(), yres()) / 1000;
1979 if (extent.compare(e, "==", tol)) {
1980 // same extent
1981 if (opt.get_filename() != "") {
1982 out = writeRaster(opt);
1983 } else {
1984 out = deepCopy();
1985 }
1986 return out;
1987 }
1988
1989
1990 if (!readStart()) {
1991 out.setError(getError());
1992 return(out);
1993 }
1994
1995 if (!out.writeStart(opt)) {
1996 readStop();
1997 return out;
1998 }
1999 out.fill(NAN);
2000 BlockSize bs = getBlockSize(opt);
2001 for (size_t i=0; i<bs.n; i++) {
2002 std::vector<double> v = readValues(bs.row[i], bs.nrows[i], 0, ncol());
2003 unsigned row1 = out.rowFromY(yFromRow(bs.row[i]));
2004 unsigned row2 = out.rowFromY(yFromRow(bs.row[i]+bs.nrows[i]-1));
2005 unsigned col1 = out.colFromX(xFromCol(0));
2006 unsigned col2 = out.colFromX(xFromCol(ncol()-1));
2007 if (!out.writeValues(v, row1, row2-row1+1, col1, col2-col1+1)) return out;
2008 }
2009 readStop();
2010 out.writeStop();
2011 return(out);
2012 }
2013
2014
2015
crop(SpatExtent e,std::string snap,SpatOptions & opt)2016 SpatRaster SpatRaster::crop(SpatExtent e, std::string snap, SpatOptions &opt) {
2017
2018 SpatRaster out = geometry(nlyr(), true);
2019
2020 if ( !e.valid() ) {
2021 out.setError("invalid extent");
2022 return out;
2023 }
2024 e.intersect(out.getExtent());
2025 if ( !e.valid() ) {
2026 out.setError("extents do not overlap");
2027 return out;
2028 }
2029
2030 out.setExtent(e, true, snap);
2031 if (!hasValues() ) {
2032 if (opt.get_filename() != "") {
2033 out.addWarning("ignoring filename argument because there are no cell values");
2034 }
2035 return(out);
2036 }
2037
2038 double xr = xres();
2039 double yr = yres();
2040 SpatExtent outext = out.getExtent();
2041 unsigned col1 = colFromX(outext.xmin + 0.5 * xr);
2042 unsigned col2 = colFromX(outext.xmax - 0.5 * xr);
2043 unsigned row1 = rowFromY(outext.ymax - 0.5 * yr);
2044 unsigned row2 = rowFromY(outext.ymin + 0.5 * yr);
2045
2046 std::vector<bool> hw = hasWindow();
2047 bool haswin = hw[0];
2048 for (size_t i=1; i<nsrc(); i++) {
2049 haswin = (haswin | hw[i]);
2050 }
2051
2052 if ((row1==0) && (row2==nrow()-1) && (col1==0) && (col2==ncol()-1) && (!haswin)) {
2053 // same extent
2054 if (opt.get_filename() != "") {
2055 out = writeRaster(opt);
2056 } else {
2057 out = deepCopy();
2058 }
2059 return out;
2060 }
2061
2062 unsigned ncols = out.ncol();
2063 if (!readStart()) {
2064 out.setError(getError());
2065 return(out);
2066 }
2067
2068 // opt.ncopies = 2;
2069 if (!out.writeStart(opt)) {
2070 readStop();
2071 return out;
2072 }
2073 std::vector<double> v;
2074 for (size_t i = 0; i < out.bs.n; i++) {
2075 v = readValues(row1+out.bs.row[i], out.bs.nrows[i], col1, ncols);
2076 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, out.ncol())) return out;
2077 }
2078 out.writeStop();
2079 readStop();
2080
2081 return(out);
2082 }
2083
flip(bool vertical,SpatOptions & opt)2084 SpatRaster SpatRaster::flip(bool vertical, SpatOptions &opt) {
2085
2086 SpatRaster out = geometry(nlyr(), true);
2087 if (!hasValues()) return out;
2088 if (!readStart()) {
2089 out.setError(getError());
2090 return(out);
2091 }
2092
2093 if (!out.writeStart(opt)) {
2094 readStop();
2095 return out;
2096 }
2097 unsigned nc = ncol();
2098 unsigned nl = nlyr();
2099
2100 if (vertical) {
2101 for (size_t i=0; i < out.bs.n; i++) {
2102 std::vector<double> b;
2103 size_t ii = out.bs.n - 1 - i;
2104 std::vector<double> a = readBlock(out.bs, ii);
2105 for (size_t j=0; j < out.nlyr(); j++) {
2106 size_t offset = j * out.bs.nrows[ii] * nc;
2107 for (size_t k=0; k < out.bs.nrows[ii]; k++) {
2108 unsigned start = offset + (out.bs.nrows[ii] - 1 - k) * nc;
2109 b.insert(b.end(), a.begin()+start, a.begin()+start+nc);
2110 }
2111 }
2112 if (!out.writeValues(b, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
2113 }
2114 } else {
2115 for (size_t i=0; i < out.bs.n; i++) {
2116 std::vector<double> b;
2117 std::vector<double> a = readBlock(out.bs, i);
2118 unsigned lyrrows = nl * out.bs.nrows[i];
2119 for (size_t j=0; j < lyrrows; j++) {
2120 unsigned start = j * nc;
2121 unsigned end = start + nc;
2122 std::vector<double> v(a.begin()+start, a.begin()+end);
2123 std::reverse(v.begin(), v.end());
2124 b.insert(b.end(), v.begin(), v.end());
2125 }
2126 if (!out.writeValues(b, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
2127 b.resize(0);
2128 }
2129 }
2130 out.writeStop();
2131 readStop();
2132 return(out);
2133 }
2134
2135
reverse(SpatOptions & opt)2136 SpatRaster SpatRaster::reverse(SpatOptions &opt) {
2137
2138 SpatRaster out = geometry(nlyr(), true);
2139 if (!hasValues()) return out;
2140 if (!readStart()) {
2141 out.setError(getError());
2142 return(out);
2143 }
2144
2145 if (!out.writeStart(opt)) {
2146 readStop();
2147 return out;
2148 }
2149 std::vector<double> b;
2150 unsigned nc = ncol();
2151 unsigned nl = nlyr();
2152
2153 for (size_t i=0; i < out.bs.n; i++) {
2154 size_t ii = out.bs.n - 1 - i;
2155 std::vector<double> a = readBlock(out.bs, ii);
2156 unsigned lyrrows = nl * out.bs.nrows[ii];
2157 for (size_t j=0; j < lyrrows; j++) {
2158 unsigned start = (lyrrows - 1 - j) * nc;
2159 unsigned end = start + nc;
2160 std::vector<double> v(a.begin()+start, a.begin()+end);
2161 std::reverse(v.begin(), v.end());
2162 b.insert(b.end(), v.begin(), v.end());
2163 }
2164 if (!out.writeValues(b, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
2165 b.resize(0);
2166 }
2167
2168 out.writeStop();
2169 readStop();
2170 return(out);
2171 }
2172
2173
shift(double x,double y,SpatOptions & opt)2174 SpatRaster SpatRaster::shift(double x, double y, SpatOptions &opt) {
2175 SpatRaster out = deepCopy();
2176 SpatExtent outext = out.getExtent();
2177 outext.xmin = outext.xmin + x;
2178 outext.xmax = outext.xmax + x;
2179 outext.ymin = outext.ymin + y;
2180 outext.ymax = outext.ymax + y;
2181 out.setExtent(outext, true);
2182 return out;
2183 }
2184
compare_origin(std::vector<double> x,double tol)2185 bool SpatRaster::compare_origin(std::vector<double> x, double tol) {
2186 std::vector<double> y = origin();
2187 if (!about_equal(x[0], y[0], xres() * tol)) return false;
2188 if (!about_equal(x[1], y[1], yres() * tol)) return false;
2189 return true;
2190 }
2191
2192
2193 /*
2194 SpatRaster SpatRasterCollection::merge(SpatOptions &opt) {
2195
2196 SpatRaster out;
2197 unsigned n = size();
2198
2199 if (n == 0) {
2200 out.setError("empty collection");
2201 return(out);
2202 }
2203 if (n == 1) {
2204 out = ds[0].deepCopy();
2205 return(out);
2206 }
2207
2208 bool any_hasvals = false;
2209 if (ds[0].hasValues()) any_hasvals = true;
2210 out = ds[0].geometry(ds[0].nlyr(), true);
2211 std::vector<double> orig = ds[0].origin();
2212 SpatExtent e = ds[0].getExtent();
2213 unsigned nl = ds[0].nlyr();
2214 for (size_t i=1; i<n; i++) {
2215 // lyrs, crs, warncrs, ext, rowcol, res
2216 if (!out.compare_geom(ds[i], false, false, false, false, false, true)) {
2217 return(out);
2218 }
2219 if (!out.compare_origin(ds[i].origin(), 0.1)) {
2220 out.setError("origin of SpatRaster " + std::to_string(i+1) + " does not match the previous SpatRaster(s)");
2221 return(out);
2222 }
2223 e.unite(ds[i].getExtent());
2224 nl = std::max(nl, ds[i].nlyr());
2225 if (ds[i].hasValues()) any_hasvals = true;
2226 }
2227 out.setExtent(e, true);
2228 out = out.geometry(nl, true);
2229 if (!any_hasvals) return out;
2230
2231 // out.setResolution(xres(), yres());
2232 if (!out.writeStart(opt)) { return out; }
2233 out.fill(NAN);
2234
2235 for (size_t i=0; i<n; i++) {
2236 SpatRaster r = ds[i];
2237 if (!r.hasValues()) continue;
2238 BlockSize bs = r.getBlockSize(opt);
2239 if (!r.readStart()) {
2240 out.setError(r.getError());
2241 return(out);
2242 }
2243
2244 for (size_t j=0; j<bs.n; j++) {
2245 std::vector<double> v = r.readValues(bs.row[j], bs.nrows[j], 0, r.ncol());
2246 unsigned row1 = out.rowFromY(r.yFromRow(bs.row[j]));
2247 unsigned row2 = out.rowFromY(r.yFromRow(bs.row[j]+bs.nrows[j]-1));
2248 unsigned col1 = out.colFromX(r.xFromCol(0));
2249 unsigned col2 = out.colFromX(r.xFromCol(r.ncol()-1));
2250 unsigned ncols = col2-col1+1;
2251 unsigned nrows = row2-row1+1;
2252 recycle(v, ncols * nrows * nl);
2253
2254 if (!out.writeValues(v, row1, nrows, col1, ncols)) return out;
2255 }
2256 r.readStop();
2257 }
2258
2259 out.writeStop();
2260 return(out);
2261 }
2262 */
2263
2264
2265
2266
merge(SpatOptions & opt)2267 SpatRaster SpatRasterCollection::merge(SpatOptions &opt) {
2268 return mosaic("first", opt);
2269 }
2270
2271
2272
2273
mosaic(std::string fun,SpatOptions & opt)2274 SpatRaster SpatRasterCollection::mosaic(std::string fun, SpatOptions &opt) {
2275
2276 SpatRaster out;
2277
2278 std::vector<std::string> f {"first", "sum", "mean", "median", "min", "max"};
2279 if (std::find(f.begin(), f.end(), fun) == f.end()) {
2280 out.setError("not a valid function");
2281 return out;
2282 }
2283
2284 unsigned n = size();
2285
2286 if (n == 0) {
2287 out.setError("empty collection");
2288 return(out);
2289 }
2290 if (n == 1) {
2291 out = ds[0].deepCopy();
2292 return(out);
2293 }
2294
2295 std::vector<bool> hvals(n);
2296 hvals[0] = ds[0].hasValues();
2297 SpatExtent e = ds[0].getExtent();
2298 unsigned nl = ds[0].nlyr();
2299 std::vector<bool> resample(n, false);
2300 for (size_t i=1; i<n; i++) {
2301 // lyrs, crs, warncrs, ext, rowcol, res
2302 if (!ds[0].compare_geom(ds[i], false, false, opt.get_tolerance(), false, false, false, true)) {
2303 out.setError(ds[0].msg.error);
2304 return(out);
2305 }
2306 e.unite(ds[i].getExtent());
2307 hvals[i] = ds[i].hasValues();
2308 nl = std::max(nl, ds[i].nlyr());
2309 }
2310 out = ds[0].geometry(nl, false);
2311 out.setExtent(e, true, "");
2312
2313 for (int i=(n-1); i>=0; i--) {
2314 if (!hvals[i]) {
2315 erase(i);
2316 }
2317 }
2318
2319 n = size();
2320 if (size() == 0) {
2321 return out;
2322 }
2323
2324 SpatExtent eout = out.getExtent();
2325 double hyr = out.yres()/2;
2326
2327 std::string warn = "";
2328 for (size_t i=0; i<n; i++) {
2329 SpatOptions topt(opt);
2330 if(!ds[i].shared_basegeom(out, 0.1, true)) {
2331 SpatRaster temp = out.crop(ds[i].getExtent(), "near", topt);
2332 std::vector<bool> hascats = ds[i].hasCategories();
2333 std::string method = hascats[0] ? "near" : "bilinear";
2334 ds[i] = ds[i].warper(temp, "", method, false, false, topt);
2335 if (ds[i].hasError()) {
2336 out.setError(ds[i].getError());
2337 return out;
2338 }
2339 warn = "rasters did not align and were resampled";
2340 }
2341 }
2342 if (warn != "") out.addWarning(warn);
2343
2344 if (!out.writeStart(opt)) { return out; }
2345 SpatOptions sopt(opt);
2346 sopt.progressbar = false;
2347 std::vector<double> v;
2348 for (size_t i=0; i < out.bs.n; i++) {
2349 eout.ymax = out.yFromRow(out.bs.row[i]) + hyr;
2350 eout.ymin = out.yFromRow(out.bs.row[i] + out.bs.nrows[i] - 1) - hyr;
2351 SpatRasterStack s;
2352 for (size_t j=0; j<n; j++) {
2353 e = ds[j].getExtent();
2354 e.intersect(eout);
2355 if ( e.valid_notequal() ) {
2356 SpatRaster r = ds[j].crop(eout, "near", sopt);
2357 //SpatExtent ec = r.getExtent();
2358 r = r.extend(eout, sopt);
2359 //SpatExtent ee = r.getExtent();
2360 if (!s.push_back(r, "", "", "", false)) {
2361 out.setError("internal error: " + s.getError());
2362 out.writeStop();
2363 return out;
2364 }
2365 }
2366 }
2367 size_t ncls = out.bs.nrows[i] * out.ncol() * nl;
2368 if (s.size() > 0) {
2369 SpatRaster r = s.summary(fun, true, sopt);
2370 if (r.hasError()) {
2371 out.setError("internal error: " + r.getError());
2372 out.writeStop();
2373 return out;
2374 }
2375 if (!r.getValuesSource(0, v)) {
2376 out.setError("internal error: " + r.getError());
2377 out.writeStop();
2378 return out;
2379 }
2380 recycle(v, ncls);
2381 } else {
2382 v = std::vector<double>(ncls, NAN);
2383 }
2384 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, out.ncol())) return out;
2385 }
2386 out.writeStop();
2387 return(out);
2388 }
2389
2390
notisnan(const std::vector<double> & x,double & n)2391 void notisnan(const std::vector<double> &x, double &n) {
2392 for (size_t i=0; i<x.size(); i++) {
2393 n += !std::isnan(x[i]);
2394 }
2395 }
2396
2397
2398
do_stats(std::vector<double> & v,std::string fun,bool narm,double & stat,double & stat2,double & n,size_t step)2399 void do_stats(std::vector<double> &v, std::string fun, bool narm, double &stat, double &stat2,double &n, size_t step) {
2400 if (v.size() == 0) return;
2401 if (fun == "sum") {
2402 if (narm && (step > 0)) {
2403 v.push_back(stat);
2404 }
2405 stat = vsum(v, narm);
2406 } else if (fun == "mean") {
2407 if (narm) {
2408 notisnan(v, n);
2409 if (step > 0) {
2410 v.push_back(stat);
2411 }
2412 } else {
2413 n += v.size();
2414 }
2415 stat = vsum(v, narm);
2416 } else if (fun == "rms") {
2417 if (narm) {
2418 notisnan(v, n);
2419 } else {
2420 n += v.size();
2421 }
2422 double s = vsum2(v, narm);
2423 if (step > 1) {
2424 std::vector<double> ss = {stat, s};
2425 stat = vsum(ss, narm);
2426 } else {
2427 stat = s;
2428 }
2429 } else if (fun == "min") {
2430 double s = vmin(v, narm);
2431 if (step > 0) {
2432 std::vector<double> ss = {stat, s};
2433 stat = vmin(ss, narm);
2434 } else {
2435 stat = s;
2436 }
2437 } else if (fun == "max") {
2438 double s = vmax(v, narm);
2439 if (step > 0) {
2440 std::vector<double> ss = {stat, s};
2441 stat = vmax(ss, narm);
2442 } else {
2443 stat = s;
2444 }
2445 } else if (fun == "range") {
2446 double sn = vmin(v, narm);
2447 double sx = vmax(v, narm);
2448 if (step > 0) {
2449 std::vector<double> ss1 = {stat, sn};
2450 stat = vmin(ss1, narm);
2451 std::vector<double> ss2 = {stat2, sx};
2452 stat2 = vmax(ss2, narm);
2453 } else {
2454 stat = sn;
2455 stat2 = sx;
2456 }
2457 } else if (fun == "sd") {
2458 if (narm) {
2459 notisnan(v, n);
2460 } else {
2461 n += v.size();
2462 }
2463 double s1 = vsum(v, narm);
2464 double s2 = vsum2(v, narm);
2465 if (step > 1) {
2466 std::vector<double> ss1 = {stat, s1};
2467 stat = vsum(ss1, narm);
2468 std::vector<double> ss2 = {stat2, s2};
2469 stat2 = vsum(ss2, narm);
2470 } else {
2471 stat = s1;
2472 stat2 = s2;
2473 }
2474 }
2475 }
2476
2477
global(std::string fun,bool narm,SpatOptions & opt)2478 SpatDataFrame SpatRaster::global(std::string fun, bool narm, SpatOptions &opt) {
2479
2480 SpatDataFrame out;
2481 std::vector<std::string> f {"sum", "mean", "min", "max", "range", "rms", "sd", "sdpop"};
2482 if (std::find(f.begin(), f.end(), fun) == f.end()) {
2483 out.setError("not a valid function");
2484 return(out);
2485 }
2486
2487 if (!hasValues()) {
2488 out.setError("SpatRaster has no values");
2489 return(out);
2490 }
2491
2492 std::string sdfun = fun;
2493 if ((fun == "sdpop") || (fun == "sd")) {
2494 fun = "sd";
2495 }
2496 std::vector<double> stats(nlyr());
2497 std::vector<double> stats2(nlyr());
2498
2499 std::vector<double> n(nlyr());
2500 if (!readStart()) {
2501 out.setError(getError());
2502 return(out);
2503 }
2504 BlockSize bs = getBlockSize(opt);
2505 for (size_t i=0; i<bs.n; i++) {
2506 std::vector<double> v = readValues(bs.row[i], bs.nrows[i], 0, ncol());
2507 unsigned off = bs.nrows[i] * ncol() ;
2508 for (size_t lyr=0; lyr<nlyr(); lyr++) {
2509 unsigned offset = lyr * off;
2510 std::vector<double> vv = { v.begin()+offset, v.begin()+offset+off };
2511 do_stats(vv, fun, narm, stats[lyr], stats2[lyr], n[lyr], i);
2512 }
2513 }
2514 readStop();
2515
2516
2517 if (fun=="mean") {
2518 for (size_t lyr=0; lyr<nlyr(); lyr++) {
2519 if (n[lyr] > 0) {
2520 stats[lyr] = stats[lyr] / n[lyr];
2521 } else {
2522 stats[lyr] = NAN;
2523 }
2524 }
2525 } else if (fun=="rms") {
2526 // rms = sqrt(sum(x^2)/(n-1))
2527 for (size_t lyr=0; lyr<nlyr(); lyr++) {
2528 if (n[lyr] > 0) {
2529 stats[lyr] = sqrt(stats[lyr] / (n[lyr]-1));
2530 } else {
2531 stats[lyr] = NAN;
2532 }
2533 }
2534 } else if (fun == "sd") {
2535 for (size_t lyr=0; lyr<nlyr(); lyr++) {
2536 if (n[lyr] > 0) {
2537 double mn = stats[lyr] / n[lyr];
2538 double mnsq = mn * mn;
2539 double mnsumsq = stats2[lyr] / n[lyr];
2540 if (sdfun == "sdpop") {
2541 stats[lyr] = sqrt(mnsumsq - mnsq);
2542 } else {
2543 stats[lyr] = sqrt((mnsumsq - mnsq) * n[lyr]/(n[lyr]-1));
2544 }
2545
2546 } else {
2547 stats[lyr] = NAN;
2548 }
2549 }
2550 }
2551 out.add_column(stats, fun);
2552 if (fun=="range") {
2553 out.add_column(stats2, "max");
2554 }
2555 return(out);
2556 }
2557
2558
2559
global_weighted_mean(SpatRaster & weights,std::string fun,bool narm,SpatOptions & opt)2560 SpatDataFrame SpatRaster::global_weighted_mean(SpatRaster &weights, std::string fun, bool narm, SpatOptions &opt) {
2561
2562 SpatDataFrame out;
2563
2564 std::vector<std::string> f {"sum", "mean"};
2565 if (std::find(f.begin(), f.end(), fun) == f.end()) {
2566 out.setError("not a valid function");
2567 return(out);
2568 }
2569
2570 if (!hasValues()) {
2571 out.setError("SpatRaster has no values");
2572 return(out);
2573 }
2574
2575 if (weights.nlyr() != 1) {
2576 out.setError("The weights raster must have 1 layer");
2577 return(out);
2578 }
2579 if (!compare_geom(weights, false, false, opt.get_tolerance(), true)) {
2580 out.setError( msg.getError() );
2581 return(out);
2582 }
2583
2584 std::vector<double> stats(nlyr());
2585 double stats2 = 0;
2586 std::vector<double> n(nlyr());
2587 std::vector<double> w(nlyr());
2588 if (!readStart()) {
2589 out.setError(getError());
2590 return(out);
2591 }
2592 if (!weights.readStart()) {
2593 out.setError(weights.getError());
2594 return(out);
2595 }
2596
2597 BlockSize bs = getBlockSize(opt);
2598 for (size_t i=0; i<bs.n; i++) {
2599 std::vector<double> v = readValues(bs.row[i], bs.nrows[i], 0, ncol());
2600 std::vector<double> wv = weights.readValues(bs.row[i], bs.nrows[i], 0, ncol());
2601
2602 unsigned off = bs.nrows[i] * ncol() ;
2603 for (size_t lyr=0; lyr<nlyr(); lyr++) {
2604 double wsum = 0;
2605 unsigned offset = lyr * off;
2606 std::vector<double> vv(v.begin()+offset, v.begin()+offset+off);
2607 for (size_t j=0; j<vv.size(); j++) {
2608 if (!std::isnan(vv[j]) && !std::isnan(wv[j])) {
2609 vv[j] *= wv[j];
2610 wsum += wv[j];
2611 } else {
2612 vv[j] = NAN;
2613 }
2614 }
2615 do_stats(vv, fun, narm, stats[lyr], stats2, n[lyr], i);
2616 w[lyr] += wsum;
2617 }
2618 }
2619 readStop();
2620 weights.readStop();
2621
2622 if (fun=="mean") {
2623 for (size_t lyr=0; lyr<nlyr(); lyr++) {
2624 if (n[lyr] > 0 && w[lyr] != 0) {
2625 stats[lyr] /= w[lyr];
2626 } else {
2627 stats[lyr] = NAN;
2628 }
2629 }
2630 out.add_column(stats, "weighted_mean");
2631 } else {
2632 out.add_column(stats, "weighted_sum");
2633 }
2634
2635 return(out);
2636 }
2637
2638
scale(std::vector<double> center,bool docenter,std::vector<double> scale,bool doscale,SpatOptions & opt)2639 SpatRaster SpatRaster::scale(std::vector<double> center, bool docenter, std::vector<double> scale, bool doscale, SpatOptions &opt) {
2640 SpatRaster out;
2641 SpatOptions opts(opt);
2642 SpatDataFrame df;
2643 if (docenter) {
2644 if (center.size() == 0) {
2645 df = global("mean", true, opts);
2646 center = df.getD(0);
2647 }
2648 if (doscale) {
2649 out = arith(center, "-", false, opts);
2650 } else {
2651 out = arith(center, "-", false, opt);
2652 }
2653 }
2654 if (doscale) {
2655 if (scale.size() == 0) {
2656 // divide by sd if centered, and the root mean square otherwise.
2657 // rms = sqrt(sum(x^2)/(n-1)); if centered rms == sd
2658 if (docenter) {
2659 df = out.global("rms", true, opts);
2660 } else {
2661 df = global("rms", true, opts);
2662 }
2663 scale = df.getD(0);
2664 }
2665 if (docenter) {
2666 out = out.arith(scale, "/", false, opt);
2667 } else {
2668 out = arith(scale, "/", false, opt);
2669 }
2670 }
2671 return out;
2672 }
2673
2674
reclass_vector(std::vector<double> & v,std::vector<std::vector<double>> rcl,unsigned doright,bool lowest,bool othNA)2675 void reclass_vector(std::vector<double> &v, std::vector<std::vector<double>> rcl, unsigned doright, bool lowest, bool othNA) {
2676
2677 size_t nc = rcl.size(); // should be 2 or 3
2678 if (nc == 2) {
2679 doright = 3; // should be 2?
2680 }
2681 bool right = false;
2682 bool leftright = false;
2683 if (doright > 1) {
2684 leftright = true;
2685 } else if (doright) {
2686 right = true;
2687 }
2688
2689 // bool hasNA = false;
2690 double NAval = NAN;
2691
2692 size_t n = v.size();
2693 unsigned nr = rcl[0].size();
2694
2695 if (nc == 1) {
2696 std::vector<double> rc = rcl[0];
2697 std::sort(rc.begin(), rc.end());
2698 if (right) { // interval closed at left and right
2699 if (lowest) {
2700 for (size_t i=0; i<n; i++) {
2701 if (std::isnan(v[i])) {
2702 v[i] = NAval;
2703 } else if ((v[i] < rc[0]) | (v[i] > rc[nr-1])) {
2704 v[i] = NAval;
2705 } else {
2706 for (size_t j=1; j<nr; j++) {
2707 if (v[i] <= rc[j]) {
2708 v[i] = j-1;
2709 break;
2710 }
2711 }
2712 }
2713 }
2714 } else {
2715 for (size_t i=0; i<n; i++) {
2716 if (std::isnan(v[i])) {
2717 v[i] = NAval;
2718 } else if ((v[i] <= rc[0]) | (v[i] > rc[nr-1])) {
2719 v[i] = NAval;
2720 } else {
2721 for (size_t j=1; j<nr; j++) {
2722 if (v[i] <= rc[j]) {
2723 v[i] = j-1;
2724 break;
2725 }
2726 }
2727 }
2728 }
2729 }
2730 } else {
2731 if (lowest) {
2732 for (size_t i=0; i<n; i++) {
2733 if (std::isnan(v[i])) {
2734 v[i] = NAval;
2735 } else if ((v[i] < rc[0]) | (v[i] > rc[nr-1])) {
2736 v[i] = NAval;
2737 } else if (v[i] == rc[nr-1]) {
2738 v[i] = nr-1;
2739 } else {
2740 for (size_t j=1; j<nr; j++) {
2741 if (v[i] < rc[j]) {
2742 v[i] = j-1;
2743 break;
2744 }
2745 }
2746 }
2747 }
2748 } else {
2749 for (size_t i=0; i<n; i++) {
2750 if (std::isnan(v[i])) {
2751 v[i] = NAval;
2752 } else if ((v[i] < rc[0]) | (v[i] >= rc[nr-1])) {
2753 v[i] = NAval;
2754 } else {
2755 for (size_t j=1; j<nr; j++) {
2756 if (v[i] < rc[j]) {
2757 v[i] = j-1;
2758 break;
2759 }
2760 }
2761 }
2762 }
2763 }
2764 }
2765
2766 // "is - becomes"
2767 } else if (nc == 2) {
2768
2769 bool hasNAN = false;
2770 double replaceNAN = NAval;
2771 for (size_t j=0; j<nr; j++) {
2772 if (std::isnan(rcl[0][j])) {
2773 hasNAN = true;
2774 replaceNAN = rcl[1][j];
2775 }
2776 }
2777 for (size_t i=0; i<n; i++) {
2778 if (std::isnan(v[i])) {
2779 if (hasNAN) {
2780 v[i] = replaceNAN;
2781 } else {
2782 v[i] = NAval;
2783 }
2784 } else {
2785 bool found = false;
2786 for (size_t j=0; j<nr; j++) {
2787 if (v[i] == rcl[0][j]) {
2788 v[i] = rcl[1][j];
2789 found = true;
2790 break;
2791 }
2792 }
2793 if ((othNA) & (!found)) {
2794 v[i] = NAval;
2795 }
2796 }
2797 }
2798
2799 // "from - to - becomes"
2800 } else {
2801
2802 bool hasNAN = false;
2803 double replaceNAN = NAval;
2804 for (size_t j=0; j<nr; j++) {
2805 if (std::isnan(rcl[0][j]) || std::isnan(rcl[1][j])) {
2806 hasNAN = true;
2807 replaceNAN = rcl[2][j];
2808 }
2809 }
2810
2811 if (leftright) { // interval closed at left and right
2812
2813 for (size_t i=0; i<n; i++) {
2814 if (std::isnan(v[i])) {
2815 if (hasNAN) {
2816 v[i] = replaceNAN;
2817 } else {
2818 v[i] = NAval;
2819 }
2820 } else {
2821 bool found = false;
2822 for (size_t j=0; j<nr; j++) {
2823 if ((v[i] >= rcl[0][j]) & (v[i] <= rcl[1][j])) {
2824 v[i] = rcl[2][j];
2825 found = true;
2826 break;
2827 }
2828 }
2829 if ((othNA) & (!found)) {
2830 v[i] = NAval;
2831 }
2832 }
2833 }
2834 } else if (right) { // interval closed at right
2835 if (lowest) { // include lowest value (left) of interval
2836
2837 double lowval = rcl[0][0];
2838 double lowres = rcl[2][0];
2839 for (size_t i=1; i<nr; i++) {
2840 if (rcl[0][i] < lowval) {
2841 lowval = rcl[0][i];
2842 lowres = rcl[2][i];
2843 }
2844 }
2845
2846 for (size_t i=0; i<n; i++) {
2847 if (std::isnan(v[i])) {
2848 if (hasNAN) {
2849 v[i] = replaceNAN;
2850 } else {
2851 v[i] = NAval;
2852 }
2853 } else if (v[i] == lowval) {
2854 v[i] = lowres;
2855 } else {
2856 bool found = false;
2857 for (size_t j=0; j<nr; j++) {
2858 if ((v[i] > rcl[0][j]) & (v[i] <= rcl[1][j])) {
2859 v[i] = rcl[2][j];
2860 found = true;
2861 break;
2862 }
2863 }
2864 if ((othNA) & (!found)) {
2865 v[i] = NAval;
2866 }
2867 }
2868 }
2869
2870 } else { // !dolowest
2871 for (size_t i=0; i<n; i++) {
2872 if (std::isnan(v[i])) {
2873 if (hasNAN) {
2874 v[i] = replaceNAN;
2875 } else {
2876 v[i] = NAval;
2877 }
2878 } else {
2879 bool found = false;
2880 for (size_t j=0; j<nr; j++) {
2881 if ((v[i] > rcl[0][j]) & (v[i] <= rcl[1][j])) {
2882 v[i] = rcl[2][j];
2883 found = true;
2884 break;
2885 }
2886 }
2887 if ((othNA) & (!found)) {
2888 v[i] = NAval;
2889 }
2890 }
2891 }
2892 }
2893
2894 } else { // !doright
2895
2896 if (lowest) { // which here means highest because right=FALSE
2897
2898 double lowval = rcl[1][0];
2899 double lowres = rcl[2][0];
2900 for (size_t i=0; i<nr; i++) {
2901 if (rcl[1][i] > lowval) {
2902 lowval = rcl[1][i];
2903 lowres = rcl[2][i];
2904 }
2905 }
2906
2907 for (size_t i=0; i<n; i++) {
2908 if (std::isnan(v[i])) {
2909 if (hasNAN) {
2910 v[i] = replaceNAN;
2911 } else {
2912 v[i] = NAval;
2913 }
2914 } else if (v[i] == lowval) {
2915 v[i] = lowres;
2916 } else {
2917 bool found = false;
2918 for (size_t j=0; j<nr; j++) {
2919 if ((v[i] >= rcl[0][j]) & (v[i] < rcl[1][j])) {
2920 v[i] = rcl[2][j];
2921 found = true;
2922 break;
2923 }
2924 }
2925 if ((othNA) & (!found)) {
2926 v[i] = NAval;
2927 }
2928 }
2929 }
2930
2931 } else { //!dolowest
2932
2933 for (size_t i=0; i<n; i++) {
2934 if (std::isnan(v[i])) {
2935 if (hasNAN) {
2936 v[i] = replaceNAN;
2937 } else {
2938 v[i] = NAval;
2939 }
2940 } else {
2941 bool found = false;
2942 for (size_t j=0; j<nr; j++) {
2943 if ((v[i] >= rcl[0][j]) & (v[i] < rcl[1][j])) {
2944 v[i] = rcl[2][j];
2945 found = true;
2946 break;
2947 }
2948 }
2949 if ((othNA) & (!found)) {
2950 v[i] = NAval;
2951 }
2952 }
2953 }
2954 }
2955 }
2956 }
2957 }
2958
2959
replaceValues(std::vector<double> from,std::vector<double> to,long nl,SpatOptions & opt)2960 SpatRaster SpatRaster::replaceValues(std::vector<double> from, std::vector<double> to, long nl, SpatOptions &opt) {
2961
2962 SpatRaster out = geometry(nl);
2963 bool multi = false;
2964 if (nl > 1) {
2965 if (nlyr() > 1) {
2966 out.setError("cannot create layer-varying replacement with multi-layer input");
2967 return out;
2968 }
2969 multi = true;
2970 }
2971
2972 if (!readStart()) {
2973 out.setError(getError());
2974 return(out);
2975 }
2976 if (!out.writeStart(opt)) {
2977 readStop();
2978 return out;
2979 }
2980
2981 if (multi) {
2982 size_t tosz = to.size() / nl;
2983 size_t nlyr = out.nlyr();
2984 for (size_t i = 0; i < out.bs.n; i++) {
2985 std::vector<double> v = readBlock(out.bs, i);
2986 size_t vs = v.size();
2987 v.reserve(vs * nlyr);
2988 for (size_t lyr = 1; lyr < nlyr; lyr++) {
2989 v.insert(v.end(), v.begin(), v.begin()+vs);
2990 }
2991 for (size_t lyr = 0; lyr < nlyr; lyr++) {
2992 std::vector<double> tolyr(to.begin()+lyr*tosz, to.begin()+(lyr+1)*tosz);
2993 recycle(tolyr, from);
2994 size_t offset = lyr*vs;
2995 for (size_t j=0; j< from.size(); j++) {
2996 if (std::isnan(from[j])) {
2997 for (size_t k=offset; k<(offset+vs); k++) {
2998 v[k] = std::isnan(v[k]) ? tolyr[j] : v[k];
2999 }
3000 } else {
3001 std::replace(v.begin()+offset, v.begin()+(offset+vs), from[j], tolyr[j]);
3002 }
3003 }
3004 }
3005 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
3006 }
3007 } else {
3008 recycle(to, from);
3009 for (size_t i = 0; i < out.bs.n; i++) {
3010 std::vector<double> v = readBlock(out.bs, i);
3011 for (size_t j=0; j< from.size(); j++) {
3012 if (std::isnan(from[j])) {
3013 for (double &d : v) d = std::isnan(d) ? to[j] : d;
3014 } else {
3015 std::replace(v.begin(), v.end(), from[j], to[j]);
3016 }
3017 }
3018 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
3019 }
3020 }
3021 readStop();
3022 out.writeStop();
3023 return(out);
3024 }
3025
3026
3027
reclassify(std::vector<std::vector<double>> rcl,unsigned right,bool lowest,bool othersNA,bool bylayer,SpatOptions & opt)3028 SpatRaster SpatRaster::reclassify(std::vector<std::vector<double>> rcl, unsigned right, bool lowest, bool othersNA, bool bylayer, SpatOptions &opt) {
3029
3030 SpatRaster out = geometry();
3031 size_t nc = rcl.size();
3032 size_t nr = rcl[0].size();
3033 size_t nl = nlyr();
3034 if (nl == 1) bylayer = false;
3035 size_t maxnc = 3 + nl * bylayer;
3036 size_t rcldim = nc;
3037
3038 if (bylayer) {
3039 if (((nc != maxnc) && (nc != (maxnc-1))) || nr < 1) {
3040 out.setError("reclass matrix is not correct. Should be nlyr(x) plus 1 or 2");
3041 return out;
3042 }
3043 rcldim = nc - (nl-1);
3044 } else {
3045 if (nc < 1 || nc > 3 || nr < 1) {
3046 out.setError("matrix must have 1, 2 or 3 columns, and at least one row");
3047 return out;
3048 }
3049 }
3050
3051 if (nc == 1) {
3052 if (nr == 1) {
3053 int breaks = rcl[0][0];
3054 if (breaks < 2) {
3055 out.setError("cannot classify with a single number that is smaller than 2");
3056 return out;
3057 }
3058 std::vector<bool> hr = hasRange();
3059 bool hasR = true;
3060 for (size_t i=0; i<hr.size(); i++) {
3061 if (!hr[i]) hasR = false;
3062 }
3063
3064 if (!hasR) {
3065 SpatOptions xopt(opt);
3066 setRange(xopt);
3067 }
3068 std::vector<double> mn = range_min();
3069 std::vector<double> mx = range_max();
3070 double mnv = vmin(mn, true);
3071 double mxv = vmax(mx, true);
3072 rcl[0] = seq_steps(mnv, mxv, breaks);
3073 }
3074
3075 if (rcl[0].size() < 256) {
3076 std::vector<std::string> s;
3077 for (size_t i=1; i<rcl[0].size(); i++) {
3078 s.push_back(double_to_string(rcl[0][i-1]) + " - " + double_to_string(rcl[0][i]));
3079 }
3080 for (size_t i=0; i<out.nlyr(); i++) {
3081 out.setLabels(i, s);
3082 }
3083 }
3084 nr = rcl[0].size();
3085 }
3086 for (size_t i=0; i<nc; i++) {
3087 if (rcl[i].size() != nr) {
3088 out.setError("matrix is not rectangular");
3089 return out;
3090 }
3091 }
3092 if (rcldim == 3) {
3093 for (size_t i=0; i<nr; i++) {
3094 if (rcl[0][i] > rcl[1][i]) {
3095 out.setError("'from' larger than 'to': (" + std::to_string(rcl[0][i]) + " - " + std::to_string(rcl[1][i]) +")");
3096 return out;
3097 }
3098 }
3099 }
3100
3101 if (!readStart()) {
3102 out.setError(getError());
3103 return(out);
3104 }
3105
3106 if (!out.writeStart(opt)) {
3107 readStop();
3108 return out;
3109 }
3110
3111 if (bylayer) {
3112 std::vector<std::vector<double>> lyrrcl(rcldim+1);
3113 for (size_t i=0; i<rcldim; i++) {
3114 lyrrcl[i] = rcl[i];
3115 }
3116 for (size_t i = 0; i < out.bs.n; i++) {
3117 unsigned off = bs.nrows[i] * ncol() ;
3118 std::vector<double> v = readBlock(out.bs, i);
3119 for (size_t lyr = 0; lyr < nl; lyr++) {
3120 unsigned offset = lyr * off;
3121 lyrrcl[rcldim] = rcl[rcldim+lyr];
3122 std::vector<double> vx(v.begin()+offset, v.begin()+offset+off);
3123 reclass_vector(vx, lyrrcl, right, lowest, othersNA);
3124 std::copy(vx.begin(), vx.end(), v.begin()+offset);
3125 }
3126 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
3127 }
3128 } else {
3129 for (size_t i = 0; i < out.bs.n; i++) {
3130 std::vector<double> v = readBlock(out.bs, i);
3131 reclass_vector(v, rcl, right, lowest, othersNA);
3132 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) return out;
3133 }
3134 }
3135
3136 readStop();
3137 out.writeStop();
3138 return(out);
3139
3140 }
3141
3142
reclassify(std::vector<double> rcl,unsigned nc,unsigned right,bool lowest,bool othersNA,bool bylayer,SpatOptions & opt)3143 SpatRaster SpatRaster::reclassify(std::vector<double> rcl, unsigned nc, unsigned right, bool lowest, bool othersNA, bool bylayer, SpatOptions &opt) {
3144
3145 SpatRaster out;
3146 if ((rcl.size() % nc) != 0) {
3147 out.setError("incorrect length of reclassify matrix");
3148 return(out);
3149 }
3150 size_t maxnc = 3 + bylayer * (nlyr() - 1);
3151 unsigned nr = rcl.size() / nc;
3152 if (nc > maxnc) {
3153 out.setError("incorrect number of columns in reclassify matrix");
3154 return(out);
3155 }
3156 std::vector< std::vector<double>> rc(nc);
3157
3158 for (size_t i=0; i<nc; i++) {
3159 rc[i] = std::vector<double>(rcl.begin()+(i*nr), rcl.begin()+(i+1)*nr);
3160 }
3161
3162 out = reclassify(rc, right, lowest, othersNA, bylayer, opt);
3163 return out;
3164 }
3165
3166
3167
clump_getRCL(std::vector<std::vector<size_t>> rcl,size_t n)3168 std::vector<std::vector<double>> clump_getRCL(std::vector<std::vector<size_t>> rcl, size_t n) {
3169 std::vector<std::vector<size_t>> rcl2(rcl[0].size());
3170 for (size_t i=0; i<rcl[0].size(); i++) {
3171 rcl2[i].push_back(rcl[0][i]);
3172 rcl2[i].push_back(rcl[1][i]);
3173 }
3174 std::sort(rcl2.begin(), rcl2.end());
3175 rcl2.erase(std::unique(rcl2.begin(), rcl2.end()), rcl2.end());
3176 std::vector<std::vector<double>> out(2);
3177 for (size_t i=0; i<rcl2.size(); i++) {
3178 out[0].push_back(rcl2[i][1]);
3179 out[1].push_back(rcl2[i][0]);
3180 }
3181 // from - to
3182 // 3 - 1
3183 // 4 - 3
3184 // becomes
3185 // 3 - 1
3186 // 4 - 1
3187 for (size_t i=1; i<out[0].size(); i++) {
3188 for (size_t j=0; j<i; j++) {
3189 if (out[0][i] == out[1][j]) {
3190 out[1][j] = out[0][i];
3191 }
3192 }
3193 }
3194
3195 std::vector<double> lost = out[0];
3196 lost.push_back(n);
3197 size_t sub = 0;
3198 for (size_t i=0; i<lost.size(); i++) {
3199 sub++;
3200 for (size_t j=lost[i]+1; j<lost[i+1]; j++) {
3201 out[0].push_back(j);
3202 out[1].push_back(j-sub);
3203 }
3204 }
3205 return out;
3206 }
3207
3208
clump_replace(std::vector<double> & v,size_t n,const std::vector<double> & d,size_t cstart,std::vector<std::vector<size_t>> & rcl)3209 void clump_replace(std::vector<double> &v, size_t n, const std::vector<double>& d, size_t cstart, std::vector<std::vector<size_t>>& rcl) {
3210 for (size_t i=0; i<n; i++) {
3211 for (size_t j=1; j<d.size(); j++) {
3212 if (v[i] == d[j]) {
3213 v[i] = d[0];
3214 }
3215 }
3216 }
3217 if (d[0] < cstart) {
3218 for (size_t j=1; j<d.size(); j++) {
3219 rcl[0].push_back(d[0]);
3220 rcl[1].push_back(d[j]);
3221 }
3222 }
3223 }
3224
3225
clump_test(std::vector<double> & d)3226 void clump_test(std::vector<double> &d) {
3227 d.erase(std::remove_if(d.begin(), d.end(),
3228 [](const double& v) { return std::isnan(v); }), d.end());
3229 std::sort(d.begin(), d.end());
3230 d.erase(std::unique(d.begin(), d.end()), d.end());
3231 }
3232
broom_clumps(std::vector<double> & v,std::vector<double> & above,const size_t & dirs,size_t & ncps,const size_t & nr,const size_t & nc,std::vector<std::vector<size_t>> & rcl)3233 void broom_clumps(std::vector<double> &v, std::vector<double>& above, const size_t &dirs, size_t &ncps, const size_t &nr, const size_t &nc, std::vector<std::vector<size_t>> &rcl) {
3234
3235 size_t nstart = ncps;
3236
3237 bool d4 = dirs == 4;
3238
3239 if ( !std::isnan(v[0]) ) { //first cell, no cell left of it
3240 if (std::isnan(above[0])) {
3241 v[0] = ncps;
3242 ncps++;
3243 } else {
3244 v[0] = above[0];
3245 }
3246 }
3247
3248 for (size_t i=1; i<nc; i++) { //first row, no row above it, use "above"
3249 if (!std::isnan(v[i])) {
3250 std::vector<double> d;
3251 if (d4) {
3252 d = {above[i], v[i-1]} ;
3253 } else {
3254 d = {above[i], above[i-1], v[i-1]} ;
3255 }
3256 clump_test(d);
3257 if (d.size() > 0) {
3258 v[i] = d[0];
3259 if (d.size() > 1) {
3260 clump_replace(v, i, d, nstart, rcl);
3261 }
3262 } else {
3263 v[i] = ncps;
3264 ncps++;
3265 }
3266 }
3267 }
3268
3269
3270 for (size_t r=1; r<nr; r++) { //other rows
3271 size_t i=r*nc;
3272 if (!std::isnan(v[i])) { // first cell
3273 if (std::isnan(v[i-nc])) {
3274 v[i] = ncps;
3275 ncps++;
3276 } else {
3277 v[i] = v[i-nc];
3278 }
3279 }
3280 for (size_t i=r*nc+1; i<((r+1)*nc); i++) { // other cells
3281 if (!std::isnan(v[i])) {
3282 std::vector<double> d;
3283 if (d4) {
3284 d = {v[i-nc], v[i-1]} ;
3285 } else {
3286 d = {v[i-nc], v[i-nc-1], v[i-1]} ;
3287 }
3288 clump_test(d);
3289 if (d.size() > 0) {
3290 v[i] = d[0];
3291 if (d.size() > 1) {
3292 clump_replace(v, i, d, nstart, rcl);
3293 }
3294 } else {
3295 v[i] = ncps;
3296 ncps++;
3297 }
3298 }
3299 }
3300 }
3301 size_t off = (nr-1) * nc;
3302 above = std::vector<double>(v.begin()+off, v.end());
3303 }
3304
3305
3306
clumps(int directions,bool zeroAsNA,SpatOptions & opt)3307 SpatRaster SpatRaster::clumps(int directions, bool zeroAsNA, SpatOptions &opt) {
3308
3309 SpatRaster out = geometry(1);
3310 if (nlyr() > 1) {
3311 SpatOptions ops(opt);
3312 std::string filename = opt.get_filename();
3313 ops.set_filenames({""});
3314 for (size_t i=0; i<nlyr(); i++) {
3315 std::vector<unsigned> lyr = {(unsigned)i};
3316 SpatRaster x = subset(lyr, ops);
3317 x = x.clumps(directions, zeroAsNA, ops);
3318 out.addSource(x);
3319 }
3320 if (filename != "") {
3321 out = out.writeRaster(opt);
3322 }
3323 return out;
3324 }
3325
3326 if (!(directions == 4 || directions == 8)) {
3327 out.setError("directions must be 4 or 8");
3328 return out;
3329 }
3330 if (!hasValues()) {
3331 out.setError("cannot compute clumps for a raster with no values");
3332 return out;
3333 }
3334
3335 std::vector<size_t> dim = {nrow(), ncol()};
3336
3337 std::string tempfile = "";
3338 std::vector<double> d, v, vv;
3339 if (!readStart()) {
3340 out.setError(getError());
3341 return(out);
3342 }
3343 std::string filename = opt.get_filename();
3344 if (filename != "") {
3345 bool overwrite = opt.get_overwrite();
3346 std::string errmsg;
3347 if (!can_write(filename, overwrite, errmsg)) {
3348 out.setError(errmsg + " (" + filename +")");
3349 return(out);
3350 }
3351 }
3352
3353 opt.set_filenames({""});
3354 if (!out.writeStart(opt)) { return out; }
3355 size_t nc = ncol();
3356 size_t ncps = 1;
3357 std::vector<double> above(nc, NAN);
3358 std::vector<std::vector<size_t>> rcl(2);
3359 for (size_t i = 0; i < out.bs.n; i++) {
3360 v = readBlock(out.bs, i);
3361 if (zeroAsNA) {
3362 std::replace(v.begin(), v.end(), 0.0, (double)NAN);
3363 }
3364 broom_clumps(v, above, directions, ncps, out.bs.nrows[i], nc, rcl);
3365 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, nc)) return out;
3366 }
3367 out.writeStop();
3368 readStop();
3369
3370 opt.set_filenames({filename});
3371 if (rcl[0].size() > 0) {
3372 std::vector<std::vector<double>> rc = clump_getRCL(rcl, ncps);
3373 out = out.reclassify(rc, 3, true, false, false, opt);
3374 } else if (filename != "") {
3375 out = out.writeRaster(opt);
3376 }
3377 return out;
3378 }
3379
3380
3381
3382