optimized geojson includes

packages/geo/GdalRaster.cpp

@@ -339,7 +339,7 @@ RasterSubset* GdalRaster::subsetAOI(OGRPolygon* poly)
 
         char* wkt;
         targetCRS.exportToWkt(&wkt);
-        subset = new RasterSubset(cols2read, rows2read, datatype, aoi_minx, aoi_maxy, cellSize, wkt, gpsTime, fileId);
+        subset = new RasterSubset(cols2read, rows2read, datatype, aoi_minx, aoi_maxy, cellSize, bbox, wkt, gpsTime, fileId);
         CPLFree(wkt);
 
         if(subset->data)
@@ -387,6 +387,117 @@ RasterSubset* GdalRaster::subsetAOI(OGRPolygon* poly)
     return subset;
 }
 
+
+/*----------------------------------------------------------------------------
+ * subsetAOI
+ *----------------------------------------------------------------------------*/
+RasterSubset* GdalRaster::subsetAOI(uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize)
+{
+    RasterSubset* subset = NULL;
+
+    /* Clear sample/subset error status */
+    ssError = SS_NO_ERRORS;
+
+    try
+    {
+        if(dset == NULL)
+            open();
+
+        if(ulx >= xsize)
+            throw RunTimeException(CRITICAL, RTE_ERROR, "Upleft pixel's x out of bounds: %u", ulx);
+
+        if(uly >= ysize)
+            throw RunTimeException(CRITICAL, RTE_ERROR, "Upleft pixel's y out of bounds: %u", uly);
+
+        if(_xsize == 0)
+        {
+            /* Read all raster columns starting at ulx */
+            _xsize = xsize - ulx;
+        }
+
+        if(_ysize == 0)
+        {
+            /* Read all raster rows starting at uly */
+            _ysize = ysize - uly;
+        }
+
+        if(ulx + _xsize > xsize)
+            throw RunTimeException(CRITICAL, RTE_ERROR, "columns out of bounds");
+
+        if(uly + _ysize > ysize)
+            throw RunTimeException(CRITICAL, RTE_ERROR, "rows out of bounds");
+
+
+        int64_t cols2read = _xsize;
+        int64_t rows2read = _ysize;
+
+        GDALDataType dtype = band->GetRasterDataType();
+        RecordObject::fieldType_t datatype = RecordObject::INVALID_FIELD;
+        switch (dtype) {
+            case GDT_Byte:      datatype = RecordObject::UINT8;     break;
+            case GDT_Int8:      datatype = RecordObject::UINT8;     break;
+            case GDT_UInt16:    datatype = RecordObject::UINT16;    break;
+            case GDT_Int16:     datatype = RecordObject::INT16;     break;
+            case GDT_UInt32:    datatype = RecordObject::UINT32;    break;
+            case GDT_Int32:     datatype = RecordObject::INT32;     break;
+            case GDT_UInt64:    datatype = RecordObject::UINT64;    break;
+            case GDT_Int64:     datatype = RecordObject::INT64;     break;
+            case GDT_Float32:   datatype = RecordObject::FLOAT;     break;
+            case GDT_Float64:   datatype = RecordObject::DOUBLE;    break;
+            default:            throw RunTimeException(CRITICAL, RTE_ERROR, "Unsupported GDT Datatype: %d", dtype);
+        }
+
+        char* wkt;
+        targetCRS.exportToWkt(&wkt);
+
+        double mapx, mapy;
+        pixel2map(ulx, uly, mapx, mapy);
+
+        subset = new RasterSubset(cols2read, rows2read, datatype, mapx, mapy, cellSize, bbox, wkt, gpsTime, fileId);
+        CPLFree(wkt);
+
+        if(subset->data)
+        {
+            int cnt = 1;
+            OGRErr err = CE_None;
+            do
+            {
+                err = band->RasterIO(GF_Read, ulx, uly, subset->cols, subset->rows, subset->data, subset->cols, subset->rows, dtype, 0, 0, NULL);
+            }
+            while(err != CE_None && cnt--);
+
+            if(err != CE_None)
+            {
+                ssError |= SS_AOI_FAILED_TO_READ_ERROR;
+                throw RunTimeException(CRITICAL, RTE_ERROR, "RasterIO call failed: %d", err);
+            }
+
+            mlog(DEBUG, "read %ld bytes (%.1fMB), pixel_ulx: %d, pixel_uly: %d, map_ulx: %.2lf, map_uly: %.2lf, cols2read: %ld, rows2read: %ld, datatype %s\n",
+                 subset->size, (float)subset->size/(1024*1024), ulx, uly, subset->map_ulx, subset->map_uly, subset->cols, subset->rows, GDALGetDataTypeName(dtype));
+        }
+        else
+        {
+            ssError |= SS_MEMPOOL_ERROR;
+            uint64_t size = cols2read * rows2read * GDALGetDataTypeSizeBytes(dtype);
+            mlog(ERROR, "RasterSubset requested memory: %lu MB, available: %lu MB, max: %lu MB", size / (1024*1024),
+                RasterSubset::poolsize / (1024*1024),
+                RasterSubset::MAX_SIZE / (1024*1024));
+        }
+    }
+    catch (const RunTimeException &e)
+    {
+        if(subset)
+        {
+            delete subset;
+            subset = NULL;
+        }
+        mlog(e.level(), "Error subsetting: %s", e.what());
+    }
+
+    return subset;
+}
+
+
 /*----------------------------------------------------------------------------
  * setCRSfromWkt
  *----------------------------------------------------------------------------*/

