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
readStart()20 bool SpatRaster::readStart() {
21
22 //if (!valid_sources(true, true)) {
23 // return false;
24 //}
25
26 for (size_t i=0; i<nsrc(); i++) {
27 if (source[i].open_read) {
28 addWarning("source already open for reading");
29 continue;
30 }
31 if (source[i].memory) {
32 source[i].open_read = true;
33 } else if (source[i].multidim) {
34 if (!readStartMulti(i)) {
35 return false;
36 }
37 } else {
38 if (!readStartGDAL(i)) {
39 return false;
40 }
41 }
42 }
43 return true;
44 }
45
readStop()46 bool SpatRaster::readStop() {
47 for (size_t i=0; i<nsrc(); i++) {
48 if (source[i].open_read) {
49 if (source[i].memory) {
50 source[i].open_read = false;
51 } else if (source[i].multidim) {
52 readStopMulti(i);
53 } else {
54 readStopGDAL(i);
55 }
56 }
57 }
58 return true;
59 }
60
61
62 // BSQ
readBlock(BlockSize bs,unsigned i)63 std::vector<double> SpatRaster::readBlock(BlockSize bs, unsigned i){
64 return readValues(bs.row[i], bs.nrows[i], 0, ncol());
65 }
66
67
68 // 2D BSQ
readBlock2(BlockSize bs,unsigned i)69 std::vector<std::vector<double>> SpatRaster::readBlock2(BlockSize bs, unsigned i) {
70 std::vector<double> x = readValues(bs.row[i], bs.nrows[i], 0, ncol());
71 std::vector<std::vector<double>> v(nlyr());
72 size_t off = bs.nrows[i] * ncol();
73 for (size_t i=0; i<nlyr(); i++) {
74 v[i] = std::vector<double>(x.begin()+(i*off), x.begin()+((i+1)*off));
75 }
76 return(v);
77 }
78
79 // BIP
readBlockIP(BlockSize bs,unsigned i)80 std::vector<double> SpatRaster::readBlockIP(BlockSize bs, unsigned i) {
81 std::vector<double> x = readValues(bs.row[i], bs.nrows[i], 0, ncol());
82 std::vector<double> v(x.size());
83 size_t off = bs.nrows[i] * ncol();
84 size_t nl = nlyr();
85 for (size_t i=0; i<nl; i++) {
86 std::vector<double> lyr = std::vector<double>(x.begin()+(i*off), x.begin()+((i+1)*off));
87 for (size_t j=0; j<off; j++){
88 size_t jj = j * nl + i;
89 v[jj] = lyr[j];
90 }
91 }
92 return(v);
93 }
94
95
readChunkMEM(std::vector<double> & out,size_t src,size_t row,size_t nrows,size_t col,size_t ncols)96 void SpatRaster::readChunkMEM(std::vector<double> &out, size_t src, size_t row, size_t nrows, size_t col, size_t ncols){
97
98 size_t nl = source[src].nlyr;
99
100 if (source[src].hasWindow) {
101 row += source[src].window.off_row;
102 col += source[src].window.off_col;
103 size_t endrow = row + nrows;
104 size_t endcol = col + ncols;
105 size_t nc = source[src].window.full_ncol;
106 double ncells = source[src].window.full_nrow * nc;
107
108 for (size_t lyr=0; lyr < nl; lyr++) {
109 size_t add = ncells * lyr;
110 for (size_t r = row; r < endrow; r++) {
111 size_t off = add + r * nc;
112 out.insert(out.end(), source[src].values.begin()+off+col, source[src].values.begin()+off+endcol);
113 }
114 }
115 /*
116 else if (source[0].window.expanded) {
117 unsigned add = ncells * lyr;
118 std::vector<double> v1(source[0].window.expand[0] * ncols, NAN);
119 out.insert(out.end(), v1.begin(), v1.end());
120 v1.resize(source[0].window.expand[1], NAN);
121 std::vector<double> v2(source[0].window.expand[2], NAN);
122 for (size_t r = wrow; r < endrow; r++) {
123 unsigned a = add + r * source[0].window.full_ncol;
124 out.insert(out.end(), v1.begin(), v1.end());
125 out.insert(out.end(), source[src].values.begin()+a+wcol, source[src].values.begin()+a+endcol);
126 out.insert(out.end(), v2.begin(), v2.end());
127 }
128 v1.resize(source[0].window.expand[3] * ncols, NAN);
129 out.insert(out.end(), v1.begin(), v1.end());
130 }
131 */
132
133 } else { // no window
134 if (row==0 && nrows==nrow() && col==0 && ncols==ncol()) {
135 out.insert(out.end(), source[src].values.begin(), source[src].values.end());
136 } else {
137 double ncells = ncell();
138 if (col==0 && ncols==ncol()) {
139 for (size_t lyr=0; lyr < nl; lyr++) {
140 size_t add = ncells * lyr;
141 size_t a = add + row * ncol();
142 size_t b = a + nrows * ncol();
143 out.insert(out.end(), source[src].values.begin()+a, source[src].values.begin()+b);
144 }
145 } else {
146 size_t endrow = row + nrows;
147 size_t endcol = col + ncols;
148 for (size_t lyr=0; lyr < nl; lyr++) {
149 size_t add = ncells * lyr;
150 for (size_t r = row; r < endrow; r++) {
151 size_t a = add + r * ncol();
152 out.insert(out.end(), source[src].values.begin()+a+col, source[src].values.begin()+a+endcol);
153 }
154 }
155 }
156 }
157 }
158 }
159
160
161
readValues(size_t row,size_t nrows,size_t col,size_t ncols)162 std::vector<double> SpatRaster::readValues(size_t row, size_t nrows, size_t col, size_t ncols){
163
164 std::vector<double> out;
165
166 if (((row + nrows) > nrow()) || ((col + ncols) > ncol())) {
167 setError("invalid rows/columns");
168 return out;
169 }
170
171
172 //row = std::min(std::max(size_t(0), row), nrow()-1);
173 //col = std::min(std::max(size_t(0), col), ncol()-1);
174 //nrows = std::max(size_t(1), std::min(nrows, nrow()-row));
175 //ncols = std::max(size_t(1), std::min(ncols, ncol()-col));
176 if ((nrows==0) | (ncols==0)) {
177 return out;
178 }
179
180 if (!hasValues()) {
181 out.resize(nrows * ncols * nlyr(), NAN);
182 addWarning("raster has no values");
183 return out; // or NAs?
184 }
185
186
187 unsigned n = nsrc();
188
189 for (size_t src=0; src<n; src++) {
190 if (source[src].memory) {
191 readChunkMEM(out, src, row, nrows, col, ncols);
192 } else {
193 // read from file
194 #ifdef useGDAL
195
196 /*
197 if (source[0].window.expanded) {
198 std::vector<double> gout;
199 readChunkGDAL(gout, src, source[0].window.off_row, nrows, source[0].window.off_col, ncols);
200
201 size_t rrow = row + source[0].window.off_row;
202 size_t rcol = col + source[0].window.off_col;
203 unsigned endrow = rrow + nrows;
204 unsigned endcol = rcol + ncols;
205 unsigned ncells = source[0].window.full_nrow * source[0].window.full_ncol;
206 unsigned nl = source[src].nlyr;
207
208 for (size_t lyr=0; lyr < nl; lyr++) {
209 unsigned add = ncells * lyr;
210 std::vector<double> v1(source[0].window.expand[0] * ncols, NAN);
211 out.insert(out.end(), v1.begin(), v1.end());
212 v1.resize(source[0].window.expand[1], NAN);
213 std::vector<double> v2(source[0].window.expand[2], NAN);
214 for (size_t r = rrow; r < endrow; r++) {
215 unsigned a = add + r * source[0].window.full_ncol;
216 out.insert(out.end(), v1.begin(), v1.end());
217 out.insert(out.end(), gout.begin()+a+rcol, gout.begin()+a+endcol);
218 out.insert(out.end(), v2.begin(), v2.end());
219 }
220 v1.resize(source[0].window.expand[3] * ncols, NAN);
221 out.insert(out.end(), v1.begin(), v1.end());
222 }
223 }
224 */
225 readChunkGDAL(out, src, row, nrows, col, ncols);
226 #endif // useGDAL
227 }
228 }
229 return out;
230 }
231
232
233
readAll()234 bool SpatRaster::readAll() {
235 if (!hasValues()) {
236 return true;
237 }
238
239 size_t row =0, col=0, nrows=nrow(), ncols=ncol();
240 readStart();
241 size_t n = nsrc();
242 for (size_t src=0; src<n; src++) {
243 if (!source[src].memory) {
244 readChunkGDAL(source[src].values, src, row, nrows, col, ncols);
245 source[src].memory = true;
246 source[src].filename = "";
247 }
248 if (src > 0) {
249 if (!source[0].combine_sources(source[src])) {
250 setError("could not combine sources");
251 return false;
252 }
253 source[src].values.resize(0);
254 }
255 }
256 readStop();
257 if (n>1) source.resize(1);
258 return true;
259 }
260
261
getValues(long lyr,SpatOptions & opt)262 std::vector<double> SpatRaster::getValues(long lyr, SpatOptions &opt) {
263
264 std::vector<double> out;
265
266 bool hw = false;
267 for (size_t i=0; i<source.size(); i++) {
268 if (source[i].hasWindow) {
269 hw = true;
270 break;
271 }
272 }
273 if (hw) {
274 if (!readStart()) return out;
275 if (lyr < 0) { // default; read all
276 out = readValues(0, nrow(), 0, ncol());
277 } else {
278 unsigned lyrnr = lyr;
279 SpatRaster sub = subset({lyrnr}, opt);
280 out = sub.readValues(0, nrow(), 0, ncol());
281 }
282 readStop();
283 return out;
284 }
285
286 if (lyr < 0) { // default; read all
287 unsigned n = nsrc();
288 for (size_t src=0; src<n; src++) {
289 if (source[src].memory) {
290 out.insert(out.end(), source[src].values.begin(), source[src].values.end());
291 } else {
292 #ifdef useGDAL
293 std::vector<double> fvals = readValuesGDAL(src, 0, nrow(), 0, ncol());
294 out.insert(out.end(), fvals.begin(), fvals.end());
295 #endif // useGDAL
296 }
297 }
298 } else { // read one lyr
299 std::vector<unsigned> sl = findLyr(lyr);
300 unsigned src=sl[0];
301 if (source[src].memory) {
302 size_t start = sl[1] * ncell();
303 out = std::vector<double>(source[src].values.begin()+start, source[src].values.begin()+start+ncell());
304 } else {
305 #ifdef useGDAL
306 out = readValuesGDAL(src, 0, nrow(), 0, ncol(), sl[1]);
307 #endif // useGDAL
308 }
309 }
310 return out;
311 }
312
313
getValuesSource(size_t src,std::vector<double> & out)314 bool SpatRaster::getValuesSource(size_t src, std::vector<double> &out) {
315
316 unsigned n = nsrc();
317 if (src > n) {
318 return false;
319 }
320
321 bool hw = false;
322 for (size_t i=0; i<source.size(); i++) {
323 if (source[i].hasWindow) {
324 hw = true;
325 break;
326 }
327 }
328 if (hw) {
329 SpatRaster sub = SpatRaster(source[src]);
330 if (!readStart()) return false;
331 out = sub.readValues(0, nrow(), 0, ncol());
332 readStop();
333 return true;
334 }
335
336
337 if (source[src].memory) {
338 out = std::vector<double>(source[src].values.begin(), source[src].values.end());
339 } else {
340 #ifdef useGDAL
341 out = readValuesGDAL(src, 0, nrow(), 0, ncol());
342 #endif // useGDAL
343 }
344 return true;
345 }
346
347