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 "spatRaster.h"
20
21 /*
22 #include "string_utils.h"
23
24 void SpatRasterSource::fsopen(std::string filename) {
25 std::string grifile = setFileExt(filename, ".gri");
26 std::ofstream fstr(grifile, std::ios::out | std::ios::binary);
27 *ofs = &fstr;
28 }
29
30 bool SpatRasterSource::fswrite(std::vector<double> &v) {
31 unsigned sz = v.size() * sizeof(double);
32 bool result = (*ofs).write(reinterpret_cast<const char*>(&v[0]), sz);
33 return result;
34 }
35
36 void SpatRasterSource::fsclose() {
37 (*ofs).close();
38 }
39 */
40
41
SpatRasterSource()42 SpatRasterSource::SpatRasterSource() {
43 open_write = false;
44 open_read = false;
45 }
46
47
48
combineSources(SpatRaster x)49 SpatRaster SpatRaster::combineSources(SpatRaster x) {
50
51 SpatRaster out = geometry();
52
53 // opt.get_tolerance()
54 if (!out.compare_geom(x, false, false, 0.1)) {
55 return out;
56 }
57
58 bool hv = hasValues();
59 if (hv != x.hasValues()) {
60 out.setError("combined sources must all have values; or none should have values");
61 return(out);
62 }
63
64 out = deepCopy();
65 // if (!hv) {
66 // out.source = x.source;
67 // } else {
68 out.source.insert(out.source.end(), x.source.begin(), x.source.end());
69 // }
70 // to make names unique
71 out.setNames(out.getNames());
72 return(out);
73 }
74
75
combine(SpatRaster x)76 void SpatRaster::combine(SpatRaster x) {
77
78
79 if (!compare_geom(x, false, false, 0.1)) {
80 return;
81 }
82
83 bool hv = hasValues();
84 if (hv != x.hasValues()) {
85 setError("combined sources must all have values; or none should have values");
86 return;
87 }
88
89 source.insert(source.end(), x.source.begin(), x.source.end());
90 setNames(getNames());
91 return;
92 }
93
addSource(SpatRaster x)94 void SpatRaster::addSource(SpatRaster x) {
95
96 if (compare_geom(x, false, false, 0.1 )) {
97 if (!hasValues()) { //or if n src == 0?
98 source = x.source;
99 } else {
100 source.insert(source.end(), x.source.begin(), x.source.end());
101 }
102 }
103 }
104
105
nsrc()106 unsigned SpatRaster::nsrc() {
107 return source.size();
108 }
109
110
sourceFromLyr(unsigned lyr)111 int SpatRaster::sourceFromLyr(unsigned lyr) {
112 if (lyr >= nlyr()) {
113 return(-1);
114 }
115 unsigned nsrc = 0;
116 unsigned nlyrs = -1;
117 for (size_t i=0; i<source.size(); i++) {
118 nlyrs += source[i].nlyr;
119 if (nlyrs >= lyr) break;
120 nsrc++;
121 }
122 return nsrc;
123 }
124
125
nlyrBySource()126 std::vector<unsigned> SpatRaster::nlyrBySource() {
127 std::vector<unsigned> lyrs(source.size());
128 for (size_t i=0; i<source.size(); i++) {
129 lyrs[i] = source[i].nlyr;
130 }
131 return lyrs;
132 }
133
lyrsBySource()134 std::vector<unsigned> SpatRaster::lyrsBySource() {
135 std::vector<unsigned> lyrs(nlyr());
136 unsigned start = 0;
137 for (size_t i=0; i<source.size(); i++) {
138 unsigned n = source[i].nlyr;
139 unsigned stop = start + n;
140 for (size_t j=start; j<stop; j++) {
141 lyrs[j] = i;
142 }
143 start = stop;
144 }
145 return lyrs;
146 }
147
148
findLyr(unsigned lyr)149 std::vector<unsigned> SpatRaster::findLyr(unsigned lyr) {
150 std::vector<unsigned> sl(2);
151 unsigned nlyrs = 0;
152 unsigned start = 0;
153 bool done = false;
154 for (size_t i=0; i<source.size(); i++) {
155 if ((nlyrs + source[i].nlyr) >= lyr) {
156 sl[0] = i;
157 for (size_t j=start; j<source[i].nlyr; j++) {
158 if ((nlyrs + j) == lyr) {
159 sl[1] = j;
160 done = true;
161 break;
162 }
163 }
164 }
165 if (done) break;
166 nlyrs += source[i].nlyr;
167 }
168 return sl;
169 }
170
getBands()171 std::vector<unsigned> SpatRaster::getBands() {
172 std::vector<unsigned> out;
173 for (size_t i=0; i<source.size(); i++) {
174 out.insert(out.end(), source[i].layers.begin(), source[i].layers.end());
175 }
176 return out;
177 }
178
179
180
181
182
183
sourcesFromLyrs(std::vector<unsigned> lyrs)184 std::vector<unsigned> SpatRaster::sourcesFromLyrs(std::vector<unsigned> lyrs) {
185 std::vector<unsigned> s(lyrs.size());
186 std::vector<unsigned> slyrs = lyrsBySource();
187 for (size_t i=0; i<lyrs.size(); i++) {
188 s[i] = slyrs[ lyrs[i] ];
189 }
190 return s;
191 }
192
193
194 /*
195 void SpatRasterSource::getValues(std::vector<double> &v, unsigned lyr) {
196 size_t nc ;
197 if (hasWindow) {
198 nc = window.full_ncol * window.full_nrow;
199 } else {
200 nc = nrow * ncol;
201 }
202 size_t start = lyr * nc;
203 v = std::vector<double>(values.begin()+start, values.begin()+start+nc);
204 }
205 */
206
207
appendValues(std::vector<double> & v,unsigned lyr)208 void SpatRasterSource::appendValues(std::vector<double> &v, unsigned lyr) {
209 size_t nc ;
210 if (hasWindow) {
211 nc = window.full_ncol * window.full_nrow;
212 } else {
213 nc = nrow * ncol;
214 }
215 size_t start = lyr * nc;
216 v.insert(v.end(), values.begin()+start, values.begin()+start+nc);
217 }
218
219
in_order()220 bool SpatRasterSource::in_order() {
221 if (memory) return true;
222 if (nlyr != nlyrfile) return false;
223 for (size_t i=0; i<layers.size(); i++) {
224 if (layers[i] != i) {
225 return false;
226 }
227 }
228 return true;
229 }
230
231
resize(unsigned n)232 void SpatRasterSource::resize(unsigned n) {
233 names.resize(n, "");
234 time.resize(n);
235 unit.resize(n);
236 depth.resize(n);
237 hasRange.resize(n, false);
238 range_min.resize(n, NAN);
239 range_max.resize(n, NAN);
240 blockcols.resize(n);
241 blockrows.resize(n);
242
243 has_scale_offset.resize(n, false);
244 scale.resize(n, 1);
245 offset.resize(n, 0);
246 hasColors.resize(n, false);
247 cols.resize(n);
248 hasCategories.resize(n, false);
249 cats.resize(n);
250 //hasAttributes.resize(n);
251 //atts.resize(n);
252 //attsIndex.resize(n);
253 nlyr = n;
254 layers.resize(n);
255 std::iota(layers.begin(), layers.end(), 0);
256 }
257
258
259 //std::vector<SpatRasterSource> SpatRasterSource::subset(std::vector<unsigned> lyrs) {
subset(std::vector<unsigned> lyrs)260 SpatRasterSource SpatRasterSource::subset(std::vector<unsigned> lyrs) {
261
262 unsigned nl = lyrs.size();
263 bool all = true;
264 if (lyrs.size() == nlyr) {
265 for (size_t i=0; i<nl; i++) {
266 if (lyrs[i] != i) {
267 all = false;
268 break;
269 }
270 }
271 } else {
272 all = false;
273 }
274 if (all) {
275 return *this ;
276 }
277
278 SpatRasterSource out;
279 if (memory) {
280 out.srs = srs;
281 out.ncol = ncol;
282 out.nrow = nrow;
283 out.nlyrfile = nlyrfile;
284 out.extent = extent;
285 out.rotated = rotated;
286 out.flipped = flipped;
287 out.hasWindow = hasWindow;
288 out.window = window;
289 out.source_name = source_name;
290 out.source_name_long = source_name_long;
291 out.timestep = timestep;
292 out.hasTime = hasTime;
293 out.memory = memory;
294 out.hasValues = hasValues;
295 //out.filename = filename;
296 //out.driver = driver;
297 //out.datatype = datatype;
298 //out.hasNAflag = hasNAflag;
299 //out.NAflag = NAflag;
300 } else {
301 //no values, deep copy is cheap
302 out = *this;
303 out.resize(0);
304 }
305
306
307 for (size_t i=0; i<nl; i++) {
308 unsigned j = lyrs[i];
309 out.names.push_back(names[j]);
310 out.time.push_back(time[j]);
311 out.depth.push_back(depth[j]);
312 out.unit.push_back(unit[j]);
313 out.hasRange.push_back(hasRange[j]);
314 out.range_min.push_back(range_min[j]);
315 out.range_max.push_back(range_max[j]);
316 out.blockrows.push_back(blockrows[j]);
317 out.blockcols.push_back(blockcols[j]);
318 out.hasColors.push_back(hasColors[j]);
319 out.cols.push_back(cols[j]);
320 out.hasCategories.push_back(hasCategories[j]);
321 out.cats.push_back(cats[j]);
322 out.has_scale_offset.push_back(has_scale_offset[j]);
323 out.scale.push_back(scale[j]);
324 out.offset.push_back(offset[j]);
325
326 if (memory) {
327 out.layers.push_back(i);
328 if (hasValues) {
329 appendValues(out.values, j);
330 }
331 } else {
332 out.layers.push_back(layers[j]);
333 }
334 }
335 out.nlyr = nl;
336 out.hasValues = hasValues;
337 return out;
338 }
339
340
validLayers(std::vector<unsigned> lyrs,unsigned nl)341 std::vector<unsigned> validLayers( std::vector<unsigned> lyrs , unsigned nl) {
342 unsigned s = lyrs.size();
343 for (size_t i=0; i<s; i++) {
344 unsigned j = s - i - 1; // start from the back
345 if (lyrs[j] >= nl) {
346 lyrs.erase(lyrs.begin() + j);
347 }
348 }
349
350 /* or
351 unsigned s = lyrs.size() - 1;
352 for (long i=s; i>=0; i--) {
353 if ((lyrs[i] < 0) | (lyrs[i] >= nl)) {
354 lyrs.erase(lyrs.begin() + i);
355 }
356 }
357 */
358 return lyrs;
359 }
360
361
subset(std::vector<unsigned> lyrs,SpatOptions & opt)362 SpatRaster SpatRaster::subset(std::vector<unsigned> lyrs, SpatOptions &opt) {
363
364 SpatRaster out = geometry(1);
365 out.source.resize(0);
366
367 unsigned oldsize = lyrs.size();
368 lyrs = validLayers(lyrs, nlyr());
369
370 if (lyrs.size() == 0) {
371 out.setError("no (valid) layer references");
372 return(out);
373 } else if (lyrs.size() != oldsize) {
374 out.addWarning("ignored " + std::to_string(oldsize - lyrs.size()) + " invalid layer reference(s)");
375 }
376
377 std::vector<unsigned> srcs = sourcesFromLyrs(lyrs);
378 unsigned ss = srcs[0];
379 std::vector<unsigned> slyr;
380 std::vector<unsigned> lyrbys = nlyrBySource();
381 // SpatRasterSource rs;
382 unsigned offset = 0;
383 for (size_t i=0; i<ss; i++) { offset += lyrbys[i]; }
384
385 for (size_t i=0; i<lyrs.size(); i++) {
386 if (srcs[i] == ss) {
387 slyr.push_back( (lyrs[i] - offset) );
388 } else {
389 out.source.push_back( source[ss].subset(slyr) );
390 ss = srcs[i];
391 offset = 0;
392 for (size_t i=0; i<ss; i++) { offset += lyrbys[i]; }
393 slyr = { lyrs[i] - offset } ;
394 }
395 }
396
397 out.source.push_back( source[ss].subset(slyr) );
398 if (opt.get_filename() != "") {
399 out = out.writeRaster(opt);
400 } //else {
401 // out.collapse();
402 //}
403
404 return out;
405 }
406
407
408
409
combine_sources(const SpatRasterSource & x)410 bool SpatRasterSource::combine_sources(const SpatRasterSource &x) {
411 if (memory & x.memory) {
412 if ((values.size() + x.values.size()) < (values.max_size()/8) ) {
413 values.insert(values.end(), x.values.begin(), x.values.end());
414 layers.resize(nlyr + x.nlyr);
415 std::iota(layers.begin(), layers.end(), 0);
416 } else {
417 return false;
418 }
419 } else if (filename == x.filename) {
420 layers.insert(layers.end(), x.layers.begin(), x.layers.end());
421 } else {
422 return false;
423 }
424 nlyr += x.nlyr;
425 names.insert(names.end(), x.names.begin(), x.names.end());
426 time.insert(time.end(), x.time.begin(), x.time.end());
427 if (!(hasTime & x.hasTime)) {
428 hasTime = false;
429 }
430 unit.insert(unit.end(), x.unit.begin(), x.unit.end());
431
432 depth.insert(depth.end(), x.depth.begin(), x.depth.end());
433 hasRange.insert(hasRange.end(), x.hasRange.begin(), x.hasRange.end());
434 range_min.insert(range_min.end(), x.range_min.begin(), x.range_min.end());
435 range_max.insert(range_max.end(), x.range_max.begin(), x.range_max.end());
436
437 blockrows.insert(blockrows.end(), x.blockrows.begin(), x.blockrows.end());
438 blockcols.insert(blockcols.end(), x.blockcols.begin(), x.blockcols.end());
439 //hasAttributes.insert(hasAttributes.end(), x.hasAttributes.begin(), x.hasAttributes.end());
440 //atts.insert(atts.end(), x.atts.begin(), x.atts.end());
441 //attsIndex.insert(attsIndex.end(), x.attsIndex.begin(), x.attsIndex.end());
442 hasCategories.insert(hasCategories.end(), x.hasCategories.begin(), x.hasCategories.end());
443 cats.insert(cats.end(), x.cats.begin(), x.cats.end());
444 hasColors.insert(hasColors.end(), x.hasColors.begin(), x.hasColors.end());
445 cols.insert(cols.end(), x.cols.begin(), x.cols.end());
446 datatype.insert(datatype.end(), x.datatype.begin(), x.datatype.end());
447 has_scale_offset.insert(has_scale_offset.end(), x.has_scale_offset.begin(), x.has_scale_offset.end());
448 scale.insert(scale.end(), x.scale.begin(), x.scale.end());
449 offset.insert(offset.end(), x.offset.begin(), x.offset.end());
450 return true;
451 }
452
453
454
combine(SpatRasterSource & x)455 bool SpatRasterSource::combine(SpatRasterSource &x) {
456 if (memory & x.memory) {
457 if ((values.size() + x.values.size()) < (values.max_size()/8) ) {
458 values.insert(values.end(), x.values.begin(), x.values.end());
459 layers.resize(nlyr + x.nlyr);
460 std::iota(layers.begin(), layers.end(), 0);
461 x.values.resize(0);
462 } else {
463 return false;
464 }
465 } else if (filename == x.filename) {
466 layers.insert(layers.end(), x.layers.begin(), x.layers.end());
467 } else {
468 return false;
469 }
470 nlyr += x.nlyr;
471 names.insert(names.end(), x.names.begin(), x.names.end());
472 time.insert(time.end(), x.time.begin(), x.time.end());
473 if (!(hasTime & x.hasTime)) {
474 hasTime = false;
475 }
476 unit.insert(unit.end(), x.unit.begin(), x.unit.end());
477 depth.insert(depth.end(), x.depth.begin(), x.depth.end());
478 hasRange.insert(hasRange.end(), x.hasRange.begin(), x.hasRange.end());
479 range_min.insert(range_min.end(), x.range_min.begin(), x.range_min.end());
480 range_max.insert(range_max.end(), x.range_max.begin(), x.range_max.end());
481 blockrows.insert(blockrows.end(), x.blockrows.begin(), x.blockrows.end());
482 blockcols.insert(blockcols.end(), x.blockcols.begin(), x.blockcols.end());
483 //hasAttributes.insert(hasAttributes.end(), x.hasAttributes.begin(), x.hasAttributes.end());
484 //atts.insert(atts.end(), x.atts.begin(), x.atts.end());
485 //attsIndex.insert(attsIndex.end(), x.attsIndex.begin(), x.attsIndex.end());
486 hasCategories.insert(hasCategories.end(), x.hasCategories.begin(), x.hasCategories.end());
487 cats.insert(cats.end(), x.cats.begin(), x.cats.end());
488 hasColors.insert(hasColors.end(), x.hasColors.begin(), x.hasColors.end());
489 cols.insert(cols.end(), x.cols.begin(), x.cols.end());
490 datatype.insert(datatype.end(), x.datatype.begin(), x.datatype.end());
491 has_scale_offset.insert(has_scale_offset.end(), x.has_scale_offset.begin(), x.has_scale_offset.end());
492 scale.insert(scale.end(), x.scale.begin(), x.scale.end());
493 offset.insert(offset.end(), x.offset.begin(), x.offset.end());
494 return true;
495 }
496
497
498
collapse()499 void SpatRaster::collapse() {
500 size_t n = nsrc();
501 if (n < 2) return;
502 std::vector<size_t> rem;
503 for (size_t i=1; i<n; i++) {
504 if (source[0].combine(source[i])) {
505 rem.push_back(i);
506 }
507 }
508 for (int i=rem.size(); i>= 0; i--) {
509 source.erase(source.begin()+i);
510 }
511 }
512
513
514
collapse_sources()515 SpatRaster SpatRaster::collapse_sources() {
516 SpatRaster out;
517 std::vector<SpatRasterSource> src;
518 SpatRasterSource s = source[0];
519 for (size_t i=1; i<nsrc(); i++) {
520 if (! s.combine_sources(source[i])) {
521 src.push_back(s);
522 s = source[i];
523 }
524 }
525 src.push_back(s);
526 out.setSources(src);
527 return out;
528 }
529
530
531