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 "gdalwarper.h"
19 #include "ogr_spatialref.h"
20 #include "gdal_alg.h"
21 #include "ogrsf_frmts.h"
22
23 #include "spatRaster.h"
24 #include "string_utils.h"
25 #include "file_utils.h"
26 #include "vecmath.h"
27
28 #include "crs.h"
29 #include "gdalio.h"
30 #include "recycle.h"
31
32
33 //#include <vector>
34 //#include "vecmath.h"
35
36
dense_extent()37 SpatVector SpatRaster::dense_extent() {
38
39 std::vector<int_64> rows, cols;
40 if (nrow() < 51) {
41 rows.resize(nrow());
42 std::iota(rows.begin(), rows.end(), 0);
43 } else {
44 rows = seq_steps((int_64) 0, (int_64) nrow()-1, 50);
45 }
46 if (ncol() < 51) {
47 cols.resize(nrow());
48 std::iota(cols.begin(), cols.end(), 0);
49 } else {
50 cols = seq_steps((int_64) 0, (int_64) ncol()-1, 50);
51 }
52
53 SpatExtent e = getExtent();
54
55 std::vector<double> xcol = xFromCol(cols) ;
56 std::vector<double> yrow = yFromRow(rows) ;
57 yrow.insert(yrow.begin(), e.ymax);
58 yrow.push_back(e.ymin);
59
60 std::vector<double> y0(xcol.size(), e.ymin);
61 std::vector<double> y1(xcol.size(), e.ymax);
62 std::vector<double> x0(yrow.size(), e.xmin);
63 std::vector<double> x1(yrow.size(), e.xmax);
64
65 std::vector<double> x = x0;
66 std::vector<double> y = yrow;
67 x.insert(x.end(), xcol.begin(), xcol.end());
68 y.insert(y.end(), y0.begin(), y0.end());
69
70 std::reverse(yrow.begin(), yrow.end());
71 std::reverse(xcol.begin(), xcol.end());
72
73 x.insert(x.end(), x1.begin(), x1.end());
74 y.insert(y.end(), yrow.begin(), yrow.end() );
75 x.insert(x.end(), xcol.begin(), xcol.end());
76 y.insert(y.end(), y1.begin(), y1.end());
77
78 x.push_back(e.xmin);
79 y.push_back(e.ymax);
80
81 SpatVector v(x, y, polygons, getSRS("wkt"));
82
83 return v;
84 }
85
86 #if GDAL_VERSION_MAJOR <= 2 && GDAL_VERSION_MINOR < 2
87
warper(SpatRaster x,std::string crs,std::string method,bool mask,bool align,SpatOptions & opt)88 SpatRaster SpatRaster::warper(SpatRaster x, std::string crs, std::string method, bool mask, bool align, SpatOptions &opt) {
89 SpatRaster out;
90 out.setError("Not supported for this old version of GDAL");
91 return(out);
92 }
93
94
95 #else
96
97
98
get_output_bounds(const GDALDatasetH & hSrcDS,std::string srccrs,const std::string dstcrs,SpatRaster & r)99 bool get_output_bounds(const GDALDatasetH &hSrcDS, std::string srccrs, const std::string dstcrs, SpatRaster &r) {
100
101 if ( hSrcDS == NULL ) {
102 r.setError("data source is NULL");
103 return false;
104 }
105
106 // Get Source coordinate system.
107 // const char *pszSrcWKT = GDALGetProjectionRef( hSrcDS );
108 const char *pszSrcWKT = srccrs.c_str();
109 if ( pszSrcWKT == NULL || strlen(pszSrcWKT) == 0 ) {
110 r.setError("data source has no WKT");
111 return false;
112 }
113
114 OGRSpatialReference* oSRS = new OGRSpatialReference;
115 std::string msg = "";
116 if (is_ogr_error(oSRS->SetFromUserInput( dstcrs.c_str() ), msg)) {
117 r.setError(msg);
118 return false;
119 };
120
121 char *pszDstWKT = NULL;
122 oSRS->exportToWkt( &pszDstWKT );
123
124 // Create a transformer that maps from source pixel/line coordinates
125 // to destination georeferenced coordinates (not destination
126 // pixel line). We do that by omitting the destination dataset
127 // handle (setting it to NULL).
128 void *hTransformArg;
129 hTransformArg =
130 GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT, FALSE, 0, 1 );
131 if (hTransformArg == NULL ) {
132 r.setError("cannot create TranformArg");
133 return false;
134 }
135 CPLFree(pszDstWKT);
136 delete oSRS;
137
138 double adfDstGeoTransform[6];
139 int nPixels=0, nLines=0;
140 CPLErr eErr = GDALSuggestedWarpOutput( hSrcDS, GDALGenImgProjTransform,
141 hTransformArg, adfDstGeoTransform, &nPixels, &nLines );
142
143 GDALDestroyGenImgProjTransformer( hTransformArg );
144 if ( eErr != CE_None ) {
145 r.setError("cannot create warp output");
146 return false;
147 }
148
149
150 r.source[0].ncol = nPixels;
151 r.source[0].nrow = nLines;
152
153 r.source[0].extent.xmin = adfDstGeoTransform[0]; /* left x */
154 /* w-e pixel resolution */
155 r.source[0].extent.xmax = r.source[0].extent.xmin + adfDstGeoTransform[1] * nPixels;
156 r.source[0].extent.ymax = adfDstGeoTransform[3]; // top y
157 r.source[0].extent.ymin = r.source[0].extent.ymax + nLines * adfDstGeoTransform[5];
158
159 r.setSRS({dstcrs});
160
161 return true;
162 }
163
164 /*
165 // Create output with same datatype as first input band.
166 GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));
167 GDALDataType eDT;
168 getGDALDataType(datatype, eDT);
169
170 // Create the output DS.
171
172 GDALDriverH hDriver = GDALGetDriverByName( driver.c_str() );
173 if ( hDriver == NULL ) {
174 msg = "empty driver";
175 return false;
176 }
177 if (driver == "MEM") {
178 hDstDS = GDALCreate( hDriver, "", nPixels, nLines, nlyrs, eDT, NULL );
179 } else {
180 hDstDS = GDALCreate( hDriver, filename.c_str(), nPixels, nLines, nlyrs, eDT, NULL );
181 }
182 if ( hDstDS == NULL ) {
183 msg = "cannot create output dataset";
184 return false;
185 }
186
187 // Write out the projection definition.
188 GDALSetProjection( hDstDS, pszDstWKT );
189 GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
190
191 // Copy the color table, if required.
192 GDALColorTableH hCT;
193 hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
194 if( hCT != NULL )
195 GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
196
197 CPLFree(pszDstWKT);
198 delete oSRS;
199
200 return true;
201 }
202 */
203
getAlgo(std::string m)204 GDALResampleAlg getAlgo(std::string m) {
205 GDALResampleAlg alg;
206 if ( m == "near" ) {
207 alg = GRA_NearestNeighbour;
208 } else if (m=="bilinear") {
209 alg = GRA_Bilinear;
210 } else if (m=="cubic") {
211 alg = GRA_Cubic;
212 } else if (m=="cubicspline") {
213 alg = GRA_CubicSpline;
214 } else if (m=="lanczos") {
215 alg = GRA_Lanczos;
216 } else if (m=="mean") {
217 alg = GRA_Average;
218 // } else if (m=="sum") {
219 // alg = GRA_Sum;
220 } else if (m=="mode") {
221 alg = GRA_Mode;
222 } else if (m=="max") {
223 alg = GRA_Max;
224 } else if (m=="min") {
225 alg = GRA_Min;
226 } else if (m=="median") {
227 alg = GRA_Med;
228 } else if (m=="q1") {
229 alg = GRA_Q1;
230 } else if (m=="q3") {
231 alg = GRA_Q3;
232 } else {
233 alg = GRA_NearestNeighbour;
234 }
235 return alg;
236 }
237
238
gdal_warper(GDALDatasetH & hSrcDS,GDALDatasetH & hDstDS,std::vector<unsigned> srcbands,std::vector<unsigned> dstbands,std::string method,std::string srccrs,std::string msg,bool verbose)239 bool gdal_warper(GDALDatasetH &hSrcDS, GDALDatasetH &hDstDS, std::vector<unsigned> srcbands, std::vector<unsigned> dstbands, std::string method, std::string srccrs, std::string msg, bool verbose) {
240
241 if (srcbands.size() != dstbands.size()) {
242 msg = "number of source bands must match number of dest bands";
243 return false;
244 }
245 int nbands = srcbands.size();
246
247 GDALResampleAlg a = getAlgo(method);
248
249 // Setup warp options.
250 GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
251 psWarpOptions->hSrcDS = hSrcDS;
252 psWarpOptions->hDstDS = hDstDS;
253 psWarpOptions->eResampleAlg = a;
254 psWarpOptions->nBandCount = nbands;
255 psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * nbands );
256 psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * nbands );
257 psWarpOptions->padfSrcNoDataReal = (double *) CPLMalloc(sizeof(double) * nbands );
258 psWarpOptions->padfDstNoDataReal = (double *) CPLMalloc(sizeof(double) * nbands );
259 psWarpOptions->padfSrcNoDataImag = (double *) CPLMalloc(sizeof(double) * nbands );
260 psWarpOptions->padfDstNoDataImag = (double *) CPLMalloc(sizeof(double) * nbands );
261
262 GDALRasterBandH hBand;
263 int hasNA;
264 for (int i=0; i<nbands; i++) {
265 psWarpOptions->panSrcBands[i] = (int) srcbands[i]+1;
266 psWarpOptions->panDstBands[i] = (int) dstbands[i]+1;
267
268 hBand = GDALGetRasterBand(hSrcDS, srcbands[i]+1);
269 double naflag = GDALGetRasterNoDataValue(hBand, &hasNA);
270 if (verbose && i == 0) {
271 #ifdef useRcpp
272 std::string hna = hasNA ? "true" : "false";
273 Rcpp::Rcout << "hasNA : " << hna << std::endl;
274 Rcpp::Rcout << "NA flag : " << naflag << std::endl;
275 #endif
276 }
277 if (hasNA) {
278 psWarpOptions->padfSrcNoDataReal[i] = naflag;
279 psWarpOptions->padfDstNoDataReal[i] = naflag;
280 hBand = GDALGetRasterBand(hDstDS, dstbands[i]+1);
281 GDALSetRasterNoDataValue(hBand, naflag);
282 } else {
283 psWarpOptions->padfSrcNoDataReal[i] = NAN;
284 psWarpOptions->padfDstNoDataReal[i] = NAN;
285 }
286 psWarpOptions->padfSrcNoDataImag[i] = 0;
287 psWarpOptions->padfDstNoDataImag[i] = 0;
288 }
289
290 //psWarpOptions->pfnProgress = GDALTermProgress;
291
292 psWarpOptions->papszWarpOptions =
293 CSLSetNameValue( psWarpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA");
294 psWarpOptions->papszWarpOptions =
295 CSLSetNameValue( psWarpOptions->papszWarpOptions, "WRITE_FLUSH", "YES");
296
297
298 psWarpOptions->pTransformerArg =
299 GDALCreateGenImgProjTransformer( hSrcDS, srccrs.c_str(),
300 hDstDS, GDALGetProjectionRef(hDstDS),
301 FALSE, 0.0, 1 );
302 psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
303
304 GDALWarpOperation oOperation;
305 if (oOperation.Initialize( psWarpOptions ) != CE_None) {
306 return false;
307 }
308 oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) );
309 GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
310 GDALDestroyWarpOptions( psWarpOptions );
311 return true;
312 }
313
314
is_valid_warp_method(const std::string & method)315 bool is_valid_warp_method(const std::string &method) {
316 std::vector<std::string> m { "near", "bilinear", "cubic", "cubicspline", "lanczos", "average", "mode", "max", "min", "med", "q1", "q3", "sum" };
317 return (std::find(m.begin(), m.end(), method) != m.end());
318 }
319
320
321 /*
322 SpatRaster SpatRaster::old_warper(SpatRaster x, std::string crs, std::string method, bool mask, SpatOptions &opt) {
323
324 SpatRaster out = x.geometry(nlyr(), false, false);
325 out.setNames(getNames());
326 if (method == "near") {
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 out.rgb = rgb;
332 out.rgblyrs = rgblyrs;
333 }
334 if (hasTime()) {
335 out.source[0].hasTime = true;
336 out.source[0].timestep = getTimeStep();
337 out.source[0].time = getTime();
338 }
339
340
341 bool use_crs = crs != "";
342 std::string filename = opt.get_filename();
343 if ((!use_crs) & (!hasValues())) {
344 if (filename != "") {
345 out.addWarning("raster has no values, not writing to file");
346 }
347 return out;
348 }
349
350 if (!is_valid_warp_method(method)) {
351 out.setError("not a valid warp method");
352 return out;
353 }
354 lrtrim(crs);
355 std::string errmsg;
356 SpatOptions mopt;
357 if (mask) {
358 mopt = opt;
359 opt = SpatOptions(opt);
360 }
361 filename = opt.get_filename();
362
363 std::string srccrs = getSRS("wkt");
364 if (srccrs == "") {
365 out.setError("input raster CRS not set");
366 return out;
367 }
368
369
370 if (filename == "") {
371 if (!canProcessInMemory(opt) || !out.canProcessInMemory(opt)) {
372 filename = tempFile(opt.get_tempdir(), ".tif");
373 }
374 } else {
375 if (!can_write(filename, opt.get_overwrite(), errmsg)) {
376 out.setError(errmsg);
377 return out;
378 }
379 }
380
381 if (opt.names.size() == out.nlyr()) {
382 out.setNames(opt.names);
383 }
384
385 if (!hasValues()) filename = ""; // for crs case
386 std::string driver = filename == "" ? "MEM" : "GTiff";
387
388 GDALDatasetH hSrcDS, hDstDS;
389 size_t ns = nsrc();
390 int bandstart = 0;
391
392 for (size_t i=0; i<ns; i++) {
393
394 if (!open_gdal(hSrcDS, i, opt)) {
395 out.setError("cannot create dataset from source");
396 return out;
397 }
398
399 // create dest source, only once
400 if (i==0) {
401 if (use_crs) {
402 if (!get_output_bounds(hSrcDS, srccrs, crs, out)) {
403 GDALClose( hSrcDS );
404 return out;
405 }
406 }
407 if (!hasValues()) {
408 GDALClose( hSrcDS );
409 return out;
410 }
411 if (!out.create_gdalDS(hDstDS, filename, driver, false, NAN, source[i].has_scale_offset, source[i].scale, source[i].offset, opt)) {
412 GDALClose( hSrcDS );
413 //GDALClose( hDstDS );
414 return out;
415 }
416 }
417 std::vector<unsigned> srcbands = source[i].layers;
418 std::vector<unsigned> dstbands(srcbands.size());
419 std::iota (dstbands.begin(), dstbands.end(), bandstart);
420 bandstart += dstbands.size();
421
422 bool success = gdal_warper(hSrcDS, hDstDS, srcbands, dstbands, method, srccrs, errmsg, opt.get_verbose());
423
424 GDALClose( hSrcDS );
425 if (!success) {
426 GDALClose( hDstDS );
427 out.setError(errmsg);
428 return out;
429 }
430 }
431
432 if (driver == "MEM") {
433 bool test = out.from_gdalMEM(hDstDS, use_crs, true);
434 GDALClose( hDstDS );
435 if (!test) {
436 out.setError("cannot do this transformation (warp)");
437 return out;
438 }
439 } else {
440 std::vector<std::string> nms = getNames();
441 for (size_t i=0; i < nlyr(); i++) {
442 GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, i+1);
443 double adfMinMax[2];
444 bool approx = ncell() > 10e+8;
445 GDALComputeRasterMinMax(hBand, approx, adfMinMax);
446 GDALSetRasterStatistics(hBand, adfMinMax[0], adfMinMax[1], NAN, NAN);
447 GDALSetDescription(hBand, nms[i].c_str());
448 }
449 GDALClose( hDstDS );
450 out = SpatRaster(filename, {-1}, {""});
451 }
452
453 if (mask) {
454 SpatVector v = dense_extent();
455 v = v.project(out.getSRS("wkt"));
456 out = out.mask(v, false, NAN, true, mopt);
457 }
458
459 return out;
460 }
461 */
462
463
464
warper(SpatRaster x,std::string crs,std::string method,bool mask,bool align,SpatOptions & opt)465 SpatRaster SpatRaster::warper(SpatRaster x, std::string crs, std::string method, bool mask, bool align, SpatOptions &opt) {
466
467
468 SpatRaster out = x.geometry(nlyr(), false, false);
469 out.setNames(getNames());
470 if (method == "near") {
471 out.source[0].hasColors = hasColors();
472 out.source[0].cols = getColors();
473 out.source[0].hasCategories = hasCategories();
474 out.source[0].cats = getCategories();
475 out.rgb = rgb;
476 out.rgblyrs = rgblyrs;
477 }
478 if (hasTime()) {
479 out.source[0].hasTime = true;
480 out.source[0].timestep = getTimeStep();
481 out.source[0].time = getTime();
482 }
483
484 bool use_crs = crs != "";
485 if (use_crs) {
486 align = false;
487 }
488 if (align) {
489 crs = out.getSRS("wkt");
490 }
491 if ((!use_crs) & (!hasValues())) {
492 std::string fname = opt.get_filename();
493 if (fname != "") {
494 out.addWarning("raster has no values, not writing to file");
495 }
496 return out;
497 }
498
499
500 if (!is_valid_warp_method(method)) {
501 out.setError("not a valid warp method");
502 return out;
503 }
504
505 std::string srccrs = getSRS("wkt");
506 if (srccrs == "") {
507 out.setError("input raster CRS not set");
508 return out;
509 }
510
511 lrtrim(crs);
512 SpatOptions sopt(opt);
513 if (use_crs || align) {
514 GDALDatasetH hSrcDS;
515 if (!open_gdal(hSrcDS, 0, true, sopt)) {
516 out.setError("cannot create dataset from source");
517 return out;
518 }
519 if (!get_output_bounds(hSrcDS, srccrs, crs, out)) {
520 GDALClose( hSrcDS );
521 out.setError("cannot get output boundaries");
522 return out;
523 }
524 GDALClose( hSrcDS );
525 }
526 if (align) {
527 SpatExtent e = out.getExtent();
528 e = x.align(e, "out");
529 out.setExtent(e, false);
530 std::vector<double> res = x.resolution();
531 out = out.setResolution(res[0], res[1]);
532 }
533 if (!hasValues()) {
534 return out;
535 }
536
537 SpatOptions mopt;
538 if (mask) {
539 mopt = opt;
540 opt = SpatOptions(opt);
541 }
542
543 opt.ncopies += 4;
544 if (!out.writeStart(opt)) {
545 return out;
546 }
547
548 std::string errmsg;
549 size_t ns = nsrc();
550 SpatExtent eout = out.getExtent();
551
552
553 std::vector<bool> has_so = source[0].has_scale_offset;
554 std::vector<double> scale = source[0].scale;
555 std::vector<double> offset = source[0].offset;
556 for (size_t i=1; i<ns; i++) {
557 has_so.insert(has_so.end(), source[0].has_scale_offset.begin(), source[0].has_scale_offset.end());
558 scale.insert(scale.end(), source[0].scale.begin(), source[0].scale.end());
559 offset.insert(offset.end(), source[0].offset.begin(), source[0].offset.end());
560 }
561
562 for (size_t i = 0; i < out.bs.n; i++) {
563 int bandstart = 0;
564 eout.ymax = out.yFromRow(out.bs.row[i]);
565 eout.ymin = out.yFromRow(out.bs.row[i] + out.bs.nrows[i]-1);
566 SpatRaster crop_out = out.crop(eout, "out", sopt);
567 GDALDatasetH hDstDS;
568
569 if (!crop_out.create_gdalDS(hDstDS, "", "MEM", false, NAN, has_so, scale, offset, sopt)) {
570 return crop_out;
571 }
572
573 for (size_t i=0; i<ns; i++) {
574
575 GDALDatasetH hSrcDS;
576
577 if (!open_gdal(hSrcDS, i, false, sopt)) {
578 out.setError("cannot create dataset from source");
579 if( hDstDS != NULL ) GDALClose( (GDALDatasetH) hDstDS );
580 return out;
581 }
582 std::vector<unsigned> srcbands = source[i].layers;
583 std::vector<unsigned> dstbands(srcbands.size());
584 std::iota (dstbands.begin(), dstbands.end(), bandstart);
585 bandstart += dstbands.size();
586
587
588 bool success = gdal_warper(hSrcDS, hDstDS, srcbands, dstbands, method, srccrs, errmsg, opt.get_verbose());
589 if( hSrcDS != NULL ) GDALClose( (GDALDatasetH) hSrcDS );
590 if (!success) {
591 if( hDstDS != NULL ) GDALClose( (GDALDatasetH) hDstDS );
592 out.setError(errmsg);
593 return out;
594 }
595 }
596
597 bool ok = crop_out.from_gdalMEM(hDstDS, false, true);
598
599 if( hDstDS != NULL ) GDALClose( (GDALDatasetH) hDstDS );
600 if (!ok) {
601 out.setError("cannot do this transformation (warp)");
602 return out;
603 }
604 std::vector<double> v = crop_out.getValues(-1, opt);
605 if (!out.writeValues(v, out.bs.row[i], out.bs.nrows[i], 0, out.ncol())) return out;
606 }
607 out.writeStop();
608
609 if (mask) {
610 SpatVector v = dense_extent();
611 v = v.project(out.getSRS("wkt"));
612 out = out.mask(v, false, NAN, true, mopt);
613 }
614
615 return out;
616 }
617
618
619
620 #endif
621
622
623
624
625
626
rectify(std::string method,SpatRaster aoi,unsigned useaoi,bool snap,SpatOptions & opt)627 SpatRaster SpatRaster::rectify(std::string method, SpatRaster aoi, unsigned useaoi, bool snap, SpatOptions &opt) {
628 SpatRaster out = geometry(0);
629
630 if (nsrc() > 1) {
631 out.setError("you can transform only one data source at a time");
632 return(out);
633 }
634 if (!source[0].rotated) {
635 out.setError("this source is not rotated");
636 return(out);
637 }
638 GDALDataset *poDataset = openGDAL(source[0].filename, GDAL_OF_RASTER | GDAL_OF_READONLY, source[0].open_ops);
639
640 if( poDataset == NULL ) {
641 setError("cannot read from " + source[0].filename);
642 return out;
643 }
644 double gt[6];
645 if( poDataset->GetGeoTransform(gt) != CE_None ) {
646 out.setError("can't get geotransform");
647 GDALClose( (GDALDatasetH) poDataset );
648 return out;
649 }
650 GDALClose( (GDALDatasetH) poDataset );
651 //SpatExtent e = getExtent();
652 //std::vector<double> x = {e.xmin, e.xmin, e.xmax, e.xmax };
653 //std::vector<double> y = {e.ymin, e.ymax, e.ymin, e.ymax };
654 double nc = ncol();
655 double nr = nrow();
656 std::vector<double> x = {0, 0, nc, nc};
657 std::vector<double> y = {0, nr, 0, nr};
658 std::vector<double> xx(4);
659 std::vector<double> yy(4);
660 for (size_t i=0; i<4; i++) {
661 xx[i] = gt[0] + x[i]*gt[1] + y[i]*gt[2];
662 yy[i] = gt[3] + x[i]*gt[4] + y[i]*gt[5];
663 }
664 double xmin = vmin(xx, TRUE);
665 double xmax = vmax(xx, TRUE);
666 double ymin = vmin(yy, TRUE);
667 double ymax = vmax(yy, TRUE);
668 SpatExtent en(xmin, xmax, ymin, ymax);
669 out = out.setResolution(gt[1], -gt[5]);
670 out.setExtent(en, false, "out");
671
672 if (useaoi == 1) { // use extent
673 en = aoi.getExtent();
674 if (snap) {
675 en = out.align(en, "near");
676 out.setExtent(en, false, "near");
677 } else {
678 out.setExtent(en, false, "");
679 }
680 } else if (useaoi == 2){ // extent and resolution
681 out = aoi.geometry(0);
682 } // else { // if (useaoi == 0) // no aoi
683
684
685 out = warper(out, "", method, false, false, opt);
686
687 return(out);
688 }
689
690
691
polygonize(bool trunc,bool values,bool narm,bool aggregate,SpatOptions & opt)692 SpatVector SpatRaster::polygonize(bool trunc, bool values, bool narm, bool aggregate, SpatOptions &opt) {
693
694 SpatVector out;
695 out.srs = source[0].srs;
696 SpatOptions topt(opt);
697
698 SpatRaster tmp;
699 if (nlyr() > 1) {
700 out.addWarning("only the first layer is polygonized when 'dissolve=TRUE'");
701 tmp = subset({0}, topt);
702 } else {
703 tmp = *this;
704 }
705
706 bool usemask = false;
707 SpatRaster mask;
708 if (narm) {
709 usemask = true;
710 SpatOptions mopt(topt);
711 mopt.set_datatype("INT1U");
712 mask = tmp.isfinite(mopt);
713 } else if (trunc) {
714 tmp = tmp.math("trunc", topt);
715 trunc = false;
716 } else if (tmp.sources_from_file()) {
717 // for NAN and INT files. Should have a check for that
718 //tmp = tmp.arith(0, "+", false, topt);
719 // riskier
720 tmp.readAll();
721 }
722 if (tmp.source[0].extset) {
723 tmp = tmp.hardCopy(topt);
724 }
725
726 GDALDatasetH rstDS;
727 if (!tmp.open_gdal(rstDS, 0, false, topt)) {
728 out.setError("cannot open dataset");
729 return out;
730 }
731
732 GDALDataset *srcDS=NULL;
733
734 #if GDAL_VERSION_MAJOR <= 2 && GDAL_VERSION_MINOR <= 2
735 srcDS = (GDALDataset *) rstDS;
736 #else
737 srcDS = srcDS->FromHandle(rstDS);
738 #endif
739
740 GDALDataset *maskDS=NULL;
741 GDALDatasetH rstMask;
742 if (usemask) {
743 if (!mask.open_gdal(rstMask, 0, false, opt)) {
744 out.setError("cannot open dataset");
745 return out;
746 }
747 #if GDAL_VERSION_MAJOR <= 2 && GDAL_VERSION_MINOR <= 2
748 maskDS = (GDALDataset *) rstMask;
749 #else
750 maskDS = srcDS->FromHandle(rstMask);
751 #endif
752 }
753
754 GDALDataset *poDS = NULL;
755 GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName( "Memory" );
756 if( poDriver == NULL ) {
757 out.setError( "cannot create output driver");
758 return out;
759 }
760 poDS = poDriver->Create("", 0, 0, 0, GDT_Unknown, NULL );
761 if( poDS == NULL ) {
762 out.setError("Creation of output dataset failed" );
763 return out;
764 }
765 std::vector<std::string> nms = getNames();
766 std::string name = nms[0];
767
768 OGRSpatialReference *SRS = NULL;
769
770 OGRLayer *poLayer;
771 poLayer = poDS->CreateLayer(name.c_str(), SRS, wkbPolygon, NULL );
772 if( poLayer == NULL ) {
773 out.setError( "Layer creation failed" );
774 return out;
775 }
776 if (SRS != NULL) SRS->Release();
777
778 OGRFieldDefn oField(name.c_str(), trunc ? OFTInteger : OFTReal);
779 if( poLayer->CreateField( &oField ) != OGRERR_NONE ) {
780 out.setError( "Creating field failed");
781 return out;
782 }
783
784 GDALRasterBand *poBand;
785 poBand = srcDS->GetRasterBand(1);
786
787 //int hasNA=1;
788 //double naflag = poBand->GetNoDataValue(&hasNA);
789
790 CPLErr err;
791 if (usemask) {
792 GDALRasterBand *maskBand;
793 maskBand = maskDS->GetRasterBand(1);
794 if (trunc) {
795 err = GDALPolygonize(poBand, maskBand, poLayer, 0, NULL, NULL, NULL);
796 } else {
797 err = GDALFPolygonize(poBand, maskBand, poLayer, 0, NULL, NULL, NULL);
798 }
799 GDALClose(maskDS);
800 } else {
801 if (trunc) {
802 err = GDALPolygonize(poBand, NULL, poLayer, 0, NULL, NULL, NULL);
803 } else {
804 err = GDALFPolygonize(poBand, NULL, poLayer, 0, NULL, NULL, NULL);
805 }
806 }
807 if (err == 4) {
808 out.setError("polygonize error");
809 return out;
810 }
811 GDALClose(srcDS);
812
813 std::vector<double> fext;
814 SpatVector fvct;
815 out.read_ogr(poDS, "", "", fext, fvct);
816 GDALClose(poDS);
817
818 if (aggregate && (out.nrow() > 0)) {
819 out = out.aggregate(name, false);
820 }
821
822 if (!values) {
823 out.df = SpatDataFrame();
824 }
825
826 return out;
827 }
828
829
830
rgb2col(size_t r,size_t g,size_t b,SpatOptions & opt)831 SpatRaster SpatRaster::rgb2col(size_t r, size_t g, size_t b, SpatOptions &opt) {
832 SpatRaster out = geometry(1);
833 if (nlyr() < 3) {
834 out.setError("need at least three layers");
835 return out;
836 }
837 size_t mxlyr = std::max(std::max(r, g), b);
838 if (nlyr() < mxlyr) {
839 out.setError("layer number for R, G, B, cannot exceed nlyr()");
840 return out;
841 }
842
843 std::string filename = opt.get_filename();
844 opt.set_datatype("INT1U");
845 std::string driver;
846 if (filename == "") {
847 if (canProcessInMemory(opt)) {
848 driver = "MEM";
849 } else {
850 filename = tempFile(opt.get_tempdir(), ".tif");
851 opt.set_filenames({filename});
852 driver = "GTiff";
853 }
854 } else {
855 driver = opt.get_filetype();
856 getGDALdriver(filename, driver);
857 if (driver == "") {
858 setError("cannot guess file type from filename");
859 return out;
860 }
861 std::string errmsg;
862 if (!can_write(filename, opt.get_overwrite(), errmsg)) {
863 out.setError(errmsg);
864 return out;
865 }
866 }
867
868 std::vector<unsigned> lyrs = {(unsigned)r, (unsigned)g, (unsigned)b};
869 SpatOptions ops(opt);
870 *this = subset(lyrs, ops);
871 *this = collapse_sources();
872 GDALDatasetH hSrcDS, hDstDS;
873 if (!open_gdal(hSrcDS, 0, false, ops)) {
874 out.setError("cannot create dataset from source");
875 return out;
876 }
877 GDALRasterBandH R = GDALGetRasterBand(hSrcDS,1);
878 GDALRasterBandH G = GDALGetRasterBand(hSrcDS,1);
879 GDALRasterBandH B = GDALGetRasterBand(hSrcDS,1);
880
881 GDALColorTableH hColorTable= GDALCreateColorTable(GPI_RGB);
882
883 if (GDALComputeMedianCutPCT(R, G, B, NULL, 256, hColorTable, NULL, NULL) != CE_None) {
884 out.setError("cannot create color table");
885 GDALClose(hSrcDS);
886 return out;
887 }
888
889 if (!out.create_gdalDS(hDstDS, filename, driver, true, 0, {false}, {0.0}, {1.0}, opt)) {
890 out.setError("cannot create new dataset");
891 GDALClose(hSrcDS);
892 return out;
893 }
894
895 GDALRasterBandH hTarget = GDALGetRasterBand(hDstDS, 1);
896 GDALSetRasterColorInterpretation(hTarget, GCI_PaletteIndex);
897 if (GDALDitherRGB2PCT(R, G, B, hTarget, hColorTable, NULL, NULL) != CE_None) {
898 out.setError("cannot set color table");
899 GDALClose(hSrcDS);
900 GDALClose(hDstDS);
901 return out;
902 }
903 GDALClose(hSrcDS);
904
905 if (driver == "MEM") {
906 if (!out.from_gdalMEM(hDstDS, false, true)) {
907 out.setError("conversion failed (mem)");
908 GDALClose(hDstDS);
909 return out;
910 }
911 SpatDataFrame cdf;
912 cdf.add_column(1, "red");
913 cdf.add_column(1, "green");
914 cdf.add_column(1, "blue");
915 cdf.add_column(1, "alpha");
916 size_t nc = GDALGetColorEntryCount(hColorTable);
917 cdf.reserve(nc);
918
919 for (size_t i=0; i<nc; i++) {
920 const GDALColorEntry* col = GDALGetColorEntry(hColorTable, i);
921 cdf.iv[0].push_back(col->c1);
922 cdf.iv[1].push_back(col->c2);
923 cdf.iv[2].push_back(col->c3);
924 cdf.iv[3].push_back(col->c4);
925 }
926 out.source[0].hasColors.resize(1);
927 out.source[0].hasColors[0] = true;
928 out.source[0].cols.resize(1);
929 out.source[0].cols[0] = cdf;
930 } else {
931 if (GDALSetRasterColorTable(hTarget, hColorTable) != CE_None) {
932 out.setError("cannot set color table");
933 GDALClose(hDstDS);
934 return out;
935 }
936 }
937 GDALClose(hDstDS);
938 if (driver != "MEM") {
939 out = SpatRaster(filename, {-1}, {""}, {});
940 }
941 return out;
942 }
943
944
sievefilter(int threshold,int connections,SpatOptions & opt)945 SpatRaster SpatRaster::sievefilter(int threshold, int connections, SpatOptions &opt) {
946 SpatRaster out;
947 std::string filename = opt.get_filename();
948 std::string driver;
949 if (filename == "") {
950 if (canProcessInMemory(opt)) {
951 driver = "MEM";
952 } else {
953 filename = tempFile(opt.get_tempdir(), ".tif");
954 opt.set_filenames({filename});
955 driver = "GTiff";
956 }
957 } else {
958 driver = opt.get_filetype();
959 getGDALdriver(filename, driver);
960 if (driver == "") {
961 setError("cannot guess file type from filename");
962 return out;
963 }
964 std::string errmsg;
965 if (!can_write(filename, opt.get_overwrite(), errmsg)) {
966 out.setError(errmsg);
967 return out;
968 }
969 }
970
971 SpatOptions ops(opt);
972 GDALDatasetH hSrcDS, hDstDS;
973 if (!open_gdal(hSrcDS, 0, false, ops)) {
974 out.setError("cannot open input dataset");
975 return out;
976 }
977
978 GDALDriverH hDriver = GDALGetDriverByName( driver.c_str() );
979 if ( hDriver == NULL ) {
980 out.setError("empty driver");
981 return out;
982 }
983 if (!out.create_gdalDS(hDstDS, filename, driver, true, 0, source[0].has_scale_offset, source[0].scale, source[0].offset, opt)) {
984 out.setError("cannot create new dataset");
985 GDALClose(hSrcDS);
986 return out;
987 }
988
989 GDALRasterBandH hSrcBand = GDALGetRasterBand(hSrcDS, 1);
990 GDALRasterBandH hTargetBand = GDALGetRasterBand(hDstDS, 1);
991
992 if (!GDALSieveFilter(hSrcBand, NULL, hTargetBand, threshold, connections, NULL, NULL, NULL)) {
993 out.setError("sieve failed");
994 GDALClose(hSrcDS);
995 GDALClose(hDstDS);
996 return out;
997 }
998
999 GDALClose(hSrcDS);
1000 if (driver == "MEM") {
1001 if (!out.from_gdalMEM(hDstDS, false, true)) {
1002 out.setError("conversion failed (mem)");
1003 GDALClose(hDstDS);
1004 return out;
1005 }
1006 } else {
1007 out = SpatRaster(filename, {-1}, {""}, {});
1008 }
1009 GDALClose(hDstDS);
1010 return out;
1011 }
1012
1013
1014 /*
1015 SpatRaster SpatRaster::fillna(int threshold, int connections, SpatOptions &opt) {
1016 CPLErr GDALFillNodata(GDALRasterBandH hTargetBand, GDALRasterBandH hMaskBand, doubledfMaxSearchDist, intbDeprecatedOption, intnSmoothingIterations, char**papszOptions, GDALProgressFuncpfnProgress, void*pProgressArg)
1017 */
1018
1019
1020