data_raster.hpp
Go to the documentation of this file.
1
31#ifndef CCGL_DATA_RASTER_H
32#define CCGL_DATA_RASTER_H
33
34// include base headers
35#include <string>
36#include <map>
37#include <set>
38#include <fstream>
39#include <iomanip>
40#include <typeinfo>
41#include <cassert>
42// include openmp if supported
43#ifdef SUPPORT_OMP
44#include <omp.h>
45#endif /* SUPPORT_OMP */
46
47#include "basic.h"
48#include "utils_filesystem.h"
49#include "utils_string.h"
50#include "utils_array.h"
51#include "utils_math.h"
52
53#include "gdal_handler.h"
54#include "db_mongoc.h"
55
56using std::map;
57using std::set;
58#ifndef HAS_VARIADIC_TEMPLATES // Not support emplace(), e.g., MSVC2010
59using std::make_pair;
60#endif
61using std::setprecision;
62
63#ifndef CONST_CHARS
64#define CONST_CHARS static const char*
65#endif
66
67namespace ccgl {
68using namespace utils_string;
69using namespace utils_filesystem;
70using namespace utils_array;
71using namespace utils_math;
72#ifdef USE_MONGODB
73using namespace db_mongoc;
74#endif /* USE_MONGODB */
75
80namespace data_raster {
81CONST_CHARS HEADER_RS_NODATA = "NODATA_VALUE";
94CONST_CHARS HEADER_INC_NODATA = "INCLUDE_NODATA";
96CONST_CHARS STATS_RS_VALIDNUM = "VALID_CELLNUMBER";
104
107CONST_CHARS PARAM_ABSTRACTION_TYPE_PHYSICAL = "PHYSICAL";
108
109
110typedef std::pair<int, int> ROW_COL;
111typedef std::pair<double, double> XY_COOR;
112
116typedef enum {
129
135string RasterDataTypeToString(int type);
136
140RasterDataType StringToRasterDataType(const string& stype);
141
145RasterDataType TypeToRasterDataType(const std::type_info& t);
146
151
152#ifdef USE_GDAL
153GDALDataType CvtToGDALDataType(RasterDataType type);
154#endif
155
160
166void CopyHeader(const STRDBL_MAP& refers, STRDBL_MAP& dst);
167
171template <typename T>
172void UpdateHeader(STRDBL_MAP& header, const string& key, T val) {
173 if (header.find(key) == header.end()) {
174#ifdef HAS_VARIADIC_TEMPLATES
175 header.emplace(key, CVT_DBL(val));
176#else
177 header.insert(make_pair(key, CVT_DBL(val)));
178#endif
179 }
180 else { header.at(key) = CVT_DBL(val); }
181}
182
187
191void UpdateStrHeader(STRING_MAP& strheader, const string& key, const string& val);
196
200void InitialStatsMap(STRDBL_MAP& stats, map<string, double*>& stats2d);
201
209template <typename T>
210bool ReadAscFile(const string& filename, STRDBL_MAP& header, T*& values) {
211 StatusMessage(("Read " + filename + "...").c_str());
212 vector<string> content_strs;
213 if (!LoadPlainTextFile(filename, content_strs)) { return false; }
214 if (content_strs.size() < 7) {
215 StatusMessage("Error: ASCII raster data requires at least 7 lines!");
216 return false;
217 }
218 double cellsize = -1.;
219 double xll = MISSINGFLOAT;
220 double yll = MISSINGFLOAT;
221 double nodata = MISSINGFLOAT;
222 bool xllcorner = false;
223 bool yllcorner = false;
224 int rows = -1;
225 int cols = -1;
226 int vidx = 0;
227 for (auto iter = content_strs.begin(); iter != content_strs.end(); ++iter) {
228 vector<string> tokens = SplitString(*iter);
229 if (tokens.empty()) { continue; }
230 bool str2value = false;
231 if (StringMatch(tokens[0], HEADER_RS_NCOLS)) {
232 if (tokens.size() < 2) { return false; }
233 cols = CVT_INT(IsInt(tokens[1], str2value));
234 if (!str2value) { return false; }
235 }
236 else if (StringMatch(tokens[0], HEADER_RS_NROWS)) {
237 if (tokens.size() < 2) { return false; }
238 rows = CVT_INT(IsInt(tokens[1], str2value));
239 if (!str2value) { return false; }
240 }
241 else if (StringMatch(tokens[0], HEADER_RS_XLL) ||
242 StringMatch(tokens[0], HEADER_RS_XLLCOR)) {
243 if (tokens.size() < 2) { return false; }
244 xll = IsDouble(tokens[1], str2value);
245 if (StringMatch(tokens[0], HEADER_RS_XLLCOR)) { xllcorner = true; }
246 if (!str2value) { return false; }
247 }
248 else if (StringMatch(tokens[0], HEADER_RS_YLL) ||
249 StringMatch(tokens[0], HEADER_RS_YLLCOR)) {
250 if (tokens.size() < 2) { return false; }
251 yll = IsDouble(tokens[1], str2value);
252 if (StringMatch(tokens[0], HEADER_RS_YLLCOR)) { yllcorner = true; }
253 if (!str2value) { return false; }
254 }
255 else if (StringMatch(tokens[0], HEADER_RS_CELLSIZE)) {
256 if (tokens.size() < 2) { return false; }
257 cellsize = IsDouble(tokens[1], str2value);
258 if (!str2value) { return false; }
259 }
260 else if (StringMatch(tokens[0], HEADER_RS_NODATA)) {
261 if (tokens.size() < 2) { return false; }
262 nodata = IsDouble(tokens[1], str2value);
263 if (!str2value) { return false; }
264 }
265 else {
266 if (rows < 0 || cols < 0) {
267 StatusMessage("Warning: NCOLS and NROWS should be defined first!");
268 }
269 if (nullptr == values) { Initialize1DArray(rows * cols, values, nodata); }
270 for (auto it_value = tokens.begin(); it_value != tokens.end(); ++it_value) {
271 double tmpv = IsDouble(*it_value, str2value);
272 if (!str2value) {
273 StatusMessage("Error: No value occurred in Raster data matrix!");
274 if (nullptr != values) { Release1DArray(values); }
275 return false;
276 }
277 if (vidx == rows * cols) {
278 StatusMessage("Error: Count of values MUST equal to rows * cols!");
279 if (nullptr != values) { Release1DArray(values); }
280 return false;
281 }
282 values[vidx++] = static_cast<T>(tmpv);
283 }
284 }
285 }
286 if (rows < 0 || cols < 0 || cellsize < 0 ||
288 StatusMessage("Error: Header information incomplete!");
289 return false;
290 }
291 // default is center, if corner, then:
292 if (xllcorner) { xll += 0.5 * cellsize; }
293 if (yllcorner) { yll += 0.5 * cellsize; }
294 UpdateHeader(header, HEADER_RS_NCOLS, cols);
295 UpdateHeader(header, HEADER_RS_NROWS, rows);
296 UpdateHeader(header, HEADER_RS_XLL, xll);
297 UpdateHeader(header, HEADER_RS_YLL, yll);
298 UpdateHeader(header, HEADER_RS_CELLSIZE, cellsize);
299 UpdateHeader(header, HEADER_RS_NODATA, nodata);
300 UpdateHeader(header, HEADER_RS_LAYERS, 1);
301 UpdateHeader(header, HEADER_RS_CELLSNUM, cols * rows);
302
303 return true;
304}
305
311bool WriteAscHeaders(const string& filename, const STRDBL_MAP& header);
312
319template <typename T>
320bool WriteSingleAsc(const string& filename, const STRDBL_MAP& header, T* values) {
321 string abs_filename = GetAbsolutePath(filename);
322 string dirname = GetPathFromFullName(abs_filename);
323 if (!PathExists(dirname)) { MakeDirectory(dirname); }
324 if (!WriteAscHeaders(abs_filename, header)) { return false; }
325 std::ofstream raster_file(abs_filename.c_str(), std::ios::app | std::ios::out);
326 if (!raster_file.is_open()) {
327 StatusMessage("Error opening file: " + abs_filename);
328 return false;
329 }
330 int rows = CVT_INT(header.at(HEADER_RS_NROWS));
331 int cols = CVT_INT(header.at(HEADER_RS_NCOLS));
332 for (int i = 0; i < rows; i++) {
333 for (int j = 0; j < cols; j++) {
334 raster_file << setprecision(6) << values[i * cols + j] << " ";
335 }
336 raster_file << endl;
337 }
338 raster_file << endl;
339 raster_file.close();
340 return true;
341}
342
343#ifdef USE_GDAL
353template <typename T>
354bool ReadRasterFileByGdal(const string& filename, STRDBL_MAP& header, T*& values,
355 RasterDataType& in_type, string& srs) {
356 StatusMessage(("Read " + filename + "...").c_str());
357 GDALRasterDSHandle po_dataset(OpenRaster(filename.c_str()));
358 // GDALDataset* po_dataset = static_cast<GDALDataset*>(GDALOpen(filename.c_str(),
359 // GA_ReadOnly));
360 if (nullptr == po_dataset) {
361 StatusMessage("Open file " + filename + " failed.");
362 return false;
363 }
364 GDALRasterBand* po_band = po_dataset->GetRasterBand(1);
365 int n_rows = po_band->GetYSize();
366 int n_cols = po_band->GetXSize();
367 int get_value_flag = false;
368 double nodata = po_band->GetNoDataValue(&get_value_flag);
369 int approx_minmax = false;
370 double minmax[2];
371 T* tmprasterdata = nullptr;
372 bool read_as_signedbyte = false;
373 unsigned char* uchar_data = nullptr;
374 signed char* char_data = nullptr; // DO NOT use char*
375 vuint16_t* uint16_data = nullptr; // 16-bit unsigned integer
376 vint16_t* int16_data = nullptr; // 16-bit signed integer
377 vuint32_t* uint32_data = nullptr; // 32-bit unsigned integer
378 vint32_t* int32_data = nullptr; // 32-bit signed integer
379 vuint64_t* uint64_data = nullptr; // 64-bit unsigned integer
380 vint64_t* int64_data = nullptr; // 64-bit signed integer
381 float* float_data = nullptr;
382 double* double_data = nullptr;
383 CPLErr result;
384 switch (po_band->GetRasterDataType()) {
385 case GDT_Byte:
386 // In GDAL, GDT_Byte is an 8-bit unsigned integer (unsigned char), ranges from 0 to 255.
387 // While in ArcGIS, both 8-bit signed char (ranges from -128 to 127) and unsigned char
388 // are supported. Both types will be recognized as GDT_Byte by GDAL.
389 //
390 // So,
391 // 1) maximum <= 127 and minimum < 0 ==> signed char
392 // 2) maximum <= 127 and minimum >= 0 and no_data_value_ < 0 ==> signed char
393 // Otherwise, unsigned char.
394 //
395 // Update (08/09/2023): GDAL>=3.7 added the support of GDT_Int8. Keep this code for compatibility!
396 //
397 po_band->ComputeRasterMinMax(approx_minmax, minmax);
398 if ((minmax[1] <= 127 && minmax[0] < 0)
399 || (minmax[1] <= 127 && minmax[0] >= 0 && (!get_value_flag || (get_value_flag && nodata < 0)))) {
400 read_as_signedbyte = true;
401 }
402 uchar_data = static_cast<unsigned char*>(CPLMalloc(sizeof(unsigned char) * n_cols * n_rows));
403 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, uchar_data,
404 n_cols, n_rows, GDT_Byte, 0, 0);
405 if (result != CE_None) {
406 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
407 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
408 return false;
409 }
410 if (read_as_signedbyte) {
411 StatusMessage("Read GDT_Byte raster as signed char!");
412 Initialize1DArray(n_rows * n_cols, char_data, uchar_data);
413 Initialize1DArray(n_rows * n_cols, tmprasterdata, char_data);
414 Release1DArray(char_data);
415 in_type = RDT_Int8;
416 } else {
417 StatusMessage("Read GDT_Byte raster as unsigned char!");
418 Initialize1DArray(n_rows * n_cols, tmprasterdata, uchar_data);
419 in_type = RDT_UInt8;
420 }
421 CPLFree(uchar_data);
422 break;
423#if GDAL_VERSION_MAJOR >= 3 && GDAL_VERSION_MINOR >= 7
424 case GDT_Int8:
425 char_data = static_cast<signed char*>(CPLMalloc(sizeof(signed char) * n_cols * n_rows));
426 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, char_data,
427 n_cols, n_rows, GDT_Int8, 0, 0);
428 if (result != CE_None) {
429 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
430 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
431 return false;
432 }
433 Initialize1DArray(n_rows * n_cols, tmprasterdata, char_data);
434 CPLFree(char_data);
435 in_type = RDT_Int8;
436 break;
437#endif
438 case GDT_UInt16:
439 uint16_data = static_cast<vuint16_t*>(CPLMalloc(sizeof(vuint16_t) * n_cols * n_rows));
440 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, uint16_data,
441 n_cols, n_rows, GDT_UInt16, 0, 0);
442 if (result != CE_None) {
443 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
444 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
445 return false;
446 }
447 Initialize1DArray(n_rows * n_cols, tmprasterdata, uint16_data);
448 CPLFree(uint16_data);
449 in_type = RDT_UInt16;
450 break;
451 case GDT_Int16:
452 int16_data = static_cast<vint16_t*>(CPLMalloc(sizeof(vint16_t) * n_cols * n_rows));
453 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, int16_data,
454 n_cols, n_rows, GDT_Int16, 0, 0);
455 if (result != CE_None) {
456 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
457 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
458 return false;
459 }
460 Initialize1DArray(n_rows * n_cols, tmprasterdata, int16_data);
461 CPLFree(int16_data);
462 in_type = RDT_Int16;
463 break;
464 case GDT_UInt32:
465 uint32_data = static_cast<vuint32_t*>(CPLMalloc(sizeof(vuint32_t) * n_cols * n_rows));
466 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, uint32_data,
467 n_cols, n_rows, GDT_UInt32, 0, 0);
468 if (result != CE_None) {
469 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
470 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
471 return false;
472 }
473 Initialize1DArray(n_rows * n_cols, tmprasterdata, uint32_data);
474 CPLFree(uint32_data);
475 in_type = RDT_UInt32;
476 break;
477 case GDT_Int32:
478 int32_data = static_cast<vint32_t*>(CPLMalloc(sizeof(vint32_t) * n_cols * n_rows));
479 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, int32_data,
480 n_cols, n_rows, GDT_Int32, 0, 0);
481 if (result != CE_None) {
482 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
483 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
484 return false;
485 }
486 Initialize1DArray(n_rows * n_cols, tmprasterdata, int32_data);
487 CPLFree(int32_data);
488 in_type = RDT_Int32;
489 break;
490#if GDAL_VERSION_MAJOR >= 3 && GDAL_VERSION_MINOR >= 5
491 case GDT_UInt64:
492 uint64_data = static_cast<vuint64_t*>(CPLMalloc(sizeof(vuint64_t) * n_cols * n_rows));
493 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, uint64_data,
494 n_cols, n_rows, GDT_UInt64, 0, 0);
495 if (result != CE_None) {
496 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
497 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
498 return false;
499 }
500 Initialize1DArray(n_rows * n_cols, tmprasterdata, uint64_data);
501 CPLFree(uint64_data);
502 in_type = RDT_UInt64;
503 break;
504 case GDT_Int64:
505 int64_data = static_cast<vint64_t*>(CPLMalloc(sizeof(vint64_t) * n_cols * n_rows));
506 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, int64_data,
507 n_cols, n_rows, GDT_Int64, 0, 0);
508 if (result != CE_None) {
509 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
510 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
511 return false;
512 }
513 Initialize1DArray(n_rows * n_cols, tmprasterdata, int64_data);
514 CPLFree(int64_data);
515 in_type = RDT_Int64;
516 break;
517#endif
518 case GDT_Float32:
519 float_data = static_cast<float*>(CPLMalloc(sizeof(float) * n_cols * n_rows));
520 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, float_data,
521 n_cols, n_rows, GDT_Float32, 0, 0);
522 if (result != CE_None) {
523 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
524 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
525 return false;
526 }
527 Initialize1DArray(n_rows * n_cols, tmprasterdata, float_data);
528 CPLFree(float_data);
529 in_type = RDT_Float;
530 break;
531 case GDT_Float64:
532 double_data = static_cast<double*>(CPLMalloc(sizeof(double) * n_cols * n_rows));
533 result = po_band->RasterIO(GF_Read, 0, 0, n_cols, n_rows, double_data,
534 n_cols, n_rows, GDT_Float64, 0, 0);
535 if (result != CE_None) {
536 StatusMessage("RaterIO trouble: " + string(CPLGetLastErrorMsg()));
537 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
538 return false;
539 }
540 Initialize1DArray(n_rows * n_cols, tmprasterdata, double_data);
541 CPLFree(double_data);
542 in_type = RDT_Double;
543 break;
544 default:
545 StatusMessage("Unexpected GDALDataType: " + string(RasterDataTypeToString(in_type)));
546 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
547 return false;
548 }
549
550 if (!get_value_flag) { // NoData value is not defined!
551 nodata = DefaultNoDataByType(in_type);
552 }
553
554 double geo_trans[6];
555 po_dataset->GetGeoTransform(geo_trans);
556 UpdateHeader(header, HEADER_RS_NCOLS, n_cols);
557 UpdateHeader(header, HEADER_RS_NROWS, n_rows);
558 UpdateHeader(header, HEADER_RS_NODATA, nodata);
559 UpdateHeader(header, HEADER_RS_CELLSIZE, geo_trans[1]);
560 UpdateHeader(header, HEADER_RS_XLL, geo_trans[0] + 0.5 * geo_trans[1]);
561 UpdateHeader(header, HEADER_RS_YLL, geo_trans[3] + (n_rows - 0.5) * geo_trans[5]);
562 UpdateHeader(header, HEADER_RS_LAYERS, 1.);
563 UpdateHeader(header, HEADER_RS_CELLSNUM, n_cols * n_rows);
564 srs = string(po_dataset->GetProjectionRef());
565
566 // GDALClose(po_dataset); // When use GDALRasterDSHandle, No need to explicitly close dataset
567
568 values = tmprasterdata;
569 return true;
570}
571
580template <typename T>
581bool WriteSingleGeotiff(const string& filename, const STRDBL_MAP& header,
582 const STRING_MAP& opts, T* values) {
583 int n_rows = CVT_INT(header.at(HEADER_RS_NROWS));
584 int n_cols = CVT_INT(header.at(HEADER_RS_NCOLS));
585 char** papsz_options = nullptr;
586 void* new_values = nullptr;
587 bool recreate_flag = true; // recreate data array because inconsistent of datatypes
588 double old_nodata = header.at(HEADER_RS_NODATA);
589 bool change_nodata = false;
590 double new_nodata = old_nodata; // change nodata when convert signed datatype to unsigned
591 bool convert_permit = true; // DO NOT allow negative cell values be converted to unsigned datatype
592 RasterDataType tmptype = TypeToRasterDataType(typeid(T));
594 if (outtype == RDT_Unknown) {
595 outtype = tmptype;
596 recreate_flag = false;
597 }
598 else if (outtype == tmptype && outtype != RDT_Int8) { // RDT_Int8 need special handling
599 recreate_flag = false;
600 }
601 else {
602 if (outtype == RDT_UInt8) { // [0, 255] --> GDT_Byte in GDAL, use unsigned char
603 new_values = static_cast<unsigned char*>(CPLMalloc(sizeof(unsigned char) * n_cols * n_rows));
604 unsigned char* values_uchar = static_cast<unsigned char*>(new_values);
605 if (old_nodata < 0 || old_nodata > UINT8_MAX) {
606 new_nodata = UINT8_MAX;
607 change_nodata = true;
608 }
609 int illegal_count = 0;
610#pragma omp parallel for reduction(+:illegal_count)
611 for (int i = 0; i < n_cols * n_rows; i++) {
612 if (FloatEqual(values[i], old_nodata) && change_nodata) {
613 values_uchar[i] = UINT8_MAX;
614 continue;
615 }
616 if (values[i] < 0 || values[i] > UINT8_MAX) illegal_count += 1;
617 values_uchar[i] = static_cast<unsigned char>(values[i]);
618 }
619 if (illegal_count > 0) convert_permit = false;
620 }
621 else if (outtype == RDT_Int8) { // [-128, 127]
622 // https://gdal.org/drivers/raster/gtiff.html
623#if GDAL_VERSION_MAJOR < 3 || (GDAL_VERSION_MAJOR == 3 && GDAL_VERSION_MINOR < 7)
624 papsz_options = CSLSetNameValue(papsz_options, "PIXELTYPE", "SIGNEDBYTE");
625#endif
626 new_values = static_cast<signed char*>(CPLMalloc(sizeof(signed char) * n_cols * n_rows));
627 signed char* values_char = static_cast<signed char*>(new_values);
628 if (old_nodata < INT8_MIN || old_nodata > INT8_MAX) {
629 new_nodata = INT8_MIN;
630 change_nodata = true;
631 }
632 int illegal_count = 0;
633#pragma omp parallel for reduction(+:illegal_count)
634 for (int i = 0; i < n_cols * n_rows; i++) {
635 if (FloatEqual(values[i], old_nodata) && change_nodata) {
636 values_char[i] = INT8_MIN;
637 continue;
638 }
639 if (values[i] < INT8_MIN || values[i] > INT8_MAX) { illegal_count += 1; }
640 values_char[i] = static_cast<signed char>(values[i]);
641 }
642 if (illegal_count > 0) convert_permit = false;
643 }
644 else if (outtype == RDT_UInt16) { // [0, 65535]
645 new_values = static_cast<vuint16_t*>(CPLMalloc(sizeof(vuint16_t) * n_cols * n_rows));
646 vuint16_t* values_uint16 = static_cast<vuint16_t*>(new_values);
647 if (old_nodata < 0 || old_nodata > UINT16_MAX) {
648 new_nodata = UINT16_MAX;
649 change_nodata = true;
650 }
651 int illegal_count = 0;
652#pragma omp parallel for reduction(+:illegal_count)
653 for (int i = 0; i < n_cols * n_rows; i++) {
654 if (FloatEqual(values[i], old_nodata) && change_nodata) {
655 values_uint16[i] = UINT16_MAX;
656 continue;
657 }
658 if (values[i] < 0 || values[i] > UINT16_MAX) illegal_count += 1;
659 values_uint16[i] = static_cast<vuint16_t>(values[i]);
660 }
661 if (illegal_count > 0) convert_permit = false;
662 }
663 else if (outtype == RDT_Int16) { // [-32768, 32767]
664 new_values = static_cast<vint16_t*>(CPLMalloc(sizeof(vint16_t) * n_cols * n_rows));
665 vint16_t* values_int16 = static_cast<vint16_t*>(new_values);
666 if (old_nodata < INT16_MIN || old_nodata > INT16_MAX) {
667 new_nodata = INT16_MIN;
668 change_nodata = true;
669 }
670 int illegal_count = 0;
671#pragma omp parallel for reduction(+:illegal_count)
672 for (int i = 0; i < n_cols * n_rows; i++) {
673 if (FloatEqual(values[i], old_nodata) && change_nodata) {
674 values_int16[i] = INT16_MIN;
675 continue;
676 }
677 if (values[i] < INT16_MIN || values[i] > INT16_MAX) illegal_count += 1;
678 values_int16[i] = static_cast<vint16_t>(values[i]);
679 }
680 if (illegal_count > 0) convert_permit = false;
681 }
682 else if (outtype == RDT_UInt32) { // [0, 4294967295]
683 new_values = static_cast<vuint32_t*>(CPLMalloc(sizeof(vuint32_t) * n_cols * n_rows));
684 vuint32_t* values_uint32 = static_cast<vuint32_t*>(new_values);
685 if (old_nodata < 0 || old_nodata > UINT32_MAX) {
686 new_nodata = UINT32_MAX;
687 change_nodata = true;
688 }
689 int illegal_count = 0;
690#pragma omp parallel for reduction(+:illegal_count)
691 for (int i = 0; i < n_cols * n_rows; i++) {
692 if (FloatEqual(values[i], old_nodata) && change_nodata) {
693 values_uint32[i] = UINT32_MAX;
694 continue;
695 }
696 if (values[i] < 0 || values[i] > UINT32_MAX) illegal_count += 1;
697 values_uint32[i] = static_cast<vuint32_t>(values[i]);
698 }
699 if (illegal_count > 0) convert_permit = false;
700 }
701 else if (outtype == RDT_Int32) { // [-2147483648, 2147483647]
702 new_values = static_cast<vint32_t*>(CPLMalloc(sizeof(vint32_t) * n_cols * n_rows));
703 vint32_t* values_int32 = static_cast<vint32_t*>(new_values);
704 if (old_nodata < INT32_MIN || old_nodata > INT32_MAX) {
705 new_nodata = INT32_MIN;
706 change_nodata = true;
707 }
708 int illegal_count = 0;
709#pragma omp parallel for reduction(+:illegal_count)
710 for (int i = 0; i < n_cols * n_rows; i++) {
711 if (FloatEqual(values[i], old_nodata) && change_nodata) {
712 values_int32[i] = INT32_MIN;
713 continue;
714 }
715 if (values[i] < INT32_MIN || values[i] > INT32_MAX) illegal_count += 1;
716 values_int32[i] = static_cast<vint32_t>(values[i]);
717 }
718 if (illegal_count > 0) convert_permit = false;
719 }
720#if GDAL_VERSION_MAJOR >= 3 && GDAL_VERSION_MINOR >=5
721 else if (outtype == RDT_UInt64) { // [0, 18446744073709551615]
722 new_values = static_cast<vuint64_t*>(CPLMalloc(sizeof(vuint64_t) * n_cols * n_rows));
723 vuint64_t* values_uint64 = static_cast<vuint64_t*>(new_values);
724 if (old_nodata < 0 || old_nodata > UINT64_MAX) {
725 new_nodata = UINT64_MAX;
726 change_nodata = true;
727 }
728 int illegal_count = 0;
729#pragma omp parallel for reduction(+:illegal_count)
730 for (int i = 0; i < n_cols * n_rows; i++) {
731 if (FloatEqual(values[i], old_nodata) && change_nodata) {
732 values_uint64[i] = UINT64_MAX;
733 continue;
734 }
735 if (values[i] < 0 || values[i] > UINT64_MAX) illegal_count += 1;
736 values_uint64[i] = static_cast<vuint64_t>(values[i]);
737 }
738 if (illegal_count > 0) convert_permit = false;
739 }
740 else if (outtype == RDT_Int64) { // [-18446744073709551615, 18446744073709551615]
741 new_values = static_cast<vint64_t*>(CPLMalloc(sizeof(vint64_t) * n_cols * n_rows));
742 vint64_t* values_int64 = static_cast<vint64_t*>(new_values);
743 if (old_nodata < INT64_MIN || old_nodata > INT64_MAX) {
744 new_nodata = INT64_MIN;
745 change_nodata = true;
746 }
747 int illegal_count = 0;
748#pragma omp parallel for reduction(+:illegal_count)
749 for (int i = 0; i < n_cols * n_rows; i++) {
750 if (FloatEqual(values[i], old_nodata) && change_nodata) {
751 values_int64[i] = INT64_MIN;
752 continue;
753 }
754 if (values[i] < INT64_MIN || values[i] > INT64_MAX) illegal_count += 1;
755 values_int64[i] = static_cast<vint64_t>(values[i]);
756 }
757 if (illegal_count > 0) convert_permit = false;
758 }
759#endif
760 else if (outtype == RDT_Float) {
761 new_values = static_cast<float*>(CPLMalloc(sizeof(float) * n_cols * n_rows));
762 float* values_float = static_cast<float*>(new_values);
763#pragma omp parallel for
764 for (int i = 0; i < n_cols * n_rows; i++) values_float[i] = static_cast<float>(values[i]);
765 }
766 else if (outtype == RDT_Double) {
767 new_values = static_cast<double*>(CPLMalloc(sizeof(double) * n_cols * n_rows));
768 double* values_double = static_cast<double*>(new_values);
769#pragma omp parallel for
770 for (int i = 0; i < n_cols * n_rows; i++) values_double[i] = static_cast<double>(values[i]);
771 }
772 }
773 if ((outtype == RDT_Unknown || (recreate_flag && nullptr == new_values))
774 || (!convert_permit)) {
775 cout << "Error: The specific raster output data type is not allowed!\n";
776 if (papsz_options != nullptr) { CSLDestroy(papsz_options); }
777 if (new_values != nullptr) { CPLFree(new_values); }
778 return false;
779 }
780 string dirname = GetPathFromFullName(filename);
781 if (!PathExists(dirname)) { MakeDirectory(dirname); }
782 GDALRasterDSHandle po_dst_ds(CreateRaster("GTiff", filename.c_str(),
783 n_cols, n_rows, 1,
784 CvtToGDALDataType(outtype), papsz_options));
785 // GDALDriver* po_driver = GetGDALDriverManager()->GetDriverByName("GTiff");
786 // if (nullptr == po_driver) { return false; }
787 // GDALDataset* po_dst_ds = po_driver->Create(filename.c_str(), n_cols, n_rows, 1,
788 // CvtToGDALDataType(outtype), papsz_options);
789 if (nullptr == po_dst_ds) { return false; }
790 GDALRasterBand* po_dst_band = po_dst_ds->GetRasterBand(1);
791 CPLErr result = CE_None;
792 if (!recreate_flag) {
793 result = po_dst_band->RasterIO(GF_Write, 0, 0, n_cols, n_rows, values,
794 n_cols, n_rows, CvtToGDALDataType(outtype), 0, 0);
795 }
796 else {
797 result = po_dst_band->RasterIO(GF_Write, 0, 0, n_cols, n_rows, new_values,
798 n_cols, n_rows, CvtToGDALDataType(outtype), 0, 0);
799 }
800 if (result != CE_None || nullptr == po_dst_band) {
801 StatusMessage("RaterIO Error: " + string(CPLGetLastErrorMsg()));
802 if (papsz_options != nullptr) { CSLDestroy(papsz_options); }
803 if (new_values != nullptr) { CPLFree(new_values); }
804 // GDALClose(po_dst_ds); // When use GDALRasterDSHandle, No need to explicitly close dataset
805 return false;
806 }
807 po_dst_band->SetNoDataValue(new_nodata);
808
809 double geo_trans[6]; // Write header information
810 geo_trans[0] = header.at(HEADER_RS_XLL) - 0.5 * header.at(HEADER_RS_CELLSIZE);
811 geo_trans[1] = header.at(HEADER_RS_CELLSIZE);
812 geo_trans[2] = 0.;
813 geo_trans[3] = header.at(HEADER_RS_YLL) + (n_rows - 0.5) * header.at(HEADER_RS_CELLSIZE);
814 geo_trans[4] = 0.;
815 geo_trans[5] = -header.at(HEADER_RS_CELLSIZE);
816 po_dst_ds->SetGeoTransform(geo_trans);
817 if (opts.find(HEADER_RS_SRS) == opts.end()) {
818 po_dst_ds->SetProjection("");
819 } else {
820 po_dst_ds->SetProjection(opts.at(HEADER_RS_SRS).c_str());
821 }
822 if (papsz_options != nullptr) { CSLDestroy(papsz_options); }
823 if (new_values != nullptr) { CPLFree(new_values); }
824 // GDALClose(po_dst_ds); // When use GDALRasterDSHandle, No need to explicitly close dataset
825
826 return true;
827}
828#endif /* USE_GDAL */
829
837template <typename T>
838bool WriteRasterToFile(const string& filename, const STRDBL_MAP& header,
839 const STRING_MAP& opts, T* values) {
840 if (GetPathFromFullName(filename).empty()) { return false; }
841 string abs_filename = GetAbsolutePath(filename);
842 string filetype = GetUpper(GetSuffix(abs_filename));
843 if (StringMatch(filetype, ASCIIExtension)) {
844 return WriteSingleAsc(abs_filename, header, values);
845 }
846#ifdef USE_GDAL
847 if (StringMatch(filetype, GTiffExtension)) {
848 return WriteSingleGeotiff(abs_filename, header, opts, values);
849 }
850 return WriteSingleGeotiff(ReplaceSuffix(abs_filename, GTiffExtension),
851 header, opts, values);
852#else
853 StatusMessage("Warning: Without GDAL, ASC file will be exported as default!");
854 return WriteSingleAsc(ReplaceSuffix(abs_filename, ASCIIExtension), header, values);
855#endif /* USE_GDAL */
856}
857
858#ifdef USE_MONGODB
868template <typename T>
869bool ReadGridFsFile(MongoGridFs* gfs, const string& filename,
870 T*& data, STRDBL_MAP& header,
871 STRING_MAP& header_str,
872 const STRING_MAP& opts /* = STRING_MAP() */) {
873 // Get stream data and metadata by file name
874 vint length;
875 char* buf = nullptr;
876 if (!gfs->GetStreamData(filename, buf, length, nullptr, &opts) ||
877 nullptr == buf) {
878 return false;
879 }
880 bson_t* bmeta = gfs->GetFileMetadata(filename, nullptr, opts);
881
882 // Retrieve raster header values
883 bson_iter_t iter; // Loop the metadata, add to `header_str` or `header`
884 if (nullptr != bmeta && bson_iter_init(&iter, bmeta)) {
885 while (bson_iter_next(&iter)) {
886 const char* key = bson_iter_key(&iter);
887 if (header.find(key) != header.end()) {
888 GetNumericFromBsonIterator(&iter, header[key]);
889 }
890 else {
891 header_str[key] = GetStringFromBsonIterator(&iter);
892 }
893 }
894 }
895 bson_destroy(bmeta); // Destroy bson of metadata immediately after use
896
897 int n_rows = CVT_INT(header.at(HEADER_RS_NROWS));
898 int n_cols = CVT_INT(header.at(HEADER_RS_NCOLS));
899 int n_lyrs = CVT_INT(header.at(HEADER_RS_LAYERS));
900 int n_cells = CVT_INT(header.at(HEADER_RS_CELLSNUM));
901 if (n_rows < 0 || n_cols < 0 || n_lyrs < 0) { // missing essential metadata
902 delete[] buf;
903 return false;
904 }
905 int value_count = n_cells * n_lyrs;
906 size_t size_dtype = length / value_count;
907
909 if (header_str.find(HEADER_RSOUT_DATATYPE) != header_str.end()) {
910 rstype = StringToRasterDataType(header_str.at(HEADER_RSOUT_DATATYPE));
911 }
912 if (rstype == RDT_Unknown) {
913 StatusMessage("Unknown data type in MongoDB GridFS!");
914 delete[] buf;
915 return false;
916 }
917
918 if (rstype == RDT_Double && size_dtype == sizeof(double)) {
919 double* data_dbl = reinterpret_cast<double*>(buf);
920 Initialize1DArray(value_count, data, data_dbl);
921 //Release1DArray(data_dbl);
922 }
923 else if (rstype == RDT_Float && size_dtype == sizeof(float)) {
924 float* data_flt = reinterpret_cast<float*>(buf);
925 Initialize1DArray(value_count, data, data_flt);
926 //Release1DArray(data_flt);
927 }
928 else if (rstype == RDT_Int32 && size_dtype == sizeof(vint32_t)) {
929 vint32_t* data_int32 = reinterpret_cast<vint32_t*>(buf);
930 Initialize1DArray(value_count, data, data_int32);
931 //Release1DArray(data_int32);
932 }
933 else if (rstype == RDT_UInt32 && size_dtype == sizeof(vuint32_t)) {
934 vuint32_t* data_uint32 = reinterpret_cast<vuint32_t*>(buf);
935 Initialize1DArray(value_count, data, data_uint32);
936 //Release1DArray(data_uint32);
937 }
938 else if (rstype == RDT_Int64 && size_dtype == sizeof(vint64_t)) {
939 vint64_t* data_int64 = reinterpret_cast<vint64_t*>(buf);
940 Initialize1DArray(value_count, data, data_int64);
941 //Release1DArray(data_int64);
942 }
943 else if (rstype == RDT_UInt64 && size_dtype == sizeof(vuint64_t)) {
944 vuint64_t* data_uint64 = reinterpret_cast<vuint64_t*>(buf);
945 Initialize1DArray(value_count, data, data_uint64);
946 //Release1DArray(data_uint64);
947 }
948 else if (rstype == RDT_Int16 && size_dtype == sizeof(vint16_t)) {
949 vint16_t* data_int16 = reinterpret_cast<vint16_t*>(buf);
950 Initialize1DArray(value_count, data, data_int16);
951 //Release1DArray(data_int16);
952 }
953 else if (rstype == RDT_UInt16 && size_dtype == sizeof(vuint16_t)) {
954 vuint16_t* data_uint16 = reinterpret_cast<vuint16_t*>(buf);
955 Initialize1DArray(value_count, data, data_uint16);
956 //Release1DArray(data_uint16);
957 }
958 else if (rstype == RDT_Int8 && size_dtype == sizeof(vint8_t)) {
959 vint8_t* data_int8 = reinterpret_cast<vint8_t*>(buf);
960 Initialize1DArray(value_count, data, data_int8);
961 //Release1DArray(data_int8);
962 }
963 else if (rstype == RDT_UInt8 && size_dtype == sizeof(vuint8_t)) {
964 vuint8_t* data_uint8 = reinterpret_cast<vuint8_t*>(buf);
965 Initialize1DArray(value_count, data, data_uint8);
966 //Release1DArray(data_uint8);
967 }
968 else {
969 StatusMessage("Unconsistent of data type and size!");
970 delete[] buf;
971 return false;
972 }
973 delete[] buf;
974 return true;
975}
976
987template <typename T>
988bool WriteStreamDataAsGridfs(MongoGridFs* gfs, const string& filename,
989 STRDBL_MAP& header, T* values, const int datalength,
990 const STRING_MAP& opts = STRING_MAP()) {
991 STRING_MAP curopts;
992 CopyStringMap(opts, curopts);
993 RasterDataType temp_type = TypeToRasterDataType(typeid(T));
994 if (curopts.find(HEADER_RSOUT_DATATYPE) == curopts.end()) {
996 }
997 gfs->RemoveFile(filename, nullptr, curopts);
998 bson_t p = BSON_INITIALIZER;
999 double intpart; // https://stackoverflow.com/a/1521682/4837280
1000 for (auto iter = header.begin(); iter != header.end(); ++iter) {
1001 if (!StringMatch(HEADER_RS_NODATA, iter->first)
1002 && std::modf(iter->second, &intpart) == 0.0) {
1003 // std::modf consider inf as an integer,
1004 // hence cannot handle -3.40282346639e+38 which is one of commonly used Nodata
1005 BSON_APPEND_INT32(&p, iter->first.c_str(), CVT_INT(iter->second));
1006 }
1007 else {
1008 BSON_APPEND_DOUBLE(&p, iter->first.c_str(), iter->second);
1009 }
1010 }
1011 // Check DATATYPE in metadata is consistent with T, transform if necessary
1012 char* buf = nullptr;
1013 vint buflength = -1;
1014 bool transformed = true;
1016 if (opt_type == temp_type) {
1017 buf = reinterpret_cast<char*>(values);
1018 buflength = datalength * sizeof(T);
1019 transformed = false;
1020 } else {
1021 if (opt_type == RDT_UInt8) {
1022 vuint8_t* newdata = nullptr;
1023 Initialize1DArray(datalength, newdata, values);
1024 buf = reinterpret_cast<char*>(newdata);
1025 buflength = datalength * sizeof(vuint8_t);
1026 } else if (opt_type == RDT_Int8) {
1027 vint8_t* newdata = nullptr;
1028 Initialize1DArray(datalength, newdata, values);
1029 buf = reinterpret_cast<char*>(newdata);
1030 buflength = datalength * sizeof(vint8_t);
1031 } else if (opt_type == RDT_UInt16) {
1032 vuint16_t* newdata = nullptr;
1033 Initialize1DArray(datalength, newdata, values);
1034 buf = reinterpret_cast<char*>(newdata);
1035 buflength = datalength * sizeof(vuint16_t);
1036 } else if (opt_type == RDT_Int16) {
1037 vint16_t* newdata = nullptr;
1038 Initialize1DArray(datalength, newdata, values);
1039 buf = reinterpret_cast<char*>(newdata);
1040 buflength = datalength * sizeof(vint16_t);
1041 } else if (opt_type == RDT_UInt32) {
1042 vuint32_t* newdata = nullptr;
1043 Initialize1DArray(datalength, newdata, values);
1044 buf = reinterpret_cast<char*>(newdata);
1045 buflength = datalength * sizeof(vuint32_t);
1046 } else if (opt_type == RDT_Int32) {
1047 vint32_t* newdata = nullptr;
1048 Initialize1DArray(datalength, newdata, values);
1049 buf = reinterpret_cast<char*>(newdata);
1050 buflength = datalength * sizeof(vint32_t);
1051 } else if (opt_type == RDT_Float) {
1052 float* newdata = nullptr;
1053 Initialize1DArray(datalength, newdata, values);
1054 buf = reinterpret_cast<char*>(newdata);
1055 buflength = datalength * sizeof(float);
1056 } else if (opt_type == RDT_Double) {
1057 double* newdata = nullptr;
1058 Initialize1DArray(datalength, newdata, values);
1059 buf = reinterpret_cast<char*>(newdata);
1060 buflength = datalength * sizeof(double);
1061 }
1062 }
1063 // Add user-specific key-values into metadata
1064 AppendStringOptionsToBson(&p, curopts);
1065
1066 int try_times = 0;
1067 bool gstatus = false;
1068 while (try_times <= 3) { // Try 3 times
1069 gstatus = gfs->WriteStreamData(filename, buf, buflength, &p);
1070 if (gstatus) { break; }
1071 SleepMs(2); // Sleep 0.002 sec and retry
1072 try_times++;
1073 }
1074 bson_destroy(&p);
1075 if (transformed) { delete[] buf; }
1076 return gstatus;
1077}
1078
1079#endif /* USE_MONGODB */
1080
1086public:
1088
1089 SubsetPositions(int srow, int erow, int scol, int ecol);
1090
1091 explicit SubsetPositions(SubsetPositions*& src, bool deep_copy = false);
1092
1094
1095 bool Initialization();
1096
1097 template <typename T>
1098 bool SetData(const int n, T* data) {
1099 if (n != n_cells) { return false; }
1100 if (nullptr == data) { return false; }
1101 if (1 != n_lyrs) { n_lyrs = 1; }
1102 if (nullptr != data_) {
1103 for (int i = 0; i < n_cells; i++) { data_[i] = CVT_DBL(data[i]); }
1104 } else {
1106 }
1107 usable = true;
1108 return true;
1109 }
1110
1111 template <typename T>
1112 bool Set2DData(const int n, const int lyr, T** data2d) {
1113 if (n != n_cells) { return false; }
1114 if (nullptr == data2d) { return false; }
1115 if (lyr != n_lyrs) { n_lyrs = lyr; }
1116 if (nullptr != data2d_) {
1117 for (int i = 0; i < n_cells; i++) {
1118 for (int j = 0; j < n_lyrs; j++) {
1119 data2d_[i][j] = CVT_DBL(data2d[i][j]);
1120 }
1121 }
1122 } else {
1124 }
1125 usable = true;
1126 return true;
1127 }
1128#ifdef USE_MONGODB
1129 bool ReadFromMongoDB(MongoGridFs* gfs, const string& fname,
1130 const STRING_MAP& opts = STRING_MAP());
1131#endif
1132 void GetHeader(double gxll, double gyll, int gnrows, double cellsize,
1133 double nodata, STRDBL_MAP& subheader);
1134 template <typename T>
1135 void Output(T nodata, vector<T*>& fulldata) {
1136 if (nullptr == data_ && nullptr == data2d_) { return; }
1137 int nrows = g_erow - g_srow + 1;
1138 int ncols = g_ecol - g_scol + 1;
1139 int fullsize = nrows * ncols;
1140 for (int ilyr = 0; ilyr < n_lyrs; ilyr++) {
1141 T* tmpdata = nullptr;
1142 Initialize1DArray(fullsize, tmpdata, nodata);
1143 for (int vi = 0; vi < n_cells; vi++) {
1144 //int j = local_pos_[vi][0] * ncols + local_pos_[vi][1];
1145 int j = local_posidx_[vi];
1146 if (n_lyrs > 1 && nullptr != data2d_) {
1147 tmpdata[j] = static_cast<T>(data2d_[vi][ilyr]);
1148 }
1149 else if (n_lyrs == 1 && nullptr != data_) {
1150 tmpdata[j] = static_cast<T>(data_[vi]);
1151 }
1152 }
1153 fulldata.emplace_back(tmpdata);
1154 }
1155 }
1156
1157 bool usable;
1164 bool alloc_;
1167 int* global_;
1168 double* data_;
1169 double** data2d_;
1170};
1171
1177template <typename T, typename MASK_T = T>
1179public:
1180 /************* Construct functions ***************/
1181
1185 explicit clsRasterData(bool is_2d = false);
1186
1190 clsRasterData(T* data, int cols, int rows, T nodata, double dx, double xll,
1191 double yll, const STRING_MAP& opts = STRING_MAP());
1192
1197 clsRasterData(T* data, int cols, int rows, T nodata, double dx, double xll,
1198 double yll, const string& srs);
1199
1203 clsRasterData(T** data2d, int cols, int rows, int nlayers, T nodata,
1204 double dx, double xll, double yll, const STRING_MAP& opts = STRING_MAP());
1205
1210 clsRasterData(T** data2d, int cols, int rows, int nlayers, T nodata,
1211 double dx, double xll, double yll, const string& srs);
1212
1224 explicit clsRasterData(const string& filename, bool calc_pos = false,
1225 clsRasterData<MASK_T>* mask = nullptr, bool use_mask_ext = true,
1226 double default_value = NODATA_VALUE, const STRING_MAP& opts = STRING_MAP());
1227
1241 static clsRasterData<T, MASK_T>* Init(const string& filename, bool calc_pos = false,
1242 clsRasterData<MASK_T>* mask = nullptr, bool use_mask_ext = true,
1243 double default_value = NODATA_VALUE, const STRING_MAP& opts = STRING_MAP());
1244
1256 explicit clsRasterData(vector<string>& filenames, bool calc_pos = false,
1257 clsRasterData<MASK_T>* mask = nullptr, bool use_mask_ext = true,
1258 double default_value = NODATA_VALUE, const STRING_MAP& opts = STRING_MAP());
1259
1275 static clsRasterData<T, MASK_T>* Init(vector<string>& filenames, bool calc_pos = false,
1276 clsRasterData<MASK_T>* mask = nullptr, bool use_mask_ext = true,
1277 double default_value = NODATA_VALUE, const STRING_MAP& opts = STRING_MAP());
1278
1282 clsRasterData(clsRasterData<MASK_T>* mask, T* values, int len,
1283 const STRING_MAP& opts = STRING_MAP());
1284
1288 clsRasterData(clsRasterData<MASK_T>* mask, T** values, int len, int lyrs,
1289 const STRING_MAP& opts = STRING_MAP());
1290
1291#ifdef USE_MONGODB
1292
1305 clsRasterData(MongoGridFs* gfs, const char* remote_filename, bool calc_pos = false,
1306 clsRasterData<MASK_T>* mask = nullptr, bool use_mask_ext = true,
1307 double default_value = NODATA_VALUE, const STRING_MAP& opts = STRING_MAP());
1308
1312 static clsRasterData<T, MASK_T>* Init(MongoGridFs* gfs, const char* remote_filename,
1313 bool calc_pos = false, clsRasterData<MASK_T>* mask = nullptr,
1314 bool use_mask_ext = true, double default_value = NODATA_VALUE,
1315 const STRING_MAP& opts = STRING_MAP());
1316#endif
1317
1332 explicit clsRasterData(clsRasterData<T, MASK_T>* another);
1333
1336
1337 /************* Read functions ***************/
1347 bool ReadFromFile(const string& filename, bool calc_pos = false, clsRasterData<MASK_T>* mask = nullptr,
1348 bool use_mask_ext = true, double default_value = NODATA_VALUE,
1349 const STRING_MAP& opts = STRING_MAP());
1354 bool ReadFromFiles(vector<string>& filenames, bool calc_pos = false, clsRasterData<MASK_T>* mask = nullptr,
1355 bool use_mask_ext = true, double default_value = NODATA_VALUE,
1356 const STRING_MAP& opts = STRING_MAP());
1357
1358#ifdef USE_MONGODB
1369 bool ReadFromMongoDB(MongoGridFs* gfs, const string& filename, bool calc_pos = false,
1370 clsRasterData<MASK_T>* mask = nullptr, bool use_mask_ext = true,
1371 double default_value = NODATA_VALUE,
1372 const STRING_MAP& opts = STRING_MAP());
1373
1374#endif /* USE_MONGODB */
1375 /************* Write functions ***************/
1376
1382 bool OutputToFile(const string& filename, bool out_origin = true);
1383
1392 bool OutputSubsetToFile(bool out_origin = false, bool out_combined = true,
1393 const string& outname = string(),
1394 const map<vint, vector<double> >& recls = map<vint, vector<double> >(),
1395 double default_value = NODATA_VALUE);
1396
1400 bool OutputAscFile(const string& filename);
1401
1402#ifdef USE_GDAL
1406 bool OutputFileByGdal(const string& filename);
1407#endif /* USE_GDAL */
1408
1409#ifdef USE_MONGODB
1410
1419 bool OutputToMongoDB(MongoGridFs* gfs, const string& filename = string(),
1420 const STRING_MAP& opts = STRING_MAP(),
1421 bool include_nodata = true,
1422 bool out_origin = true);
1423
1439 bool OutputSubsetToMongoDB(MongoGridFs* gfs, const string& filename = string(),
1440 const STRING_MAP& opts = STRING_MAP(),
1441 bool include_nodata = true,
1442 bool out_origin = false, bool out_combined = true,
1443 const map<vint, vector<double> >& recls = map<vint, vector<double> >(),
1444 double default_value = NODATA_VALUE);
1445
1446#endif /* USE_MONGODB */
1447
1448 /************************************************************************/
1449 /* Set information functions */
1450 /************************************************************************/
1451
1452 void SetHeader(const STRDBL_MAP& refers);
1453
1455 void SetCoreName(const string& name) { core_name_ = name; }
1456
1460 void SetValue(int row, int col, T value, int lyr = 1);
1461
1466 bool SetCalcPositions();
1467
1471 bool SetPositions(int len, int** pdata);
1472 bool SetPositions(int len, int* pdata);
1473
1478 bool SetUseMaskExt();
1479
1481 void SetDataType(RasterDataType const type) { rs_type_ = type; }
1482
1484 void SetDataType(const string& strtype) { rs_type_ = StringToRasterDataType(strtype); }
1485
1487 void SetOutDataType(RasterDataType type);
1488
1490 void SetOutDataType(const string& strtype);
1491
1493 void SetDefaultValue(const double defaultv) { default_value_ = defaultv; }
1494
1498 bool BuildSubSet(map<int, int> groups = map<int, int>());
1499
1503 bool ReleaseSubset();
1504
1508 bool RebuildSubSet(map<int, int> groups = map<int, int>());
1509
1510 /************************************************************************/
1511 /* Get information functions */
1512 /************************************************************************/
1513
1518 void CalculateStatistics();
1519
1524 void UpdateStatistics();
1525
1529 void ReleaseStatsMap2D();
1530
1538 double GetStatistics(string sindex, int lyr = 1);
1539
1547 void GetStatistics(string sindex, int* lyrnum, double** values);
1548
1553 double GetAverage(const int lyr = 1) { return GetStatistics(STATS_RS_MEAN, lyr); }
1554
1559 double GetMaximum(const int lyr = 1) { return GetStatistics(STATS_RS_MAX, lyr); }
1560
1565 double GetMinimum(const int lyr = 1) { return GetStatistics(STATS_RS_MIN, lyr); }
1566
1571 double GetStd(const int lyr = 1) { return GetStatistics(STATS_RS_STD, lyr); }
1572
1578 double GetRange(const int lyr = 1) { return GetStatistics(STATS_RS_RANGE, lyr); }
1579
1585 void GetAverage(int* lyrnum, double** values) { GetStatistics(STATS_RS_MEAN, lyrnum, values); }
1586
1591 void GetMaximum(int* lyrnum, double** values) { GetStatistics(STATS_RS_MAX, lyrnum, values); }
1592
1597 void GetMinimum(int* lyrnum, double** values) { GetStatistics(STATS_RS_MIN, lyrnum, values); }
1598
1603 void GetStd(int* lyrnum, double** values) { GetStatistics(STATS_RS_STD, lyrnum, values); }
1604
1609 void GetRange(int* lyrnum, double** values) { GetStatistics(STATS_RS_RANGE, lyrnum, values); }
1610
1615 void GetValidNumber(int* lyrnum, double** values) { GetStatistics(STATS_RS_VALIDNUM, lyrnum, values); }
1616
1622 int GetValidNumber(const int lyr = 1) { return CVT_INT(GetStatistics(STATS_RS_VALIDNUM, lyr)); }
1623
1624 int GetCellNumber() const { return n_cells_; }
1626 int GetDataLength() const { return n_cells_ < 0 || n_lyrs_ < 0 ? -1 : n_cells_ * n_lyrs_; }
1627 int GetCols() const { return CVT_INT(headers_.at(HEADER_RS_NCOLS)); }
1628 int GetRows() const { return CVT_INT(headers_.at(HEADER_RS_NROWS)); }
1629 double GetCellWidth() const { return headers_.at(HEADER_RS_CELLSIZE); }
1630
1632 double GetXllCenter() const {
1633 if (headers_.find(HEADER_RS_XLL) != headers_.end()) {
1634 return headers_.at(HEADER_RS_XLL);
1635 }
1636 if (headers_.find(HEADER_RS_XLLCOR) != headers_.end()) {
1637 return headers_.at(HEADER_RS_XLLCOR) + 0.5 * GetCellWidth();
1638 }
1639 return NODATA_VALUE;
1640 }
1641
1643 double GetYllCenter() const {
1644 if (headers_.find(HEADER_RS_YLL) != headers_.end()) {
1645 return headers_.at(HEADER_RS_YLL);
1646 }
1647 if (headers_.find(HEADER_RS_YLLCOR) != headers_.end()) {
1648 return headers_.at(HEADER_RS_YLLCOR) + 0.5 * GetCellWidth();
1649 }
1650 return NODATA_VALUE;
1651 }
1652
1653 int GetLayers() const { return n_lyrs_; }
1654 RasterDataType GetDataType() const { return rs_type_; }
1655 RasterDataType GetOutDataType() const { return rs_type_out_; }
1656 T GetNoDataValue() const { return no_data_value_; }
1657 double GetDefaultValue() const { return default_value_; }
1658
1664 int GetPosition(int row, int col);
1665
1667 int GetPosition(float x, float y);
1668
1670 int GetPosition(double x, double y);
1671
1673 map<int, SubsetPositions*>& GetSubset() { return subset_; }
1674
1678 bool GetRasterData(int* n_cells, T** data);
1679
1684 bool Get2DRasterData(int* n_cells, int* n_lyrs, T*** data);
1685
1687 const STRDBL_MAP& GetRasterHeader() const { return headers_; }
1688
1690 const STRDBL_MAP& GetStatistics() const { return stats_; }
1691
1693 const map<string, double *>& GetStatistics2D() const { return stats_2d_; };
1694
1696 string GetFilePath() const { return full_path_; }
1697
1699 string GetCoreName() const { return core_name_; }
1700
1706 void GetRasterPositionData(int* datalength, int*** positiondata);
1707
1708 void GetRasterPositionData(int* datalength, int** positiondata);
1709
1710 T* GetRasterDataPointer() const { return raster_; }
1711 int** GetRasterPositionDataPointer() const { return pos_data_; }
1712 int* GetRasterPositionIndexPointer() const { return pos_idx_; }
1713 T** Get2DRasterDataPointer() const { return raster_2d_; }
1714 const char* GetSrs();
1715 string GetSrsString();
1716 string GetOption(const char* key);
1717 const STRING_MAP& GetOptions() const { return options_; }
1718
1723 T GetValueByIndex(int cell_index, int lyr = 1);
1724
1730 void GetValueByIndex(int cell_index, T*& values);
1731
1736 T GetValue(int row, int col, int lyr = 1);
1737
1744 void GetValue(int row, int col, T*& values);
1745
1750 bool IsNoData(const int row, const int col, const int lyr = 1) {
1751 return FloatEqual(GetValue(row, col, lyr), no_data_value_);
1752 }
1753
1754 bool Is2DRaster() const { return is_2draster; }
1755 bool PositionsCalculated() const { return calc_pos_; }
1756 bool PositionsAllocated() const { return store_pos_; }
1757 bool MaskExtented() const { return use_mask_ext_; }
1758 bool StatisticsCalculated() const { return stats_calculated_; }
1759 bool Initialized() const { return initialized_; }
1760
1765 if ((!is_2draster && nullptr != raster_) || // Valid 1D raster
1766 (is_2draster && nullptr != raster_2d_)) // Valid 2D raster
1767 { return true; }
1768 StatusMessage("Error: Please initialize the raster object first.");
1769 return false;
1770 }
1771
1775 bool ValidateRowCol(const int row, const int col) {
1776 if ((row < 0 || row >= GetRows()) || (col < 0 || col >= GetCols())) {
1777 StatusMessage("The row must between 0 and " + ValueToString(GetRows() - 1) +
1778 ", and the col must between 0 and " + ValueToString(GetCols() - 1));
1779 return false;
1780 }
1781 return true;
1782 }
1783
1787 bool ValidateLayer(const int lyr) {
1788 if (lyr <= 0 || lyr > n_lyrs_) {
1789 StatusMessage("The layer must be 1 ");
1790 if (n_lyrs_ > 1) StatusMessage(" or between 1 and " +
1792 return false;
1793 }
1794 return true;
1795 }
1796
1800 bool ValidateIndex(const int idx) {
1801 if (idx < 0 || idx >= n_cells_) {
1802 StatusMessage("The index must between 0 and " + utils_string::ValueToString(n_cells_ - 1));
1803 return false;
1804 }
1805 return true;
1806 }
1807
1809 string GetFullFileName() const { return full_path_; }
1810
1812 clsRasterData<MASK_T>* GetMask() const { return mask_; }
1813
1817 void Copy(clsRasterData<T, MASK_T>* orgraster);
1818
1822 void ReplaceNoData(T replacedv);
1823
1827 void Reclassify(map<int, T> reclass_map);
1828
1829 /************* Utility functions ***************/
1830
1838 XY_COOR GetCoordinateByRowCol(int row, int col);
1839
1848 ROW_COL GetPositionByCoordinate(double x, double y, STRDBL_MAP* header = nullptr);
1849
1850private:
1854 void InitializeRasterClass(bool is_2d = false);
1855
1859 void InitializeReadFunction(const string& filename, bool calc_pos = false, clsRasterData<MASK_T>* mask = nullptr,
1860 bool use_mask_ext = true, double default_value = NODATA_VALUE,
1861 const STRING_MAP& opts = STRING_MAP());
1862
1874 bool ConstructFromSingleFile(const string& filename, bool calc_pos = false, clsRasterData<MASK_T>* mask = nullptr,
1875 bool use_mask_ext = true, double default_value = NODATA_VALUE,
1876 const STRING_MAP& opts = STRING_MAP());
1877
1886 int MaskAndCalculateValidPosition();
1887
1893 void CalculateValidPositionsFromGridData();
1894
1898 bool PrepareCombSubsetData(T**values, int* datalen, int* datalyrs,
1899 bool out_origin = false, bool include_nodata = true,
1900 const map<vint, vector<double> >&recls = map<vint, vector<double> >(),
1901 double default_value = NODATA_VALUE);
1902
1906 bool PrepareSubsetData(int sub_id, SubsetPositions* sub,
1907 T** values, int* datalen, int* datalyrs,
1908 bool out_origin = false, bool include_nodata = true,
1909 const map<vint, vector<double> >& recls = map<vint, vector<double> >(),
1910 double default_value = NODATA_VALUE);
1911
1915 bool OutputFullsizeToFiles(T* fullsizedata, int fsize, int datalyrs,
1916 const string& fullfilename, const STRDBL_MAP& header,
1917 const STRING_MAP& opts);
1918
1928 void AddOtherLayerRasterData(int row, int col, int cellidx, int lyr,
1929 STRDBL_MAP lyrheader, T* lyrdata);
1930
1934 void CheckDefaultValue() {
1935 if (FloatEqual(default_value_, NODATA_VALUE) && !FloatEqual(no_data_value_, NODATA_VALUE)) {
1936 default_value_ = CVT_DBL(no_data_value_);
1937 }
1938 }
1939
1943 clsRasterData& operator=(const clsRasterData&) { }
1944
1951 int n_cells_;
1953 int n_lyrs_;
1955 RasterDataType rs_type_;
1957 RasterDataType rs_type_out_;
1959 T no_data_value_;
1961 double default_value_;
1963 string full_path_;
1965 string core_name_;
1967 T* raster_;
1969 T** raster_2d_;
1971 int** pos_data_;
1973 int* pos_idx_;
1975 STRING_MAP options_;
1977 STRDBL_MAP headers_;
1979 STRDBL_MAP stats_;
1981 map<string, double *> stats_2d_;
1983 clsRasterData<MASK_T>* mask_;
1985 map<int, SubsetPositions*> subset_;
1987 bool initialized_;
1989 bool is_2draster;
1991 bool calc_pos_;
1993 bool store_pos_;
1995 bool use_mask_ext_;
1997 bool stats_calculated_;
1998};
1999
2000/******** Define common used raster types **************/
2001#ifndef IntRaster
2002#define IntRaster clsRasterData<int>
2003#endif
2004#ifndef FltRaster
2005#define FltRaster clsRasterData<float>
2006#endif
2007#ifndef DblRaster
2008#define DblRaster clsRasterData<double>
2009#endif
2010#ifndef IntIntRaster
2011#define IntIntRaster clsRasterData<int, int>
2012#endif
2013#ifndef FltIntRaster
2014#define FltIntRaster clsRasterData<float, int>
2015#endif
2016#ifndef DblIntRaster
2017#define DblIntRaster clsRasterData<double, int>
2018#endif
2019
2020/*******************************************************/
2021/************* Implementation Code Begin ***************/
2022/*******************************************************/
2023
2024/************* Construct functions ***************/
2025
2026template <typename T, typename MASK_T>
2027void clsRasterData<T, MASK_T>::InitializeRasterClass(bool is_2d /* = false */) {
2028 full_path_ = "";
2029 core_name_ = "";
2030 n_cells_ = -1;
2031 rs_type_ = RDT_Unknown;
2032 rs_type_out_ = RDT_Unknown;
2033 no_data_value_ = DefaultNoDataByType(TypeToRasterDataType(typeid(T)));
2034 //no_data_value_ = static_cast<T>(NODATA_VALUE); // Be careful of unsigned data type!
2035 default_value_ = NODATA_VALUE;
2036 raster_ = nullptr;
2037 pos_data_ = nullptr;
2038 pos_idx_ = nullptr;
2039 mask_ = nullptr;
2040 subset_ = map<int, SubsetPositions*>();
2041 n_lyrs_ = -1;
2042 is_2draster = is_2d;
2043 raster_2d_ = nullptr;
2044 calc_pos_ = false;
2045 store_pos_ = false;
2046 use_mask_ext_ = false;
2047 stats_calculated_ = false;
2048 headers_ = InitialHeader();
2049 options_ = InitialStrHeader();
2050 InitialStatsMap(stats_, stats_2d_);
2051 initialized_ = true;
2052}
2053
2054template <typename T, typename MASK_T>
2055void clsRasterData<T, MASK_T>::InitializeReadFunction(const string& filename, const bool calc_pos /* = false */,
2056 clsRasterData<MASK_T>* mask /* = nullptr */,
2057 const bool use_mask_ext /* = true */,
2058 double default_value /* NODATA_VALUE */,
2059 const STRING_MAP& opts /* = STRING_MAP() */) {
2060 if (!initialized_) { InitializeRasterClass(); }
2061 full_path_ = filename;
2062 mask_ = mask;
2063 calc_pos_ = calc_pos;
2064 use_mask_ext_ = use_mask_ext;
2065 if (nullptr == mask_) { use_mask_ext_ = false; }
2066 default_value_ = default_value;
2067 rs_type_out_ = RasterDataTypeInOptionals(opts);
2068 core_name_ = GetCoreFileName(full_path_);
2069 CopyStringMap(opts, options_);
2070}
2071
2072template <typename T, typename MASK_T>
2074 InitializeRasterClass(is_2d);
2075}
2076
2077template <typename T, typename MASK_T>
2078clsRasterData<T, MASK_T>::clsRasterData(T* data, int cols, int rows, T nodata,
2079 double dx, double xll, double yll,
2080 const string& srs) {
2082 UpdateStrHeader(opts, HEADER_RS_SRS, srs);
2083 clsRasterData(data, cols, rows, nodata, dx, xll, yll, opts);
2084}
2085
2086template <typename T, typename MASK_T>
2087clsRasterData<T, MASK_T>::clsRasterData(T* data, const int cols, const int rows, T nodata,
2088 const double dx, const double xll, const double yll,
2089 const STRING_MAP& opts /* = STRING_MAP() */) {
2090 InitializeRasterClass(false);
2091 rs_type_out_ = RasterDataTypeInOptionals(opts);
2092 raster_ = data;
2093 no_data_value_ = nodata;
2094 CopyStringMap(opts, options_);
2095 n_cells_ = cols * rows;
2096 n_lyrs_ = 1;
2097 UpdateHeader(headers_, HEADER_RS_NCOLS, cols);
2098 UpdateHeader(headers_, HEADER_RS_NROWS, rows);
2099 UpdateHeader(headers_, HEADER_RS_XLL, xll);
2100 UpdateHeader(headers_, HEADER_RS_YLL, yll);
2101 UpdateHeader(headers_, HEADER_RS_CELLSIZE, dx);
2102 UpdateHeader(headers_, HEADER_RS_NODATA, nodata);
2103 UpdateHeader(headers_, HEADER_RS_LAYERS, 1);
2104 UpdateHeader(headers_, HEADER_RS_CELLSNUM, n_cells_);
2105}
2106
2107template <typename T, typename MASK_T>
2108clsRasterData<T, MASK_T>::clsRasterData(T** data2d, const int cols, const int rows, const int nlayers,
2109 T nodata, const double dx, const double xll, const double yll,
2110 const string& srs) {
2112 UpdateStrHeader(opts, HEADER_RS_SRS, srs);
2113 clsRasterData(data2d, cols, rows, nlayers, nodata, dx, xll, yll, opts);
2114}
2115
2116template <typename T, typename MASK_T>
2117clsRasterData<T, MASK_T>::clsRasterData(T** data2d, const int cols, const int rows, const int nlayers,
2118 T nodata, const double dx, const double xll, const double yll,
2119 const STRING_MAP& opts /* = STRING_MAP() */) {
2120 InitializeRasterClass(true);
2121 rs_type_out_ = RasterDataTypeInOptionals(opts);
2122 raster_2d_ = data2d;
2123 no_data_value_ = nodata;
2124 CopyStringMap(opts, options_);
2125 n_cells_ = cols * rows;
2126 n_lyrs_ = nlayers;
2127 headers_[HEADER_RS_NCOLS] = cols;
2128 headers_[HEADER_RS_NROWS] = rows;
2129 headers_[HEADER_RS_XLL] = xll;
2130 headers_[HEADER_RS_YLL] = yll;
2131 headers_[HEADER_RS_CELLSIZE] = dx;
2132 headers_[HEADER_RS_NODATA] = no_data_value_;
2133 headers_[HEADER_RS_LAYERS] = n_lyrs_;
2134 headers_[HEADER_RS_CELLSNUM] = n_cells_;
2135}
2136
2137template <typename T, typename MASK_T>
2138clsRasterData<T, MASK_T>::clsRasterData(const string& filename, const bool calc_pos /* = false */,
2139 clsRasterData<MASK_T>* mask /* = nullptr */,
2140 const bool use_mask_ext /* = true */,
2141 double default_value /* NODATA_VALUE */,
2142 const STRING_MAP& opts /* = STRING_MAP() */) {
2143 InitializeRasterClass();
2144 rs_type_out_ = RasterDataTypeInOptionals(opts);
2145 ReadFromFile(filename, calc_pos, mask, use_mask_ext, default_value, opts);
2146}
2147
2148template <typename T, typename MASK_T>
2150 const bool calc_pos /* = false */,
2151 clsRasterData<MASK_T>* mask /* = nullptr */,
2152 const bool use_mask_ext /* = true */,
2153 double default_value /* = NODATA_VALUE */,
2154 const STRING_MAP& opts /* = STRING_MAP() */) {
2155 if (!FileExists(filename)) { return nullptr; }
2158 if (!rs->ReadFromFile(filename, calc_pos, mask, use_mask_ext, default_value, opts)) {
2159 delete rs;
2160 return nullptr;
2161 }
2162 return rs;
2163}
2164
2165template <typename T, typename MASK_T>
2167 const bool calc_pos /* = false */,
2168 clsRasterData<MASK_T>* mask /* = nullptr */,
2169 const bool use_mask_ext /* = true */,
2170 double default_value /* = NODATA_VALUE */,
2171 const STRING_MAP& opts /* = STRING_MAP() */) {
2172 if (!FilesExist(filenames)) {
2173 StatusMessage("Please make sure all file path existed!");
2174 return;
2175 }
2176 if (filenames.size() == 1) {
2177 ReadFromFile(filenames[0], calc_pos, mask,
2178 use_mask_ext, default_value, opts);
2179 }
2180 ReadFromFiles(filenames, calc_pos, mask,
2181 use_mask_ext, default_value, opts);
2182}
2183
2184template <typename T, typename MASK_T>
2186 bool calc_pos /* = false */,
2187 clsRasterData<MASK_T>* mask /* = nullptr */,
2188 bool use_mask_ext /* = true */,
2189 double default_value /* = NODATA_VALUE */,
2190 const STRING_MAP& opts /* = STRING_MAP() */) {
2191 if (!FilesExist(filenames)) {
2192 StatusMessage("Please make sure all file path existed!");
2193 return nullptr;
2194 }
2195 if (filenames.size() == 1) {
2196 return clsRasterData<T, MASK_T>::Init(filenames[0], calc_pos, mask, use_mask_ext,
2197 default_value, opts);
2198 }
2201 if (!rs2d->ReadFromFiles(filenames, calc_pos, mask, use_mask_ext, default_value, opts)) {
2202 delete rs2d;
2203 return nullptr;
2204 }
2205 return rs2d;
2206}
2207
2208template <typename T, typename MASK_T>
2210 const STRING_MAP& opts /* = STRING_MAP() */) {
2211 InitializeRasterClass(false);
2212 rs_type_out_ = RasterDataTypeInOptionals(opts);
2213 calc_pos_ = false;
2214 mask_ = mask;
2215 use_mask_ext_ = true;
2216 n_lyrs_ = 1;
2217 mask->GetRasterPositionData(&n_cells_, &pos_data_);
2218 mask->GetRasterPositionData(&n_cells_, &pos_idx_);
2219 if (n_cells_ != len) {
2220 StatusMessage("Input data length MUST EQUALS TO valid cell's number of mask!");
2221 initialized_ = false;
2222 return;
2223 }
2224 Initialize1DArray(n_cells_, raster_, values); // DO NOT ASSIGN ARRAY DIRECTLY!
2225 default_value_ = mask_->GetDefaultValue();
2226 CopyHeader(mask_->GetRasterHeader(), headers_);
2227 no_data_value_ = static_cast<T>(mask_->GetNoDataValue());
2228 UpdateStrHeader(options_, HEADER_RS_SRS, mask_->GetSrsString());
2229 CopyStringMap(opts, options_);
2230}
2231
2232template <typename T, typename MASK_T>
2234 const int lyrs, const STRING_MAP& opts /* = STRING_MAP() */) {
2235 InitializeRasterClass(true);
2236 calc_pos_ = false;
2237 mask_ = mask;
2238 use_mask_ext_ = true;
2239 n_lyrs_ = lyrs;
2240 mask->GetRasterPositionData(&n_cells_, &pos_data_);
2241 mask->GetRasterPositionData(&n_cells_, &pos_idx_);
2242 if (n_cells_ != len) {
2243 StatusMessage("Input data length MUST EQUALS TO valid cell's number of mask!");
2244 initialized_ = false;
2245 return;
2246 }
2247 Initialize2DArray(n_cells_, n_lyrs_, raster_2d_, values); // DO NOT ASSIGN ARRAY DIRECTLY!
2248 CopyHeader(mask_->GetRasterHeader(), headers_);
2249 no_data_value_ = static_cast<T>(mask_->GetNoDataValue());
2250 UpdateHeader(headers_, HEADER_RS_LAYERS, n_lyrs_);
2251 CopyStringMap(opts, options_);
2252 UpdateStrHeader(options_, HEADER_RS_SRS, mask_->GetSrsString());
2253}
2254
2255#ifdef USE_MONGODB
2256
2257template <typename T, typename MASK_T>
2259 const bool calc_pos /* = false */,
2260 clsRasterData<MASK_T>* mask /* = nullptr */,
2261 const bool use_mask_ext /* = true */,
2262 double default_value /* NODATA_VALUE */,
2263 const STRING_MAP& opts /* = STRING_MAP() */) {
2264 InitializeRasterClass();
2265 rs_type_out_ = RasterDataTypeInOptionals(opts);
2266 ReadFromMongoDB(gfs, remote_filename, calc_pos, mask, use_mask_ext, default_value, opts);
2267}
2268
2269template <typename T, typename MASK_T>
2271 const bool calc_pos /* = false */,
2272 clsRasterData<MASK_T>* mask /* = nullptr */,
2273 const bool use_mask_ext /* = true */,
2274 double default_value /* NODATA_VALUE */,
2275 const STRING_MAP& opts /* = STRING_MAP() */) {
2278 if (!rs_mongo->ReadFromMongoDB(gfs, remote_filename, calc_pos, mask, use_mask_ext, default_value, opts)) {
2279 delete rs_mongo;
2280 return nullptr;
2281 }
2282 return rs_mongo;
2283}
2284
2285#endif /* USE_MONGODB */
2286
2287template <typename T, typename MASK_T>
2288bool clsRasterData<T, MASK_T>::ConstructFromSingleFile(const string& filename,
2289 const bool calc_pos /* = false */,
2290 clsRasterData<MASK_T>* mask /* = nullptr */,
2291 const bool use_mask_ext /* = true */,
2292 double default_value /* NODATA_VALUE */,
2293 const STRING_MAP& opts /* = STRING_MAP() */) {
2294 InitializeReadFunction(filename, calc_pos, mask, use_mask_ext, default_value, opts);
2295
2296 bool readflag = false;
2297 string srs = string();
2298 if (StringMatch(GetUpper(GetSuffix(filename)), ASCIIExtension)) {
2299 readflag = ReadAscFile(full_path_, headers_, raster_);
2300 } else {
2301#ifdef USE_GDAL
2302 readflag = ReadRasterFileByGdal(full_path_, headers_, raster_, rs_type_, srs);
2303#else
2304 StatusMessage("Warning: Only ASC format is supported without GDAL!");
2305 return false;
2306#endif /* USE_GDAL */
2307 }
2308 // After read raster data from ASCII file or GeoTiff.
2309 no_data_value_ = static_cast<T>(headers_.at(HEADER_RS_NODATA));
2310 UpdateStrHeader(options_, HEADER_RS_SRS, srs);
2311 // if not specified, set output data type the same as input
2312 if (rs_type_out_ == 0) {
2313 rs_type_out_ = rs_type_;
2315 }
2316 CheckDefaultValue();
2317 if (readflag) {
2318 if (n_lyrs_ < 0) { n_lyrs_ = 1; }
2319 if (is_2draster) { is_2draster = false; }
2320 return MaskAndCalculateValidPosition() >= 0;
2321 }
2322 return false;
2323}
2324
2325template <typename T, typename MASK_T>
2327 if (!core_name_.empty()) { StatusMessage(("Release raster: " + core_name_).c_str()); }
2328 if (nullptr != raster_) { Release1DArray(raster_); }
2329 if (nullptr != pos_data_ && store_pos_) { Release2DArray(pos_data_); }
2330 if (nullptr != pos_idx_ && store_pos_) { Release1DArray(pos_idx_); }
2331 if (nullptr != raster_2d_ && is_2draster) { Release2DArray(raster_2d_); }
2332 if (is_2draster && stats_calculated_) { ReleaseStatsMap2D(); }
2333 ReleaseSubset();
2334}
2335
2336template <typename T, typename MASK_T>
2337bool clsRasterData<T, MASK_T>::BuildSubSet(map<int, int> groups /* = map<int, int>() */) {
2338 if (!ValidateRasterData()) { return false; }
2339 if (nullptr == pos_data_ || nullptr == pos_idx_) {
2340 if (!SetCalcPositions()) { return false; }
2341 }
2342 if (!subset_.empty()) { return true; }
2343
2344 int global_ncols = GetCols();
2345 map<int, vector<int> > global_idx;
2346 for (int vi = 0; vi < n_cells_; vi++) {
2347 T curv = GetValueByIndex(vi); // compatible with 2D Raster
2348 if (FloatEqual(curv, no_data_value_)) { continue; }
2349 int groupv = CVT_INT(curv); // By default, group value is original raster value
2350 if (!groups.empty() && groups.find(groupv) != groups.end()) {
2351 groupv = groups.at(groupv); // original raster value --> specified group ID
2352 }
2353 //int currow = pos_data_[vi][0];
2354 //int curcol = pos_data_[vi][1];
2355 int currow = pos_idx_[vi] / global_ncols;
2356 int curcol = pos_idx_[vi] % global_ncols;
2357 if (subset_.find(groupv) == subset_.end()) {
2358#ifdef HAS_VARIADIC_TEMPLATES
2359 subset_.emplace(groupv, new SubsetPositions(currow, currow, curcol, curcol));
2360#else
2361 subset_.insert(make_pair(groupv, new SubsetPositions(currow, currow, curcol, curcol)));
2362#endif
2363 }
2364 SubsetPositions*& cursubset = subset_.at(groupv);
2365 if (currow > cursubset->g_erow) { cursubset->g_erow = currow; }
2366 if (currow < cursubset->g_srow) { cursubset->g_srow = currow; }
2367 if (curcol > cursubset->g_ecol) { cursubset->g_ecol = curcol; }
2368 if (curcol < cursubset->g_scol) { cursubset->g_scol = curcol; }
2369 if (global_idx.find(groupv) == global_idx.end()) {
2370#ifdef HAS_VARIADIC_TEMPLATES
2371 global_idx.emplace(groupv, vector<int>());
2372#else
2373 global_idx.insert(make_pair(groupv, vector<int>()));
2374#endif
2375 }
2376 global_idx[groupv].emplace_back(vi);
2377 cursubset->n_cells += 1;
2378 }
2379 for (auto it = subset_.begin(); it != subset_.end(); ++it) {
2380 Initialize1DArray(it->second->n_cells, it->second->global_, 0);
2381 for (auto it2 = global_idx[it->first].begin(); it2 != global_idx[it->first].end(); ++it2) {
2382 it->second->global_[it2 - global_idx[it->first].begin()] = *it2;
2383 }
2384 }
2385 for (auto it = subset_.begin(); it != subset_.end(); ++it) {
2386 Initialize2DArray(it->second->n_cells, 2, it->second->local_pos_, -1);
2387 Initialize1DArray(it->second->n_cells, it->second->local_posidx_, -1);
2388
2389 int nrows = it->second->g_erow - it->second->g_srow + 1;
2390 int local_ncols = it->second->g_ecol - it->second->g_scol + 1;
2391 for (int gidx = 0; gidx < it->second->n_cells; gidx++) {
2392 //it->second->local_pos_[gidx][0] = pos_data_[it->second->global_[gidx]][0] - it->second->g_srow;
2393 //it->second->local_pos_[gidx][1] = pos_data_[it->second->global_[gidx]][1] - it->second->g_scol;
2394 int local_row = pos_idx_[it->second->global_[gidx]] / global_ncols - it->second->g_srow;
2395 int local_col = pos_idx_[it->second->global_[gidx]] % global_ncols - it->second->g_scol;
2396 it->second->local_pos_[gidx][0] = local_row;
2397 it->second->local_pos_[gidx][1] = local_col;
2398 it->second->local_posidx_[gidx] = local_row * local_ncols + local_col;
2399 }
2400 it->second->alloc_ = true;
2401 }
2402 return true;
2403}
2404
2405template <typename T, typename MASK_T>
2407 if (!subset_.empty()) {
2408 for (auto it = subset_.begin(); it != subset_.end(); ++it) {
2409 delete it->second;
2410 it->second = nullptr;
2411 }
2412 subset_.clear();
2413 }
2414 return true;
2415}
2416
2417template <typename T, typename MASK_T>
2418bool clsRasterData<T, MASK_T>::RebuildSubSet(const map<int, int> groups /* = map<int, int>() */) {
2419 return ReleaseSubset() && BuildSubSet(groups);
2420}
2421
2422/************* Set information functions ***************/
2423
2424template <typename T, typename MASK_T>
2426 rs_type_out_ = type;
2428}
2429
2430template <typename T, typename MASK_T>
2431void clsRasterData<T, MASK_T>::SetOutDataType(const string& strtype) {
2432 rs_type_out_ = StringToRasterDataType(strtype);
2434}
2435
2436/************* Get information functions ***************/
2437template <typename T, typename MASK_T>
2439 if (stats_calculated_) { return; }
2440 if (stats_.empty() || stats_2d_.empty()) { InitialStatsMap(stats_, stats_2d_); }
2441 if (is_2draster && nullptr != raster_2d_) {
2442 double** derivedvs = nullptr;
2443 BasicStatistics(raster_2d_, n_cells_, n_lyrs_, &derivedvs, no_data_value_);
2444 stats_2d_.at(STATS_RS_VALIDNUM) = derivedvs[0];
2445 stats_2d_.at(STATS_RS_MEAN) = derivedvs[1];
2446 stats_2d_.at(STATS_RS_MAX) = derivedvs[2];
2447 stats_2d_.at(STATS_RS_MIN) = derivedvs[3];
2448 stats_2d_.at(STATS_RS_STD) = derivedvs[4];
2449 stats_2d_.at(STATS_RS_RANGE) = derivedvs[5];
2450 // 1D array elements of derivedvs will be released by the destructor: releaseStatsMap2D()
2451 Release1DArray(derivedvs);
2452 } else {
2453 double* derivedv = nullptr;
2454 BasicStatistics(raster_, n_cells_, &derivedv, no_data_value_);
2455 stats_.at(STATS_RS_VALIDNUM) = derivedv[0];
2456 stats_.at(STATS_RS_MEAN) = derivedv[1];
2457 stats_.at(STATS_RS_MAX) = derivedv[2];
2458 stats_.at(STATS_RS_MIN) = derivedv[3];
2459 stats_.at(STATS_RS_STD) = derivedv[4];
2460 stats_.at(STATS_RS_RANGE) = derivedv[5];
2461 Release1DArray(derivedv);
2462 }
2463 stats_calculated_ = true;
2464}
2465
2466template <typename T, typename MASK_T>
2468 for (auto it = stats_2d_.begin(); it != stats_2d_.end(); ++it) {
2469 if (nullptr != it->second) {
2470 Release1DArray(it->second);
2471 }
2472 }
2473 stats_2d_.clear();
2474}
2475
2476template <typename T, typename MASK_T>
2478 if (is_2draster && stats_calculated_) ReleaseStatsMap2D();
2479 stats_calculated_ = false;
2480 CalculateStatistics();
2481}
2482
2483template <typename T, typename MASK_T>
2484double clsRasterData<T, MASK_T>::GetStatistics(string sindex, const int lyr /* = 1 */) {
2485 sindex = GetUpper(sindex);
2486 if (!ValidateRasterData() || !ValidateLayer(lyr)) {
2487 StatusMessage("No available raster statistics!");
2488 return CVT_DBL(default_value_);
2489 }
2490 if (is_2draster && nullptr != raster_2d_) {
2491 // for 2D raster data
2492 auto it = stats_2d_.find(sindex);
2493 if (it != stats_2d_.end()) {
2494 if (nullptr == it->second) {
2495 stats_calculated_ = false;
2496 CalculateStatistics();
2497 }
2498 return stats_2d_.at(sindex)[lyr - 1];
2499 }
2500 StatusMessage("WARNING: " + ValueToString(sindex) + " is not supported currently.");
2501 return CVT_DBL(default_value_);
2502 }
2503 // Else, for 1D raster data
2504 auto it = stats_.find(sindex);
2505 if (it != stats_.end()) {
2506 if (FloatEqual(it->second, CVT_DBL(NODATA_VALUE)) || !stats_calculated_) {
2507 CalculateStatistics();
2508 }
2509 return stats_.at(sindex);
2510 }
2511 StatusMessage("WARNING: " + ValueToString(sindex) + " is not supported currently.");
2512 return CVT_DBL(default_value_);
2513}
2514
2515template <typename T, typename MASK_T>
2516void clsRasterData<T, MASK_T>::GetStatistics(string sindex, int* lyrnum, double** values) {
2517 if (!is_2draster || nullptr != raster_2d_) {
2518 StatusMessage("Please initialize the raster object first.");
2519 *values = nullptr;
2520 return;
2521 }
2522 sindex = GetUpper(sindex);
2523 *lyrnum = n_lyrs_;
2524 auto it = stats_2d_.find(sindex);
2525 if (!Is2DRaster()) {
2526 *values = nullptr;
2527 return;
2528 }
2529 if (it == stats_2d_.end()) {
2530 *values = nullptr;
2531 StatusMessage("WARNING: " + ValueToString(sindex) + " is not supported currently.");
2532 return;
2533 }
2534 if (nullptr == it->second || !stats_calculated_) {
2535 CalculateStatistics();
2536 }
2537 *values = it->second;
2538}
2539
2540template <typename T, typename MASK_T>
2541int clsRasterData<T, MASK_T>::GetPosition(const int row, const int col) {
2542 if (!ValidateRasterData() || !ValidateRowCol(row, col)) {
2543 return -2; // means error occurred!
2544 }
2545 int pos_idx = GetCols() * row + col;
2546 if (!calc_pos_ || nullptr == pos_data_ || nullptr == pos_idx_) {
2547 return pos_idx;
2548 }
2549// previous low-efficiency code
2550// for (int i = 0; i < n_cells_; i++) {
2551// if (row == pos_data_[i][0] && col == pos_data_[i][1]) {
2552// return i;
2553// }
2554// }
2555 // Use binary search method, refers to https://leetcode.cn/problems/binary-search
2556 //int search(vector<int>& nums, int target) {
2557 int left = 0;
2558 int right = n_cells_ - 1; // assumes pos_idx belongs to pos_idx_[left, right]
2559 while (left <= right) { // when left==right,[left, right] still works,so use <=
2560 int middle = left + ((right - left) / 2); // equals to (left + right)/2
2561 if (pos_idx_[middle] > pos_idx) {
2562 right = middle - 1; // pos_idx belongs to [left, middle - 1]
2563 } else if (pos_idx_[middle] < pos_idx) {
2564 left = middle + 1; // pos_idx belongs to [middle + 1, right]
2565 } else { // pos_idx_[middle] == pos_idx
2566 return middle;
2567 }
2568 }
2569 return -1; // means the location of the raster data or mask data is NODATA
2570}
2571
2572template <typename T, typename MASK_T>
2573int clsRasterData<T, MASK_T>::GetPosition(const float x, const float y) {
2574 return GetPosition(CVT_DBL(x), CVT_DBL(y));
2575}
2576
2577template <typename T, typename MASK_T>
2578int clsRasterData<T, MASK_T>::GetPosition(const double x, const double y) {
2579 if (!initialized_) return -2;
2580 double xll_center = GetXllCenter();
2581 double yll_center = GetYllCenter();
2582 double dx = GetCellWidth();
2583 double dy = GetCellWidth();
2584 int n_rows = GetRows();
2585 int n_cols = GetCols();
2586
2587 if (FloatEqual(xll_center, CVT_DBL(NODATA_VALUE))
2588 || FloatEqual(yll_center, CVT_DBL(NODATA_VALUE))
2590 || n_rows < 0 || n_cols < 0) {
2591 StatusMessage("No available header information!");
2592 return -2;
2593 }
2594
2595 double x_min = xll_center - dx / 2.;
2596 double x_max = x_min + dx * n_cols;
2597 if (x > x_max || x < x_min) {
2598 return -2;
2599 }
2600
2601 double y_min = yll_center - dy / 2.;
2602 double y_max = y_min + dy * n_rows;
2603 if (y > y_max || y < y_min) {
2604 return -2;
2605 }
2606
2607 int n_row = CVT_INT((y_max - y) / dy); //calculate from ymax
2608 int n_col = CVT_INT((x - x_min) / dx); //calculate from xmin
2609
2610 return GetPosition(n_row, n_col);
2611}
2612
2613template <typename T, typename MASK_T>
2614bool clsRasterData<T, MASK_T>::GetRasterData(int* n_cells, T** data) {
2615 if (ValidateRasterData() && !is_2draster) {
2616 *n_cells = n_cells_;
2617 *data = raster_;
2618 return true;
2619 }
2620 *n_cells = -1;
2621 *data = nullptr;
2622 return false;
2623}
2624
2625template <typename T, typename MASK_T>
2626bool clsRasterData<T, MASK_T>::Get2DRasterData(int* n_cells, int* n_lyrs, T*** data) {
2627 if (ValidateRasterData() && is_2draster) {
2628 *n_cells = n_cells_;
2629 *n_lyrs = n_lyrs_;
2630 *data = raster_2d_;
2631 return true;
2632 }
2633 *n_cells = -1;
2634 *n_lyrs = -1;
2635 *data = nullptr;
2636 return false;
2637}
2638
2639template <typename T, typename MASK_T>
2640void clsRasterData<T, MASK_T>::GetRasterPositionData(int* datalength, int*** positiondata) {
2641 if (nullptr != pos_data_) {
2642 *datalength = n_cells_;
2643 *positiondata = pos_data_;
2644 } else {
2645 // reCalculate position data
2646 if (!ValidateRasterData()) {
2647 *datalength = -1;
2648 *positiondata = nullptr;
2649 return;
2650 }
2651 CalculateValidPositionsFromGridData();
2652 *datalength = n_cells_;
2653 *positiondata = pos_data_;
2654 }
2655}
2656
2657template <typename T, typename MASK_T>
2658void clsRasterData<T, MASK_T>::GetRasterPositionData(int* datalength, int** positiondata) {
2659 if (nullptr != pos_idx_) {
2660 *datalength = n_cells_;
2661 *positiondata = pos_idx_;
2662 } else {
2663 // reCalculate position data
2664 if (!ValidateRasterData()) {
2665 *datalength = -1;
2666 *positiondata = nullptr;
2667 return;
2668 }
2669 CalculateValidPositionsFromGridData();
2670 *datalength = n_cells_;
2671 *positiondata = pos_idx_;
2672 }
2673}
2674
2675template <typename T, typename MASK_T>
2677 return GetSrsString().c_str();
2678}
2679
2680template <typename T, typename MASK_T>
2682 return GetOption(HEADER_RS_SRS);
2683}
2684
2685template <typename T, typename MASK_T>
2687 if (options_.find(key) == options_.end()) {
2688 StatusMessage((string(key) + " is not existed in the options of the raster data!").c_str());
2689 return "";
2690 }
2691 return options_.at(key);
2692}
2693
2694template <typename T, typename MASK_T>
2695T clsRasterData<T, MASK_T>::GetValueByIndex(const int cell_index, const int lyr /* = 1 */) {
2696 if (!ValidateRasterData() || !ValidateIndex(cell_index) || !ValidateLayer(lyr)) {
2697 return no_data_value_;
2698 }
2699 if (is_2draster) {
2700 return raster_2d_[cell_index][lyr - 1];
2701 }
2702 return raster_[cell_index];
2703}
2704
2705template <typename T, typename MASK_T>
2706void clsRasterData<T, MASK_T>::GetValueByIndex(const int cell_index, T*& values) {
2707 if (!ValidateRasterData() || !ValidateIndex(cell_index)) {
2708 if (nullptr != values) { Release1DArray(values); }
2709 values = nullptr;
2710 return;
2711 }
2712 if (nullptr == values) {
2713 Initialize1DArray(n_lyrs_, values, no_data_value_);
2714 }
2715 if (is_2draster) {
2716 for (int i = 0; i < n_lyrs_; i++) {
2717 values[i] = raster_2d_[cell_index][i];
2718 }
2719 } else {
2720 values[0] = raster_[cell_index];
2721 }
2722}
2723
2724template <typename T, typename MASK_T>
2725T clsRasterData<T, MASK_T>::GetValue(const int row, const int col, const int lyr /* = 1 */) {
2726 if (!ValidateRasterData() || !ValidateRowCol(row, col) || !ValidateLayer(lyr)) {
2727 return no_data_value_;
2728 }
2729 // get index according to position data if possible
2730 if (calc_pos_ && (nullptr != pos_data_ || nullptr != pos_idx_)) {
2731 int valid_cell_index = GetPosition(row, col);
2732 if (valid_cell_index < 0) { return no_data_value_; }// error or NODATA
2733 return GetValueByIndex(valid_cell_index, lyr);
2734 }
2735 // get data directly from row and col
2736 if (is_2draster) { return raster_2d_[row * GetCols() + col][lyr - 1]; }
2737 return raster_[row * GetCols() + col];
2738}
2739
2740template <typename T, typename MASK_T>
2741void clsRasterData<T, MASK_T>::GetValue(const int row, const int col, T*& values) {
2742 if (!ValidateRasterData() || !ValidateRowCol(row, col)) {
2743 if (nullptr != values) { Release1DArray(values); }
2744 values = nullptr;
2745 return;
2746 }
2747 if (nullptr == values) {
2748 Initialize1DArray(n_lyrs_, values, no_data_value_);
2749 }
2750 if (calc_pos_ && (nullptr != pos_data_ || nullptr != pos_idx_)) {
2751 int valid_cell_index = GetPosition(row, col);
2752 if (valid_cell_index == -1) {
2753 for (int i = 0; i < n_lyrs_; i++) {
2754 values[i] = no_data_value_; // NODATA
2755 }
2756 } else {
2757 GetValueByIndex(valid_cell_index, values);
2758 }
2759 } else {
2760 // get data directly from row and col
2761 if (is_2draster) {
2762 for (int i = 0; i < n_lyrs_; i++) {
2763 values[i] = raster_2d_[row * GetCols() + col][i];
2764 }
2765 } else {
2766 values[0] = raster_[row * GetCols() + col];
2767 }
2768 }
2769}
2770
2771template <typename T, typename MASK_T>
2773 CopyHeader(refers, headers_);
2774 // Update header related variables
2775 auto it = headers_.find(HEADER_RS_CELLSNUM);
2776 if (it != headers_.end()) { n_cells_ = CVT_INT(it->second); }
2777 it = headers_.find(HEADER_RS_LAYERS);
2778 if (it != headers_.end()) { n_lyrs_ = CVT_INT(it->second); }
2779 it = headers_.find(HEADER_RS_NODATA);
2780 if (it != headers_.end()) { no_data_value_ = static_cast<T>(it->second); }
2781}
2782
2783
2784template <typename T, typename MASK_T>
2785void clsRasterData<T, MASK_T>::SetValue(const int row, const int col, T value, const int lyr /* = 1 */) {
2786 if (!ValidateRasterData() || !ValidateRowCol(row, col) || !ValidateLayer(lyr)) {
2787 StatusMessage("Set value failed!");
2788 return;
2789 }
2790 int idx = GetPosition(row, col);
2791 if (idx == -1) {
2792 // the origin value is NODATA, and positions of valid values are calculated
2793 StatusMessage("Current version do not support to setting value to NoDATA location!");
2794 } else {
2795 if (is_2draster) {
2796 raster_2d_[idx][lyr - 1] = value;
2797 } else {
2798 raster_[idx] = value;
2799 }
2800 }
2801}
2802
2803template <typename T, typename MASK_T>
2805 if (!ValidateRasterData()) { return false; }
2806 if (calc_pos_ && (nullptr != pos_data_ || nullptr != pos_idx_)) {
2807 // already set as True, no need to recalculate.
2808 return false;
2809 }
2810 calc_pos_ = true;
2811 return MaskAndCalculateValidPosition() >= 0;
2812}
2813
2814template <typename T, typename MASK_T>
2815bool clsRasterData<T, MASK_T>::SetPositions(int len, int** pdata) {
2816 if (nullptr != pos_data_) {
2817 if (len != n_cells_) { return false; } // cannot change origin n_cells_
2818 Release2DArray(pos_data_);
2819 }
2820 pos_data_ = pdata;
2821 calc_pos_ = true;
2822 store_pos_ = false;
2823 return true;
2824}
2825
2826template <typename T, typename MASK_T>
2827bool clsRasterData<T, MASK_T>::SetPositions(int len, int* pdata) {
2828 if (nullptr != pos_idx_) {
2829 if (len != n_cells_) { return false; } // cannot change origin n_cells_
2830 Release1DArray(pos_idx_);
2831 }
2832 pos_idx_ = pdata;
2833 calc_pos_ = true;
2834 store_pos_ = false;
2835 return true;
2836}
2837
2838template <typename T, typename MASK_T>
2840 if (!ValidateRasterData()) return false;
2841 if (nullptr == mask_) return false;
2842 if (use_mask_ext_) return false; // already set as True, no need to recalculate.
2843 use_mask_ext_ = true;
2844 return MaskAndCalculateValidPosition() >= 0;
2845}
2846
2847/************* Output to file functions ***************/
2848
2849template <typename T, typename MASK_T>
2850bool clsRasterData<T, MASK_T>::OutputToFile(const string& filename, const bool out_origin /* = true */) {
2851 if (GetPathFromFullName(filename).empty()) { return false; }
2852 string abs_filename = GetAbsolutePath(filename);
2853 if (!ValidateRasterData()) { return false; }
2854 if (!out_origin) { return OutputSubsetToFile(false, true, filename); }
2855 string filetype = GetUpper(GetSuffix(abs_filename));
2856 if (StringMatch(filetype, ASCIIExtension)) {
2857 return OutputAscFile(abs_filename);
2858 }
2859#ifdef USE_GDAL
2860 if (StringMatch(filetype, GTiffExtension)) {
2861 return OutputFileByGdal(abs_filename);
2862 }
2863 return OutputFileByGdal(ReplaceSuffix(abs_filename, GTiffExtension));
2864#else
2865 StatusMessage("Warning: Without GDAL, ASC file will be exported as default!");
2866 return OutputAscFile(ReplaceSuffix(abs_filename, ASCIIExtension));
2867#endif /* USE_GDAL */
2868}
2869
2870template <typename T, typename MASK_T>
2871bool clsRasterData<T, MASK_T>::OutputSubsetToFile(const bool out_origin /* = false */,
2872 const bool out_combined /* = true */,
2873 const string& outname /* = string() */,
2874 const map<vint, vector<double> >& recls /* map() */,
2875 const double default_value /* = NODATA_VALUE */) {
2876 if (!ValidateRasterData()) { return false; }
2877 if (subset_.empty()) { return false; }
2878 string outpathact = outname.empty() ? full_path_ : outname;
2879 // output combined subsets
2880 bool out_comb = out_combined;
2881 if (out_comb) {
2882 T* data1d = nullptr;
2883 int sublyrs;
2884 int sublen;
2885 if (!PrepareCombSubsetData(&data1d, &sublen, &sublyrs,
2886 out_origin, true, recls, default_value)) {
2887 return false;
2888 }
2889 STRDBL_MAP tmpheader;
2890 CopyHeader(headers_, tmpheader);
2891 UpdateHeader(tmpheader, HEADER_RS_LAYERS, sublyrs);
2892 UpdateHeader(tmpheader, HEADER_RS_CELLSNUM, sublen / sublyrs);
2893 bool combflag = OutputFullsizeToFiles(data1d, sublen / sublyrs, sublyrs,
2894 PrefixCoreFileName(outpathact, 0),
2895 tmpheader, options_);
2896 Release1DArray(data1d);
2897 return combflag;
2898 }
2899 STRDBL_MAP subheader;
2900 bool flag = true;
2901 for (auto it = subset_.begin(); it != subset_.end(); ++it) {
2902 if (!it->second->usable) { continue; }
2903 it->second->GetHeader(GetXllCenter(), GetYllCenter(), GetRows(),
2904 GetCellWidth(), CVT_DBL(no_data_value_), subheader);
2905 T* tmpdata1d = nullptr;
2906 int tmpdatalen;
2907 int tmplyrs;
2908 if (!PrepareSubsetData(it->first, it->second, &tmpdata1d, &tmpdatalen, &tmplyrs,
2909 out_origin, true, recls, default_value)) {
2910 continue;
2911 }
2912 UpdateHeader(subheader, HEADER_RS_LAYERS, tmplyrs);
2913 UpdateHeader(subheader, HEADER_RS_CELLSNUM, tmpdatalen / tmplyrs);
2914 flag = flag && OutputFullsizeToFiles(tmpdata1d, tmpdatalen / tmplyrs, tmplyrs,
2915 PrefixCoreFileName(outpathact, it->first),
2916 subheader, options_);
2917 Release1DArray(tmpdata1d);
2918 }
2919 return flag;
2920}
2921
2922template <typename T, typename MASK_T>
2923bool clsRasterData<T, MASK_T>::PrepareCombSubsetData(T** values, int* datalen, int* datalyrs,
2924 bool out_origin /* false */, bool include_nodata /* true */,
2925 const map<vint, vector<double> >& recls /* map() */,
2926 double default_value /* NODATA_VALUE*/) {
2927 if (subset_.empty()) { return false; }
2928 T* data1d = nullptr; // Both raster 1D and 2D data can be combined as 1D array
2929 int lyr_recls = -1;
2930 if (!recls.empty()) {
2931 for (auto it = recls.begin(); it != recls.end(); ++it) {
2932 if (lyr_recls < 0) { lyr_recls = CVT_INT(it->second.size()); }
2933 else if (lyr_recls != it->second.size()) {
2934 StatusMessage("Error: Reclassification layer count MUST be consistent!");
2935 return false;
2936 }
2937 }
2938 }
2939 int lyrs_subset = -1;
2940 for (auto it = subset_.begin(); it != subset_.end(); ++it) {
2941 if (!it->second->usable) { continue; }
2942 if (!(nullptr != it->second->data_ // Only if all subset have data_
2943 || nullptr != it->second->data2d_ // or data2d_,
2944 || !recls.empty())) { // or reclassification map specified
2945 return false;
2946 }
2947 if (lyrs_subset < 0) { lyrs_subset = it->second->n_lyrs; }
2948 else if (it->second->n_lyrs != lyrs_subset) {
2949 StatusMessage("Error: Subset's layer count MUST be consistent!");
2950 return false;
2951 }
2952 }
2953 int lyrs = -1;
2954 if (out_origin) {
2955 lyrs = n_lyrs_;
2956 if (lyr_recls > 0) { lyrs = lyr_recls; }
2957 } else {
2958 lyrs = lyrs_subset;
2959 if (lyrs < 0 && lyr_recls > 0) {
2960 lyrs = lyr_recls;
2961 }
2962 }
2963 if (lyrs < 0) {
2964 StatusMessage("Error: Cannot determine valid layer count!");
2965 return false;
2966 }
2967
2968 int gnrows = GetRows();
2969 int gncols = GetCols();
2970 if (FloatEqual(default_value, NODATA_VALUE)) {
2971 default_value = default_value_;
2972 }
2973 int gncells = include_nodata ? gnrows * gncols : n_cells_;
2974 int data_length = gncells * lyrs;
2975 Initialize1DArray(data_length, data1d, no_data_value_);
2976 for (auto it = subset_.begin(); it != subset_.end(); ++it) {
2977 bool use_defaultv_directly = false;
2978 if (!it->second->usable) {
2979 if (!FloatEqual(no_data_value_, default_value)) {
2980 use_defaultv_directly = true;
2981 } else {
2982 continue;
2983 }
2984 }
2985 for (int vi = 0; vi < it->second->n_cells; vi++) {
2986 for (int ilyr = 0; ilyr < lyrs; ilyr++) {
2987 int gidx = it->second->global_[vi];
2988 //int tmpr = pos_data_[gidx][0];
2989 //int tmpc = pos_data_[gidx][1];
2990 //int tmprc = tmpr * gncols + tmpc;
2991 int tmprc = pos_idx_[gidx];
2992 if (!include_nodata) { tmprc = gidx; }
2993 if (!recls.empty()) { // first priority
2994 double uniqe_value = default_value;
2995 int recls_key = it->first; // default reclassification key is subset's ID
2996 if (out_origin) {
2997 if (nullptr != raster_) { recls_key = CVT_INT(raster_[gidx]); }
2998 else if (nullptr != raster_2d_) {
2999 recls_key = CVT_INT(raster_2d_[gidx][ilyr]);
3000 }
3001 }
3002 if (recls.count(recls_key) > 0 && recls.at(recls_key).size() > ilyr) {
3003 uniqe_value = recls.at(recls_key).at(ilyr);
3004 }
3005 data1d[tmprc * lyrs + ilyr] = static_cast<T>(uniqe_value);
3006 }
3007 else if (use_defaultv_directly) {
3008 data1d[tmprc * lyrs + ilyr] = static_cast<T>(default_value);
3009 }
3010 else if (!out_origin && lyrs > 1 && nullptr != it->second->data2d_) { // raster 2D
3011 data1d[tmprc * lyrs + ilyr] = static_cast<T>(it->second->data2d_[vi][ilyr]);
3012 }
3013 else if (!out_origin && lyrs == 1 && nullptr != it->second->data_) { // raster 1D
3014 data1d[tmprc * lyrs + ilyr] = static_cast<T>(it->second->data_[vi]);
3015 }
3016 else { // Will not happen
3017 StatusMessage("Error: No subset or reclassification map can be output!");
3018 if (nullptr != data1d) { Release1DArray(data1d); }
3019 return false;
3020 }
3021 }
3022 }
3023 }
3024 *values = data1d;
3025 *datalen = data_length;
3026 *datalyrs = lyrs;
3027 return true;
3028}
3029
3030template <typename T, typename MASK_T>
3031bool clsRasterData<T, MASK_T>::PrepareSubsetData(const int sub_id, SubsetPositions* sub,
3032 T** values, int* datalen, int* datalyrs,
3033 bool out_origin /* false */, bool include_nodata /* true */,
3034 const map<vint, vector<double> >& recls /* map() */,
3035 double default_value /* NODATA_VALUE*/) {
3036 if (nullptr == sub) { return false; }
3037 T* data1d = nullptr; // Both raster 1D and 2D data can be combined as 1D array
3038 int lyr_recls = -1;
3039 if (!recls.empty()) {
3040 for (auto it = recls.begin(); it != recls.end(); ++it) {
3041 if (lyr_recls < 0) { lyr_recls = CVT_INT(it->second.size()); }
3042 else if (lyr_recls != it->second.size()) {
3043 StatusMessage("Error: Reclassification layer count MUST be consistent!");
3044 return false;
3045 }
3046 }
3047 }
3048 int lyrs = - 1;
3049 if (out_origin) {
3050 lyrs = n_lyrs_;
3051 if (lyr_recls > 0) { lyrs = lyr_recls; }
3052 } else {
3053 lyrs = sub->n_lyrs;
3054 }
3055 if (lyrs < 0) {
3056 StatusMessage("Error: Cannot determine valid layer count!");
3057 return false;
3058 }
3059 int nrows = sub->g_erow - sub->g_srow + 1;
3060 int ncols = sub->g_ecol - sub->g_scol + 1;
3061 if (FloatEqual(default_value, NODATA_VALUE)) {
3062 default_value = default_value_;
3063 }
3064 int ncells = include_nodata ? nrows * ncols : sub->n_cells;
3065 int data_length = ncells * lyrs;
3066 Initialize1DArray(data_length, data1d, no_data_value_);
3067 for (int vi = 0; vi < sub->n_cells; vi++) {
3068 for (int ilyr = 0; ilyr < lyrs; ilyr++) {
3069 //int j = sub->local_pos_[vi][0] * ncols + sub->local_pos_[vi][1];
3070 int j = sub->local_posidx_[vi];
3071 int gidx = sub->global_[vi];
3072 if (!include_nodata) { j = vi; }
3073 if (!recls.empty()) { // first priority
3074 double uniqe_value = default_value;
3075 int recls_key = sub_id;
3076 if (out_origin) {
3077 if (nullptr != raster_) { recls_key = CVT_INT(raster_[gidx]); }
3078 else if (nullptr != raster_2d_) {
3079 recls_key = CVT_INT(raster_2d_[gidx][ilyr]);
3080 }
3081 }
3082 if (recls.count(recls_key) > 0 && recls.at(recls_key).size() > ilyr) {
3083 uniqe_value = recls.at(recls_key).at(ilyr);
3084 }
3085 data1d[j * lyrs + ilyr] = static_cast<T>(uniqe_value);
3086 }
3087 else if (out_origin && lyrs > 1 && nullptr != raster_2d_) {
3088 data1d[j * lyrs + ilyr] = raster_2d_[gidx][ilyr];
3089 }
3090 else if (out_origin && lyrs == 1 && nullptr != raster_) {
3091 data1d[j * lyrs + ilyr] = raster_[gidx];
3092 }
3093 else if (!out_origin && lyrs > 1 && nullptr != sub->data2d_) { // raster 2D
3094 data1d[j * lyrs + ilyr] = static_cast<T>(sub->data2d_[vi][ilyr]);
3095 }
3096 else if (!out_origin && lyrs == 1 && nullptr != sub->data_) { // raster 1D
3097 data1d[j * lyrs + ilyr] = static_cast<T>(sub->data_[vi]);
3098 }
3099 else { // Will not happen
3100 StatusMessage("Error: No subset or reclassification map can be output!");
3101 if (nullptr != data1d) { Release1DArray(data1d); }
3102 return false;
3103 }
3104 }
3105 }
3106 *values = data1d;
3107 *datalen = data_length;
3108 *datalyrs = lyrs;
3109 return true;
3110}
3111
3112template <typename T, typename MASK_T>
3113bool clsRasterData<T, MASK_T>::OutputFullsizeToFiles(T* fullsizedata, const int fsize, const int datalyrs,
3114 const string& fullfilename,
3115 const STRDBL_MAP& header, const STRING_MAP& opts) {
3116 if (nullptr == fullsizedata) { return false; }
3117 if (datalyrs < 1) { return false; }
3118 if (datalyrs == 1) {
3119 return WriteRasterToFile(fullfilename, header, opts, fullsizedata);
3120 }
3121 T* tmpdata1d = nullptr;
3122 Initialize1DArray(fsize, tmpdata1d, NODATA_VALUE);
3123 bool flag = true;
3124 for (int ilyr = 0; ilyr < datalyrs; ilyr++) {
3125 for (int gi = 0; gi < fsize; gi++) {
3126 tmpdata1d[gi] = fullsizedata[gi * datalyrs + ilyr];
3127 }
3128 flag = flag && WriteRasterToFile(AppendCoreFileName(fullfilename, ilyr + 1),
3129 header, opts, tmpdata1d);
3130 }
3131 if (nullptr != tmpdata1d) { Release1DArray(tmpdata1d); }
3132 return flag;
3133}
3134
3135template <typename T, typename MASK_T>
3136bool clsRasterData<T, MASK_T>::OutputAscFile(const string& filename) {
3137 string abs_filename = GetAbsolutePath(filename);
3138 // Is there need to calculate valid position index?
3139 int count;
3140 int** position = nullptr;
3141 int* position_idx = nullptr;
3142 bool outputdirectly = true;
3143 if ((nullptr != pos_data_ || nullptr != pos_idx_)) {
3144 GetRasterPositionData(&count, &position);
3145 GetRasterPositionData(&count, &position_idx);
3146 outputdirectly = false;
3147 assert(nullptr != position);
3148 assert(nullptr != position_idx);
3149 }
3150 // Begin to write raster data
3151 int rows = CVT_INT(headers_.at(HEADER_RS_NROWS));
3152 int cols = CVT_INT(headers_.at(HEADER_RS_NCOLS));
3153 if (is_2draster) { // 3.1 2D raster data
3154 string pre_path = GetPathFromFullName(abs_filename);
3155 if (StringMatch(pre_path, "")) { return false; }
3156 for (int lyr = 0; lyr < n_lyrs_; lyr++) {
3157 string tmpfilename = AppendCoreFileName(abs_filename, itoa(CVT_VINT(lyr)));
3158 if (!WriteAscHeaders(tmpfilename, headers_)) { return false; }
3159 std::ofstream raster_file(tmpfilename.c_str(), std::ios::app | std::ios::out);
3160 if (!raster_file.is_open()) {
3161 StatusMessage("Error opening file: " + tmpfilename);
3162 return false;
3163 }
3164 int index = 0;
3165 for (int i = 0; i < rows; ++i) {
3166 for (int j = 0; j < cols; ++j) {
3167 if (outputdirectly) {
3168 index = i * cols + j;
3169 raster_file << setprecision(6) << raster_2d_[index][lyr] << " ";
3170 continue;
3171 }
3172 if (index < n_cells_ && (position[index][0] == i && position[index][1] == j)) {
3173 raster_file << setprecision(6) << raster_2d_[index][lyr] << " ";
3174 index++;
3175 } else { raster_file << setprecision(6) << NODATA_VALUE << " "; }
3176 }
3177 raster_file << endl;
3178 }
3179 raster_file.close();
3180 }
3181 } else { // 1D raster data
3182 if (!WriteAscHeaders(abs_filename, headers_)) { return false; }
3183 std::ofstream raster_file(filename.c_str(), std::ios::app | std::ios::out);
3184 if (!raster_file.is_open()) {
3185 StatusMessage("Error opening file: " + abs_filename);
3186 return false;
3187 }
3188 int index = 0;
3189 for (int i = 0; i < rows; ++i) {
3190 for (int j = 0; j < cols; ++j) {
3191 if (outputdirectly) {
3192 index = i * cols + j;
3193 raster_file << setprecision(6) << raster_[index] << " ";
3194 continue;
3195 }
3196 if (index < n_cells_) {
3197 if (position[index][0] == i && position[index][1] == j) {
3198 raster_file << setprecision(6) << raster_[index] << " ";
3199 index++;
3200 } else { raster_file << setprecision(6) << no_data_value_ << " "; }
3201 } else { raster_file << setprecision(6) << no_data_value_ << " "; }
3202 }
3203 raster_file << endl;
3204 }
3205 raster_file.close();
3206 }
3207 position = nullptr;
3208 return true;
3209}
3210
3211#ifdef USE_GDAL
3212template <typename T, typename MASK_T>
3214 string abs_filename = GetAbsolutePath(filename);
3215 bool outputdirectly = (nullptr == pos_data_ || nullptr == pos_idx_);
3216 int n_rows = CVT_INT(headers_.at(HEADER_RS_NROWS));
3217 int n_cols = CVT_INT(headers_.at(HEADER_RS_NCOLS));
3218 bool outflag = false;
3219 T* data_1d = nullptr;
3220 if (is_2draster) {
3221 string pre_path = GetPathFromFullName(abs_filename);
3222 if (StringMatch(pre_path, "")) { return false; }
3223 for (int lyr = 0; lyr < n_lyrs_; lyr++) {
3224 string tmpfilename = AppendCoreFileName(abs_filename, lyr + 1);
3225 if (outputdirectly) {
3226 if (nullptr == data_1d) {
3227 Initialize1DArray(n_rows * n_cols, data_1d, no_data_value_);
3228 }
3229 for (int gi = 0; gi < n_rows * n_cols; gi++) {
3230 data_1d[gi] = raster_2d_[gi][lyr];
3231 }
3232 } else {
3233 if (nullptr == data_1d) {
3234 Initialize1DArray(n_rows * n_cols, data_1d, no_data_value_);
3235 }
3236 for (int vi = 0; vi < n_cells_; vi++) {
3237 //data_1d[pos_data_[vi][0] * n_cols + pos_data_[vi][1]] = raster_2d_[vi][lyr];
3238 data_1d[pos_idx_[vi]] = raster_2d_[vi][lyr];
3239 }
3240 }
3241 outflag = WriteSingleGeotiff(tmpfilename, headers_, options_, data_1d);
3242 if (!outflag) {
3243 Release1DArray(data_1d);;
3244 return false;
3245 }
3246 }
3247 Release1DArray(data_1d);
3248 } else {
3249 if (outputdirectly) {
3250 outflag = WriteSingleGeotiff(abs_filename, headers_, options_, raster_);
3251 } else {
3252 Initialize1DArray(n_rows * n_cols, data_1d, no_data_value_);
3253 for (int vi = 0; vi < n_cells_; vi++) {
3254 //data_1d[pos_data_[vi][0] * n_cols + pos_data_[vi][1]] = raster_[vi];
3255 data_1d[pos_idx_[vi]] = raster_[vi];
3256 }
3257 outflag = WriteSingleGeotiff(abs_filename, headers_, options_, data_1d);
3258 Release1DArray(data_1d);
3259 }
3260 if (!outflag) { return false; }
3261 }
3262 return true;
3263}
3264#endif /* USE_GDAL */
3265
3266#ifdef USE_MONGODB
3267template <typename T, typename MASK_T>
3268bool clsRasterData<T, MASK_T>::OutputToMongoDB(MongoGridFs* gfs, const string& filename /* = string() */,
3269 const STRING_MAP& opts /* = STRING_MAP() */,
3270 bool include_nodata /* = true */,
3271 bool out_origin /* true */) {
3272 if (nullptr == gfs) { return false; }
3273 if (!out_origin) { // Output subset's data
3274 return OutputSubsetToMongoDB(gfs, filename, opts, include_nodata, false, true);
3275 }
3276 CopyStringMap(opts, options_); // Update metadata
3277 // Added by ljzhu, for compatible with yjwang's code. But, can this key-value be passed by the opts argument?
3278 UpdateStringMapIfNotExist(options_, HEADER_RS_PARAM_ABSTRACTION_TYPE, PARAM_ABSTRACTION_TYPE_PHYSICAL);
3279 if (options_.find(HEADER_RSOUT_DATATYPE) == options_.end()
3280 || StringMatch("Unknown", options_.at(HEADER_RSOUT_DATATYPE))) {
3283 }
3284 // Check if we can output directly: 1) pos_data_ is not NULL and include_nodata is false;
3285 // 2) pos_data_ is NULL and include_nodata is true.
3286 bool outputdirectly = true; // output directly or create new full size array
3287 int cnt;
3288 int** pos = nullptr;
3289 if ((nullptr != pos_data_ || nullptr != pos_idx_) && include_nodata) {
3290 outputdirectly = false;
3291 GetRasterPositionData(&cnt, &pos);
3292 }
3293 if ((nullptr == pos_data_ || nullptr == pos_idx_) && !include_nodata) {
3294 SetCalcPositions();
3295 GetRasterPositionData(&cnt, &pos);
3296 }
3297 int n_rows = GetRows();
3298 int n_cols = GetCols();
3299 int n_fullsize = n_rows * n_cols;
3300 if (include_nodata) {
3301 UpdateStringMap(options_, HEADER_INC_NODATA, "TRUE");
3302 } else {
3303 UpdateStringMap(options_, HEADER_INC_NODATA, "FALSE");
3304 }
3305 // 2. Get raster data
3306 T* data_1d = nullptr;
3307 T no_data_value = GetNoDataValue();
3308 int datalength;
3309 string core_name = filename.empty() ? core_name_ : filename;
3310 if (is_2draster) { // 2.1 2D raster data
3311 if (outputdirectly) {
3312 data_1d = raster_2d_[0]; // refers to Initialize2DArray() for why we can do this assignment
3313 datalength = n_lyrs_ * n_cells_; // can be n_lyrs_ * (n_rows*n_cols or valid cells' number)
3314 } else {
3315 datalength = n_lyrs_ * n_fullsize;
3316 Initialize1DArray(datalength, data_1d, no_data_value);
3317 for (int idx = 0; idx < n_cells_; idx++) {
3318 int rowcol_index = pos[idx][0] * n_cols + pos[idx][1];
3319 for (int k = 0; k < n_lyrs_; k++) {
3320 data_1d[n_lyrs_ * rowcol_index + k] = raster_2d_[idx][k];
3321 }
3322 }
3323 }
3324 } else { // 3.2 1D raster data
3325 if (outputdirectly) {
3326 data_1d = raster_;
3327 datalength = n_cells_; // can be n_rows*n_cols or valid cells' number
3328 } else {
3329 datalength = n_fullsize;
3330 Initialize1DArray(datalength, data_1d, no_data_value);
3331 for (int idx = 0; idx < n_cells_; idx++) {
3332 data_1d[pos[idx][0] * n_cols + pos[idx][1]] = raster_[idx];
3333 }
3334 }
3335 }
3336 STRDBL_MAP tmpheader;
3337 CopyHeader(headers_, tmpheader);
3338 if (include_nodata) { UpdateHeader(tmpheader, HEADER_RS_CELLSNUM, n_fullsize); }
3339 bool saved = WriteStreamDataAsGridfs(gfs, core_name, tmpheader,
3340 data_1d, datalength, options_);
3341 if (outputdirectly) {
3342 data_1d = nullptr;
3343 } else {
3344 Release1DArray(data_1d);
3345 }
3346 return saved;
3347}
3348
3349template <typename T, typename MASK_T>
3351 const string& filename /* string() */,
3352 const STRING_MAP& opts /* STRING_MAP() */,
3353 bool include_nodata /* true */,
3354 bool out_origin /* false */,
3355 bool out_combined /* true */,
3356 const map<vint, vector<double> >& recls /* map()*/,
3357 double default_value /* = NODATA_VALUE */) {
3358 if (!ValidateRasterData()) { return false; }
3359 if (nullptr == gfs) { return false; }
3360 if (subset_.empty()) { return false; }
3361 CopyStringMap(opts, options_); // Update metadata
3362 // Added by ljzhu, for compatible with yjwang's code. But, can this key-value be passed by the opts argument?
3363 UpdateStringMapIfNotExist(options_, HEADER_RS_PARAM_ABSTRACTION_TYPE, PARAM_ABSTRACTION_TYPE_PHYSICAL);
3364 if (include_nodata) {
3365 UpdateStringMap(options_, HEADER_INC_NODATA, "TRUE");
3366 }
3367 else {
3368 UpdateStringMap(options_, HEADER_INC_NODATA, "FALSE");
3369 }
3370 if (options_.find(HEADER_RSOUT_DATATYPE) == options_.end()
3371 || StringMatch("Unknown", options_.at(HEADER_RSOUT_DATATYPE))) {
3374 }
3375
3376
3377 int grows = GetRows();
3378 string outnameact = filename.empty() ? core_name_ : filename;
3379 bool out_comb = out_combined;
3380 if (out_comb) {
3381 T* data1d = nullptr;
3382 int sublyrs;
3383 int sublen;
3384 bool flag = PrepareCombSubsetData(&data1d, &sublen, &sublyrs,
3385 out_origin, include_nodata, recls, default_value);
3386 STRDBL_MAP tmpheader;
3387 CopyHeader(headers_, tmpheader);
3388 UpdateHeader(tmpheader, HEADER_RS_LAYERS, sublyrs);
3389 UpdateHeader(tmpheader, HEADER_RS_CELLSNUM, sublen / sublyrs);
3390 flag = flag && WriteStreamDataAsGridfs(gfs, "0_" + outnameact,
3391 tmpheader, data1d, sublen, options_);
3392 Release1DArray(data1d);
3393 return flag;
3394 }
3395 // output each subset
3396 STRDBL_MAP subheader;
3397 bool flag = true;
3398 for (auto it = subset_.begin(); it != subset_.end(); ++it) {
3399 if (!it->second->usable) { continue; }
3400 it->second->GetHeader(GetXllCenter(), GetYllCenter(), grows,
3401 GetCellWidth(), CVT_DBL(no_data_value_), subheader);
3402 T* tmpdata1d = nullptr;
3403 int tmpdatalen;
3404 int tmplyrs;
3405 if (!PrepareSubsetData(it->first, it->second, &tmpdata1d, &tmpdatalen, &tmplyrs,
3406 out_origin, include_nodata, recls, default_value)) {
3407 continue;
3408 }
3409 UpdateHeader(subheader, HEADER_RS_LAYERS, tmplyrs);
3410 UpdateHeader(subheader, HEADER_RS_CELLSNUM, tmpdatalen / tmplyrs);
3411 string tmpfname = itoa(CVT_VINT(it->first)) + "_" + outnameact;
3412 flag = flag && WriteStreamDataAsGridfs(gfs, tmpfname, subheader, tmpdata1d, tmpdatalen, options_);
3413 if (nullptr != tmpdata1d) { Release1DArray(tmpdata1d); }
3414 }
3415 return true;
3416}
3417
3418#endif /* USE_MONGODB */
3419
3420/************* Read functions ***************/
3421template <typename T, typename MASK_T>
3422bool clsRasterData<T, MASK_T>::ReadFromFile(const string& filename, const bool calc_pos /* = false */,
3423 clsRasterData<MASK_T>* mask /* = nullptr */,
3424 const bool use_mask_ext /* = true */,
3425 double default_value /* NODATA_VALUE */,
3426 const STRING_MAP& opts /* = STRING_MAP() */) {
3427 if (!FileExists(filename)) { return false; }
3428 if (!initialized_) { InitializeRasterClass(); }
3429 rs_type_out_ = RasterDataTypeInOptionals(opts);
3430 return ConstructFromSingleFile(filename, calc_pos, mask, use_mask_ext, default_value, opts);
3431}
3432
3433template <typename T, typename MASK_T>
3434bool clsRasterData<T, MASK_T>::ReadFromFiles(vector<string>& filenames, const bool calc_pos /* = false */,
3435 clsRasterData<MASK_T>* mask /* = nullptr */,
3436 const bool use_mask_ext /* = true */,
3437 double default_value /* NODATA_VALUE */,
3438 const STRING_MAP& opts /* = STRING_MAP() */) {
3439 if (!FilesExist(filenames)) { return false; }
3440 n_lyrs_ = CVT_INT(filenames.size());
3441 if (n_lyrs_ == 1) { return ReadFromFile(filenames[0], calc_pos, mask, use_mask_ext, default_value, opts); }
3442 if (!initialized_) { InitializeRasterClass(true); }
3443 rs_type_out_ = RasterDataTypeInOptionals(opts);
3444 // 1. firstly, take the first layer as the main input, to calculate position index or extract by mask.
3445 if (!ConstructFromSingleFile(filenames[0], calc_pos, mask, use_mask_ext, default_value, opts)) {
3446 return false;
3447 }
3448 // 2. change corename and filepath template which format is: `<file dir>/CoreName_%d.<suffix>`
3449 // support "corename.tif or corename_1.tif", "corename_2.tif", and "corename_3.tif", etc.
3450 string::size_type last_underline = core_name_.find_last_of('_');
3451 if (last_underline == string::npos) {
3452 last_underline = core_name_.length();
3453 }
3454 core_name_ = core_name_.substr(0, last_underline);
3455 full_path_ = GetPathFromFullName(filenames[0]) + core_name_ + "_%d." + GetSuffix(filenames[0]);
3456 // 3. initialize raster_2d_ and read the other layers according to position data if stated,
3457 // or just read by row and col
3458 Initialize2DArray(n_cells_, n_lyrs_, raster_2d_, no_data_value_);
3459#pragma omp parallel for
3460 for (int i = 0; i < n_cells_; i++) {
3461 raster_2d_[i][0] = raster_[i];
3462 }
3463 Release1DArray(raster_);
3464 int rows = GetRows();
3465 int cols = GetCols();
3466 // take the first layer as mask, and use_mask_ext is true, and no need to calculate position data
3467 for (int fileidx = 1; fileidx < CVT_INT(filenames.size()); fileidx++) {
3468 STRDBL_MAP tmpheader;
3469 T* tmplyrdata = nullptr;
3470 string curfilename = filenames[fileidx];
3471 if (StringMatch(GetUpper(GetSuffix(curfilename)), string(ASCIIExtension))) {
3472 ReadAscFile(curfilename, tmpheader, tmplyrdata);
3473 }
3474 else {
3475#ifdef USE_GDAL
3476 RasterDataType tmpintype;
3477 string tmpsrs;
3478 ReadRasterFileByGdal(curfilename, tmpheader, tmplyrdata, tmpintype, tmpsrs);
3479#else
3480 StatusMessage("Warning: Only ASC format is supported without GDAL!");
3481 return false;
3482#endif
3483 }
3484 if (nullptr != pos_data_ || nullptr != pos_idx_) {
3485#pragma omp parallel for
3486 for (int i = 0; i < n_cells_; ++i) {
3487 //int tmp_row = pos_data_[i][0]; // row * ncols + col,
3488 //int tmp_col = pos_data_[i][1];
3489 int tmp_row = pos_idx_[i] / cols;
3490 int tmp_col = pos_idx_[i] % cols;
3491 AddOtherLayerRasterData(tmp_row, tmp_col, i, fileidx, tmpheader, tmplyrdata);
3492 }
3493 }
3494 else {
3495#pragma omp parallel for
3496 for (int i = 0; i < rows; i++) {
3497 for (int j = 0; j < cols; j++) {
3498 int cellidx = i * cols + j;
3499 if (FloatEqual(no_data_value_, raster_2d_[cellidx][0])) {
3500 raster_2d_[cellidx][fileidx] = no_data_value_;
3501 continue;
3502 }
3503 AddOtherLayerRasterData(i, j, cellidx, fileidx, tmpheader, tmplyrdata);
3504 }
3505 }
3506 }
3507 Release1DArray(tmplyrdata);
3508 }
3509 if(!is_2draster) is_2draster = true;
3510 UpdateHeader(headers_, HEADER_RS_LAYERS, n_lyrs_); // repair layers count in headers
3511 return true;
3512}
3513
3514#ifdef USE_MONGODB
3515
3516template <typename T, typename MASK_T>
3518 const string& filename,
3519 const bool calc_pos /* = false */,
3520 clsRasterData<MASK_T>* mask /* = nullptr */,
3521 const bool use_mask_ext /* = true */,
3522 double default_value /* NODATA_VALUE */,
3523 const STRING_MAP& opts /* = STRING_MAP() */) {
3524 InitializeReadFunction(filename, calc_pos, mask, use_mask_ext, default_value, opts);
3525 T* dbdata = nullptr;
3526 STRDBL_MAP header_dbl = InitialHeader();
3527 STRING_MAP header_str = InitialStrHeader();
3528 STRING_MAP opts_upd;
3529 CopyStringMap(opts, opts_upd);
3530 UpdateStringMapIfNotExist(opts_upd, HEADER_RS_PARAM_ABSTRACTION_TYPE, PARAM_ABSTRACTION_TYPE_PHYSICAL);
3532 if (opts.empty() || opts.count(HEADER_INC_NODATA) < 1) { // GFS file include nodata by default
3533 UpdateStringMap(opts_upd, HEADER_INC_NODATA, "TRUE");
3534 }
3535 //if(is_conceptual) {
3536 // UpdateStringMap(opts_upd, HEADER_INC_NODATA, "FALSE");
3537 //}
3538 if (!ReadGridFsFile(gfs, filename, dbdata, header_dbl, header_str, opts_upd)) { return false; }
3539
3540 if (headers_.at(HEADER_RS_NROWS) > 0 && headers_.at(HEADER_RS_NCOLS) > 0
3541 && headers_.at(HEADER_RS_LAYERS) > 0 && headers_.at(HEADER_RS_CELLSNUM) > 0) {
3542 // means the raster is preassigned header information, and this function is used for read data only
3543 if (!FloatEqual(headers_.at(HEADER_RS_NROWS), header_dbl.at(HEADER_RS_NROWS))
3544 || !FloatEqual(headers_.at(HEADER_RS_NCOLS), header_dbl.at(HEADER_RS_NCOLS))
3545 || !FloatEqual(headers_.at(HEADER_RS_LAYERS), header_dbl.at(HEADER_RS_LAYERS))
3546 || !FloatEqual(headers_.at(HEADER_RS_CELLSNUM), header_dbl.at(HEADER_RS_CELLSNUM))) {
3547 Release1DArray(dbdata);
3548 return false;
3549 }
3550 } else {
3551 CopyHeader(header_dbl, headers_);
3552 }
3553 CopyStringMap(header_str, options_);
3554
3555 bool include_nodata = !StringMatch(options_.at(HEADER_INC_NODATA), "FALSE");
3556
3557 int n_rows = CVT_INT(headers_.at(HEADER_RS_NROWS));
3558 int n_cols = CVT_INT(headers_.at(HEADER_RS_NCOLS));
3559 n_lyrs_ = CVT_INT(headers_.at(HEADER_RS_LAYERS));
3560 // if (n_rows < 0 || n_cols < 0 || n_lyrs_ < 0) { return false; } // Needless
3561 int fullsize = n_rows * n_cols;
3562 no_data_value_ = static_cast<T>(headers_.at(HEADER_RS_NODATA));
3563 n_cells_ = CVT_INT(headers_.at(HEADER_RS_CELLSNUM));
3564
3565 if (include_nodata && n_cells_ != fullsize) { return false; }
3566 if (n_cells_ != fullsize) { calc_pos_ = true; }
3567
3568 // check the valid values count and determine whether we can read directly.
3569 bool mask_pos_subset = true;
3570 if (nullptr != mask_ && calc_pos_ && use_mask_ext_ && n_cells_ == mask_->GetValidNumber()) {
3571 store_pos_ = false;
3572 mask_->GetRasterPositionData(&n_cells_, &pos_data_);
3573 mask_->GetRasterPositionData(&n_cells_, &pos_idx_);
3574 if (!mask->GetSubset().empty()) {
3575 map<int, SubsetPositions*>& mask_subset = mask_->GetSubset();
3576 for (auto it = mask_subset.begin(); it != mask_subset.end(); ++it) {
3577 SubsetPositions* tmp = new SubsetPositions(it->second, true);
3578 tmp->n_lyrs = n_lyrs_;
3579#ifdef HAS_VARIADIC_TEMPLATES
3580 subset_.emplace(it->first, tmp);
3581#else
3582 subset_.insert(make_pair(it->first, tmp));
3583#endif
3584 }
3585 }
3586 mask_pos_subset = false;
3587 }
3588
3589 if (!include_nodata && (nullptr == pos_data_ || nullptr == pos_idx_)) { return false; }
3590
3591 if (n_lyrs_ == 1) {
3592 is_2draster = false;
3593 if (include_nodata && !mask_pos_subset) {
3594 Initialize1DArray(n_cells_, raster_, no_data_value_);
3595#pragma omp parallel for
3596 for (int i = 0; i < n_cells_; i++) {
3597 //int tmpidx = pos_data_[i][0] * n_cols + pos_data_[i][1];
3598 int tmpidx = pos_idx_[i];
3599 raster_[i] = static_cast<T>(dbdata[tmpidx]);
3600 }
3601 Release1DArray(dbdata);
3602 } else {
3603 raster_ = dbdata;
3604 }
3605 } else {
3606 Initialize2DArray(n_cells_, n_lyrs_, raster_2d_, no_data_value_);
3607#pragma omp parallel for
3608 for (int i = 0; i < n_cells_; i++) {
3609 int tmpidx = i;
3610 if (include_nodata && !mask_pos_subset) {
3611 //tmpidx = pos_data_[i][0] * n_cols + pos_data_[i][1];
3612 tmpidx = pos_idx_[i];
3613 }
3614 for (int j = 0; j < n_lyrs_; j++) {
3615 int idx = tmpidx * n_lyrs_ + j;
3616 raster_2d_[i][j] = dbdata[idx];
3617 }
3618 }
3619 is_2draster = true;
3620 Release1DArray(dbdata);
3621 }
3622 CheckDefaultValue();
3623 if (include_nodata && mask_pos_subset) {
3624 return MaskAndCalculateValidPosition() >= 0;
3625 }
3626 return true;
3627}
3628
3629#endif /* USE_MONGODB */
3630
3631template <typename T, typename MASK_T>
3632void clsRasterData<T, MASK_T>::AddOtherLayerRasterData(const int row, const int col,
3633 const int cellidx, const int lyr,
3634 STRDBL_MAP lyrheader, T* lyrdata) {
3635 int tmpcols = CVT_INT(lyrheader.at(HEADER_RS_NCOLS));
3636 double tmpnodata = lyrheader.at(HEADER_RS_NODATA);
3637 XY_COOR tmp_xy = GetCoordinateByRowCol(row, col);
3638 // get current raster layer's value by XY
3639 ROW_COL tmp_pos = GetPositionByCoordinate(tmp_xy.first, tmp_xy.second, &lyrheader);
3640 if (tmp_pos.first == -1 || tmp_pos.second == -1) {
3641 raster_2d_[cellidx][lyr] = no_data_value_;
3642 } else {
3643 T tmpvalue = lyrdata[tmp_pos.first * tmpcols + tmp_pos.second];
3644 raster_2d_[cellidx][lyr] = FloatEqual(tmpvalue, tmpnodata) ? no_data_value_ : tmpvalue;
3645 }
3646}
3647
3648template <typename T, typename MASK_T>
3650 InitializeRasterClass();
3651 Copy(another);
3652}
3653
3654template <typename T, typename MASK_T>
3656 // Release current data
3657 if (is_2draster && nullptr != raster_2d_ && n_cells_ > 0) {
3658 Release2DArray(raster_2d_);
3659 }
3660 if (!is_2draster && nullptr != raster_) {
3661 Release1DArray(raster_);
3662 }
3663 if (nullptr != pos_data_) {
3664 Release2DArray(pos_data_);
3665 }
3666 if (nullptr != pos_idx_) {
3667 Release1DArray(pos_idx_);
3668 }
3669 if (stats_calculated_) {
3670 ReleaseStatsMap2D();
3671 stats_calculated_ = false;
3672 }
3673 if (!subset_.empty()) {
3674 ReleaseSubset();
3675 }
3676 // Initialize now Raster and copy data
3677 InitializeReadFunction(orgraster->GetFilePath(), orgraster->PositionsCalculated(),
3678 orgraster->GetMask(), orgraster->MaskExtented(),
3679 orgraster->GetDefaultValue(), orgraster->GetOptions());
3680 n_cells_ = orgraster->GetCellNumber();
3681 n_lyrs_ = orgraster->GetLayers();
3682 rs_type_ = orgraster->GetDataType();
3683 rs_type_out_ = orgraster->GetOutDataType();
3684 no_data_value_ = orgraster->GetNoDataValue();
3685 if (orgraster->Is2DRaster()) {
3686 is_2draster = true;
3687 Initialize2DArray(n_cells_, n_lyrs_, raster_2d_,
3688 orgraster->Get2DRasterDataPointer());
3689 } else {
3690 Initialize1DArray(n_cells_, raster_, orgraster->GetRasterDataPointer());
3691 }
3692 if (calc_pos_) {
3693 store_pos_ = true;
3694 Initialize2DArray(n_cells_, 2, pos_data_, orgraster->GetRasterPositionDataPointer());
3695 Initialize1DArray(n_cells_, pos_idx_, orgraster->GetRasterPositionIndexPointer());
3696 }
3697 stats_calculated_ = orgraster->StatisticsCalculated();
3698 if (stats_calculated_) {
3699 if (is_2draster) {
3700 map<string, double *> stats2D = orgraster->GetStatistics2D();
3701 for (auto iter = stats2D.begin(); iter != stats2D.end(); ++iter) {
3702 double* tmpstatvalues = nullptr;
3703 Initialize1DArray(n_lyrs_, tmpstatvalues, iter->second);
3704 stats_2d_.at(iter->first) = tmpstatvalues;
3705 }
3706 } else {
3707 STRDBL_MAP stats = orgraster->GetStatistics();
3708 for (auto iter = stats.begin(); iter != stats.end(); ++iter) {
3709 stats_[iter->first] = iter->second;
3710 }
3711 }
3712 }
3713 CopyHeader(orgraster->GetRasterHeader(), headers_);
3714 // deep copy subset
3715 if (!orgraster->GetSubset().empty()) {
3716 for (auto it = orgraster->GetSubset().begin(); it != orgraster->GetSubset().end(); ++it) {
3717 SubsetPositions* tmp = new SubsetPositions(it->second, true);
3718#ifdef HAS_VARIADIC_TEMPLATES
3719 subset_.emplace(it->first, tmp);
3720#else
3721 subset_.insert(make_pair(it->first, tmp));
3722#endif
3723 }
3724 }
3725}
3726
3727template <typename T, typename MASK_T>
3729#pragma omp parallel for
3730 for (int i = 0; i < n_cells_; i++) {
3731 for (int lyr = 0; lyr < n_lyrs_; lyr++) {
3732 bool flag = is_2draster && nullptr != raster_2d_
3733 ? FloatEqual(raster_2d_[i][lyr], no_data_value_)
3734 : FloatEqual(raster_[i], no_data_value_);
3735 if (!flag) { continue; }
3736 if (is_2draster && nullptr != raster_2d_) {
3737 raster_2d_[i][lyr] = replacedv;
3738 } else if (nullptr != raster_) {
3739 raster_[i] = replacedv;
3740 }
3741 }
3742 }
3743 no_data_value_ = replacedv;
3744 default_value_ = CVT_DBL(replacedv);
3745 UpdateHeader(headers_, HEADER_RS_NODATA, replacedv);
3746}
3747
3748template <typename T, typename MASK_T>
3749void clsRasterData<T, MASK_T>::Reclassify(const map<int, T> reclass_map) {
3750 map<int, T> recls;
3751 for(auto it = reclass_map.begin(); it != reclass_map.end(); ++it) {
3752#ifdef HAS_VARIADIC_TEMPLATES
3753 recls.emplace(it->first, it->second);
3754#else
3755 recls.insert(make_pair(it->first, it->second));
3756#endif
3757 }
3758 for (int i = 0; i < n_cells_; i++) {
3759 for (int lyr = 0; lyr < n_lyrs_; lyr++) {
3760 T curv = is_2draster && nullptr != raster_2d_
3761 ? CVT_INT(raster_2d_[i][lyr])
3762 : CVT_INT(raster_[i]);
3763 if (recls.count(curv) < 0) {
3764#ifdef HAS_VARIADIC_TEMPLATES
3765 recls.emplace(curv, no_data_value_);
3766#else
3767 recls.insert(make_pair(curv, no_data_value_));
3768#endif
3769 }
3770 if (is_2draster && nullptr != raster_2d_) {
3771 raster_2d_[i][lyr] = recls.at(curv);
3772 } else if (nullptr != raster_) {
3773 raster_[i] = recls.at(curv);
3774 }
3775 }
3776 }
3777}
3778
3779/************* Utility functions ***************/
3780
3781template <typename T, typename MASK_T>
3783 return XY_COOR(GetXllCenter() + col * GetCellWidth(),
3784 GetYllCenter() + (GetRows() - row - 1) * GetCellWidth());
3785}
3786
3787template <typename T, typename MASK_T>
3788ROW_COL clsRasterData<T, MASK_T>::GetPositionByCoordinate(const double x, const double y,
3789 STRDBL_MAP* header /* = nullptr */) {
3790 if (nullptr == header) { header = &headers_; }
3791 double xll_center = (*header).at(HEADER_RS_XLL);
3792 double yll_center = (*header).at(HEADER_RS_YLL);
3793 double dx = (*header).at(HEADER_RS_CELLSIZE);
3794 double dy = dx;
3795 int n_rows = CVT_INT((*header).at(HEADER_RS_NROWS));
3796 int n_cols = CVT_INT((*header).at(HEADER_RS_NCOLS));
3797
3798 double x_min = xll_center - dx / 2.;
3799 double x_max = x_min + dx * n_cols;
3800
3801 double y_min = yll_center - dy / 2.;
3802 double y_max = y_min + dy * n_rows;
3803 if ((x > x_max || x < xll_center) || (y > y_max || y < yll_center)) {
3804 return ROW_COL(-1, -1);
3805 }
3806 return ROW_COL(CVT_INT((y_max - y) / dy), CVT_INT((x - x_min) / dx));
3807}
3808
3809template <typename T, typename MASK_T>
3811 vector<T> values; // store 1st layer for both Rater1D and Raster2D
3812 vector<vector<T> > values_2d; // store layer 2~n
3813 vector<int> pos_rows;
3814 vector<int> pos_cols;
3815 int nrows = CVT_INT(headers_.at(HEADER_RS_NROWS));
3816 int ncols = CVT_INT(headers_.at(HEADER_RS_NCOLS));
3817 // get all valid values (i.e., exclude NODATA_VALUE)
3818 for (int i = 0; i < nrows; ++i) {
3819 for (int j = 0; j < ncols; ++j) {
3820 int idx = i * ncols + j;
3821 T tmp_value;
3822 if (is_2draster) {
3823 tmp_value = raster_2d_[idx][0];
3824 } else {
3825 tmp_value = raster_[idx];
3826 }
3827 if (FloatEqual(tmp_value, static_cast<T>(no_data_value_))) continue;
3828 values.emplace_back(tmp_value);
3829 if (is_2draster && n_lyrs_ > 1) {
3830 vector<T> tmpv(n_lyrs_ - 1);
3831 for (int lyr = 1; lyr < n_lyrs_; lyr++) {
3832 tmpv[lyr - 1] = raster_2d_[idx][lyr];
3833 }
3834 values_2d.emplace_back(tmpv);
3835 }
3836 pos_rows.emplace_back(i);
3837 pos_cols.emplace_back(j);
3838 }
3839 }
3840 vector<T>(values).swap(values);
3841 if (is_2draster && n_lyrs_ > 1) {
3842 vector<vector<T> >(values_2d).swap(values_2d);
3843 }
3844 vector<int>(pos_rows).swap(pos_rows);
3845 vector<int>(pos_cols).swap(pos_cols);
3846 // reCreate raster data array
3847 n_cells_ = CVT_INT(values.size());
3848 UpdateHeader(headers_, HEADER_RS_CELLSNUM, n_cells_);
3849 if (is_2draster) {
3850 Release2DArray(raster_2d_);
3851 Initialize2DArray(n_cells_, n_lyrs_, raster_2d_, no_data_value_);
3852 } else {
3853 Release1DArray(raster_);
3854 Initialize1DArray(n_cells_, raster_, no_data_value_);
3855 }
3856 // pos_data_ is nullptr till now.
3857 Initialize2DArray(n_cells_, 2, pos_data_, 0);
3858 Initialize1DArray(n_cells_, pos_idx_, 0);
3859 store_pos_ = true;
3860#pragma omp parallel for
3861 for (int i = 0; i < n_cells_; ++i) {
3862 if (is_2draster) {
3863 raster_2d_[i][0] = values.at(i);
3864 if (n_lyrs_ > 1) {
3865 for (int lyr = 1; lyr < n_lyrs_; lyr++) {
3866 raster_2d_[i][lyr] = values_2d[i][lyr - 1];
3867 }
3868 }
3869 } else {
3870 raster_[i] = values.at(i);
3871 }
3872 pos_data_[i][0] = pos_rows.at(i);
3873 pos_data_[i][1] = pos_cols.at(i);
3874 pos_idx_[i] = pos_rows.at(i) * ncols + pos_cols.at(i);
3875 }
3876 calc_pos_ = true;
3877}
3878
3879template <typename T, typename MASK_T>
3880int clsRasterData<T, MASK_T>::MaskAndCalculateValidPosition() {
3881 int old_fullsize = GetRows() * GetCols();
3882 if (nullptr == mask_) {
3883 if (calc_pos_) {
3884 if (nullptr == pos_data_ || nullptr == pos_idx_) {
3885 CalculateValidPositionsFromGridData();
3886 return 1;
3887 }
3888 return 0;
3889 }
3890 n_cells_ = old_fullsize;
3891 UpdateHeader(headers_, HEADER_RS_CELLSNUM, n_cells_);
3892 // do nothing
3893 return 0;
3894 }
3895 // Use mask data
3896 // 1. Get new values, positions, and subsets (if exist) according to Mask's position data
3897 int mask_ncells;
3898 int** valid_pos = nullptr;
3899 int mask_rows = mask_->GetRows();
3900 int mask_cols = mask_->GetCols();
3901 // Get the position data from mask
3902 if (!mask_->PositionsCalculated()) { mask_->SetCalcPositions(); }
3903 mask_->GetRasterPositionData(&mask_ncells, &valid_pos);
3904 // Masked raster data have the same size with mask's valid positions
3905 vector<T> values(mask_ncells); // store layer 1 data
3906 vector<vector<T> > values_2d(mask_ncells); // store layer 2~n data (excluding the first layer)
3907 vector<int> pos_rows(mask_ncells); // position rows in mask
3908 vector<int> pos_cols(mask_ncells); // position cols in mask
3909 // calculate the intersection extent between mask and the raster data
3910 int max_row = -1;
3911 int min_row = mask_rows;
3912 int max_col = -1;
3913 int min_col = mask_cols;
3914 int masked_count = 0; // position matched count
3915 int matched_count = 0; // valid value matched count
3916 // Get the valid data according to coordinate
3917 for (int i = 0; i < mask_ncells; i++) {
3918 int tmp_row = valid_pos[i][0];
3919 int tmp_col = valid_pos[i][1];
3920 XY_COOR tmp_xy = mask_->GetCoordinateByRowCol(tmp_row, tmp_col);
3921 ROW_COL tmp_pos = GetPositionByCoordinate(tmp_xy.first, tmp_xy.second);
3922 T tmp_value;
3923 if (tmp_pos.first == -1 || tmp_pos.second == -1) {
3924 tmp_value = no_data_value_; // location exceeds the extent of raster data
3925 if (is_2draster && n_lyrs_ > 1) {
3926 vector<T> tmp_values(n_lyrs_ - 1);
3927 for (int lyr = 1; lyr < n_lyrs_; lyr++) { tmp_values[lyr - 1] = no_data_value_; }
3928 values_2d[i] = tmp_values;
3929 }
3930 values[i] = tmp_value;
3931 pos_rows[i] = tmp_row;
3932 pos_cols[i] = tmp_col;
3933 continue;
3934 }
3935 tmp_value = GetValue(tmp_pos.first, tmp_pos.second, 1);
3936 if (FloatEqual(tmp_value, no_data_value_)) {
3937 if (!FloatEqual(default_value_, no_data_value_)) {
3938 tmp_value = static_cast<T>(default_value_);
3939 SetValue(tmp_pos.first, tmp_pos.second, tmp_value);
3940 }
3941 } else { // the intersect extents dependent on the valid raster values
3942 matched_count++;
3943 if (max_row < tmp_row) max_row = tmp_row;
3944 if (min_row > tmp_row) min_row = tmp_row;
3945 if (max_col < tmp_col) max_col = tmp_col;
3946 if (min_col > tmp_col) min_col = tmp_col;
3947 }
3948 if (is_2draster && n_lyrs_ > 1) {
3949 vector<T> tmp_values(n_lyrs_ - 1);
3950 for (int lyr = 1; lyr < n_lyrs_; lyr++) {
3951 tmp_values[lyr - 1] = GetValue(tmp_pos.first, tmp_pos.second, lyr + 1);
3952 if (FloatEqual(tmp_values[lyr - 1], no_data_value_)
3953 && !FloatEqual(default_value_, no_data_value_)) {
3954 tmp_values[lyr - 1] = static_cast<T>(default_value_);
3955 SetValue(tmp_pos.first, tmp_pos.second, tmp_values[lyr - 1], lyr - 1);
3956 }
3957 }
3958 values_2d[i] = tmp_values;
3959 }
3960 values[i] = tmp_value;
3961 pos_rows[i] = tmp_row;
3962 pos_cols[i] = tmp_col;
3963 masked_count++;
3964 }
3965 if (masked_count == 0) { return -1; }
3966 n_cells_ = masked_count;
3967
3968 // Priority Copy header of mask data, and update NoData, SRS, and Layers' count
3969 CopyHeader(mask_->GetRasterHeader(), headers_);
3970 UpdateHeader(headers_, HEADER_RS_NODATA, no_data_value_);
3971 UpdateStrHeader(options_, HEADER_RS_SRS, mask_->GetSrsString());
3972 UpdateHeader(headers_, HEADER_RS_LAYERS, n_lyrs_);
3973
3974 // Priority DEEP Copy subset of mask data
3975 map<int, SubsetPositions*>& mask_subset = mask_->GetSubset();
3976 bool mask_has_subset = !mask_subset.empty(); // if mask data has subsets
3977 if (mask_has_subset) {
3978 ReleaseSubset();
3979 for (auto it = mask_subset.begin(); it != mask_subset.end(); ++it) {
3980 SubsetPositions* tmp = new SubsetPositions(it->second, true);
3981 tmp->n_lyrs = n_lyrs_;
3982#ifdef HAS_VARIADIC_TEMPLATES
3983 subset_.emplace(it->first, tmp);
3984#else
3985 subset_.insert(make_pair(it->first, tmp));
3986#endif
3987 }
3988 }
3989
3990 // Set several flags
3991 // Flag 1: If masked cells' count equals valid positions' count of the mask layer
3992 bool match_exactly = matched_count == mask_ncells;
3993 // Flag2 : Is the valid grid extent same as the mask data, which means the mask
3994 // is within the extent of the raster data
3995 bool within_ext = true;
3996 int new_rows = max_row - min_row + 1;
3997 int new_cols = max_col - min_col + 1;
3998 if (new_rows != mask_rows || new_cols != mask_cols) {
3999 within_ext = false;
4000 } else { // within_ext is true means
4001 // masked cells' count equals valid positions' count of the mask layer
4002 assert(masked_count == mask_ncells);
4003 }
4004 bool upd_header_rowcol = false; // need to update row and col counts in header info?
4005 bool recalc_pos = false; // need to recalculate valid positions? If true, set mask to nullptr!
4006 bool upd_header_valid_num = false; // need to update valid cell count in header info?
4007 bool store_fullsize = false; // need to store full size raster data after masked?
4008 bool recalc_subset = false; // need to recalculate subset? TODO, need more unittests
4009
4010 // Although this if-then may redundancy, I think this still necessary for clearly thinking.
4011 if (within_ext) {
4012 upd_header_rowcol = false;
4013 if (use_mask_ext_) { // use_mask_ext_ has the highest priority
4014 if (calc_pos_) {
4015 recalc_pos = !match_exactly;
4016 upd_header_valid_num = recalc_pos;
4017 store_fullsize = false;
4018 //recalc_subset = mask_has_subset && !match_exactly;
4019 } else { // calc_pos_ = false
4020 recalc_pos = false;
4021 upd_header_valid_num = false;
4022 store_fullsize = false;
4023 calc_pos_ = true; // use_mask_ext_ is priority
4024 //recalc_subset = mask_has_subset && !match_exactly;
4025 }
4026 } else { // use_mask_ext_ is false
4027 if (calc_pos_) {
4028 recalc_pos = !match_exactly;
4029 upd_header_valid_num = recalc_pos;
4030 store_fullsize = false;
4031 //recalc_subset = mask_has_subset && !match_exactly;
4032 } else { // calc_pos_ = false
4033 recalc_pos = false;
4034 upd_header_valid_num = true;
4035 store_fullsize = false;
4036 calc_pos_ = true; // use_mask_ext_ is priority
4037 //recalc_subset = mask_has_subset && !match_exactly;
4038 }
4039 }
4040 } else { // within_ext is false
4041 if (use_mask_ext_) { // use_mask_ext_ has the highest priority
4042 upd_header_rowcol = false;
4043 if (calc_pos_) {
4044 recalc_pos = !match_exactly;
4045 upd_header_valid_num = recalc_pos;
4046 store_fullsize = false;
4047 //recalc_subset = mask_has_subset && !match_exactly;
4048 } else { // calc_pos_ = false
4049 recalc_pos = false;
4050 upd_header_valid_num = false;
4051 store_fullsize = false;
4052 calc_pos_ = true; // use_mask_ext_ is priority
4053 //recalc_subset = mask_has_subset && !match_exactly;
4054 }
4055 } else { // use_mask_ext_ is false
4056 upd_header_rowcol = true;
4057 upd_header_valid_num = true;
4058 if (calc_pos_) {
4059 recalc_pos = true;
4060 store_fullsize = false;
4061 //recalc_subset = mask_has_subset && !match_exactly;
4062 } else { // calc_pos_ = false
4063 recalc_pos = false;
4064 store_fullsize = true;
4065 //recalc_subset = mask_has_subset && !match_exactly;
4066 }
4067 }
4068 }
4069 if (mask_has_subset && !match_exactly && upd_header_rowcol) {
4070 recalc_subset = true;
4071 }
4072 // Update row and col counts in header information firstly.
4073 if (upd_header_rowcol) {
4074 UpdateHeader(headers_, HEADER_RS_NROWS, new_rows);
4075 UpdateHeader(headers_, HEADER_RS_NCOLS, new_cols);
4076 headers_.at(HEADER_RS_XLL) += min_col * mask_->GetCellWidth();
4077 headers_.at(HEADER_RS_YLL) += (mask_rows - max_row - 1) * mask_->GetCellWidth();
4078 headers_.at(HEADER_RS_CELLSIZE) = mask_->GetCellWidth();
4079 }
4080
4081 // ReCalculate valid position
4082 if (recalc_pos) {
4083 // clean redundant values (i.e., NODATA)
4084 auto rit = pos_rows.begin();
4085 auto cit = pos_cols.begin();
4086 auto vit = values.begin();
4087 auto data2dit = values_2d.begin();
4088 for (auto it = values.begin(); it != values.end();) {
4089 size_t idx = distance(vit, it);
4090 int tmpr = pos_rows.at(idx);
4091 int tmpc = pos_cols.at(idx);
4092 if (tmpr > max_row || tmpr < min_row || tmpc > max_col || tmpc < min_col
4093 || FloatEqual(static_cast<T>(*it), no_data_value_)) {
4094 it = values.erase(it);
4095 if (is_2draster && n_lyrs_ > 1) {
4096 values_2d.erase(data2dit + idx);
4097 data2dit = values_2d.begin();
4098 }
4099 pos_cols.erase(cit + idx);
4100 pos_rows.erase(rit + idx);
4101 rit = pos_rows.begin(); // reset the iterators
4102 cit = pos_cols.begin();
4103 vit = values.begin();
4104 } else {
4105 pos_rows[idx] -= min_row; // get new column and row number
4106 pos_cols[idx] -= min_col;
4107 ++it;
4108 }
4109 }
4110 vector<T>(values).swap(values);
4111 if (is_2draster && n_lyrs_ > 1) { vector<vector<T> >(values_2d).swap(values_2d); }
4112 vector<int>(pos_rows).swap(pos_rows);
4113 vector<int>(pos_cols).swap(pos_cols);
4114
4115 store_pos_ = true;
4116 n_cells_ = CVT_INT(values.size());
4117 if (nullptr != pos_data_) { Release2DArray(pos_data_); }
4118 if (nullptr != pos_idx_) { Release1DArray(pos_idx_); }
4119 Initialize2DArray(n_cells_, 2, pos_data_, 0);
4120 Initialize1DArray(n_cells_, pos_idx_, 0);
4121 for (size_t k = 0; k < pos_rows.size(); ++k) {
4122 pos_data_[k][0] = pos_rows.at(k);
4123 pos_data_[k][1] = pos_cols.at(k);
4124 if (upd_header_rowcol) {
4125 pos_idx_[k] = pos_rows.at(k) * new_cols + pos_cols.at(k);
4126 } else {
4127 pos_idx_[k] = pos_rows.at(k) * mask_cols + pos_cols.at(k);
4128 }
4129 }
4130 } else {
4131 if (nullptr != pos_data_) { Release2DArray(pos_data_); }
4132 if (nullptr != pos_idx_) { Release1DArray(pos_idx_); }
4133 if (calc_pos_) {
4134 mask_->GetRasterPositionData(&n_cells_, &pos_data_);
4135 mask_->GetRasterPositionData(&n_cells_, &pos_idx_);
4136 }
4137 store_pos_ = false;
4138 }
4139
4140 // Release the original raster values, and create new
4141 // raster array and positions data array (if necessary)
4142 assert(ValidateRasterData());
4143 int ncols = CVT_INT(headers_.at(HEADER_RS_NCOLS));
4144 int nrows = CVT_INT(headers_.at(HEADER_RS_NROWS));
4145 if (store_fullsize) {
4146 n_cells_ = ncols * nrows;
4147 }
4148 bool release_origin = true;
4149 if (store_fullsize && old_fullsize == n_cells_) {
4150 release_origin = false;
4151 }
4152 if (release_origin) {
4153 if (is_2draster && nullptr != raster_2d_) { // multiple layers
4154 Release2DArray(raster_2d_);
4155 Initialize2DArray(n_cells_, n_lyrs_, raster_2d_, no_data_value_);
4156 } else { // single layer
4157 Release1DArray(raster_);
4158 Initialize1DArray(n_cells_, raster_, no_data_value_);
4159 }
4160 // Loop the masked raster values
4161 size_t synthesis_idx = 0;
4162 for (size_t k = 0; k < pos_rows.size(); ++k) {
4163 synthesis_idx = k;
4164 int tmpr = pos_rows.at(k);
4165 int tmpc = pos_cols.at(k);
4166 if (store_fullsize) {
4167 if (tmpr > max_row || tmpr < min_row || tmpc > max_col || tmpc < min_col) {
4168 continue;
4169 }
4170 synthesis_idx = (tmpr - min_row) * ncols + tmpc - min_col;
4171 if (synthesis_idx > n_cells_ - 1) {
4172 continue; // error may occurred!
4173 }
4174 }
4175 if (is_2draster) { // multiple layers
4176 raster_2d_[synthesis_idx][0] = values.at(k);
4177 if (n_lyrs_ > 1) {
4178 for (int lyr = 1; lyr < n_lyrs_; lyr++) {
4179 raster_2d_[synthesis_idx][lyr] = values_2d[k][lyr - 1];
4180 }
4181 }
4182 } else { // single layer
4183 raster_[synthesis_idx] = values.at(k);
4184 }
4185 }
4186 }
4187
4188 if (upd_header_valid_num) {
4189 UpdateHeader(headers_, HEADER_RS_CELLSNUM, n_cells_);
4190 }
4191
4192 if (mask_has_subset) { // check former assigned mask's subset
4193 for (auto it = subset_.begin(); it != subset_.end();) {
4194 int count = 0;
4195 int srow = nrows;
4196 int erow = -1;
4197 int scol = ncols;
4198 int ecol = -1;
4199 vector<int> globalpos;
4200 for (int i = 0; i < it->second->n_cells; i++) {
4201 int gi = it->second->global_[i];
4202 int tmprow = valid_pos[gi][0];
4203 int tmpcol = valid_pos[gi][1];
4204 XY_COOR tmpxy = mask_->GetCoordinateByRowCol(tmprow, tmpcol);
4205 if (GetPosition(tmpxy.first, tmpxy.second) < 0) { continue; }
4206 ROW_COL tmppos = GetPositionByCoordinate(tmpxy.first, tmpxy.second);
4207 if (IsNoData(tmppos.first, tmppos.second)) { continue; }
4208 if (tmppos.first < srow) { srow = tmppos.first; }
4209 if (tmppos.first > erow) { erow = tmppos.first; }
4210 if (tmppos.second < scol) { scol = tmppos.second; }
4211 if (tmppos.second > ecol) { ecol = tmppos.second; }
4212 globalpos.emplace_back(GetPosition(tmppos.first, tmppos.second));
4213 count++;
4214 }
4215 if (count == 0) {
4216 delete it->second;
4217 subset_.erase(it++);
4218 continue;
4219 }
4220 if (recalc_subset) {
4221 it->second->g_srow = srow;
4222 it->second->g_erow = erow;
4223 it->second->g_scol = scol;
4224 it->second->g_ecol = ecol;
4225 it->second->n_cells = count;
4226 it->second->n_lyrs = n_lyrs_;
4227 int local_ncols = ecol - scol + 1;
4228 if (it->second->alloc_) {
4229 Release1DArray(it->second->global_);
4230 Release2DArray(it->second->local_pos_);
4231 Release1DArray(it->second->local_posidx_);
4232 }
4233 it->second->global_ = nullptr; // not affect mask's subset
4234 it->second->local_pos_ = nullptr;
4235 it->second->local_posidx_ = nullptr;
4236 it->second->alloc_ = true;
4237 Initialize1DArray(count, it->second->global_, -1);
4238 Initialize2DArray(count, 2, it->second->local_pos_, -1);
4239 Initialize1DArray(count, it->second->local_posidx_, -1);
4240 for (int ii = 0; ii < count; ii++) {
4241 it->second->global_[ii] = globalpos[ii];
4242 int local_row = -1;
4243 int local_col = -1;
4244 if (nullptr == pos_data_ || nullptr == pos_idx_) {
4245 local_row = globalpos[ii] / ncols - it->second->g_srow;
4246 local_col = globalpos[ii] % ncols - it->second->g_scol;
4247 } else {
4248 //local_row = pos_data_[globalpos[ii]][0] - it->second->g_srow;
4249 //local_col = pos_data_[globalpos[ii]][1] - it->second->g_scol;
4250 local_row = pos_idx_[globalpos[ii]] / ncols - it->second->g_srow;
4251 local_col = pos_idx_[globalpos[ii]] % ncols - it->second->g_scol;
4252 }
4253 it->second->local_pos_[ii][0] = local_row;
4254 it->second->local_pos_[ii][1] = local_col;
4255 it->second->local_posidx_[ii] = local_row * local_ncols + local_col;
4256 }
4257 it->second->data_ = nullptr;
4258 it->second->data2d_ = nullptr;
4259 it->second->alloc_ = true;
4260 }
4261 ++it;
4262 }
4263 }
4264 if (store_fullsize || recalc_pos) {
4265 mask_ = nullptr;
4266 }
4267 return 2; // all situations that use mask data
4268}
4269
4270} // namespace data_raster
4271} // namespace ccgl
4272#endif /* CCGL_DATA_RASTER_H */
Basic definitions.
#define CVT_INT(param)
A reference to the postfix of executable file for RELWITHDEBINFO mode.
Definition: basic.h:325
#define MISSINGFLOAT
Missing float value.
Definition: basic.h:250
#define CVT_DBL(param)
Convert to double double
Definition: basic.h:331
#define CVT_VINT(param)
Convert to 8-byte (64-bit) signed integer vint
Definition: basic.h:340
#define NODATA_VALUE
Global utility definitions.
Definition: basic.h:245
Base class for classes that cannot be copied.
Definition: basic.h:385
Subset positions of raster data.
Definition: data_raster.hpp:1085
int * global_
global position index
Definition: data_raster.hpp:1167
int * local_posidx_
local position index
Definition: data_raster.hpp:1166
double ** data2d_
valid 2d data array
Definition: data_raster.hpp:1169
int n_cells
valid cell count
Definition: data_raster.hpp:1158
int g_erow
end row in global data
Definition: data_raster.hpp:1161
bool alloc_
local_pos_ and global_ are allocated?
Definition: data_raster.hpp:1164
bool usable
flag for usable subset data
Definition: data_raster.hpp:1157
int g_ecol
end col in global data
Definition: data_raster.hpp:1163
int ** local_pos_
local position data
Definition: data_raster.hpp:1165
double * data_
valid data array
Definition: data_raster.hpp:1168
int g_srow
start row in global data
Definition: data_raster.hpp:1160
int n_lyrs
layer count
Definition: data_raster.hpp:1159
int g_scol
start col in global data
Definition: data_raster.hpp:1162
Raster data (1D and 2D) I/O class Support I/O among ASCII file, TIFF, and MongoDB database.
Definition: data_raster.hpp:1178
const STRDBL_MAP & GetRasterHeader() const
Get raster header information.
Definition: data_raster.hpp:1687
clsRasterData(bool is_2d=false)
Constructor an empty clsRasterData instance for 1D or 2D raster.
Definition: data_raster.hpp:2073
bool OutputFileByGdal(const string &filename)
Write 1D or 2D raster data into TIFF file by GDAL.
Definition: data_raster.hpp:3213
RasterDataType GetOutDataType() const
Get data type of source.
Definition: data_raster.hpp:1655
XY_COOR GetCoordinateByRowCol(int row, int col)
Calculate XY coordinates by given row and col number.
Definition: data_raster.hpp:3782
bool SetUseMaskExt()
Set the flag of use_mask_ext_ to true and recalculate positions and stored raster data if necessary.
Definition: data_raster.hpp:2839
~clsRasterData()
Destructor.
Definition: data_raster.hpp:2326
bool ReleaseSubset()
Release subsets.
Definition: data_raster.hpp:2406
bool MaskExtented() const
position data is allocated or a pointer
Definition: data_raster.hpp:1757
string GetFullFileName() const
Get full filename.
Definition: data_raster.hpp:1809
double GetCellWidth() const
Get row number.
Definition: data_raster.hpp:1629
int ** GetRasterPositionDataPointer() const
Get pointer of raster 1D data.
Definition: data_raster.hpp:1711
bool SetPositions(int len, int **pdata)
Set valid positions data, without mask raster layer.
Definition: data_raster.hpp:2815
int * GetRasterPositionIndexPointer() const
Get pointer of position data.
Definition: data_raster.hpp:1712
double GetYllCenter() const
Get Y coordinate of left lower center of raster data.
Definition: data_raster.hpp:1643
string GetCoreName() const
Get core name.
Definition: data_raster.hpp:1699
bool OutputAscFile(const string &filename)
Write 1D or 2D raster data into ASC file(s)
Definition: data_raster.hpp:3136
bool ReadFromMongoDB(MongoGridFs *gfs, const string &filename, bool calc_pos=false, clsRasterData< MASK_T > *mask=nullptr, bool use_mask_ext=true, double default_value=NODATA_VALUE, const STRING_MAP &opts=STRING_MAP())
Read raster data from MongoDB.
Definition: data_raster.hpp:3517
int GetDataLength() const
Get the first dimension size.
Definition: data_raster.hpp:1626
bool GetRasterData(int *n_cells, T **data)
Get raster data, include valid cell number and data.
Definition: data_raster.hpp:2614
bool StatisticsCalculated() const
Use mask extent or not.
Definition: data_raster.hpp:1758
string GetFilePath() const
Get full path name.
Definition: data_raster.hpp:1696
double GetRange(const int lyr=1)
Get the range of raster data.
Definition: data_raster.hpp:1578
string GetSrsString()
Get the spatial reference (char*)
Definition: data_raster.hpp:2681
void SetDataType(RasterDataType const type)
Set raster data type.
Definition: data_raster.hpp:1481
static clsRasterData< T, MASK_T > * Init(const string &filename, bool calc_pos=false, clsRasterData< MASK_T > *mask=nullptr, bool use_mask_ext=true, double default_value=NODATA_VALUE, const STRING_MAP &opts=STRING_MAP())
Validation check before the construction of clsRasterData, i.e., the raster file is existed and the d...
Definition: data_raster.hpp:2149
double GetStatistics(string sindex, int lyr=1)
Get basic statistics value Mean, Max, Min, STD, Range, etc.
Definition: data_raster.hpp:2484
void GetStd(int *lyrnum, double **values)
Get the standard derivation of 2D raster data.
Definition: data_raster.hpp:1603
clsRasterData< MASK_T > * GetMask() const
Get mask data pointer.
Definition: data_raster.hpp:1812
void SetDefaultValue(const double defaultv)
Set default value.
Definition: data_raster.hpp:1493
double GetMinimum(const int lyr=1)
Get the minimum of raster data.
Definition: data_raster.hpp:1565
bool ValidateRowCol(const int row, const int col)
Validate the input row and col.
Definition: data_raster.hpp:1775
int GetPosition(int row, int col)
Get default value.
Definition: data_raster.hpp:2541
bool OutputToMongoDB(MongoGridFs *gfs, const string &filename=string(), const STRING_MAP &opts=STRING_MAP(), bool include_nodata=true, bool out_origin=true)
Write the whole raster data to MongoDB, name format: corename__LyrNum.
Definition: data_raster.hpp:3268
const STRDBL_MAP & GetStatistics() const
Get raster statistics information.
Definition: data_raster.hpp:1690
void SetValue(int row, int col, T value, int lyr=1)
Set value to the given position and layer.
Definition: data_raster.hpp:2785
T ** Get2DRasterDataPointer() const
Get pointer of position data.
Definition: data_raster.hpp:1713
void ReleaseStatsMap2D()
Release statistics map of 2D raster data.
Definition: data_raster.hpp:2467
bool ValidateLayer(const int lyr)
Validate the input layer number.
Definition: data_raster.hpp:1787
ROW_COL GetPositionByCoordinate(double x, double y, STRDBL_MAP *header=nullptr)
Calculate position by given coordinate.
Definition: data_raster.hpp:3788
bool PositionsAllocated() const
Calculate positions or not.
Definition: data_raster.hpp:1756
void Reclassify(map< int, T > reclass_map)
classify raster
Definition: data_raster.hpp:3749
const map< string, double * > & GetStatistics2D() const
Get raster statistics information of 2D raster.
Definition: data_raster.hpp:1693
double GetAverage(const int lyr=1)
Get the average of raster data.
Definition: data_raster.hpp:1553
T GetValueByIndex(int cell_index, int lyr=1)
Get options.
Definition: data_raster.hpp:2695
void GetAverage(int *lyrnum, double **values)
Get the average of 2D raster data.
Definition: data_raster.hpp:1585
void SetDataType(const string &strtype)
Set raster data type.
Definition: data_raster.hpp:1484
bool RebuildSubSet(map< int, int > groups=map< int, int >())
Re-build subsets by given groups (cell value->group ID) or by discrete values.
Definition: data_raster.hpp:2418
const char * GetSrs()
Get pointer of raster 2D data.
Definition: data_raster.hpp:2676
int GetValidNumber(const int lyr=1)
Get the non-NoDATA cells number of the given raster layer data.
Definition: data_raster.hpp:1622
void UpdateStatistics()
Force to update basic statistics values Mean, Max, Min, STD, Range, etc.
Definition: data_raster.hpp:2477
int GetRows() const
Get column number.
Definition: data_raster.hpp:1628
void ReplaceNoData(T replacedv)
Replace NoData value by the given value.
Definition: data_raster.hpp:3728
void SetOutDataType(RasterDataType type)
Set output raster data type.
Definition: data_raster.hpp:2425
T GetNoDataValue() const
Get data type of output.
Definition: data_raster.hpp:1656
map< int, SubsetPositions * > & GetSubset()
Get subset.
Definition: data_raster.hpp:1673
RasterDataType GetDataType() const
Get layer number.
Definition: data_raster.hpp:1654
bool ReadFromFiles(vector< string > &filenames, bool calc_pos=false, clsRasterData< MASK_T > *mask=nullptr, bool use_mask_ext=true, double default_value=NODATA_VALUE, const STRING_MAP &opts=STRING_MAP())
Read raster data from two or more files, mask data is optional.
Definition: data_raster.hpp:3434
bool OutputToFile(const string &filename, bool out_origin=true)
Write the whole raster to raster file, if 2D raster, output name will be corename_LyrNum.
Definition: data_raster.hpp:2850
void GetValidNumber(int *lyrnum, double **values)
Get the non-NoDATA number of 2D raster data.
Definition: data_raster.hpp:1615
double GetMaximum(const int lyr=1)
Get the maximum of raster data.
Definition: data_raster.hpp:1559
void Copy(clsRasterData< T, MASK_T > *orgraster)
Copy clsRasterData object.
Definition: data_raster.hpp:3655
bool OutputSubsetToFile(bool out_origin=false, bool out_combined=true, const string &outname=string(), const map< vint, vector< double > > &recls=map< vint, vector< double > >(), double default_value=NODATA_VALUE)
Write one or more raster's subset to files, default name format: corename_SubID_LyrNum.
Definition: data_raster.hpp:2871
bool ValidateIndex(const int idx)
Validate the input index.
Definition: data_raster.hpp:1800
double GetXllCenter() const
Get cell size.
Definition: data_raster.hpp:1632
void GetRange(int *lyrnum, double **values)
Get the range of 2D raster data.
Definition: data_raster.hpp:1609
double GetDefaultValue() const
Get NoDATA value.
Definition: data_raster.hpp:1657
const STRING_MAP & GetOptions() const
Get option by key, including the spatial reference by "SRS".
Definition: data_raster.hpp:1717
bool Get2DRasterData(int *n_cells, int *n_lyrs, T ***data)
Get 2D raster data, include valid cell number of each layer, layer number, and data.
Definition: data_raster.hpp:2626
T GetValue(int row, int col, int lyr=1)
Get raster data via row and col The default lyr is 1, which means the 1D raster data,...
Definition: data_raster.hpp:2725
string GetOption(const char *key)
Get the spatial reference (string)
Definition: data_raster.hpp:2686
void CalculateStatistics()
Calculate basic statistics values in one time Mean, Max, Min, STD, Range, etc.
Definition: data_raster.hpp:2438
bool IsNoData(const int row, const int col, const int lyr=1)
Check if the raster data is NoData via row and col The default lyr is 1, which means the 1D raster da...
Definition: data_raster.hpp:1750
void GetMinimum(int *lyrnum, double **values)
Get the minimum of 2D raster data.
Definition: data_raster.hpp:1597
bool PositionsCalculated() const
Is 2D raster data?
Definition: data_raster.hpp:1755
bool Initialized() const
Basic statistics calculated?
Definition: data_raster.hpp:1759
bool ValidateRasterData()
Instance of clsRasterData initialized?
Definition: data_raster.hpp:1764
bool OutputSubsetToMongoDB(MongoGridFs *gfs, const string &filename=string(), const STRING_MAP &opts=STRING_MAP(), bool include_nodata=true, bool out_origin=false, bool out_combined=true, const map< vint, vector< double > > &recls=map< vint, vector< double > >(), double default_value=NODATA_VALUE)
Write one or more raster's subset to MongoDB,.
Definition: data_raster.hpp:3350
double GetStd(const int lyr=1)
Get the stand derivation of raster data.
Definition: data_raster.hpp:1571
void SetCoreName(const string &name)
Set new core file name.
Definition: data_raster.hpp:1455
void GetMaximum(int *lyrnum, double **values)
Get the maximum of 2D raster data.
Definition: data_raster.hpp:1591
bool SetCalcPositions()
Set the flag of calc_pos_ to true and recalculate positions and stored raster data if necessary.
Definition: data_raster.hpp:2804
bool BuildSubSet(map< int, int > groups=map< int, int >())
Build subsets by given groups (cell value->group ID) or by discrete values (default)
Definition: data_raster.hpp:2337
bool ReadFromFile(const string &filename, bool calc_pos=false, clsRasterData< MASK_T > *mask=nullptr, bool use_mask_ext=true, double default_value=NODATA_VALUE, const STRING_MAP &opts=STRING_MAP())
Read raster data from file, mask data is optional.
Definition: data_raster.hpp:3422
void GetRasterPositionData(int *datalength, int ***positiondata)
Get position index data and the data length.
Definition: data_raster.hpp:2640
A simple wrapper of the class of MongoDB database mongoc_gridfs_t.
Definition: db_mongoc.h:141
bool RemoveFile(string const &gfilename, mongoc_gridfs_t *gfs=NULL, STRING_MAP opts=STRING_MAP())
Remove GridFS all matching files and their data chunks.
Definition: db_mongoc.cpp:338
bool GetStreamData(string const &gfilename, char *&databuf, vint &datalength, mongoc_gridfs_t *gfs=NULL, const STRING_MAP *opts=nullptr)
Get stream data of a given GridFS file name.
Definition: db_mongoc.cpp:416
bool WriteStreamData(const string &gfilename, char *&buf, vint length, const bson_t *p, mongoc_gridfs_t *gfs=NULL)
Write stream data to a GridFS file.
Definition: db_mongoc.cpp:447
bson_t * GetFileMetadata(string const &gfilename, mongoc_gridfs_t *gfs=NULL, STRING_MAP opts=STRING_MAP())
Get metadata of a given GridFS file name, remember to destory bson_t after use.
Definition: db_mongoc.cpp:396
#define CONST_CHARS
const string
Definition: data_raster.hpp:64
Simple wrappers of the API of MongoDB C driver mongo-c-driver, see MongoDB C Driver for more informat...
Simple wrappers of GDAL API that compatible with all versions.
void UpdateStrHeader(STRING_MAP &strheader, const string &key, const string &val)
Update header information in string.
Definition: data_raster.cpp:155
CONST_CHARS HEADER_INC_NODATA
Desired output data type of raster.
Definition: data_raster.hpp:94
bool WriteStreamDataAsGridfs(MongoGridFs *gfs, const string &filename, STRDBL_MAP &header, T *values, const int datalength, const STRING_MAP &opts=STRING_MAP())
Write array data (both valid and full-sized raster data) as GridFS file.
Definition: data_raster.hpp:988
bool ReadRasterFileByGdal(const string &filename, STRDBL_MAP &header, T *&values, RasterDataType &in_type, string &srs)
Read single raster file by GDAL.
Definition: data_raster.hpp:354
bool WriteSingleGeotiff(const string &filename, const STRDBL_MAP &header, const STRING_MAP &opts, T *values)
Write single geotiff file If the file exists, delete it first.
Definition: data_raster.hpp:581
STRING_MAP InitialStrHeader()
Initialize header information in string.
Definition: data_raster.cpp:137
bool WriteRasterToFile(const string &filename, const STRDBL_MAP &header, const STRING_MAP &opts, T *values)
Write single raster file, if the file exists, delete it first.
Definition: data_raster.hpp:838
CONST_CHARS STATS_RS_MEAN
Valid cell number.
Definition: data_raster.hpp:97
CONST_CHARS HEADER_RS_CELLSNUM
Layers number.
Definition: data_raster.hpp:90
CONST_CHARS HEADER_RS_PARAM_ABSTRACTION_TYPE
GeoTIFF extension.
Definition: data_raster.hpp:105
CONST_CHARS HEADER_RS_YLL
X coordinate value of left low center.
Definition: data_raster.hpp:83
CONST_CHARS STATS_RS_VALIDNUM
Mask layer's name if only store valid values.
Definition: data_raster.hpp:96
RasterDataType
Coordinate pair.
Definition: data_raster.hpp:116
@ RDT_UInt16
GDT_UInt16.
Definition: data_raster.hpp:120
@ RDT_Int32
GDT_Int32.
Definition: data_raster.hpp:123
@ RDT_UInt32
GDT_UInt32.
Definition: data_raster.hpp:122
@ RDT_Int64
GDT_Int64, GDAL>=3.5.
Definition: data_raster.hpp:125
@ RDT_Float
GDT_Float32.
Definition: data_raster.hpp:126
@ RDT_UInt64
GDT_UInt64, GDAL>=3.5.
Definition: data_raster.hpp:124
@ RDT_Double
GDT_Float64.
Definition: data_raster.hpp:127
@ RDT_Int16
GDT_Int16.
Definition: data_raster.hpp:121
@ RDT_Unknown
GDT_Unknown.
Definition: data_raster.hpp:117
@ RDT_Int8
GDT_Int8, GDAL>=3.7.
Definition: data_raster.hpp:119
@ RDT_UInt8
GDT_Byte.
Definition: data_raster.hpp:118
void CopyHeader(const STRDBL_MAP &refers, STRDBL_MAP &dst)
Copy header information from one to another.
Definition: data_raster.cpp:121
bool WriteSingleAsc(const string &filename, const STRDBL_MAP &header, T *values)
Write raster data as a single ASC file.
Definition: data_raster.hpp:320
CONST_CHARS GTiffExtension
ASCII extension.
Definition: data_raster.hpp:103
STRDBL_MAP InitialHeader()
Initialize header information in double.
Definition: data_raster.cpp:105
double DefaultNoDataByType(const RasterDataType type)
Default NoData value by data type.
Definition: data_raster.cpp:62
CONST_CHARS HEADER_RS_DATATYPE
SRS.
Definition: data_raster.hpp:92
CONST_CHARS STATS_RS_MIN
Mean value.
Definition: data_raster.hpp:98
RasterDataType TypeToRasterDataType(const std::type_info &t)
Convert C++ data type to RasterDataType.
Definition: data_raster.cpp:48
CONST_CHARS HEADER_RS_LAYERS
Cell size (length)
Definition: data_raster.hpp:89
CONST_CHARS STATS_RS_RANGE
Standard derivation value.
Definition: data_raster.hpp:101
CONST_CHARS ASCIIExtension
Range value.
Definition: data_raster.hpp:102
CONST_CHARS HEADER_RSOUT_DATATYPE
Data type of original raster.
Definition: data_raster.hpp:93
CONST_CHARS HEADER_RS_SRS
Number of the first layer's valid cells.
Definition: data_raster.hpp:91
void InitialStatsMap(STRDBL_MAP &stats, map< string, double * > &stats2d)
Initialize statistics values for 1D and 2D raster data.
Definition: data_raster.cpp:173
bool WriteAscHeaders(const string &filename, const STRDBL_MAP &header)
Write raster header information into a ASC file.
Definition: data_raster.cpp:189
CONST_CHARS STATS_RS_STD
Maximum value.
Definition: data_raster.hpp:100
CONST_CHARS HEADER_RS_CELLSIZE
Column number.
Definition: data_raster.hpp:88
RasterDataType StringToRasterDataType(const string &stype)
Convert string to RasterDataType.
Definition: data_raster.cpp:32
string RasterDataTypeToString(const int type)
Common functions independent to clsRasterData.
Definition: data_raster.cpp:16
CONST_CHARS HEADER_RS_YLLCOR
X coordinate value of left low center.
Definition: data_raster.hpp:85
CONST_CHARS PARAM_ABSTRACTION_TYPE_CONEPTUAL
spatial parameter type, physical or conceptual
Definition: data_raster.hpp:106
CONST_CHARS HEADER_MASK_NAME
Include nodata ("TRUE") or not ("FALSE"), for DB only.
Definition: data_raster.hpp:95
CONST_CHARS STATS_RS_MAX
Minimum value.
Definition: data_raster.hpp:99
std::pair< double, double > XY_COOR
Row and Col pair.
Definition: data_raster.hpp:111
void UpdateHeader(STRDBL_MAP &header, const string &key, T val)
Update value in header information.
Definition: data_raster.hpp:172
bool ReadGridFsFile(MongoGridFs *gfs, const string &filename, T *&data, STRDBL_MAP &header, STRING_MAP &header_str, const STRING_MAP &opts)
Read GridFs file from MongoDB.
Definition: data_raster.hpp:869
CONST_CHARS HEADER_RS_NCOLS
Rows number.
Definition: data_raster.hpp:87
CONST_CHARS HEADER_RS_XLLCOR
Y coordinate value of left low center.
Definition: data_raster.hpp:84
CONST_CHARS HEADER_RS_NROWS
Y coordinate value of left low center.
Definition: data_raster.hpp:86
bool ReadAscFile(const string &filename, STRDBL_MAP &header, T *&values)
Read raster data from ASC file, the simply usage.
Definition: data_raster.hpp:210
RasterDataType RasterDataTypeInOptionals(const STRING_MAP &opts)
Get output raster data type from optional inputs.
Definition: data_raster.cpp:166
CONST_CHARS HEADER_RS_XLL
NoData value.
Definition: data_raster.hpp:82
void AppendStringOptionsToBson(bson_t *bson_opts, const STRING_MAP &opts, const string &prefix)
Append options to bson_t
Definition: db_mongoc.cpp:486
string GetStringFromBsonIterator(bson_iter_t *iter)
Get String from the iterator (bson_iter_t) of bson_t
Definition: db_mongoc.cpp:508
bool GetNumericFromBsonIterator(bson_iter_t *iter, T &numericvalue)
Get numeric value from the iterator (bson_iter_t) of bson_taccording to a given key.
Definition: db_mongoc.h:192
void Release1DArray(T *&data)
Release DT_Array1D data.
Definition: utils_array.h:476
bool Initialize1DArray(int row, T *&data, INI_T init_value)
Initialize DT_Array1D data.
Definition: utils_array.h:344
bool Initialize2DArray(int row, int col, T **&data, INI_T init_value)
Initialize DT_Array2D data.
Definition: utils_array.h:394
void Release2DArray(T **&data)
Release DT_Array2D data.
Definition: utils_array.h:484
bool FileExists(string const &filename)
Return a flag indicating if the given file exists.
Definition: utils_filesystem.cpp:23
string GetSuffix(string const &full_filename)
Return the suffix of a given file's path without dot, e.g., "tif", "asc".
Definition: utils_filesystem.cpp:312
string GetAbsolutePath(string const &full_filename)
Return the absolute file path from a given file path.
Definition: utils_filesystem.cpp:284
bool PathExists(string const &fullpath)
Return a flag indicating if the given path (directory or file) exists.
Definition: utils_filesystem.cpp:43
string AppendCoreFileName(string const &full_filename, string const &endstr, char deli)
Append a given string to the core filename.
Definition: utils_filesystem.cpp:329
string GetPathFromFullName(string const &full_filename)
Get Path From full file path string.
Definition: utils_filesystem.cpp:351
string ReplaceSuffix(string const &full_filename, string const &new_suffix)
Replace the suffix by a given suffix.
Definition: utils_filesystem.cpp:321
string GetCoreFileName(string const &full_filename)
Return the file name from a given file's path.
Definition: utils_filesystem.cpp:297
string PrefixCoreFileName(string const &full_filename, string const &prestr, char deli)
Add a prefix to the core filename.
Definition: utils_filesystem.cpp:340
bool FilesExist(vector< string > &filenames)
Return a flag indicating if given files exist.
Definition: utils_filesystem.cpp:35
bool MakeDirectory(const string &dirpath)
Make directory if not exists.
Definition: utils_filesystem.cpp:140
bool LoadPlainTextFile(const string &filepath, vector< string > &content_strs)
Load short plain text file as string vector, ignore comments begin with '#' and empty lines.
Definition: utils_filesystem.cpp:373
void BasicStatistics(const T *values, int num, double **derivedvalues, T exclude=static_cast< T >(NODATA_VALUE))
calculate basic statistics at one time_funcs
Definition: utils_math.h:416
bool FloatEqual(T1 v1, T2 v2)
Whether v1 is equal to v2.
Definition: utils_math.h:145
string ValueToString(const T &val)
Convert value to string.
Definition: utils_string.h:88
void CopyStringMap(const STRING_MAP &in_opts, STRING_MAP &out_opts)
Copy string map.
Definition: utils_string.cpp:78
void UpdateStringMap(STRING_MAP &opts, const string &key, const string &value)
Add or modify element in a string map.
Definition: utils_string.cpp:94
string itoa(vint number)
Convert a signed integer to a string.
Definition: utils_string.cpp:162
vint IsInt(const string &num_str, bool &success)
Check if a string is an signed integer, if ture, return the converted integer.
Definition: utils_string.cpp:260
string GetUpper(const string &str)
Get Uppercase of given string.
Definition: utils_string.cpp:16
vector< string > SplitString(const string &item, const char delimiter)
Splits the given string based on the given delimiter.
Definition: utils_string.cpp:35
bool StringMatch(const char *a, const char *b)
Match char ignore cases.
Definition: utils_string.cpp:61
double IsDouble(const string &num_str, bool &success)
Convert an Ansi string to 64-bits floating point number.
Definition: utils_string.cpp:316
Common Cross-platform Geographic Library (CCGL)
Definition: basic.cpp:17
std::map< string, string > STRING_MAP
Map of string key and string value.
Definition: basic.h:349
void StatusMessage(const char *msg)
Print status messages for Debug.
Definition: basic.cpp:136
std::map< string, double > STRDBL_MAP
Map of string key and double value.
Definition: basic.h:352
void SleepMs(const int millisecs)
Sleep milliseconds.
Definition: basic.h:491
Template functions to initialize and release arrays.
File system related functions in CCGL.
Useful math equations in CCGL.
Handling string related issues in CCGL.