diff --git a/R/RcppExports.R b/R/RcppExports.R index ef94095..e29b99c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -93,8 +93,8 @@ gc_create_view <- function(v) { .Call('_gdalcubes_gc_create_view', PACKAGE = 'gdalcubes', v) } -gc_create_image_collection_cube <- function(pin, chunk_sizes, mask, v = NULL) { - .Call('_gdalcubes_gc_create_image_collection_cube', PACKAGE = 'gdalcubes', pin, chunk_sizes, mask, v) +gc_create_image_collection_cube <- function(pin, chunk_sizes, mask, strict = TRUE, v = NULL) { + .Call('_gdalcubes_gc_create_image_collection_cube', PACKAGE = 'gdalcubes', pin, chunk_sizes, mask, strict, v) } gc_create_ncdf_cube <- function(path, chunk_sizes, auto_unpack) { @@ -209,8 +209,8 @@ gc_create_stream_cube <- function(pin, cmd) { .Call('_gdalcubes_gc_create_stream_cube', PACKAGE = 'gdalcubes', pin, cmd) } -gc_create_simple_cube <- function(files, datetime_values, bands, band_names, dx, dy, chunk_sizes) { - .Call('_gdalcubes_gc_create_simple_cube', PACKAGE = 'gdalcubes', files, datetime_values, bands, band_names, dx, dy, chunk_sizes) +gc_create_simple_cube <- function(files, datetime_values, bands, band_names, dx, dy, chunk_sizes, strict = TRUE) { + .Call('_gdalcubes_gc_create_simple_cube', PACKAGE = 'gdalcubes', files, datetime_values, bands, band_names, dx, dy, chunk_sizes, strict) } gc_create_fill_time_cube <- function(pin, method) { diff --git a/R/cube.R b/R/cube.R index 6f1430a..1b56789 100644 --- a/R/cube.R +++ b/R/cube.R @@ -7,6 +7,7 @@ #' @param view A data cube view defining the shape (spatiotemporal extent, resolution, and spatial reference), if missing, a default overview is used #' @param mask mask pixels of images based on band values, see \code{\link{image_mask}} #' @param chunking length-3 vector or a function returning a vector of length 3, defining the size of data cube chunks in the order time, y, x. +#' @param incomplete_ok logical, if TRUE (the default), chunks will ignore IO failures and simply use as much images as possible, otherwise the result will contain empty chunks if IO errors or similar occur. #' @return A proxy data cube object #' @details #' The following steps will be performed when the data cube is requested to read data of a chunk: @@ -38,7 +39,7 @@ #' #' @note This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result. #' @export -raster_cube <- function(image_collection, view, mask=NULL, chunking=.pkgenv$default_chunksize) { +raster_cube <- function(image_collection, view, mask=NULL, chunking=.pkgenv$default_chunksize, incomplete_ok = TRUE) { stopifnot(is.image_collection(image_collection)) if (is.function(chunking)) { @@ -60,10 +61,10 @@ raster_cube <- function(image_collection, view, mask=NULL, chunking=.pkgenv$defa x = NULL if (!missing(view)) { stopifnot(is.cube_view(view)) - x = gc_create_image_collection_cube(image_collection, as.integer(chunking), mask, view) + x = gc_create_image_collection_cube(image_collection, as.integer(chunking), mask, !incomplete_ok, view) } else { - x = gc_create_image_collection_cube(image_collection, as.integer(chunking), mask) + x = gc_create_image_collection_cube(image_collection, as.integer(chunking), mask, !incomplete_ok) } class(x) <- c("image_collection_cube", "cube", "xptr") return(x) @@ -100,6 +101,7 @@ raster_cube <- function(image_collection, view, mask=NULL, chunking=.pkgenv$defa #' @param chunking vector of length 3 defining the size of data cube chunks in the order time, y, x. #' @param dx optional target pixel size in x direction, by default (NULL) the original or highest resolution of images is used #' @param dy optional target pixel size in y direction, by default (NULL) the original or highest resolution of images is used +#' @param incomplete_ok logical, if TRUE (the default), chunks will ignore IO failures and simply use as much images as possible, otherwise the result will contain empty chunks if IO errors or similar occur. #' @return A proxy data cube object #' @examples #' # toy example, repeating the same image as a daily time series @@ -121,7 +123,7 @@ raster_cube <- function(image_collection, view, mask=NULL, chunking=.pkgenv$defa #' #' @note This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result. #' @export -stack_cube <- function(x, datetime_values, bands = NULL, band_names = NULL, chunking = c(1, 256, 256), dx=NULL, dy=NULL) { +stack_cube <- function(x, datetime_values, bands = NULL, band_names = NULL, chunking = c(1, 256, 256), dx=NULL, dy=NULL, incomplete_ok = TRUE) { if (length(datetime_values) != length(x)) { @@ -157,7 +159,7 @@ stack_cube <- function(x, datetime_values, bands = NULL, band_names = NULL, chun dy = -1.0 } - x = gc_create_simple_cube(x, datetime_values, bands, band_names, dx, dy, as.integer(chunking)) + x = gc_create_simple_cube(x, datetime_values, bands, band_names, dx, dy, as.integer(chunking), !incomplete_ok) class(x) <- c("simple_cube", "cube", "xptr") return(x) } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 27119e5..2dfafe3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -265,16 +265,17 @@ BEGIN_RCPP END_RCPP } // gc_create_image_collection_cube -SEXP gc_create_image_collection_cube(SEXP pin, Rcpp::IntegerVector chunk_sizes, SEXP mask, SEXP v); -RcppExport SEXP _gdalcubes_gc_create_image_collection_cube(SEXP pinSEXP, SEXP chunk_sizesSEXP, SEXP maskSEXP, SEXP vSEXP) { +SEXP gc_create_image_collection_cube(SEXP pin, Rcpp::IntegerVector chunk_sizes, SEXP mask, bool strict, SEXP v); +RcppExport SEXP _gdalcubes_gc_create_image_collection_cube(SEXP pinSEXP, SEXP chunk_sizesSEXP, SEXP maskSEXP, SEXP strictSEXP, SEXP vSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< SEXP >::type pin(pinSEXP); Rcpp::traits::input_parameter< Rcpp::IntegerVector >::type chunk_sizes(chunk_sizesSEXP); Rcpp::traits::input_parameter< SEXP >::type mask(maskSEXP); + Rcpp::traits::input_parameter< bool >::type strict(strictSEXP); Rcpp::traits::input_parameter< SEXP >::type v(vSEXP); - rcpp_result_gen = Rcpp::wrap(gc_create_image_collection_cube(pin, chunk_sizes, mask, v)); + rcpp_result_gen = Rcpp::wrap(gc_create_image_collection_cube(pin, chunk_sizes, mask, strict, v)); return rcpp_result_gen; END_RCPP } @@ -656,8 +657,8 @@ BEGIN_RCPP END_RCPP } // gc_create_simple_cube -SEXP gc_create_simple_cube(std::vector files, std::vector datetime_values, std::vector bands, std::vector band_names, double dx, double dy, Rcpp::IntegerVector chunk_sizes); -RcppExport SEXP _gdalcubes_gc_create_simple_cube(SEXP filesSEXP, SEXP datetime_valuesSEXP, SEXP bandsSEXP, SEXP band_namesSEXP, SEXP dxSEXP, SEXP dySEXP, SEXP chunk_sizesSEXP) { +SEXP gc_create_simple_cube(std::vector files, std::vector datetime_values, std::vector bands, std::vector band_names, double dx, double dy, Rcpp::IntegerVector chunk_sizes, bool strict); +RcppExport SEXP _gdalcubes_gc_create_simple_cube(SEXP filesSEXP, SEXP datetime_valuesSEXP, SEXP bandsSEXP, SEXP band_namesSEXP, SEXP dxSEXP, SEXP dySEXP, SEXP chunk_sizesSEXP, SEXP strictSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -668,7 +669,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double >::type dx(dxSEXP); Rcpp::traits::input_parameter< double >::type dy(dySEXP); Rcpp::traits::input_parameter< Rcpp::IntegerVector >::type chunk_sizes(chunk_sizesSEXP); - rcpp_result_gen = Rcpp::wrap(gc_create_simple_cube(files, datetime_values, bands, band_names, dx, dy, chunk_sizes)); + Rcpp::traits::input_parameter< bool >::type strict(strictSEXP); + rcpp_result_gen = Rcpp::wrap(gc_create_simple_cube(files, datetime_values, bands, band_names, dx, dy, chunk_sizes, strict)); return rcpp_result_gen; END_RCPP } @@ -875,7 +877,7 @@ static const R_CallMethodDef CallEntries[] = { {"_gdalcubes_gc_add_images", (DL_FUNC) &_gdalcubes_gc_add_images, 4}, {"_gdalcubes_gc_list_collection_formats", (DL_FUNC) &_gdalcubes_gc_list_collection_formats, 0}, {"_gdalcubes_gc_create_view", (DL_FUNC) &_gdalcubes_gc_create_view, 1}, - {"_gdalcubes_gc_create_image_collection_cube", (DL_FUNC) &_gdalcubes_gc_create_image_collection_cube, 4}, + {"_gdalcubes_gc_create_image_collection_cube", (DL_FUNC) &_gdalcubes_gc_create_image_collection_cube, 5}, {"_gdalcubes_gc_create_ncdf_cube", (DL_FUNC) &_gdalcubes_gc_create_ncdf_cube, 3}, {"_gdalcubes_gc_create_dummy_cube", (DL_FUNC) &_gdalcubes_gc_create_dummy_cube, 4}, {"_gdalcubes_gc_create_empty_cube", (DL_FUNC) &_gdalcubes_gc_create_empty_cube, 3}, @@ -904,7 +906,7 @@ static const R_CallMethodDef CallEntries[] = { {"_gdalcubes_gc_write_chunks_ncdf", (DL_FUNC) &_gdalcubes_gc_write_chunks_ncdf, 4}, {"_gdalcubes_gc_write_tif", (DL_FUNC) &_gdalcubes_gc_write_tif, 8}, {"_gdalcubes_gc_create_stream_cube", (DL_FUNC) &_gdalcubes_gc_create_stream_cube, 2}, - {"_gdalcubes_gc_create_simple_cube", (DL_FUNC) &_gdalcubes_gc_create_simple_cube, 7}, + {"_gdalcubes_gc_create_simple_cube", (DL_FUNC) &_gdalcubes_gc_create_simple_cube, 8}, {"_gdalcubes_gc_create_fill_time_cube", (DL_FUNC) &_gdalcubes_gc_create_fill_time_cube, 2}, {"_gdalcubes_gc_create_aggregate_time_cube", (DL_FUNC) &_gdalcubes_gc_create_aggregate_time_cube, 4}, {"_gdalcubes_gc_create_aggregate_space_cube", (DL_FUNC) &_gdalcubes_gc_create_aggregate_space_cube, 5}, diff --git a/src/gdalcubes.cpp b/src/gdalcubes.cpp index 11c0f2c..8eee658 100644 --- a/src/gdalcubes.cpp +++ b/src/gdalcubes.cpp @@ -844,7 +844,7 @@ SEXP gc_create_view(SEXP v) { // [[Rcpp::export]] -SEXP gc_create_image_collection_cube(SEXP pin, Rcpp::IntegerVector chunk_sizes, SEXP mask, SEXP v = R_NilValue) { +SEXP gc_create_image_collection_cube(SEXP pin, Rcpp::IntegerVector chunk_sizes, SEXP mask,bool strict = true, SEXP v = R_NilValue) { try { Rcpp::XPtr> aa = Rcpp::as>>(pin); @@ -858,6 +858,7 @@ SEXP gc_create_image_collection_cube(SEXP pin, Rcpp::IntegerVector chunk_sizes, x = new std::shared_ptr( image_collection_cube::create(*aa, cv)); } (*x)->set_chunk_size(chunk_sizes[0], chunk_sizes[1], chunk_sizes[2]); + (*x)->set_strict(strict); if (mask != R_NilValue) { @@ -1464,12 +1465,13 @@ SEXP gc_create_stream_cube(SEXP pin, std::string cmd) { // [[Rcpp::export]] SEXP gc_create_simple_cube(std::vector files,std::vector datetime_values, std::vector bands, std::vector band_names, - double dx, double dy, Rcpp::IntegerVector chunk_sizes) { + double dx, double dy, Rcpp::IntegerVector chunk_sizes, bool strict = true) { try { std::shared_ptr* x = new std::shared_ptr( simple_cube::create(files, datetime_values, bands, band_names, dx, dy)); (*x)->set_chunk_size(chunk_sizes[0], chunk_sizes[1], chunk_sizes[2]); + (*x)->set_strict(strict); Rcpp::XPtr< std::shared_ptr > p(x, true) ; return p; } diff --git a/src/gdalcubes/src/aggregate_space.cpp b/src/gdalcubes/src/aggregate_space.cpp index 8f83ca3..275442f 100644 --- a/src/gdalcubes/src/aggregate_space.cpp +++ b/src/gdalcubes/src/aggregate_space.cpp @@ -382,6 +382,15 @@ std::shared_ptr aggregate_space_cube::read_chunk(chunkid_t id) { chunkid_t input_chunk_id = _in_cube->chunk_id_from_coords({ccoords[0], ch_y, ch_x}); std::shared_ptr in_chunk = _in_cube->read_chunk(input_chunk_id); + // propagate chunk status + if (in_chunk->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (in_chunk->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + + if (in_chunk->empty()) { continue; } diff --git a/src/gdalcubes/src/aggregate_time.cpp b/src/gdalcubes/src/aggregate_time.cpp index a9a74f9..f2e34aa 100644 --- a/src/gdalcubes/src/aggregate_time.cpp +++ b/src/gdalcubes/src/aggregate_time.cpp @@ -397,6 +397,15 @@ std::shared_ptr aggregate_time_cube::read_chunk(chunkid_t id) { } std::shared_ptr in_chunk = chunk_cache[cur_in_chunk]; + + // propagate chunk status + if (in_chunk->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (in_chunk->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + if (in_chunk->empty()) { continue; } diff --git a/src/gdalcubes/src/apply_pixel.cpp b/src/gdalcubes/src/apply_pixel.cpp index 27b4aef..3a6cd8c 100644 --- a/src/gdalcubes/src/apply_pixel.cpp +++ b/src/gdalcubes/src/apply_pixel.cpp @@ -37,8 +37,9 @@ std::shared_ptr apply_pixel_cube::read_chunk(chunkid_t id) { return std::make_shared(); // chunk is outside of the view, we don't need to read anything. std::shared_ptr out = std::make_shared(); - std::shared_ptr in = _in_cube->read_chunk(id); + + out->set_status(in->status()); // propagate chunk status if (in->empty()) { return out; } diff --git a/src/gdalcubes/src/crop.cpp b/src/gdalcubes/src/crop.cpp index 5a497ee..a3ed756 100644 --- a/src/gdalcubes/src/crop.cpp +++ b/src/gdalcubes/src/crop.cpp @@ -58,6 +58,14 @@ std::shared_ptr crop_cube::read_chunk(chunkid_t id) { chunkid_t input_chunk_id = _in_cube->chunk_id_from_coords({ch_t, ch_y, ch_x}); std::shared_ptr in_chunk = _in_cube->read_chunk(input_chunk_id); + + // propagate chunk status + if (in_chunk->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (in_chunk->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } if (in_chunk->empty()) { continue; } diff --git a/src/gdalcubes/src/cube.cpp b/src/gdalcubes/src/cube.cpp index fb8fd18..12639d4 100644 --- a/src/gdalcubes/src/cube.cpp +++ b/src/gdalcubes/src/cube.cpp @@ -45,73 +45,6 @@ SOFTWARE. namespace gdalcubes { -void cube::write_chunks_gtiff(std::string dir, std::shared_ptr p) { - if (!filesystem::exists(dir)) { - filesystem::mkdir_recursive(dir); - } - - if (!filesystem::is_directory(dir)) { - throw std::string("ERROR in cube::write_chunks_gtiff(): output is not a directory."); - } - - GDALDriver *gtiff_driver = (GDALDriver *)GDALGetDriverByName("GTiff"); - if (gtiff_driver == NULL) { - throw std::string("ERROR: cannot find GDAL driver for GTiff."); - } - - if (!_st_ref->has_regular_space()) { - throw std::string("ERROR: GeoTIFF export currently does not support irregular spatial dimensions"); - } - - // NOTE: the following will only work as long as all cube st reference types with regular spatial dimensions inherit from cube_stref_regular class - std::shared_ptr stref = std::dynamic_pointer_cast(_st_ref); - - std::shared_ptr prg = config::instance()->get_default_progress_bar()->get(); - prg->set(0); // explicitly set to zero to show progress bar immediately - - std::function, std::mutex &)> f = [this, dir, prg, gtiff_driver, stref](chunkid_t id, std::shared_ptr dat, std::mutex &m) { - bounds_st cextent = this->bounds_from_chunk(id); // implemented in derived classes - double affine[6]; - affine[0] = cextent.s.left; - affine[3] = cextent.s.top; - affine[1] = stref->dx(); - affine[5] = -stref->dy(); - affine[2] = 0.0; - affine[4] = 0.0; - - CPLStringList out_co; - //out_co.AddNameValue("TILED", "YES"); - //out_co.AddNameValue("BLOCKXSIZE", "256"); - // out_co.AddNameValue("BLOCKYSIZE", "256"); - - for (uint16_t ib = 0; ib < dat->count_bands(); ++ib) { - for (uint32_t it = 0; it < dat->size()[1]; ++it) { - std::string out_file = filesystem::join(dir, (std::to_string(id) + "_" + std::to_string(ib) + "_" + std::to_string(it) + ".tif")); - - GDALDataset *gdal_out = gtiff_driver->Create(out_file.c_str(), dat->size()[3], dat->size()[2], 1, GDT_Float64, out_co.List()); - CPLErr res = gdal_out->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, dat->size()[3], dat->size()[2], ((double *)dat->buf()) + (ib * dat->size()[1] * dat->size()[2] * dat->size()[3] + it * dat->size()[2] * dat->size()[3]), dat->size()[3], dat->size()[2], GDT_Float64, 0, 0, NULL); - if (res != CE_None) { - GCBS_WARN("RasterIO (write) failed for band " + _bands.get(ib).name); - } - gdal_out->GetRasterBand(1)->SetNoDataValue(std::stod(_bands.get(ib).no_data_value)); - char *wkt_out; - OGRSpatialReference srs_out; - srs_out.SetFromUserInput(_st_ref->srs().c_str()); - srs_out.exportToWkt(&wkt_out); - - GDALSetProjection(gdal_out, wkt_out); - GDALSetGeoTransform(gdal_out, affine); - CPLFree(wkt_out); - - GDALClose(gdal_out); - } - } - prg->increment((double)1 / (double)this->count_chunks()); - }; - - p->apply(shared_from_this(), f); - prg->finalize(); -} void cube::write_tif_collection(std::string dir, std::string prefix, bool overviews, bool cog, @@ -437,301 +370,6 @@ void cube::write_tif_collection(std::string dir, std::string prefix, prg->finalize(); } -void cube::write_png_collection(std::string dir, std::string prefix, std::vector bands, std::vector zlim_min, - std::vector zlim_max, double gamma, std::vector na_color, bool na_transparent, - std::map creation_options, bool drop_empty_slices, std::shared_ptr p) { - if (zlim_min.empty()) { - zlim_min = {0}; - } - if (zlim_max.empty()) { - zlim_max = {255}; - } - - if (size_bands() > 1 && bands.empty()) { - if (size_bands() == 3) { - bands.push_back(_bands.get(0).name); - bands.push_back(_bands.get(1).name); - bands.push_back(_bands.get(2).name); - if (zlim_min.size() == 1) { - zlim_min.push_back(zlim_min[0]); - zlim_min.push_back(zlim_min[0]); - } - if (zlim_max.size() == 1) { - zlim_max.push_back(zlim_max[0]); - zlim_max.push_back(zlim_max[0]); - } - GCBS_WARN("Input data cube has three bands, using as (R,G,B) = (" + _bands.get(0).name + "," + _bands.get(1).name + "," + _bands.get(2).name + ")"); - } else { - GCBS_WARN("Input data cube has more than one band, using " + _bands.get(0).name + " to produce a grayscale image"); - bands.push_back(_bands.get(0).name); - } - } - if (bands.size() != 1 && bands.size() != 3) { - throw std::string("ERROR in cube::write_png_collection(): either one (grayscale) or three (RGB) bands expected, got " + std::to_string(bands.size())); - } - assert(bands.size() == 1 || bands.size() == 3); - - bool grayscale_as_rgb = false; - if (na_transparent) { - na_color.clear(); - } else { - if (!(na_color.empty() || na_color.size() == 1 || na_color.size() == 3)) { - throw std::string("ERROR in cube::write_png_collection(): either one (grayscale) or three (RGB) values expected for na_color, got " + std::to_string(na_color.size())); - } - if (na_color.size() == 3 && bands.size() == 1) { - grayscale_as_rgb = true; - bands.push_back(bands[0]); - bands.push_back(bands[0]); - if (zlim_min.size() == 1) { - zlim_min.push_back(zlim_min[0]); - zlim_min.push_back(zlim_min[0]); - } - if (zlim_max.size() == 1) { - zlim_max.push_back(zlim_max[0]); - zlim_max.push_back(zlim_max[0]); - } - } else if (na_color.size() == 1 && bands.size() == 3) { - na_color.clear(); - GCBS_WARN("Expected three (RGB) no data color; no data color will be ignored"); - } - } - assert(zlim_min.size() == bands.size()); - assert(zlim_max.size() == bands.size()); - - GDALDriver *png_driver = (GDALDriver *)GDALGetDriverByName("PNG"); - if (png_driver == NULL) { - throw std::string("ERROR: cannot find GDAL driver for PNG."); - } - // GDALDriver *mem_driver = (GDALDriver *)GDALGetDriverByName("MEM"); - // if (png_driver == NULL) { - // throw std::string("ERROR: cannot find GDAL driver for MEM."); - // } - GDALDriver *gtiff_driver = (GDALDriver *)GDALGetDriverByName("GTiff"); - if (png_driver == NULL) { - throw std::string("ERROR: cannot find GDAL driver for GTiff."); - } - - if (!filesystem::exists(dir)) { - filesystem::mkdir_recursive(dir); - } - if (!filesystem::is_directory(dir)) { - throw std::string("ERROR in cube::write_png_collection(): invalid output directory."); - } - - std::shared_ptr prg = config::instance()->get_default_progress_bar()->get(); - prg->set(0); // explicitly set to zero to show progress bar immediately - - CPLStringList out_co; - for (auto it = creation_options.begin(); it != creation_options.end(); ++it) { - std::string key = it->first; - std::transform(key.begin(), key.end(), key.begin(), (int (*)(int))std::toupper); - out_co.AddNameValue(it->first.c_str(), it->second.c_str()); - } - - std::vector band_nums; - for (uint16_t i = 0; i < bands.size(); ++i) { - if (!_bands.has(bands[i])) { - throw std::string("Data cube has no band with name '" + bands[i] + "'"); - } - band_nums.push_back(_bands.get_index(bands[i])); - } - uint8_t nbands = band_nums.size(); - int add_alpha = na_transparent ? 1 : 0; - - /* GTiff out */ - // avoid parallel RasterIO calls writing to the same file - std::map mtx; // time_slice_index -> mutex - - CPLStringList out_co_gtiff; - out_co_gtiff.AddNameValue("TILED", "YES"); - out_co_gtiff.AddNameValue("BLOCKXSIZE", std::to_string(chunk_size()[2]).c_str()); - out_co_gtiff.AddNameValue("BLOCKYSIZE", std::to_string(chunk_size()[1]).c_str()); - - std::string tempdir = filesystem::join(filesystem::get_tempdir(), utils::generate_unique_filename()); - filesystem::mkdir_recursive(tempdir); - - // create all temporary tif datasets - for (uint32_t it = 0; it < size_t(); ++it) { - std::string name = filesystem::join(tempdir, std::to_string(it) + ".tif"); - - GDALDataset *gdal_out = gtiff_driver->Create(name.c_str(), size_x(), size_y(), nbands + add_alpha, GDT_Byte, out_co.List()); - char *wkt_out; - OGRSpatialReference srs_out; - srs_out.SetFromUserInput(_st_ref->srs().c_str()); - srs_out.exportToWkt(&wkt_out); - GDALSetProjection(gdal_out, wkt_out); - - double affine[6]; - affine[0] = st_reference()->left(); - affine[3] = st_reference()->top(); - affine[1] = _st_ref->dx(); - affine[5] = -_st_ref->dy(); - affine[2] = 0.0; - affine[4] = 0.0; - GDALSetGeoTransform(gdal_out, affine); - CPLFree(wkt_out); - - if (nbands == 1) { - gdal_out->GetRasterBand(1)->SetColorInterpretation(GCI_GrayIndex); - } else { - gdal_out->GetRasterBand(1)->SetColorInterpretation(GCI_RedBand); - gdal_out->GetRasterBand(2)->SetColorInterpretation(GCI_GreenBand); - gdal_out->GetRasterBand(3)->SetColorInterpretation(GCI_BlueBand); - } - if (add_alpha == 1) { - gdal_out->GetRasterBand(nbands + 1)->SetColorInterpretation(GCI_AlphaBand); - } - - // gdal_out->GetRasterBand(ib + 1)->SetNoDataValue(packing.nodata[0]); - // gdal_out->GetRasterBand(ib + 1)->SetOffset(packing.offset[0]); - // gdal_out->GetRasterBand(ib + 1)->SetScale(packing.scale[0]); - // gdal_out->GetRasterBand(ib + 1)->Fill(packing.nodata[0]); - - GDALClose((GDALDatasetH)gdal_out); - } - - std::function, std::mutex &)> f = [this, dir, prg, &mtx, &tempdir, &nbands, &add_alpha, &na_color, &gamma, &band_nums, &zlim_min, &zlim_max, &grayscale_as_rgb](chunkid_t id, std::shared_ptr dat, std::mutex &m) { - if (!dat->empty()) { - for (uint32_t it = 0; it < dat->size()[1]; ++it) { - uint32_t cur_t_index = chunk_limits(id).low[0] + it; - std::string name = filesystem::join(tempdir, std::to_string(cur_t_index) + ".tif"); - mtx[cur_t_index].lock(); - GDALDataset *gdal_out = (GDALDataset *)GDALOpen(name.c_str(), GA_Update); - if (!gdal_out) { - GCBS_WARN("GDAL failed to open " + name); - mtx[cur_t_index].unlock(); - continue; - } - - uint8_t *buf_mask = nullptr; - if (add_alpha == 1) buf_mask = (uint8_t *)std::malloc(sizeof(uint8_t) * dat->size()[2] * dat->size()[3]); - - if (grayscale_as_rgb) { - uint8_t *buf = (uint8_t *)std::malloc(sizeof(uint8_t) * dat->size()[2] * dat->size()[3]); - for (uint16_t ib = 0; ib < nbands; ++ib) { - uint16_t b = band_nums[ib]; - for (uint32_t i = 0; i < dat->size()[2] * dat->size()[3]; ++i) { - double v = ((double *)(dat->buf()))[b * dat->size()[1] * dat->size()[2] * dat->size()[3] + - it * dat->size()[2] * dat->size()[3] + i]; - if (std::isnan(v)) { - if (add_alpha == 1) { - buf_mask[i] = 0; - } else if (!na_color.empty()) { - v = na_color[ib]; - } else { - v = 0; // TODO: require either na_color not empty or na transparent == true - } - } else { - if (add_alpha == 1) { - buf_mask[i] = 255; - } - v = ((v - zlim_min[ib]) / (zlim_max[ib] - zlim_min[ib])); - v = std::round(std::pow(v, gamma) * 255); - v = std::min(std::max(v, 0.0), 255.0); - } - buf[i] = v; - } - - CPLErr res = gdal_out->GetRasterBand(ib + 1)->RasterIO(GF_Write, chunk_limits(id).low[2], chunk_limits(id).low[1], dat->size()[3], dat->size()[2], - buf, dat->size()[3], dat->size()[2], GDT_Byte, 0, 0, NULL); - if (res != CE_None) { - GCBS_WARN("RasterIO (write) failed for " + name); - break; - } - } - std::free(buf); - } else { // modify band values in-place - for (uint16_t ib = 0; ib < nbands; ++ib) { - uint16_t b = band_nums[ib]; - for (uint32_t i = 0; i < dat->size()[2] * dat->size()[3]; ++i) { - double &v = ((double *)(dat->buf()))[b * dat->size()[1] * dat->size()[2] * dat->size()[3] + - it * dat->size()[2] * dat->size()[3] + i]; - if (std::isnan(v)) { - if (add_alpha == 1) { - buf_mask[i] = 0; - } else if (!na_color.empty()) { - v = na_color[ib]; - } else { - v = 0; // TODO: require either na_color not empty or na transparent == true - } - } else { - if (add_alpha == 1) { - buf_mask[i] = 255; - } - v = ((v - zlim_min[ib]) / (zlim_max[ib] - zlim_min[ib])); - v = std::round(std::pow(v, gamma) * 255); - v = std::min(std::max(v, 0.0), 255.0); - } - } - - CPLErr res = gdal_out->GetRasterBand(ib + 1)->RasterIO(GF_Write, chunk_limits(id).low[2], size_y() - chunk_limits(id).high[1] - 1, dat->size()[3], dat->size()[2], - ((double *)dat->buf()) + (b * dat->size()[1] * dat->size()[2] * dat->size()[3] + it * dat->size()[2] * dat->size()[3]), - dat->size()[3], dat->size()[2], GDT_Float64, 0, 0, NULL); - if (res != CE_None) { - GCBS_WARN("RasterIO (write) failed for " + name); - break; - } - } - } - if (add_alpha == 1) { - CPLErr res = gdal_out->GetRasterBand(nbands + 1)->RasterIO(GF_Write, chunk_limits(id).low[2], size_y() - chunk_limits(id).high[1] - 1, dat->size()[3], dat->size()[2], buf_mask, dat->size()[3], dat->size()[2], GDT_Byte, 0, 0, NULL); - if (res != CE_None) { - GCBS_WARN("RasterIO (write) failed for chunk " + std::string(gdal_out->GetDescription())); - } - std::free(buf_mask); - } - GDALClose(gdal_out); - mtx[cur_t_index].unlock(); - } - } - prg->increment((double)1 / (double)this->count_chunks()); - }; - p->apply(shared_from_this(), f); - - for (uint32_t it = 0; it < size_t(); ++it) { - std::string name = filesystem::join(tempdir, std::to_string(it) + ".tif"); - - CPLStringList translate_args; - translate_args.AddString("-of"); - translate_args.AddString("PNG"); - for (auto it = creation_options.begin(); it != creation_options.end(); ++it) { - translate_args.AddString("-co"); - std::string v = it->first + "=" + it->second; - translate_args.AddString(v.c_str()); - } - - GDALTranslateOptions *trans_options = GDALTranslateOptionsNew(translate_args.List(), NULL); - if (trans_options == NULL) { - GCBS_WARN("Cannot create gdal_translate options."); - continue; - } - GDALDataset *dataset = (GDALDataset *)GDALOpen(name.c_str(), GA_ReadOnly); - if (!dataset) { - GCBS_WARN("Cannot open GDAL dataset '" + name + "'."); - GDALTranslateOptionsFree(trans_options); - continue; - } - std::string outfile = filesystem::join(dir, prefix + _st_ref->datetime_at_index(it).to_string() + ".png"); - if (!filesystem::exists(dir)) { - filesystem::mkdir(dir); - } - GDALDatasetH out = GDALTranslate(outfile.c_str(), (GDALDatasetH)dataset, trans_options, NULL); - if (!out) { - GCBS_WARN("Cannot translate GDAL dataset '" + name + "'."); - GDALClose((GDALDatasetH)dataset); - GDALTranslateOptionsFree(trans_options); - } - GDALClose(out); - GDALTranslateOptionsFree(trans_options); - - GDALClose((GDALDatasetH)dataset); - filesystem::remove(name); - } - filesystem::remove(tempdir); - - prg->set(1.0); - prg->finalize(); -} void cube::write_netcdf_file(std::string path, uint8_t compression_level, bool with_VRT, bool write_bounds, packed_export packing, bool drop_empty_slices, std::shared_ptr p) { @@ -898,6 +536,11 @@ void cube::write_netcdf_file(std::string path, uint8_t compression_level, bool w nc_def_var(ncout, "x_bnds", NC_DOUBLE, 2, d_xbnds, &v_xbnds); } + int d_chunkstatus; + nc_def_dim(ncout, "chunks", count_chunks(), &d_chunkstatus); + int v_chunkstatus; + nc_def_var(ncout, "chunk_status", NC_INT, 1, &d_chunkstatus, &v_chunkstatus); + std::string att_source = "gdalcubes " + std::to_string(GDALCUBES_VERSION_MAJOR) + "." + std::to_string(GDALCUBES_VERSION_MINOR) + "." + std::to_string(GDALCUBES_VERSION_PATCH); nc_put_att_text(ncout, NC_GLOBAL, "Conventions", std::strlen("CF-1.6"), "CF-1.6"); @@ -1073,9 +716,12 @@ void cube::write_netcdf_file(std::string path, uint8_t compression_level, bool w if (dim_x_bnds) std::free(dim_x_bnds); } - std::function, std::mutex &)> f = [this, op, prg, &v_bands, ncout, &packing](chunkid_t id, std::shared_ptr dat, std::mutex &m) { + std::function, std::mutex &)> f = [this, op, prg, &v_bands, ncout, &packing, &v_chunkstatus](chunkid_t id, std::shared_ptr dat, std::mutex &m) { // TODO: check if it is OK to simply not write anything to netCDF or if we need to fill dat explicity with no data values, check also for packed output + int nc_chunk_status = (int)dat->status(); + std::size_t nc_chunk_id = std::size_t(id); + nc_put_var1_int(ncout, v_chunkstatus, &nc_chunk_id, &nc_chunk_status); if (!dat->empty()) { chunk_size_btyx csize = dat->size(); bounds_nd climits = chunk_limits(id); @@ -1923,6 +1569,9 @@ void chunk_data::write_ncdf(std::string path, uint8_t compression_level, bool fo std::string att_source = "gdalcubes " + std::to_string(GDALCUBES_VERSION_MAJOR) + "." + std::to_string(GDALCUBES_VERSION_MINOR) + "." + std::to_string(GDALCUBES_VERSION_PATCH); nc_put_att_text(ncout, NC_GLOBAL, "source", std::strlen(att_source.c_str()), att_source.c_str()); + + int nc_chunk_status = (int)_status; + nc_put_att(ncout, NC_GLOBAL, "chunk_status", NC_INT, 1, &nc_chunk_status); int d_all[] = {d_b, d_t, d_y, d_x}; int v; @@ -1954,6 +1603,15 @@ void chunk_data::read_ncdf(std::string path) { return; } + int s = 0; + if (nc_get_att_int(ncfile, NC_GLOBAL, "chunk_status", &s)) { + set_status(static_cast(s)); + } + else { + GCBS_DEBUG("NetCDF chunk file does not contain chunk status data. "); + set_status(chunk_data::chunk_status::UNKNOWN); + } + int dim_id_x = -1; int dim_id_y = -1; int dim_id_t = -1; diff --git a/src/gdalcubes/src/cube.h b/src/gdalcubes/src/cube.h index a474cb5..901a356 100644 --- a/src/gdalcubes/src/cube.h +++ b/src/gdalcubes/src/cube.h @@ -317,10 +317,18 @@ class band_collection { */ class chunk_data { public: + + enum class chunk_status { + OK = 0, + ERROR = 1, + INCOMPLETE = 2, + UNKNOWN = 128 + }; + /** * @brief Default constructor that creates an empty chunk */ - chunk_data() : _buf(nullptr), _size({{0, 0, 0, 0}}) {} + chunk_data() : _buf(nullptr), _size({{0, 0, 0, 0}}), _status(chunk_status::OK) {} ~chunk_data() { if (_buf && _size[0] * _size[1] * _size[2] * _size[3] > 0) std::free(_buf); @@ -427,9 +435,19 @@ class chunk_data { */ void read_ncdf(std::string path); + + void set_status(chunk_status s) { + _status = s; + } + + chunk_status status() { + return _status; + } + private: void *_buf; chunk_size_btyx _size; + chunk_status _status; // error flags }; /** @@ -777,14 +795,6 @@ class cube : public std::enable_shared_from_this { */ std::shared_ptr read_window(std::array lower, std::array upper); - /** - * @brief Write a data cube as a set of GeoTIFF files under a given directory - * - * @param dir directory where to store the files - * @param p chunk processor instance, defaults to the global configuration - */ - void write_chunks_gtiff(std::string dir, - std::shared_ptr p = config::instance()->get_default_chunk_processor()); /** * @brief Writes a data cube as a collection of GeoTIFF files @@ -845,38 +855,6 @@ class cube : public std::enable_shared_from_this { void write_single_chunk_netcdf(chunkid_t id, std::string path, uint8_t compression_level = 0); - /** - * @brief Writes a data cube as a collection of PNG files - * - * Each time slice of the cube will be written to one GeoTIFF file. - * @param dir output path - * @param prefix name prefix for created files - * @param bands vector with band names; must contain either no, one, or three elements; If empty, the first band is used to produce a grayscale image - * @param zlim_min minimum color value per channel (if vector has size 1 but three bands are provided, the same value is used for all bands) - * @param zlim_max maximum color value per channel (if vector has size 1 but three bands are provided, the same value is used for all bands) - * @param gamma gamma correction factor to be applied - * @param na_color = {} - * @param na_transparent if true, the resulting image will cintain an alpha channel, where NA values will be assigned alpha = 0% - * @param creation_options key value pairs passed to GDAL as PNG creation options (see https://gdal.org/drivers/raster/png.html) - * @param drop_empty_slices if true, empty time slices will be skipped; not yet implemented - * @param p chunk processor instance, defaults to the global configuration - * - * @note argument `drop_empty_slices` is not yet implemented. - * - * @note This function can be used to create single-band grayscale or three-bands color RGB plots. - * If bands is empty, the first band of the data cube will be used to produce a single-band grayscale image. If - * three bands are provided, they are interpreted in the order R,G,B. Depsite the number of provided bands, - * the color scale defined by zlim_min and zlim_max may be provided as single element vectors in which case the values are used for all color channels, or as - * three element vectors with separate scales for each RGB channel. If na_color is not empty, it can be used to set the color values of NA pixels either as a single gray value or as three separate RGB values. - * In the case of producing a single band grayscale image but three separate vlaues as na_color, the output image will contain three channels with gray values (R=G=B) for non-NA pixels and the na_color for NA values. - * Notice that if na_transparent is true, na_color is ignored. - * Further creation options can be used, e.g., to control the compression level of PNG files. - */ - void write_png_collection(std::string dir, std::string prefix = "", - std::vector bands = {}, std::vector zlim_min = {0}, std::vector zlim_max = {255}, double gamma = 1, std::vector na_color = {}, bool na_transparent = true, - std::map creation_options = std::map(), - bool drop_empty_slices = false, - std::shared_ptr p = config::instance()->get_default_chunk_processor()); /** * Get the cube's bands diff --git a/src/gdalcubes/src/cube_factory.cpp b/src/gdalcubes/src/cube_factory.cpp index a29368e..c6275e0 100644 --- a/src/gdalcubes/src/cube_factory.cpp +++ b/src/gdalcubes/src/cube_factory.cpp @@ -247,6 +247,9 @@ void cube_factory::register_default() { dy = j["dy"].number_value(); } auto x = simple_cube::create(files, datetime, bands, band_names, dx, dy); + if (!j["strict"].is_null()) { + x->set_strict(j["strict"].bool_value()); + } if (!j["chunk_size"].is_null()) { if (j["chunk_size"].array_items().size() == 3) { x->set_chunk_size(j["chunk_size"][0].int_value(), j["chunk_size"][1].int_value(), j["chunk_size"][2].int_value()); @@ -342,7 +345,9 @@ void cube_factory::register_default() { cube_view v = cube_view::read_json_string(j["view"].dump()); auto x = image_collection_cube::create(j["file"].string_value(), v); x->set_chunk_size(j["chunk_size"][0].int_value(), j["chunk_size"][1].int_value(), j["chunk_size"][2].int_value()); - + if (!j["strict"].is_null()) { + x->set_strict(j["strict"].bool_value()); + } if (!j["mask"].is_null()) { if (j["mask"]["mask_type"].is_null()) { GCBS_WARN("ERROR in cube_generators[\"image_collection\"](): missing mask type, mask will be ignored"); diff --git a/src/gdalcubes/src/fill_time.cpp b/src/gdalcubes/src/fill_time.cpp index 8be76fa..5338f76 100644 --- a/src/gdalcubes/src/fill_time.cpp +++ b/src/gdalcubes/src/fill_time.cpp @@ -26,7 +26,10 @@ std::shared_ptr fill_time_cube::read_chunk(chunkid_t id) { //std::vector> r_chunks; std::unordered_map> in_chunks; - in_chunks.insert(std::pair>(id, _in_cube->read_chunk(id))); + auto ic = _in_cube->read_chunk(id); + out->set_status(ic->status()); // propagate chunk status + + in_chunks.insert(std::pair>(id, ic)); if (in_chunks[id]->empty()) { // if input chunk is empty, fill with NANs in_chunks[id]->size(size_btyx); @@ -57,7 +60,15 @@ std::shared_ptr fill_time_cube::read_chunk(chunkid_t id) { while (prev_chunk >= 0 && !found) { // load chunk (only if needed) if (in_chunks.find(prev_chunk) == in_chunks.end()) { - in_chunks.insert(std::pair>(prev_chunk, _in_cube->read_chunk(prev_chunk))); + auto ic = _in_cube->read_chunk(prev_chunk); + // propagate chunk status + if (ic->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (ic->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + in_chunks.insert(std::pair>(prev_chunk, ic)); } if (!in_chunks[prev_chunk]->empty()) { prev_t = _in_cube->chunk_size()[0] - 1; @@ -87,6 +98,14 @@ std::shared_ptr fill_time_cube::read_chunk(chunkid_t id) { while (next_chunk < (int32_t)_in_cube->count_chunks() && !found) { // load chunk (only if needed) if (in_chunks.find(next_chunk) == in_chunks.end()) { + auto ic = _in_cube->read_chunk(next_chunk); + // propagate chunk status + if (ic->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (ic->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } in_chunks.insert(std::pair>(next_chunk, _in_cube->read_chunk(next_chunk))); } if (!in_chunks[next_chunk]->empty()) { diff --git a/src/gdalcubes/src/filter_geom.cpp b/src/gdalcubes/src/filter_geom.cpp index 7a85f8c..57cbd44 100644 --- a/src/gdalcubes/src/filter_geom.cpp +++ b/src/gdalcubes/src/filter_geom.cpp @@ -215,6 +215,14 @@ std::shared_ptr filter_geom_cube::read_chunk(chunkid_t id) { } std::shared_ptr in = _in_cube->read_chunk(id); + // propagate chunk status + if (in->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (in->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + if (in->empty()) { GDALClose(in_ogr_dataset); return out; diff --git a/src/gdalcubes/src/filter_pixel.cpp b/src/gdalcubes/src/filter_pixel.cpp index d4013fb..72c3618 100644 --- a/src/gdalcubes/src/filter_pixel.cpp +++ b/src/gdalcubes/src/filter_pixel.cpp @@ -39,6 +39,8 @@ std::shared_ptr filter_pixel_cube::read_chunk(chunkid_t id) { std::shared_ptr out = std::make_shared(); std::shared_ptr in = _in_cube->read_chunk(id); + + out->set_status(in->status()); // propagate chunk status if (in->empty()) { return out; } @@ -98,7 +100,9 @@ std::shared_ptr filter_pixel_cube::read_chunk(chunkid_t id) { } // check if chunk is completely NAN and if yes, return empty chunk if (out->all_nan()) { + chunk_data::chunk_status s = out->status(); out = std::make_shared(); + out->set_status(s); } return out; diff --git a/src/gdalcubes/src/gdalcubes.cpp b/src/gdalcubes/src/gdalcubes.cpp index 6996a06..05c6fc8 100644 --- a/src/gdalcubes/src/gdalcubes.cpp +++ b/src/gdalcubes/src/gdalcubes.cpp @@ -37,9 +37,6 @@ #include "image_collection_cube.h" #include "image_collection_ops.h" #include "stream.h" -#ifndef GDALCUBES_NO_SWARM -#include "swarm.h" -#endif #include "utils.h" using namespace gdalcubes; diff --git a/src/gdalcubes/src/gdalcubes.h b/src/gdalcubes/src/gdalcubes.h index c221eea..570687f 100644 --- a/src/gdalcubes/src/gdalcubes.h +++ b/src/gdalcubes/src/gdalcubes.h @@ -58,7 +58,6 @@ #include "stream_reduce_space.h" #include "stream_reduce_time.h" #include "utils.h" -#include "vector_queries.h" #include "window_space.h" #include "window_time.h" diff --git a/src/gdalcubes/src/image_collection_cube.cpp b/src/gdalcubes/src/image_collection_cube.cpp index 41f328c..6b1c6ff 100644 --- a/src/gdalcubes/src/image_collection_cube.cpp +++ b/src/gdalcubes/src/image_collection_cube.cpp @@ -34,16 +34,16 @@ namespace gdalcubes { -image_collection_cube::image_collection_cube(std::shared_ptr ic, cube_view v) : cube(std::make_shared(v)), _collection(ic), _input_bands(), _mask(nullptr), _mask_band("") { load_bands(); } -image_collection_cube::image_collection_cube(std::string icfile, cube_view v) : cube(std::make_shared(v)), _collection(std::make_shared(icfile)), _input_bands(), _mask(nullptr), _mask_band("") { load_bands(); } -image_collection_cube::image_collection_cube(std::shared_ptr ic, std::string vfile) : cube(std::make_shared(cube_view::read_json(vfile))), _collection(ic), _input_bands(), _mask(nullptr), _mask_band("") { load_bands(); } -image_collection_cube::image_collection_cube(std::string icfile, std::string vfile) : cube(std::make_shared(cube_view::read_json(vfile))), _collection(std::make_shared(icfile)), _input_bands(), _mask(nullptr), _mask_band("") { load_bands(); } -image_collection_cube::image_collection_cube(std::shared_ptr ic) : cube(), _collection(ic), _input_bands(), _mask(nullptr), _mask_band("") { +image_collection_cube::image_collection_cube(std::shared_ptr ic, cube_view v) : cube(std::make_shared(v)), _collection(ic), _input_bands(), _mask(nullptr), _mask_band(""), _strict(true) { load_bands(); } +image_collection_cube::image_collection_cube(std::string icfile, cube_view v) : cube(std::make_shared(v)), _collection(std::make_shared(icfile)), _input_bands(), _mask(nullptr), _mask_band(""), _strict(true) { load_bands(); } +image_collection_cube::image_collection_cube(std::shared_ptr ic, std::string vfile) : cube(std::make_shared(cube_view::read_json(vfile))), _collection(ic), _input_bands(), _mask(nullptr), _mask_band(""), _strict(true) { load_bands(); } +image_collection_cube::image_collection_cube(std::string icfile, std::string vfile) : cube(std::make_shared(cube_view::read_json(vfile))), _collection(std::make_shared(icfile)), _input_bands(), _mask(nullptr), _mask_band(""), _strict(true) { load_bands(); } +image_collection_cube::image_collection_cube(std::shared_ptr ic) : cube(), _collection(ic), _input_bands(), _mask(nullptr), _mask_band(""), _strict(true) { st_reference(std::make_shared(image_collection_cube::default_view(_collection))); load_bands(); } -image_collection_cube::image_collection_cube(std::string icfile) : cube(), _collection(std::make_shared(icfile)), _input_bands(), _mask(nullptr), _mask_band("") { +image_collection_cube::image_collection_cube(std::string icfile) : cube(), _collection(std::make_shared(icfile)), _input_bands(), _mask(nullptr), _mask_band(""), _strict(true) { st_reference(std::make_shared(image_collection_cube::default_view(_collection))); load_bands(); } @@ -317,7 +317,7 @@ std::shared_ptr image_collection_cube::read_chunk(chunkid_t id) { std::shared_ptr out = std::make_shared(); if (id >= count_chunks()) { // chunk is outside of the cube, we don't need to read anything. - GCBS_WARN("Chunk id " + std::to_string(id) + " is out of range"); + GCBS_DEBUG("Chunk id " + std::to_string(id) + " is out of range"); return out; } @@ -377,6 +377,7 @@ std::shared_ptr image_collection_cube::read_chunk(chunkid_t id) { } uint32_t i = 0; + uint32_t count_success = 0; // count successful image reads while (i < datasets.size()) { std::pair mask_dataset_band; mask_dataset_band.first = ""; @@ -422,12 +423,22 @@ std::shared_ptr image_collection_cube::read_chunk(chunkid_t id) { std::string bandsel_vrt_name = ""; GDALDataset *g = (GDALDataset *)GDALOpen(it->first.c_str(), GA_ReadOnly); if (!g) { - GCBS_WARN("GDAL cannot open '" + it->first + "', image will be ignored."); + GCBS_WARN("GDAL could not open '" + it->first + "': ERROR no " + std::to_string(CPLGetLastErrorNo()) + ":" + CPLGetLastErrorMsg()); + if (_strict) { + if (img_buf) std::free(img_buf); + if (mask_buf) std::free(mask_buf); + delete agg; + out = std::make_shared(); + out->set_status(chunk_data::chunk_status::ERROR); + return out; + } + GCBS_WARN("Dataset '" + it->first + "' will be ignored."); + out->set_status(chunk_data::chunk_status::INCOMPLETE); continue; } if (g->GetRasterCount() == 0) { - GCBS_WARN("GDAL dataset'" + it->first + "' does not contain any raster bands and will be ignored."); + GCBS_DEBUG("GDAL dataset'" + it->first + "' does not contain any raster bands and will be ignored."); continue; } @@ -495,6 +506,20 @@ std::shared_ptr image_collection_cube::read_chunk(chunkid_t id) { cextent.s.top, cextent.s.bottom, size_btyx[3], size_btyx[2], resampling::to_string(view()->resampling_method()), nodata_value_list); } + if (!gdal_out) { + GCBS_WARN("GDAL could not warp '" + it->first + "': ERROR no " + std::to_string(CPLGetLastErrorNo()) + ":" + CPLGetLastErrorMsg()); + if (_strict) { + if (img_buf) std::free(img_buf); + if (mask_buf) std::free(mask_buf); + delete agg; + out = std::make_shared(); + out->set_status(chunk_data::chunk_status::ERROR); + return out; + } + GCBS_WARN("Dataset '" + it->first + "' will be ignored."); + out->set_status(chunk_data::chunk_status::INCOMPLETE); + continue; + } // For each band, call RasterIO to read and copy data to the right position in the buffers for (uint16_t b = 0; b < it->second.size(); ++b) { @@ -511,7 +536,18 @@ std::shared_ptr image_collection_cube::read_chunk(chunkid_t id) { res = gdal_out->GetRasterBand(std::get<1>(it->second[b]))->RasterIO(GF_Read, 0, 0, size_btyx[3], size_btyx[2], ((double *)img_buf) + b_internal * size_btyx[2] * size_btyx[3], size_btyx[3], size_btyx[2], GDT_Float64, 0, 0, NULL); } if (res != CE_None) { - GCBS_WARN("RasterIO (read) failed for " + std::string(gdal_out->GetDescription())); + GCBS_WARN("RasterIO (read) failed for '" + std::string(gdal_out->GetDescription()) + "': ERROR no " + std::to_string(CPLGetLastErrorNo()) + ":" + CPLGetLastErrorMsg()); + if (_strict) { + if (img_buf) std::free(img_buf); + if (mask_buf) std::free(mask_buf); + delete agg; + out = std::make_shared(); + out->set_status(chunk_data::chunk_status::ERROR); + return out; + } + GCBS_WARN("Dataset '" + it->first + "' will be ignored."); + out->set_status(chunk_data::chunk_status::INCOMPLETE); + continue; } } if (!bandsel_vrt_name.empty()) { @@ -532,7 +568,15 @@ std::shared_ptr image_collection_cube::read_chunk(chunkid_t id) { GDALDataset *bandsel_vrt = nullptr; GDALDataset *g = (GDALDataset *)GDALOpen(mask_dataset_band.first.c_str(), GA_ReadOnly); if (!g) { - GCBS_WARN("GDAL cannot open '" + mask_dataset_band.first + "', mask will be ignored"); + GCBS_WARN("GDAL could not open '" + mask_dataset_band.first + "': ERROR no " + std::to_string(CPLGetLastErrorNo()) + ":" + CPLGetLastErrorMsg()); + if (_strict) { + if (img_buf) std::free(img_buf); + if (mask_buf) std::free(mask_buf); + delete agg; + return std::make_shared(); + } + GCBS_WARN("Mask dataset '" + mask_dataset_band.first + "' will be ignored."); + continue; } else { // If input dataset has more bands than requested @@ -573,9 +617,35 @@ std::shared_ptr image_collection_cube::read_chunk(chunkid_t id) { cextent.s.top, cextent.s.bottom, size_btyx[3], size_btyx[2], "near", std::vector()); } + if (!gdal_out) { + GCBS_WARN("GDAL could not warp '" + mask_dataset_band.first + "': ERROR no " + std::to_string(CPLGetLastErrorNo()) + ":" + CPLGetLastErrorMsg()); + if (_strict) { + if (img_buf) std::free(img_buf); + if (mask_buf) std::free(mask_buf); + delete agg; + out = std::make_shared(); + out->set_status(chunk_data::chunk_status::ERROR); + return out; + } + GCBS_WARN("Mask dataset '" + mask_dataset_band.first + "' will be ignored."); + out->set_status(chunk_data::chunk_status::INCOMPLETE); + continue; + } CPLErr res = gdal_out->GetRasterBand(mask_dataset_band.second)->RasterIO(GF_Read, 0, 0, size_btyx[3], size_btyx[2], mask_buf, size_btyx[3], size_btyx[2], GDT_Float64, 0, 0, NULL); + if (res != CE_None) { - GCBS_WARN("RasterIO (read) failed for " + std::string(gdal_out->GetDescription())); + GCBS_WARN("RasterIO (read) failed for '" + std::string(gdal_out->GetDescription()) + "': ERROR no " + std::to_string(CPLGetLastErrorNo()) + ":" + CPLGetLastErrorMsg()); + if (_strict) { + if (img_buf) std::free(img_buf); + if (mask_buf) std::free(mask_buf); + delete agg; + out = std::make_shared(); + out->set_status(chunk_data::chunk_status::ERROR); + return out; + } + GCBS_WARN("Mask dataset '" + mask_dataset_band.first + "' will be ignored."); + out->set_status(chunk_data::chunk_status::INCOMPLETE); + continue; } GDALClose(gdal_out); _mask->apply((double *)mask_buf, (double *)img_buf, size_btyx[0], size_btyx[2], size_btyx[3]); @@ -585,6 +655,13 @@ std::shared_ptr image_collection_cube::read_chunk(chunkid_t id) { // feed the aggregator agg->update(out->buf(), img_buf, itime); + + count_success++; + } + + + if (out->status() == chunk_data::chunk_status::INCOMPLETE && count_success == 0) { + out->set_status(chunk_data::chunk_status::ERROR); } agg->finalize(out->buf()); diff --git a/src/gdalcubes/src/image_collection_cube.h b/src/gdalcubes/src/image_collection_cube.h index faed513..835fdc7 100644 --- a/src/gdalcubes/src/image_collection_cube.h +++ b/src/gdalcubes/src/image_collection_cube.h @@ -286,6 +286,10 @@ class image_collection_cube : public cube { _chunk_size = {t, y, x}; } + void set_strict(bool s) { + _strict = s; + } + json11::Json make_constructible_json() override { if (_collection->is_temporary()) { throw std::string("ERROR in image_collection_cube::make_constructible_json(): image collection is temporary, please export as file using write() first."); @@ -300,6 +304,7 @@ class image_collection_cube : public cube { out["mask"] = _mask->as_json(); out["mask_band"] = _mask_band; } + out["strict"] = _strict; return out; } @@ -314,6 +319,8 @@ class image_collection_cube : public cube { std::shared_ptr _mask; std::string _mask_band; + + bool _strict; }; } // namespace gdalcubes diff --git a/src/gdalcubes/src/join_bands.cpp b/src/gdalcubes/src/join_bands.cpp index 378f7e0..80020cc 100644 --- a/src/gdalcubes/src/join_bands.cpp +++ b/src/gdalcubes/src/join_bands.cpp @@ -54,6 +54,14 @@ std::shared_ptr join_bands_cube::read_chunk(chunkid_t id) { bool allempty = true; for (uint16_t i = 0; i < _in.size(); ++i) { std::shared_ptr dat = _in[i]->read_chunk(id); + // propagate chunk status + if (dat->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (dat->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + if (!dat->empty()) { allempty = false; std::memcpy(((double *)out->buf()) + offset, ((double *)dat->buf()), dat->size()[0] * dat->size()[1] * dat->size()[2] * dat->size()[3] * sizeof(double)); @@ -61,8 +69,10 @@ std::shared_ptr join_bands_cube::read_chunk(chunkid_t id) { offset += _in[i]->size_bands() * size_tyx[0] * size_tyx[1] * size_tyx[2]; } if (allempty) { + auto s = out->status(); // propagated empty chunk if all input chunks are emtpy out = std::make_shared(); + out->set_status(s); } return out; } diff --git a/src/gdalcubes/src/ncdf_cube.cpp b/src/gdalcubes/src/ncdf_cube.cpp index 6ffd828..fbad72e 100644 --- a/src/gdalcubes/src/ncdf_cube.cpp +++ b/src/gdalcubes/src/ncdf_cube.cpp @@ -460,7 +460,7 @@ std::shared_ptr ncdf_cube::read_chunk(chunkid_t id) { std::shared_ptr out = std::make_shared(); if (id >= count_chunks()) { // chunk is outside of the cube, we don't need to read anything. - GCBS_WARN("Chunk id " + std::to_string(id) + " is out of range"); + GCBS_DEBUG("Chunk id " + std::to_string(id) + " is out of range"); return out; } @@ -491,8 +491,19 @@ std::shared_ptr ncdf_cube::read_chunk(chunkid_t id) { throw std::string("Failed to open netCDF file '" + _path + "'; nc_open() returned " + std::to_string(retval)); } + int varid = -1; + if (nc_inq_varid(ncfile, "chunk_status", &varid) == NC_NOERR) { + int s = 0; + std::size_t nc_chunk_id = std::size_t(id); + nc_get_var1_int(ncfile, varid, &nc_chunk_id, &s); + out->set_status(static_cast(s)); + } + else { + GCBS_DEBUG("NetCDF input file does not contain chunk status data. "); + out->set_status(chunk_data::chunk_status::UNKNOWN); + } + for (uint16_t i = 0; i < size_btyx[0]; ++i) { - int varid = -1; if (nc_inq_varid(ncfile, _bands.get(i).name.c_str(), &varid) == NC_NOERR) { nc_get_vara_double(ncfile, varid, startp, countp, (&((double *)(out->buf()))[i * size_btyx[1] * size_btyx[2] * size_btyx[3]])); } else { @@ -527,7 +538,9 @@ std::shared_ptr ncdf_cube::read_chunk(chunkid_t id) { // check if chunk is completely NAN and if yes, return empty chunk if (out->all_nan()) { + auto s = out->status(); out = std::make_shared(); + out->set_status(s); } return out; } diff --git a/src/gdalcubes/src/reduce_space.cpp b/src/gdalcubes/src/reduce_space.cpp index bca394a..bbef997 100644 --- a/src/gdalcubes/src/reduce_space.cpp +++ b/src/gdalcubes/src/reduce_space.cpp @@ -396,6 +396,15 @@ std::shared_ptr reduce_space_cube::read_chunk(chunkid_t id) { bool initialized = false; // lazy initialization after the first non-empty chunk for (chunkid_t i = id * _in_cube->count_chunks_x() * _in_cube->count_chunks_y(); i < (id + 1) * _in_cube->count_chunks_x() * _in_cube->count_chunks_y(); ++i) { std::shared_ptr x = _in_cube->read_chunk(i); + + // propagate chunk status + if (x->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (x->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + if (!x->empty()) { if (!initialized) { // Fill buffers with NAN @@ -416,7 +425,9 @@ std::shared_ptr reduce_space_cube::read_chunk(chunkid_t id) { } } if (empty) { + auto s = out->status(); out = std::make_shared(); + out->set_status(s); } else { for (uint16_t i = 0; i < _reducer_bands.size(); ++i) { diff --git a/src/gdalcubes/src/reduce_time.cpp b/src/gdalcubes/src/reduce_time.cpp index 9d0dc53..eb30596 100644 --- a/src/gdalcubes/src/reduce_time.cpp +++ b/src/gdalcubes/src/reduce_time.cpp @@ -580,6 +580,15 @@ std::shared_ptr reduce_time_cube::read_chunk(chunkid_t id) { bool initialized = false; // lazy initialization after the first non-empty chunk for (chunkid_t i = id; i < _in_cube->count_chunks(); i += _in_cube->count_chunks_x() * _in_cube->count_chunks_y()) { std::shared_ptr x = _in_cube->read_chunk(i); + + // propagate chunk status + if (x->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (x->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + if (!x->empty()) { if (!initialized) { // Fill buffers with NAN @@ -600,7 +609,9 @@ std::shared_ptr reduce_time_cube::read_chunk(chunkid_t id) { } } if (empty) { + auto s = out->status(); out = std::make_shared(); + out->set_status(s); } else { for (uint16_t i = 0; i < _reducer_bands.size(); ++i) { diff --git a/src/gdalcubes/src/select_bands.cpp b/src/gdalcubes/src/select_bands.cpp index 9242165..933d864 100644 --- a/src/gdalcubes/src/select_bands.cpp +++ b/src/gdalcubes/src/select_bands.cpp @@ -37,12 +37,17 @@ std::shared_ptr select_bands_cube::read_chunk(chunkid_t id) { } std::shared_ptr in = _in_cube->read_chunk(id); + + // propagate chunk status + + std::shared_ptr out = std::make_shared(); + out->set_status(in->status()); // propagate chunk status if (in->empty()) { - return std::make_shared(); + return out; } // Fill buffers accordingly - std::shared_ptr out = std::make_shared(); + out->size({_bands.count(), in->size()[1], in->size()[2], in->size()[3]}); out->buf(std::calloc(_bands.count() * in->size()[1] * in->size()[2] * in->size()[3], sizeof(double))); diff --git a/src/gdalcubes/src/select_time.cpp b/src/gdalcubes/src/select_time.cpp index 229a3d1..b206bb6 100644 --- a/src/gdalcubes/src/select_time.cpp +++ b/src/gdalcubes/src/select_time.cpp @@ -43,6 +43,14 @@ std::shared_ptr select_time_cube::read_chunk(chunkid_t id) { } } + // propagate chunk status + if (in_chunk->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (in_chunk->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + // for all bands if (!in_chunk->empty()) { for (uint16_t ib = 0; ib < size_btyx[0]; ++ib) { diff --git a/src/gdalcubes/src/simple_cube.cpp b/src/gdalcubes/src/simple_cube.cpp index d4c0056..b289bb5 100644 --- a/src/gdalcubes/src/simple_cube.cpp +++ b/src/gdalcubes/src/simple_cube.cpp @@ -30,7 +30,7 @@ namespace gdalcubes { simple_cube::simple_cube(std::vector files, std::vector datetime_values, std::vector bands, - std::vector band_names, double dx, double dy) : cube(), _in_files(files), _in_datetime(datetime_values), _in_bands(bands), _in_band_names(band_names), _in_dx(dx), _in_dy(dy), _orig_bands(), _band_selection() { + std::vector band_names, double dx, double dy) : cube(), _in_files(files), _in_datetime(datetime_values), _in_bands(bands), _in_band_names(band_names), _in_dx(dx), _in_dy(dy), _strict(true), _orig_bands(), _band_selection() { if (files.size() != datetime_values.size()) { GCBS_ERROR("Number of files differs from number of provided datetime_values values"); throw std::string("Number of files differs from number of provided datetime_values values"); @@ -231,7 +231,7 @@ std::shared_ptr simple_cube::read_chunk(chunkid_t id) { std::shared_ptr out = std::make_shared(); if (id >= count_chunks()) { // chunk is outside of the cube, we don't need to read anything. - GCBS_WARN("Chunk id " + std::to_string(id) + " is out of range"); + GCBS_DEBUG("Chunk id " + std::to_string(id) + " is out of range"); return out; } @@ -261,6 +261,8 @@ std::shared_ptr simple_cube::read_chunk(chunkid_t id) { // + + uint32_t count_success = 0; // count successful image reads for (uint32_t it = 0; it < size_btyx[1]; ++it) { datetime dt = _st_ref->datetime_at_index(climits.low[0] + it); auto iter = _index.find(dt); @@ -278,7 +280,15 @@ std::shared_ptr simple_cube::read_chunk(chunkid_t id) { GDALDataset *dataset = (GDALDataset *)GDALOpen(gdal_file.c_str(), GA_ReadOnly); if (!dataset) { - GCBS_DEBUG("GDAL failed to open '" + gdal_file + "'"); + GCBS_WARN("GDAL could not open '" + gdal_file + "': ERROR no " + std::to_string(CPLGetLastErrorNo()) + ":" + CPLGetLastErrorMsg()); + if (_strict) { + GDALClose(dataset); + out = std::make_shared(); + out->set_status(chunk_data::chunk_status::ERROR); + return out; + } + GCBS_WARN("Dataset '" + gdal_file + "' will be ignored."); + out->set_status(chunk_data::chunk_status::INCOMPLETE); continue; } @@ -337,11 +347,24 @@ std::shared_ptr simple_cube::read_chunk(chunkid_t id) { it * size_btyx[2] * size_btyx[3], size_btyx[3], size_btyx[2], GDT_Float64, 0, 0, NULL); if (res != CE_None) { - GCBS_WARN("RasterIO (read) failed for " + gdal_file); + GCBS_WARN("RasterIO (read) failed for '" + gdal_file + "': ERROR no " + std::to_string(CPLGetLastErrorNo()) + ":" + CPLGetLastErrorMsg()); + if (_strict) { + GDALClose(dataset); + out = std::make_shared(); + out->set_status(chunk_data::chunk_status::ERROR); + return out; + } + GCBS_WARN("Dataset '" + gdal_file + "' will be ignored."); + out->set_status(chunk_data::chunk_status::INCOMPLETE); + continue; } } GDALClose(dataset); } + count_success++; + } + if (out->status() == chunk_data::chunk_status::INCOMPLETE && count_success == 0) { + out->set_status(chunk_data::chunk_status::ERROR); } return out; } diff --git a/src/gdalcubes/src/simple_cube.h b/src/gdalcubes/src/simple_cube.h index a6429d3..97d5e30 100644 --- a/src/gdalcubes/src/simple_cube.h +++ b/src/gdalcubes/src/simple_cube.h @@ -98,6 +98,10 @@ class simple_cube : public cube { _chunk_size = {t, y, x}; } + void set_strict(bool s) { + _strict = s; + } + json11::Json make_constructible_json() override { json11::Json::object out; out["cube_type"] = "simple_cube"; @@ -133,6 +137,7 @@ class simple_cube : public cube { out["dx"] = _in_dx; out["dy"] = _in_dy; + out["strict"] = _strict; return out; } @@ -145,6 +150,7 @@ class simple_cube : public cube { std::vector _in_band_names; double _in_dx; double _in_dy; + bool _strict; std::map>> _index; band_collection _orig_bands; diff --git a/src/gdalcubes/src/slice_space.cpp b/src/gdalcubes/src/slice_space.cpp index 7f188d3..c0a7883 100644 --- a/src/gdalcubes/src/slice_space.cpp +++ b/src/gdalcubes/src/slice_space.cpp @@ -18,6 +18,7 @@ std::shared_ptr slice_space_cube::read_chunk(chunkid_t id) { chunkid_t input_chunk_id = _in_cube->chunk_id_from_coords(input_chunk_coords); std::shared_ptr in_chunk = _in_cube->read_chunk(input_chunk_id); + out->set_status(in_chunk->status()); // propagate chunk status if (in_chunk && !in_chunk->empty()) { out->size(size_btyx); // Fill buffers accordingly diff --git a/src/gdalcubes/src/slice_time.cpp b/src/gdalcubes/src/slice_time.cpp index ca44495..75500b6 100644 --- a/src/gdalcubes/src/slice_time.cpp +++ b/src/gdalcubes/src/slice_time.cpp @@ -18,6 +18,7 @@ std::shared_ptr slice_time_cube::read_chunk(chunkid_t id) { chunkid_t input_chunk_id = _in_cube->chunk_id_from_coords(input_chunk_coords); std::shared_ptr in_chunk = _in_cube->read_chunk(input_chunk_id); + out->set_status(in_chunk->status()); // propagate chunk status if (in_chunk && !in_chunk->empty()) { out->size(size_btyx); // Fill buffers accordingly diff --git a/src/gdalcubes/src/stream_apply_pixel.cpp b/src/gdalcubes/src/stream_apply_pixel.cpp index 840d73e..9747176 100644 --- a/src/gdalcubes/src/stream_apply_pixel.cpp +++ b/src/gdalcubes/src/stream_apply_pixel.cpp @@ -14,6 +14,8 @@ std::shared_ptr stream_apply_pixel_cube::read_chunk(chunkid_t id) { std::shared_ptr inbuf = _in_cube->read_chunk(id); + + out->set_status(inbuf->status()); // propagate chunk status // check whether input chunk is empty and if yes, avoid computations if (inbuf->empty()) { return out; diff --git a/src/gdalcubes/src/stream_apply_time.cpp b/src/gdalcubes/src/stream_apply_time.cpp index 394bc04..cbe6b46 100644 --- a/src/gdalcubes/src/stream_apply_time.cpp +++ b/src/gdalcubes/src/stream_apply_time.cpp @@ -30,6 +30,13 @@ std::shared_ptr stream_apply_time_cube::read_chunk(chunkid_t id) { for (chunkid_t i = id; i < _in_cube->count_chunks(); i += _in_cube->count_chunks_x() * _in_cube->count_chunks_y()) { std::shared_ptr x = _in_cube->read_chunk(i); + // propagate chunk status + if (x->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (x->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } if (!x->empty()) { if (!initialized) { // Fill buffers with NAN @@ -65,7 +72,10 @@ std::shared_ptr stream_apply_time_cube::read_chunk(chunkid_t id) { } // check if inbuf is completely empty and if yes, avoid streaming at all and return empty chunk if (empty) { - return std::make_shared(); + auto s = out->status(); + out = std::make_shared(); + out->set_status(s); + return out; } // generate in and out filename std::string f_in = filesystem::join(config::instance()->get_streaming_dir(), diff --git a/src/gdalcubes/src/stream_reduce_space.cpp b/src/gdalcubes/src/stream_reduce_space.cpp index 7c0f167..aff2720 100644 --- a/src/gdalcubes/src/stream_reduce_space.cpp +++ b/src/gdalcubes/src/stream_reduce_space.cpp @@ -32,6 +32,13 @@ std::shared_ptr stream_reduce_space_cube::read_chunk(chunkid_t id) { uint32_t in_chunk_id = id * nchunks_in_space + i; auto in_chunk_coords = _in_cube->chunk_coords_from_id(in_chunk_id); std::shared_ptr x = _in_cube->read_chunk(in_chunk_id); + // propagate chunk status + if (x->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (x->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } if (!x->empty()) { if (!initialized) { @@ -70,7 +77,10 @@ std::shared_ptr stream_reduce_space_cube::read_chunk(chunkid_t id) { } // check if inbuf is completely empty and if yes, avoid streaming at all and return empty chunk if (empty) { - return std::make_shared(); + auto s = out->status(); + out = std::make_shared(); + out->set_status(s); + return out; } // generate in and out filename diff --git a/src/gdalcubes/src/stream_reduce_time.cpp b/src/gdalcubes/src/stream_reduce_time.cpp index eb10b96..824f1ba 100644 --- a/src/gdalcubes/src/stream_reduce_time.cpp +++ b/src/gdalcubes/src/stream_reduce_time.cpp @@ -30,6 +30,14 @@ std::shared_ptr stream_reduce_time_cube::read_chunk(chunkid_t id) { for (chunkid_t i = id; i < _in_cube->count_chunks(); i += _in_cube->count_chunks_x() * _in_cube->count_chunks_y()) { std::shared_ptr x = _in_cube->read_chunk(i); + // propagate chunk status + if (x->status() == chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::ERROR); + } + else if (x->status() == chunk_data::chunk_status::INCOMPLETE && out->status() != chunk_data::chunk_status::ERROR) { + out->set_status(chunk_data::chunk_status::INCOMPLETE); + } + if (!x->empty()) { if (!initialized) { @@ -66,7 +74,10 @@ std::shared_ptr stream_reduce_time_cube::read_chunk(chunkid_t id) { } // check if inbuf is completely empty and if yes, avoid streaming at all and return empty chunk if (empty) { - return std::make_shared(); + auto s = out->status(); + out = std::make_shared(); + out->set_status(s); + return out; } // generate in and out filename diff --git a/src/gdalcubes/src/swarm.cpp b/src/gdalcubes/src/swarm.cpp deleted file mode 100644 index 0397362..0000000 --- a/src/gdalcubes/src/swarm.cpp +++ /dev/null @@ -1,301 +0,0 @@ -/* - MIT License - - Copyright (c) 2019 Marius Appel - - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be included in all - copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ - -#ifndef GDALCUBES_NO_SWARM - -#include "swarm.h" - -#include -#include - -namespace gdalcubes { - -size_t post_file_read_callback(char *buffer, size_t size, size_t nitems, void *userdata) { - std::ifstream *is = ((std::ifstream *)userdata); - is->read(buffer, size * nitems); // or use readsome() ? - return is->gcount(); -} - -std::shared_ptr gdalcubes_swarm::from_txtfile(std::string path) { - std::ifstream file(path); - std::vector urllist; - std::string line; - while (std::getline(file, line)) { - // TODO: check if valid URL - if (!line.empty()) urllist.push_back(line); - } - - file.close(); - return gdalcubes_swarm::from_urls(urllist); -} - -void gdalcubes_swarm::post_file(std::string path, uint16_t server_index) { - if (_server_handles[server_index]) { - curl_easy_reset(_server_handles[server_index]); - - // Step 1: send a HEAD HTTP request to check whether the file already exists on the server - curl_easy_setopt(_server_handles[server_index], CURLOPT_URL, (_server_uris[server_index] + "/file" + "?name=" + path + "&size=" + std::to_string(filesystem::file_size(path))).c_str()); - curl_easy_setopt(_server_handles[server_index], CURLOPT_CUSTOMREQUEST, "HEAD"); - curl_easy_setopt(_server_handles[server_index], CURLOPT_VERBOSE, config::instance()->get_swarm_curl_verbose() ? 1L : 0L); - - CURLcode res = curl_easy_perform(_server_handles[server_index]); - long response_code; - if (res != CURLE_OK) { - throw std::string("ERROR in gdalcubes_swarm::post_file(): HEAD /file?name='" + path + "' to '" + _server_uris[server_index] + "' failed"); - } else { - curl_easy_getinfo(_server_handles[server_index], CURLINFO_RESPONSE_CODE, &response_code); - } - - if (response_code == 200) { - // already exists on server - return; - } else if (response_code == 204 || response_code == 409) { - // file does not exist or has different size - - // Step 2: if not upload - - std::ifstream is(path, std::ifstream::in | std::ifstream::binary); - if (!is.is_open()) { - throw std::string("ERROR in gdalcubes_swarm::push_execution_context(): cannot open file '" + path + "'"); - } - - curl_easy_setopt(_server_handles[server_index], CURLOPT_URL, (_server_uris[server_index] + "/file" + "?name=" + path).c_str()); - - //struct curl_slist *header = NULL; - //header = curl_slist_append(header, "Content-Type: application/octet-stream"); - // curl_easy_setopt(curl, CURLOPT_HTTPHEADER, header); - - curl_easy_setopt(_server_handles[server_index], CURLOPT_UPLOAD, 1L); - curl_easy_setopt(_server_handles[server_index], CURLOPT_POST, 1L); - curl_easy_setopt(_server_handles[server_index], CURLOPT_CUSTOMREQUEST, "POST"); // for whatever reason all requests are PUT if not set - - // This will most likely work only as long as curl_easy_perform is synchronous, otherwise &is might become invalid - curl_easy_setopt(_server_handles[server_index], CURLOPT_READDATA, &is); - curl_easy_setopt(_server_handles[server_index], CURLOPT_READFUNCTION, &post_file_read_callback); - - curl_easy_setopt(_server_handles[server_index], CURLOPT_INFILESIZE_LARGE, (curl_off_t)filesystem::file_size(path)); - - curl_easy_setopt(_server_handles[server_index], CURLOPT_VERBOSE, config::instance()->get_swarm_curl_verbose() ? 1L : 0L); - - CURLcode res = curl_easy_perform(_server_handles[server_index]); - if (res != CURLE_OK) { - if (is.is_open()) is.close(); - throw std::string("ERROR in gdalcubes_swarm::post_file(): uploading '" + path + "' to '" + _server_uris[server_index] + "' failed"); - } - if (is.is_open()) is.close(); - } else { - throw std::string("ERROR in gdalcubes_swarm::post_file(): HEAD /file?name='" + path + "' to '" + _server_uris[server_index] + "' returned HTTP code " + std::to_string(response_code)); - } - } -} - -void gdalcubes_swarm::push_execution_context(bool recursive) { - std::string p = filesystem::get_working_dir(); - - std::vector file_list; - if (recursive) { - filesystem::iterate_directory_recursive(p, [&file_list](const std::string &p) { - if (filesystem::is_regular_file(p)) { - file_list.push_back(p); - } - }); - } else { - filesystem::iterate_directory(p, [&file_list](const std::string &p) { - if (filesystem::is_regular_file(p)) { - file_list.push_back(p); - } - }); - } - - // TODO: parallel upload - for (auto it_f = file_list.begin(); it_f != file_list.end(); ++it_f) { - for (uint16_t is = 0; is < _server_uris.size(); ++is) { - post_file(*it_f, is); - } - } -} - -size_t post_cube_read_callback(char *buffer, size_t size, size_t nitems, void *userdata) { - std::istringstream *is = (std::istringstream *)userdata; - is->readsome(buffer, size * nitems); // or use readsome() ? - return is->gcount(); -} - -size_t post_cube_write_callback(char *buffer, size_t size, size_t nitems, void *userdata) { - std::vector *x = (std::vector *)userdata; - x->insert(x->end(), buffer, buffer + (size * nitems)); - return size * nitems; -} - -uint32_t gdalcubes_swarm::post_cube(std::string json, uint16_t server_index) { - std::istringstream is(json); - if (_server_handles[server_index]) { - curl_easy_reset(_server_handles[server_index]); - curl_easy_setopt(_server_handles[server_index], CURLOPT_URL, (_server_uris[server_index] + "/cube").c_str()); - - // TODO: Disable Expect: 100-Continue header? - - curl_easy_setopt(_server_handles[server_index], CURLOPT_UPLOAD, 1L); - curl_easy_setopt(_server_handles[server_index], CURLOPT_POST, 1L); - curl_easy_setopt(_server_handles[server_index], CURLOPT_CUSTOMREQUEST, "POST"); // for whatever reason all requests are PUT if not set - - struct curl_slist *header = NULL; - header = curl_slist_append(header, "Content-Type: application/json"); - header = curl_slist_append(header, "Expect:"); - curl_easy_setopt(_server_handles[server_index], CURLOPT_HTTPHEADER, header); - - curl_easy_setopt(_server_handles[server_index], CURLOPT_READFUNCTION, &post_cube_read_callback); - - // This will most likely work only as long as curl_easy_perform is synchronous, otherwise &is might become invalid - curl_easy_setopt(_server_handles[server_index], CURLOPT_READDATA, &is); - - //curl_easy_setopt(_server_handles[server_index], CURLOPT_INFILESIZE, json.size() + 1); - - curl_easy_setopt(_server_handles[server_index], CURLOPT_WRITEFUNCTION, &post_cube_write_callback); - - // This will most likely work only as long as curl_easy_perform is synchronous, otherwise &response_body_bytes might become invalid - std::vector response_body_bytes; - curl_easy_setopt(_server_handles[server_index], CURLOPT_WRITEDATA, &response_body_bytes); - - curl_easy_setopt(_server_handles[server_index], CURLOPT_VERBOSE, config::instance()->get_swarm_curl_verbose() ? 1L : 0L); - - CURLcode res = curl_easy_perform(_server_handles[server_index]); - if (res != CURLE_OK) { - throw std::string("ERROR in gdalcubes_swarm::post_cube(): POST /cube to '" + _server_uris[server_index] + "' failed"); - } - std::string response_body(response_body_bytes.begin(), response_body_bytes.end()); - return std::stoi(response_body); - } else { - throw std::string("ERROR in gdalcubes_swarm::post_cube(): no connection with given server index available"); - } -} - -void gdalcubes_swarm::push_cube(std::shared_ptr c) { - _cube = c; // Does this require a mutex? - std::string json = c->make_constructible_json().dump(); - - _cube_ids.clear(); - - // TODO: parallelize requests - for (uint16_t is = 0; is < _server_uris.size(); ++is) { - _cube_ids.push_back(post_cube(json, is)); - } -} - -void gdalcubes_swarm::post_start(uint32_t chunk_id, uint16_t server_index) { - if (_server_handles[server_index]) { - // TODO: URLencode? - curl_easy_reset(_server_handles[server_index]); - curl_easy_setopt(_server_handles[server_index], CURLOPT_URL, (_server_uris[server_index] + "/cube/" + std::to_string(_cube_ids[server_index]) + "/" + std::to_string(chunk_id) + "/start").c_str()); - curl_easy_setopt(_server_handles[server_index], CURLOPT_POST, 1L); - curl_easy_setopt(_server_handles[server_index], CURLOPT_CUSTOMREQUEST, "POST"); - curl_easy_setopt(_server_handles[server_index], CURLOPT_VERBOSE, config::instance()->get_swarm_curl_verbose() ? 1L : 0L); - - CURLcode res = curl_easy_perform(_server_handles[server_index]); - if (res != CURLE_OK) { - throw std::string("ERROR in gdalcubes_swarm::post_start(): POST /cube/{cube_id}/{chunk_id}/start to '" + _server_uris[server_index] + "' failed"); - } - return; - } -} - -size_t get_download_callback(char *buffer, size_t size, size_t nitems, void *userdata) { - std::vector *x = (std::vector *)userdata; - x->insert(x->end(), buffer, buffer + (size * nitems)); - return size * nitems; -} - -std::shared_ptr gdalcubes_swarm::get_download(uint32_t chunk_id, uint16_t server_index) { - if (_server_handles[server_index]) { - // TODO: URLencode? - curl_easy_reset(_server_handles[server_index]); - curl_easy_setopt(_server_handles[server_index], CURLOPT_URL, (_server_uris[server_index] + "/cube/" + std::to_string(_cube_ids[server_index]) + "/" + std::to_string(chunk_id) + "/download").c_str()); - curl_easy_setopt(_server_handles[server_index], CURLOPT_HTTPGET, 1L); - curl_easy_setopt(_server_handles[server_index], CURLOPT_CUSTOMREQUEST, "GET"); - curl_easy_setopt(_server_handles[server_index], CURLOPT_WRITEFUNCTION, &get_download_callback); - - // This will most likely work only as long as curl_easy_perform is synchronous, otherwise &response_body_bytes might become invalid - std::vector response_body_bytes; - curl_easy_setopt(_server_handles[server_index], CURLOPT_WRITEDATA, &response_body_bytes); - - curl_easy_setopt(_server_handles[server_index], CURLOPT_VERBOSE, config::instance()->get_swarm_curl_verbose() ? 1L : 0L); - - CURLcode res = curl_easy_perform(_server_handles[server_index]); - if (res != CURLE_OK) { - throw std::string("ERROR in gdalcubes_swarm::get_download(): GET /cube/{cube_id}/{chunk_id}/download to '" + _server_uris[server_index] + "' failed"); - } - - // construct chunk_data from byte vector - std::shared_ptr out = std::make_shared(); - std::array size = {((uint32_t *)response_body_bytes.data())[0], ((uint32_t *)response_body_bytes.data())[1], ((uint32_t *)response_body_bytes.data())[2], ((uint32_t *)response_body_bytes.data())[3]}; - out->size(size); - if (size[0] * size[1] * size[2] * size[3] > 0) { - out->buf(std::malloc(out->size()[0] * out->size()[1] * out->size()[2] * out->size()[3] * sizeof(double))); - std::copy(response_body_bytes.begin() + sizeof(std::array), response_body_bytes.end(), - (char *)out->buf()); - } - - return out; - } else { - throw std::string("ERROR in gdalcubes_swarm::get_download(): no connection with given server index available"); - } -} - -void gdalcubes_swarm::apply(std::shared_ptr c, std::function, std::mutex &)> f) { - uint32_t nthreads = 1; - // try whether default chunk processor is multithread and read number of threads if successful - if (std::dynamic_pointer_cast(config::instance()->get_default_chunk_processor())) { - nthreads = std::dynamic_pointer_cast(config::instance()->get_default_chunk_processor())->get_threads(); - } - - push_execution_context(false); - push_cube(c); - - // TODO: what if servers fail? - std::map> chunk_distr; - for (uint32_t i = 0; i < c->count_chunks(); ++i) { - post_start(i, i % _server_uris.size()); // start computations on remote servers - chunk_distr[i % _server_uris.size()].push_back(i); - } - std::mutex mutex; - std::vector workers; - for (uint16_t it = 0; it < nthreads; ++it) { - workers.push_back(std::thread([this, c, &chunk_distr, &f, it, &mutex](void) { - // One remote server is processed by a single thread for now, changing this will require to using the multi interface of curl - for (uint32_t iserver = it; iserver < _server_uris.size(); iserver += _server_uris.size()) { - for (uint32_t ichunk = 0; ichunk < chunk_distr[iserver].size(); ++ichunk) { - std::shared_ptr dat = get_download(chunk_distr[iserver][ichunk], iserver); - f(chunk_distr[iserver][ichunk], dat, mutex); - } - } - })); - } - for (uint16_t it = 0; it < nthreads; ++it) { - workers[it].join(); - } -} - -} // namespace gdalcubes - -#endif // GDALCUBES_NO_SWARM diff --git a/src/gdalcubes/src/swarm.h b/src/gdalcubes/src/swarm.h deleted file mode 100644 index 805aad6..0000000 --- a/src/gdalcubes/src/swarm.h +++ /dev/null @@ -1,98 +0,0 @@ -/* - MIT License - - Copyright (c) 2019 Marius Appel - - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be included in all - copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ - -#ifndef GDALCUBES_NO_SWARM - -#ifndef SWARM_H -#define SWARM_H - -#include - -#include "cube.h" - -namespace gdalcubes { - -/** - * @brief Chunk processor implementation for distributed processing by connecting to gdalcubes_server instances - * - * This class connects the several gdalcubes_server instances in order to distribute read_chunk() operations - * - * @todo implement add / remove method for workers - */ -class gdalcubes_swarm : public chunk_processor { - public: - gdalcubes_swarm(std::vector urls) : _cube(nullptr), _cube_ids(), _server_handles(), _server_uris(urls), _nthreads(1) { - for (uint16_t i = 0; i < _server_uris.size(); ++i) - _server_handles.push_back(curl_easy_init()); - } - - ~gdalcubes_swarm() { - for (uint16_t i = 0; i < _server_uris.size(); ++i) - curl_easy_cleanup(_server_handles[i]); - } - - inline static std::shared_ptr from_urls(std::vector urls) { return std::make_shared(urls); } - - // Create from txt file where each line is a uri - static std::shared_ptr from_txtfile(std::string path); - - // upload execution context to all servers - void push_execution_context(bool recursive = false); - - // create cube on all servers - void push_cube(std::shared_ptr c); - - // Mimic cube::apply with distributed calls to cube::read_chunk() - void apply(std::shared_ptr c, std::function, std::mutex &)> f) override; - - /** - * @copydoc chunk_processor::max_threads - */ - uint32_t max_threads() override { - return _nthreads; - } - inline void set_threads(uint16_t threads) { _nthreads = threads; } - - private: - void post_file(std::string path, uint16_t server_index); - - void post_start(uint32_t chunk_id, uint16_t server_index); - std::shared_ptr get_download(uint32_t chunk_id, uint16_t server_index); - - uint32_t post_cube(std::string json, uint16_t server_index); - - std::shared_ptr _cube; - std::vector _cube_ids; // IDs of the cube for each server - - std::vector _server_handles; - std::vector _server_uris; - - uint16_t _nthreads; -}; - -} // namespace gdalcubes - -#endif //SWARM_H - -#endif // GDALCUBES_NO_SWARM diff --git a/src/gdalcubes/src/vector_queries.cpp b/src/gdalcubes/src/vector_queries.cpp deleted file mode 100644 index 11e9759..0000000 --- a/src/gdalcubes/src/vector_queries.cpp +++ /dev/null @@ -1,1170 +0,0 @@ -/* - MIT License - - Copyright (c) 2019 Marius Appel - - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be included in all - copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ - -#include "vector_queries.h" - -#include -#include - -#include -#include -#include - -namespace gdalcubes { - -std::vector> vector_queries::query_points(std::shared_ptr cube, std::vector x, - std::vector y, std::vector t, - std::string srs) { - // make sure that x, y, and t have same size - if (x.size() != y.size() || y.size() != t.size()) { - GCBS_ERROR("Point coordinate vectors x, y, t must have identical length"); - throw std::string("Point coordinate vectors x, y, t must have identical length"); - } - - if (x.empty()) { - GCBS_ERROR("Point coordinate vectors x, y, t must have length > 0"); - throw std::string("Point coordinate vectors x, y, t must have length > 0"); - } - - if (!cube) { - GCBS_ERROR("Invalid data cube pointer"); - throw std::string("Invalid data cube pointer"); - } - - uint32_t nthreads = config::instance()->get_default_chunk_processor()->max_threads(); - - std::shared_ptr prg = config::instance()->get_default_progress_bar()->get(); - prg->set(0); // explicitly set to zero to show progress bar immediately - - // coordinate transformation - if (cube->st_reference()->srs() != srs) { - OGRSpatialReference srs_in; - OGRSpatialReference srs_out; - srs_in.SetFromUserInput(srs.c_str()); - srs_out.SetFromUserInput(cube->st_reference()->srs().c_str()); - - if (!srs_in.IsSame(&srs_out)) { - std::vector workers_transform; - uint32_t n = (uint32_t)std::ceil(double(x.size()) / double(nthreads)); // points per thread - - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers_transform.push_back(std::thread([&cube, &srs, &srs_in, &srs_out, &x, &y, ithread, n](void) { - OGRCoordinateTransformation *coord_transform = OGRCreateCoordinateTransformation(&srs_in, &srs_out); - - int begin = ithread * n; - int end = std::min(uint32_t(ithread * n + n), uint32_t(x.size())); - int count = end - begin; - - // change coordinates in place, should be safe because vectors don't change their sizes - if (count > 0) { - if (coord_transform == nullptr || !coord_transform->Transform(count, x.data() + begin, y.data() + begin)) { - throw std::string("ERROR: coordinate transformation failed (from " + - cube->st_reference()->srs() + " to " + srs + ")."); - } - } - OCTDestroyCoordinateTransformation(coord_transform); - })); - } - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers_transform[ithread].join(); - } - } - } - - // TODO: possible without additional copy? - std::vector it; // array indexes - // ix.resize(x.size()); - // iy.resize(x.size()); - it.resize(x.size()); - - std::map> chunk_index; - - std::vector workers_preprocess; - std::mutex mtx; - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers_preprocess.push_back(std::thread([&mtx, &cube, &x, &y, &t, &it, &chunk_index, ithread, nthreads](void) { - for (uint32_t i = ithread; i < x.size(); i += nthreads) { - coords_st st; - - st.s.x = x[i]; - st.s.y = y[i]; - - // array coordinates - x[i] = (x[i] - cube->st_reference()->left()) / cube->st_reference()->dx(); - //iy.push_back(cube->st_reference()->ny() - 1 - ((y[i] - cube->st_reference()->bottom()) / cube->st_reference()->dy())); // top 0 - y[i] = (y[i] - cube->st_reference()->bottom()) / cube->st_reference()->dy(); - - datetime dt = datetime::from_string(t[i]); - if (dt.unit() > cube->st_reference()->dt().dt_unit) { - dt.unit(cube->st_reference()->dt().dt_unit); - GCBS_WARN("date / time of query point has coarser granularity than the data cube; converting '" + t[i] + "' -> '" + dt.to_string() + "'"); - } else { - dt.unit(cube->st_reference()->dt().dt_unit); - } - it[i] = cube->st_reference()->index_at_datetime(dt); - - if (it[i] < 0 || it[i] >= cube->size_t() || - x[i] < 0 || x[i] >= cube->size_x() || - y[i] < 0 || y[i] >= cube->size_y()) { // if point is outside of the cube - continue; - } - st.t = dt; - chunkid_t c = cube->find_chunk_that_contains(st); - - mtx.lock(); - chunk_index[c].push_back(i); - mtx.unlock(); - } - })); - } - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers_preprocess[ithread].join(); - } - - std::vector> out; - out.resize(cube->bands().count()); - for (uint16_t ib = 0; ib < out.size(); ++ib) { - out[ib].resize(x.size(), NAN); - } - - std::vector chunks; // vector of keys in chunk_index - for (auto iter = chunk_index.begin(); iter != chunk_index.end(); ++iter) { - chunks.push_back(iter->first); - } - - std::vector workers; - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers.push_back(std::thread([&prg, &cube, &out, &chunk_index, &chunks, &x, &it, &y, ithread, nthreads](void) { - for (uint32_t ic = ithread; ic < chunks.size(); ic += nthreads) { - try { - if (chunks[ic] < cube->count_chunks()) { // if chunk exists - std::shared_ptr dat = cube->read_chunk(chunks[ic]); - if (!dat->empty()) { // if chunk is not empty - // iterate over all query points within the current chunk - for (uint32_t i = 0; i < chunk_index[chunks[ic]].size(); ++i) { - double ixc = x[chunk_index[chunks[ic]][i]]; - double iyc = y[chunk_index[chunks[ic]][i]]; - double itc = it[chunk_index[chunks[ic]][i]]; - - int iix = ((int)std::floor(ixc)) % cube->chunk_size()[2]; - int iiy = dat->size()[2] - 1 - (((int)std::floor(iyc)) % cube->chunk_size()[1]); - int iit = ((int)std::floor(itc)) % cube->chunk_size()[0]; - - // check to prevent out of bounds faults - if (iix < 0 || uint32_t(iix) >= dat->size()[3]) continue; - if (iiy < 0 || uint32_t(iiy) >= dat->size()[2]) continue; - if (iit < 0 || uint32_t(iit) >= dat->size()[1]) continue; - - for (uint16_t ib = 0; ib < out.size(); ++ib) { - out[ib][chunk_index[chunks[ic]][i]] = ((double *)dat->buf())[ib * dat->size()[1] * dat->size()[2] * dat->size()[3] + iit * dat->size()[2] * dat->size()[3] + iiy * dat->size()[3] + iix]; - } - } - } - } - prg->increment((double)1 / (double)chunks.size()); - } catch (std::string s) { - GCBS_ERROR(s); - continue; - } catch (...) { - GCBS_ERROR("unexpected exception while processing chunk " + std::to_string(chunks[ic])); - continue; - } - } - })); - } - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers[ithread].join(); - } - prg->finalize(); - - return out; -} - -std::vector>> vector_queries::query_timeseries(std::shared_ptr cube, - std::vector x, - std::vector y, - std::string srs) { - if (x.size() != y.size()) { - GCBS_ERROR("Point coordinate vectors x, y must have identical length"); - throw std::string("Point coordinate vectors x, y must have identical length"); - } - - if (x.empty()) { - GCBS_ERROR("Point coordinate vectors x, y must have length > 0"); - throw std::string("Point coordinate vectors x, y must have length > 0"); - } - - if (!cube) { - GCBS_ERROR("Invalid data cube pointer"); - throw std::string("Invalid data cube pointer"); - } - - uint32_t nthreads = config::instance()->get_default_chunk_processor()->max_threads(); - - std::shared_ptr prg = config::instance()->get_default_progress_bar()->get(); - prg->set(0); // explicitly set to zero to show progress bar immediately - - // coordinate transformation - if (cube->st_reference()->srs() != srs) { - OGRSpatialReference srs_in; - OGRSpatialReference srs_out; - srs_in.SetFromUserInput(cube->st_reference()->srs().c_str()); - srs_out.SetFromUserInput(srs.c_str()); - - if (!srs_in.IsSame(&srs_out)) { - std::vector workers_transform; - uint32_t n = (uint32_t)std::ceil(double(x.size()) / double(nthreads)); // points per thread - - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers_transform.push_back(std::thread([&cube, &srs, &srs_in, &srs_out, &x, &y, ithread, n](void) { - OGRCoordinateTransformation *coord_transform = OGRCreateCoordinateTransformation(&srs_in, &srs_out); - - int begin = ithread * n; - int end = std::min(uint32_t(ithread * n + n), uint32_t(x.size())); - int count = end - begin; - - // change coordinates in place, should be safe because vectors don't change their sizes - if (count > 0) { - if (coord_transform == nullptr || !coord_transform->Transform(count, x.data() + begin, y.data() + begin)) { - throw std::string("ERROR: coordinate transformation failed (from " + - cube->st_reference()->srs() + " to " + srs + ")."); - } - } - OCTDestroyCoordinateTransformation(coord_transform); - })); - } - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers_transform[ithread].join(); - } - } - } - - // TODO: possible without additional copy? - std::vector ipoints; // array indexes - // ix.resize(x.size()); - // iy.resize(x.size()); - ipoints.resize(x.size()); - - std::map> chunk_index; - std::vector workers_preprocess; - std::mutex mtx; - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers_preprocess.push_back(std::thread([&mtx, &cube, &x, &y, &chunk_index, ithread, nthreads](void) { - for (uint32_t i = ithread; i < x.size(); i += nthreads) { - coords_st st; - - st.s.x = x[i]; - st.s.y = y[i]; - - // array coordinates - double xarr = (x[i] - cube->st_reference()->left()) / cube->st_reference()->dx(); - //iy.push_back(cube->st_reference()->ny() - 1 - ((y[i] - cube->st_reference()->bottom()) / cube->st_reference()->dy())); // top 0 - double yarr = (y[i] - cube->st_reference()->bottom()) / cube->st_reference()->dy(); - - if (xarr < 0 || xarr >= cube->size_x() || - yarr < 0 || yarr >= cube->size_y()) { // if point is outside of the cube - continue; - } - uint32_t cx = xarr / cube->chunk_size()[2]; - uint32_t cy = yarr / cube->chunk_size()[1]; - chunkid_t c = cube->chunk_id_from_coords({0, cy, cx}); - // st.t = cube->st_reference()->t0(); - // chunkid_t c = cube->find_chunk_that_contains(st); - // - mtx.lock(); - chunk_index[c].push_back(i); - mtx.unlock(); - } - })); - } - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers_preprocess[ithread].join(); - } - - std::vector>> out; - out.resize(cube->bands().count()); - for (uint16_t ib = 0; ib < out.size(); ++ib) { - out[ib].resize(cube->size_t()); - for (uint32_t it = 0; it < out[ib].size(); ++it) { - out[ib][it].resize(x.size(), NAN); - } - } - - std::vector chunks; // vector of keys in chunk_index - for (auto iter = chunk_index.begin(); iter != chunk_index.end(); ++iter) { - chunks.push_back(iter->first); - } - - std::vector workers; - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers.push_back(std::thread([&prg, &cube, &out, &chunk_index, &chunks, &x, &y, ithread, nthreads](void) { - for (uint32_t ic = ithread; ic < chunks.size(); ic += nthreads) { - try { - for (uint32_t ct = 0; ct < cube->count_chunks_t(); ct++) { - chunkid_t cur_chunk = chunks[ic] + ct * (cube->count_chunks_x() * cube->count_chunks_y()); - if (cur_chunk < cube->count_chunks()) { // if chunk exists - uint32_t nt_in_chunk = cube->chunk_size(cur_chunk)[0]; - std::shared_ptr dat = cube->read_chunk(cur_chunk); - if (!dat->empty()) { // if chunk is not empty - // iterate over all query points within the current chunk - for (uint32_t i = 0; i < chunk_index[chunks[ic]].size(); ++i) { - double ixc = x[chunk_index[chunks[ic]][i]]; - double iyc = y[chunk_index[chunks[ic]][i]]; - - int iix = (ixc - cube->bounds_from_chunk(cur_chunk).s.left) / - cube->st_reference()->dx(); - int iiy = (cube->bounds_from_chunk(cur_chunk).s.top - iyc) / - cube->st_reference()->dy(); - - // check to prevent out of bounds faults - if (iix < 0 || uint32_t(iix) >= dat->size()[3]) continue; - if (iiy < 0 || uint32_t(iiy) >= dat->size()[2]) continue; - - for (uint16_t ib = 0; ib < out.size(); ++ib) { - for (uint32_t it = 0; it < nt_in_chunk; ++it) { - out[ib][ct * cube->chunk_size()[0] + it][chunk_index[chunks[ic]][i]] = ((double *)dat->buf())[ib * dat->size()[1] * dat->size()[2] * dat->size()[3] + it * dat->size()[2] * dat->size()[3] + iiy * dat->size()[3] + iix]; - } - } - } - } - } - prg->increment((double)1 / ((double)chunks.size() * (double)cube->count_chunks_t())); - } - } catch (std::string s) { - GCBS_ERROR(s); - continue; - } catch (...) { - GCBS_ERROR("unexpected exception while processing chunk " + std::to_string(chunks[ic])); - continue; - } - } - })); - } - for (uint32_t ithread = 0; ithread < nthreads; ++ithread) { - workers[ithread].join(); - } - prg->finalize(); - - return out; -} - -struct zonal_statistics_func { - zonal_statistics_func() : _nfeatures(0), _nt(0){}; - virtual ~zonal_statistics_func(){}; - - virtual void init(uint32_t nfeatures, uint32_t nt) { - _nfeatures = nfeatures; - _nt = nt; - }; - virtual void update(double x, uint32_t ifeature, uint32_t it) = 0; - virtual std::shared_ptr> finalize() = 0; - - protected: - uint32_t _nfeatures; - uint32_t _nt; -}; - -struct zonal_statistics_count : public zonal_statistics_func { - void init(uint32_t nfeatures, uint32_t nt) override { - zonal_statistics_func::init(nfeatures, nt); - _x = std::make_shared>(); - _x->resize(_nt * _nfeatures, 0); - } - - void update(double x, uint32_t ifeature, uint32_t it) override { - if (std::isfinite(x)) (*_x)[ifeature * _nt + it] += 1; - } - - std::shared_ptr> finalize() override { - return _x; - } - - std::shared_ptr> _x; -}; - -struct zonal_statistics_sum : public zonal_statistics_func { - void init(uint32_t nfeatures, uint32_t nt) override { - zonal_statistics_func::init(nfeatures, nt); - _x = std::make_shared>(); - _x->resize(_nt * _nfeatures, NAN); - } - - void update(double x, uint32_t ifeature, uint32_t it) override { - if (std::isfinite(x)) { - if (std::isnan((*_x)[ifeature * _nt + it])) { - (*_x)[ifeature * _nt + it] = x; - } else { - (*_x)[ifeature * _nt + it] += x; - } - } - } - - std::shared_ptr> finalize() override { - return _x; - } - - std::shared_ptr> _x; -}; - -struct zonal_statistics_prod : public zonal_statistics_func { - void init(uint32_t nfeatures, uint32_t nt) override { - zonal_statistics_func::init(nfeatures, nt); - _x = std::make_shared>(); - _x->resize(_nt * _nfeatures, NAN); - } - - void update(double x, uint32_t ifeature, uint32_t it) override { - if (std::isfinite(x)) { - if (std::isnan((*_x)[ifeature * _nt + it])) { - (*_x)[ifeature * _nt + it] = x; - } else { - (*_x)[ifeature * _nt + it] *= x; - } - } - } - - std::shared_ptr> finalize() override { - return _x; - } - - std::shared_ptr> _x; -}; - -struct zonal_statistics_mean : public zonal_statistics_func { - void init(uint32_t nfeatures, uint32_t nt) override { - zonal_statistics_func::init(nfeatures, nt); - _x = std::make_shared>(); - _x->resize(_nt * _nfeatures, 0); - _n.resize(_nt * _nfeatures, 0); - } - - void update(double x, uint32_t ifeature, uint32_t it) override { - if (std::isfinite(x)) { - (*_x)[ifeature * _nt + it] += x; - _n[ifeature * _nt + it]++; - } - } - - std::shared_ptr> finalize() override { - for (uint32_t i = 0; i < _nfeatures * _nt; ++i) { - (*_x)[i] = (*_x)[i] / double(_n[i]); - } - return _x; - } - - std::shared_ptr> _x; - std::vector _n; -}; - -struct zonal_statistics_min : public zonal_statistics_func { - void init(uint32_t nfeatures, uint32_t nt) override { - zonal_statistics_func::init(nfeatures, nt); - _x = std::make_shared>(); - _x->resize(_nt * _nfeatures, std::numeric_limits::max()); - } - void update(double x, uint32_t ifeature, uint32_t it) override { - if (std::isfinite(x)) { - (*_x)[ifeature * _nt + it] = std::min(x, (*_x)[ifeature * _nt + it]); - } - } - std::shared_ptr> finalize() override { - for (uint32_t i = 0; i < _nfeatures * _nt; ++i) { - if ((*_x)[i] == std::numeric_limits::max()) { - (*_x)[i] = NAN; - } - } - return _x; - } - std::shared_ptr> _x; -}; - -struct zonal_statistics_max : public zonal_statistics_func { - void init(uint32_t nfeatures, uint32_t nt) override { - zonal_statistics_func::init(nfeatures, nt); - _x = std::make_shared>(); - _x->resize(_nt * _nfeatures, std::numeric_limits::lowest()); - } - void update(double x, uint32_t ifeature, uint32_t it) override { - if (std::isfinite(x)) { - (*_x)[ifeature * _nt + it] = std::max(x, (*_x)[ifeature * _nt + it]); - } - } - std::shared_ptr> finalize() override { - for (uint32_t i = 0; i < _nfeatures * _nt; ++i) { - if ((*_x)[i] == std::numeric_limits::lowest()) { - (*_x)[i] = NAN; - } - } - return _x; - } - std::shared_ptr> _x; -}; - -struct zonal_statistics_median : public zonal_statistics_func { - void init(uint32_t nfeatures, uint32_t nt) override { - zonal_statistics_func::init(nfeatures, nt); - _values.resize(_nt * _nfeatures); - } - - void update(double x, uint32_t ifeature, uint32_t it) override { - if (std::isfinite(x)) { - _values[ifeature * _nt + it].push_back(x); - } - } - - std::shared_ptr> finalize() override { - std::shared_ptr> out = std::make_shared>(); - out->resize(_nt * _nfeatures); - for (uint32_t i = 0; i < _nfeatures * _nt; ++i) { - std::sort(_values[i].begin(), _values[i].end()); - if (_values[i].size() == 0) { - (*out)[i] = NAN; - } else if (_values[i].size() % 2 == 1) { - (*out)[i] = _values[i][_values[i].size() / 2]; - } else { - (*out)[i] = (_values[i][_values[i].size() / 2] + _values[i][_values[i].size() / 2 - 1]) / ((double)2); - } - } - return out; - } - - std::vector> _values; -}; - -struct zonal_statistics_var : public zonal_statistics_func { - void init(uint32_t nfeatures, uint32_t nt) override { - zonal_statistics_func::init(nfeatures, nt); - _n.resize(_nt * _nfeatures, 0); - _cur_mean.resize(_nt * _nfeatures, 0); - _cur_M2 = std::make_shared>(); - _cur_M2->resize(_nt * _nfeatures, 0); - } - - void update(double x, uint32_t ifeature, uint32_t it) override { - if (std::isfinite(x)) { - _n[ifeature * _nt + it]++; - double delta = x - _cur_mean[ifeature * _nt + it]; - _cur_mean[ifeature * _nt + it] += delta / _n[ifeature * _nt + it]; - double delta2 = x - _cur_mean[ifeature * _nt + it]; - (*_cur_M2)[ifeature * _nt + it] += delta * delta2; - } - } - - std::shared_ptr> finalize() override { - for (uint32_t i = 0; i < _nfeatures * _nt; ++i) { - if (_n[i] < 2) { - (*_cur_M2)[i] = NAN; - } else { - (*_cur_M2)[i] = (*_cur_M2)[i] / double(_n[i]); - } - } - return _cur_M2; - } - - std::vector _n; - std::vector _cur_mean; - std::shared_ptr> _cur_M2; -}; - -struct zonal_statistics_sd : public zonal_statistics_var { - std::shared_ptr> finalize() override { - for (uint32_t i = 0; i < _nfeatures * _nt; ++i) { - if (_n[i] < 2) { - (*_cur_M2)[i] = NAN; - } else { - (*_cur_M2)[i] = std::sqrt((*_cur_M2)[i] / double(_n[i])); - } - } - return _cur_M2; - } -}; - -void vector_queries::zonal_statistics(std::shared_ptr cube, std::string ogr_dataset, - std::vector> agg_band_functions, - std::string out_path, bool overwrite_if_exists, std::string ogr_layer) { - if (!OGRGeometryFactory::haveGEOS()) { - GCBS_ERROR("Missing GEOS support in GDAL installation"); - throw std::string("Missing GEOS support in GDAL installation"); - } - - if (filesystem::exists(out_path)) { - if (!overwrite_if_exists) { - GCBS_ERROR("Output file already exists; please select a different name or set argument overwrite_if_exists = true"); - throw std::string("Output file already exists; please select a different name or set argument overwrite_if_exists = true"); - } - } - if (filesystem::filename(out_path).empty()) { - GCBS_ERROR("Invalid output path, missing a filename"); - throw std::string("Invalid output path, missing a filename"); - } - - std::vector band_index; - std::vector agg_func_names; - std::vector()>> agg_func_creators; - for (uint16_t i = 0; i < agg_band_functions.size(); ++i) { - if (!cube->bands().has(agg_band_functions[i].second)) { - GCBS_WARN("Data cube has no band '" + agg_band_functions[i].second + "', statistics on this band will be ignored"); - } else { - if (agg_band_functions[i].first == "min") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_min()); }); - } else if (agg_band_functions[i].first == "max") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_max()); }); - } else if (agg_band_functions[i].first == "count") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_count()); }); - } else if (agg_band_functions[i].first == "sum") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_sum()); }); - } else if (agg_band_functions[i].first == "prod") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_prod()); }); - } else if (agg_band_functions[i].first == "mean") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_mean()); }); - } else if (agg_band_functions[i].first == "median") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_median()); }); - } else if (agg_band_functions[i].first == "sd") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_sd()); }); - } else if (agg_band_functions[i].first == "var") { - agg_func_creators.push_back([]() { return std::unique_ptr(new zonal_statistics_var()); }); - } else { - GCBS_WARN("There is no aggregation function '" + agg_band_functions[i].first + "', related summary statistics will be ignored."); - continue; - } - band_index.push_back(cube->bands().get_index(agg_band_functions[i].second)); - agg_func_names.push_back(agg_band_functions[i].first); - } - } - - if (agg_func_names.empty() || band_index.empty()) { - GCBS_ERROR("No valid summary statistics functions given"); - return; - } - - // open input OGR dataset - GDALDataset *in_ogr_dataset; - in_ogr_dataset = (GDALDataset *)GDALOpenEx(ogr_dataset.c_str(), GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, - NULL); - if (in_ogr_dataset == NULL) { - GCBS_ERROR("failed to open '" + ogr_dataset + "'"); - throw std::string("failed to open '" + ogr_dataset + "'"); - } - - OGRLayer *layer; - if (in_ogr_dataset->GetLayerCount() > 1) { - if (ogr_layer.empty()) { - GCBS_WARN("input OGR dataset has multiple layers, using the first."); - layer = in_ogr_dataset->GetLayer(0); - } else { - layer = in_ogr_dataset->GetLayerByName(ogr_layer.c_str()); - } - } else { - if (ogr_layer.empty()) { - layer = in_ogr_dataset->GetLayer(0); - } else { - layer = in_ogr_dataset->GetLayerByName(ogr_layer.c_str()); - } - } - if (layer == NULL) { - GCBS_ERROR("invalid OGR layer"); - throw std::string("invalid OGR layer"); - } - ogr_layer = layer->GetName(); - - // If layer has more than one geometry column, only the first will be used. - // Warn if there are more geometry columns - if (layer->GetLayerDefn()->GetGeomFieldCount() > 1) { - std::string geom_field_name = layer->GetLayerDefn()->GetGeomFieldDefn(0)->GetNameRef(); - GCBS_WARN("Found more than one geometry field for input features, using only the first ('" + geom_field_name + "'"); - } - - // Make sure that geometry column has type Polygon / MultiPolygon - OGRwkbGeometryType geom_type = layer->GetGeomType(); - if (geom_type != wkbPolygon && geom_type != wkbPolygonM && geom_type != wkbPolygon25D && - geom_type != wkbPolygonZM && geom_type != wkbMultiPolygon) { - GCBS_ERROR("Zonal statistics requires Polygon or MultiPolygon input geometries"); - GDALClose(in_ogr_dataset); - // TODO: do we have to clean up more things here? - throw std::string("Zonal statistics requires Polygon or MultiPolygon input geometries"); - } - - // check that cube and ogr dataset have same spatial reference system. - OGRSpatialReference srs_cube = cube->st_reference()->srs_ogr(); - srs_cube.AutoIdentifyEPSG(); - OGRSpatialReference srs_features = *(layer->GetSpatialRef()); - srs_features.AutoIdentifyEPSG(); - - bool in_ogr_was_transformed = false; - if (!srs_cube.IsSame(&srs_features)) { - // Transform features to srs of cube - - CPLStringList translate_args; - translate_args.AddString("-f"); - translate_args.AddString("GPKG"); - translate_args.AddString("-t_srs"); - translate_args.AddString(cube->st_reference()->srs().c_str()); - GDALVectorTranslateOptions *opts = GDALVectorTranslateOptionsNew(translate_args.List(), NULL); - if (opts == NULL) { - GDALVectorTranslateOptionsFree(opts); - throw std::string("ERROR in vector_queries::zonal_statistics(): cannot create ogr2ogr options."); - } - // remember filename of temporary copy with transformed input features - ogr_dataset = filesystem::join(filesystem::get_tempdir(), utils::generate_unique_filename() + ".gpkg"); - GDALDatasetH temp = GDALVectorTranslate(ogr_dataset.c_str(), NULL, 1, (GDALDatasetH *)&in_ogr_dataset, opts, NULL); - if (!temp) { - GDALClose(in_ogr_dataset); - GCBS_ERROR("Failed to transform input feature to data cube SRS"); - throw std::string("Failed to transform input feature to data cube SRS"); - } - GDALClose(in_ogr_dataset); - in_ogr_dataset = (GDALDataset *)temp; - layer = in_ogr_dataset->GetLayerByName(ogr_layer.c_str()); - GDALVectorTranslateOptionsFree(opts); - in_ogr_was_transformed = true; - } - - // Helper data structure: map: chunk index -> vector of FIDs of all features intersecting the chunk - std::map> features_in_chunk; - - // Check assumption: Input dataset has FIDs - std::string fid_column = layer->GetFIDColumn(); - if (fid_column.empty()) { - GCBS_ERROR("Input feature dataset must have FIDs"); - throw std::string("Input feature dataset must have FIDs"); - } - - if (!layer->TestCapability(OLCRandomRead)) { - GCBS_WARN("Input feature layer does not support efficient random reads; computations may take considerably longer."); - } - - // start progress bar - std::shared_ptr prg = config::instance()->get_default_progress_bar()->get(); - prg->set(0); // explicitly set to zero to show progress bar immediately - - // TODO: use SetSpatialFilterRect() ??? - - for (uint32_t cx = 0; cx < cube->count_chunks_x(); ++cx) { - for (uint32_t cy = 0; cy < cube->count_chunks_y(); ++cy) { - chunkid_t id = cube->chunk_id_from_coords({0, cy, cx}); - bounds_2d sextent = cube->bounds_from_chunk(id).s; - OGRPolygon pp; - OGRLinearRing a; - - a.addPoint(sextent.left, sextent.bottom); - a.addPoint(sextent.right, sextent.bottom); - a.addPoint(sextent.right, sextent.top); - a.addPoint(sextent.left, sextent.top); - a.addPoint(sextent.left, sextent.bottom); - pp.addRing(&a); - pp.assignSpatialReference(&srs_cube); - - // iterate over all features - layer->ResetReading(); - OGRFeature *cur_feature; - while ((cur_feature = layer->GetNextFeature()) != NULL) { - OGRGeometry *cur_geometry = cur_feature->GetGeometryRef(); - if (cur_geometry != NULL) { - if (cur_geometry->Intersects(&pp)) { - features_in_chunk[id].push_back(cur_feature->GetFID()); - } - } - OGRFeature::DestroyFeature(cur_feature); - } - } - } - - // Helper data structure: map: FID -> internal integer index - std::unordered_map FID_of_index; - std::unordered_map index_of_FID; - - uint32_t nfeatures = layer->GetFeatureCount(); - layer->ResetReading(); - OGRFeature *cur_feature; - int32_t cur_index = 0; - while ((cur_feature = layer->GetNextFeature()) != NULL) { - FID_of_index[cur_index] = cur_feature->GetFID(); - index_of_FID[cur_feature->GetFID()] = cur_index; - cur_index++; - OGRFeature::DestroyFeature(cur_feature); - } - layer->ResetReading(); - - GDALDriver *gpkg_driver = GetGDALDriverManager()->GetDriverByName("GPKG"); - if (gpkg_driver == NULL) { - GCBS_ERROR("OGR GeoPackage driver not found"); - throw std::string("OGR GeoPackage driver not found"); - } - - std::string output_file = out_path; - GDALDataset *gpkg_out = gpkg_driver->Create(output_file.c_str(), 0, 0, 0, GDT_Unknown, NULL); - if (gpkg_out == NULL) { - GCBS_ERROR("Creation of GPKG file '" + output_file + "' failed"); - throw std::string("Creation of GPKG file '" + output_file + "' failed"); - } - - // Copy geometries from input dataset - OGRLayer *geom_layer_out = gpkg_out->CreateLayer("geom", &srs_features, geom_type, NULL); - if (geom_layer_out == NULL) { - GCBS_ERROR("Failed to create output layer in '" + output_file + "'"); - throw std::string("Failed to create output layer in '" + output_file + "'"); - } - for (uint32_t ifeature = 0; ifeature < nfeatures; ++ifeature) { - OGRFeature *geom_feature_out = OGRFeature::CreateFeature(geom_layer_out->GetLayerDefn()); - int32_t fid = FID_of_index[ifeature]; - OGRFeature *in_feature = layer->GetFeature(FID_of_index[ifeature]); - // TODO: if in_feature != NULL? - geom_feature_out->SetGeometry(in_feature->GetGeometryRef()); - geom_feature_out->SetFID(fid); - if (geom_layer_out->CreateFeature(geom_feature_out) != OGRERR_NONE) { - GCBS_ERROR("Failed to create output feature with FID '" + std::to_string(FID_of_index[ifeature]) + "' in '" + output_file + "'"); - throw std::string("Failed to create output feature with FID '" + std::to_string(FID_of_index[ifeature]) + "' in '" + output_file + "'"); - } - OGRFeature::DestroyFeature(geom_feature_out); - OGRFeature::DestroyFeature(in_feature); - } - GDALClose(gpkg_out); - GDALClose(in_ogr_dataset); - - uint16_t nthreads = config::instance()->get_default_chunk_processor()->max_threads(); - - std::mutex mutex; - std::vector workers; - std::vector out_temp_files; - for (uint16_t ithread = 0; ithread < nthreads; ++ithread) { - workers.push_back(std::thread([ithread, nthreads, &cube, &agg_func_names, &agg_func_creators, nfeatures, &features_in_chunk, &fid_column, &band_index, &output_file, &index_of_FID, FID_of_index, &mutex, &prg, &out_temp_files, &ogr_dataset, &ogr_layer](void) { - GDALDriver *gpkg_driver = GetGDALDriverManager()->GetDriverByName("GPKG"); - // if (gpkg_driver == NULL) { - // GCBS_ERROR("OGR GeoPackage driver not found"); - // throw std::string("OGR GeoPackage driver not found"); - // } - - GDALDataset *in_ogr_dataset = (GDALDataset *)GDALOpenEx(ogr_dataset.c_str(), GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, - NULL); - // if (in_ogr_dataset == NULL) { - // GCBS_ERROR("failed to open '" + ogr_dataset + "'"); - // throw std::string("failed to open '" + ogr_dataset + "'"); - // } - - OGRLayer *layer = in_ogr_dataset->GetLayerByName(ogr_layer.c_str()); - // if (layer == NULL) { - // GCBS_ERROR("invalid OGR layer"); - // throw std::string("invalid OGR layer"); - // } - - for (uint32_t ct = ithread; ct < cube->count_chunks_t(); ct += nthreads) { - std::string output_file_cur = filesystem::join(filesystem::get_tempdir(), utils::generate_unique_filename(8, "zs_" + std::to_string(ithread) + "_", ".gpkg")); - GDALDataset *gpkg_out = gpkg_driver->Create(output_file_cur.c_str(), 0, 0, 0, GDT_Unknown, NULL); - if (gpkg_out == NULL) { - GCBS_ERROR("Creation of GPKG file '" + output_file_cur + "' failed"); - throw std::string("Creation of GPKG file '" + output_file_cur + "' failed"); - } - - // initialize per geometry + time aggregators - uint32_t nt = cube->chunk_size(cube->chunk_id_from_coords({ct, 0, 0}))[0]; - std::vector> pixel_aggregators; - for (uint16_t i = 0; i < agg_func_names.size(); ++i) { - pixel_aggregators.push_back(agg_func_creators[i]()); - pixel_aggregators[i]->init(nfeatures, nt); - } - - for (uint32_t cx = 0; cx < cube->count_chunks_x(); ++cx) { - for (uint32_t cy = 0; cy < cube->count_chunks_y(); ++cy) { - chunkid_t id = cube->chunk_id_from_coords({ct, cy, cx}); - chunkid_t id_spatial = cube->chunk_id_from_coords({0, cy, cx}); - - bounds_st chunk_bounds = cube->bounds_from_chunk(id); - - // if chunk intersects with at least one geometry - if (features_in_chunk.count(id_spatial) > 0) { - // read chunk - std::shared_ptr chunk = cube->read_chunk(id); - - if (chunk->empty()) { - continue; - } - - // iterate over all features that intersect spatially with current chunk - for (uint32_t ifeature = 0; ifeature < features_in_chunk[id_spatial].size(); ++ifeature) { - int32_t fid = features_in_chunk[id_spatial][ifeature]; - - // filter chunk pixels by bounding box of current feature - OGRFeature *feature = layer->GetFeature(fid); - OGREnvelope feature_bbox; - feature->GetGeometryRef()->getEnvelope(&feature_bbox); - - int32_t x_start = std::max((int32_t)std::floor((feature_bbox.MinX - chunk_bounds.s.left) / - cube->st_reference()->dx()), - 0); - int32_t x_end = std::min((int32_t)std::ceil((feature_bbox.MaxX - chunk_bounds.s.left) / - cube->st_reference()->dx()), - (int32_t)chunk->size()[3] - 1); - - bool outside = false; - if (x_end < 0 || x_start > int32_t(chunk->size()[3] - 1)) { - outside = true; - } - - int32_t y_start = std::max((int32_t)std::floor((chunk_bounds.s.top - feature_bbox.MaxY) / - cube->st_reference()->dy()), - 0); - int32_t y_end = std::min((int32_t)std::ceil((chunk_bounds.s.top - feature_bbox.MinY) / - cube->st_reference()->dy()), - (int32_t)chunk->size()[2] - 1); - - if (y_end < 0 || y_start > int32_t(chunk->size()[2] - 1)) { - outside = true; - } - - if (!outside) { - // rasterize - - CPLStringList rasterize_args; - rasterize_args.AddString("-burn"); - rasterize_args.AddString("1"); - rasterize_args.AddString("-ot"); - rasterize_args.AddString("Byte"); - rasterize_args.AddString("-of"); - rasterize_args.AddString("MEM"); - rasterize_args.AddString("-init"); - rasterize_args.AddString("0"); - rasterize_args.AddString("-tr"); - rasterize_args.AddString(std::to_string(cube->st_reference()->dx()).c_str()); - rasterize_args.AddString(std::to_string(cube->st_reference()->dy()).c_str()); - rasterize_args.AddString("-te"); - rasterize_args.AddString(std::to_string(chunk_bounds.s.left + x_start * - cube->st_reference()->dx()) - .c_str()); // xmin - rasterize_args.AddString(std::to_string(chunk_bounds.s.top - (y_end + 1) * - cube->st_reference()->dy()) - .c_str()); // ymin - rasterize_args.AddString(std::to_string(chunk_bounds.s.left + (x_end + 1) * - cube->st_reference()->dx()) - .c_str()); // xmax - rasterize_args.AddString(std::to_string(chunk_bounds.s.top - y_start * - cube->st_reference()->dy()) - .c_str()); // ymax - rasterize_args.AddString("-where"); - std::string where = fid_column + "=" + std::to_string(fid); - rasterize_args.AddString(where.c_str()); - rasterize_args.AddString("-l"); - rasterize_args.AddString(layer->GetName()); - - // log gdal_rasterize call - // std::stringstream ss; - // ss << "Running gdal_rasterize "; - // for (uint16_t iws = 0; iws < rasterize_args.size(); ++iws) { - // ss << rasterize_args[iws] << " "; - // } - // ss << ogr_dataset; - // GCBS_TRACE(ss.str()); - - GDALRasterizeOptions *rasterize_opts = GDALRasterizeOptionsNew(rasterize_args.List(), NULL); - if (rasterize_opts == NULL) { - GDALRasterizeOptionsFree(rasterize_opts); - throw std::string("ERROR in vector_queries::zonal_statistics(): cannot create gdal_rasterize options."); - } - - int err = 0; - GDALDataset *gdal_rasterized = (GDALDataset *)GDALRasterize("", NULL, (GDALDatasetH)in_ogr_dataset, rasterize_opts, &err); - if (gdal_rasterized == NULL) { - GCBS_ERROR("gdal_rasterize failed for feature with FID " + std::to_string(fid)); - } - GDALRasterizeOptionsFree(rasterize_opts); - - uint8_t *geom_mask = (uint8_t *)std::malloc(sizeof(uint8_t) * (x_end - x_start + 1) * (y_end - y_start + 1)); - if (gdal_rasterized->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, x_end - x_start + 1, y_end - y_start + 1, geom_mask, x_end - x_start + 1, y_end - y_start + 1, GDT_Byte, 0, 0, NULL) != CE_None) { - GCBS_ERROR("RasterIO failed for feature with FID " + std::to_string(fid)); - } - - for (int32_t iy = y_start; iy <= y_end; ++iy) { - for (int32_t ix = x_start; ix <= x_end; ++ix) { - // if mask is 1 - if (geom_mask[(iy - y_start) * (x_end - x_start + 1) + ix - x_start] == 1) { - for (uint32_t it = 0; it < chunk->size()[1]; ++it) { - for (uint16_t ifield = 0; ifield < agg_func_names.size(); ++ifield) { - uint16_t b_index = band_index[ifield]; - double v = ((double *)chunk->buf())[b_index * chunk->size()[1] * chunk->size()[2] * chunk->size()[3] + - it * chunk->size()[2] * chunk->size()[3] + - iy * chunk->size()[3] + - ix]; - pixel_aggregators[ifield]->update(v, index_of_FID[fid], it); - } - } - } - } - } - std::free(geom_mask); - GDALClose(gdal_rasterized); - } - OGRFeature::DestroyFeature(feature); - } - } - } - } - - std::vector>> res; - for (uint16_t iband = 0; iband < agg_func_names.size(); ++iband) { - res.push_back(pixel_aggregators[iband]->finalize()); - } - - // write output layers - for (uint32_t it = 0; it < nt; ++it) { - std::string layer_name = "attr_" + cube->st_reference()->datetime_at_index((it + cube->chunk_limits({ct, 0, 0}).low[0])).to_string(); - - OGRLayer *cur_attr_layer_out = gpkg_out->CreateLayer(layer_name.c_str(), NULL, wkbNone, NULL); - if (cur_attr_layer_out == NULL) { - GCBS_ERROR("Failed to create output layer in '" + output_file + "'"); - throw std::string("Failed to create output layer in '" + output_file + "'"); - } - - for (uint16_t ifield = 0; ifield < agg_func_names.size(); ++ifield) { - std::string field_name = cube->bands().get(band_index[ifield]).name + "_" + agg_func_names[ifield]; - OGRFieldDefn oField(field_name.c_str(), OFTReal); - // TODO: set precision? set_nullable? - - if (cur_attr_layer_out->CreateField(&oField) != OGRERR_NONE) { - GCBS_ERROR("Failed to create output field '" + field_name + "' in '" + output_file + "'"); - throw std::string("Failed to create output field '" + field_name + "' in '" + output_file + "'"); - } - } - - if (cur_attr_layer_out->StartTransaction() != OGRERR_NONE) { - GCBS_WARN("Failed to start transaction on attribute table '" + layer_name + "'"); - } - for (uint32_t ifeature = 0; ifeature < nfeatures; ++ifeature) { - OGRFeature *cur_feature_out = OGRFeature::CreateFeature(cur_attr_layer_out->GetLayerDefn()); - for (uint16_t ifield = 0; ifield < agg_func_names.size(); ++ifield) { - double v = (*(res[ifield]))[ifeature * nt + it]; - cur_feature_out->SetField(ifield, v); - } - cur_feature_out->SetFID(FID_of_index.at(ifeature)); - if (cur_attr_layer_out->CreateFeature(cur_feature_out) != OGRERR_NONE) { - GCBS_ERROR("Failed to create output feature with FID '" + std::to_string(FID_of_index.at(ifeature)) + "' in '" + output_file + "'"); - throw std::string("Failed to create output feature with FID '" + std::to_string(FID_of_index.at(ifeature)) + "' in '" + output_file + "'"); - } - OGRFeature::DestroyFeature(cur_feature_out); - } - if (cur_attr_layer_out->CommitTransaction() != OGRERR_NONE) { - GCBS_WARN("Failed to commit transaction on attribute table '" + layer_name + "'"); - } - - prg->increment(double(1) / cube->size_t()); - } - mutex.lock(); - out_temp_files.push_back(output_file_cur); - mutex.unlock(); - GDALClose(gpkg_out); - } - GDALClose(in_ogr_dataset); - })); - } - for (uint16_t ithread = 0; ithread < nthreads; ++ithread) { - workers[ithread].join(); - } - - // Combine layers with ogr2ogr - CPLStringList ogr2ogr_args; - ogr2ogr_args.AddString("-update"); - ogr2ogr_args.AddString("-gt"); - ogr2ogr_args.AddString("65536"); - - GDALVectorTranslateOptions *ogr2ogr_opts = GDALVectorTranslateOptionsNew(ogr2ogr_args.List(), NULL); - if (ogr2ogr_opts == NULL) { - GDALVectorTranslateOptionsFree(ogr2ogr_opts); - throw std::string("ERROR in vector_queries::zonal_statistics(): cannot create ogr2ogr options."); - } - - for (uint32_t i = 0; i < out_temp_files.size(); ++i) { - GDALDataset *temp_in = (GDALDataset *)GDALOpenEx(out_temp_files[i].c_str(), GDAL_OF_VECTOR | GDAL_OF_UPDATE, NULL, NULL, NULL); - gpkg_out = (GDALDataset *)GDALVectorTranslate(output_file.c_str(), NULL, 1, (GDALDatasetH *)&temp_in, ogr2ogr_opts, NULL); - if (gpkg_out == NULL) { - GCBS_ERROR("ogr2ogr failed"); - } - GDALClose(gpkg_out); - GDALClose(temp_in); - - filesystem::remove(out_temp_files[i]); - } - GDALVectorTranslateOptionsFree(ogr2ogr_opts); - - gpkg_out = (GDALDataset *)GDALOpenEx(output_file.c_str(), GDAL_OF_VECTOR | GDAL_OF_UPDATE, NULL, NULL, NULL); - if (gpkg_out == NULL) { - GCBS_ERROR("Opening output GPKG file '" + output_file + "' failed"); - throw std::string("Opening output GPKG file '" + output_file + "' failed"); - } - - // create spatial views see https://gdal.org/drivers/vector/gpkg.html#spatial-views - const char *gpkg_has_column_md = gpkg_driver->GetMetadataItem("SQLITE_HAS_COLUMN_METADATA", NULL); - - int32_t srs_auth_code = 0; - const char *authcode_char = srs_features.GetAuthorityCode(NULL); - std::string authcode_str; - if (std::strlen(authcode_char) > 0) { - authcode_str = authcode_char; - } - if (!authcode_str.empty()) { - srs_auth_code = std::atoi(authcode_str.c_str()); - } - - if (gpkg_has_column_md == NULL || std::strcmp(gpkg_has_column_md, "YES") != 0) { - GCBS_WARN("GeoPackage OGR driver does not support SQLite column metadata; skipping creation of spatial views"); - } else { - if (srs_auth_code == 0) { - GCBS_WARN("Failed to identify authority code of spatial reference system; creating spatial views with missing SRS"); - } - - std::string field_names_str; - for (uint16_t ifield = 0; ifield < agg_func_names.size() - 1; ++ifield) { - field_names_str += (cube->bands().get(band_index[ifield]).name + "_" + agg_func_names[ifield]) + ","; - } - field_names_str += (cube->bands().get(band_index[agg_func_names.size() - 1]).name + "_" + agg_func_names[agg_func_names.size() - 1]); - - for (uint32_t it = 0; it < cube->size_t(); ++it) { - std::string layer_name = "attr_" + cube->st_reference()->datetime_at_index(it).to_string(); - std::string view_name = "map_" + cube->st_reference()->datetime_at_index(it).to_string(); - - std::string query; - query = "CREATE VIEW '" + view_name + "' AS SELECT geom.fid AS OGC_FID, geom.geom, " + field_names_str + " FROM geom JOIN '" + layer_name + - "' ON geom.fid = '" + layer_name + "'.fid;"; - gpkg_out->ExecuteSQL(query.c_str(), NULL, NULL); - query = "INSERT INTO gpkg_contents (table_name, identifier, data_type, srs_id) VALUES ( '" + view_name + - "', '" + view_name + "', 'features'," + std::to_string(srs_auth_code) + - ")"; - - gpkg_out->ExecuteSQL(query.c_str(), NULL, NULL); - query = "INSERT INTO gpkg_geometry_columns (table_name, column_name, geometry_type_name, srs_id, z, m) values ('" + - view_name + "', 'geom', 'GEOMETRY', " + std::to_string(srs_auth_code) + - ", 0, 0)"; - gpkg_out->ExecuteSQL(query.c_str(), NULL, NULL); - } - // fix geometry type column in gpkg_geometry_columns table - std::string query = "UPDATE gpkg_geometry_columns SET geometry_type_name = (SELECT geometry_type_name FROM gpkg_geometry_columns WHERE table_name = \"geom\")"; - gpkg_out->ExecuteSQL(query.c_str(), NULL, NULL); - //GDALClose(gpkg_out); - } - - GDALClose(gpkg_out); - if (in_ogr_was_transformed) { // if temporary copy (due to coordinate transformation needed) - filesystem::remove(ogr_dataset); - } - prg->finalize(); -} - -} // namespace gdalcubes diff --git a/src/gdalcubes/src/vector_queries.h b/src/gdalcubes/src/vector_queries.h deleted file mode 100644 index 61c2f48..0000000 --- a/src/gdalcubes/src/vector_queries.h +++ /dev/null @@ -1,86 +0,0 @@ -/* - MIT License - - Copyright (c) 2019 Marius Appel - - Permission is hereby granted, free of charge, to any person obtaining a copy - of this software and associated documentation files (the "Software"), to deal - in the Software without restriction, including without limitation the rights - to use, copy, modify, merge, publish, distribute, sublicense, and/or sell - copies of the Software, and to permit persons to whom the Software is - furnished to do so, subject to the following conditions: - - The above copyright notice and this permission notice shall be included in all - copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE - AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER - LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, - OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ -#ifndef VECTOR_H -#define VECTOR_H - -#include "cube.h" - -namespace gdalcubes { - -class vector_queries { - public: - /** - * Query values of a data cube at irregular spatiotmeporal points. - * - * @brief This function extracts values of data cube cells at provided spatiotemporal points. - * The result will contain data from all bands of the data cube, but will not contain the input point coordinates. - * - * @param cube data cube - * @param x x coordinates of query points - * @param y y coordinates of query points - * @param t date/time of query points - * @param srs spatial reference system of spatial point coordinates - * @return data frame (vector of double vectors) where first vector represents columns (bands) - */ - static std::vector> query_points(std::shared_ptr cube, std::vector x, std::vector y, std::vector t, std::string srs); - - /** - * Query complete time series of a datacube at irregular spatial points - * - * @brief This function extracts time series of data cube cells at provided spatial points. - * The result will contain data from all bands of the data cube, but will not contain the input point coordinates. - * - * @param cube data cube - * @param x x coordinates of query points - * @param y y coordinates of query points - * @param t date/time of query points - * @param srs spatial reference system of spatial point coordinates - * @return vector of data frames (vector of double vectors) where first vector represents bands, second vector represents - * data frame columns (time) and the third vector represents the points (data frame rows) - */ - static std::vector>> query_timeseries(std::shared_ptr cube, std::vector x, std::vector y, std::string srs); - - /** - * Query summary statistics of a data cube over spatial polygons - * - * The function produces a single geopackage file with one layer "geom" containing the geometries and feature IDs, and other layers containing resulting summary statistics - * per time slice of the data cube. These layers are named "attr_%DATETIME%" and contain only attribute values and feature IDs. Additional layers with names "map_%DATETIME%" - * are spatial views of these time slices by joining attribute layers with the geometry layer in a SQLite database view. - * - * Available aggregation functions currently include "min", "max", "mean", "median", "sum", "prod", "count", "var" and "sd". - * - * - * @param cube input data cube - * @param ogr_dataset input OGR dataset identifier with polygon geometries - * @param agg_band_functions vector of aggregation functions, band pairs representing combinations of available summary statistics functions (first element) and data cube bands (second element), e.g. {"min", "band1"}. - * @param out_path path where the resulting GeoPackage file will be written to - * @param overwrite_if_exists overwrite output file if already existing - * @param ogr_layer defines from which layer geometries are taken, if the ogr_dataset has multiple layers - */ - static void zonal_statistics(std::shared_ptr cube, std::string ogr_dataset, std::vector> agg_band_functions, std::string out_path, bool overwrite_if_exists = false, std::string ogr_layer = ""); -}; - -} // namespace gdalcubes - -#endif //VECTOR_H diff --git a/src/gdalcubes/src/warp.cpp b/src/gdalcubes/src/warp.cpp index 07191d0..dfc74eb 100644 --- a/src/gdalcubes/src/warp.cpp +++ b/src/gdalcubes/src/warp.cpp @@ -490,13 +490,13 @@ GDALDataset *gdalwarp_client::warp_complex(GDALDataset *in, std::string s_srs, s // see https://github.com/OSGeo/gdal/blob/64cf9b4e889c93e34177237665fe842186d1f581/alg/gdaltransformer.cpp#L1274C5-L1275C67 CPLStringList trnsfrm_opts; - // TODO: Do we need to check if t_srs and/or s_src are empty? trnsfrm_opts.AddNameValue("SRC_SRS", s_srs.c_str()); trnsfrm_opts.AddNameValue("GEOLOC_USE_TEMP_DATASETS", "NO"); GDALGenImgProjTransformInfo *psInfo = static_cast(GDALCreateGenImgProjTransformer2(in, out, trnsfrm_opts.List())); - //GDALGenImgProjTransformInfo *psInfo = static_cast(GDALCreateGenImgProjTransformer(in, s_srs.c_str(), out, NULL, TRUE, 0.0, 1 )); - - // TODO DO we need to check if psInfo is NULL? + if (!psInfo) { + GCBS_ERROR("Cannot find coordinate transformation from input image to target data cube"); + throw std::string("Cannot find coordinate transformation from input image to target data cube"); + } // Overwrite reprojection transform if needed, to benefit from cache diff --git a/src/gdalcubes/src/window_space.cpp b/src/gdalcubes/src/window_space.cpp index 385a897..3969751 100644 --- a/src/gdalcubes/src/window_space.cpp +++ b/src/gdalcubes/src/window_space.cpp @@ -219,7 +219,7 @@ std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { // 2. Apply padding (if needed) if (_pad.mode != padding::MODE::NONE) { - double v; + double v = NAN; if (_pad.mode == padding::MODE::CONSTANT) { v = _pad.constant_value; } @@ -586,7 +586,7 @@ std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { for (uint32_t ix = 0; ix < size_btyx[3]; ++ix) { // apply kernel - double sum = NAN; // TODO: complete.cases only? + double sum = 0.0; // TODO: complete.cases only? for (int16_t ky=0; ky < _win_size_y; ++ky) { for (int16_t kx=0; kx < _win_size_x; ++kx) { double v = ((double*)(cwin->buf()))[ @@ -595,12 +595,13 @@ std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { (iy+ky) * cwin->size()[3] + (ix + kx) ]; - if (std::isnan(sum)) { - sum = _kernel[ky * _win_size_x + kx] * v; - } - else if (std::isfinite(v)) + if (std::isfinite(v)) { // sum += _kernel[ky * _win_size_x + kx] * v; - + } + else { + sum = NAN; + break; + } } }