diff --git a/.zenodo.json b/.zenodo.json index dc6c0adeaf0661332a03fe95bebd30ece5021e5a..b540a52fc46970be7cf18602aed572f940e3c33b 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -166,6 +166,11 @@ "affiliation": "United States Geological Survey, Astro Geology Science Center", "name": "Fergason, Robin" }, + { + "affiliation": "Italian National Institute for Astrophysics (INAF), Istituto di Astrofisica e Planetologia Spaziali (IAPS), Rome, Italy", + "name": "Frigeri, Alessandro", + "orcid": "0000-0002-9140-3977" + }, { "name": "Gaddis, Lisa" }, diff --git a/AUTHORS.rst b/AUTHORS.rst index d564245b80844f30b30ba1a4dd4690e9c7cabf10..e3330cd117ef81bac5be392b6638d352bd0bfb68 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -176,4 +176,4 @@ Integrated Software for Imagers and Spectrometers Contributors ----- This list was generated from the .zenodo.json file by running the -isis/scripts/zenodo_to_authors.py Python program. \ No newline at end of file +isis/scripts/zenodo_to_authors.py Python program. diff --git a/CHANGELOG.md b/CHANGELOG.md index efde84a228212a4959f163d4e642d7dfba1eec5b..e6861eaa1cdebe551336e36cc7a9c39d2c6ac205 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -35,8 +35,16 @@ release. ## [Unreleased] -### Fixed +### Added +- Added TOVECT output parameter which generate a geospatial CSV file with a VRT metadata sidecar file [#5571](https://github.com/DOI-USGS/ISIS3/issues/5571) +- Added Vectorize to ProcessGroundPolygon library +- Added gtest files for the app and unit test +### Changed +- Refactored the pixel2map app +- Updated pixel2map documentation + +### Fixed - Fixed a bug in kaguyasp2isis that doesn't work for data with a detached label. ## [8.3.0] - 2024-08-16 diff --git a/isis/src/base/apps/pixel2map/main.cpp b/isis/src/base/apps/pixel2map/main.cpp index 825a383f159d033d322a157ca3f408b1e644ecd8..ef3e5a1765336334107453c6dde22ed7c65c815f 100755 --- a/isis/src/base/apps/pixel2map/main.cpp +++ b/isis/src/base/apps/pixel2map/main.cpp @@ -1,364 +1,34 @@ #define GUIHELPERS - #include "Isis.h" -#include <QDebug> -#include <QList> -#include <QPointF> -#include <QString> - -#include "Camera.h" -#include "Cube.h" -#include "CubeAttribute.h" -#include "FileList.h" -#include "IException.h" -#include "PixelFOV.h" -#include "ProcessByBrick.h" -#include "ProcessGroundPolygons.h" -#include "ProcessRubberSheet.h" -#include "ProjectionFactory.h" -#include "Pvl.h" -#include "PvlGroup.h" -#include "Target.h" +#include "Application.h" #include "pixel2map.h" +#include "UserInterface.h" + using namespace std; using namespace Isis; void PrintMap(); -void rasterizePixel(Isis::Buffer &in); -map <QString, void *> GuiHelpers() { +// Check for docs for GUI helpers -> KEEP and move it to pixel2map.cpp + std::map <QString, void *> GuiHelpers() { map <QString, void *> helper; helper ["PrintMap"] = (void *) PrintMap; return helper; } -// Global variables -ProcessGroundPolygons g_processGroundPolygons; -Camera *g_incam; -int g_numIFOVs = 0; void IsisMain() { - - g_incam = NULL; - - // Get the map projection file provided by the user UserInterface &ui = Application::GetUserInterface(); - Pvl userMap; - userMap.read(ui.GetFileName("MAP")); - PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse); - - FileList list; - if (ui.GetString("FROMTYPE") == "FROM") { - // GetAsString will capture the entire string, including attributes - list.push_back(FileName(ui.GetAsString("FROM"))); - } - else { - try { - list.read(ui.GetFileName("FROMLIST")); - } - catch (IException &e) { - throw IException(e); - } - } - - if (ui.GetString("FOVRANGE") == "INSTANTANEOUS") { - g_numIFOVs = 1; - } - else { - g_numIFOVs = ui.GetInteger("NUMIFOV"); - } - - double newminlat, newmaxlat, newminlon, newmaxlon; - double minlat = 90; - double maxlat = -90; - double minlon = 360.0; - double maxlon = 0; - PvlGroup camGrp; - QString lastBandString; - - // Get the combined lat/lon range for all input cubes - int bands = 1; - for (int i = 0; i < list.size(); i++) { - - // Open the input cube and get the camera - CubeAttributeInput atts0(list[i]); - Cube icube; - if(atts0.bands().size() != 0) { - vector<QString> lame = atts0.bands(); - icube.setVirtualBands(lame); - } - icube.open( list[i].toString() ); - bands = icube.bandCount(); - g_incam = icube.camera(); - - // Make sure it is not the sky - if (g_incam->target()->isSky()) { - QString msg = "The image [" + list[i].toString() + - "] is targeting the sky, use skymap instead."; - throw IException(IException::User, msg, _FILEINFO_); - } - - // Make sure all the bands for all the files match - if (i >= 1 && atts0.bandsString() != lastBandString) { - QString msg = "The Band numbers for all the files do not match."; - throw IException(IException::User, msg, _FILEINFO_); - } - else { - lastBandString = atts0.bandsString(); - } - - // Get the mapping group - Pvl camMap; - g_incam->BasicMapping(camMap); - camGrp = camMap.findGroup("Mapping"); - - g_incam->GroundRange(newminlat, newmaxlat, newminlon, newmaxlon, userMap); - //set min lat/lon - if (newminlat < minlat) { - minlat = newminlat; - } - if (newminlon < minlon) { - minlon = newminlon; - } - - //set max lat/lon - if (newmaxlat > maxlat) { - maxlat = newmaxlat; - } - if (newmaxlon > maxlon) { - maxlon = newmaxlon; - } - - // The camera will be deleted when the Cube is deleted so NULL g_incam - g_incam = NULL; - } //end for list.size - - camGrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace); - camGrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace); - camGrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace); - camGrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace); - - - // We want to delete the keywords we just added if the user wants the range - // out of the mapfile, otherwise they will replace any keywords not in the - // mapfile - if (ui.GetString("DEFAULTRANGE") == "MAP") { - camGrp.deleteKeyword("MinimumLatitude"); - camGrp.deleteKeyword("MaximumLatitude"); - camGrp.deleteKeyword("MinimumLongitude"); - camGrp.deleteKeyword("MaximumLongitude"); - } - - // Otherwise, remove the keywords from the map file so the camera keywords - // will be propogated correctly - else { - while (userGrp.hasKeyword("MinimumLatitude")) { - userGrp.deleteKeyword("MinimumLatitude"); - } - while (userGrp.hasKeyword("MinimumLongitude")) { - userGrp.deleteKeyword("MinimumLongitude"); - } - while (userGrp.hasKeyword("MaximumLatitude")) { - userGrp.deleteKeyword("MaximumLatitude"); - } - while (userGrp.hasKeyword("MaximumLongitude")) { - userGrp.deleteKeyword("MaximumLongitude"); - } - } - - // If the user decided to enter a ground range then override - if (ui.WasEntered("MINLON")) { - userGrp.addKeyword(PvlKeyword("MinimumLongitude", - toString(ui.GetDouble("MINLON"))), Pvl::Replace); - } - - if (ui.WasEntered("MAXLON")) { - userGrp.addKeyword(PvlKeyword("MaximumLongitude", - toString(ui.GetDouble("MAXLON"))), Pvl::Replace); - } - - if (ui.WasEntered("MINLAT")) { - userGrp.addKeyword(PvlKeyword("MinimumLatitude", - toString(ui.GetDouble("MINLAT"))), Pvl::Replace); - } - - if (ui.WasEntered("MAXLAT")) { - userGrp.addKeyword(PvlKeyword("MaximumLatitude", - toString(ui.GetDouble("MAXLAT"))), Pvl::Replace); - } - - // If they want the res. from the mapfile, delete it from the camera so - // nothing gets overriden - if (ui.GetString("PIXRES") == "MAP") { - camGrp.deleteKeyword("PixelResolution"); - } - // Otherwise, delete any resolution keywords from the mapfile so the camera - // info is propogated over - else if (ui.GetString("PIXRES") == "CAMERA") { - if (userGrp.hasKeyword("Scale")) { - userGrp.deleteKeyword("Scale"); - } - if (userGrp.hasKeyword("PixelResolution")) { - userGrp.deleteKeyword("PixelResolution"); - } - } - - // Copy any defaults that are not in the user map from the camera map file - for (int k = 0; k < camGrp.keywords(); k++) { - if (!userGrp.hasKeyword(camGrp[k].name())) { - userGrp += camGrp[k]; - } - } - - // If the user decided to enter a resolution then override - if (ui.GetString("PIXRES") == "MPP") { - userGrp.addKeyword(PvlKeyword("PixelResolution", - toString(ui.GetDouble("RESOLUTION"))), - Pvl::Replace); - if (userGrp.hasKeyword("Scale")) { - userGrp.deleteKeyword("Scale"); - } - } - else if (ui.GetString("PIXRES") == "PPD") { - userGrp.addKeyword(PvlKeyword("Scale", - toString(ui.GetDouble("RESOLUTION"))), - Pvl::Replace); - if (userGrp.hasKeyword("PixelResolution")) { - userGrp.deleteKeyword("PixelResolution"); - } - } - - // See if the user want us to handle the longitude seam - if (ui.GetString("DEFAULTRANGE") == "CAMERA" || ui.GetString("DEFAULTRANGE") == "MINIMIZE") { - - // Open the last cube and use its camera - CubeAttributeInput atts0( list.back() ); - Cube icube; - if(atts0.bands().size() != 0) { - vector<QString> lame = atts0.bands(); - icube.setVirtualBands(lame); - } - icube.open( list.back().toString() ); - g_incam = icube.camera(); - - if (g_incam->IntersectsLongitudeDomain(userMap)) { - if (ui.GetString("LONSEAM") == "AUTO") { - if ((int) userGrp["LongitudeDomain"] == 360) { - userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(180)), - Pvl::Replace); - if (g_incam->IntersectsLongitudeDomain(userMap)) { - // Its looks like a global image so switch back to the - // users preference - userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(360)), - Pvl::Replace); - } - } - else { - userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(360)), - Pvl::Replace); - if (g_incam->IntersectsLongitudeDomain(userMap)) { - // Its looks like a global image so switch back to the - // users preference - userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(180)), - Pvl::Replace); - } - } - // Make the target info match the new longitude domain - double minlat, maxlat, minlon, maxlon; - g_incam->GroundRange(minlat, maxlat, minlon, maxlon, userMap); - if (!ui.WasEntered("MINLAT")) { - userGrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace); - } - if (!ui.WasEntered("MAXLAT")) { - userGrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace); - } - if (!ui.WasEntered("MINLON")) { - userGrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace); - } - if (!ui.WasEntered("MAXLON")) { - userGrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace); - } - } - - else if (ui.GetString("LONSEAM") == "ERROR") { - QString msg = "The image [" + ui.GetCubeName("FROM") + "] crosses the " + - "longitude seam"; - throw IException(IException::User, msg, _FILEINFO_); - } - } - - // The camera will be deleted when the Cube is deleted so NULL g_incam - g_incam = NULL; - } - - - Pvl pvl; - pvl.addGroup(userGrp); - - // If there is only one input cube, we need to attach an AlphaCube to the outputs. - if (list.size() == 1) { - Cube parent(list[0].toString()); - if (!parent.hasGroup("AlphaCube")) { - PvlGroup alpha("AlphaCube"); - alpha += PvlKeyword("AlphaSamples", toString(parent.sampleCount())); - alpha += PvlKeyword("AlphaLines", toString(parent.lineCount())); - alpha += PvlKeyword("AlphaStartingSample", toString(0.5)); - alpha += PvlKeyword("AlphaStartingLine", toString(0.5)); - alpha += PvlKeyword("AlphaEndingSample", toString(parent.sampleCount() + 0.5)); - alpha += PvlKeyword("AlphaEndingLine", toString(parent.lineCount() + 0.5)); - alpha += PvlKeyword("BetaSamples", toString(parent.sampleCount())); - alpha += PvlKeyword("BetaLines", toString(parent.lineCount())); - pvl.addGroup(alpha); - } - } - - g_processGroundPolygons.SetStatCubes("TO", pvl, bands); - - bool useCenter = true; - if (ui.GetString("METHOD") == "CENTER") { - useCenter = true; - } - else if (ui.GetString("METHOD") == " WHOLEPIXEL") { - useCenter = false; - } - - g_processGroundPolygons.SetIntersectAlgorithm(useCenter); - - for (int f = 0; f < list.size(); f++) { - - Cube cube(list[f].toString(), "r"); - // Loop through the input cube and get the all pixels values for all bands - ProcessByBrick processBrick; - processBrick.Progress()->SetText("Working on file: " + list[f].toString()); - processBrick.SetBrickSize(1, 1, bands); - // Recall list[f] is a FileName, which stores the attributes - CubeAttributeInput atts0(list[f]); - Cube *icube = processBrick.SetInputCube(list[f].toString(), atts0, 0); - g_incam = icube->camera(); - - processBrick.StartProcess(rasterizePixel); - processBrick.EndProcess(); - } + Pvl appLog; + pixel2map(ui, &appLog); - // When there is only one input cube, we want to propagate IsisCube labels to output cubes - if (list.size() == 1) { - // Note that polygons and original labels are not propagated - g_processGroundPolygons.PropagateLabels(list[0].toString()); - // Tell Process which tables we want to propagate - QList<QString> tablesToPropagate; - tablesToPropagate << "InstrumentPointing" << "InstrumentPosition" << "BodyRotation" - << "SunPosition"; - g_processGroundPolygons.PropagateTables(list[0].toString(), tablesToPropagate); + if( ui.WasEntered("TO") && ui.IsInteractive() ) { + Application::GuiLog(appLog); } - g_processGroundPolygons.EndProcess(); - - // WARNING: rasterizePixel() method alters the current state of the camera. - // If any code is added after this point, you must call setImage to return - // to original camera state before rasterization. - + } @@ -366,56 +36,17 @@ void IsisMain() { * Helper function to print out mapfile to session log */ void PrintMap() { - UserInterface &ui = Application::GetUserInterface(); +//removed in the refactoring process +UserInterface &ui = Application::GetUserInterface(); // Get mapping group from map file - Pvl userMap; - userMap.read(ui.GetFileName("MAP")); + Pvl userMap(ui.GetFileName("MAP")); PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse); //Write map file out to the log +// Isis::Application::GuiLog(userGrp); -} - +} // PrintMap -/** - * This method uses the ProcessGroundPolygons object to rasterize each - * pixel in the given buffer. - * - * @param in Input ProcessByBrick buffer. - */ -void rasterizePixel(Isis::Buffer &in) { - - std::vector<double>lat, lon; - std::vector<double>dns; - - for (int i = 0; i < in.size(); i++) { - dns.push_back(in[i]); - } - - int l = in.Line(); - int s = in.Sample(); - - // TODO: This needs to be done for each band for band dependant instruments - // Note: This can slow this down a lot - - // Get the IFOVs in lat/lon space - PixelFOV fov; - QList< QList< QPointF > > fovVertices = fov.latLonVertices(*g_incam, l, s, g_numIFOVs); + - // loop through each ifov list - for (int ifov = 0; ifov < fovVertices.size(); ifov++) { - // we need at least 3 vertices for a polygon - if (fovVertices[ifov].size() > 3) { - // Get lat/lon for each vertex of the ifov - for (int point = 0; point < fovVertices[ifov].size(); point++) { - lat.push_back(fovVertices[ifov][point].x()); - lon.push_back(fovVertices[ifov][point].y()); - } - // rasterize this ifov and clear vectors for next ifov - g_processGroundPolygons.Rasterize(lat, lon, dns); - lat.clear(); - lon.clear(); - } - } -} diff --git a/isis/src/base/apps/pixel2map/pixel2map.cpp b/isis/src/base/apps/pixel2map/pixel2map.cpp new file mode 100644 index 0000000000000000000000000000000000000000..70e0b4dc19bb6cd860cfc918aaa270bd8d0310cf --- /dev/null +++ b/isis/src/base/apps/pixel2map/pixel2map.cpp @@ -0,0 +1,500 @@ +#include <cmath> + +#include "pixel2map.h" + +#include "Process.h" + +using namespace std; +using namespace Isis; + +namespace Isis { + + static void rasterizePixel(Isis::Buffer &in); + static void vectorizePixel(Isis::Buffer &in); + + // Global variables + static Isis::ProcessGroundPolygons g_processGroundPolygons; + static Isis::Camera *g_incam; + static int g_numIFOVs = 0; + static int g_vectorOut = 0; + + static int csamples; + static int clines; + + static QString vectOut; + static QString outvect; + + static std::ofstream fout_csv; + static QString ogc_SRS; + + void pixel2map(UserInterface &ui, Pvl *appLog){ + + Pvl userMap(ui.GetFileName("MAP")); + PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse); + + FileList list; + if (ui.GetString("FROMTYPE") == "FROM") { + // GetAsString will capture the entire string, including attributes + list.push_back(FileName(ui.GetAsString("FROM"))); + } + else { + try { + list.read(ui.GetFileName("FROMLIST")); + } + catch (IException &e) { + throw IException(e); + } + } + + if (ui.GetString("FOVRANGE") == "INSTANTANEOUS") { + g_numIFOVs = 1; + } + else { + g_numIFOVs = ui.GetInteger("NUMIFOV"); + } + + double newminlat, newmaxlat, newminlon, newmaxlon; + double minlat = 90; + double maxlat = -90; + double minlon = 360.0; + double maxlon = 0; + PvlGroup camGrp; + QString lastBandString; + + // Get the combined lat/lon range for all input cubes + int bands = 1; + for (int i = 0; i < list.size(); i++) { + + // Open the input cube and get the camera + CubeAttributeInput atts0(list[i]); + Cube icube; + if(atts0.bands().size() != 0) { + vector<QString> lame = atts0.bands(); + icube.setVirtualBands(lame); + } + icube.open( list[i].toString() ); + bands = icube.bandCount(); + g_incam = icube.camera(); + + csamples = g_incam->Samples(); + clines = g_incam->Lines(); + + // Preparing the OGC IAU code for the vector data being generated + Spice spi(icube); + ogc_SRS = "IAU:" + Isis::toString(spi.target()->naifBodyCode()) + "00"; + + // Make sure it is not the sky + if (g_incam->target()->isSky()) { + QString msg = "The image [" + list[i].toString() + + "] is targeting the sky, use skymap instead."; + throw IException(IException::User, msg, _FILEINFO_); + } + + // Make sure all the bands for all the files match + if (i >= 1 && atts0.bandsString() != lastBandString) { + QString msg = "The Band numbers for all the files do not match."; + throw IException(IException::User, msg, _FILEINFO_); + } + else { + lastBandString = atts0.bandsString(); + } + + // Get the mapping group + Isis::Pvl camMap; + g_incam->BasicMapping(camMap); + camGrp = camMap.findGroup("Mapping"); + + g_incam->GroundRange(newminlat, newmaxlat, newminlon, newmaxlon, userMap); + //set min lat/lon + if (newminlat < minlat) { + minlat = newminlat; + } + if (newminlon < minlon) { + minlon = newminlon; + } + + //set max lat/lon + if (newmaxlat > maxlat) { + maxlat = newmaxlat; + } + if (newmaxlon > maxlon) { + maxlon = newmaxlon; + } + + // The camera will be deleted when the Cube is deleted so NULL g_incam + g_incam = NULL; + } //end for list.size + + camGrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace); + camGrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace); + camGrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace); + camGrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace); + + + // We want to delete the keywords we just added if the user wants the range + // out of the mapfile, otherwise they will replace any keywords not in the + // mapfile + if (ui.GetString("DEFAULTRANGE") == "MAP") { + camGrp.deleteKeyword("MinimumLatitude"); + camGrp.deleteKeyword("MaximumLatitude"); + camGrp.deleteKeyword("MinimumLongitude"); + camGrp.deleteKeyword("MaximumLongitude"); + } + + // Otherwise, remove the keywords from the map file so the camera keywords + // will be propogated correctly + else { + while (userGrp.hasKeyword("MinimumLatitude")) { + userGrp.deleteKeyword("MinimumLatitude"); + } + while (userGrp.hasKeyword("MinimumLongitude")) { + userGrp.deleteKeyword("MinimumLongitude"); + } + while (userGrp.hasKeyword("MaximumLatitude")) { + userGrp.deleteKeyword("MaximumLatitude"); + } + while (userGrp.hasKeyword("MaximumLongitude")) { + userGrp.deleteKeyword("MaximumLongitude"); + } + } + + // If the user decided to enter a ground range then override + if (ui.WasEntered("MINLON")) { + userGrp.addKeyword(PvlKeyword("MinimumLongitude", + toString(ui.GetDouble("MINLON"))), Pvl::Replace); + } + + if (ui.WasEntered("MAXLON")) { + userGrp.addKeyword(PvlKeyword("MaximumLongitude", + toString(ui.GetDouble("MAXLON"))), Pvl::Replace); + } + + if (ui.WasEntered("MINLAT")) { + userGrp.addKeyword(PvlKeyword("MinimumLatitude", + toString(ui.GetDouble("MINLAT"))), Pvl::Replace); + } + + if (ui.WasEntered("MAXLAT")) { + userGrp.addKeyword(PvlKeyword("MaximumLatitude", + toString(ui.GetDouble("MAXLAT"))), Pvl::Replace); + } + + // If they want the res. from the mapfile, delete it from the camera so + // nothing gets overriden + if (ui.GetString("PIXRES") == "MAP") { + camGrp.deleteKeyword("PixelResolution"); + } + // Otherwise, delete any resolution keywords from the mapfile so the camera + // info is propogated over + else if (ui.GetString("PIXRES") == "CAMERA") { + if (userGrp.hasKeyword("Scale")) { + userGrp.deleteKeyword("Scale"); + } + if (userGrp.hasKeyword("PixelResolution")) { + userGrp.deleteKeyword("PixelResolution"); + } + } + + // Copy any defaults that are not in the user map from the camera map file + for (int k = 0; k < camGrp.keywords(); k++) { + if (!userGrp.hasKeyword(camGrp[k].name())) { + userGrp += camGrp[k]; + } + } + + // If the user decided to enter a resolution then override + if (ui.GetString("PIXRES") == "MPP") { + userGrp.addKeyword(PvlKeyword("PixelResolution", + toString(ui.GetDouble("RESOLUTION"))), + Pvl::Replace); + if (userGrp.hasKeyword("Scale")) { + userGrp.deleteKeyword("Scale"); + } + } + else if (ui.GetString("PIXRES") == "PPD") { + userGrp.addKeyword(PvlKeyword("Scale", + toString(ui.GetDouble("RESOLUTION"))), + Pvl::Replace); + if (userGrp.hasKeyword("PixelResolution")) { + userGrp.deleteKeyword("PixelResolution"); + } + } + + // See if the user want us to handle the longitude seam + if (ui.GetString("DEFAULTRANGE") == "CAMERA" || ui.GetString("DEFAULTRANGE") == "MINIMIZE") { + + // Open the last cube and use its camera + CubeAttributeInput atts0( list.back() ); + Cube icube; + if(atts0.bands().size() != 0) { + vector<QString> lame = atts0.bands(); + icube.setVirtualBands(lame); + } + icube.open( list.back().toString() ); + g_incam = icube.camera(); + + if (g_incam->IntersectsLongitudeDomain(userMap)) { + if (ui.GetString("LONSEAM") == "AUTO") { + if ((int) userGrp["LongitudeDomain"] == 360) { + userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(180)), + Pvl::Replace); + if (g_incam->IntersectsLongitudeDomain(userMap)) { + // Its looks like a global image so switch back to the + // users preference + userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(360)), + Pvl::Replace); + } + } + else { + userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(360)), + Pvl::Replace); + if (g_incam->IntersectsLongitudeDomain(userMap)) { + // Its looks like a global image so switch back to the + // users preference + userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(180)), + Pvl::Replace); + } + } + // Make the target info match the new longitude domain + double minlat, maxlat, minlon, maxlon; + g_incam->GroundRange(minlat, maxlat, minlon, maxlon, userMap); + if (!ui.WasEntered("MINLAT")) { + userGrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace); + } + if (!ui.WasEntered("MAXLAT")) { + userGrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace); + } + if (!ui.WasEntered("MINLON")) { + userGrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace); + } + if (!ui.WasEntered("MAXLON")) { + userGrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace); + } + } + + else if (ui.GetString("LONSEAM") == "ERROR") { + QString msg = "The image [" + ui.GetCubeName("FROM") + "] crosses the " + + "longitude seam"; + throw IException(IException::User, msg, _FILEINFO_); + } + } + + // The camera will be deleted when the Cube is deleted so NULL g_incam + g_incam = NULL; + } + + Pvl pvl; + pvl.addGroup(userGrp); + + // If there is only one input cube, we need to attach an AlphaCube to the outputs. + if (list.size() == 1) { + Cube parent(list[0].toString()); + if (!parent.hasGroup("AlphaCube")) { + PvlGroup alpha("AlphaCube"); + alpha += PvlKeyword("AlphaSamples", toString(parent.sampleCount())); + alpha += PvlKeyword("AlphaLines", toString(parent.lineCount())); + alpha += PvlKeyword("AlphaStartingSample", toString(0.5)); + alpha += PvlKeyword("AlphaStartingLine", toString(0.5)); + alpha += PvlKeyword("AlphaEndingSample", toString(parent.sampleCount() + 0.5)); + alpha += PvlKeyword("AlphaEndingLine", toString(parent.lineCount() + 0.5)); + alpha += PvlKeyword("BetaSamples", toString(parent.sampleCount())); + alpha += PvlKeyword("BetaLines", toString(parent.lineCount())); + pvl.addGroup(alpha); + } + } + + if ( ui.WasEntered("TO") ) { + g_processGroundPolygons.SetStatCubes("TO", pvl, bands); + } + + if ( ui.WasEntered("TOVECT") ) { + g_vectorOut = 1; + outvect = ui.GetFileName("TOVECT"); + } + + bool useCenter = true; + if (ui.GetString("METHOD") == "CENTER") { + useCenter = true; + } + else if (ui.GetString("METHOD") == " WHOLEPIXEL") { + useCenter = false; + } + + g_processGroundPolygons.SetIntersectAlgorithm(useCenter); + + + for (int f = 0; f < list.size(); f++) { + + Cube cube(list[f].toString(), "r"); + // Loop through the input cube and get the all pixels values for all bands + ProcessByBrick processBrick; + processBrick.Progress()->SetText("Working on file: " + list[f].toString()); + processBrick.SetBrickSize(1, 1, bands); + // Recall list[f] is a FileName, which stores the attributes + CubeAttributeInput atts0(list[f]); + Cube *icube = processBrick.SetInputCube(list[f].toString(), atts0, 0); + g_incam = icube->camera(); + + if ( ui.WasEntered("TO") ) { + processBrick.StartProcess(rasterizePixel); + processBrick.EndProcess(); + } + if ( ui.WasEntered("TOVECT") ) { + + ofstream fout_vrt; + + QString outFileNameVRT = FileName( outvect.toLatin1().data() ).removeExtension().addExtension("vrt").expanded(); + QString layer_name = FileName( outvect.toLatin1().data() ).baseName(); + QString csv_fname = FileName( outvect.toLatin1().data() ).name(); + + fout_vrt.open( outFileNameVRT.toLatin1().data() ); + + fout_vrt << "<OGRVRTDataSource>" << endl; + fout_vrt << " <OGRVRTLayer name=\""<< layer_name.toLatin1().data() << "\"> " << endl; + fout_vrt << " <SrcDataSource relativeToVRT=\"1\">" << csv_fname.toLatin1().data() << "</SrcDataSource>" << endl; + fout_vrt << " <GeometryType>wkbPolygon</GeometryType>" << endl; + fout_vrt << " <LayerSRS>"<< ogc_SRS.toLatin1().data() << "</LayerSRS>" << endl; + fout_vrt << " <Field name=\"sample\" src=\"sampleno\" type=\"Integer\"/> "<< endl; + fout_vrt << " <Field name=\"line\" src=\"lineno\" type=\"Integer\"/> "<< endl; + fout_vrt << " <Field name=\"pixelvalue\" src=\"pixelvalue\" type=\"Real\"/>" << endl; + fout_vrt << " <GeometryField encoding=\"WKT\" field=\"geom\" reportSrcColumn=\"FALSE\" />" << endl; + fout_vrt << " </OGRVRTLayer>" << endl; + fout_vrt << "</OGRVRTDataSource>" << endl; + + fout_vrt.close(); + + // write the header + fout_csv.open( outvect.toLatin1().data() ); + fout_csv << "sampleno" << "," << "lineno" << "," << "pixelvalue" << "," << "geom" << endl; + fout_csv.close(); + // re-open in append mode + fout_csv.open( outvect.toLatin1().data(), std::ios_base::app ); + //fout_csv << std::setprecision(7); // see isis2ascii + processBrick.StartProcess(vectorizePixel); + processBrick.EndProcess(); + fout_csv.close(); + + } // IF TOVECT + + } + + if ( ui.WasEntered("TO") ) { + // When there is only one input cube, we want to propagate IsisCube labels to output cubes + if (list.size() == 1) { + // Note that polygons and original labels are not propagated + g_processGroundPolygons.PropagateLabels(list[0].toString()); + // Tell Process which tables we want to propagate + QList<QString> tablesToPropagate; + tablesToPropagate << "InstrumentPointing" << "InstrumentPosition" << "BodyRotation" + << "SunPosition"; + g_processGroundPolygons.PropagateTables(list[0].toString(), tablesToPropagate); + } + } + g_processGroundPolygons.EndProcess(); + + // WARNING: rasterizePixel() method alters the current state of the camera. + // If any code is added after this point, you must call setImage to return + // to original camera state before rasterization. + + } // void pixel2map + + /** + * This method uses the ProcessGroundPolygons object to rasterize each + * pixel in the given buffer. + * + * @param in Input ProcessByBrick buffer. + */ + void rasterizePixel(Isis::Buffer &in) { + + std::vector<double>lat, lon; + std::vector<double>dns; + + for (int i = 0; i < in.size(); i++) { + dns.push_back(in[i]); + } + + int l = in.Line(); + int s = in.Sample(); + + // TODO: This needs to be done for each band for band dependant instruments + // Note: This can slow this down a lot + + // Get the IFOVs in lat/lon space + PixelFOV fov; + QList< QList< QPointF > > fovVertices = fov.latLonVertices(*g_incam, l, s, g_numIFOVs); + + // loop through each ifov list + for (int ifov = 0; ifov < fovVertices.size(); ifov++) { + // we need at least 3 vertices for a polygon + if (fovVertices[ifov].size() > 3) { + // Get lat/lon for each vertex of the ifov + for (int point = 0; point < fovVertices[ifov].size(); point++) { + lat.push_back(fovVertices[ifov][point].x()); + lon.push_back(fovVertices[ifov][point].y()); + } + // rasterize this ifov and clear vectors for next ifov + g_processGroundPolygons.Rasterize(lat, lon, dns); + lat.clear(); + lon.clear(); + } + } + } // rasterizePixel + + /** + * This method uses the ProcessGroundPolygons object to vectorize each + * pixel in the given buffer. + * + * @param in Input ProcessByBrick buffer. + */ + void vectorizePixel(Isis::Buffer &in) { + + std::vector<double>lat, lon; + std::vector<double>dns; + geos::geom::Geometry* GndPixel; + + // Setup the WKT writer + geos::io::WKTWriter *wkt = new geos::io::WKTWriter(); + + // With one band at input, size will be 1 + for (int i = 0; i < in.size(); i++) { + // adds the element at the end of the vector + dns.push_back(in[i]); + } + + int l = in.Line(); + int s = in.Sample(); + + // Get the IFOVs in lat/lon space + PixelFOV fov; + QList< QList< QPointF > > fovVertices = fov.latLonVertices(*g_incam, l, s, g_numIFOVs); + + + // loop through each ifov list + for (int ifov = 0; ifov < fovVertices.size(); ifov++) { + // we need at least 3 vertices for a polygon + if (fovVertices[ifov].size() > 3) { + // Get lat/lon for each vertex of the ifov + for (int point = 0; point < fovVertices[ifov].size(); point++) { + lat.push_back(fovVertices[ifov][point].x()); + lon.push_back(fovVertices[ifov][point].y()); + } + + // vectorize this ifov and clear coordinate vectors for next ifov + GndPixel = g_processGroundPolygons.Vectorize(lat, lon); + + QString v = QString::number(dns[0]); + + // Generate vectors of only non null pixel + if ( dns[0] != Isis::Null ) { + fout_csv << s << "," << l << "," << dns[0] << ",\"" << wkt->write(GndPixel) << "\"" << endl; + } + + lat.clear(); + lon.clear(); + } + + } + + } // vectorizePixel +} // namespace Isis + diff --git a/isis/src/base/apps/pixel2map/pixel2map.h b/isis/src/base/apps/pixel2map/pixel2map.h index 9e616cf2aedfc4a13a4adb8b2ee1dc75ee4f3d9c..93e15cc75eef5494957dc7165ef2ca83d7a16907 100755 --- a/isis/src/base/apps/pixel2map/pixel2map.h +++ b/isis/src/base/apps/pixel2map/pixel2map.h @@ -1,29 +1,32 @@ #ifndef pixel2map_h #define pixel2map_h +#include <QDebug> +#include <QList> +#include <QPointF> +#include <QString> + +#include "Camera.h" +#include "Cube.h" +#include "FileList.h" +#include "PixelFOV.h" +#include "ProcessByBrick.h" +#include "ProcessGroundPolygons.h" +#include "Pvl.h" +#include "Target.h" + #include "Transform.h" -/** - * @author 2008-02-13 Stacy Alley - * - * @internal - * @history 2013-07-30 Stuart C. Sides & Tracie Sucharski Renamed from vim2map - */ -class pixel2map : public Isis::Transform { - private: - Isis::Camera *p_incam; - Isis::Projection *p_outmap; +#include "geos/geom/Geometry.h" +#include "geos/io/WKTWriter.h" - public: - // constructor - pixel2map(const int inputSamples, const int inputLines, Isis::Camera *incam, - const int outputSamples, const int outputLines, Isis::Projection *outmap, - bool trim); +#include "UserInterface.h" - // destructor - ~pixel2map() {}; +namespace Isis { + +extern void pixel2map(UserInterface &ui, Pvl *log=nullptr); -}; +} #endif diff --git a/isis/src/base/apps/pixel2map/pixel2map.xml b/isis/src/base/apps/pixel2map/pixel2map.xml index 0cce8b7b57b583d21270667c51208a55b38063fc..bbefbb8efaf5c070c5dcd79d176f940f2c6c93f8 100755 --- a/isis/src/base/apps/pixel2map/pixel2map.xml +++ b/isis/src/base/apps/pixel2map/pixel2map.xml @@ -323,6 +323,9 @@ <change name="Ian Humphrey" date="2017-05-23"> Updated documentation and added examples. Fixes #4537. </change> + <change name="Alessandro Frigeri" date="2024-08-22"> + Added Vector output + </change> </history> <groups> @@ -418,10 +421,13 @@ </filter> </parameter> + + <parameter name="TO"> <type>cube</type> <fileMode>output</fileMode> - <pixelType>real</pixelType> + <pixelType>real</pixelType> + <default><item>Null</item></default> <brief> Newly mapped cube </brief> @@ -431,10 +437,36 @@ filename with '-count-' appended. The DN values in this output cube indicate how many input pixels were included in the average calculation to create the output pixel. </description> + <exclusions> + <item>TOVECT</item> + </exclusions> <filter> *.cub </filter> </parameter> + + <parameter name="TOVECT"> + <type>filename</type> + <fileMode>output</fileMode> + <default><item>output.csv</item></default> + <brief> + Geospatial CSV/BVRT Vector File + </brief> + <description> + The output file is a comma-delimited file (.csv). An XML VRT sidecar file with geospatial metadata will be generated automatically. The vector pixel will include the values of the first band of the input cube. The coordinate are spherical, unprojected. This file is interoperable with <i>Open Gis Consortium (OGC)</i> compatible geospatial tools and desktop GIS. + </description> + <exclusions> + <item>TO</item> + <item>MAP</item> + <item>FROMLIST</item> + </exclusions> + <filter> + *.csv + </filter> + </parameter> + + + </group> <group name="Output Average"> diff --git a/isis/src/base/objs/ProcessGroundPolygons/ProcessGroundPolygons.cpp b/isis/src/base/objs/ProcessGroundPolygons/ProcessGroundPolygons.cpp index 4f5a1907f635e957ba557286603ae2e77bbe976e..60c03752ecac17ac21c5dfc6701887ba99021579 100644 --- a/isis/src/base/objs/ProcessGroundPolygons/ProcessGroundPolygons.cpp +++ b/isis/src/base/objs/ProcessGroundPolygons/ProcessGroundPolygons.cpp @@ -10,6 +10,7 @@ find files of those names at the top level of this repository. **/ #include "geos/geom/CoordinateSequence.h" #include "geos/geom/LineString.h" +#include "geos/io/WKTWriter.h" #include "Application.h" #include "BoxcarCachingAlgorithm.h" @@ -104,6 +105,100 @@ namespace Isis { } + + +/** + * This method gets called from the application with the lat/lon + * vertices of a polygon. + * + * it returns a geos::geom::Geometry* representation of the + * geometry + * + * @param lat + * @param lon + * + */ +geos::geom::Geometry* ProcessGroundPolygons::Vectorize(std::vector<double> &lat, + std::vector<double> &lon) + //std::vector<double> &values) + { + + + + // Decide if we need to split the poly on the 360 boundry + // Yes, we can do this better. The lat and lon vectors should be passed in as polygons, but + // for now this is reusing older code. + bool crosses = false; + for (unsigned int i = 0; i < lon.size() - 1; i++) { + if (fabs(lon[i] - lon[i+1]) > 180.0) { + crosses = true; + break; + } + } + + + if (crosses) { + // Make a polygon from the lat/lon vectors and split it on 360 + geos::geom::CoordinateSequence pts; + for (unsigned int i = 0; i < lat.size(); i++) { + pts.add(geos::geom::Coordinate(lon[i], lat[i])); + } + pts.add(geos::geom::Coordinate(lon[0], lat[0])); + + geos::geom::Polygon *crossingPoly = Isis::globalFactory->createPolygon( + globalFactory->createLinearRing(pts)).release(); + + geos::geom::MultiPolygon *splitPoly = NULL; + try { + splitPoly = PolygonTools::SplitPolygonOn360(crossingPoly); + } + // Ignore any pixel footprints that could not be split. This should only be pixels + // that contain the pole. + catch (IException &) { + // See leading comment + } + + delete crossingPoly; + + if (splitPoly != NULL) { + // Process the polygons in the split multipolygon as if we were still using the lat/lon vectors + for (unsigned int g = 0; g < splitPoly->getNumGeometries(); ++g) { + const geos::geom::Polygon *poly = + dynamic_cast<const geos::geom::Polygon *>(splitPoly->getGeometryN(g)); + + geos::geom::CoordinateSequence *llcoords = + poly->getExteriorRing()->getCoordinates().release(); + + // Move each coordinate in the exterior ring of this lat/lon polygon to vectors + // Ignore any holes in the polygon + std::vector<double> tlat; + std::vector<double> tlon; + for (unsigned int cord = 0; cord < llcoords->getSize() - 1; ++cord) { + tlon.push_back(llcoords->getAt(cord).x); + tlat.push_back(llcoords->getAt(cord).y); + } + + } + return splitPoly; + } + } + else { // if does not crosses + + // Make a polygon from the lat/lon vectors and split it on 360 + geos::geom::CoordinateSequence pts; + for (unsigned int i = 0; i < lat.size(); i++) { + pts.add(geos::geom::Coordinate(lon[i], lat[i])); + } + pts.add(geos::geom::Coordinate(lon[0], lat[0])); + + geos::geom::Polygon *singlePoly = Isis::globalFactory->createPolygon( + globalFactory->createLinearRing(pts)).release(); + + return singlePoly; + } + +} + /** * This method gets called from the application with the lat/lon * verticies of a polygon along with the band number and the @@ -234,7 +329,7 @@ namespace Isis { /** * This is a method that is called directly from the * application. Using the "TO" parameter we also create a - * count cube name. The we call the overloaded SetStatCubes + * count cube name. Then we call the overloaded SetStatCubes * method above. * * @param parameter diff --git a/isis/src/base/objs/ProcessGroundPolygons/ProcessGroundPolygons.h b/isis/src/base/objs/ProcessGroundPolygons/ProcessGroundPolygons.h index 3ccb449308b67c4fceeccbf86439fe89e2581369..44ef4adf3d07ae3adaefe4fac21e094827741e32 100644 --- a/isis/src/base/objs/ProcessGroundPolygons/ProcessGroundPolygons.h +++ b/isis/src/base/objs/ProcessGroundPolygons/ProcessGroundPolygons.h @@ -73,6 +73,10 @@ namespace Isis { void Rasterize(std::vector<double> &lat, std::vector<double> &lon, int &band, double &value); + + geos::geom::Geometry* Vectorize (std::vector<double> &lat, + std::vector<double> &lon); + //std::vector<double> &values); void EndProcess(); void Finalize(); diff --git a/isis/tests/CameraFixtures.cpp b/isis/tests/CameraFixtures.cpp index 716c40add3ebc02bfe4453c6352a504d75c55f91..eb8fa8c8270f27bb4eb07cde60a51f59cd59d79c 100644 --- a/isis/tests/CameraFixtures.cpp +++ b/isis/tests/CameraFixtures.cpp @@ -988,4 +988,4 @@ namespace Isis { testCube.reset(); } -} \ No newline at end of file +} diff --git a/isis/tests/FunctionalTestsPixel2map.cpp b/isis/tests/FunctionalTestsPixel2map.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1862fffed0fc052f8dc75dd0c2c723e0b544985a --- /dev/null +++ b/isis/tests/FunctionalTestsPixel2map.cpp @@ -0,0 +1,76 @@ +#include <QTextStream> +#include <QStringList> +#include <QTemporaryFile> + +#include "Fixtures.h" +#include "Camera.h" +#include "CameraFixtures.h" +#include "Pvl.h" +#include "PvlGroup.h" +#include "TestUtilities.h" +#include "XmlToPvlTranslationManager.h" +#include "CSVReader.h" +#include "LineManager.h" +#include "Histogram.h" + +#include "UserInterface.h" + +#include "pixel2map.h" + +#include "gtest/gtest.h" + +using namespace Isis; + +static QString PIXEL2MAP_XML = FileName("$ISISROOT/bin/xml/pixel2map.xml").expanded(); + +TEST_F(MroCtxCube, FunctionalTestPixel2mapVector) { + + QString csvFileName = tempDir.path() + "/vect.csv"; + QString vrtFileName = tempDir.path() + "/vect.vrt"; + + QFile csvFile( csvFileName ); + QFile vrtFile( vrtFileName ); + + QVector<QString> args = {"TOVECT="+csvFile.fileName(), "FROM="+ testCube->fileName() }; + + UserInterface options(PIXEL2MAP_XML, args); + + try { + pixel2map(options); + } + catch (IException &e) { + FAIL() << "Unable to open image: " << e.what() << std::endl; + } + + // pre-TEST: Check that we have no null values in the test cube + std::unique_ptr<Histogram> hist (testCube->histogram()); + + EXPECT_EQ(hist->ValidPixels(), testCube->sampleCount()*testCube->lineCount()); + EXPECT_EQ(hist->NullPixels(), 0); + + + // TEST 1a: Check we have both csv and vrt output files + FileName csvFileOut( csvFileName ); + EXPECT_TRUE(csvFileOut.fileExists()); + FileName vrtFileOut( vrtFileName ); + EXPECT_TRUE(vrtFileOut.fileExists()); + + // TEST 1b: Check the output csv file header + CSVReader::CSVAxis csvLine; + CSVReader csvout = CSVReader( csvFileName , false, 0, ','); + + csvLine = csvout.getRow(0); + + EXPECT_PRED_FORMAT2(AssertQStringsEqual, csvLine[0], "sampleno"); + EXPECT_PRED_FORMAT2(AssertQStringsEqual, csvLine[1], "lineno"); + EXPECT_PRED_FORMAT2(AssertQStringsEqual, csvLine[2], "pixelvalue"); + EXPECT_PRED_FORMAT2(AssertQStringsEqual, csvLine[3], "geom"); + + // The number of lines must be equal to the number of pixels + EXPECT_EQ(testCube->sampleCount()*testCube->lineCount(), csvout.rows()-1); + +} + + + + diff --git a/isis/tests/UnitTestProcessGroundPolygons.cpp b/isis/tests/UnitTestProcessGroundPolygons.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c4672b4405fbc6b5ea5a71067631594354160540 --- /dev/null +++ b/isis/tests/UnitTestProcessGroundPolygons.cpp @@ -0,0 +1,55 @@ +#include "ImagePolygon.h" +#include "LineManager.h" +#include "PolygonTools.h" +#include "Preference.h" +#include "ProgramLauncher.h" +#include "Pvl.h" +#include "geos/geom/Point.h" +#include "geos/geom/MultiPolygon.h" +#include "geos/geom/CoordinateSequence.h" + +#include <gtest/gtest.h> + +#include "CameraFixtures.h" +#include "TempFixtures.h" + +#include "ProcessGroundPolygons.h" + +using namespace Isis; + +// Default test to check integrity of the geometry +TEST( UnitTestVectorize, Default ) { + + Isis::ProcessGroundPolygons g_processGroundPolygons; + geos::geom::Geometry* groundpixel; + + std::vector<double> lons = {255.645377, 256.146301, 256.146301, 255.645377}; + std::vector<double> lats = {9.928429, 9.928429, 10.434929, 10.434929}; + + groundpixel = g_processGroundPolygons.Vectorize(lats, lons); + + // Check that the geometry is valid + EXPECT_EQ( groundpixel->isValid(), 1 ); + EXPECT_EQ( groundpixel->isPolygonal(), 1 ); + EXPECT_EQ( groundpixel->getGeometryType(), "Polygon"); + EXPECT_EQ( groundpixel->getNumGeometries(), 1); + EXPECT_EQ( groundpixel->getNumPoints(), lons.size() + 1 ); + +}; + +TEST( UnitTestVectorize, Crosses360 ) { + + Isis::ProcessGroundPolygons g_processGroundPolygons; + geos::geom::Geometry* groundpixel; + + std::vector<double> lons = {359.0, 1.0, 1.0, 359.0, 359.0}; + std::vector<double> lats = { 0.0, 0.0, 1.0, 1.0, 0.0}; + + groundpixel = g_processGroundPolygons.Vectorize(lats, lons); + + // Check that the geometry is valid + EXPECT_EQ( groundpixel->isValid(), 1 ); + EXPECT_EQ( groundpixel->isPolygonal(), 1 ); + EXPECT_EQ( groundpixel->getGeometryType(), "MultiPolygon"); + EXPECT_EQ( groundpixel->getNumGeometries(), 2); +}; \ No newline at end of file