packages/geo/GdalRaster.h

@@ -106,6 +106,7 @@ class GdalRaster
         void               open           (void);
         RasterSample*      samplePOI      (OGRPoint* poi);
         RasterSubset*      subsetAOI      (OGRPolygon* poly);
+        RasterSubset*      subsetAOI      (uint32_t ulx, uint32_t uly, uint32_t _xsize, uint32_t _ysize);
         const std::string& getFileName    (void) const { return fileName;}
         int                getRows        (void) const { return ysize; }
         int                getCols        (void) const { return xsize; }

packages/geo/GeoJsonRaster.cpp

@@ -81,7 +81,7 @@ GeoJsonRaster* GeoJsonRaster::create (lua_State* L, int index)
 
     /* Get cellsize */
     lua_getfield(L, index, CELLSIZE_KEY);
-    double _cellsize = getLuaFloat(L, -1);
+    double cellsize = getLuaFloat(L, -1);
     lua_pop(L, 1);
 
     /* Get Geo Parameters */
@@ -91,7 +91,7 @@ GeoJsonRaster* GeoJsonRaster::create (lua_State* L, int index)
     lua_pop(L, 1);
 
     /* Create GeoJsonRaster */
-    return new GeoJsonRaster(L, _parms, geojstr, _cellsize);
+    return new GeoJsonRaster(L, _parms, geojstr, cellsize);
 }
 
 
@@ -100,18 +100,29 @@ GeoJsonRaster* GeoJsonRaster::create (lua_State* L, int index)
  *----------------------------------------------------------------------------*/
 bool GeoJsonRaster::includes(double lon, double lat, double height)
 {
-    std::vector<RasterSample*> slist;
+    std::ignore = height;
+    bool pixel_on = false;
 
-    OGRPoint poi(lon, lat, height);
-    getSamples(&poi, 0, slist);
-    int sampleCnt = slist.size();
+    /*
+     * Skip transforming POI since geojsons should be in geographic coordinates.
+     * Raster created from geojson is also in geo.
+     *
+     * Don't need a mutex, multiple threads can read the same data.
+     */
 
-    if( sampleCnt == 0 ) return false; // no need to delete anything
-    if( sampleCnt > 1  ) mlog(ERROR, "Multiple samples returned for lon: %.2lf, lat: %.2lf, using first sample", lon, lat);
+    const GdalRaster::bbox_t& bbox = subset->bbox;
 
-    bool pixel_on = static_cast<int>(slist[0]->value) == RASTER_PIXEL_ON;
+    if((lon >= bbox.lon_min) && (lon <= bbox.lon_max) &&
+       (lat >= bbox.lat_min) && (lat <= bbox.lat_max))
+    {
+        uint32_t row = (bbox.lat_max - lat) / subset->cellsize;
+        uint32_t col = (lon - bbox.lon_min) / subset->cellsize;
 
-    for(auto sample: slist) delete sample;
+        if((row < subset->rows) && (col < subset->cols))
+        {
+            pixel_on = rawPixel(row, col);
+        }
+    }
 
     return pixel_on;
 }
@@ -121,6 +132,7 @@ bool GeoJsonRaster::includes(double lon, double lat, double height)
  *----------------------------------------------------------------------------*/
 GeoJsonRaster::~GeoJsonRaster(void)
 {
+    delete subset;
     VSIUnlink(rasterFileName.c_str());
 }
 
@@ -132,8 +144,9 @@ GeoJsonRaster::~GeoJsonRaster(void)
 /*----------------------------------------------------------------------------
  * Constructor
  *----------------------------------------------------------------------------*/
