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 <functional>
19
20 #include "spatRasterMultiple.h"
21 #include "distance.h"
22 #include "vecmath.h"
23
rowColToCell(unsigned ncols,unsigned row,unsigned col)24 double rowColToCell(unsigned ncols, unsigned row, unsigned col) {
25 return row * ncols + col;
26 }
27
28
fourCellsFromXY(const std::vector<double> & x,const std::vector<double> & y)29 std::vector<double> SpatRaster::fourCellsFromXY(const std::vector<double> &x, const std::vector<double> &y) {
30
31 size_t n = x.size();
32 SpatExtent e = getExtent();
33 std::vector<double> out;
34 out.reserve(4*n);
35
36 double xmin = e.xmin;
37 double xmax = e.xmax;
38 double xr = xres();
39 double ymin = e.ymin;
40 double ymax = e.ymax;
41 double yr = yres();
42 double nc = ncol();
43 int_64 mxr = nrow()-1;
44 int_64 mxc = ncol()-1;
45 int_64 r1, r2, c1, c2;
46 std::vector<double> bad = {NAN, NAN, NAN, NAN};
47 for (size_t i = 0; i < n; i++) {
48 if (y[i] < ymin || y[i] > ymax || x[i] < xmin || x[i] > xmax) {
49 out.insert(out.end(), bad.begin(), bad.end());
50 continue;
51 }
52 if (y[i] == ymin) {
53 r1 = mxr;
54 r2 = mxr;
55 } else {
56 double p = (ymax - y[i]) / yr;
57 r1 = trunc(p);
58 if ((p - r1) > 0.5) {
59 r2 = r1 == mxr ? mxr : r1 + 1;
60 } else {
61 r2 = r1;
62 r1 = r1 == 0 ? 0 : r1 - 1;
63 }
64 }
65 if (x[i] == xmax) {
66 c1 = mxc;
67 c2 = mxc;
68 } else {
69 double p = (x[i] - xmin) / xr;
70 c1 = trunc(p);
71 if ((p - c1) > 0.5) {
72 c2 = c1 == mxc ? mxc : c1 + 1;
73 } else {
74 c2 = c1;
75 c1 = c2 == 0 ? 0 : c2 - 1;
76 }
77 }
78 out.push_back(r1 * nc + c1);
79 out.push_back(r1 * nc + c2);
80 out.push_back(r2 * nc + c1);
81 out.push_back(r2 * nc + c2);
82 }
83 return out;
84 }
85
86 /*
87 double linearInt(const double& d, const double& x, const double& x1, const double& x2, const double& v1, const double& v2) {
88 double result = (v2 * (x - x1) + v1 * (x2 - x)) / d;
89 return result;
90 }
91 */
92
93 /* ok but cannot handle NA
94 double bilinearInt(const double& x, const double& y, const double& x1, const double& x2, const double& y1, const double& y2, const double& v11, const double& v21, const double& v12, const double& v22) {
95 double d = x2-x1;
96 double h1 = linearInt(d, x, x1, x2, v11, v21);
97 double h2 = linearInt(d, x, x1, x2, v12, v22);
98 d = y2-y1;
99 double v = linearInt(d, y, y1, y2, h1, h2);
100 return v;
101 }
102 */
103
104 /*
105 double bilinearIntold(const double& x, const double& y,
106 const double& x1, const double& x2, const double& y1, const double& y2,
107 const double& v11, const double& v21, const double& v12, const double& v22) {
108 double d = x2-x1;
109 double h1=NAN;
110 double h2=NAN;
111 if (!std::isnan(v11) && !std::isnan(v21)) {
112 h1 = linearInt(d, x, x1, x2, v11, v21);
113 } else if (!std::isnan(v11)) {
114 h1 = v11;
115 } else if (!std::isnan(v21)) {
116 h1 = v21;
117 }
118
119 if (!std::isnan(v12) && !std::isnan(v22)) {
120 h2 = linearInt(d, x, x1, x2, v12, v22);
121 } else if (!std::isnan(v12)) {
122 h2 = v12;
123 } else if (!std::isnan(v22)) {
124 h2 = v22;
125 }
126 if (!std::isnan(h1) && !std::isnan(h2)) {
127 d = y2-y1;
128 double v = linearInt(d, y, y1, y2, h1, h2);
129 return v;
130 } else if (!std::isnan(h1)) {
131 return h1;
132 } else if (!std::isnan(h2)) {
133 return h2;
134 }
135 return NAN;
136 }
137 */
138
139
bilinearInt(const double & x,const double & y,const double & x1,const double & x2,const double & y1,const double & y2,double & v11,double & v12,double & v21,double & v22,bool weights)140 std::vector<double> bilinearInt(const double& x, const double& y,
141 const double& x1, const double& x2, const double& y1, const double& y2,
142 double& v11, double& v12, double& v21, double& v22, bool weights) {
143
144 bool n1 = std::isnan(v11);
145 bool n2 = std::isnan(v12);
146 bool n3 = std::isnan(v21);
147 bool n4 = std::isnan(v22);
148 if (std::isnan(x) || std::isnan(y) || (n1 && n2 && n3 && n4)) {
149 if (weights) {
150 std::vector<double> out(4, NAN);
151 return out;
152 }
153 std::vector<double> out(1, NAN);
154 return out;
155 }
156 double dx = (x2 - x1);
157 bool intx = dx > 0;
158 double dy = (y1 - y2);
159 bool inty = dy > 0;
160 double w11, w12, w21, w22;
161 if (weights) {
162 v11 = 1;
163 v12 = 1;
164 v21 = 1;
165 v22 = 1;
166 }
167
168 if (intx & inty) {
169 double d = dx * dy;
170 if (!(n1 || n2)) {
171 w11 = v11 * ((x2 - x) * (y - y2)) / d;
172 w12 = v12 * ((x - x1) * (y - y2)) / d;
173 } else {
174 w11 = n1 ? 0.0 : 0.5 * v11;
175 w12 = n2 ? 0.0 : 0.5 * v12;
176 }
177 if (!(n3 || n4)) {
178 w21 = v21 * ((x2 - x) * (y1 - y)) / d;
179 w22 = v22 * ((x - x1) * (y1 - y)) / d;
180 } else {
181 w21 = n3 ? 0.0 : 0.5 * v21;
182 w22 = n4 ? 0.0 : 0.5 * v22;
183 }
184 } else if (intx) {
185 if (!(n1 || n2)) {
186 w11 = v11 * (x2 - x) / dx;
187 w12 = v12 * (x - x1) / dx;
188 } else {
189 w11 = n1 ? 0.0 : v11;
190 w12 = n2 ? 0.0 : v12;
191 }
192 w21 = 0;
193 w22 = 0;
194 } else if (inty) {
195 if (!(n1 || n3)) {
196 w11 = v11 * (y - y2) / dy;
197 w21 = v21 * (y1 - y) / dy;
198 } else {
199 w11 = n1 ? 0.0 : v11;
200 w21 = n3 ? 0.0 : v21;
201 }
202 w12 = 0;
203 w22 = 0;
204 } else {
205 w11 = v11;
206 w21 = 0;
207 w12 = 0;
208 w22 = 0;
209 }
210
211 if (weights) {
212 std::vector<double> out = { w11, w12, w21, w22 };
213 return out;
214 }
215 std::vector<double> out = { w11 + w12 + w21 + w22 };
216 return out;
217 }
218
219
220 /*
221 double distInt(double d, double pd1, double pd2, double v1, double v2) {
222 double result = (v2 * pd1 + v1 * pd2) / d;
223 return result;
224 }
225
226
227 double bilinear_geo(double x, double y, double x1, double x2, double y1, double y2, double halfyres, std::vector<double> vv) {
228
229 double a = 6378137.0;
230 double f = 1/298.257223563;
231 double hy = y1 - halfyres;
232 double d = distance_lonlat(x1, hy, x2, hy, a, f);
233
234 std::vector<double> dist(4);
235 double pd1 = distance_lonlat(x, hy, x1, hy, a, f);
236 double pd2 = distance_lonlat(x, hy, x2, hy, a, f);
237 double h1 = distInt(d, pd1, pd2, vv[0], vv[1]);
238 double h2 = distInt(d, pd1, pd2, vv[2], vv[3]);
239 d = y2 - y1;
240 double v = linearInt(d, y, y1, y2, h1, h2);
241
242 return v;
243 }
244 */
245
246
bilinearValues(const std::vector<double> & x,const std::vector<double> & y)247 std::vector<std::vector<double>> SpatRaster::bilinearValues(const std::vector<double> &x, const std::vector<double> &y) {
248
249
250 std::vector<double> four = fourCellsFromXY(x, y);
251 std::vector<std::vector<double>> xy = xyFromCell(four);
252 std::vector<std::vector<double>> v = extractCell(four);
253 size_t n = x.size();
254 // double halfyres = yres()/2;
255 std::vector<std::vector<double>> res(nlyr(), std::vector<double>(n));
256 /* if (lonlat) {
257 for (size_t i=0; i<n; i++) {
258 size_t ii = i * 4;
259 for (size_t j=0; j<nlyr(); j++) {
260 std::vector<double> vv(v[j].begin()+ii, v[j].begin()+ii+4);
261 if (i==0) {
262 std::cout << x[i] << " "<< y[i] << "\n";
263 std::cout << xy[0][ii] << " " << xy[0][ii+1] << " " << xy[1][ii] << " " << xy[1][ii+3] << "\n";
264 std::cout << v[j][ii] << " " << v[j][ii+1] << " " << v[j][ii+2] << " " << v[j][ii+3] << "\n";
265 }
266 res[j][i] = bilinear_geo(x[i], y[i], xy[0][ii], xy[0][ii+1], xy[1][ii], xy[1][ii+3], halfyres, vv);
267
268 }
269 }
270 } else {
271
272 */
273 for (size_t i=0; i<n; i++) {
274 size_t ii = i * 4;
275 for (size_t j=0; j<nlyr(); j++) {
276 /* if (i==0) {
277 std::cout << x[i] << " "<< y[i] << "\n";
278 std::cout << xy[0][ii] << " " << xy[0][ii+1] << " " << xy[1][ii] << " " << xy[1][ii+3] << "\n";
279 std::cout << v[j][ii] << " " << v[j][ii+1] << " " << v[j][ii+2] << " " << v[j][ii+3] << "\n";
280 }
281 */
282 std::vector<double> value = bilinearInt(x[i], y[i], xy[0][ii], xy[0][ii+1], xy[1][ii], xy[1][ii+3], v[j][ii], v[j][ii+1], v[j][ii+2], v[j][ii+3], false);
283 res[j][i] = value[0];
284 }
285 }
286 // }
287 return res;
288 }
289
290
bilinearCells(const std::vector<double> & x,const std::vector<double> & y)291 std::vector<double> SpatRaster::bilinearCells(const std::vector<double> &x, const std::vector<double> &y) {
292 std::vector<double> four = fourCellsFromXY(x, y);
293 std::vector<std::vector<double>> xy = xyFromCell(four);
294 std::vector<std::vector<double>> v = extractCell(four);
295 size_t n = x.size();
296 std::vector<double> res;
297 std::vector<double> w;
298 for (size_t i=0; i<n; i++) {
299 size_t ii = i * 4;
300 size_t j=0;
301 std::vector<double> w = bilinearInt(x[i], y[i], xy[0][ii], xy[0][ii+1], xy[1][ii], xy[1][ii+3], v[j][ii], v[j][ii+1], v[j][ii+2], v[j][ii+3], true);
302 res.insert(res.end(), four.begin()+ii, four.begin()+ii+4);
303 res.insert(res.end(), w.begin(), w.end());
304 }
305 return res;
306 }
307
308
309
bilinear(const std::vector<double> & v,const std::vector<double> & e,const double & dxdy,const double & x,const double & y)310 double bilinear(const std::vector<double> &v, const std::vector<double> &e, const double &dxdy, const double &x, const double &y) {
311 // values
312 // v[0] v[1]
313 // v[2] v[3]
314
315 // coordinates
316 // e[3] (ymax)
317 // (xmin)e[0] e[1] (xmax)
318 // e[2] (ymin)
319
320 double dx1 = x - e[0];
321 double dx2 = e[1] - x;
322 double dy1 = y - e[2];
323 double dy2 = e[3] - y;
324 return (v[2] * dx2 * dy2 + v[3] * dx1 * dy2 + v[0] * dx2 * dy1 + v[1] * dx1 * dy1) / dxdy;
325 }
326
327
328
line_cells(SpatGeom & g)329 std::vector<double> SpatRaster::line_cells(SpatGeom& g) {
330
331 unsigned nrows = nrow();
332 unsigned ncols = ncol();
333 SpatExtent extent = getExtent();
334
335 double xmin = extent.xmin;
336 double ymax = extent.ymax;
337 double rx = xres();
338 double ry = yres();
339 std::vector<double> out;
340
341 unsigned np = g.size();
342 for (size_t prt=0; prt<np; prt++) {
343 SpatPart p = g.getPart(prt);
344 double miny = vmin(p.y, true);
345 double maxy = vmax(p.y, true);
346
347 double minrow = rowFromY(miny);
348 double maxrow = rowFromY(maxy);
349 if (minrow > nrows || maxrow < 0) {
350 return(out);
351 }
352 size_t startrow = minrow < 0 ? 0 : minrow;
353 size_t endrow = maxrow >= nrows ? (nrows-1) : maxrow;
354 unsigned n = p.x.size();
355 out.reserve(2*(startrow-endrow+1));
356
357 for (size_t row=startrow; row<endrow; row++) {
358 double y = ymax - (row+0.5) * ry;
359 unsigned rowcell = ncols * row;
360 for (size_t i=1; i<n; i++) {
361 size_t j = i-1;
362 if (((p.y[i] < y) && (p.y[j] >= y)) || ((p.y[j] < y) && (p.y[i] >= y))) {
363 double col = ((p.x[i] - xmin + (y-p.y[i])/(p.y[j]-p.y[i]) * (p.x[j]-p.x[i])) + 0.5 * rx ) / rx;
364 if ((col >= 0) & (col < ncols)) {
365 out.push_back(rowcell + col);
366 }
367 }
368 }
369 }
370 }
371 return(out);
372 }
373
374
375
376
polygon_cells(SpatGeom & g)377 std::vector<double> SpatRaster::polygon_cells(SpatGeom& g) {
378
379 // does not deal with holes yet.
380
381 unsigned nrows = nrow();
382 unsigned ncols = ncol();
383 SpatExtent extent = getExtent();
384
385 double xmin = extent.xmin;
386 double ymax = extent.ymax;
387 double rx = xres();
388 double ry = yres();
389 std::vector<double> out;
390
391 unsigned np = g.size();
392 for (size_t prt=0; prt<np; prt++) {
393
394 SpatPart p = g.getPart(prt);
395 double miny = vmin(p.y, true);
396 double maxy = vmax(p.y, true);
397 double minrow = rowFromY(miny);
398 double maxrow = rowFromY(maxy);
399 if (minrow > nrows || maxrow < 0) {
400 return(out);
401 }
402 size_t startrow = minrow < 0 ? 0 : minrow;
403 size_t endrow = maxrow >= nrows ? (nrows-1) : maxrow;
404 unsigned n = p.x.size();
405 out.reserve(5*(startrow-endrow+1));
406
407 std::vector<unsigned> nCol(n);
408 for (size_t row=0; row<nrows; row++) {
409 double y = ymax - (row+0.5) * ry;
410 // find nodes.
411 unsigned nodes = 0;
412 size_t j = n-1;
413 for (size_t i=0; i<n; i++) {
414 if (((p.y[i] < y) && (p.y[j] >= y)) || ((p.y[j] < y) && (p.y[i] >= y))) {
415 // nCol[nodes++]=(int) (((pX[i] - xmin + (y-pY[i])/(pY[j]-pY[i]) * (pX[j]-pX[i])) + 0.5 * rx ) / rx);
416 double nds = ((p.x[i] - xmin + (y-p.y[i])/(p.y[j]-p.y[i]) * (p.x[j]-p.x[i])) + 0.5 * rx ) / rx;
417 nds = nds < 0 ? 0 : nds;
418 nds = nds > ncols ? ncols : nds;
419 nCol[nodes] = (unsigned) nds;
420 nodes++;
421 }
422 j = i;
423 }
424
425 // now remove the holes?
426
427 std::sort(nCol.begin(), nCol.begin()+nodes);
428 unsigned rowcell = ncols * row;
429
430 // fill cells between node pairs.
431 for (size_t i=0; i < nodes; i+=2) {
432 if (nCol[i+1] > 0 && nCol[i] < ncols) { // surely should be >= 0?
433 for (size_t col = nCol[i]; col < nCol[i+1]; col++) {
434 out.push_back(col + rowcell);
435 }
436 }
437 }
438 }
439 }
440 return(out);
441 }
442
443
444 /*
445 idw
446 bool lonlat = could_be_lonlat();
447 //bool globalLonLat = is_global_lonlat();
448 //size_t n = x.size();
449
450 if (method == "idw") {
451 std::function<std::vector<double>(std::vector<double>&,std::vector<double>&,double,double)> distFun;
452 // std::vector<double> distance_plane(std::vector<double> &x1, std::vector<double> &y1, std::vector<double> &x2, std::vector<double> &y2);
453 if (lonlat) {
454 distFun = distance_lonlat_vd;
455 } else {
456 distFun = distance_plane_vd;
457 }
458 */
459 /*
460 cxy = xyFromCell(cells);
461 d = distFun(cxy[0], cxy[1], x[i], y[i]);
462 v = extractCell(cells);
463
464 double a=0, b=0;
465 for (size_t j=0; j<4; j++) {
466 a += v[j] * d[j];
467 b += d[j];
468 }
469 out[i] = a / b;
470 */
471
472
473 // <layer<values>>
extractXY(const std::vector<double> & x,const std::vector<double> & y,const std::string & method,const bool & cells)474 std::vector<std::vector<double>> SpatRaster::extractXY(const std::vector<double> &x, const std::vector<double> &y, const std::string & method, const bool &cells) {
475
476 unsigned nl = nlyr();
477 unsigned np = x.size();
478 if (!hasValues()) {
479 std::vector<std::vector<double>> out(nl+cells, std::vector<double>(np, NAN));
480 return out;
481 }
482 std::vector<std::vector<double>> out;
483
484 if (method == "bilinear") {
485 out = bilinearValues(x, y);
486 if (cells) {
487 std::vector<double> cell = cellFromXY(x, y);
488 out.push_back(cell);
489 }
490 } else {
491 std::vector<double> cell = cellFromXY(x, y);
492 out = extractCell(cell);
493 if (cells) {
494 out.push_back(cell);
495 }
496 }
497
498 return out;
499 }
500
501
extractXYFlat(const std::vector<double> & x,const std::vector<double> & y,const std::string & method,const bool & cells)502 std::vector<double> SpatRaster::extractXYFlat(const std::vector<double> &x, const std::vector<double> &y, const std::string & method, const bool &cells) {
503
504 // <layer<values>>
505 std::vector<std::vector<double>> e = extractXY(x, y, method, cells);
506 std::vector<double> out = e[0];
507 for (size_t i=1; i<e.size(); i++) {
508 out.insert(out.end(), e[i].begin(), e[i].end());
509 }
510 return out;
511 }
512
513 /*
514 std::vector<double> SpatRaster::extractXYFlat(const std::vector<double> &x, const std::vector<double> &y, const std::string & method, const bool &cells) {
515
516 unsigned nl = nlyr();
517 unsigned np = x.size();
518 if (!hasValues()) {
519 std::vector<double> out(nl * np, NAN);
520 return out;
521 }
522 std::vector<double> out;
523 if (method == "bilinear") {
524 std::vector<std::vector<double>> bil = bilinearValues(x, y);
525 if (cells) {
526 std::vector<double> cell = cellFromXY(x, y);
527 bil.push_back(cell);
528 }
529 out = flatten(bil);
530 } else {
531 std::vector<double> cell = cellFromXY(x, y);
532 if (cells) {
533 std::vector<std::vector<double>> xout;
534 xout = extractCell(cell);
535 xout.push_back(cell);
536 out = flatten(xout);
537 } else {
538 out = extractCellFlat(cell);
539 }
540 }
541 return out;
542 }
543 */
544
545
546 // <geom<layer<values>>>
extractVector(SpatVector v,bool touches,std::string method,bool cells,bool xy,bool weights,bool exact,SpatOptions & opt)547 std::vector<std::vector<std::vector<double>>> SpatRaster::extractVector(SpatVector v, bool touches, std::string method, bool cells, bool xy, bool weights, bool exact, SpatOptions &opt) {
548
549 std::string gtype = v.type();
550 if (gtype != "polygons") weights = false;
551
552 if (exact) weights = false;
553 unsigned nl = nlyr();
554 unsigned ng = v.size();
555 std::vector<std::vector<std::vector<double>>> out(ng, std::vector<std::vector<double>>(nl + cells + 2*xy + (weights || exact)));
556
557 if (!hasValues()) {
558 setError("raster has no value");
559 return out;
560 }
561 /*
562 #if GDAL_VERSION_MAJOR < 3
563 if (weights) {
564 setError("extract with weights not supported for your GDAL version");
565 return out;
566 }
567 #endif
568 */
569 std::vector<std::vector<double>> srcout;
570 if (gtype == "points") {
571 if (method != "bilinear") method = "simple";
572 SpatDataFrame vd = v.getGeometryDF();
573 if (vd.nrow() == ng) { // single point geometry
574 std::vector<double> x = vd.getD(0);
575 std::vector<double> y = vd.getD(1);
576 srcout = extractXY(x, y, method, cells);
577 for (size_t i=0; i<ng; i++) {
578 for (size_t j=0; j<nl; j++) {
579 out[i][j].push_back( srcout[j][i] );
580 }
581 if (cells) {
582 out[i][nl].push_back( srcout[nl][i] );
583 }
584 if (xy) {
585 out[i][nl+cells].push_back(x[i]);
586 out[i][nl+cells+1].push_back(y[i]);
587 }
588 }
589 } else { //multipoint
590 Rcpp::Rcout << "multipoint" << std::endl;
591 for (size_t i=0; i<ng; i++) {
592 SpatVector vv = v.subset_rows(i);
593 SpatDataFrame vd = vv.getGeometryDF();
594 std::vector<double> x = vd.getD(0);
595 std::vector<double> y = vd.getD(1);
596 //srcout = extractXY(x, y, method, cells);
597 Rcpp::Rcout << srcout.size() << " " << srcout[0].size() << std::endl;
598
599 /*
600 for (size_t j=0; j<nl; j++) {
601 out[i][j] = srcout[j];
602 }
603 if (cells) {
604 out[i][nl] = srcout[nl];
605 }
606 if (xy) {
607 out[i][nl+cells] = x;
608 out[i][nl+cells+1] = y;
609 }
610 */
611 }
612 }
613 } else {
614 SpatRaster r = geometry(1);
615 //SpatOptions opt;
616 //std::vector<double> feats(1, 1) ;
617 for (size_t i=0; i<ng; i++) {
618 SpatGeom g = v.getGeom(i);
619 SpatVector p(g);
620 p.srs = v.srs;
621 std::vector<double> cell, wgt;
622 if (weights) {
623 rasterizeCellsWeights(cell, wgt, p, opt);
624 } else if (exact) {
625 rasterizeCellsExact(cell, wgt, p, opt);
626 } else {
627 cell = rasterizeCells(p, touches, opt);
628 }
629 srcout = extractCell(cell);
630 for (size_t j=0; j<nl; j++) {
631 out[i][j] = srcout[j];
632 }
633 if (cells) {
634 out[i][nl] = cell;
635 }
636 if (xy) {
637 std::vector<std::vector<double>> crds = xyFromCell(cell);
638 out[i][nl+cells] = crds[0];
639 out[i][nl+cells+1] = crds[1];
640 }
641 if (weights || exact) {
642 out[i][nl + cells + 2*xy] = wgt;
643 }
644 }
645 }
646 return out;
647 }
648
649
extractVectorFlat(SpatVector v,bool touches,std::string method,bool cells,bool xy,bool weights,bool exact,SpatOptions & opt)650 std::vector<double> SpatRaster::extractVectorFlat(SpatVector v, bool touches, std::string method, bool cells, bool xy, bool weights, bool exact, SpatOptions &opt) {
651
652 std::vector<double> flat;
653 std::string gtype = v.type();
654 if (gtype != "polygons") weights = false;
655
656 unsigned nl = nlyr();
657 unsigned ng = v.size();
658
659 if (!hasValues()) {
660 setError("raster has no value");
661 return flat;
662 }
663 /*
664 #if GDAL_VERSION_MAJOR < 3
665 if (weights) {
666 setError("extract with weights not supported for your GDAL version");
667 return out;
668 }
669 #endif
670 */
671
672 std::vector<std::vector<std::vector<double>>> out;
673 if (gtype != "points") {
674 out.resize(ng, std::vector<std::vector<double>>(nl + cells + 2*xy + (weights||exact)));
675 }
676 std::vector<std::vector<double>> srcout;
677 if (gtype == "points") {
678 if (method != "bilinear") method = "simple";
679 SpatDataFrame vd = v.getGeometryDF();
680 //if (vd.nrow() == ng) { // single point geometry
681 std::vector<double> x = vd.getD(0);
682 std::vector<double> y = vd.getD(1);
683 if (!cells & !xy & !weights) {
684 return( extractXYFlat(x, y, method, cells));
685 } else {
686 srcout = extractXY(x, y, method, cells);
687 nl += cells;
688 flat.reserve(ng * nl);
689 for (size_t i=0; i<ng; i++) {
690 //flat.push_back( i+1 );//no id for points
691 for (size_t j=0; j<nl; j++) {
692 flat.push_back( srcout[j][i] );
693 }
694 if (xy) {
695 flat.push_back(x[i]);
696 flat.push_back(y[i]);
697 }
698 }
699 }
700 return flat;
701 /*
702 } else { // multipoint
703 Rcpp::Rcout << "multipoint" << std::endl;
704 std::vector<double> x = vd.getD(0);
705 std::vector<double> y = vd.getD(1);
706 if (!cells & !xy & !weights) {
707 return( extractXYFlat(x, y, method, cells));
708 }
709 for (size_t i=0; i<ng; i++) {
710 SpatVector vv = v.subset_rows(i);
711 SpatDataFrame vd = vv.getGeometryDF();
712 std::vector<double> x = vd.getD(0);
713 std::vector<double> y = vd.getD(1);
714 srcout = extractXY(x, y, method, cells);
715 Rcpp::Rcout << srcout.size() << " " << srcout[0].size();
716
717 out.push_back(srcout);
718
719 if (cells) {
720 out[i][nl] = srcout[nl];
721 }
722 if (xy) {
723 out[i][nl+cells] = x;
724 out[i][nl+cells+1] = y;
725 }
726
727
728 }
729
730 }
731 */
732 } else {
733 SpatRaster r = geometry(1);
734 //std::vector<double> feats(1, 1) ;
735 for (size_t i=0; i<ng; i++) {
736 SpatGeom g = v.getGeom(i);
737 SpatVector p(g);
738 p.srs = v.srs;
739 std::vector<double> cell, wgt;
740 if (weights) {
741 rasterizeCellsWeights(cell, wgt, p, opt);
742 } else if (exact) {
743 rasterizeCellsExact(cell, wgt, p, opt);
744 } else {
745 cell = rasterizeCells(p, touches, opt);
746 }
747 srcout = extractCell(cell);
748 for (size_t j=0; j<nl; j++) {
749 out[i][j] = srcout[j];
750 }
751 if (cells) {
752 out[i][nl] = cell;
753 }
754 if (xy) {
755 std::vector<std::vector<double>> crds = xyFromCell(cell);
756 out[i][nl+cells] = crds[0];
757 out[i][nl+cells+1] = crds[1];
758 }
759 if (weights || exact) {
760 out[i][nl + cells + 2*xy] = wgt;
761 }
762 }
763 }
764
765 size_t fsize = 0;
766 for (size_t i=0; i<out.size(); i++) { // geoms
767 fsize += (out[i].size()+1) * nl;
768 }
769 flat.reserve(fsize);
770
771 for (size_t i=0; i<out.size(); i++) { // geoms
772 for (size_t j=0; j<out[i][0].size(); j++) { // cells
773 flat.push_back(i+1);
774 for (size_t k=0; k<out[i].size(); k++) { // layers
775 flat.push_back(out[i][k][j]);
776 }
777 }
778 }
779 return flat;
780 }
781
782
783
extractCell(std::vector<double> & cell)784 std::vector<std::vector<double>> SpatRaster::extractCell(std::vector<double> &cell) {
785
786 std::vector<double> wcell;
787 std::vector<std::vector<int_64>> rc, wrc;
788 rc = rowColFromCell(cell);
789
790 size_t n = cell.size();
791 std::vector<std::vector<double>> out(nlyr(), std::vector<double>(n, NAN));
792 if (!hasValues()) return out;
793
794 unsigned ns = nsrc();
795 unsigned lyr = 0;
796 size_t nc;
797 for (size_t src=0; src<ns; src++) {
798 unsigned slyrs = source[src].layers.size();
799 bool win = source[src].hasWindow;
800 if (win) {
801 nc = source[src].window.full_ncol * source[src].window.full_nrow;
802 wrc = rc;
803 wcell.reserve(cell.size());
804 for (size_t i=0; i<cell.size(); i++) {
805 if ((wrc[0][i] < 0) || (wrc[1][i] <0)) {
806 wcell.push_back( NAN );
807 } else {
808 wrc[0][i] += source[src].window.off_row;
809 wrc[1][i] += source[src].window.off_col;
810 wcell.push_back( wrc[0][i] * source[src].window.full_ncol + wrc[1][i] );
811 }
812 }
813 } else {
814 nc = ncell();
815 }
816 if (source[src].memory) {
817 for (size_t i=0; i<slyrs; i++) {
818 size_t j = i * nc;
819 if (win) {
820 for (size_t k=0; k<n; k++) {
821 if (!is_NA(wcell[k]) && wcell[k] >= 0 && wcell[k] < nc) {
822 out[lyr][k] = source[src].values[j + wcell[k]];
823 }
824 }
825 } else {
826 for (size_t k=0; k<n; k++) {
827 if (!is_NA(cell[k]) && cell[k] >= 0 && cell[k] < nc) {
828 out[lyr][k] = source[src].values[j + cell[k]];
829 }
830 }
831 }
832 lyr++;
833 }
834 } else {
835 std::vector<std::vector<double>> srcout;
836 //if (source[0].driver == "raster") {
837 // srcout = readCellsBinary(src, cell);
838 //} else {
839 #ifdef useGDAL
840 if (win) {
841 srcout = readRowColGDAL(src, wrc[0], wrc[1]);
842 } else {
843 srcout = readRowColGDAL(src, rc[0], rc[1]);
844 }
845 #endif
846 if (hasError()) return out;
847 //}
848 for (size_t i=0; i<slyrs; i++) {
849 out[lyr] = srcout[i];
850 lyr++;
851 }
852 }
853 }
854 return out;
855 }
856
857
858
859
860
extractCellFlat(std::vector<double> & cell)861 std::vector<double> SpatRaster::extractCellFlat(std::vector<double> &cell) {
862
863 std::vector<double> wcell;
864 std::vector<std::vector<int_64>> rc, wrc;
865 rc = rowColFromCell(cell);
866
867 size_t n = cell.size();
868 std::vector<double> out(nlyr() * n, NAN);
869
870 unsigned ns = nsrc();
871 unsigned lyr = 0;
872 size_t nc;
873 size_t off = 0;
874 for (size_t src=0; src<ns; src++) {
875 unsigned slyrs = source[src].layers.size();
876 bool win = source[src].hasWindow;
877 if (win) {
878 nc = source[src].window.full_ncol * source[src].window.full_nrow;
879 wrc = rc;
880 wcell.reserve(cell.size());
881 for (size_t i=0; i<cell.size(); i++) {
882 if ((wrc[0][i] < 0) || (wrc[1][i] <0)) {
883 wcell.push_back( NAN );
884 } else {
885 wrc[0][i] += source[src].window.off_row;
886 wrc[1][i] += source[src].window.off_col;
887 wcell.push_back( wrc[0][i] * source[src].window.full_ncol + wrc[1][i] );
888 }
889 }
890 } else {
891 nc = ncell();
892 }
893
894 if (source[src].memory) {
895 for (size_t i=0; i<slyrs; i++) {
896 size_t j = i * nc;
897 size_t off2 = off + i*n;
898 if (win) {
899 for (size_t k=0; k<n; k++) {
900 if (!is_NA(wcell[k]) && wcell[k] >= 0 && wcell[k] < nc) {
901 out[off2+k] = source[src].values[j + wcell[k]] ;
902 }
903 }
904 } else {
905 for (size_t k=0; k<n; k++) {
906 if (!is_NA(cell[k]) && cell[k] >= 0 && cell[k] < nc) {
907 out[off2+k] = source[src].values[j + cell[k]];
908 }
909 }
910 }
911 lyr++;
912 }
913 } else {
914 //if (source[0].driver == "raster") {
915 // srcout = readCellsBinary(src, cell);
916 //} else {
917 #ifdef useGDAL
918 std::vector<double> g;
919 if (win) {
920 g = readRowColGDALFlat(src, wrc[0], wrc[1]);
921 } else {
922 g = readRowColGDALFlat(src, rc[0], rc[1]);
923 }
924 for (size_t i=0; i<slyrs; i++) {
925 size_t j = i * n;
926 size_t off2 = off + j;
927 for (size_t k=0; k<n; k++) {
928 out[off2+k] = g[j + k];
929 }
930 }
931 #endif
932 if (hasError()) {
933 return out;
934 }
935 }
936 off += source[src].nlyr;
937 }
938 return out;
939 }
940
941
vectCells(SpatVector v,bool touches,std::string method,bool weights,bool exact,SpatOptions & opt)942 std::vector<double> SpatRaster::vectCells(SpatVector v, bool touches, std::string method, bool weights, bool exact, SpatOptions &opt) {
943
944 std::string gtype = v.type();
945 if (gtype != "polygons") weights = false;
946 std::vector<double> out, cells, wghts;
947 if (gtype == "points") {
948 SpatDataFrame vd = v.getGeometryDF();
949 std::vector<long> id = vd.getI(0);
950 if (method == "bilinear") {
951 return bilinearCells(vd.getD(0), vd.getD(1));
952 } else {
953 return cellFromXY(vd.getD(0), vd.getD(1));
954 //cells = cellFromXY(vd.getD(0), vd.getD(1));
955 //out.insert(out.end(), id.begin(), id.end());
956 //out.insert(out.end(), cells.begin(), cells.end());
957 }
958 } else {
959 unsigned ng = v.size();
960 SpatRaster r = geometry(1);
961 std::vector<double> feats(1, 1) ;
962 for (size_t i=0; i<ng; i++) {
963 SpatGeom g = v.getGeom(i);
964 SpatVector p(g);
965 p.srs = v.srs;
966 if (weights) {
967 std::vector<double> cnr, wght;
968 rasterizeCellsWeights(cnr, wght, p, opt);
969 std::vector<double> id(cnr.size(), i);
970 out.insert(out.end(), id.begin(), id.end());
971 cells.insert(cells.end(), cnr.begin(), cnr.end());
972 wghts.insert(wghts.end(), wght.begin(), wght.end());
973 } else if (exact) {
974 std::vector<double> cnr, wght;
975 rasterizeCellsExact(cnr, wght, p, opt);
976 std::vector<double> id(cnr.size(), i);
977 out.insert(out.end(), id.begin(), id.end());
978 cells.insert(cells.end(), cnr.begin(), cnr.end());
979 wghts.insert(wghts.end(), wght.begin(), wght.end());
980 } else {
981 std::vector<double> geomc = rasterizeCells(p, touches, opt);
982 std::vector<double> id(geomc.size(), i);
983 out.insert(out.end(), id.begin(), id.end());
984 cells.insert(cells.end(), geomc.begin(), geomc.end());
985 }
986 }
987 if (weights || exact) {
988 out.insert(out.end(), cells.begin(), cells.end());
989 out.insert(out.end(), wghts.begin(), wghts.end());
990 } else {
991 out.insert(out.end(), cells.begin(), cells.end());
992 }
993 }
994 return out;
995 }
996
997
extCells(SpatExtent ext)998 std::vector<double> SpatRaster::extCells(SpatExtent ext) {
999
1000 std::vector<double> out;
1001 ext = align(ext, "near");
1002 ext.intersect(getExtent());
1003 if (!ext.valid()) {
1004 return(out);
1005 }
1006 double resx = xres() / 2;
1007 double resy = yres() / 2;
1008 std::vector<double> e = ext.asVector();
1009 e[0] += resx;
1010 e[1] -= resx;
1011 e[2] += resy;
1012 e[3] -= resy;
1013 std::vector<double> ex = {e[0], e[1]};
1014 std::vector<double> ey = {e[3], e[2]};
1015 std::vector<int_64> r = rowFromY(ey);
1016 std::vector<int_64> c = colFromX(ex);
1017 int_64 nc = ncol();
1018 out.reserve((r[1]-r[0]) * (c[1]-c[0]));
1019 for (int_64 i=r[0]; i <= r[1]; i++) {
1020 for (int_64 j=c[0]; j <= c[1]; j++) {
1021 out.push_back(i*nc+j);
1022 }
1023 }
1024 return out;
1025 }
1026
1027
1028
extractXY(std::vector<double> & x,std::vector<double> & y,std::string method)1029 std::vector<std::vector<std::vector<double>>> SpatRasterStack::extractXY(std::vector<double> &x, std::vector<double> &y, std::string method) {
1030 unsigned ns = nsds();
1031 std::vector<std::vector<std::vector<double>>> out(ns);
1032 bool cells = false;
1033 for (size_t i=0; i<ns; i++) {
1034 SpatRaster r = getsds(i);
1035 out[i] = r.extractXY(x, y, method, cells);
1036 }
1037 return out;
1038 }
1039
extractCell(std::vector<double> & cell)1040 std::vector<std::vector<std::vector<double>>> SpatRasterStack::extractCell(std::vector<double> &cell) {
1041 unsigned ns = nsds();
1042 std::vector<std::vector<std::vector<double>>> out(ns);
1043 for (size_t i=0; i<ns; i++) {
1044 SpatRaster r = getsds(i);
1045 out[i] = r.extractCell(cell);
1046 }
1047 return out;
1048 }
1049
1050
1051 // this is rather inefficient (repeated rasterization)
extractVector(SpatVector v,bool touches,std::string method,SpatOptions & opt)1052 std::vector<std::vector<std::vector<std::vector<double>>>> SpatRasterStack::extractVector(SpatVector v, bool touches, std::string method, SpatOptions &opt) {
1053 unsigned ns = nsds();
1054 std::vector<std::vector<std::vector<std::vector<double>>>> out(ns);
1055 for (size_t i=0; i<ns; i++) {
1056 SpatRaster r = getsds(i);
1057 out[i] = r.extractVector(v, touches, method, false, false, false, false, opt);
1058 }
1059 return out;
1060 }
1061
1062
1063
1064 /*
1065 } else if (method == "oldbilinear") {
1066
1067 // this is much too slow
1068 SpatExtent extent = getExtent();
1069
1070 double ymax = extent.ymax;
1071 double xmin = extent.xmin;
1072 double yrs = yres();
1073 double xrs = xres();
1074
1075 //SpatOptions opt;
1076 SpatRaster g = geometry();
1077 std::vector<unsigned> f = {2,2};
1078 SpatRaster gd = g.disaggregate(f, opt);
1079
1080 double dyrs = gd.yres();
1081 double dxrs = gd.xres();
1082 std::vector<double> d, cells(4);
1083
1084 std::vector<std::vector<double> > cxy;
1085
1086 std::vector<double> rc(4);
1087 unsigned nr = nrow();
1088 unsigned nc = ncol();
1089 unsigned mnr = nr-1;
1090 unsigned mnc = nc-1;
1091
1092 // needs row-wise adjustment for lonlat
1093 double dxdy = xres() * yres();
1094
1095 for (size_t i=0; i<n; i++) {
1096 long row_d = ((ymax - y[i]) / dyrs);
1097 long col_d = ((x[i] - xmin) / dxrs);
1098 unsigned rq = row_d % 2;
1099 unsigned cq = col_d % 2;
1100
1101 double row1 = std::floor((ymax - y[i]) / yrs);
1102 double col1 = std::floor((x[i] - xmin) / xrs);
1103 if ((row1 < 0) || (row1 > mnr)) { continue; }
1104
1105 double row2 = (rq == 0) ? row1-1 : row1+1;
1106 row2 = row2 < 0 ? row1+1 : row2==nr ? row1-1 : row2;
1107 double col2;
1108 if (globalLonLat) {
1109 if ((col1 < -1) | (col1 > nc)) { continue; }
1110 col1 = col1 < 0 ? mnc : col1 > mnc ? 0 : col1;
1111 col2 = (cq == 0) ? col1-1 : col1 + 1;
1112 col2 = col2 < 0 ? mnc : col2 > mnc ? 0 : col2;
1113 } else {
1114 if ((col1 < 0) | (col1 > mnc)) { continue; }
1115 col2 = (cq == 0) ? col1-1 : col1 + 1;
1116 col2 = col2 < 0 ? col1+1 : col2 == nc ? col1-1 : col2;
1117 }
1118 cells[0] = nc * row1 + col1;
1119 cells[1] = nc * row1 + col2;
1120 cells[2] = nc * row2 + col1;
1121 cells[3] = nc * row2 + col2;
1122 std::sort(cells.begin(), cells.end());
1123 std::vector<std::vector<double>> xy = xyFromCell(cells);
1124 std::vector<std::vector<double>> v = extractCell(cells);
1125 std::vector<double> e = {xy[0][0], xy[0][1], xy[1][2], xy[1][0]};
1126 for (size_t j=0; j<nl; j++) {
1127 out[j][i] = bilinear(v[j], e, dxdy, x[i], y[i]);
1128 }
1129 }
1130 }
1131 */
1132
1133
1134 /*
1135 std::vector<double> SpatRaster::extractCell(std::vector<double> &cell) {
1136
1137 unsigned n = cell.size();
1138 unsigned nc = ncell();
1139
1140 std::vector<double> out;
1141 if (!hasValues()) {
1142 out = std::vector<double>(n * nlyr(), NAN)
1143 return out;
1144 }
1145
1146 unsigned ns = nsrc();
1147 for (size_t src=0; src<ns; src++) {
1148 unsigned slyrs = source[src].layers.size();
1149 std::vector<double> srcout;
1150 if (source[src].memory) {
1151 srcout = std::vector<double>(n * slyrs, NAN)
1152 std::vector<size_t> off1(slyrs);
1153 std::vector<size_t> off2(slyrs);
1154 for (size_t i=0; i<slyrs; i++) {
1155 off1[i] = i * n;
1156 off2[i] = i * nc;
1157 }
1158 for (size_t i=0; i<n; i++) {
1159 if (!is_NA(cell[i]) && cell[i] >= 0 && cell[i] < nc) {
1160 for (size_t j=0; j<slyrs; j++) {
1161 out[off1[j]+i] = source[src].values[off2[j] + cell[i]];
1162 }
1163 }
1164 }
1165 } else {
1166 #ifdef useGDAL
1167 std::vector<std::vector<int_64>> rc = rowColFromCell(cell);
1168 srcout = readRowColGDAL(src, rc[0], rc[1]);
1169 #endif
1170 if (hasError()) return out;
1171 //}
1172 }
1173 out.insert(out.end(), srcout.begin(), srcout.end());
1174 }
1175 return out;
1176 }
1177 */
1178