diff --git a/CMakeLists.txt b/CMakeLists.txt index 05bcccb2d..baeb22823 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -52,6 +52,7 @@ if(CMAKE_BUILD_TYPE MATCHES "Debug") "-misc-non-private-member-variables-in-classes," "-misc-include-cleaner," "-misc-use-anonymous-namespace," + "-misc-no-recursion," "hicpp-*," "-hicpp-special-member-functions," @@ -139,6 +140,8 @@ if(CMAKE_BUILD_TYPE MATCHES "Debug") "--suppress=constParameterReference:*/BathyRefractionCorrector.cpp" # in run() param extent cannot be const "--suppress=constParameterPointer:*/CcsdsPayloadDispatch.cpp" # Not trivial to fix, would have to change DispachObject class as well. "--suppress=knownConditionTrueFalse" # -1 < 0 etc + "--suppress=constVariableReference:*/RasterObject.cpp" # List [] const issue + "--suppress=constVariableReference:*/GeoIndexedRaster.cpp" # List [] const issue ) endif() diff --git a/datasets/gebco/package/GebcoBathyRaster.cpp b/datasets/gebco/package/GebcoBathyRaster.cpp index dbad2cb83..48da88f03 100644 --- a/datasets/gebco/package/GebcoBathyRaster.cpp +++ b/datasets/gebco/package/GebcoBathyRaster.cpp @@ -70,9 +70,10 @@ GebcoBathyRaster::~GebcoBathyRaster(void) = default; /*---------------------------------------------------------------------------- * getIndexFile *----------------------------------------------------------------------------*/ -void GebcoBathyRaster::getIndexFile(const OGRGeometry* geo, std::string& file) +void GebcoBathyRaster::getIndexFile(const OGRGeometry* geo, std::string& file, const std::vector* points) { static_cast(geo); + static_cast(points); file = filePath + "/" + indexFile; mlog(DEBUG, "Using %s", file.c_str()); } @@ -83,17 +84,17 @@ void GebcoBathyRaster::getIndexFile(const OGRGeometry* geo, std::string& file) *----------------------------------------------------------------------------*/ bool GebcoBathyRaster::findRasters(finder_t* finder) { - const OGRGeometry* geo = finder->geo; - const uint32_t start_indx = finder->range.start_indx; - const uint32_t end_indx = finder->range.end_indx; + const std::vector* flist = finder->featuresList; + const OGRGeometry* geo = finder->geo; + const uint32_t start = 0; + const uint32_t end = flist->size(); try { - for(uint32_t i = start_indx; i < end_indx; i++) + for(uint32_t i = start; i < end; i++) { - OGRFeature* feature = featuresList[i]; + OGRFeature* feature = flist->at(i); OGRGeometry *rastergeo = feature->GetGeometryRef(); - CHECKPTR(geo); if (!rastergeo->Intersects(geo)) continue; diff --git a/datasets/gebco/package/GebcoBathyRaster.h b/datasets/gebco/package/GebcoBathyRaster.h index 12b3cb9f2..b49ae0dcc 100644 --- a/datasets/gebco/package/GebcoBathyRaster.h +++ b/datasets/gebco/package/GebcoBathyRaster.h @@ -71,7 +71,7 @@ class GebcoBathyRaster: public GeoIndexedRaster GebcoBathyRaster (lua_State* L, GeoParms* _parms); ~GebcoBathyRaster (void) override; - void getIndexFile (const OGRGeometry* geo, std::string& file) final; + void getIndexFile (const OGRGeometry* geo, std::string& file, const std::vector* points) final; bool findRasters (finder_t* finder) final; /*-------------------------------------------------------------------- diff --git a/datasets/icesat2/package/MeritRaster.cpp b/datasets/icesat2/package/MeritRaster.cpp index f9e098c0b..90cf8e18f 100644 --- a/datasets/icesat2/package/MeritRaster.cpp +++ b/datasets/icesat2/package/MeritRaster.cpp @@ -117,6 +117,8 @@ uint32_t MeritRaster::getSamples (const MathLib::point_3d_t& point, int64_t gps, (void)param; (void)gps; + lockSampling(); + /* Determine Upper Left Coordinates */ int left_lon = ((int)floor(point.x / 5.0)) * 5; int upper_lat = ((int)ceil(point.y / 5.0)) * 5; @@ -203,6 +205,8 @@ uint32_t MeritRaster::getSamples (const MathLib::point_3d_t& point, int64_t gps, mlog(ERROR, "Failed to sample dataset: %s", dataset.c_str()); } + unlockSampling(); + return SS_NO_ERRORS; } diff --git a/datasets/landsat/package/LandsatHlsRaster.cpp b/datasets/landsat/package/LandsatHlsRaster.cpp index 7af241d12..59aa305fa 100644 --- a/datasets/landsat/package/LandsatHlsRaster.cpp +++ b/datasets/landsat/package/LandsatHlsRaster.cpp @@ -100,17 +100,15 @@ LandsatHlsRaster::LandsatHlsRaster(lua_State *L, GeoParms* _parms): VSIFCloseL(fp); /* Create dictionary of bands and algo names to process */ - bool returnBandSample; for(int i = 0; i < _parms->bands->length; i++) { const char* name = (*_parms->bands)[i].c_str(); if( isValidL8Band(name) || isValidS2Band(name) || isValidAlgoName(name)) { - if(!bandsDict.find(name, &returnBandSample)) - { - returnBandSample = true; - bandsDict.add(name, returnBandSample); - } + /* Add band to dictionary of bands but don't override if already there */ + auto it = bandsDict.find(name); + if(it == bandsDict.end()) + bandsDict[name] = true; if(strcasecmp((*_parms->bands)[i].c_str(), "NDSI") == 0) ndsi = true; if(strcasecmp((*_parms->bands)[i].c_str(), "NDVI") == 0) ndvi = true; @@ -123,24 +121,22 @@ LandsatHlsRaster::LandsatHlsRaster(lua_State *L, GeoParms* _parms): { for (unsigned i=0; iflags_file) { + /* Add band to dictionary of bands but don't override if already there */ const char* band = "Fmask"; - if(!bandsDict.find(band, &returnBandSample)) - { - returnBandSample = false; - bandsDict.add(band, returnBandSample); - } + auto it = bandsDict.find(band); + if(it == bandsDict.end()) + bandsDict[band] = false; } } @@ -155,9 +151,10 @@ LandsatHlsRaster::~LandsatHlsRaster(void) /*---------------------------------------------------------------------------- * getIndexFile *----------------------------------------------------------------------------*/ -void LandsatHlsRaster::getIndexFile(const OGRGeometry* geo, std::string& file) +void LandsatHlsRaster::getIndexFile(const OGRGeometry* geo, std::string& file, const std::vector* points) { static_cast(geo); + static_cast(points); file = indexFile; mlog(DEBUG, "Using %s", file.c_str()); } @@ -168,17 +165,17 @@ void LandsatHlsRaster::getIndexFile(const OGRGeometry* geo, std::string& file) *----------------------------------------------------------------------------*/ bool LandsatHlsRaster::findRasters(finder_t* finder) { - const OGRGeometry* geo = finder->geo; - const uint32_t start_indx = finder->range.start_indx; - const uint32_t end_indx = finder->range.end_indx; + const std::vector* flist = finder->featuresList; + const OGRGeometry* geo = finder->geo; + const uint32_t start = 0; + const uint32_t end = flist->size(); try { - for(uint32_t i = start_indx; i < end_indx; i++) + for(uint32_t i = start; i < end; i++) { - OGRFeature* feature = featuresList[i]; + OGRFeature* feature = flist->at(i); OGRGeometry *rastergeo = feature->GetGeometryRef(); - CHECKPTR(geo); if (!rastergeo->Intersects(geo)) continue; @@ -189,18 +186,13 @@ bool LandsatHlsRaster::findRasters(finder_t* finder) rgroup->gpsTime = getGmtDate(feature, DATE_TAG, rgroup->gmtDate); /* Find each requested band in the index file */ - bool val; - - bandsDictMutex.lock(); - const char* bandName = bandsDict.first(&val); - while(bandName != NULL) + for(auto it = bandsDict.begin(); it != bandsDict.end(); it++) { - /* skip algo names (NDIS, etc) */ - if(isValidAlgoName(bandName)) - { - bandName = bandsDict.next(&val); + const char* bandName = it->first.c_str(); + + /* skip algo names (NDIS, etc) and invalid bands (user error) */ + if(!isValidL8Band(bandName) && !isValidS2Band(bandName) && !isValidAlgoName(bandName)) continue; - } const char* fname = feature->GetFieldAsString(bandName); if(fname && strlen(fname) > 0) @@ -227,10 +219,7 @@ bool LandsatHlsRaster::findRasters(finder_t* finder) rgroup->infovect.push_back(rinfo); } } - bandName = bandsDict.next(&val); } - bandsDictMutex.unlock(); - mlog(DEBUG, "Added group: %s with %ld rasters", rgroup->id.c_str(), rgroup->infovect.size()); finder->rasterGroups.push_back(rgroup); } @@ -245,11 +234,54 @@ bool LandsatHlsRaster::findRasters(finder_t* finder) } -/*---------------------------------------------------------------------------- - * getGroupSamples - *----------------------------------------------------------------------------*/ -void LandsatHlsRaster::getGroupSamples (const rasters_group_t* rgroup, List& slist, uint32_t flags) +/****************************************************************************** + * PRIVATE METHODS + ******************************************************************************/ + +bool LandsatHlsRaster::validateBand(band_type_t type, const char* bandName) { + const char** tags; + int tagsCnt = 0; + + if(bandName == NULL) return false; + + if(type == SENTINEL2) + { + tags = S2_bands; + tagsCnt = S2_bandCnt; + } + else if(type == LANDSAT8) + { + tags = L8_bands; + tagsCnt = L8_bandCnt; + } + else if(type == ALGOBAND) + { + tags = ALGO_bands; + tagsCnt = ALGO_bandCnt; + } + else if(type == ALGONAME) + { + tags = ALGO_names; + tagsCnt = ALGO_nameCnt; + } + + for(int i = 0; i < tagsCnt; i++) + { + const char* tag = tags[i]; + if( strncasecmp(bandName, tag, strlen(tag)) == 0) + return true; + } + + return false; +} + + +uint32_t LandsatHlsRaster::_getGroupSamples(sample_mode_t mode, const rasters_group_t* rgroup, + List* slist, uint32_t flags, uint32_t pointIndx) +{ + uint32_t errors = SS_NO_ERRORS; + /* Which group is it? Landsat8 or Sentinel2 */ bool isL8 = false; bool isS2 = false; @@ -272,46 +304,108 @@ void LandsatHlsRaster::getGroupSamples (const rasters_group_t* rgroup, Listinfovect) + if(mode == SERIAL) { - const char* key = rinfo.fileName.c_str(); - cacheitem_t* item; - if(cache.find(key, &item)) + for(const auto& rinfo : rgroup->infovect) { - RasterSample* sample = item->sample; - if(sample) + const char* key = rinfo.fileName.c_str(); + cacheitem_t* item; + if(cache.find(key, &item)) { - sample->flags = flags; - - /* Is this band's sample to be returned to the user? */ - bool returnBandSample = false; - const char* bandName = rinfo.tag.c_str(); - if(bandsDict.find(bandName, &returnBandSample)) + RasterSample* sample = item->sample; + if(sample) { - if(returnBandSample) + sample->flags = flags; + + /* Is this band's sample to be returned to the user? */ + const char* bandName = rinfo.tag.c_str(); + auto it = bandsDict.find(bandName); + if(it != bandsDict.end()) { - slist.add(sample); - item->sample = NULL; + const bool returnBandSample = it->second; + if(returnBandSample) + { + slist->add(sample); + item->sample = NULL; + } } - } - /* green and red bands are the same for L8 and S2 */ - if(rinfo.tag == "B03") green = sample->value; - if(rinfo.tag == "B04") red = sample->value; + /* green and red bands are the same for L8 and S2 */ + if(rinfo.tag == "B03") green = sample->value; + if(rinfo.tag == "B04") red = sample->value; - if(isL8) - { - if(rinfo.tag == "B05") nir08 = sample->value; - if(rinfo.tag == "B06") swir16 = sample->value; + if(isL8) + { + if(rinfo.tag == "B05") nir08 = sample->value; + if(rinfo.tag == "B06") swir16 = sample->value; + } + else /* Must be Sentinel2 */ + { + if(rinfo.tag == "B8A") nir08 = sample->value; + if(rinfo.tag == "B11") swir16 = sample->value; + } } - else /* Must be Sentinel2 */ + } + } + } + else if(mode == BATCH) + { + for(const auto& rinfo : rgroup->infovect) + { + /* This is the unique raster we are looking for, it cannot be NULL */ + unique_raster_t* ur = rinfo.uraster; + assert(ur); + + /* Get the sample for this point from unique raster */ + for(const point_sample_t& ps : ur->pointSamples) + { + if(ps.pointInfo.index == pointIndx) { - if(rinfo.tag == "B8A") nir08 = sample->value; - if(rinfo.tag == "B11") swir16 = sample->value; + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ + if(ps.sample == NULL) break; + + /* Is this band's sample to be returned to the user? */ + const char* bandName = rinfo.tag.c_str(); + + auto it = bandsDict.find(bandName); + if(it != bandsDict.end()) + { + const bool returnBandSample = it->second; + if(returnBandSample) + { + /* Create a copy of the sample and add it to the list */ + RasterSample* sample = new RasterSample(*ps.sample); + + /* Set flags for this sample */ + sample->flags = flags; + slist->add(sample); + } + } + errors |= ps.ssErrors; + + /* green and red bands are the same for L8 and S2 */ + if(rinfo.tag == "B03") green = ps.sample->value; + if(rinfo.tag == "B04") red = ps.sample->value; + + if(isL8) + { + if(rinfo.tag == "B05") nir08 = ps.sample->value; + if(rinfo.tag == "B06") swir16 = ps.sample->value; + } + else /* Must be Sentinel2 */ + { + if(rinfo.tag == "B8A") nir08 = ps.sample->value; + if(rinfo.tag == "B11") swir16 = ps.sample->value; + } + break; } } } } + else + { + throw RunTimeException(DEBUG, RTE_ERROR, "Invalid sample mode"); + } const double groupTime = rgroup->gpsTime / 1000; const std::string groupName = rgroup->id + " {\"algo\": \""; @@ -323,7 +417,7 @@ void LandsatHlsRaster::getGroupSamples (const rasters_group_t* rgroup, Listvalue = (green - swir16) / (green + swir16); else sample->value = invalid; - slist.add(sample); + slist->add(sample); } if(ndvi) @@ -332,7 +426,7 @@ void LandsatHlsRaster::getGroupSamples (const rasters_group_t* rgroup, Listvalue = (nir08 - red) / (nir08 + red); else sample->value = invalid; - slist.add(sample); + slist->add(sample); } if(ndwi) @@ -341,64 +435,8 @@ void LandsatHlsRaster::getGroupSamples (const rasters_group_t* rgroup, Listvalue = (nir08 - swir16) / (nir08 + swir16); else sample->value = invalid; - slist.add(sample); + slist->add(sample); } -} - - -/*---------------------------------------------------------------------------- - * getMaxBatchThreads - *----------------------------------------------------------------------------*/ -uint32_t LandsatHlsRaster::getMaxBatchThreads(void) -{ - /* - * The average number of strips for a point is ~10 - * Limit the number of batch threads to 8. - */ - return 8; + return errors; } - - - -/****************************************************************************** - * PRIVATE METHODS - ******************************************************************************/ - -bool LandsatHlsRaster::validateBand(band_type_t type, const char* bandName) -{ - const char** tags; - int tagsCnt = 0; - - if(bandName == NULL) return false; - - if(type == SENTINEL2) - { - tags = S2_bands; - tagsCnt = S2_bandCnt; - } - else if(type == LANDSAT8) - { - tags = L8_bands; - tagsCnt = L8_bandCnt; - } - else if(type == ALGOBAND) - { - tags = ALGO_bands; - tagsCnt = ALGO_bandCnt; - } - else if(type == ALGONAME) - { - tags = ALGO_names; - tagsCnt = ALGO_nameCnt; - } - - for(int i = 0; i < tagsCnt; i++) - { - const char* tag = tags[i]; - if( strncasecmp(bandName, tag, strlen(tag)) == 0) - return true; - } - - return false; -} \ No newline at end of file diff --git a/datasets/landsat/package/LandsatHlsRaster.h b/datasets/landsat/package/LandsatHlsRaster.h index 876125aee..cb2d9317a 100644 --- a/datasets/landsat/package/LandsatHlsRaster.h +++ b/datasets/landsat/package/LandsatHlsRaster.h @@ -75,21 +75,38 @@ class LandsatHlsRaster: public GeoIndexedRaster protected: + /*-------------------------------------------------------------------- + * Types + *--------------------------------------------------------------------*/ + + typedef enum { + BATCH, + SERIAL, + } sample_mode_t; + /*-------------------------------------------------------------------- * Methods *--------------------------------------------------------------------*/ - LandsatHlsRaster (lua_State* L, GeoParms* _parms); - ~LandsatHlsRaster (void) override; + LandsatHlsRaster (lua_State* L, GeoParms* _parms); + ~LandsatHlsRaster (void) override; + + void getIndexFile (const OGRGeometry* geo, std::string& file, const std::vector* points) final; + bool findRasters (finder_t* finder) final; + + void getGroupSamples (const rasters_group_t* rgroup, List& slist, uint32_t flags) final + { _getGroupSamples(SERIAL, rgroup, &slist, flags);} - void getIndexFile (const OGRGeometry* geo, std::string& file) final; - bool findRasters (finder_t* finder) final; - void getGroupSamples (const rasters_group_t* rgroup, List& slist, uint32_t flags) final; - uint32_t getMaxBatchThreads (void) final; + uint32_t getBatchGroupSamples(const rasters_group_t* rgroup, List* slist, uint32_t flags, uint32_t pointIndx) final + { return _getGroupSamples(BATCH, rgroup, slist, flags, pointIndx);} private: + /*-------------------------------------------------------------------- + * Methods + *--------------------------------------------------------------------*/ + static bool validateBand (band_type_t type, const char* bandName); static bool isValidL8Band (const char* bandName) {return validateBand(LANDSAT8, bandName);} @@ -97,18 +114,21 @@ class LandsatHlsRaster: public GeoIndexedRaster static bool isValidAlgoBand (const char* bandName) {return validateBand(ALGOBAND, bandName);} static bool isValidAlgoName (const char* bandName) {return validateBand(ALGONAME, bandName);} + uint32_t _getGroupSamples(sample_mode_t mode, const rasters_group_t* rgroup, + List* slist, uint32_t flags, uint32_t pointIndx=0); + /*-------------------------------------------------------------------- * Data *--------------------------------------------------------------------*/ std::string filePath; std::string indexFile; - Mutex bandsDictMutex; - Dictionary bandsDict; + std::unordered_map bandsDict; bool ndsi; bool ndvi; bool ndwi; + }; #endif /* __landsat_hls_raster__ */ diff --git a/datasets/landsat/selftests/landsat_reader.lua b/datasets/landsat/selftests/landsat_reader.lua index 483ae8fdb..fadb019c6 100644 --- a/datasets/landsat/selftests/landsat_reader.lua +++ b/datasets/landsat/selftests/landsat_reader.lua @@ -187,7 +187,7 @@ local tstrOverride = "2000:2:4:23:3:0" local expectedGroup = "T01UCS.2021001T225941.v2.0" print(string.format("\n--------------------------------\nTest: %s Temporal Filter Override: closest_time=%s\n--------------------------------", demType, tstrOverride, tstr)) -dem = geo.raster(geo.parms({ asset = demType, algorithm = "NearestNeighbour", radius = 0, with_flags=true, closest_time=tstr, bands = {"NDSI", "NDVI", "NDWI"}, catalog = contents })) +dem = geo.raster(geo.parms({ asset = demType, algorithm = "NearestNeighbour", radius = 0, with_flags=true, closest_time=tstr, bands = {"NDSI", "NDVI", "NDWI"}, catalog = contents })) runner.check(dem ~= nil) tbl, err = dem:sample(lon, lat, height, tstrOverride) runner.check(err == 0) diff --git a/datasets/landsat/selftests/plugin_unittest.lua b/datasets/landsat/selftests/plugin_unittest.lua new file mode 100644 index 000000000..cac7dc3ba --- /dev/null +++ b/datasets/landsat/selftests/plugin_unittest.lua @@ -0,0 +1,57 @@ +local runner = require("test_executive") +local console = require("console") +local asset = require("asset") +local assets = asset.loaddir() +local json = require("json") +local _,td = runner.srcscript() + +-- console.monitor:config(core.LOG, core.DEBUG) +-- sys.setlvl(core.LOG, core.DEBUG) + +-- Check If Present -- +if not core.UNITTEST then return end + +-- Setup -- +local assets = asset.loaddir() + +local script_parms = {earthdata="https://data.lpdaac.earthdatacloud.nasa.gov/s3credentials", identity="lpdaac-cloud"} +local earthdata_auth_script = core.script("earth_data_auth", json.encode(script_parms)) +while not aws.csget("lpdaac-cloud") do + print("Waiting to authenticate to LPDAAC...") + sys.wait(1) +end + + +local geojsonfile = td.."../data/hls_trimmed.geojson" +local f = io.open(geojsonfile, "r") +local contents = f:read("*all") +f:close() + +-- Unit Test -- + +local lon = -179.0 +local lat = 51.0 + +local lon_incr = 0.01 +local lat_incr = 0.00 +local pointCount = 100 + +print(string.format("\n-------------------------------------------------\nLandsat Plugin test (NDVI)\n-------------------------------------------------")) +local demType = "landsat-hls" +local dem = geo.raster(geo.parms({ asset = demType, algorithm = "NearestNeighbour", radius = 0, bands = {"NDVI"}, catalog = contents })) +runner.check(dem ~= nil) + +ut = geo.ut_sample(dem) +runner.check(ut ~= nil) +status = ut:test(lon, lat, lon_incr, lat_incr, pointCount) +runner.check(status, "Failed sampling test") +ut = nil + + + +-- Clean Up -- + +-- Report Results -- + +runner.report() + diff --git a/datasets/pgc/package/PgcDemStripsRaster.cpp b/datasets/pgc/package/PgcDemStripsRaster.cpp index ecc8ce109..53956529a 100644 --- a/datasets/pgc/package/PgcDemStripsRaster.cpp +++ b/datasets/pgc/package/PgcDemStripsRaster.cpp @@ -34,6 +34,8 @@ ******************************************************************************/ #include "PgcDemStripsRaster.h" +#include + /****************************************************************************** * PROTECTED METHODS @@ -59,7 +61,14 @@ PgcDemStripsRaster::PgcDemStripsRaster(lua_State *L, GeoParms* _parms, const cha /*---------------------------------------------------------------------------- * Destructor *----------------------------------------------------------------------------*/ -PgcDemStripsRaster::~PgcDemStripsRaster(void) = default; +PgcDemStripsRaster::~PgcDemStripsRaster(void) +{ + /* Remove combined geojson file */ + if(!combinedGeoJSON.empty()) + { + VSIUnlink(combinedGeoJSON.c_str()); + } +} /*---------------------------------------------------------------------------- @@ -82,137 +91,90 @@ bool PgcDemStripsRaster::getFeatureDate(const OGRFeature* feature, TimeLib::gmt_ /*---------------------------------------------------------------------------- - * openGeoIndex + * getIndexFile *----------------------------------------------------------------------------*/ -bool PgcDemStripsRaster::openGeoIndex(const OGRGeometry* geo) +void PgcDemStripsRaster::getIndexFile(const OGRGeometry* geo, std::string& file, const std::vector* points) { - /* For point call parent class */ - if(GdalRaster::ispoint(geo)) - return GeoIndexedRaster::openGeoIndex(geo); + const OGRPolygon* poly = NULL; - /* - * Create a list of minx, miny 1° x 1° geocell points contained in AOI - * For each point get geojson file - * Open file, get list of features. - */ - const OGRPolygon* poly = geo->toPolygon(); - OGREnvelope env; - poly->getEnvelope(&env); - - const double minx = floor(env.MinX); - const double miny = floor(env.MinY); - const double maxx = ceil(env.MaxX); - const double maxy = ceil(env.MaxY); + if(geo == NULL && points == NULL) + { + mlog(ERROR, "Both geo and points are NULL"); + ssErrors |= SS_INDEX_FILE_ERROR; + return; + } - /* Create poly geometry for all index files */ - geoIndexPoly = GdalRaster::makeRectangle(minx, miny, maxx, maxy); + /* Determine if we have a point */ + if(geo && GdalRaster::ispoint(geo)) + { + /* Only one index file for a point from one of the geocells */ + const OGRPoint* poi = geo->toPoint(); + _getIndexFile(poi->getX(), poi->getY(), file); + return; + } - emptyFeaturesList(); + /* Vector holding all geojson files from all geocells */ + std::vector files; - for(long ix = minx; ix < maxx; ix++ ) + /* Determine if we have a polygon */ + if(geo && GdalRaster::ispoly(geo)) { - for(long iy = miny; iy < maxy; iy++) - { - std::string newFile; - _getIndexFile(ix, iy, newFile); - - GDALDataset* dset = NULL; - try - { - /* Open new vector data set*/ - dset = static_cast(GDALOpenEx(newFile.c_str(), GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, NULL)); - if(dset == NULL) - { - /* If index file for this ix, iy does not exist, continue */ - mlog(DEBUG, "Failed to open geojson index file: %s:", newFile.c_str()); - continue; - } + OGREnvelope env; + poly = geo->toPolygon(); + poly->getEnvelope(&env); - OGRLayer* layer = dset->GetLayer(0); - CHECKPTR(layer); + const double minx = floor(env.MinX); + const double miny = floor(env.MinY); + const double maxx = ceil(env.MaxX); + const double maxy = ceil(env.MaxY); - /* - * Clone all features and store them for performance/speed of feature lookup - */ - layer->ResetReading(); - while(OGRFeature* feature = layer->GetNextFeature()) + for(long ix = minx; ix < maxx; ix++) + { + for(long iy = miny; iy < maxy; iy++) + { + std::string newFile; + _getIndexFile(ix, iy, newFile); + if(!newFile.empty()) { - /* Temporal filter */ - TimeLib::gmt_time_t gmtDate; - if(parms->filter_time && getFeatureDate(feature, gmtDate)) - { - /* Check if feature is in time range */ - if(!TimeLib::gmtinrange(gmtDate, parms->start_time, parms->stop_time)) - { - OGRFeature::DestroyFeature(feature); - continue; - } - } - - /* Clone feature and store it */ - OGRFeature* fp = feature->Clone(); - featuresList.push_back(fp); - OGRFeature::DestroyFeature(feature); + files.push_back(newFile); } - - mlog(DEBUG, "Loaded %lu index file features/rasters from: %s", featuresList.size(), newFile.c_str()); - GDALClose((GDALDatasetH)dset); - } - catch(const RunTimeException& e) - { - /* If geocell does not have a geojson index file, ignore it and don't count it as error */ - if(dset) GDALClose((GDALDatasetH)dset); } } + mlog(INFO, "Found %ld geojson files in polygon", files.size()); } - if(featuresList.empty()) + /* If we don't have a polygon but we have points get all files for the points */ + if(poly == NULL && points != NULL) { - /* All geocells were 'empty' */ - geoIndexPoly.empty(); - ssError |= SS_INDEX_FILE_ERROR; - return false; + for(const auto& p : *points) + { + std::string newFile; + _getIndexFile(p.point.x, p.point.y, newFile); + if(!newFile.empty()) + { + files.push_back(newFile); + } + } + mlog(INFO, "Found %zu geojson files with %zu points", files.size(), points->size()); } - return true; -} + /* Remove any duplicate files */ + std::sort(files.begin(), files.end()); + files.erase(std::unique(files.begin(), files.end()), files.end()); -/*---------------------------------------------------------------------------- - * getIndexFile - *----------------------------------------------------------------------------*/ -void PgcDemStripsRaster::getIndexFile(const OGRGeometry* geo, std::string& file) -{ - if(GdalRaster::ispoint(geo)) + /* Combine all geojson files into a single file stored in vsimem */ + if(!combinedGeoJSON.empty()) { - const OGRPoint* poi = geo->toPoint(); - _getIndexFile(poi->getX(), poi->getY(), file); + /* Remove previous combined geojson file */ + VSIUnlink(combinedGeoJSON.c_str()); } -} - -/*---------------------------------------------------------------------------- - * getMaxBatchThreads - *----------------------------------------------------------------------------*/ -uint32_t PgcDemStripsRaster::getMaxBatchThreads(void) -{ - /* - * Typically, the average number of strips for a point ranges between 10 to 20, - * but in some areas, the number can exceed 100. To avoid overwhelming the system - * in these high-density areas we limit the number of batch threads to 1 - * . - */ - uint32_t numThreads = 1; - /* - * If we are filtering by closest time or using POI time, we only process - * the data strip raster and the quality mask raster. In this scenario, - * we allow the default number of threads to maximize performance. - */ - if(parms->filter_closest_time || parms->use_poi_time) + combinedGeoJSON = "/vsimem/" + GdalRaster::getUUID() + "_combined.geojson"; + if(combineGeoJSONFiles(files)) { - numThreads = MAX_BATCH_THREADS; + /* Set the combined geojson file as the index file */ + file = combinedGeoJSON; } - - return numThreads; } /*---------------------------------------------------------------------------- @@ -220,10 +182,10 @@ uint32_t PgcDemStripsRaster::getMaxBatchThreads(void) *----------------------------------------------------------------------------*/ bool PgcDemStripsRaster::findRasters(finder_t* finder) { - - const OGRGeometry* geo = finder->geo; - const uint32_t start_indx = finder->range.start_indx; - const uint32_t end_indx = finder->range.end_indx; + const std::vector* flist = finder->featuresList; + const OGRGeometry* geo = finder->geo; + const uint32_t start = 0; + const uint32_t end = flist->size(); /* * Find rasters and their dates. @@ -239,11 +201,10 @@ bool PgcDemStripsRaster::findRasters(finder_t* finder) */ try { - for(uint32_t i = start_indx; i < end_indx; i++) + for(uint32_t i = start; i < end; i++) { - OGRFeature* feature = featuresList[i]; + OGRFeature* feature = flist->at(i); OGRGeometry* rastergeo = feature->GetGeometryRef(); - CHECKPTR(geo); if (!rastergeo->Intersects(geo)) continue; @@ -362,3 +323,107 @@ void PgcDemStripsRaster::_getIndexFile(double lon, double lat, std::string& file mlog(DEBUG, "Using %s", file.c_str()); } + +/*---------------------------------------------------------------------------- + * combineGeoJSONFiles + *----------------------------------------------------------------------------*/ +bool PgcDemStripsRaster::combineGeoJSONFiles(const std::vector& inputFiles) +{ + /* Create an in-memory data source for the output */ + GDALDriver* memDriver = GetGDALDriverManager()->GetDriverByName("Memory"); + GDALDataset* memDataset = memDriver->Create("memory", 0, 0, 0, GDT_Unknown, NULL); + + if (memDataset == NULL) + { + mlog(ERROR, "Failed to create in-memory dataset."); + return false; + } + + OGRLayer* combinedLayer = NULL; + + for (const auto& infile : inputFiles) + { + GDALDataset* inputDataset = static_cast(GDALOpenEx(infile.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL)); + if (inputDataset == NULL) + { + mlog(DEBUG, "Failed to open input file: %s", infile.c_str()); + continue; + } + + /* Assuming that each file has a single layer */ + OGRLayer* inputLayer = inputDataset->GetLayer(0); + if (inputLayer == NULL) + { + mlog(ERROR, "No layer found in file: %s", infile.c_str()); + GDALClose(inputDataset); + continue; + } + + /* If this is the first file, create the combined layer in the memory dataset */ + if (combinedLayer == NULL) + { + combinedLayer = memDataset->CreateLayer(inputLayer->GetName(), inputLayer->GetSpatialRef(), inputLayer->GetGeomType(), NULL); + if (combinedLayer == NULL) + { + mlog(ERROR, "Failed to create combined layer in memory dataset."); + GDALClose(inputDataset); + GDALClose(memDataset); + return false; + } + + /* Copy the fields from the input layer to the combined layer */ + OGRFeatureDefn* inputFeatureDefn = inputLayer->GetLayerDefn(); + for (int i = 0; i < inputFeatureDefn->GetFieldCount(); i++) + combinedLayer->CreateField(inputFeatureDefn->GetFieldDefn(i)); + } + + /* Copy features from the input layer to the combined layer */ + inputLayer->ResetReading(); + OGRFeature* inputFeature = NULL; + while ((inputFeature = inputLayer->GetNextFeature()) != NULL) + { + OGRErr err = OGRERR_NONE; + OGRFeature* combinedFeature = OGRFeature::CreateFeature(combinedLayer->GetLayerDefn()); + err |= combinedFeature->SetFrom(inputFeature); + err |= combinedLayer->CreateFeature(combinedFeature); + if(err != OGRERR_NONE) + { + mlog(ERROR, "Failed to copy feature from input layer to combined layer."); + OGRFeature::DestroyFeature(inputFeature); + OGRFeature::DestroyFeature(combinedFeature); + GDALClose(inputDataset); + GDALClose(memDataset); + return false; + } + OGRFeature::DestroyFeature(inputFeature); + OGRFeature::DestroyFeature(combinedFeature); + } + + GDALClose(inputDataset); + } + + /* Write the combined layer to a GeoJSON file in the /vsimem filesystem */ + GDALDriver* jsonDriver = GetGDALDriverManager()->GetDriverByName("GeoJSON"); + if (jsonDriver == NULL) + { + mlog(ERROR, "GeoJSON driver not available."); + GDALClose(memDataset); + return false; + } + + GDALDataset* vsiDataset = jsonDriver->CreateCopy(combinedGeoJSON.c_str(), memDataset, FALSE, NULL, NULL, NULL); + if (vsiDataset == NULL) + { + mlog(ERROR, "Failed to create GeoJSON in /vsimem."); + } + else + { + mlog(DEBUG, "GeoJSON successfully created: %s", combinedGeoJSON.c_str()); + } + + /* Cleanup */ + GDALClose(vsiDataset); + GDALClose(memDataset); + + return true; +} diff --git a/datasets/pgc/package/PgcDemStripsRaster.h b/datasets/pgc/package/PgcDemStripsRaster.h index 7264acc63..744d5f141 100644 --- a/datasets/pgc/package/PgcDemStripsRaster.h +++ b/datasets/pgc/package/PgcDemStripsRaster.h @@ -53,13 +53,12 @@ class PgcDemStripsRaster: public GeoIndexedRaster PgcDemStripsRaster (lua_State* L, GeoParms* _parms, const char* dem_name, const char* geo_suffix, GdalRaster::overrideCRS_t cb); ~PgcDemStripsRaster (void) override; bool getFeatureDate (const OGRFeature* feature, TimeLib::gmt_time_t& gmtDate) final; - bool openGeoIndex (const OGRGeometry* geo) final; - void getIndexFile (const OGRGeometry* geo, std::string& file) final; + void getIndexFile (const OGRGeometry* geo, std::string& file, const std::vector* points) final; bool findRasters (finder_t* finder) final; - uint32_t getMaxBatchThreads (void) final; private: void _getIndexFile (double lon, double lat, std::string& file); + bool combineGeoJSONFiles (const std::vector& inputFiles); /*-------------------------------------------------------------------- * Data @@ -67,6 +66,8 @@ class PgcDemStripsRaster: public GeoIndexedRaster std::string filePath; std::string demName; std::string path2geocells; + std::string combinedGeoJSON; + }; #endif /* __pgcdem_strips_raster__ */ diff --git a/datasets/pgc/selftests/plugin_unittest.lua b/datasets/pgc/selftests/plugin_unittest.lua index 9b08c4371..54f80d353 100644 --- a/datasets/pgc/selftests/plugin_unittest.lua +++ b/datasets/pgc/selftests/plugin_unittest.lua @@ -13,7 +13,6 @@ if not core.UNITTEST then return end -- Unit Test -- - print('\n------------------\nTest01 RasterSubset::rema-mosaic \n------------------') local demType = "rema-mosaic" @@ -37,6 +36,38 @@ runner.check(status, "Failed subset test") ut = nil +lon = -178.0 +lat = 51.7 +lon_incr = 0.0001 +lat_incr = 0.0001 +pointCount = 10 + +print('\n------------------\nTest03 RasterSampler::arcticdem-mosaic\n------------------') +demType = "arcticdem-mosaic" +dem = geo.raster(geo.parms({asset=demType, algorithm="NearestNeighbour", radius=0, with_flags=false})) +runner.check(dem ~= nil) + +ut = geo.ut_sample(dem) +runner.check(ut ~= nil) +status = ut:test(lon, lat, lon_incr, lat_incr, pointCount) +runner.check(status, "Failed sampling test") +ut = nil + + +pointCount = 2 +lon_incr = 2.0 +print('\n------------------\nTest04 RasterSampler::arcticdem-strips\n------------------') +demType = "arcticdem-strips" +dem = geo.raster(geo.parms({asset=demType, algorithm="NearestNeighbour", radius=0, with_flags=true})) +runner.check(dem ~= nil) + +ut = geo.ut_sample(dem) +runner.check(ut ~= nil) +status = ut:test(lon, lat, lon_incr, lat_incr, pointCount) +runner.check(status, "Failed sampling test") +ut = nil + + -- Clean Up -- -- Report Results -- diff --git a/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp b/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp index 053b5a42a..4be50ac5e 100644 --- a/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp +++ b/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.cpp @@ -83,9 +83,10 @@ Usgs3dep1meterDemRaster::~Usgs3dep1meterDemRaster(void) /*---------------------------------------------------------------------------- * getIndexFile *----------------------------------------------------------------------------*/ -void Usgs3dep1meterDemRaster::getIndexFile(const OGRGeometry* geo, std::string& file) +void Usgs3dep1meterDemRaster::getIndexFile(const OGRGeometry* geo, std::string& file, const std::vector* points) { static_cast(geo); + static_cast(points); file = indexFile; mlog(DEBUG, "Using %s", file.c_str()); } @@ -96,18 +97,18 @@ void Usgs3dep1meterDemRaster::getIndexFile(const OGRGeometry* geo, std::string& *----------------------------------------------------------------------------*/ bool Usgs3dep1meterDemRaster::findRasters(finder_t* finder) { - const OGRGeometry* geo = finder->geo; - const uint32_t start_indx = finder->range.start_indx; - const uint32_t end_indx = finder->range.end_indx; + const std::vector* flist = finder->featuresList; + const OGRGeometry* geo = finder->geo; + const uint32_t start = 0; + const uint32_t end = flist->size(); try { /* Linearly search through feature list */ - for(uint32_t i = start_indx; i < end_indx; i++) + for(uint32_t i = start; i < end; i++) { - OGRFeature* feature = featuresList[i]; + OGRFeature* feature = flist->at(i); OGRGeometry *rastergeo = feature->GetGeometryRef(); - CHECKPTR(geo); if (!rastergeo->Intersects(geo)) continue; diff --git a/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.h b/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.h index 3322de3ac..8d923302e 100644 --- a/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.h +++ b/datasets/usgs3dep/package/Usgs3dep1meterDemRaster.h @@ -72,7 +72,7 @@ class Usgs3dep1meterDemRaster: public GeoIndexedRaster Usgs3dep1meterDemRaster (lua_State* L, GeoParms* _parms); ~Usgs3dep1meterDemRaster (void) override; - void getIndexFile (const OGRGeometry* geo, std::string& file) final; + void getIndexFile (const OGRGeometry* geo, std::string& file, const std::vector* points) final; bool findRasters (finder_t* finder) final; static OGRErr overrideTargetCRS(OGRSpatialReference& target); diff --git a/packages/arrow/ArrowSampler.cpp b/packages/arrow/ArrowSampler.cpp index da938459b..c6751ca1a 100644 --- a/packages/arrow/ArrowSampler.cpp +++ b/packages/arrow/ArrowSampler.cpp @@ -145,42 +145,10 @@ ArrowSampler::BatchSampler::BatchSampler(const char* _rkey, RasterObject* _robj, *----------------------------------------------------------------------------*/ ArrowSampler::BatchSampler::~BatchSampler(void) { - clearSamples(); delete [] rkey; robj->releaseLuaObject(); } -/*---------------------------------------------------------------------------- - * clearSamples - *----------------------------------------------------------------------------*/ -void ArrowSampler::BatchSampler::clearSamples(void) -{ - for(ArrowSampler::sample_list_t* slist : samples) - { - delete slist; - } - samples.clear(); -} - -/*---------------------------------------------------------------------------- - * Reader Constructor - *----------------------------------------------------------------------------*/ -ArrowSampler::Reader::Reader(RasterObject* _robj, ArrowSampler* _obj) : - robj(_robj), - obj(_obj), - range({0, 0}) -{ -} - -/*---------------------------------------------------------------------------- - * Reader Destructor - *----------------------------------------------------------------------------*/ -ArrowSampler::Reader::~Reader(void) -{ - delete robj; /* This is locally created RasterObject, not lua created */ -} - - /*---------------------------------------------------------------------------- * mainThread *----------------------------------------------------------------------------*/ @@ -197,7 +165,9 @@ void* ArrowSampler::mainThread(void* parm) { if(s->active) { - batchSampling(sampler); + const double t = TimeLib::latchtime(); + sampler->robj->getSamples(sampler->obj->points, sampler->samples); + mlog(INFO, "getSamples time: %.3lf", TimeLib::latchtime() - t); /* batchSampling can take minutes, check active again */ if(s->active) @@ -205,7 +175,7 @@ void* ArrowSampler::mainThread(void* parm) } /* Release since not needed anymore */ - sampler->clearSamples(); + sampler->samples.clear(); sampler->file_ids.clear(); } @@ -227,9 +197,8 @@ void* ArrowSampler::mainThread(void* parm) } catch(const RunTimeException& e) { - mlog(e.level(), "Error creating output file: %s", e.what()); + mlog(e.level(), "Error creating output file, PARQUET_ARROW reported: %s", e.what()); stop_trace(INFO, trace_id); - throw; } /* Signal Completion */ @@ -355,13 +324,11 @@ void ArrowSampler::Delete(void) active = false; delete mainPid; - parms->releaseLuaObject(); - for(batch_sampler_t* sampler : batchSamplers) + { + sampler->robj->stopSampling(); delete sampler; - - for(point_info_t* pinfo : points) - delete pinfo; + } delete [] dataFile; delete [] metadataFile; @@ -369,195 +336,6 @@ void ArrowSampler::Delete(void) delete [] outputMetadataPath; delete outQ; delete impl; -} - -/*---------------------------------------------------------------------------- - * getReadersRange - *----------------------------------------------------------------------------*/ -void ArrowSampler::getReadersRange(std::vector& ranges, uint32_t maxNumThreads) -{ - /* - * If points are geographically dispersed and fall into different data blocks of a raster, - * the initial read operation from the AWS S3 bucket can take approximately one second due - * to network latency and data retrieval time. Subsequent reads from the same data blocks - * are significantly faster due to caching mechanisms. - * - * The worst-case scenario occurs when points are not located within the same data block, - * leading to multiple time-consuming read operations. - * - * To optimize performance and balance the overhead of creating new RasterObjects and - * managing multiple threads, a threshold of 5 seconds (minPointsPerThread) is used. This value - * determines when to initiate multiple threads for parallel processing. By doing so, - * we aim to enhance efficiency and reduce overall processing time. - */ - const uint32_t minPointsPerThread = 5; - - /* Determine how many reader threads to use and index range for each */ - if(points.size() <= minPointsPerThread) - { - ranges.emplace_back(reader_range_t{0, static_cast(points.size())}); - return; - } - - uint32_t numThreads = std::min(maxNumThreads, static_cast(points.size()) / minPointsPerThread); - - /* Ensure at least two threads if points.size() > minPointsPerThread */ - if(numThreads == 1 && maxNumThreads > 1) - { - numThreads = 2; - } - - const uint32_t pointsPerThread = points.size() / numThreads; - uint32_t remainingPoints = points.size() % numThreads; - - uint32_t start = 0; - for(uint32_t i = 0; i < numThreads; i++) - { - const uint32_t end = start + pointsPerThread + (remainingPoints > 0 ? 1 : 0); - ranges.emplace_back(reader_range_t{start, end}); - - start = end; - if(remainingPoints > 0) - { - remainingPoints--; - } - } -} - -/*---------------------------------------------------------------------------- - * batchSampling - *----------------------------------------------------------------------------*/ -void ArrowSampler::batchSampling(batch_sampler_t* sampler) -{ - /* Get maximum number of batch processing threads allowed */ - const uint32_t maxNumThreads = sampler->robj->getMaxBatchThreads(); - - /* Get readers ranges */ - std::vector ranges; - sampler->obj->getReadersRange(ranges, maxNumThreads); - - for(uint32_t i = 0; i < ranges.size(); i++) - { - const reader_range_t& r = ranges[i]; - mlog(DEBUG, "%s: ragne-%u: %u to %u\n", sampler->rkey, i, r.start_indx, r.end_indx); - } - - const uint32_t numThreads = ranges.size(); - - if(numThreads == 1) - { - /* Single thread, read all samples in this thread using user RasterObject */ - readSamples(sampler->robj, ranges[0].start_indx, ranges[0].end_indx, sampler->obj, sampler->samples); - } - else - { - /* Start reader threads */ - std::vector pids; - std::vector readers; - for(uint32_t i = 0; i < numThreads; i++) - { - /* Create RasterObject for each reader. - * These are local objects and will be deleted in the reader destructor. - * User RasterObject is not used for sampling. It is used for acumulating samples from all readers. - */ - RasterObject* _robj = RasterObject::cppCreate(sampler->robj); - reader_t* reader = new reader_t(_robj, sampler->obj); - reader->range = ranges[i]; - readers.push_back(reader); - Thread* pid = new Thread(readerThread, reader); - pids.push_back(pid); - } - /* Wait for all reader threads to finish */ - for(Thread* pid : pids) - { - delete pid; - } - - /* Copy samples lists (slist pointers only) from each reader. */ - for(const reader_t* reader : readers) - { - for(sample_list_t* slist : reader->samples) - { - for(int32_t i = 0; i < slist->length(); i++) - { - /* NOTE: sample.fileId is an index of the file name in the reader's file dictionary. - * we need to convert it to the index in the batch sampler's dictionary (user's RasterObject dict). - */ - RasterSample* sample = slist->get(i); - - /* Find the file name for the sample id in reader's dictionary */ - const char* name = reader->robj->fileDictGetFile(sample->fileId); - - /* Use user's RasterObject dictionary to store the file names. */ - const uint64_t id = sampler->robj->fileDictAdd(name); - - /* Update the sample file id */ - sample->fileId = id; - } - - /* slist pointer is now in two samples vectors, one in the reader and one in the batch sampler - * batch sampler destructor will delete it - */ - sampler->samples.push_back(slist); - } - - delete reader; - } - } -} - -/*---------------------------------------------------------------------------- - * readerThread - *----------------------------------------------------------------------------*/ -void* ArrowSampler::readerThread(void* parm) -{ - reader_t* reader = static_cast(parm); - readSamples(reader->robj, - reader->range.start_indx, - reader->range.end_indx, - reader->obj, - reader->samples); - - /* Exit Thread */ - return NULL; -} - -/*---------------------------------------------------------------------------- - * readSamples - *----------------------------------------------------------------------------*/ -void* ArrowSampler::readSamples(RasterObject* robj, uint32_t start_indx, uint32_t end_indx, - ArrowSampler* obj, std::vector& samples) -{ - for(uint32_t i = start_indx; i < end_indx; i++) - { - if(!obj->active) break; // early exit if lua object is being destroyed - - point_info_t* pinfo = obj->points[i]; - - const MathLib::point_3d_t point = {pinfo->x, pinfo->y, 0.0}; - const double gps = robj->usePOItime() ? pinfo->gps : 0.0; - - sample_list_t* slist = new sample_list_t; - bool listvalid = true; - const uint32_t err = robj->getSamples(point, gps, *slist, NULL); - - if(err & SS_THREADS_LIMIT_ERROR) - { - listvalid = false; - mlog(CRITICAL, "Too many rasters to sample"); - } - - if(!listvalid) - { - /* Clear the list but don't delete it, empty slist indicates no samples for this point */ - slist->clear(); - } - - /* Add sample list */ - samples.push_back(slist); - } - - /* Exit Thread */ - return NULL; + parms->releaseLuaObject(); } \ No newline at end of file diff --git a/packages/arrow/ArrowSampler.h b/packages/arrow/ArrowSampler.h index 7ac11ced0..6d69305a8 100644 --- a/packages/arrow/ArrowSampler.h +++ b/packages/arrow/ArrowSampler.h @@ -75,51 +75,28 @@ class ArrowSampler: public LuaObject /*-------------------------------------------------------------------- * Types *--------------------------------------------------------------------*/ + typedef RasterObject::sample_list_t sample_list_t; + typedef RasterObject::point_info_t point_info_t; + typedef struct { const char* rkey; RasterObject* robj; } raster_info_t; - typedef struct - { - double x; - double y; - double gps; - } point_info_t; - - typedef List sample_list_t; typedef struct BatchSampler { const char* rkey; RasterObject* robj; ArrowSampler* obj; - std::vector samples; + List samples; std::set file_ids; std::vector> filemap; explicit BatchSampler (const char* _rkey, RasterObject* _robj, ArrowSampler* _obj); ~BatchSampler (void); - void clearSamples (void); } batch_sampler_t; - typedef struct - { - uint32_t start_indx; - uint32_t end_indx; - } reader_range_t; - - typedef struct Reader - { - RasterObject* robj; - ArrowSampler* obj; - reader_range_t range; - std::vector samples; - - explicit Reader (RasterObject* _robj, ArrowSampler* _obj); - ~Reader (void); - } reader_t; - /*-------------------------------------------------------------------- * Methods *--------------------------------------------------------------------*/ @@ -146,7 +123,7 @@ class ArrowSampler: public LuaObject Thread* mainPid; ArrowParms* parms; Publisher* outQ; - std::vector points; + std::vector points; std::vector batchSamplers; ArrowSamplerImpl* impl; const char* dataFile; // used locally to build parquet file @@ -163,12 +140,7 @@ class ArrowSampler: public LuaObject const char* outq_name, const std::vector& rasters); ~ArrowSampler (void) override; void Delete (void); - void getReadersRange (std::vector& ranges, uint32_t maxNumThreads); static void* mainThread (void* parm); - static void batchSampling (batch_sampler_t* sampler); - static void* readerThread (void* parm); - static void* readSamples (RasterObject* robj, uint32_t start_indx, uint32_t end_indx, - ArrowSampler* obj, std::vector& samples); }; #endif /* __parquet_sampler__*/ diff --git a/packages/arrow/ArrowSamplerImpl.cpp b/packages/arrow/ArrowSamplerImpl.cpp index 218c03a57..6d44aa44f 100644 --- a/packages/arrow/ArrowSamplerImpl.cpp +++ b/packages/arrow/ArrowSamplerImpl.cpp @@ -87,7 +87,7 @@ ArrowSamplerImpl::~ArrowSamplerImpl (void) /*---------------------------------------------------------------------------- * processInputFile *----------------------------------------------------------------------------*/ -void ArrowSamplerImpl::processInputFile(const char* file_path, std::vector& points) +void ArrowSamplerImpl::processInputFile(const char* file_path, std::vector& points) { try { @@ -276,7 +276,7 @@ void ArrowSamplerImpl::getMetadata(void) /*---------------------------------------------------------------------------- * getPoints *----------------------------------------------------------------------------*/ -void ArrowSamplerImpl::getPoints(std::vector& points) +void ArrowSamplerImpl::getPoints(std::vector& points) { if(asGeo) getGeoPoints(points); @@ -294,7 +294,7 @@ void ArrowSamplerImpl::getPoints(std::vector& point /* Update gps time for each point */ for(int64_t i = 0; i < time_column->length(); i++) - points[i]->gps = time_column->Value(i); + points[i].gps = time_column->Value(i); } else mlog(DEBUG, "Time column not found."); } @@ -302,7 +302,7 @@ void ArrowSamplerImpl::getPoints(std::vector& point /*---------------------------------------------------------------------------- * getXYPoints *----------------------------------------------------------------------------*/ -void ArrowSamplerImpl::getXYPoints(std::vector& points) +void ArrowSamplerImpl::getXYPoints(std::vector& points) { const std::vector columnNames = {xKey, yKey}; @@ -316,22 +316,21 @@ void ArrowSamplerImpl::getXYPoints(std::vector& poi auto x_column = std::static_pointer_cast(table->column(x_column_index)->chunk(0)); auto y_column = std::static_pointer_cast(table->column(y_column_index)->chunk(0)); - /* x and y columns are the same longth */ + /* x and y columns are the same length */ for(int64_t i = 0; i < x_column->length(); i++) { const double x = x_column->Value(i); const double y = y_column->Value(i); - ArrowSampler::point_info_t* pinfo = new ArrowSampler::point_info_t({x, y, 0.0}); - points.push_back(pinfo); + points.emplace_back(point_info_t({{x, y, 0.0}, 0.0})); } - mlog(DEBUG, "Read %ld points from file", points.size()); + mlog(DEBUG, "Read %zu points from file", points.size()); } /*---------------------------------------------------------------------------- * getGeoPoints *----------------------------------------------------------------------------*/ -void ArrowSamplerImpl::getGeoPoints(std::vector& points) +void ArrowSamplerImpl::getGeoPoints(std::vector& points) { const char* geocol = "geometry"; const std::vector columnNames = {geocol}; @@ -351,10 +350,10 @@ void ArrowSamplerImpl::getGeoPoints(std::vector& po { const std::string wkb_data = binary_array->GetString(i); /* Get WKB data as string (binary data) */ const ArrowCommon::wkbpoint_t point = convertWKBToPoint(wkb_data); - ArrowSampler::point_info_t* pinfo = new ArrowSampler::point_info_t({point.x, point.y, 0.0}); - points.push_back(pinfo); + + points.emplace_back(point_info_t({{point.x, point.y, 0.0}, 0.0})); } - mlog(INFO, "Read %ld geo points from file", points.size()); + mlog(INFO, "Read %zu geo points from file", points.size()); } /*---------------------------------------------------------------------------- @@ -441,7 +440,7 @@ std::shared_ptr ArrowSamplerImpl::addNewColumns(const std::shared_ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampler) { auto* pool = arrow::default_memory_pool(); - RasterObject* robj = sampler->robj; + const RasterObject* robj = sampler->robj; /* Create list builders for the new columns */ arrow::ListBuilder value_list_builder(pool, std::make_shared()); @@ -482,8 +481,10 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl std::shared_ptr count_list_array, min_list_array, max_list_array, mean_list_array, median_list_array, stdev_list_array, mad_list_array; /* Iterate over each sample in a vector of lists of samples */ - for(ArrowSampler::sample_list_t* slist : sampler->samples) + for(int i = 0; i < sampler->samples.length(); i++) { + sample_list_t* slist = sampler->samples[i]; + /* Start new lists */ PARQUET_THROW_NOT_OK(value_list_builder.Append()); PARQUET_THROW_NOT_OK(time_list_builder.Append()); @@ -507,9 +508,9 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl * If slist is empty the column will contain an empty list * to keep the number of rows consistent with the other columns */ - for(int i = 0; i < slist->length(); i++) + for(int j = 0; j < slist->length(); j++) { - const RasterSample* sample = slist->get(i); + const RasterSample* sample = slist->get(j); /* Append the value to the value list */ PARQUET_THROW_NOT_OK(value_builder->Append(sample->value)); @@ -624,7 +625,7 @@ bool ArrowSamplerImpl::makeColumnsWithLists(ArrowSampler::batch_sampler_t* sampl bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* sampler) { auto* pool = arrow::default_memory_pool(); - RasterObject* robj = sampler->robj; + const RasterObject* robj = sampler->robj; /* Create builders for the new columns */ arrow::DoubleBuilder value_builder(pool); @@ -647,8 +648,9 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s RasterSample fakeSample(0.0, 0); fakeSample.value = std::nan(""); - for(ArrowSampler::sample_list_t* slist : sampler->samples) + for(int i = 0; i < sampler->samples.length(); i++) { + sample_list_t* slist = sampler->samples[i]; RasterSample* sample; if(!slist->empty()) @@ -771,7 +773,7 @@ bool ArrowSamplerImpl::makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* s /*---------------------------------------------------------------------------- * getFirstValidSample *----------------------------------------------------------------------------*/ -RasterSample* ArrowSamplerImpl::getFirstValidSample(ArrowSampler::sample_list_t* slist) +RasterSample* ArrowSamplerImpl::getFirstValidSample(sample_list_t* slist) { for(int i = 0; i < slist->length(); i++) { diff --git a/packages/arrow/ArrowSamplerImpl.h b/packages/arrow/ArrowSamplerImpl.h index 3f13bd61b..a8da9b1fa 100644 --- a/packages/arrow/ArrowSamplerImpl.h +++ b/packages/arrow/ArrowSamplerImpl.h @@ -56,6 +56,9 @@ class ArrowSamplerImpl * Types *--------------------------------------------------------------------*/ + typedef RasterObject::sample_list_t sample_list_t; + typedef RasterObject::point_info_t point_info_t; + /*-------------------------------------------------------------------- * Methods *--------------------------------------------------------------------*/ @@ -63,7 +66,7 @@ class ArrowSamplerImpl explicit ArrowSamplerImpl (ArrowSampler* _sampler); ~ArrowSamplerImpl (void); - void processInputFile (const char* file_path, std::vector& points); + void processInputFile (const char* file_path, std::vector& points); bool processSamples (ArrowSampler::batch_sampler_t* sampler); void createOutpuFiles (void); @@ -99,14 +102,14 @@ class ArrowSamplerImpl *--------------------------------------------------------------------*/ void getMetadata (void); - void getPoints (std::vector& points); - void getXYPoints (std::vector& points); - void getGeoPoints (std::vector& points); + void getPoints (std::vector& points); + void getXYPoints (std::vector& points); + void getGeoPoints (std::vector& points); std::shared_ptr inputFileToTable (const std::vector& columnNames = {}); std::shared_ptr addNewColumns (const std::shared_ptr& table); bool makeColumnsWithLists (ArrowSampler::batch_sampler_t* sampler); bool makeColumnsWithOneSample(ArrowSampler::batch_sampler_t* sampler); - static RasterSample* getFirstValidSample (ArrowSampler::sample_list_t* slist); + static RasterSample* getFirstValidSample (sample_list_t* slist); static void tableToParquet (const std::shared_ptr& table, const char* file_path); static void tableToCsv (const std::shared_ptr& table, diff --git a/packages/geo/CMakeLists.txt b/packages/geo/CMakeLists.txt index 11157a67c..2928b3c44 100644 --- a/packages/geo/CMakeLists.txt +++ b/packages/geo/CMakeLists.txt @@ -36,6 +36,7 @@ if (GDAL_FOUND AND PROJ_FOUND AND TIFF_FOUND) ${CMAKE_CURRENT_LIST_DIR}/GeoParms.cpp ${CMAKE_CURRENT_LIST_DIR}/GeoLib.cpp $<$:${CMAKE_CURRENT_LIST_DIR}/UT_RasterSubset.cpp> + $<$:${CMAKE_CURRENT_LIST_DIR}/UT_RasterSample.cpp> ) target_include_directories (slideruleLib @@ -59,6 +60,7 @@ if (GDAL_FOUND AND PROJ_FOUND AND TIFF_FOUND) ${CMAKE_CURRENT_LIST_DIR}/GeoParms.h ${CMAKE_CURRENT_LIST_DIR}/GeoLib.h $<$:${CMAKE_CURRENT_LIST_DIR}/UT_RasterSubset.h> + $<$:${CMAKE_CURRENT_LIST_DIR}/UT_RasterSample.h> DESTINATION ${INCDIR} ) diff --git a/packages/geo/GeoIndexedRaster.cpp b/packages/geo/GeoIndexedRaster.cpp index f6e2dc923..2cb3f4434 100644 --- a/packages/geo/GeoIndexedRaster.cpp +++ b/packages/geo/GeoIndexedRaster.cpp @@ -59,14 +59,14 @@ const char* GeoIndexedRaster::DATE_TAG = "datetime"; /*---------------------------------------------------------------------------- * Reader Constructor *----------------------------------------------------------------------------*/ -GeoIndexedRaster::Reader::Reader (GeoIndexedRaster* raster): - obj(raster), +GeoIndexedRaster::Reader::Reader (GeoIndexedRaster* _obj): + obj(_obj), geo(NULL), entry(NULL), sync(NUM_SYNC_SIGNALS), run(true) { - thread = new Thread(GeoIndexedRaster::readerThread, this); + thread = new Thread(readerThread, this); } /*---------------------------------------------------------------------------- @@ -90,14 +90,10 @@ GeoIndexedRaster::Reader::~Reader (void) /*---------------------------------------------------------------------------- * Finder Constructor *----------------------------------------------------------------------------*/ -GeoIndexedRaster::Finder::Finder (GeoIndexedRaster* raster): - obj(raster), - geo(NULL), - range{0, 0}, - sync(NUM_SYNC_SIGNALS), - run(true) +GeoIndexedRaster::Finder::Finder (OGRGeometry* _geo, std::vector* _featuresList): + geo(_geo), + featuresList(_featuresList) { - thread = new Thread(GeoIndexedRaster::finderThread, this); } /*---------------------------------------------------------------------------- @@ -105,19 +101,32 @@ GeoIndexedRaster::Finder::Finder (GeoIndexedRaster* raster): *----------------------------------------------------------------------------*/ GeoIndexedRaster::Finder::~Finder (void) { - sync.lock(); - { - run = false; /* Set run flag to false */ - sync.signal(DATA_TO_SAMPLE, Cond::NOTIFY_ONE); - } - sync.unlock(); - - delete thread; /* delete thread waits on thread to join */ - /* geometry geo is cloned not 'newed' on GDAL heap. Use this call to free it */ OGRGeometryFactory::destroyGeometry(geo); } +/*---------------------------------------------------------------------------- + * UnionMaker Constructor + *----------------------------------------------------------------------------*/ +GeoIndexedRaster::UnionMaker::UnionMaker(GeoIndexedRaster* _obj, const std::vector* _points): + obj(_obj), + pointsRange({0, 0}), + points(_points), + unionPolygon(NULL), + stats({0.0, 0.0}) +{ +} + +/*---------------------------------------------------------------------------- + * GroupsFinder Constructor + *----------------------------------------------------------------------------*/ +GeoIndexedRaster::GroupsFinder::GroupsFinder(GeoIndexedRaster* _obj, const std::vector* _points): + obj(_obj), + pointsRange({0, 0}), + points(_points) +{ +} + /*---------------------------------------------------------------------------- * init *----------------------------------------------------------------------------*/ @@ -140,15 +149,16 @@ uint32_t GeoIndexedRaster::getSamples(const MathLib::point_3d_t& point, int64_t { static_cast(param); - samplingMutex.lock(); + lockSampling(); try { - ssError = SS_NO_ERRORS; + GroupOrdering groupList; + OGRPoint ogr_point(point.x, point.y, point.z); - OGRPoint ogr_point(point.x, point.y, point.z); + ssErrors = SS_NO_ERRORS; /* Sample Rasters */ - if(sample(&ogr_point, gps)) + if(sample(&ogr_point, gps, &groupList)) { /* Populate Return Vector of Samples (slist) */ const GroupOrdering::Iterator iter(groupList); @@ -187,9 +197,125 @@ uint32_t GeoIndexedRaster::getSamples(const MathLib::point_3d_t& point, int64_t } key = cache.next(&item); } - samplingMutex.unlock(); + unlockSampling(); - return ssError; + return ssErrors; +} + +/*---------------------------------------------------------------------------- + * getSamples + *----------------------------------------------------------------------------*/ +uint32_t GeoIndexedRaster::getSamples(const std::vector& points, List& sllist, void* param) +{ + static_cast(param); + + lockSampling(); + + perfStats.clear(); + + /* Vector of points and their associated raster groups */ + std::vector pointsGroups; + + /* Vector of rasters and all points they contain */ + std::vector uniqueRasters; + + try + { + ssErrors = SS_NO_ERRORS; + + /* Open the index file for all points */ + openGeoIndex(NULL, &points); + + /* Rasters to points map */ + raster_points_map_t rasterToPointsMap; + + /* For all points from the caller, create a vector of raster group lists */ + if(!findAllGroups(&points, pointsGroups, rasterToPointsMap)) + { + throw RunTimeException(CRITICAL, RTE_ERROR, "Error creating groups"); + } + + /* For all points from the caller, create a vector of unique rasters */ + if(!findUniqueRasters(uniqueRasters, pointsGroups, rasterToPointsMap)) + { + throw RunTimeException(CRITICAL, RTE_ERROR, "Error finding unique rasters"); + } + + /* For all unique rasters, sample them */ + if(!sampleUniqueRasters(uniqueRasters)) + { + throw RunTimeException(CRITICAL, RTE_ERROR, "Error sampling unique rasters"); + } + + /* + * Populate sllist with samples + */ + if(isSampling()) + { + mlog(DEBUG, "Populating sllist with samples"); + for(uint32_t pointIndx = 0; pointIndx < pointsGroups.size(); pointIndx++) + { + const point_groups_t& pg = pointsGroups[pointIndx]; + const GroupOrdering::Iterator iter(*pg.groupList); + + /* Allocate a new sample list for groupList */ + sample_list_t* slist = new sample_list_t(); + + for(int i = 0; i < iter.length; i++) + { + const rasters_group_t* rgroup = iter[i].value; + uint32_t flags = 0; + + /* Get flags value for this group of rasters */ + if(parms->flags_file) + flags = getBatchGroupFlags(rgroup, pointIndx); + + ssErrors |= getBatchGroupSamples(rgroup, slist, flags, pointIndx); + } + + sllist.add(slist); + } + } + else + { + /* Sampling has been stopped, return empty sllist */ + sllist.clear(); + } + } + catch (const RunTimeException &e) + { + mlog(e.level(), "Error getting samples: %s", e.what()); + } + + /* Clean up pointsGroups */ + for(const point_groups_t& pg : pointsGroups) + delete pg.groupList; + + /* Clean up unique rasters */ + for(unique_raster_t* ur : uniqueRasters) + { + for(const point_sample_t& ps : ur->pointSamples) + { + /* Delete samples which have not been returned (quality masks, etc) */ + if(!ps.sampleReturned) + { + delete ps.sample; + } + } + + delete ur; + } + + unlockSampling(); + + /* Print performance stats */ + mlog(INFO, "Performance Stats:"); + mlog(INFO, "spatialFilter: %12.3lf", perfStats.spatialFilterTime); + mlog(INFO, "findingRasters:%12.3lf", perfStats.findRastersTime); + mlog(INFO, "findingUnique: %12.3lf", perfStats.findUniqueRastersTime); + mlog(INFO, "sampling: %12.3lf", perfStats.samplesTime); + + return ssErrors; } /*---------------------------------------------------------------------------- @@ -199,15 +325,15 @@ uint32_t GeoIndexedRaster::getSubsets(const MathLib::extent_t& extent, int64_t g { static_cast(param); - samplingMutex.lock(); + lockSampling(); try { - ssError = SS_NO_ERRORS; - - OGRPolygon poly = GdalRaster::makeRectangle(extent.ll.x, extent.ll.y, extent.ur.x, extent.ur.y); + GroupOrdering groupList; + OGRPolygon poly = GdalRaster::makeRectangle(extent.ll.x, extent.ll.y, extent.ur.x, extent.ur.y); + ssErrors = SS_NO_ERRORS; /* Sample Subsets */ - if(sample(&poly, gps)) + if(sample(&poly, gps, &groupList)) { /* Populate Return Vector of Subsets (slist) */ const GroupOrdering::Iterator iter(groupList); @@ -222,9 +348,9 @@ uint32_t GeoIndexedRaster::getSubsets(const MathLib::extent_t& extent, int64_t g { mlog(e.level(), "Error subsetting raster: %s", e.what()); } - samplingMutex.unlock(); + unlockSampling(); - return ssError; + return ssErrors; } /*---------------------------------------------------------------------------- @@ -232,7 +358,6 @@ uint32_t GeoIndexedRaster::getSubsets(const MathLib::extent_t& extent, int64_t g *----------------------------------------------------------------------------*/ GeoIndexedRaster::~GeoIndexedRaster(void) { - delete [] findersRange; emptyFeaturesList(); } @@ -240,15 +365,40 @@ GeoIndexedRaster::~GeoIndexedRaster(void) * PROTECTED METHODS ******************************************************************************/ +/*---------------------------------------------------------------------------- + * BatchReader Constructor + *----------------------------------------------------------------------------*/ +GeoIndexedRaster::BatchReader::BatchReader (GeoIndexedRaster* _obj): + obj(_obj), + uraster(NULL), + sync(NUM_SYNC_SIGNALS), + run(true) +{ + thread = new Thread(GeoIndexedRaster::batchReaderThread, this); +} + +/*---------------------------------------------------------------------------- + * BatchReader Destructor + *----------------------------------------------------------------------------*/ +GeoIndexedRaster::BatchReader::~BatchReader (void) +{ + sync.lock(); + { + run = false; /* Set run flag to false */ + sync.signal(DATA_TO_SAMPLE, Cond::NOTIFY_ONE); + } + sync.unlock(); + + delete thread; /* delete thread waits on thread to join */ +} + /*---------------------------------------------------------------------------- * Constructor *----------------------------------------------------------------------------*/ GeoIndexedRaster::GeoIndexedRaster(lua_State *L, GeoParms* _parms, GdalRaster::overrideCRS_t cb): RasterObject (L, _parms), cache (MAX_READER_THREADS), - ssError (SS_NO_ERRORS), - numFinders (0), - findersRange (NULL), + ssErrors (SS_NO_ERRORS), crscb (cb), bbox {0, 0, 0, 0}, rows (0), @@ -264,35 +414,122 @@ GeoIndexedRaster::GeoIndexedRaster(lua_State *L, GeoParms* _parms, GdalRaster::o /* Mark index file bbox/extent poly as empty */ geoIndexPoly.empty(); +} + + +/*---------------------------------------------------------------------------- + * getBatchGroupSamples + *----------------------------------------------------------------------------*/ +uint32_t GeoIndexedRaster::getBatchGroupSamples(const rasters_group_t* rgroup, List* slist, uint32_t flags, uint32_t pointIndx) +{ + uint32_t errors = SS_NO_ERRORS; + + for(const auto& rinfo: rgroup->infovect) + { + if(strcmp(VALUE_TAG, rinfo.tag.c_str()) != 0) continue; + + /* This is the unique raster we are looking for, it cannot be NULL */ + unique_raster_t* ur = rinfo.uraster; + assert(ur); + + /* Get the sample for this point from unique raster */ + for(point_sample_t& ps : ur->pointSamples) + { + if(ps.pointInfo.index == pointIndx) + { + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ + if(ps.sample == NULL) break; + + RasterSample* s; - /* Create finder threads used to find rasters intersecting with point/polygon */ - createFinderThreads(); + if(!ps.sampleReturned) + { + ps.sampleReturned = true; + s = ps.sample; + } + else + { + /* Sample has already been returned, must create a copy */ + s = new RasterSample(*ps.sample); + } + + /* Set flags for this sample, add it to the list */ + s->flags = flags; + slist->add(s); + errors |= ps.ssErrors; + + /* + * There is only one raster with VALUE_TAG and one with FLAGS_TAG in a group. + * If group has other type rasters the dataset must override the getGroupFlags method. + */ + return errors; + } + } + } + + return errors; } /*---------------------------------------------------------------------------- - * getGroupSamples + * getBatchGroupFlags *----------------------------------------------------------------------------*/ -void GeoIndexedRaster::getGroupSamples(const rasters_group_t* rgroup, List& slist, uint32_t flags) +uint32_t GeoIndexedRaster::getBatchGroupFlags(const rasters_group_t* rgroup, uint32_t pointIndx) { for(const auto& rinfo: rgroup->infovect) { - if(strcmp(VALUE_TAG, rinfo.tag.c_str()) == 0) + if(strcmp(FLAGS_TAG, rinfo.tag.c_str()) != 0) continue; + + /* This is the unique raster we are looking for, it cannot be NULL */ + unique_raster_t* ur = rinfo.uraster; + assert(ur); + + /* Get the sample for this point from unique raster */ + for(const point_sample_t& ps : ur->pointSamples) { - const char* key = rinfo.fileName.c_str(); - cacheitem_t* item; - if(cache.find(key, &item)) + if(ps.pointInfo.index == pointIndx) { - RasterSample* _sample = item->sample; - if(_sample) + /* sample can be NULL if raster read failed, (e.g. point out of bounds) */ + if(ps.sample) { - _sample->flags = flags; - slist.add(_sample); - item->sample = NULL; + return ps.sample->value; } + } + } + } + + return 0; +} + - /* Get sampling/subset error status */ - ssError |= item->raster->getSSerror(); +/*---------------------------------------------------------------------------- + * getGroupSamples + *----------------------------------------------------------------------------*/ +void GeoIndexedRaster::getGroupSamples(const rasters_group_t* rgroup, List& slist, uint32_t flags) +{ + for(const auto& rinfo: rgroup->infovect) + { + if(strcmp(VALUE_TAG, rinfo.tag.c_str()) != 0) continue; + + const char* key = rinfo.fileName.c_str(); + cacheitem_t* item; + if(cache.find(key, &item)) + { + RasterSample* _sample = item->sample; + if(_sample) + { + _sample->flags = flags; + slist.add(_sample); + item->sample = NULL; } + + /* Get sampling/subset error status */ + ssErrors |= item->raster->getSSerror(); + + /* + * There is only one raster with VALUE_TAG and one with FLAGS_TAG in a group. + * If group has other type rasters the dataset must override the getGroupFlags method. + */ + break; } } } @@ -317,7 +554,7 @@ void GeoIndexedRaster::getGroupSubsets(const rasters_group_t* rgroup, Listraster->getSSerror(); + ssErrors |= item->raster->getSSerror(); } } } @@ -328,27 +565,23 @@ void GeoIndexedRaster::getGroupSubsets(const rasters_group_t* rgroup, Listinfovect) { - if(strcmp(FLAGS_TAG, rinfo.tag.c_str()) == 0) + if(strcmp(FLAGS_TAG, rinfo.tag.c_str()) != 0) continue; + + cacheitem_t* item; + const char* key = rinfo.fileName.c_str(); + if(cache.find(key, &item)) { - cacheitem_t* item; - const char* key = rinfo.fileName.c_str(); - if(cache.find(key, &item)) + const RasterSample* _sample = item->sample; + if(_sample) { - const RasterSample* _sample = item->sample; - if(_sample) - { - flags = _sample->value; - } + return _sample->value; } - break; } } - return flags; + return 0; } /*---------------------------------------------------------------------------- @@ -410,10 +643,10 @@ bool GeoIndexedRaster::getFeatureDate(const OGRFeature* feature, TimeLib::gmt_ti /*---------------------------------------------------------------------------- * openGeoIndex *----------------------------------------------------------------------------*/ -bool GeoIndexedRaster::openGeoIndex(const OGRGeometry* geo) +bool GeoIndexedRaster::openGeoIndex(const OGRGeometry* geo, const std::vector* points) { std::string newFile; - getIndexFile(geo, newFile); + getIndexFile(geo, newFile, points); /* Trying to open the same file? */ if(!featuresList.empty() && newFile == indexFile) @@ -437,6 +670,12 @@ bool GeoIndexedRaster::openGeoIndex(const OGRGeometry* geo) OGRLayer* layer = dset->GetLayer(0); CHECKPTR(layer); + /* If caller provided points, use spatial filter to reduce number of features */ + if(points) + { + applySpatialFilter(layer, points); + } + /* * Clone features and store them for performance/speed of feature lookup */ @@ -480,13 +719,13 @@ bool GeoIndexedRaster::openGeoIndex(const OGRGeometry* geo) } GDALClose((GDALDatasetH)dset); - mlog(DEBUG, "Loaded %lu raster index file", featuresList.size()); + mlog(DEBUG, "Loaded %ld features from: %s", featuresList.size(), newFile.c_str()); } catch (const RunTimeException &e) { if(dset) GDALClose((GDALDatasetH)dset); emptyFeaturesList(); - ssError |= SS_INDEX_FILE_ERROR; + ssErrors |= SS_INDEX_FILE_ERROR; return false; } @@ -540,23 +779,32 @@ void GeoIndexedRaster::sampleRasters(OGRGeometry* geo) /*---------------------------------------------------------------------------- * sample *----------------------------------------------------------------------------*/ -bool GeoIndexedRaster::sample(OGRGeometry* geo, int64_t gps) +bool GeoIndexedRaster::sample(OGRGeometry* geo, int64_t gps, GroupOrdering* groupList) { - /* For AOI always open new index file, for POI it depends... */ - const bool openNewFile = GdalRaster::ispoly(geo) || geoIndexPoly.IsEmpty() || !geoIndexPoly.Contains(geo); + const bool openNewFile = geoIndexPoly.IsEmpty() || !geoIndexPoly.Contains(geo); if(openNewFile) { - if(!openGeoIndex(geo)) + if(!openGeoIndex(geo, NULL)) return false; + } - setFindersRange(); + /* Find rasters that intersect with the geometry */ + finder_t finder(geo->clone(), &featuresList); + if(!findRasters(&finder)) + return false; + + /* Copy finder's raster groups to groupList */ + for(uint32_t j = 0; j < finder.rasterGroups.size(); j++) + { + rasters_group_t* rgroup = finder.rasterGroups[j]; + groupList->add(groupList->length(), rgroup); } - if(!_findRasters(geo)) return false; - if(!filterRasters(gps)) return false; + if(!filterRasters(gps, groupList)) + return false; uint32_t rasters2sample = 0; - if(!updateCache(rasters2sample)) + if(!updateCache(rasters2sample, groupList)) return false; /* Create additional reader threads if needed */ @@ -650,39 +898,6 @@ int GeoIndexedRaster::luaBoundingBox(lua_State *L) return returnLuaStatus(L, status, num_ret); } -/*---------------------------------------------------------------------------- - * finderThread - *----------------------------------------------------------------------------*/ -void* GeoIndexedRaster::finderThread(void *param) -{ - finder_t *finder = static_cast(param); - - while(finder->run) - { - finder->sync.lock(); - { - while((finder->geo == NULL) && finder->run) - finder->sync.wait(DATA_TO_SAMPLE, SYS_TIMEOUT); - } - finder->sync.unlock(); - - if(finder->geo) - { - finder->obj->findRasters(finder); - - finder->sync.lock(); - { - OGRGeometryFactory::destroyGeometry(finder->geo); - finder->geo = NULL; - finder->sync.signal(DATA_SAMPLED, Cond::NOTIFY_ONE); - } - finder->sync.unlock(); - } - } - - return NULL; -} - /*---------------------------------------------------------------------------- * luaCellSize - :cell() --> cell size *----------------------------------------------------------------------------*/ @@ -773,22 +988,259 @@ void* GeoIndexedRaster::readerThread(void *param) } /*---------------------------------------------------------------------------- - * createFinderThreads + * batchReaderThread + *----------------------------------------------------------------------------*/ +void* GeoIndexedRaster::batchReaderThread(void *param) +{ + batch_reader_t *breader = static_cast(param); + + while(breader->run) + { + breader->sync.lock(); + { + /* Wait for raster to work on */ + while((breader->uraster == NULL) && breader->run) + breader->sync.wait(DATA_TO_SAMPLE, SYS_TIMEOUT); + } + breader->sync.unlock(); + + if(breader->uraster != NULL) + { + unique_raster_t* ur = breader->uraster; + GdalRaster* raster = new GdalRaster(breader->obj->parms, + ur->rinfo.fileName, + static_cast(ur->gpsTime / 1000), + ur->fileId, + ur->rinfo.dataIsElevation, + breader->obj->crscb); + + /* Sample all points for this raster */ + for(point_sample_t& ps : ur->pointSamples) + { + ps.sample = raster->samplePOI(&ps.pointInfo.point); + ps.ssErrors |= raster->getSSerror(); + } + + delete raster; + breader->sync.lock(); + { + breader->uraster = NULL; /* Done with this raster and all points */ + breader->sync.signal(DATA_SAMPLED, Cond::NOTIFY_ONE); + } + breader->sync.unlock(); + } + } + + return NULL; +} + +/*---------------------------------------------------------------------------- + * unionThread + *----------------------------------------------------------------------------*/ +void* GeoIndexedRaster::unionThread(void* param) +{ + union_maker_t* um = static_cast(param); + + /* Create an empty geometry collection to hold all points */ + OGRGeometryCollection geometryCollection; + + const double distance = 0.01; /* Aproxiately 1km at equator */ + const double tolerance = distance; /* Tolerance for simplification */ + + /* + * Create geometry collection from the bounding boxes around each point. + * NOTE: Using buffered points is significantly more computationally + * expensive than bounding boxes (10x slower). + */ + mlog(DEBUG, "Creating collection of bboxes from %d points, range: %d - %d", + um->pointsRange.end - um->pointsRange.start, um->pointsRange.start, um->pointsRange.end); + + double startTime = TimeLib::latchtime(); + + const uint32_t start = um->pointsRange.start; + const uint32_t end = um->pointsRange.end; + + for(uint32_t i = start; i < end; i++) + { + const point_info_t& point_info = um->points->at(i); + const double lon = point_info.point.x; + const double lat = point_info.point.y; + + /* Create a linear ring representing the bounding box */ + OGRLinearRing ring; + ring.addPoint(lon - distance, lat - distance); /* Lower-left corner */ + ring.addPoint(lon + distance, lat - distance); /* Lower-right corner */ + ring.addPoint(lon + distance, lat + distance); /* Upper-right corner */ + ring.addPoint(lon - distance, lat + distance); /* Upper-left corner */ + ring.closeRings(); + + /* Create a polygon from the ring */ + OGRPolygon* bboxPoly = new OGRPolygon(); + bboxPoly->addRing(&ring); + + /* Add the polygon to the geometry collection */ + geometryCollection.addGeometry(bboxPoly); + + /* Clean up the bounding box */ + OGRGeometryFactory::destroyGeometry(bboxPoly); + } + um->stats.points2polyTime = TimeLib::latchtime() - startTime; + mlog(DEBUG, "Creating collection took %.3lf seconds", um->stats.points2polyTime); + + /* + * Union the bounding boxes in batches to reduce computational complexity. + * NOTE: Union call fails with large number of geometries - must use batch union. + */ + const int batchSize = 60; /* Sweet spot for performance based on my testing */ + OGRGeometry* unionPolygon = NULL; + + mlog(DEBUG, "Unioning point geometries using batch size: %d", batchSize); + startTime = TimeLib::latchtime(); + for(int i = 0; i < geometryCollection.getNumGeometries(); i += batchSize) + { + if(!um->obj->isSampling()) + { + mlog(WARNING, "Sampling has been stopped, exiting union thread"); + OGRGeometryFactory::destroyGeometry(unionPolygon); + return NULL; + } + + OGRGeometry* batchUnion = NULL; + for(int j = i; j < std::min(i + batchSize, geometryCollection.getNumGeometries()); ++j) + { + OGRGeometry* geo = geometryCollection.getGeometryRef(j); + if(batchUnion == NULL) + { + batchUnion = geo->clone(); + } + else + { + OGRGeometry* newUnion = batchUnion->Union(geo); + OGRGeometryFactory::destroyGeometry(batchUnion); + batchUnion = newUnion; + } + } + + /* Combine the batch union with the overall union polygon */ + if(batchUnion) + { + /* Simplify the batch union to reduce complexity */ + OGRGeometry* simplifiedBatchUnion = batchUnion->Simplify(tolerance); + OGRGeometryFactory::destroyGeometry(batchUnion); + batchUnion = simplifiedBatchUnion; + + if(unionPolygon == NULL) + { + /* Initial union polygon */ + unionPolygon = batchUnion; + } + else + { + OGRGeometry* newUnionPolygon = unionPolygon->Union(batchUnion); + OGRGeometryFactory::destroyGeometry(unionPolygon); + unionPolygon = newUnionPolygon; + OGRGeometryFactory::destroyGeometry(batchUnion); + } + } + } + + if(unionPolygon == NULL) + { + mlog(ERROR, "Unioning geometries failed"); + um->stats.unioningTime = TimeLib::latchtime() - startTime; + return NULL; + } + + /* Simplify the final union polygon to reduce complexity */ + OGRGeometry* simplifiedPolygon = unionPolygon->Simplify(tolerance); + OGRGeometryFactory::destroyGeometry(unionPolygon); + unionPolygon = simplifiedPolygon; + if(simplifiedPolygon == NULL) + { + mlog(ERROR, "Failed to simplify polygon"); + } + + um->stats.unioningTime = TimeLib::latchtime() - startTime; + mlog(DEBUG, "Unioning all geometries took %.3lf seconds", um->stats.unioningTime); + + + /* Set the unioned polygon in the union maker object */ + um->unionPolygon = unionPolygon; + + /* Exit Thread */ + return NULL; +} + +/*---------------------------------------------------------------------------- + * groupsFinderThread *----------------------------------------------------------------------------*/ -bool GeoIndexedRaster::createFinderThreads(void) +void* GeoIndexedRaster::groupsFinderThread(void *param) { - /* Finder threads are created in the constructor, this call should not fail */ + groups_finder_t* gf = static_cast(param); + + /* Find groups of rasters for each point in range. + * NOTE: cannot filter rasters here, must be done in one thread + */ + const uint32_t start = gf->pointsRange.start; + const uint32_t end = gf->pointsRange.end; - for(int i = 0; i < MAX_FINDER_THREADS; i++) + /* Clone features list for this thread (OGRFeature objects are not thread save, must copy/clone) */ + const uint32_t size = gf->obj->featuresList.size(); + std::vector localFeaturesList; + for(uint32_t j = 0; j < size; j++) { - Finder* f = new Finder(this); - finders.add(f); + localFeaturesList.push_back(gf->obj->featuresList[j]->Clone()); } - /* Array of ranges for each thread */ - findersRange = new finder_range_t[MAX_FINDER_THREADS]; + for(uint32_t i = start; i < end; i++) + { + if(!gf->obj->isSampling()) + { + mlog(WARNING, "Sampling has been stopped, exiting groups finder thread"); + break; + } - return finders.length() == MAX_FINDER_THREADS; + const point_info_t& pinfo = gf->points->at(i); + OGRPoint* ogr_point = new OGRPoint(pinfo.point.x, pinfo.point.y, pinfo.point.z); + + GroupOrdering* groupList = new GroupOrdering(); + + /* Set finder for the whole range of features */ + Finder finder(ogr_point, &localFeaturesList); + + /* Find rasters intersecting with ogr_point */ + gf->obj->findRasters(&finder); + + /* Filter rasters based on gps time */ + const int64_t gps = gf->obj->usePOItime() ? pinfo.gps : 0.0; + gf->obj->filterRasters(gps, groupList); + + /* Copy from finder to pointGroups */ + for(rasters_group_t* rgroup : finder.rasterGroups) + { + groupList->add(groupList->length(), rgroup); + } + gf->pointsGroups.emplace_back(point_groups_t{ {*ogr_point, i}, groupList }); + + /* Add raster file names from this groupList to raster to points map */ + const GroupOrdering::Iterator iter(*groupList); + for(int j = 0; j < iter.length; j++) + { + const rasters_group_t* rgroup = iter[j].value; + for(const raster_info_t& rinfo : rgroup->infovect) + { + gf->rasterToPointsMap[rinfo.fileName].insert(i); + } + } + } + + /* Cleanup cloned features */ + for(OGRFeature* feature : localFeaturesList) + { + OGRFeature::DestroyFeature(feature); + } + + return NULL; } /*---------------------------------------------------------------------------- @@ -813,17 +1265,46 @@ bool GeoIndexedRaster::createReaderThreads(uint32_t rasters2sample) } catch (const RunTimeException &e) { - ssError |= SS_RESOURCE_LIMIT_ERROR; + ssErrors |= SS_RESOURCE_LIMIT_ERROR; mlog(CRITICAL, "Failed to create reader threads, needed: %d, created: %d", newThreadsCnt, readers.length() - threadsNow); } return readers.length() == threadsNeeded; } +/*---------------------------------------------------------------------------- + * createBatchReaderThreads + *----------------------------------------------------------------------------*/ +bool GeoIndexedRaster::createBatchReaderThreads(uint32_t rasters2sample) +{ + const int threadsNeeded = rasters2sample; + const int threadsNow = batchReaders.length(); + const int newThreadsCnt = threadsNeeded - threadsNow; + + if(threadsNeeded <= threadsNow) + return true; + + try + { + for(int i = 0; i < newThreadsCnt; i++) + { + BatchReader* r = new BatchReader(this); + batchReaders.add(r); + } + } + catch (const RunTimeException &e) + { + ssErrors |= SS_RESOURCE_LIMIT_ERROR; + mlog(CRITICAL, "Failed to create batch reader threads"); + } + + return batchReaders.length() == threadsNeeded; +} + /*---------------------------------------------------------------------------- * updateCache *----------------------------------------------------------------------------*/ -bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample) +bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample, const GroupOrdering* groupList) { /* Mark all items in cache as not enabled */ { @@ -837,7 +1318,7 @@ bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample) } /* Cache contains items/rasters from previous sample run */ - const GroupOrdering::Iterator group_iter(groupList); + const GroupOrdering::Iterator group_iter(*groupList); for(int i = 0; i < group_iter.length; i++) { const rasters_group_t* rgroup = group_iter[i].value; @@ -896,7 +1377,7 @@ bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample) /* Check for max limit of concurent reading raster threads */ if(rasters2sample > MAX_READER_THREADS) { - ssError |= SS_THREADS_LIMIT_ERROR; + ssErrors |= SS_THREADS_LIMIT_ERROR; mlog(ERROR, "Too many rasters to read: %d, max allowed: %d", cache.length(), MAX_READER_THREADS); return false; } @@ -907,12 +1388,12 @@ bool GeoIndexedRaster::updateCache(uint32_t& rasters2sample) /*---------------------------------------------------------------------------- * filterRasters *----------------------------------------------------------------------------*/ -bool GeoIndexedRaster::filterRasters(int64_t gps) +bool GeoIndexedRaster::filterRasters(int64_t gps, GroupOrdering* groupList) { /* NOTE: temporal filter is applied in openGeoIndex() */ if(parms->url_substring || parms->filter_doy_range) { - const GroupOrdering::Iterator group_iter(groupList); + const GroupOrdering::Iterator group_iter(*groupList); for(int i = 0; i < group_iter.length; i++) { const rasters_group_t* rgroup = group_iter[i].value; @@ -955,7 +1436,7 @@ bool GeoIndexedRaster::filterRasters(int64_t gps) if(removeGroup) { - groupList.remove(group_iter[i].key); + groupList->remove(group_iter[i].key); } } } @@ -978,7 +1459,7 @@ bool GeoIndexedRaster::filterRasters(int64_t gps) int64_t minDelta = abs(std::numeric_limits::max() - closestGps); /* Find raster group with the closest time */ - const GroupOrdering::Iterator group_iter(groupList); + const GroupOrdering::Iterator group_iter(*groupList); for(int i = 0; i < group_iter.length; i++) { const rasters_group_t* rgroup = group_iter[i].value; @@ -998,104 +1479,413 @@ bool GeoIndexedRaster::filterRasters(int64_t gps) if(delta > minDelta) { - groupList.remove(group_iter[i].key); + groupList->remove(group_iter[i].key); } } } - return (!groupList.empty()); + return (!groupList->empty()); } /*---------------------------------------------------------------------------- - * setFindersRange + * getConvexHull *----------------------------------------------------------------------------*/ -void GeoIndexedRaster::setFindersRange(void) +OGRGeometry* GeoIndexedRaster::getConvexHull(const std::vector* points) { - const uint32_t minFeaturesPerThread = MIN_FEATURES_PER_FINDER_THREAD; - const uint32_t features = featuresList.size(); + if(points->empty()) + return NULL; + + /* Create an empty geometry collection to hold all points */ + OGRGeometryCollection geometryCollection; + + mlog(INFO, "Creating convex hull from %zu points", points->size()); + + /* Collect all points into a geometry collection */ + for(const point_info_t& point_info : *points) + { + const double lon = point_info.point.x; + const double lat = point_info.point.y; + + OGRPoint* point = new OGRPoint(lon, lat); + geometryCollection.addGeometryDirectly(point); + } - /* Determine how many finder threads to use and index range for each */ - if(features <= minFeaturesPerThread) + /* Create a convex hull that wraps around all the points */ + OGRGeometry* convexHull = geometryCollection.ConvexHull(); + if(convexHull == NULL) { - numFinders = 1; - findersRange[0].start_indx = 0; - findersRange[0].end_indx = features; - return; + mlog(ERROR, "Failed to create a convex hull around points."); } - numFinders = std::min(static_cast(MAX_FINDER_THREADS), features / minFeaturesPerThread); + return convexHull; +} - /* Ensure at least two threads if features > minFeaturesPerThread */ - if(numFinders == 1) +/*---------------------------------------------------------------------------- + * getBufferedPoints + *----------------------------------------------------------------------------*/ +OGRGeometry* GeoIndexedRaster::getBufferedPoints(const std::vector* points) +{ + if(points->empty()) + return NULL; + + /* Start geometry union threads */ + std::vector pids; + std::vector unionMakers; + + const uint32_t numMaxThreads = std::thread::hardware_concurrency(); + const uint32_t minPointsPerThread = 100; + + std::vector pointsRanges; + getThreadsRanges(pointsRanges, points->size(), minPointsPerThread, numMaxThreads); + mlog(DEBUG, "Using %ld groupfinder threads to search for %ld points", pointsRanges.size(), points->size()); + + const uint32_t numThreads = pointsRanges.size(); + const double startTime = TimeLib::latchtime(); + + for(uint32_t i = 0; i < numThreads; i++) + { + union_maker_t* um = new union_maker_t(this, points); + um->pointsRange = pointsRanges[i]; + unionMakers.push_back(um); + Thread* pid = new Thread(unionThread, um); + pids.push_back(pid); + } + + /* Wait for all union maker threads to finish */ + for(Thread* pid : pids) { - numFinders = 2; + delete pid; } - const uint32_t featuresPerThread = features / numFinders; - uint32_t remainingFeatures = features % numFinders; + /* Combine results from all union maker threads */ + OGRGeometry* unionPolygon = NULL; - uint32_t start = 0; - for(uint32_t i = 0; i < numFinders; i++) + for(union_maker_t* um : unionMakers) { - findersRange[i].start_indx = start; - findersRange[i].end_indx = start + featuresPerThread + (remainingFeatures > 0 ? 1 : 0); - start = findersRange[i].end_indx; - if(remainingFeatures > 0) + if(um->unionPolygon == NULL) + { + mlog(WARNING, "Union polygon is NULL"); + continue; + } + + if(unionPolygon == NULL) + { + unionPolygon = um->unionPolygon; + } + else + { + OGRGeometry* newUnionPolygon = unionPolygon->Union(um->unionPolygon); + OGRGeometryFactory::destroyGeometry(unionPolygon); + unionPolygon = newUnionPolygon; + OGRGeometryFactory::destroyGeometry(um->unionPolygon); + } + + const double tolerance = 0.0005; /* Simplification tolerance */ + OGRGeometry* simplifiedPolygon = unionPolygon->Simplify(tolerance); + OGRGeometryFactory::destroyGeometry(unionPolygon); + unionPolygon = simplifiedPolygon; + if(simplifiedPolygon == NULL) { - remainingFeatures--; + mlog(ERROR, "Failed to simplify polygon"); } } + + const double elapsedTime = TimeLib::latchtime() - startTime; + mlog(DEBUG, "Unioning point geometries took %.3lf", elapsedTime); + + /* Log stats and cleanup */ + for(uint32_t i = 0; i < unionMakers.size(); i++) + { + union_maker_t* um = unionMakers[i]; + + mlog(DEBUG, "Thread[%d]: %10d - %-10d p2poly: %.3lf, unioning: %.3lf", + i, um->pointsRange.start, um->pointsRange.end, um->stats.points2polyTime, um->stats.unioningTime); + + delete um; + } + + return unionPolygon; } +/*---------------------------------------------------------------------------- + * applySpatialFilter + *----------------------------------------------------------------------------*/ +void GeoIndexedRaster::applySpatialFilter(OGRLayer* layer, const std::vector* points) +{ + mlog(INFO, "Features before spatial filter: %lld", layer->GetFeatureCount()); + + const double startTime = TimeLib::latchtime(); + + OGRGeometry* filter = getBufferedPoints(points); + if(filter != NULL) + { + layer->SetSpatialFilter(filter); + OGRGeometryFactory::destroyGeometry(filter); + } + else + { + mlog(ERROR, "Failed to create union polygon for spatial filter"); + } + perfStats.spatialFilterTime = TimeLib::latchtime() - startTime; + + mlog(INFO, "Features after spatial filter: %lld", layer->GetFeatureCount()); + mlog(INFO, "Spatial filter time: %.3lf", perfStats.spatialFilterTime); +} /*---------------------------------------------------------------------------- - * _findRasters + * findAllGroups *----------------------------------------------------------------------------*/ -bool GeoIndexedRaster::_findRasters(OGRGeometry* geo) +bool GeoIndexedRaster::findAllGroups(const std::vector* points, + std::vector& pointsGroups, + raster_points_map_t& rasterToPointsMap) { - groupList.clear(); + bool status = false; + const double startTime = TimeLib::latchtime(); - /* Start finder threads to find rasters intersecting with point/polygon */ - uint32_t signaledFinders = 0; - for(uint32_t i = 0; i < numFinders; i++) + try { - Finder* finder = finders[i]; - finder->sync.lock(); + mlog(DEBUG, "Finding rasters groups for all points"); + + /* Start rasters groups finder threads */ + std::vector pids; + std::vector rgroupFinders; + + const uint32_t numMaxThreads = std::thread::hardware_concurrency(); + const uint32_t minPointsPerThread = 100; + + std::vector pointsRanges; + getThreadsRanges(pointsRanges, points->size(), minPointsPerThread, numMaxThreads); + const uint32_t numThreads = pointsRanges.size(); + + for(uint32_t i = 0; i < numThreads; i++) + { + GroupsFinder* gf = new GroupsFinder(this, points); + gf->pointsRange = pointsRanges[i]; + rgroupFinders.push_back(gf); + Thread* pid = new Thread(groupsFinderThread, gf); + pids.push_back(pid); + } + + /* Wait for all groups finder threads to finish */ + for(Thread* pid : pids) { - OGRGeometryFactory::destroyGeometry(finder->geo); - finder->geo = geo->clone(); - finder->range = findersRange[i]; - finder->rasterGroups.clear(); - finder->sync.signal(DATA_TO_SAMPLE, Cond::NOTIFY_ONE); - signaledFinders++; + delete pid; } - finder->sync.unlock(); + + /* Merge the pointGroups for each thread */ + for(GroupsFinder* gf : rgroupFinders) + { + pointsGroups.insert(pointsGroups.end(), gf->pointsGroups.begin(), gf->pointsGroups.end()); + + /* Merge the rasterToPointsMap for each thread */ + for(const raster_points_map_t::value_type& pair : gf->rasterToPointsMap) + { + rasterToPointsMap[pair.first].insert(pair.second.begin(), pair.second.end()); + } + + delete gf; + } + + /* Verify that the number of points groups is the same as the number of points */ + if(pointsGroups.size() != points->size()) + { + mlog(ERROR, "Number of points groups: %zu does not match number of points: %zu", pointsGroups.size(), points->size()); + throw RunTimeException(CRITICAL, RTE_ERROR, "Number of points groups does not match number of points"); + } + + status = true; } + catch (const RunTimeException &e) + { + mlog(e.level(), "Error creating groups: %s", e.what()); + } + + perfStats.findRastersTime = TimeLib::latchtime() - startTime; + return status; +} + + +/*---------------------------------------------------------------------------- + * findUniqueRasters + *----------------------------------------------------------------------------*/ +bool GeoIndexedRaster::findUniqueRasters(std::vector& uniqueRasters, + const std::vector& pointsGroups, + raster_points_map_t& rasterToPointsMap) +{ + bool status = false; + const double startTime = TimeLib::latchtime(); - /* Wait for finder threads to finish searching for rasters */ - for(uint32_t i = 0; i < signaledFinders; i++) + try { - finder_t* finder = finders[i]; - finder->sync.lock(); + /* Create vector of unique rasters. */ + mlog(DEBUG, "Finding unique rasters"); + for(const point_groups_t& pg : pointsGroups) + { + const GroupOrdering::Iterator iter(*pg.groupList); + for(int64_t i = 0; i < iter.length; i++) + { + rasters_group_t* rgroup = iter[i].value; + for(raster_info_t& rinfo : rgroup->infovect) + { + /* Is this raster already in the list of unique rasters? */ + bool addNewRaster = true; + for(unique_raster_t* ur : uniqueRasters) + { + if(ur->rinfo.fileName == rinfo.fileName) + { + /* already in unique rasters list, set pointer in rinfo to this raster */ + rinfo.uraster = ur; + addNewRaster = false; + break; + } + } + + if(addNewRaster) + { + unique_raster_t* ur = new unique_raster_t(); + ur->rinfo = rinfo; + ur->gpsTime = rgroup->gpsTime; + ur->fileId = fileDictAdd(rinfo.fileName); + uniqueRasters.push_back(ur); + + /* Set pointer in rinfo to this new unique raster */ + rinfo.uraster = ur; + } + } + } + } + + /* + * For each unique raster, find the points that belong to it + */ + mlog(DEBUG, "Finding points for unique rasters"); + for(unique_raster_t* ur : uniqueRasters) { - while(finder->geo != NULL) - finder->sync.wait(DATA_SAMPLED, SYS_TIMEOUT); + const std::string& rasterName = ur->rinfo.fileName; + + auto it = rasterToPointsMap.find(rasterName); + if(it != rasterToPointsMap.end()) + { + for(const uint32_t pointIndx : it->second) + { + const point_groups_t& pg = pointsGroups[pointIndx]; + ur->pointSamples.push_back({ pg.pointInfo, NULL, false, SS_NO_ERRORS }); + } + } } - finder->sync.unlock(); + + status = true; } + catch(const RunTimeException& e) + { + mlog(e.level(), "Error creating groups: %s", e.what()); + } + + perfStats.findUniqueRastersTime = TimeLib::latchtime() - startTime; + return status; +} + - /* Combine results from all finder threads */ - for(uint32_t i = 0; i < numFinders; i++) +/*---------------------------------------------------------------------------- + * sampleUniqueRasters + *----------------------------------------------------------------------------*/ +bool GeoIndexedRaster::sampleUniqueRasters(const std::vector& uniqueRasters) +{ + bool status = false; + const double startTime = TimeLib::latchtime(); + + try { - const Finder* finder = finders[i]; - for(uint32_t j = 0; j < finder->rasterGroups.size(); j++) + /* Testing has shown that 20 threads performs twice as fast on a 8 core system than 50 or 100 threads. */ + const uint32_t maxThreads = 20; + + /* Create batch reader threads */ + const uint32_t numRasters = uniqueRasters.size(); + createBatchReaderThreads(std::min(maxThreads, numRasters)); + + const uint32_t numThreads = batchReaders.length(); + mlog(DEBUG, "Sampling %u rasters with %u threads", numRasters, numThreads); + + /* Sample unique rasters utilizing numThreads */ + uint32_t currentRaster = 0; + + while(currentRaster < numRasters) { - rasters_group_t* rgroup = finder->rasterGroups[j]; - groupList.add(groupList.length(), rgroup); + /* Calculate how many rasters we can process in this batch */ + const uint32_t batchSize = std::min(numThreads, numRasters - currentRaster); + + /* Keep track of how many threads have been assigned work */ + uint32_t activeReaders = 0; + + /* Assign rasters to batch readers as soon as they are free */ + while(currentRaster < numRasters || activeReaders > 0) + { + for(uint32_t i = 0; i < batchSize; i++) + { + BatchReader* breader = batchReaders[i]; + + breader->sync.lock(); + { + /* If this thread is done with its previous raster, assign a new one */ + if(breader->uraster == NULL && currentRaster < numRasters) + { + breader->uraster = uniqueRasters[currentRaster++]; + breader->sync.signal(DATA_TO_SAMPLE, Cond::NOTIFY_ONE); + activeReaders++; + } + } + breader->sync.unlock(); + + if(!isSampling()) + { + /* Sampling has been stopped, stop assigning new rasters */ + activeReaders = 0; + currentRaster = numRasters; + break; + } + + /* Check if the current breader has completed its work */ + breader->sync.lock(); + { + if(breader->uraster == NULL && activeReaders > 0) + { + /* Mark one reader as free */ + activeReaders--; + } + } + breader->sync.unlock(); + } + + /* Short wait before checking again to avoid busy waiting */ + std::this_thread::sleep_for(std::chrono::milliseconds(10)); + } } + + /* Wait for all batch readers to finish sampling */ + for(int32_t i = 0; i < batchReaders.length(); i++) + { + BatchReader* breader = batchReaders[i]; + + breader->sync.lock(); + { + while(breader->uraster != NULL) + breader->sync.wait(DATA_SAMPLED, SYS_TIMEOUT); + } + breader->sync.unlock(); + } + + status = true; + } + catch(const RunTimeException& e) + { + mlog(e.level(), "Error creating groups: %s", e.what()); } - return (!groupList.empty()); + perfStats.samplesTime = TimeLib::latchtime() - startTime; + return status; } + + + diff --git a/packages/geo/GeoIndexedRaster.h b/packages/geo/GeoIndexedRaster.h index 0b36ef0eb..5bb86a5e5 100644 --- a/packages/geo/GeoIndexedRaster.h +++ b/packages/geo/GeoIndexedRaster.h @@ -39,6 +39,8 @@ #include "GdalRaster.h" #include "RasterObject.h" #include "Ordering.h" +#include +#include /****************************************************************************** @@ -56,9 +58,6 @@ class GeoIndexedRaster: public RasterObject static const int MAX_CACHE_SIZE = 20; static const int MAX_READER_THREADS = 200; - static const int MAX_FINDER_THREADS = 16; - static const int MIN_FEATURES_PER_FINDER_THREAD = 20; - static const char* FLAGS_TAG; static const char* VALUE_TAG; static const char* DATE_TAG; @@ -67,14 +66,36 @@ class GeoIndexedRaster: public RasterObject * Typedefs *--------------------------------------------------------------------*/ + typedef struct { + OGRPoint point; + int64_t index; + } ogr_point_info_t; + + typedef struct { + ogr_point_info_t pointInfo; + RasterSample* sample; + bool sampleReturned; + uint32_t ssErrors; + } point_sample_t; + + struct unique_raster_t; typedef struct RasterInfo { - bool dataIsElevation; - std::string tag; - std::string fileName; + bool dataIsElevation; + std::string tag; + std::string fileName; + unique_raster_t* uraster; // only used for batch reading - RasterInfo(void): dataIsElevation(false) {} + RasterInfo(void): dataIsElevation(false), uraster(NULL) {} } raster_info_t; + typedef struct unique_raster_t { + raster_info_t rinfo; + double gpsTime; + uint64_t fileId; + std::vector pointSamples; + } unique_raster_t; + + typedef struct RaserGroup { std::string id; std::vector infovect; @@ -84,6 +105,19 @@ class GeoIndexedRaster: public RasterObject RaserGroup(void): gmtDate{0,0,0,0,0,0}, gpsTime(0) {} } rasters_group_t; + typedef Ordering GroupOrdering; + + typedef struct BatchReader { + GeoIndexedRaster* obj; + unique_raster_t* uraster; + Thread* thread; + Cond sync; + bool run; + explicit BatchReader(GeoIndexedRaster* _obj); + ~BatchReader(void); + } batch_reader_t; + + typedef struct CacheItem { bool enabled; RasterSample* sample; @@ -99,27 +133,53 @@ class GeoIndexedRaster: public RasterObject cacheitem_t* entry; Cond sync; bool run; - explicit Reader(GeoIndexedRaster* raster); + explicit Reader(GeoIndexedRaster* _obj); ~Reader(void); } reader_t; - typedef struct { - uint32_t start_indx; - uint32_t end_indx; - } finder_range_t; - + typedef RasterObject::range_t range_t; typedef struct Finder { - GeoIndexedRaster* obj; OGRGeometry* geo; - finder_range_t range; + std::vector* featuresList; std::vector rasterGroups; - Thread* thread; - Cond sync; - bool run; - explicit Finder(GeoIndexedRaster* raster); - ~Finder(void); + explicit Finder(OGRGeometry* geo, std::vector* _featuresList); + ~Finder(void); } finder_t; + + typedef struct { + double points2polyTime; + double unioningTime; + } union_maker_stats_t; + + typedef struct UnionMaker { + GeoIndexedRaster* obj; + range_t pointsRange; + const std::vector* points; + OGRGeometry* unionPolygon; + union_maker_stats_t stats; + + explicit UnionMaker (GeoIndexedRaster* _obj, const std::vector* _points); + } union_maker_t; + + typedef struct { + ogr_point_info_t pointInfo; + GroupOrdering* groupList; + } point_groups_t; + + /* Typedef for the global map (raster file name -> set of unique point IDs) */ + typedef std::unordered_map> raster_points_map_t; + + typedef struct GroupsFinder { + GeoIndexedRaster* obj; + range_t pointsRange; + const std::vector* points; + std::vector pointsGroups; + raster_points_map_t rasterToPointsMap; + + explicit GroupsFinder (GeoIndexedRaster* _obj, const std::vector* _points); + } groups_finder_t; + /*-------------------------------------------------------------------- * Methods *--------------------------------------------------------------------*/ @@ -127,48 +187,62 @@ class GeoIndexedRaster: public RasterObject static void init (void); static void deinit (void); uint32_t getSamples (const MathLib::point_3d_t& point, int64_t gps, List& slist, void* param=NULL) final; + uint32_t getSamples (const std::vector& points, List& sllist, void* param=NULL) final; uint32_t getSubsets (const MathLib::extent_t& extent, int64_t gps, List& slist, void* param=NULL) final; ~GeoIndexedRaster (void) override; protected: /*-------------------------------------------------------------------- - * Methods + * Types *--------------------------------------------------------------------*/ - GeoIndexedRaster (lua_State* L, GeoParms* _parms, GdalRaster::overrideCRS_t cb=NULL); - virtual void getGroupSamples (const rasters_group_t* rgroup, List& slist, uint32_t flags); - virtual void getGroupSubsets (const rasters_group_t* rgroup, List& slist); - uint32_t getGroupFlags (const rasters_group_t* rgroup); - static double getGmtDate (const OGRFeature* feature, const char* field, TimeLib::gmt_time_t& gmtDate); - virtual bool openGeoIndex (const OGRGeometry* geo); - virtual bool getFeatureDate (const OGRFeature* feature, TimeLib::gmt_time_t& gmtDate); - virtual void getIndexFile (const OGRGeometry* geo, std::string& file) = 0; - virtual bool findRasters (finder_t* finder) = 0; - void sampleRasters (OGRGeometry* geo); - bool sample (OGRGeometry* geo, int64_t gps); - void emptyFeaturesList (void); + typedef Dictionary CacheDictionary; /*-------------------------------------------------------------------- - * Types + * Methods *--------------------------------------------------------------------*/ - typedef Dictionary CacheDictionary; - typedef Ordering GroupOrdering; + GeoIndexedRaster (lua_State* L, GeoParms* _parms, GdalRaster::overrideCRS_t cb=NULL); + virtual uint32_t getBatchGroupSamples (const rasters_group_t* rgroup, List* slist, uint32_t flags, uint32_t pointIndx); + static uint32_t getBatchGroupFlags (const rasters_group_t* rgroup, uint32_t pointIndx); + + virtual void getGroupSamples (const rasters_group_t* rgroup, List& slist, uint32_t flags); + virtual void getGroupSubsets (const rasters_group_t* rgroup, List& slist); + uint32_t getGroupFlags (const rasters_group_t* rgroup); + + static double getGmtDate (const OGRFeature* feature, const char* field, TimeLib::gmt_time_t& gmtDate); + bool openGeoIndex (const OGRGeometry* geo, const std::vector* points); + virtual bool getFeatureDate (const OGRFeature* feature, TimeLib::gmt_time_t& gmtDate); + virtual void getIndexFile (const OGRGeometry* geo, std::string& file, const std::vector* points=NULL) = 0; + virtual bool findRasters (finder_t* finder) = 0; + void sampleRasters (OGRGeometry* geo); + bool sample (OGRGeometry* geo, int64_t gps, GroupOrdering* groupList); + void emptyFeaturesList (void); /*-------------------------------------------------------------------- * Data *--------------------------------------------------------------------*/ - Mutex samplingMutex; - GroupOrdering groupList; - CacheDictionary cache; - vector featuresList; - OGRPolygon geoIndexPoly; - uint32_t ssError; + CacheDictionary cache; + uint32_t ssErrors; private: + /*-------------------------------------------------------------------- + * Types + *--------------------------------------------------------------------*/ + + typedef struct PerfStats { + double spatialFilterTime; + double findRastersTime; + double findUniqueRastersTime; + double samplesTime; + + PerfStats(void) : spatialFilterTime(0), findRastersTime(0), findUniqueRastersTime(0), samplesTime(0) {} + void clear(void) { spatialFilterTime = 0; findRastersTime = 0; findUniqueRastersTime = 0; samplesTime = 0;} + } perf_stats_t; + /*-------------------------------------------------------------------- * Constants *--------------------------------------------------------------------*/ @@ -181,10 +255,12 @@ class GeoIndexedRaster: public RasterObject * Data *--------------------------------------------------------------------*/ - uint32_t numFinders; - finder_range_t* findersRange; - List finders; + std::vector featuresList; + OGRPolygon geoIndexPoly; + List readers; + List batchReaders; + perf_stats_t perfStats; GdalRaster::overrideCRS_t crscb; std::string indexFile; @@ -196,19 +272,33 @@ class GeoIndexedRaster: public RasterObject * Methods *--------------------------------------------------------------------*/ - static int luaDimensions (lua_State* L); - static int luaBoundingBox (lua_State* L); - static int luaCellSize (lua_State* L); + static int luaDimensions (lua_State* L); + static int luaBoundingBox (lua_State* L); + static int luaCellSize (lua_State* L); - static void* finderThread (void *param); - static void* readerThread (void *param); + static void* readerThread (void *param); + static void* batchReaderThread (void *param); + static void* unionThread (void* param); + static void* groupsFinderThread (void *param); - bool createFinderThreads (void); bool createReaderThreads (uint32_t rasters2sample); - bool updateCache (uint32_t& rasters2sample); - bool filterRasters (int64_t gps); - void setFindersRange (void); - bool _findRasters (OGRGeometry* geo); + bool createBatchReaderThreads(uint32_t rasters2sample); + + bool updateCache (uint32_t& rasters2sample, const GroupOrdering* groupList); + bool filterRasters (int64_t gps, GroupOrdering* groupList); + static OGRGeometry* getConvexHull (const std::vector* points); + OGRGeometry* getBufferedPoints (const std::vector* points); + void applySpatialFilter (OGRLayer* layer, const std::vector* points); + + bool findAllGroups (const std::vector* points, + std::vector& pointsGroups, + raster_points_map_t& rasterToPointsMap); + + bool findUniqueRasters (std::vector& uniqueRasters, + const std::vector& pointsGroups, + raster_points_map_t& rasterToPointsMap); + + bool sampleUniqueRasters (const std::vector& uniqueRasters); }; #endif /* __geo_indexed_raster__ */ diff --git a/packages/geo/GeoParms.cpp b/packages/geo/GeoParms.cpp index 16c152763..60f295893 100644 --- a/packages/geo/GeoParms.cpp +++ b/packages/geo/GeoParms.cpp @@ -398,6 +398,61 @@ const char* GeoParms::tojson(void) const return StringLib::duplicate(buffer.GetString()); } +/*---------------------------------------------------------------------------- + * sserror2str + *----------------------------------------------------------------------------*/ +std::string GeoParms::sserror2str(uint32_t error) +{ + std::string errorStr; + + if(error == SS_NO_ERRORS) + { + errorStr = "SS_NO_ERRORS"; + return errorStr; + } + + if(error & SS_THREADS_LIMIT_ERROR) + { + errorStr += "SS_THREADS_LIMIT_ERROR, "; + } + if(error & SS_MEMPOOL_ERROR) + { + errorStr += "SS_MEMPOOL_ERROR, "; + } + if(error & SS_OUT_OF_BOUNDS_ERROR) + { + errorStr += "SS_OUT_OF_BOUNDS_ERROR, "; + } + if(error & SS_READ_ERROR) + { + errorStr += "SS_READ_ERROR, "; + } + if(error & SS_WRITE_ERROR) + { + errorStr += "SS_WRITE_ERROR, "; + } + if(error & SS_SUBRASTER_ERROR) + { + errorStr += "SS_SUBRASTER_ERROR, "; + } + if(error & SS_INDEX_FILE_ERROR) + { + errorStr += "SS_INDEX_FILE_ERROR, "; + } + if(error & SS_RESOURCE_LIMIT_ERROR) + { + errorStr += "SS_RESOURCE_LIMIT_ERROR, "; + } + + /* Remove the last ", " if it exists */ + if (errorStr.size() >= 2 && errorStr[errorStr.size() - 2] == ',') + { + errorStr.resize(errorStr.size() - 2); + } + return errorStr; +} + + /*---------------------------------------------------------------------------- * cleanup *----------------------------------------------------------------------------*/ diff --git a/packages/geo/GeoParms.h b/packages/geo/GeoParms.h index 9350b5c89..4156f06bd 100644 --- a/packages/geo/GeoParms.h +++ b/packages/geo/GeoParms.h @@ -155,6 +155,8 @@ class GeoParms: public LuaObject ~GeoParms (void) override; const char* tojson (void) const override; + static std::string sserror2str(uint32_t error); + private: /*-------------------------------------------------------------------- diff --git a/packages/geo/GeoRaster.cpp b/packages/geo/GeoRaster.cpp index a60cede5b..0c9af9c2e 100644 --- a/packages/geo/GeoRaster.cpp +++ b/packages/geo/GeoRaster.cpp @@ -83,7 +83,7 @@ uint32_t GeoRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, Li static_cast(gps); static_cast(param); - samplingMutex.lock(); + lockSampling(); try { OGRPoint ogr_point(point.x, point.y, point.z); @@ -94,7 +94,7 @@ uint32_t GeoRaster::getSamples(const MathLib::point_3d_t& point, int64_t gps, Li { mlog(e.level(), "Error getting samples: %s", e.what()); } - samplingMutex.unlock(); + unlockSampling(); return raster.getSSerror(); } @@ -107,7 +107,7 @@ uint32_t GeoRaster::getSubsets(const MathLib::extent_t& extent, int64_t gps, Lis static_cast(gps); static_cast(param); - samplingMutex.lock(); + lockSampling(); /* Enable multi-threaded decompression in Gtiff driver */ CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); @@ -146,7 +146,7 @@ uint32_t GeoRaster::getSubsets(const MathLib::extent_t& extent, int64_t gps, Lis /* Disable multi-threaded decompression in Gtiff driver */ CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "1"); - samplingMutex.unlock(); + unlockSampling(); return raster.getSSerror(); } @@ -159,7 +159,7 @@ uint8_t* GeoRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t xsize, uint32 static_cast(param); uint8_t* data = NULL; - samplingMutex.lock(); + lockSampling(); /* Enable multi-threaded decompression in Gtiff driver */ CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); @@ -169,7 +169,7 @@ uint8_t* GeoRaster::getPixels(uint32_t ulx, uint32_t uly, uint32_t xsize, uint32 /* Disable multi-threaded decompression in Gtiff driver */ CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "1"); - samplingMutex.unlock(); + unlockSampling(); return data; } diff --git a/packages/geo/GeoRaster.h b/packages/geo/GeoRaster.h index 340b2a86d..5210e6e29 100644 --- a/packages/geo/GeoRaster.h +++ b/packages/geo/GeoRaster.h @@ -90,7 +90,6 @@ class GeoRaster: public RasterObject * Data *--------------------------------------------------------------------*/ - Mutex samplingMutex; GdalRaster raster; /*-------------------------------------------------------------------- diff --git a/packages/geo/RasterObject.cpp b/packages/geo/RasterObject.cpp index 335df530a..8eddeaa44 100644 --- a/packages/geo/RasterObject.cpp +++ b/packages/geo/RasterObject.cpp @@ -180,6 +180,122 @@ bool RasterObject::registerRaster (const char* _name, factory_f create) return status; } +/*---------------------------------------------------------------------------- + * getSamples + *----------------------------------------------------------------------------*/ +uint32_t RasterObject::getSamples(const std::vector& points, List& sllist, void* param) +{ + static_cast(param); + uint32_t ssErrors = SS_NO_ERRORS; + + samplingMut.lock(); + try + { + /* Get maximum number of batch processing threads allowed */ + const uint32_t maxNumThreads = getMaxBatchThreads(); + + /* Get readers ranges */ + std::vector ranges; + getThreadsRanges(ranges, points.size(), 5, maxNumThreads); + + for(uint32_t i = 0; i < ranges.size(); i++) + { + const range_t& range = ranges[i]; + mlog(DEBUG, "ragne-%u: %u to %u\n", i, range.start, range.end); + } + + const uint32_t numThreads = ranges.size(); + mlog(INFO, "Number of reader threads: %u", numThreads); + + if(numThreads == 1) + { + /* Single thread, read all samples in one thread using this RasterObject */ + std::vector samples; + ssErrors = readSamples(this, ranges[0], points, samples); + for(sample_list_t* slist : samples) + { + sllist.add(slist); + } + } + else + { + /* Start reader threads */ + std::vector pids; + + for(uint32_t i = 0; i < numThreads; i++) + { + /* Create a RasterObject for each reader thread. + * These objects are local and will be deleted in the reader destructor. + * The user's (this) RasterObject is not directly used for sampling; it is used to accumulate samples from all readers. + */ + RasterObject* _robj = RasterObject::cppCreate(this); + reader_t* reader = new reader_t(_robj, points); + reader->range = ranges[i]; + readersMut.lock(); + { + readers.push_back(reader); + } + readersMut.unlock(); + Thread* pid = new Thread(readerThread, reader); + pids.push_back(pid); + } + + /* Wait for all reader threads to finish */ + for(Thread* pid : pids) + { + delete pid; + } + + /* Copy samples lists (slist pointers only) from each reader. */ + for(const reader_t* reader : readers) + { + /* Acumulate errors from all reader threads */ + ssErrors |= reader->ssErrors; + + for(sample_list_t* slist : reader->samples) + { + for(int32_t i = 0; i < slist->length(); i++) + { + /* NOTE: sample.fileId is an index of the file name in the reader's file dictionary. + * we need to convert it to the index in the batch sampler's dictionary (user's RasterObject dict). + */ + RasterSample* sample = slist->get(i); + + /* Find the file name for the sample id in reader's dictionary */ + const char* name = reader->robj->fileDictGetFile(sample->fileId); + + /* Use user's RasterObject dictionary to store the file names. */ + const uint64_t id = fileDictAdd(name); + + /* Update the sample file id */ + sample->fileId = id; + } + + sllist.add(slist); + } + } + + /* Clear raders */ + readersMut.lock(); + { + for(const reader_t* reader : readers) + delete reader; + + readers.clear(); + } + readersMut.unlock(); + + } + } + catch (const RunTimeException &e) + { + mlog(e.level(), "Error getting samples: %s", e.what()); + } + samplingMut.unlock(); + + return ssErrors; +} + /*---------------------------------------------------------------------------- * getPixels *----------------------------------------------------------------------------*/ @@ -199,9 +315,10 @@ uint8_t* RasterObject::getPixels(uint32_t ulx, uint32_t uly, uint32_t xsize, uin uint32_t RasterObject::getMaxBatchThreads(void) { /* Maximum number of batch threads. - * Each batch thread creates multiple raster reading threads. + * Each batch thread may create multiple raster reading threads. */ - return MAX_BATCH_THREADS; + const uint32_t maxThreads = 16; + return std::min(std::thread::hardware_concurrency(), maxThreads); } /*---------------------------------------------------------------------------- @@ -213,6 +330,17 @@ RasterObject::~RasterObject(void) parms->releaseLuaObject(); } +void RasterObject::stopSampling(void) +{ + sampling = false; + readersMut.lock(); + { + for(const reader_t* reader : readers) + reader->robj->stopSampling(); + } + readersMut.unlock(); +} + /*---------------------------------------------------------------------------- * fileDictAdd *----------------------------------------------------------------------------*/ @@ -243,6 +371,46 @@ const char* RasterObject::fileDictGetFile (uint64_t fileId) return NULL; } +/*---------------------------------------------------------------------------- + * getThreadsRanges + *----------------------------------------------------------------------------*/ +void RasterObject::getThreadsRanges(std::vector& ranges, uint32_t num, + uint32_t minPerThread, uint32_t maxNumThreads) +{ + ranges.clear(); + + /* Determine how many threads to use */ + if(num <= minPerThread) + { + ranges.emplace_back(range_t{0, num}); + return; + } + + uint32_t numThreads = std::min(maxNumThreads, num / minPerThread); + + /* Ensure at least two threads if num > minPerThread */ + if(numThreads == 1 && maxNumThreads > 1) + { + numThreads = 2; + } + + const uint32_t pointsPerThread = num / numThreads; + uint32_t remainingPoints = num % numThreads; + + uint32_t start = 0; + for(uint32_t i = 0; i < numThreads; i++) + { + const uint32_t end = start + pointsPerThread + (remainingPoints > 0 ? 1 : 0); + ranges.emplace_back(range_t{start, end}); + + start = end; + if(remainingPoints > 0) + { + remainingPoints--; + } + } +} + /****************************************************************************** * PROTECTED METHODS ******************************************************************************/ @@ -252,7 +420,8 @@ const char* RasterObject::fileDictGetFile (uint64_t fileId) *----------------------------------------------------------------------------*/ RasterObject::RasterObject(lua_State *L, GeoParms* _parms): LuaObject(L, OBJECT_TYPE, LUA_META_NAME, LUA_META_TABLE), - parms(_parms) + parms(_parms), + sampling(true) { /* Add Lua Functions */ LuaEngine::setAttrFunc(L, "sample", luaSamples); @@ -416,6 +585,27 @@ int RasterObject::luaSubsets(lua_State *L) * PRIVATE METHODS ******************************************************************************/ +/*---------------------------------------------------------------------------- + * Reader Constructor + *----------------------------------------------------------------------------*/ +RasterObject::Reader::Reader(RasterObject* _robj, const std::vector& _points) : + robj(_robj), + range({0, 0}), + points(_points), + ssErrors(SS_NO_ERRORS) +{ +} + +/*---------------------------------------------------------------------------- + * Reader Destructor + *----------------------------------------------------------------------------*/ +RasterObject::Reader::~Reader(void) +{ + delete robj; /* This is locally created RasterObject, not lua created */ +} + + + /*---------------------------------------------------------------------------- * slist2table *----------------------------------------------------------------------------*/ @@ -466,4 +656,65 @@ int RasterObject::slist2table(const List& slist, uint32_t errors, else mlog(DEBUG, "No subsets read"); return num_ret; -} \ No newline at end of file +} + +/*---------------------------------------------------------------------------- + * readerThread + *----------------------------------------------------------------------------*/ +void* RasterObject::readerThread(void* parm) +{ + reader_t* reader = static_cast(parm); + reader->ssErrors = readSamples(reader->robj, reader->range, reader->points, reader->samples); + + /* Exit Thread */ + return NULL; +} + +/*---------------------------------------------------------------------------- + * readSamples + *----------------------------------------------------------------------------*/ +uint32_t RasterObject::readSamples(RasterObject* robj, const range_t& range, + const std::vector& points, + std::vector& samples) +{ + uint32_t ssErrors = SS_NO_ERRORS; + + for(uint32_t i = range.start; i < range.end; i++) + { + if(!robj->sampling) + { + mlog(DEBUG, "Sampling stopped"); + samples.clear(); + break; + } + + const RasterObject::point_info_t& pinfo = points[i]; + const MathLib::point_3d_t& point = pinfo.point; + const int64_t gps = robj->usePOItime() ? pinfo.gps : 0.0; + + sample_list_t* slist = new sample_list_t; + bool listvalid = true; + const uint32_t err = robj->getSamples(point, gps, *slist, NULL); + + /* Acumulate errors from all getSamples calls */ + ssErrors |= err; + + if(err & SS_THREADS_LIMIT_ERROR) + { + listvalid = false; + mlog(CRITICAL, "Too many rasters to sample"); + } + + if(!listvalid) + { + /* Clear the list but don't delete it, empty slist indicates no samples for this point */ + slist->clear(); + } + + /* Add sample list */ + samples.push_back(slist); + } + + return ssErrors; +} + diff --git a/packages/geo/RasterObject.h b/packages/geo/RasterObject.h index b429000a3..073404944 100644 --- a/packages/geo/RasterObject.h +++ b/packages/geo/RasterObject.h @@ -54,7 +54,6 @@ class RasterObject: public LuaObject /*-------------------------------------------------------------------- * Constants *--------------------------------------------------------------------*/ - static const int MAX_BATCH_THREADS = 16; static const char* OBJECT_TYPE; static const char* LUA_META_NAME; static const struct luaL_Reg LUA_META_TABLE[]; @@ -65,10 +64,36 @@ class RasterObject: public LuaObject typedef RasterObject* (*factory_f) (lua_State* L, GeoParms* _parms); - typedef struct { + typedef struct + { factory_f create; } factory_t; + typedef struct + { + MathLib::point_3d_t point; + double gps; + } point_info_t; + + typedef struct + { + uint32_t start; + uint32_t end; + } range_t; + + typedef List sample_list_t; + typedef struct Reader + { + RasterObject* robj; + range_t range; + const std::vector& points; + std::vector samples; + uint32_t ssErrors; + + explicit Reader (RasterObject* _robj, const std::vector& _points); + ~Reader (void); + } reader_t; + /*-------------------------------------------------------------------- * Methods *--------------------------------------------------------------------*/ @@ -80,22 +105,23 @@ class RasterObject: public LuaObject static RasterObject* cppCreate (const RasterObject* obj); static bool registerRaster (const char* _name, factory_f create); virtual uint32_t getSamples (const MathLib::point_3d_t& point, int64_t gps, List& slist, void* param=NULL) = 0; + virtual uint32_t getSamples (const std::vector& points, List& sllist, void* param=NULL); virtual uint32_t getSubsets (const MathLib::extent_t& extent, int64_t gps, List& slist, void* param=NULL) = 0; virtual uint8_t* getPixels (uint32_t ulx, uint32_t uly, uint32_t xsize=0, uint32_t ysize=0, void* param=NULL); virtual uint32_t getMaxBatchThreads(void); ~RasterObject (void) override; - bool hasZonalStats (void) + bool hasZonalStats (void) const { return parms->zonal_stats; } - bool hasFlags (void) + bool hasFlags (void) const { return parms->flags_file; } - bool usePOItime(void) + bool usePOItime(void) const { return parms->use_poi_time; } @@ -105,8 +131,22 @@ class RasterObject: public LuaObject return fileDict; } + void lockSampling(void) + { + samplingMut.lock(); + } + + void unlockSampling(void) + { + samplingMut.unlock(); + } + + void stopSampling (void); + bool isSampling (void) {return sampling;}; uint64_t fileDictAdd (const std::string& fileName); const char* fileDictGetFile (uint64_t fileId); + static void getThreadsRanges(std::vector& ranges, uint32_t num, + uint32_t minPerThread, uint32_t maxNumThreads); protected: @@ -131,7 +171,12 @@ class RasterObject: public LuaObject * Methods *--------------------------------------------------------------------*/ - static int slist2table(const List& slist, uint32_t errors, lua_State* L); + static int slist2table (const List& slist, uint32_t errors, lua_State* L); + + static void* readerThread (void* parm); + static uint32_t readSamples (RasterObject* robj, const range_t& range, + const std::vector& points, std::vector& samples); + /*-------------------------------------------------------------------- * Data @@ -140,6 +185,11 @@ class RasterObject: public LuaObject static Mutex factoryMut; static Dictionary factories; Dictionary fileDict; + + Mutex readersMut; + Mutex samplingMut; + std::atomic sampling; /* Used by batch getSample code */ + std::vector readers; }; #endif /* __raster_object__ */ diff --git a/packages/geo/RasterSample.h b/packages/geo/RasterSample.h index d077c9bf2..c6b09d876 100644 --- a/packages/geo/RasterSample.h +++ b/packages/geo/RasterSample.h @@ -78,6 +78,13 @@ class RasterSample { } + RasterSample(const RasterSample& sample): RasterSample(sample.time, sample.fileId, sample.verticalShift) + { + value = sample.value; + flags = sample.flags; + stats = sample.stats; + } + std::string toString(void) const { char buffer[1024]; diff --git a/packages/geo/UT_RasterSample.cpp b/packages/geo/UT_RasterSample.cpp new file mode 100644 index 000000000..ae334feb7 --- /dev/null +++ b/packages/geo/UT_RasterSample.cpp @@ -0,0 +1,274 @@ +/* + * Copyright (c) 2021, University of Washington + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the University of Washington nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF WASHINGTON AND CONTRIBUTORS + * “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF WASHINGTON OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/****************************************************************************** + * INCLUDES + ******************************************************************************/ + +#include "core.h" +#include "GdalRaster.h" +#include "UT_RasterSample.h" + +/****************************************************************************** + * STATIC DATA + ******************************************************************************/ + +const char* UT_RasterSample::OBJECT_TYPE = "UT_RasterSample"; + +const char* UT_RasterSample::LUA_META_NAME = "UT_RasterSample"; +const struct luaL_Reg UT_RasterSample::LUA_META_TABLE[] = { + {"test", luaSampleTest}, + {NULL, NULL} +}; + +/****************************************************************************** + * PUBLIC METHODS + ******************************************************************************/ + +/*---------------------------------------------------------------------------- + * luaCreate - :test() + *----------------------------------------------------------------------------*/ +int UT_RasterSample::luaCreate (lua_State* L) +{ + RasterObject* _raster = NULL; + try + { + /* Get Parameters */ + _raster = dynamic_cast(getLuaObject(L, 1, RasterObject::OBJECT_TYPE)); + + return createLuaObject(L, new UT_RasterSample(L, _raster)); + } + catch(const RunTimeException& e) + { + if(_raster) _raster->releaseLuaObject(); + mlog(e.level(), "Error creating %s: %s", LUA_META_NAME, e.what()); + return returnLuaStatus(L, false); + } +} + +/****************************************************************************** + * PRIVATE METHODS + *******************************************************************************/ + +/*---------------------------------------------------------------------------- + * Constructor + *----------------------------------------------------------------------------*/ +UT_RasterSample::UT_RasterSample (lua_State* L, RasterObject* _raster): + LuaObject(L, OBJECT_TYPE, LUA_META_NAME, LUA_META_TABLE), + raster(_raster) +{ + assert(raster); +} + +/*---------------------------------------------------------------------------- + * Destructor - + *----------------------------------------------------------------------------*/ +UT_RasterSample::~UT_RasterSample(void) +{ + raster->releaseLuaObject(); +} + + +/*---------------------------------------------------------------------------- + * luaSampleTest + *----------------------------------------------------------------------------*/ +int UT_RasterSample::luaSampleTest(lua_State* L) +{ + bool status = false; + uint32_t errors = 0; + + UT_RasterSample* lua_obj = NULL; + + try + { + lua_obj = dynamic_cast(getLuaSelf(L, 1)); + + /* Get Coordinates */ + double lon = getLuaFloat(L, 2); + double lat = getLuaFloat(L, 3); + + /* Get Increments */ + const double lonIncr = getLuaFloat(L, 4); + const double latIncr = getLuaFloat(L, 5); + + /* Get Count */ + const long pointsCnt = getLuaInteger(L, 6); + + /* Create of points to sample */ + std::vector points2sample; + for(long i = 0; i < pointsCnt; i++) + { + RasterObject::point_info_t point; + point.point.x = lon; + point.point.y = lat; + point.point.z = 0.0; + point.gps = 0.0; + + points2sample.push_back(point); + + lon += lonIncr; + lat += latIncr; + } + + print2term("Points to sample: %ld\n", pointsCnt); + print2term("Starting at (%.4lf, %.4lf), incrementing by (%+.4lf, %+.4lf)\n", lon, lat, lonIncr, latIncr); + + /* Get samples using serial method */ + print2term("Getting samples using serial method\n"); + List serial_sllist; + const double serialStartTime = TimeLib::latchtime(); + for(uint32_t i = 0; i < points2sample.size(); i++) + { + RasterObject::sample_list_t* slist = new RasterObject::sample_list_t(); + lua_obj->raster->getSamples(points2sample[i].point, 0, *slist, NULL); + + /* Add to list */ + serial_sllist.add(slist); + } + const double serialStopTime = TimeLib::latchtime(); + + /* Get samples using batch method */ + print2term("Getting samples using batch method\n"); + List batch_sllist; + const double batchStartTime = TimeLib::latchtime(); + lua_obj->raster->getSamples(points2sample, batch_sllist, NULL); + const double batchStopTime = TimeLib::latchtime(); + + /* Print times */ + const double serialTime = serialStopTime - serialStartTime; + const double batchTime = batchStopTime - batchStartTime; + print2term("Serial time: %lf\n", serialTime); + print2term("Batch time: %lf\n", batchTime); + + /* Compare times */ + if (batchTime > 0) + { + const double performanceDifference = ((serialTime - batchTime) / batchTime) * 100.0; + print2term("Performance difference: %.2lf%%\n", performanceDifference); + } + + /* Compare results */ + uint32_t validStreamSamples = 0; + uint32_t validBatchSamples = 0; + uint32_t nanStreamSamples = 0; + uint32_t nanBatchSamples = 0; + + print2term("Comparing samples\n"); + if(serial_sllist.length() != batch_sllist.length()) + { + print2term("Number of samples differ, serial: %d, batch: %d\n", serial_sllist.length(), batch_sllist.length()); + errors++; + } + + for(int i = 0; i < serial_sllist.length(); i++) + { + RasterObject::sample_list_t* serial_slist = serial_sllist[i]; + RasterObject::sample_list_t* batch_slist = batch_sllist[i]; + + if(serial_slist->length() != batch_slist->length()) + { + print2term("Number of samples differ, serial: %d, batch: %d\n", serial_slist->length(), batch_slist->length()); + errors++; + } + + for(int j = 0; j < serial_slist->length(); j++) + { + RasterSample* serial = (*serial_slist)[j]; + RasterSample* batch = (*batch_slist)[j]; + + const char* serialName = lua_obj->raster->fileDictGetFile(serial->fileId); + const char* batchName = lua_obj->raster->fileDictGetFile(batch->fileId); + + if (strcmp(serialName, batchName) != 0) + { + print2term("Files differ:\n"); + print2term("Serial: %s\n", serialName); + print2term("Batch: %s\n", batchName); + errors++; + } + + /* Compare as whole seconds */ + if(static_cast(serial->time) != static_cast(batch->time)) + { + print2term("Time differ: %lf != %lf\n", serial->time, batch->time); + errors++; + } + + /* Nodata sample values are returned as NaN, cannot compare two NaNs */ + const bool dontCompare = std::isnan(serial->value) && std::isnan(batch->value); + if(!dontCompare && serial->value != batch->value) + { + print2term("Value differ: %lf != %lf\n", serial->value, batch->value); + errors++; + } + + if(std::isnan(serial->value)) nanStreamSamples++; + else validStreamSamples++; + + if(std::isnan(batch->value)) nanBatchSamples++; + else validBatchSamples++; + + if(serial->stats.mean != batch->stats.mean) + { + print2term("Mean differ: %lf != %lf\n", serial->stats.mean, batch->stats.mean); + errors++; + } + + if(serial->stats.stdev != batch->stats.stdev) + { + print2term("Stdev differ: %lf != %lf\n", serial->stats.stdev, batch->stats.stdev); + errors++; + } + + if(serial->stats.mad != batch->stats.mad) + { + print2term("Mad differ: %lf != %lf\n", serial->stats.mad, batch->stats.mad); + errors++; + } + } + } + + print2term("Stream samples, valid: %u, nan: %u\n", validStreamSamples, nanStreamSamples); + print2term("Batch samples, valid: %u, nan: %u\n", validBatchSamples, nanBatchSamples); + } + catch(const RunTimeException& e) + { + mlog(e.level(), "Error creating %s: %s", LUA_META_NAME, e.what()); + errors++; + } + + if(errors == 0) status = true; + else status = false; + + /* Return Status */ + return returnLuaStatus(L, status); +} + diff --git a/packages/geo/UT_RasterSample.h b/packages/geo/UT_RasterSample.h new file mode 100644 index 000000000..76561a038 --- /dev/null +++ b/packages/geo/UT_RasterSample.h @@ -0,0 +1,89 @@ +/* + * Copyright (c) 2021, University of Washington + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the University of Washington nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF WASHINGTON AND CONTRIBUTORS + * “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY OF WASHINGTON OR + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef __ut_raster_sample__ +#define __ut_raster_sample__ + +/****************************************************************************** + * INCLUDES + ******************************************************************************/ + +#include "OsApi.h" +#include "LuaObject.h" +#include "RasterObject.h" + +/****************************************************************************** + * ATL03 READER UNIT TEST CLASS + ******************************************************************************/ + +class UT_RasterSample: public LuaObject +{ + public: + + /*-------------------------------------------------------------------- + * Types + *--------------------------------------------------------------------*/ + + /*-------------------------------------------------------------------- + * Constants + *--------------------------------------------------------------------*/ + + static const char* OBJECT_TYPE; + + static const char* LUA_META_NAME; + static const struct luaL_Reg LUA_META_TABLE[]; + + /*-------------------------------------------------------------------- + * Methods + *--------------------------------------------------------------------*/ + + static int luaCreate (lua_State* L); + + private: + + + /*-------------------------------------------------------------------- + * Methods + *--------------------------------------------------------------------*/ + + explicit UT_RasterSample (lua_State* L, RasterObject* _raster); + ~UT_RasterSample (void) override; + + static int luaSampleTest (lua_State* L); + + /*-------------------------------------------------------------------- + * Data + *--------------------------------------------------------------------*/ + + RasterObject* raster; +}; + +#endif /* __ut_raster_subset__*/ diff --git a/packages/geo/geo.cpp b/packages/geo/geo.cpp index bf93aa3df..c45954cdc 100644 --- a/packages/geo/geo.cpp +++ b/packages/geo/geo.cpp @@ -37,6 +37,7 @@ #ifdef __unittesting__ #include "UT_RasterSubset.h" +#include "UT_RasterSample.h" #endif #include "GeoRaster.h" @@ -254,7 +255,8 @@ int geo_open (lua_State* L) {"calcutm", GeoLib::luaCalcUTM}, {"tiff", GeoLib::TIFFImage::luaCreate}, #ifdef __unittesting__ - {"ut_subset", UT_RasterSubset::luaCreate}, + {"ut_subset", UT_RasterSubset::luaCreate}, + {"ut_sample", UT_RasterSample::luaCreate}, #endif {NULL, NULL} }; diff --git a/scripts/selftests/test_runner.lua b/scripts/selftests/test_runner.lua index 1f7e00789..4e85c42db 100644 --- a/scripts/selftests/test_runner.lua +++ b/scripts/selftests/test_runner.lua @@ -102,6 +102,7 @@ end -- Run Landsat Plugin Self Tests if __landsat__ and incloud then local landsat_td = td .. "../../datasets/landsat/selftests/" + runner.script(landsat_td .. "plugin_unittest.lua") runner.script(landsat_td .. "landsat_reader.lua") end diff --git a/scripts/systests/parquet_sampler_perf_test.lua b/scripts/systests/parquet_sampler_perf_test.lua index 072a8541c..f284d2ffc 100644 --- a/scripts/systests/parquet_sampler_perf_test.lua +++ b/scripts/systests/parquet_sampler_perf_test.lua @@ -6,7 +6,9 @@ local td = runner.rootdir(arg[0]) local outq_name = "outq-luatest" -local in_parquet = '/data/3dep/wrzesien_snow.parquet' +local in_parquet = '/data/3dep/wrzesien_snow_64k.parquet' +-- local in_parquet = '/data/3dep/wrzesien_snow_525k.parquet' +-- local in_parquet = '/data/3dep/wrzesien_snow_2618k.parquet' -- Indicates local file system (no s3 or client) local prefix = "file://" @@ -34,8 +36,7 @@ function getFileSize(filePath) return fileSize end -local robj_3dep = geo.raster(geo.parms({asset="usgs3dep-1meter-dem", algorithm="NearestNeighbour", radius=0, catalog = contents, use_poi_time=true})) --- local robj_3dep = geo.raster(geo.parms({asset="usgs3dep-1meter-dem", algorithm="NearestNeighbour", radius=0, catalog = contents, use_poi_time=true, single_stop=true, })) +local robj_3dep = geo.raster(geo.parms({asset="usgs3dep-1meter-dem", algorithm="NearestNeighbour", radius=0, catalog = contents, use_poi_time=false})) runner.check(robj_3dep ~= nil) print('\n--------------------------------------\nTest01: input/output parquet (x, y)\n--------------------------------------') diff --git a/targets/slideruleearth-aws/Makefile b/targets/slideruleearth-aws/Makefile index c64332bf5..730bd1bb0 100644 --- a/targets/slideruleearth-aws/Makefile +++ b/targets/slideruleearth-aws/Makefile @@ -313,8 +313,8 @@ distclean: ## fully remove all non-version controlled files and directories ######################## TEST ?= $(ROOT)/scripts/selftests/test_runner.lua cloud -MOSAICS_PERFORMANCE_TEST ?= $(ROOT)/plugins/pgc/systests/arcticdem_mosaic_perf.lua -STRIPS_PERFORMANCE_TEST ?= $(ROOT)/plugins/pgc/systests/arcticdem_strips_perf.lua +MOSAICS_PERFORMANCE_TEST ?= $(ROOT)/datasets/pgc/systests/arcticdem_mosaic_perf.lua +STRIPS_PERFORMANCE_TEST ?= $(ROOT)/datasets/pgc/systests/arcticdem_strips_perf.lua SUBSET_PERFORMANCE_TEST ?= $(ROOT)/scripts/systests/subset_perf_test.lua selftest: ## run the self test on the server code