-GeoJsonRaster::GeoJsonRaster(lua_State* L, GeoParms* _parms, const char* geojstr, double _cellsize):
- GeoRaster(L, _parms, std::string("/vsimem/" + GdalRaster::getUUID() + ".tif"), TimeLib::gpstime(), false /* not elevation*/)
+GeoJsonRaster::GeoJsonRaster(lua_State* L, GeoParms* _parms, const char* geojstr, double cellsize):
+ GeoRaster(L, _parms, std::string("/vsimem/" + GdalRaster::getUUID() + ".tif"), TimeLib::gpstime(), false /* not elevation*/),
+ subset(NULL)
 {
     bool rasterCreated = false;
     GDALDataset* rasterDset = NULL;
@@ -144,8 +157,8 @@ GeoJsonRaster::GeoJsonRaster(lua_State* L, GeoParms* _parms, const char* geojstr
     if (geojstr == NULL)
         throw RunTimeException(CRITICAL, RTE_ERROR, "Invalid file pointer (NULL)");
 
-    if (_cellsize <= 0.0)
-        throw RunTimeException(CRITICAL, RTE_ERROR, "Invalid cellSize: %.2lf:", _cellsize);
+    if (cellsize <= 0.0)
+        throw RunTimeException(CRITICAL, RTE_ERROR, "Invalid cellSize: %.2lf:", cellsize);
 
     try
     {
@@ -164,9 +177,8 @@ GeoJsonRaster::GeoJsonRaster(lua_State* L, GeoParms* _parms, const char* geojstr
         OGRErr ogrerr = srcLayer->GetExtent(&e);
         CHECK_GDALERR(ogrerr);
 
-        double cellsize = _cellsize;
-        int cols = int((e.MaxX - e.MinX) / cellsize);
-        int rows = int((e.MaxY - e.MinY) / cellsize);
+        int cols = static_cast<int>((e.MaxX - e.MinX) / cellsize);
+        int rows = static_cast<int>((e.MaxY - e.MinY) / cellsize);
 
         char **options = NULL;
         options = CSLSetNameValue(options, "COMPRESS", "DEFLATE");
@@ -216,6 +228,11 @@ GeoJsonRaster::GeoJsonRaster(lua_State* L, GeoParms* _parms, const char* geojstr
         GDALClose((GDALDatasetH)rasterDset);
         rasterDset = NULL;
         rasterCreated = true;
+
+        /* Subset newly created raster, get all pixels */
+        subset = getSubset();
+        CHECKPTR(subset);
+        CHECKPTR(subset->data);
     }
     catch(const RunTimeException& e)
     {

packages/geo/GeoJsonRaster.h

@@ -68,21 +68,31 @@ class GeoJsonRaster: public GeoRaster
         bool                  includes       (double lon, double lat, double height=0);
         virtual              ~GeoJsonRaster  (void);
 
+        /*--------------------------------------------------------------------
+         * Inline Methods
+         *--------------------------------------------------------------------*/
+
+        bool rawPixel (const uint32_t row, const uint32_t col)
+        {
+            return static_cast<int>(subset->data[(row * subset->cols) + col]) == RASTER_PIXEL_ON;
+        }
+
     protected:
 
         /*--------------------------------------------------------------------
          * Methods
          *--------------------------------------------------------------------*/
 
-        GeoJsonRaster(lua_State* L, GeoParms* _parms, const char* geojstr, double _cellsize);
+        GeoJsonRaster(lua_State* L, GeoParms* _parms, const char* geojstr, double cellsize);
 
     private:
 
         /*--------------------------------------------------------------------
          * Data
          *--------------------------------------------------------------------*/
 
-         std::string rasterFileName;
+         std::string   rasterFileName;
+         RasterSubset* subset;
 };
 
 #endif  /* __geojson_raster__ */

packages/geo/GeoRaster.cpp

@@ -63,7 +63,7 @@ uint32_t GeoRaster::getSamples(OGRGeometry* geo, int64_t gps, std::vector<Raster
 }
 
 /*----------------------------------------------------------------------------
- * getSubset
+ * getSubsets
  *----------------------------------------------------------------------------*/
 uint32_t GeoRaster::getSubsets(OGRGeometry* geo, int64_t gps, std::vector<RasterSubset*>& slist, void* param)
 {
@@ -95,6 +95,38 @@ uint32_t GeoRaster::getSubsets(OGRGeometry* geo, int64_t gps, std::vector<Raster
 }
 
 
+/*----------------------------------------------------------------------------
+ * getSubset
+ *----------------------------------------------------------------------------*/
+RasterSubset* GeoRaster::getSubset(uint32_t ulx, uint32_t uly, uint32_t xsize, uint32_t ysize, void* param)
+{
+    std::ignore = param;
+    RasterSubset* subset = NULL;
+
+    samplingMutex.lock();
+
+    /* Enable multi-threaded decompression in Gtiff driver */
+    CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
+
+    try
+    {
+        subset = raster.subsetAOI(ulx, uly, xsize, ysize);
+    }
+    catch (const RunTimeException &e)
+    {
+        mlog(e.level(), "Error subsetting raster: %s", e.what());
+        delete subset;
+        subset = NULL;
+    }
+
+    /* Disable multi-threaded decompression in Gtiff driver */
+    CPLSetThreadLocalConfigOption("GDAL_NUM_THREADS", "1");
+
+    samplingMutex.unlock();
+
+    return subset;
+}
+
 /*----------------------------------------------------------------------------
  * Destructor
  *----------------------------------------------------------------------------*/

packages/geo/GeoRaster.h

@@ -51,9 +51,10 @@ class GeoRaster: public RasterObject
          * Methods
          *--------------------------------------------------------------------*/
 
-        virtual ~GeoRaster  (void);
-        uint32_t getSamples (OGRGeometry* geo, int64_t gps, std::vector<RasterSample*>& slist, void* param=NULL) final;
-        uint32_t getSubsets (OGRGeometry* geo, int64_t gps, std::vector<RasterSubset*>& slist, void* param=NULL) final;
+        virtual      ~GeoRaster  (void);
+        uint32_t      getSamples (OGRGeometry* geo, int64_t gps, std::vector<RasterSample*>& slist, void* param=NULL) final;
+        uint32_t      getSubsets (OGRGeometry* geo, int64_t gps, std::vector<RasterSubset*>& slist, void* param=NULL) final;
+        RasterSubset* getSubset  (uint32_t ulx=0, uint32_t uly=0, uint32_t xsize=0, uint32_t ysize=0, void* param=NULL);
 
     protected:
 

packages/geo/RasterObject.cpp

@@ -368,6 +368,10 @@ int RasterObject::luaSubset(lua_State *L)
                 LuaEngine::setAttrNum(L, "ulx", subset->map_ulx);
                 LuaEngine::setAttrNum(L, "uly", subset->map_uly);
                 LuaEngine::setAttrNum(L, "cellsize", subset->cellsize);
+                LuaEngine::setAttrNum(L, "bbox.lonmin", subset->bbox.lon_min);
+                LuaEngine::setAttrNum(L, "bbox.latmin", subset->bbox.lat_min);
+                LuaEngine::setAttrNum(L, "bbox.lonmax", subset->bbox.lon_max);
+                LuaEngine::setAttrNum(L, "bbox.latmax", subset->bbox.lat_max);
                 LuaEngine::setAttrStr(L, "wkt", subset->wkt.c_str());
                 lua_rawseti(L, -2, i+1);
             }

packages/geo/RasterSubset.cpp

@@ -52,7 +52,7 @@ Mutex RasterSubset::mutex;
  * Constructor
  *----------------------------------------------------------------------------*/
 RasterSubset::RasterSubset(uint32_t _cols, uint32_t _rows, RecordObject::fieldType_t _datatype,
-                           double ulx, double uly, double _cellsize, const char* _wkt,
+                           double ulx, double uly, double _cellsize, GeoParms::bbox_t& _bbox, const char* _wkt,
                            double _time, double _fileId):
     data(NULL),
     size(0),
@@ -62,6 +62,7 @@ RasterSubset::RasterSubset(uint32_t _cols, uint32_t _rows, RecordObject::fieldTy
     map_ulx(ulx),
     map_uly(uly),
     cellsize(_cellsize),
+    bbox(_bbox),
     wkt(_wkt),
     time(_time),
     fileId(_fileId)

packages/geo/RasterSubset.h

@@ -38,6 +38,7 @@
 
 #include "OsApi.h"
 #include "RecordObject.h"
+#include "GeoParms.h"
 
 /******************************************************************************
  * RASTER SUBSET CLASS
@@ -55,6 +56,7 @@ class RasterSubset
         double                      map_ulx;
         double                      map_uly;
         double                      cellsize;
+        GeoParms::bbox_t            bbox;
         std::string                 wkt;
         double                      time;     // gps seconds
         uint64_t                    fileId;
@@ -66,7 +68,7 @@ class RasterSubset
         static Mutex                mutex;
 
         RasterSubset(uint32_t _cols, uint32_t _rows, RecordObject::fieldType_t _datatype,
-                     double ulx, double uly, double _cellsize, const char* _wkt,
+                     double ulx, double uly, double _cellsize, GeoParms::bbox_t& _bbox, const char* _wkt,
                      double _time, double _fileId);
         ~RasterSubset(void);
 };