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 "spatRaster.h"
19 #include "file_utils.h"
20 #include "string_utils.h"
21 #include "math_utils.h"
22
23
24
writeValuesMem(std::vector<double> & vals,size_t startrow,size_t nrows,size_t startcol,size_t ncols)25 bool SpatRaster::writeValuesMem(std::vector<double> &vals, size_t startrow, size_t nrows, size_t startcol, size_t ncols) {
26
27 if (vals.size() == size()) {
28 source[0].values = vals;
29 return true;
30 }
31
32 // this seems rather inefficient for most cases
33 // where values could simply be appended
34 if (source[0].values.size() == 0) { // && startrow != 0 && startcol != 0) {
35 source[0].values = std::vector<double>(size(), NAN);
36 }
37
38 size_t nc = ncell();
39 size_t chunk = nrows * ncols;
40
41 //complete rows
42 if (startcol==0 && ncols==ncol()) {
43 for (size_t i=0; i<nlyr(); i++) {
44 size_t off1 = i * chunk;
45 size_t off2 = startrow * ncols + i * nc;
46 std::copy( vals.begin()+off1, vals.begin()+off1+chunk, source[0].values.begin()+off2 );
47 }
48
49 // block writing
50 } else {
51 for (size_t i=0; i<nlyr(); i++) {
52 unsigned off = i*chunk;
53 for (size_t r=0; r<nrows; r++) {
54 size_t start1 = r * ncols + off;
55 size_t start2 = (startrow+r)*ncol() + i*nc + startcol;
56 std::copy(vals.begin()+start1, vals.begin()+start1+ncols, source[0].values.begin()+start2);
57 }
58 }
59 }
60 return true;
61 }
62
63
64
fill(double x)65 void SpatRaster::fill(double x) {
66 if (source[0].driver == "gdal") {
67 #ifdef useGDAL
68 fillValuesGDAL(x);
69 #endif
70 } else {
71 source[0].values.resize(size(), x);
72 }
73
74 }
75
76
77
isSource(std::string filename)78 bool SpatRaster::isSource(std::string filename) {
79 std::vector<std::string> ff = filenames();
80 for (size_t i=0; i<ff.size(); i++) {
81 if (ff[i] == filename) {
82 return true;
83 }
84 }
85 return false;
86 }
87
88 /*
89 #include <experimental/filesystem>
90 bool SpatRaster::differentFilenames(std::vector<std::string> outf) {
91 std::vector<std::string> inf = filenames();
92 for (size_t i=0; i<inf.size(); i++) {
93 if (inf[i] == "") continue;
94 std::experimental::filesystem::path pin = inf[i];
95 for (size_t j=0; j<outf.size(); j++) {
96 std::experimental::filesystem::path pout = outf[i];
97 if (pin.compare(pout) == 0) return false;
98 }
99 }
100 return true;
101 }
102 */
103
differentFilenames(std::vector<std::string> outf,bool & duplicates,bool & empty)104 bool SpatRaster::differentFilenames(std::vector<std::string> outf, bool &duplicates, bool &empty) {
105 std::vector<std::string> inf = filenames();
106 duplicates = false;
107 empty = false;
108 for (size_t j=0; j<outf.size(); j++) {
109 if (outf[j] == "") {
110 empty = true;
111 return false;
112 }
113 }
114
115 for (size_t i=0; i<inf.size(); i++) {
116 if (inf[i] == "") continue;
117 #ifdef _WIN32
118 lowercase(inf[i]);
119 #endif
120 for (size_t j=0; j<outf.size(); j++) {
121 #ifdef _WIN32
122 lowercase(outf[j]);
123 #endif
124 if (inf[i] == outf[j]) return false;
125 }
126 }
127
128 size_t n = outf.size();
129 outf.erase(std::unique(outf.begin(), outf.end()), outf.end());
130 if (n > outf.size()) {
131 duplicates = true;
132 return false;
133 }
134 return true;
135 }
136
137
138
writeRaster(SpatOptions & opt)139 SpatRaster SpatRaster::writeRaster(SpatOptions &opt) {
140
141 // here we could check if we can simple make a copy if
142 // a) the SpatRaster is backed by a file
143 // b) there are no write options
144
145 SpatRaster out = geometry(nlyr(), true);
146 if (!hasValues()) {
147 out.setError("there are no cell values");
148 return out;
149 }
150
151 // recursive writing of layers
152 std::vector<std::string> fnames = opt.get_filenames();
153 bool dups, empty;
154 if (!differentFilenames(fnames, dups, empty)) {
155 if (dups) {
156 out.setError("duplicate filenames");
157 } else if (empty) {
158 out.setError("empty filename");
159 } else {
160 out.setError("source and target filename cannot be the same");
161 }
162 return(out);
163 }
164
165 size_t nl = nlyr();
166 if (fnames.size() > 1) {
167 if (fnames.size() != nl) {
168 out.setError("the number of filenames should either be one, or equal to the number of layers");
169 return out;
170 } else {
171 bool overwrite = opt.get_overwrite();
172 std::string errmsg;
173 for (size_t i=0; i<nl; i++) {
174 if (!can_write(fnames[i], overwrite, errmsg)) {
175 out.setError(errmsg + " (" + fnames[i] +")");
176 return(out);
177 }
178 }
179 for (unsigned i=0; i<nl; i++) {
180 opt.set_filenames({fnames[i]});
181 SpatRaster out = subset({i}, opt);
182 if (out.hasError()) {
183 return out;
184 }
185 fnames[i] = out.source[0].filename;
186 }
187 SpatRaster out(fnames, {-1}, {""}, false, {}, {});
188 return out;
189 }
190 }
191
192 if (!readStart()) {
193 out.setError(getError());
194 return(out);
195 }
196
197 opt.ncopies = 2;
198 if (!out.writeStart(opt)) {
199 readStop();
200 return out;
201 }
202 for (size_t i=0; i<out.bs.n; i++) {
203 std::vector<double> v = readBlock(out.bs, i);
204 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, ncol())) {
205 readStop();
206 out.writeStop();
207 return out;
208 }
209 }
210 out.writeStop();
211 //if (!out.writeStopGDAL()) {
212 // out.setError("cannot close file");
213 //}
214 readStop();
215 return out;
216 }
217
218
219
220
writeStart(SpatOptions & opt)221 bool SpatRaster::writeStart(SpatOptions &opt) {
222
223 if (opt.names.size() == nlyr()) {
224 setNames(opt.names);
225 }
226
227 std::vector<std::string> fnames = opt.get_filenames();
228 if (fnames.size() > 1) {
229 addWarning("only the first filename supplied is used");
230 }
231 std::string filename = fnames[0];
232 if (filename == "") {
233 if (!canProcessInMemory(opt)) {
234 std::string extension = ".tif";
235 filename = tempFile(opt.get_tempdir(), extension);
236 opt.set_filenames({filename});
237 opt.gdal_options = {"COMPRESS=NONE"};
238 }
239 }
240
241 if (filename != "") {
242 // open GDAL filestream
243 #ifdef useGDAL
244 if (! writeStartGDAL(opt) ) {
245 return false;
246 }
247 #else
248 setError("GDAL is not available");
249 return false;
250 #endif
251 }
252 if (source[0].open_write) {
253 addWarning("file was already open");
254 }
255 source[0].open_write = true;
256 source[0].filename = filename;
257 bs = getBlockSize(opt);
258 #ifdef useRcpp
259 if (opt.verbose) {
260 std::vector<double> mems = mem_needs(opt);
261 double gb = 1073741824 / 8;
262 //{memneed, memavail, frac, csize, inmem} ;
263 // << "max vect size : " << roundn(mems.max_size() / gb, 2) << " GB" << std::endl;
264 Rcpp::Rcout<< "memory avail. : " << roundn(mems[1] / gb, 2) << " GB" << std::endl;
265 Rcpp::Rcout<< "memory allow. : " << roundn(mems[2] * mems[1] / gb, 2) << " GB" << std::endl;
266 Rcpp::Rcout<< "memory needed : " << roundn(mems[0] / gb, 3) << " GB" << " (" << opt.ncopies << " copies)" << std::endl;
267 std::string inmem = mems[4] < 0.5 ? "false" : "true";
268 Rcpp::Rcout<< "in memory : " << inmem << std::endl;
269 Rcpp::Rcout<< "block size : " << mems[3] << " rows" << std::endl;
270 Rcpp::Rcout<< "n blocks : " << bs.n << std::endl;
271 Rcpp::Rcout<< "pb : " << opt.show_progress(bs.n) << std::endl;
272 Rcpp::Rcout<< std::endl;
273 }
274
275 if (opt.progressbar) {
276 unsigned long steps = bs.n+2;
277 pbar = new Progress(steps, opt.show_progress(bs.n));
278 pbar->increment();
279 progressbar = true;
280 } else {
281 progressbar = false;
282 }
283 #endif
284 return true;
285 }
286
287
288
writeValues(std::vector<double> & vals,size_t startrow,size_t nrows,size_t startcol,size_t ncols)289 bool SpatRaster::writeValues(std::vector<double> &vals, size_t startrow, size_t nrows, size_t startcol, size_t ncols) {
290 bool success = true;
291
292 if (!source[0].open_write) {
293 setError("cannot write (no open file)");
294 return false;
295 }
296
297 if ((startrow + nrows) > nrow()) {
298 setError("incorrect start and/or nrows value");
299 return false;
300 }
301
302 if (source[0].driver == "gdal") {
303 #ifdef useGDAL
304
305 success = writeValuesGDAL(vals, startrow, nrows, startcol, ncols);
306 #else
307 setError("GDAL is not available");
308 return false;
309 #endif
310 } else {
311 success = writeValuesMem(vals, startrow, nrows, startcol, ncols);
312 }
313
314 #ifdef useRcpp
315 if (progressbar) {
316 if (Progress::check_abort()) {
317 pbar->cleanup();
318 setError("aborted");
319 return(false);
320 }
321 pbar->increment();
322 }
323 #endif
324 return success;
325 }
326
327
328
329
330 template <typename T>
flatten(const std::vector<std::vector<T>> & v)331 std::vector<T> flatten(const std::vector<std::vector<T>>& v) {
332 std::size_t total_size = 0;
333 for (const auto& sub : v)
334 total_size += sub.size();
335 std::vector<T> result;
336 result.reserve(total_size);
337 for (const auto& sub : v)
338 result.insert(result.end(), sub.begin(), sub.end());
339 return result;
340 }
341
342
writeValues2(std::vector<std::vector<double>> & vals,size_t startrow,size_t nrows,size_t startcol,size_t ncols)343 bool SpatRaster::writeValues2(std::vector<std::vector<double>> &vals, size_t startrow, size_t nrows, size_t startcol, size_t ncols) {
344 std::vector<double> vv = flatten(vals);
345 return writeValues(vv, startrow, nrows, startcol, ncols);
346 }
347
writeStop()348 bool SpatRaster::writeStop(){
349 if (!source[0].open_write) {
350 setError("cannot close a file that is not open");
351 return false;
352 }
353 source[0].open_write = false;
354 bool success = true;
355 source[0].memory = false;
356 if (source[0].driver=="gdal") {
357 #ifdef useGDAL
358 success = writeStopGDAL();
359 //source[0].hasValues = true;
360 #else
361 setError("GDAL is not available");
362 return false;
363 #endif
364 } else {
365 source[0].setRange();
366 //source[0].driver = "memory";
367 source[0].memory = true;
368 if (source[0].values.size() > 0) {
369 source[0].hasValues = true;
370 }
371 }
372
373 #ifdef useRcpp
374 if (progressbar) {
375 if (Progress::check_abort()) {
376 pbar->cleanup();
377 setError("aborted");
378 return(false);
379 }
380 pbar->increment();
381 delete pbar;
382 }
383 #endif
384
385 return success;
386 }
387
388 //bool SpatRaster::replaceValues(std::vector<double> cells, std::vector<double> _values, int ncols) {
389 //}
390
setValues(std::vector<double> & v,SpatOptions & opt)391 bool SpatRaster::setValues(std::vector<double> &v, SpatOptions &opt) {
392 SpatRaster g = geometry();
393 SpatRasterSource s = g.source[0];
394 s.hasValues = true;
395 s.memory = true;
396 s.names = getNames();
397 s.driver = "memory";
398 setSource(s);
399
400 if (v.size() < g.size()) {
401 *this = init(v, opt);
402 return (!hasError());
403 } else if (v.size() == g.size()) {
404 /*
405 if (!canProcessInMemory(opt)) {
406 // this should be chunked to avoid the copy
407 // but this may still help in some cases
408 source[0].values = v;
409 std::string f = opt.get_filename();
410 if (f == "") {
411 std::string filename = tempFile(opt.get_tempdir(), ".tif");
412 opt.set_filenames({filename});
413 }
414 *this = writeRaster(opt);
415 return true;
416 } else {
417 */
418 source[0].values = v;
419 // }
420 } else {
421
422 setError("incorrect number of values");
423 return false;
424 }
425 source[0].setRange();
426 return true;
427 }
428
setRange(SpatOptions & opt)429 void SpatRaster::setRange(SpatOptions &opt) {
430
431 for (size_t i=0; i<nsrc(); i++) {
432 if (source[i].hasRange[0]) continue;
433
434 if (source[i].memory) {
435 source[i].setRange();
436 } else {
437 SpatRaster r(source[i]);
438 SpatDataFrame x = r.global("range", true, opt);
439 source[i].range_min = x.getD(0);
440 source[i].range_max = x.getD(1);
441 source[i].hasRange = std::vector<bool>(source[i].hasRange.size(), true);
442 }
443 }
444 }
445
setRange()446 void SpatRasterSource::setRange() {
447 double vmin, vmax;
448 size_t nc = ncol * nrow;
449 size_t start;
450 range_min.resize(nlyr);
451 range_max.resize(nlyr);
452 hasRange.resize(nlyr);
453 if (values.size() == nc * nlyr) {
454 for (size_t i=0; i<nlyr; i++) {
455 start = nc * i;
456 minmax(values.begin()+start, values.begin()+start+nc, vmin, vmax);
457 range_min[i] = vmin;
458 range_max[i] = vmax;
459 hasRange[i] = true;
460 }
461 }
462 }
463
464
465
466
467
468
469 /*
470 bool SpatRaster::writeStartFs(std::string filename, std::string format, std::string datatype, bool overwrite, fstream& f) {
471
472 lrtrim(filename);
473 if (filename == "") {
474 if (!canProcessInMemory(4)) {
475 filename = "random_file_name.grd";
476 }
477 }
478
479 if (filename == "") {
480 source.driver = {"memory"};
481
482 } else {
483 // if (!overwrite) check if file exists
484 string ext = getFileExt(filename);
485 lowercase(ext);
486 if (ext == ".grd") {
487 source.driver = {"raster"};
488 string fname = setFileExt(filename, ".gri");
489 f.open(fname, ios::out | ios::binary);
490 (*fs).open(fname, ios::out | ios::binary);
491 fs = &f;
492 } else {
493 source.driver = {"gdal"} ;
494 // open GDAL filestream
495 }
496 }
497
498 source.filename = {filename};
499 bs = getBlockSize(4);
500 return true;
501 }
502 */
503
504
505