diff --git a/CMakeLists.txt b/CMakeLists.txt index 46009a88..f59f240d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -227,6 +227,16 @@ if(C_API) add_definitions(-DC_API) endif(C_API) +# Enable phone home and message server startup +option(PHONE_HOME_QUERIES "Enable Phone Home Queries" ON) + +if (PHONE_HOME_QUERIES) + add_definitions(-DPHONE_HOME_QUERIES_ENABLED) + message(STATUS "Phone Home Queries are enabled.") +else() + message(STATUS "Phone Home Queries are disabled.") +endif() + # OpenFOAM related options option(NINJAFOAM "Enable OpenFOAM solver" ON) if(NINJAFOAM) @@ -421,6 +431,7 @@ set(share_files date_time_zonespec.csv map.htm thredds.csv surface_data.zip + srtm_region.geojson us_srtm_region.dbf us_srtm_region.prj us_srtm_region.shp diff --git a/CONTRIBUTORS b/CONTRIBUTORS index 77376d92..3a63f6e5 100644 --- a/CONTRIBUTORS +++ b/CONTRIBUTORS @@ -27,3 +27,6 @@ Nicolás Chiaraviglio Peter Mehlitz Tanner Finney Tomàs Margalef +Hannah Gedlaman +Nicholas Kim +Rui Zhang diff --git a/src/cli/cmake_cli.cpp b/src/cli/cmake_cli.cpp index 68c91c47..d9e39b8e 100644 --- a/src/cli/cmake_cli.cpp +++ b/src/cli/cmake_cli.cpp @@ -51,7 +51,7 @@ int main(int argc, char *argv[]) { CPLSetConfigOption( "NINJA_DISABLE_CALL_HOME", "ON" ); - NinjaInitialize(); + NinjaInitialize("cli"); //set a different key for the CLI since this is where users are likely to abuse access if(CPLGetConfigOption("NINJA_CLI_SRTM_API_KEY", NULL) != NULL) { diff --git a/src/gui/WidgetDownloadDEM.cpp b/src/gui/WidgetDownloadDEM.cpp index f7c13187..4da7acaa 100644 --- a/src/gui/WidgetDownloadDEM.cpp +++ b/src/gui/WidgetDownloadDEM.cpp @@ -302,7 +302,7 @@ void WidgetDownloadDEM::updateProgress() if(result < 0) { - progressBar->setLabelText("The surface data download failed. \nThis normally happens when either the data source doesn't cover your region or the server that provides the surface data is down or under high usage. \nPlease try again later or try a different data source."); + progressBar->setLabelText("The surface data download failed. \nThis can happen when either the data source doesn't cover your region or the server that provides the surface data is down or under high usage. \nPlease try again later or try a different data source."); progressBar->setRange(0,1); progressBar->setValue( 0 ); progressBar->setCancelButtonText("Close"); diff --git a/src/gui/cmake_gui.cpp b/src/gui/cmake_gui.cpp index 7ff641cb..f82f6d4d 100644 --- a/src/gui/cmake_gui.cpp +++ b/src/gui/cmake_gui.cpp @@ -67,7 +67,7 @@ int main(int argc, char *argv[]) return result; } - NinjaInitialize(); + NinjaInitialize("gui"); QApplication app(argc, argv); diff --git a/src/gui/mainWindow.cpp b/src/gui/mainWindow.cpp index e89042a2..4b3c23aa 100644 --- a/src/gui/mainWindow.cpp +++ b/src/gui/mainWindow.cpp @@ -84,6 +84,10 @@ mainWindow::mainWindow(QWidget *parent) meshCellSize = 200.0; +#ifdef PHONE_HOME_QUERIES_ENABLED + checkMessages(); +#endif + QString v(NINJA_VERSION_STRING); v = "Welcome to WindNinja " + v; @@ -113,32 +117,30 @@ mainWindow::mainWindow(QWidget *parent) /* ** Check for version updates, or messages from the server. */ +#ifdef PHONE_HOME_QUERIES_ENABLED void mainWindow::checkMessages(void) { - QMessageBox mbox; - char **papszMsg = NinjaCheckVersion(); - if (papszMsg != NULL) { - const char *pszVers = - CSLFetchNameValueDef(papszMsg, "VERSION", NINJA_VERSION_STRING); - if (strcmp(pszVers, NINJA_VERSION_STRING) > 0) { - mbox.setText("A new version of WindNinja is available: " + - QString(pszVers)); - mbox.exec(); - } - char **papszUserMsg = CSLFetchNameValueMultiple(papszMsg, "MESSAGE"); - for (int i = 0; i < CSLCount(papszUserMsg); i++) { - mbox.setText(QString(papszUserMsg[i])); - mbox.exec(); - } - CSLDestroy(papszUserMsg); - if (CSLFetchNameValue(papszMsg, "ABORT") != NULL) { + QMessageBox mbox; + char *papszMsg = NinjaQueryServerMessages(true); + if (papszMsg != NULL) { + if (strcmp(papszMsg, "TRUE\n") == 0) { mbox.setText("There is a fatal flaw in Windninja, it must close."); mbox.exec(); abort(); - } } - CSLDestroy(papszMsg); + + else { + char *papszMsg = NinjaQueryServerMessages(false); + if (papszMsg != NULL) { + mbox.setText(papszMsg); + + mbox.exec(); + } + } + } } +#endif + bool mainWindow::okToContinue() { if(okToContinueCheck) @@ -1223,8 +1225,6 @@ void mainWindow::aboutWindNinja() "Kyle Shannon
" \ "Natalie Wagenbrenner
" \ "Bret Butler
" \ - "Levi Malott
" \ - "Cody Posey

"); aboutText.append("

Missoula Fire Sciences Laboratory
"); aboutText.append("Rocky Mountain Research Station
"); aboutText.append("USDA Forest Service
"); diff --git a/src/gui/mainWindow.h b/src/gui/mainWindow.h index 16434d5f..647bb27d 100644 --- a/src/gui/mainWindow.h +++ b/src/gui/mainWindow.h @@ -145,7 +145,10 @@ class mainWindow : public QMainWindow void mainDiurnalChanged(bool dC); public slots: + +#ifdef PHONE_HOME_QUERIES_ENABLED void checkMessages(); +#endif #ifdef NINJAFOAM void openExistingCase(); void updateFileInputForCase(const char* file); diff --git a/src/ninja/KmlVector.cpp b/src/ninja/KmlVector.cpp index a9da41be..65549923 100644 --- a/src/ninja/KmlVector.cpp +++ b/src/ninja/KmlVector.cpp @@ -43,6 +43,9 @@ KmlVector::KmlVector() wxModelName = ""; coordTransform = NULL; turbulenceFlag = false; + colMaxFlag = false; + colMax_colHeightAGL = -1.0; + colMax_colHeightAGL_units = lengthUnits::meters; } KmlVector::~KmlVector() @@ -77,6 +80,14 @@ void KmlVector::setTurbulenceGrid(AsciiGrid &turb, velocityUnits::eVeloc turbulence = turb; } +void KmlVector::setColMaxGrid(AsciiGrid &columnMax, velocityUnits::eVelocityUnits units, const double colHeightAGL, const lengthUnits::eLengthUnits colHeightAGL_units) +{ + speedUnits = units; + colMax = columnMax; + + colMax_colHeightAGL = colHeightAGL; + colMax_colHeightAGL_units = colHeightAGL_units; +} #ifdef FRICTION_VELOCITY void KmlVector::setUstarGrid(AsciiGrid &ust) @@ -280,6 +291,14 @@ bool KmlVector::writeKml(std::string cScheme, bool vector_scaling) writeTurbulence(fout); VSIFPrintfL(fout, ""); } + + if(colMaxFlag) + { + VSIFPrintfL(fout, ""); + VSIFPrintfL(fout, "\n\tColumn Max Velocity Fluctuations\n"); + writeColMax(fout); + VSIFPrintfL(fout, ""); + } #ifdef FRICTION_VELOCITY if(ustarFlag ==true) @@ -366,6 +385,14 @@ bool KmlVector::writeKml(egoogSpeedScaling scaling, string cScheme,bool vector_s VSIFPrintfL(fout, ""); } + if(colMaxFlag) + { + VSIFPrintfL(fout, ""); + VSIFPrintfL(fout, "\n\tColumn Max Velocity Fluctuations\n"); + writeColMax(fout); + VSIFPrintfL(fout, ""); + } + #ifdef FRICTION_VELOCITY ustarFlag = false; @@ -1181,30 +1208,10 @@ bool KmlVector::writeScreenOverlayDateTimeLegendWxModelRun(VSILFILE *fileOut) bool KmlVector::writeTurbulence(VSILFILE *fileOut) { - double xPoint, yPoint; - double xCenter, yCenter; - double left_x, right_x, lower_y, upper_y; - double u = 0; - double cSize; - int nR; - int nC; - double upper, lower, upper_mid, lower_mid, mid; - std::string icon; - - turbulence_png = "turbulence_png.png"; - - cSize = turbulence.get_cellSize(); - nR = turbulence.get_nRows(); - nC = turbulence.get_nCols(); + turbulence_png = "turbulence.png"; - lower = turbulence.get_minValue(); - upper = turbulence.get_maxValue(); - lower_mid = lower + (turbulence.get_maxValue() - turbulence.get_minValue())/4; - upper_mid = upper - (turbulence.get_maxValue() - turbulence.get_minValue())/4; - mid = upper_mid - (turbulence.get_maxValue() - turbulence.get_minValue())/4; - - //---------------make single png for overlay------------------ - std::string outFilename = "turbulence_png.png"; + //---------------make single png for overlay------------------ + std::string outFilename = turbulence_png; std::string scalarLegendFilename = "turbulence_legend"; std::string legendTitle = "Speed Fluctuation"; std::string legendUnits = ""; @@ -1227,19 +1234,56 @@ bool KmlVector::writeTurbulence(VSILFILE *fileOut) break; } bool writeLegend = TRUE; + bool keepTiff = TRUE; // always true for turbulence outputs, input.writeTurbulence was already triggered to get here - turbulence.ascii2png(outFilename, legendTitle, legendUnits, scalarLegendFilename, writeLegend); + turbulence.ascii2png(outFilename, legendTitle, legendUnits, scalarLegendFilename, writeLegend, keepTiff); - turbulence.get_cellPosition(0, 0, &xCenter, &yCenter); //sw - left_x = xCenter - cSize/2; //west - lower_y = yCenter - cSize/2; //south - turbulence.get_cellPosition(nR-1, nC-1, &xCenter, &yCenter); //ne - right_x = xCenter + cSize/2; //east - upper_y = yCenter + cSize/2; //north + if ( keepTiff == TRUE ) + { + std::string base_outFilename = outFilename; + int pos = outFilename.find_last_of(".png"); + if(pos != -1) + base_outFilename = outFilename.substr(0, pos - 4 + 1); // .png is 4 letters back, + 1 to go from digit Id to a count + std::string ascii_tiff_filename = base_outFilename + ".tif"; + + std::string kmz_baseFilename = kmzFile; + pos = kmzFile.find_last_of(".kmz"); + if(pos != -1) + kmz_baseFilename = kmzFile.substr(0, pos - 4 + 1); // .kmz is 4 letters back, + 1 to go from digit Id to a count + std::string kmz_tiff_filename = kmz_baseFilename + "__" + base_outFilename + ".tif"; + + //// now move the keepTiff file to the updated filename + CPLMoveFile( kmz_tiff_filename.c_str(), ascii_tiff_filename.c_str() ); + } - coordTransform->Transform(1, &right_x, &upper_y); - coordTransform->Transform(1, &left_x, &lower_y); - coordTransform->Transform(1, &xCenter, &yCenter); + double left_x = turbulence.get_xllCorner(); + double lower_y = turbulence.get_yllCorner(); + double right_x = turbulence.get_xDimension()+turbulence.get_xllCorner(); + double upper_y = turbulence.get_yDimension()+turbulence.get_yllCorner(); + double xll = left_x; + double yll = lower_y; + double xul = left_x; + double yul = upper_y; + double xur = right_x; + double yur = upper_y; + double xlr = right_x; + double ylr = lower_y; + + coordTransform->Transform(1, &xll, &yll); + coordTransform->Transform(1, &xul, &yul); + coordTransform->Transform(1, &xur, &yur); + coordTransform->Transform(1, &xlr, &ylr); + + double north = std::max(yul,yur); + double south = std::min(yll,ylr); + // calculating for east and west gets more complicated for the rare case that it crosses between -180 and 180 degrees + double east = std::max(xlr,xur); + double west = std::min(xll,xul); + // check if crosses between -180 and 180 degrees, swap min for max or max for min if swapped dirs around circle + if ( std::max(xlr,xur) - std::min(xlr,xur) > 180 ) + east = std::min(xlr,xur); + if ( std::max(xll,xul) - std::min(xll,xul) > 180 ) + west = std::max(xll,xul); int pos; std::string shortName; @@ -1262,21 +1306,21 @@ bool KmlVector::writeTurbulence(VSILFILE *fileOut) //VSIFPrintfL(fileOut, "\n\t\t88ffffff"); VSIFPrintfL(fileOut, "\n\t"); - VSIFPrintfL(fileOut, "\n\t\t%s", shortName.c_str()); //turbulence_png.png + VSIFPrintfL(fileOut, "\n\t\t%s", shortName.c_str()); //turbulence.png VSIFPrintfL(fileOut, "\n\t"); VSIFPrintfL(fileOut, "\n\t"); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", upper_y); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", lower_y); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", right_x); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", left_x); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", north); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", south); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", east); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", west); VSIFPrintfL(fileOut, "\n\t\t0"); VSIFPrintfL(fileOut, "\n\t"); VSIFPrintfL(fileOut, "\n\n"); //---add legend---------------------------------------------- - turbulence_legend = "./turbulence_legend"; + turbulence_legend = CPLSPrintf("./%s",scalarLegendFilename.c_str()); //int pos; //std::string shortName; pos = turbulence_legend.find_last_of('\\'); @@ -1302,50 +1346,189 @@ bool KmlVector::writeTurbulence(VSILFILE *fileOut) return true; } -#ifdef FRICTION_VELOCITY -bool KmlVector::writeUstar(FILE *fileOut) +bool KmlVector::writeColMax(VSILFILE *fileOut) { - double xPoint, yPoint; - double xCenter, yCenter; - double left_x, right_x, lower_y, upper_y; - double u = 0; - double cSize; - int nR; - int nC; - double upper, lower, upper_mid, lower_mid, mid; - std::string icon; + colMax_png = CPLSPrintf("colMax_%ld%sColHeightAGL.png", (long) (colMax_colHeightAGL+0.5), lengthUnits::getString(colMax_colHeightAGL_units).c_str()); + + //---------------make single png for overlay------------------ + std::string outFilename = colMax_png; + std::string scalarLegendFilename = "colMax_legend"; + std::string legendTitle = "Speed Fluctuation"; + std::string legendUnits = ""; + switch(speedUnits) + { + case velocityUnits::metersPerSecond: // m/s + legendUnits = "(m/s)"; + break; + case velocityUnits::milesPerHour: // mph + legendUnits = "(mph)"; + break; + case velocityUnits::kilometersPerHour: // kph + legendUnits = "(kph)"; + break; + case velocityUnits::knots: // kts + legendUnits = "(knots)"; + break; + default: // default is mph + legendUnits = "(mph)"; + break; + } + bool writeLegend = TRUE; + bool keepTiff = TRUE; // always true for turbulence outputs, input.writeTurbulence was already triggered to get here + + colMax.ascii2png(outFilename, legendTitle, legendUnits, scalarLegendFilename, writeLegend, keepTiff); + + if ( keepTiff == TRUE ) + { + std::string base_outFilename = outFilename; + int pos = outFilename.find_last_of(".png"); + if(pos != -1) + base_outFilename = outFilename.substr(0, pos - 4 + 1); // .png is 4 letters back, + 1 to go from digit Id to a count + std::string ascii_tiff_filename = base_outFilename + ".tif"; + + std::string kmz_baseFilename = kmzFile; + pos = kmzFile.find_last_of(".kmz"); + if(pos != -1) + kmz_baseFilename = kmzFile.substr(0, pos - 4 + 1); // .kmz is 4 letters back, + 1 to go from digit Id to a count + std::string kmz_tiff_filename = kmz_baseFilename + "__" + base_outFilename + ".tif"; + + //// now move the keepTiff file to the updated filename + CPLMoveFile( kmz_tiff_filename.c_str(), ascii_tiff_filename.c_str() ); + } + + double left_x = colMax.get_xllCorner(); + double lower_y = colMax.get_yllCorner(); + double right_x = colMax.get_xDimension()+colMax.get_xllCorner(); + double upper_y = colMax.get_yDimension()+colMax.get_yllCorner(); + double xll = left_x; + double yll = lower_y; + double xul = left_x; + double yul = upper_y; + double xur = right_x; + double yur = upper_y; + double xlr = right_x; + double ylr = lower_y; + + coordTransform->Transform(1, &xll, &yll); + coordTransform->Transform(1, &xul, &yul); + coordTransform->Transform(1, &xur, &yur); + coordTransform->Transform(1, &xlr, &ylr); + + double north = std::max(yul,yur); + double south = std::min(yll,ylr); + // calculating for east and west gets more complicated for the rare case that it crosses between -180 and 180 degrees + double east = std::max(xlr,xur); + double west = std::min(xll,xul); + // check if crosses between -180 and 180 degrees, swap min for max or max for min if swapped dirs around circle + if ( std::max(xlr,xur) - std::min(xlr,xur) > 180 ) + east = std::min(xlr,xur); + if ( std::max(xll,xul) - std::min(xll,xul) > 180 ) + west = std::max(xll,xul); + + int pos; + std::string shortName; + pos = colMax_png.find_last_of('\\'); + if(pos == -1) + pos = colMax_png.find_last_of('/'); + + shortName = colMax_png.substr(pos + 1, colMax_png.size()); - ustar_png = "ustar_png.png"; + VSIFPrintfL(fileOut, ""); + VSIFPrintfL(fileOut, "\n\tColumn Max Velocity Fluctuations"); + VSIFPrintfL(fileOut, "\n\t"); + VSIFPrintfL(fileOut, "\n\t\t"); + VSIFPrintfL(fileOut, "\n\t\t\t2"); + VSIFPrintfL(fileOut, "\n\t\t"); + VSIFPrintfL(fileOut, "\n\t"); + + VSIFPrintfL(fileOut, "\n\t0"); + VSIFPrintfL(fileOut, "\n\tclampToGround"); + //VSIFPrintfL(fileOut, "\n\t\t88ffffff"); + + VSIFPrintfL(fileOut, "\n\t"); + VSIFPrintfL(fileOut, "\n\t\t%s", shortName.c_str()); //colMax.png + VSIFPrintfL(fileOut, "\n\t"); + + VSIFPrintfL(fileOut, "\n\t"); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", north); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", south); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", east); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", west); + VSIFPrintfL(fileOut, "\n\t\t0"); + VSIFPrintfL(fileOut, "\n\t"); + + VSIFPrintfL(fileOut, "\n\n"); + + //---add legend---------------------------------------------- + colMax_legend = CPLSPrintf("./%s",scalarLegendFilename.c_str()); + //int pos; + //std::string shortName; + pos = colMax_legend.find_last_of('\\'); + if(pos == -1) + pos = colMax_legend.find_last_of('/'); + + shortName = colMax_legend.substr(pos + 1, colMax_legend.size()); + + VSIFPrintfL(fileOut, ""); + VSIFPrintfL(fileOut, "\nLegend"); + VSIFPrintfL(fileOut, "\n1"); + VSIFPrintfL(fileOut, "\n9bffffff"); + VSIFPrintfL(fileOut, "\n"); + VSIFPrintfL(fileOut, "\n"); + VSIFPrintfL(fileOut, "\n%s", shortName.c_str()); + VSIFPrintfL(fileOut, "\n"); + VSIFPrintfL(fileOut, "\n"); + VSIFPrintfL(fileOut, "\n"); + VSIFPrintfL(fileOut, "\n"); + VSIFPrintfL(fileOut, "\n"); + VSIFPrintfL(fileOut, "\n\n"); - cSize = ustar.get_cellSize(); - nR = ustar.get_nRows(); - nC = ustar.get_nCols(); + return true; +} - lower = ustar.get_minValue(); - upper = ustar.get_maxValue(); - lower_mid = lower + (ustar.get_maxValue() - ustar.get_minValue())/4; - upper_mid = upper - (ustar.get_maxValue() - ustar.get_minValue())/4; - mid = upper_mid - (ustar.get_maxValue() - ustar.get_minValue())/4; +#ifdef FRICTION_VELOCITY +bool KmlVector::writeUstar(FILE *fileOut) +{ + ustar_png = "ustar.png"; - //---------------make single png for overlay------------------ - std::string outFilename = "ustar_png.png"; + //---------------make single png for overlay------------------ + std::string outFilename = ustar_png; std::string scalarLegendFilename = "ustar_legend"; std::string legendTitle = "Friction Velocity"; std::string legendUnits = "(m/s)"; bool writeLegend = TRUE; - - ustar.ascii2png(outFilename, legendTitle, legendUnits, scalarLegendFilename, writeLegend); - - ustar.get_cellPosition(0, 0, &xCenter, &yCenter); //sw - left_x = xCenter - cSize/2; //west - lower_y = yCenter - cSize/2; //south - ustar.get_cellPosition(nR-1, nC-1, &xCenter, &yCenter); //ne - right_x = xCenter + cSize/2; //east - upper_y = yCenter + cSize/2; //north - - coordTransform->Transform(1, &right_x, &upper_y); - coordTransform->Transform(1, &left_x, &lower_y); - coordTransform->Transform(1, &xCenter, &yCenter); + bool keepTiff = FALSE; + + ustar.ascii2png(outFilename, legendTitle, legendUnits, scalarLegendFilename, writeLegend, keepTiff); + + double left_x = turbulence.get_xllCorner(); + double lower_y = turbulence.get_yllCorner(); + double right_x = turbulence.get_xDimension()+turbulence.get_xllCorner(); + double upper_y = turbulence.get_yDimension()+turbulence.get_yllCorner(); + double xll = left_x; + double yll = lower_y; + double xul = left_x; + double yul = upper_y; + double xur = right_x; + double yur = upper_y; + double xlr = right_x; + double ylr = lower_y; + + coordTransform->Transform(1, &xll, &yll); + coordTransform->Transform(1, &xul, &yul); + coordTransform->Transform(1, &xur, &yur); + coordTransform->Transform(1, &xlr, &ylr); + + double north = std::max(yul,yur); + double south = std::min(yll,ylr); + // calculating for east and west gets more complicated for the rare case that it crosses between -180 and 180 degrees + double east = std::max(xlr,xur); + double west = std::min(xll,xul); + // check if crosses between -180 and 180 degrees, swap min for max or max for min if swapped dirs around circle + if ( std::max(xlr,xur) - std::min(xlr,xur) > 180 ) + east = std::min(xlr,xur); + if ( std::max(xll,xul) - std::min(xll,xul) > 180 ) + west = std::max(xll,xul); int pos; std::string shortName; @@ -1368,21 +1551,21 @@ bool KmlVector::writeUstar(FILE *fileOut) //VSIFPrintfL(fileOut, "\n\t\t88ffffff"); VSIFPrintfL(fileOut, "\n\t"); - VSIFPrintfL(fileOut, "\n\t\t%s", shortName.c_str()); //ustar_png.png + VSIFPrintfL(fileOut, "\n\t\t%s", shortName.c_str()); //ustar.png VSIFPrintfL(fileOut, "\n\t"); VSIFPrintfL(fileOut, "\n\t"); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", upper_y); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", lower_y); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", right_x); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", left_x); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", north); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", south); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", east); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", west); VSIFPrintfL(fileOut, "\n\t\t0"); VSIFPrintfL(fileOut, "\n\t"); VSIFPrintfL(fileOut, "\n\n"); //---add legend---------------------------------------------- - ustar_legend = "./ustar_legend"; + ustar_legend = CPLSPrintf("./%s",scalarLegendFilename.c_str()); //int pos; //std::string shortName; pos = ustar_legend.find_last_of('\\'); @@ -1412,48 +1595,46 @@ bool KmlVector::writeUstar(FILE *fileOut) #ifdef EMISSIONS bool KmlVector::writeDust(FILE *fileOut) { - double xPoint, yPoint; - double xCenter, yCenter; - double left_x, right_x, lower_y, upper_y; - double u = 0; - double cSize; - int nR; - int nC; - double upper, lower, upper_mid, lower_mid, mid; - std::string icon; - - cSize = dust.get_cellSize(); - nR = dust.get_nRows(); - nC = dust.get_nCols(); - - lower = dust.get_minValue(); - upper = dust.get_maxValue(); - lower_mid = lower + (dust.get_maxValue() - dust.get_minValue())/4; - upper_mid = upper - (dust.get_maxValue() - dust.get_minValue())/4; - mid = upper_mid - (dust.get_maxValue() - dust.get_minValue())/4; + dust_png = "dust.png"; //---------------make png for overlay------------------ - std::string outFilename = "dust_png.png"; + std::string outFilename = dust_png; std::string scalarLegendFilename = "dust_legend"; std::string legendTitle = "PM10"; std::string legendUnits = "(mg/m2/s)"; bool writeLegend = true; - - dust.ascii2png(outFilename, legendTitle, legendUnits, scalarLegendFilename, writeLegend); - - dust_png = "dust_png.png"; - - dust.get_cellPosition(0, 0, &xCenter, &yCenter); //sw - left_x = xCenter - cSize/2; //west - lower_y = yCenter - cSize/2; //south - dust.get_cellPosition(nR-1, nC-1, &xCenter, &yCenter); //ne - right_x = xCenter + cSize/2; //east - upper_y = yCenter + cSize/2; //north - - - coordTransform->Transform(1, &right_x, &upper_y); - coordTransform->Transform(1, &left_x, &lower_y); - coordTransform->Transform(1, &xCenter, &yCenter); + bool keepTiff = FALSE; + + dust.ascii2png(outFilename, legendTitle, legendUnits, scalarLegendFilename, writeLegend, keepTiff); + + double left_x = turbulence.get_xllCorner(); + double lower_y = turbulence.get_yllCorner(); + double right_x = turbulence.get_xDimension()+turbulence.get_xllCorner(); + double upper_y = turbulence.get_yDimension()+turbulence.get_yllCorner(); + double xll = left_x; + double yll = lower_y; + double xul = left_x; + double yul = upper_y; + double xur = right_x; + double yur = upper_y; + double xlr = right_x; + double ylr = lower_y; + + coordTransform->Transform(1, &xll, &yll); + coordTransform->Transform(1, &xul, &yul); + coordTransform->Transform(1, &xur, &yur); + coordTransform->Transform(1, &xlr, &ylr); + + double north = std::max(yul,yur); + double south = std::min(yll,ylr); + // calculating for east and west gets more complicated for the rare case that it crosses between -180 and 180 degrees + double east = std::max(xlr,xur); + double west = std::min(xll,xul); + // check if crosses between -180 and 180 degrees, swap min for max or max for min if swapped dirs around circle + if ( std::max(xlr,xur) - std::min(xlr,xur) > 180 ) + east = std::min(xlr,xur); + if ( std::max(xll,xul) - std::min(xll,xul) > 180 ) + west = std::max(xll,xul); std::string shortName; int pos; @@ -1475,21 +1656,22 @@ bool KmlVector::writeDust(FILE *fileOut) //VSIFPrintfL(fileOut, "\n\t\t88ffffff"); VSIFPrintfL(fileOut, "\n\t"); - VSIFPrintfL(fileOut, "\n\t\t%s", shortName.c_str()); //dust_png.png + VSIFPrintfL(fileOut, "\n\t\t%s", shortName.c_str()); //dust.png VSIFPrintfL(fileOut, "\n\t"); VSIFPrintfL(fileOut, "\n\t"); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", upper_y); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", lower_y); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", right_x); - VSIFPrintfL(fileOut, "\n\t\t%.10lf", left_x); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", north); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", south); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", east); + VSIFPrintfL(fileOut, "\n\t\t%.10lf", west); VSIFPrintfL(fileOut, "\n\t\t0"); VSIFPrintfL(fileOut, "\n\t"); VSIFPrintfL(fileOut, "\n\n"); //---add legend---------------------------------------------- - dust_legend = "dust_legend"; + dust_legend = scalarLegendFilename; + //dust_legend = CPLSPrintf("./%s",scalarLegendFilename.c_str()); pos = dust_legend.find_last_of('\\'); if(pos == -1) pos = dust_legend.find_last_of('/'); @@ -1721,6 +1903,14 @@ bool KmlVector::makeKmz() filesInZip.push_back(getShortName(turbulence_png)); filesInZip.push_back(getShortName(turbulence_legend)); } + + if(colMaxFlag) + { + filesToZip.push_back(colMax_png); + filesToZip.push_back(colMax_legend); + filesInZip.push_back(getShortName(colMax_png)); + filesInZip.push_back(getShortName(colMax_legend)); + } #ifdef FRICTION_VELOCITY if(ustarFlag==1) @@ -1799,6 +1989,13 @@ bool KmlVector::removeKmlFile() if(turbulence_legend.c_str() !="") VSIUnlink(turbulence_legend.c_str()); + if(colMax_png.c_str() != ""){ + VSIUnlink(colMax_png.c_str()); + VSIUnlink((colMax_png + ".aux.xml").c_str()); + } + if(colMax_legend.c_str() !="") + VSIUnlink(colMax_legend.c_str()); + #ifdef FRICTION_VELOCITY if(ustar_png.c_str() != ""){ VSIUnlink(ustar_png.c_str()); diff --git a/src/ninja/KmlVector.h b/src/ninja/KmlVector.h index 9e790bc8..b82c0ad7 100644 --- a/src/ninja/KmlVector.h +++ b/src/ninja/KmlVector.h @@ -86,6 +86,7 @@ class KmlVector AsciiGrid spd, dir; AsciiGrid turbulence; + AsciiGrid colMax; #ifdef FRICTION_VELOCITY AsciiGrid ustar; #endif @@ -123,6 +124,7 @@ class KmlVector bool writeVectors(VSILFILE *fileOut); bool writeTurbulence(VSILFILE *fileOut); + bool writeColMax(VSILFILE *fileOut); #ifdef FRICTION_VELOCITY bool writeUstar(FILE *fileOut); #endif @@ -142,6 +144,10 @@ class KmlVector void setDirGrid(AsciiGrid &d); void setTurbulenceGrid(AsciiGrid &turb, velocityUnits::eVelocityUnits units); void setTurbulenceFlag(bool inputTurbulenceFlag){turbulenceFlag = inputTurbulenceFlag;} + double colMax_colHeightAGL; + lengthUnits::eLengthUnits colMax_colHeightAGL_units; + void setColMaxGrid(AsciiGrid &columnMax, velocityUnits::eVelocityUnits units, const double colHeightAGL, const lengthUnits::eLengthUnits colHeightAGL_units); + void setColMaxFlag(bool inputColMaxFlag){colMaxFlag = inputColMaxFlag;} #ifdef FRICTION_VELOCITY void setUstarGrid(AsciiGrid &ust); void setUstarFlag(bool inputUstarFlag){ustarFlag = inputUstarFlag;} @@ -178,8 +184,11 @@ class KmlVector bool turbulenceFlag; std::string turbulence_tiff; std::string turbulence_png; - std::string turbulence_legend; - std::string ustar_legend; + std::string turbulence_legend; + bool colMaxFlag; + std::string colMax_tiff; + std::string colMax_png; + std::string colMax_legend; #ifdef FRICTION_VELOCITY bool ustarFlag; std::string ustar_tiff; diff --git a/src/ninja/WindNinjaInputs.cpp b/src/ninja/WindNinjaInputs.cpp index 5fff9386..28807d4f 100644 --- a/src/ninja/WindNinjaInputs.cpp +++ b/src/ninja/WindNinjaInputs.cpp @@ -156,6 +156,8 @@ WindNinjaInputs::WindNinjaInputs() foamVelocityGrid = -1.0; foamAngleGrid = -1.0; writeTurbulence = false; + colMax_colHeightAGL = 300.0; // default value of 300 m + colMax_colHeightAGL_units = lengthUnits::meters; #endif outputPointsFilename = "!set"; @@ -267,6 +269,8 @@ WindNinjaInputs &WindNinjaInputs::operator=(const WindNinjaInputs &rhs) foamVelocityGrid = rhs.foamVelocityGrid; foamAngleGrid = rhs.foamAngleGrid; writeTurbulence = rhs.writeTurbulence; + colMax_colHeightAGL = rhs.colMax_colHeightAGL; + colMax_colHeightAGL_units = rhs.colMax_colHeightAGL_units; #endif outputPointsFilename = rhs.outputPointsFilename; inputPointsFilename = rhs.inputPointsFilename; diff --git a/src/ninja/WindNinjaInputs.h b/src/ninja/WindNinjaInputs.h index fbc3751b..e340ad6d 100644 --- a/src/ninja/WindNinjaInputs.h +++ b/src/ninja/WindNinjaInputs.h @@ -304,6 +304,8 @@ struct WindNinjaInputs AsciiGrid foamVelocityGrid; //output velocity grid from ninjafoam AsciiGrid foamAngleGrid; //output angle grid from ninjafoam bool writeTurbulence; + double colMax_colHeightAGL; + lengthUnits::eLengthUnits colMax_colHeightAGL_units; #endif }; diff --git a/src/ninja/ascii_grid.cpp b/src/ninja/ascii_grid.cpp index 7aba5566..f2043b9d 100644 --- a/src/ninja/ascii_grid.cpp +++ b/src/ninja/ascii_grid.cpp @@ -1611,7 +1611,7 @@ GDALDatasetH AsciiGrid::ascii2GDAL() return hDS; } - + /** * Create a png for an ascii grid. @@ -1629,14 +1629,28 @@ void AsciiGrid::ascii2png(std::string outFilename, std::string legendTitle, std::string legendUnits, std::string scalarLegendFilename, - bool writeLegend) -{ + bool writeLegend, bool keepTiff) +{ + // png driver doesn't support create(), only createCopy() + // so create the image as a tif first, then use createCopy() to create the png from that tif + + std::string base_outFilename = outFilename; + int pos = outFilename.find_last_of(".png"); + //std::cout << "outFilename = " << outFilename.c_str() << std::endl; + //std::cout << "pos = " << pos << std::endl; + if(pos != -1) + base_outFilename = outFilename.substr(0, pos - 4 + 1); // .png is 4 letters back, + 1 to go from digit Id to a count + //std::cout << "base_outFilename = " << base_outFilename.c_str() << std::endl; + std::string tiff_utm_fileout = base_outFilename + "_utm.tif"; + std::string tiff_latlon_fileout = base_outFilename + "_latlon.tif"; + std::string rawTiff_utm_fileout = base_outFilename + "_raw_utm.tif"; + std::string rawTiff_latlon_fileout = base_outFilename + ".tif"; + GDALDataset *poDS; GDALDriver *tiffDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); char** papszOptions = NULL; - std::string tempFileout = "temp_fileout"; - poDS = tiffDriver->Create(tempFileout.c_str(), get_nCols(), get_nRows(), 1, + poDS = tiffDriver->Create(tiff_utm_fileout.c_str(), get_nCols(), get_nRows(), 1, GDT_Byte, papszOptions); double adfGeoTransform[6] = {get_xllCorner(), get_cellSize(), 0, @@ -1656,33 +1670,27 @@ void AsciiGrid::ascii2png(std::string outFilename, /* -------------------------------------------------------------------- */ AsciiGridscaledDataGrid(*this); - double translation = 0; - double scalingFactor = 1; - // convert nodata values to 0 (0 is transparent channel in color table) - for(int i=0;i(); //if want to show *real* 0 values in image - } - if(scaledDataGrid(i,j) == scaledDataGrid.get_noDataValue()) - { - scaledDataGrid(i,j) = 0; - } - } - }//scaledDataGrid.write_Grid("scaled_datagrid", 2); - // need min value (without 0s) later to make legend + //// set a patch of no data vals for debugging + //for(int i=0;i::max(); for(int i=0; i::ascii2png(std::string outFilename, } } - if(raw_minValue<0) - { - translation = fabs(raw_minValue)+1.1; //+1.1 since 0 = transparent in color table - } - else if(raw_minValue<2 && raw_minValue!=0) - { - translation = raw_minValue+1.1; - } - else if(raw_minValue!=0) //since 0 is transparent channel - { - translation = -raw_minValue-1.1; - } - if(raw_maxValue>255 || raw_maxValue-raw_minValue<5) - { - scalingFactor = 254/(raw_maxValue+translation); - } + // didn't see much difference when varying the range, the important thing is that scaled values are between 0 and 255 + // however, brk2 and brk4 divide may divide more evenly int wise into numbers/range divisible by 2 and 5. + // oops, need to reserve idx = 0 for the transparent channel in the color table, looks like need to manually convert noData vals to 0 + // if really want to show *real* 0 values in image, the noData vals, I guess for debugging?, set the noData vals to minIdx instead of to 0 + int idxRangeMin = 1; + int idxRangeMax = 255; + // numVals = idxRangeMax - idxRangeMin + 1, numBins = numVals - 1, so numBins = idxRangeMax - idxRangeMin + int numBins = idxRangeMax - idxRangeMin; + double binWidth = (raw_maxValue - raw_minValue)/double(numBins); for(int i=0;i raw_maxValue ){ + scaledDataGrid(i,j) = idxRangeMax; + } else { + // int( val + 0.5 ) here is to make int() behave like round() + int binIdx = int( (scaledDataGrid(i,j) - raw_minValue)/binWidth + 0.5 ) + idxRangeMin; + scaledDataGrid(i,j) = binIdx; } } } //scaledDataGrid.write_Grid("scaled_datagrid_again", 2); @@ -1773,13 +1778,13 @@ void AsciiGrid::ascii2png(std::string outFilename, } } + // need min value without 0s aka without no data vals, to make legend double _minValue = std::numeric_limits::max(); for(int i=nYSize-1;i>=0;i--) { for (int j=0;j::ascii2png(std::string outFilename, brk1 = _minValue; brk2 = 0.2*(_maxValue-_minValue)+_minValue; brk4 = _maxValue; - brk3 = (brk4+brk2)/2; + brk3 = (brk4+brk2)/2.0; poCT->SetColorEntry(brk0, &white); poCT->SetColorEntry(brk1, &blue); @@ -1814,7 +1819,7 @@ void AsciiGrid::ascii2png(std::string outFilename, const GDALColorEntry psEndColor3 = red; GDALCreateColorRamp(poCT, nStartIndex, &psStartColor3, nEndIndex, &psEndColor3); - //poBand->SetColorTable(poCT); + poBand->SetColorTable(poCT); /* -------------------------------------------------------------------- */ /* Create the optional legend */ @@ -1824,34 +1829,32 @@ void AsciiGrid::ascii2png(std::string outFilename, { //make bitmap int legendWidth = 180; - int legendHeight = int(legendWidth / 0.75); - BMP legend; + int legendHeight = int(legendWidth / 0.75); // increase by a factor of 4/3, rounded down. For legendWidth of 180, this comes out to be 240 + BMP legend; std::string legendStrings[6]; ostringstream os; - double maxxx = get_maxValue(); - double minnn = raw_minValue; double _brk0 = 0; double _brk1 = raw_minValue; - double _brk2 = 0.25*(raw_maxValue-raw_minValue)+raw_minValue; + double _brk2 = 0.20*(raw_maxValue-raw_minValue)+raw_minValue; double _brk4 = raw_maxValue; - double _brk3 = (_brk4+_brk2)/2; + double _brk3 = (_brk4+_brk2)/2.0; - for(int i = 0;i < 5;i++) + for(int labelIdx = 0; labelIdx < 5; labelIdx++) { os << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2); - if(i == 0) - os << maxxx; - else if (i==1) + if(labelIdx == 0) + os << (double)_brk4; + else if(labelIdx == 1) os << (double)_brk3; - else if(i == 2) + else if(labelIdx == 2) os << (double)_brk2; - else if(i == 3) - os << minnn; - else if(i == 4) - os << "0.00"; - legendStrings[i] = os.str(); + else if(labelIdx == 3) + os << (double)_brk1; + else if(labelIdx == 4) + os << "0.00"; // not currently used, but just in case leave it as a stored value + legendStrings[labelIdx] = os.str(); os.str(""); } @@ -1871,11 +1874,11 @@ void AsciiGrid::ascii2png(std::string outFilename, } //for white text - RGBApixel w; - w.Red = 255; - w.Green = 255; - w.Blue = 255; - w.Alpha = 0; + RGBApixel white; + white.Red = 255; + white.Green = 255; + white.Blue = 255; + white.Alpha = 0; RGBApixel colors[10]; @@ -1930,72 +1933,69 @@ void AsciiGrid::ascii2png(std::string outFilename, colors[9].Blue = 255; colors[9].Alpha = 0; - int arrowLength = 30; //pixels; - int textHeight = 12; - int titleTextHeight = int(1.2 * textHeight); - int titleX, titleY; + int cbarWidth = 40; // pixels; + int textHeight = 12; // cbarLabelTextHeight expected/desired value, DrawLine() seems to preserve the value, but DrawArc() in PrintString() seems to randomly pad up to +1 above and below, in addition to always padding +1 above and below to the value + int titleTextHeight = int(1.2 * textHeight); // increase by a factor, rounded down. For textHeight of 12, this comes out to be int(14.4) = 14. Also is expected/desired value, but the exact value isn't as sensitive as it is for textHeight + + double leftMarginPad = 0.05; // percent of total image width, the empty space to the left of the title box and cbar box regions + double topMarginPad = 0.05; // percent of total image height, the empty space to the top of the title box + + int titleX = legendWidth * leftMarginPad; // top left x pixel position for the title box. Note the int rounds it down. For a legendWidth of 180 and a leftMarginPad of 0.05, this comes out to be 9 + int titleY = legendHeight * topMarginPad; // top left y pixel position for the title box. Note the int rounds it down. For a legendHeight of 240 and a topMarginPad of 0.05, this comes out to be 12 + PrintString(legend, legendTitle.c_str() , titleX, titleY, titleTextHeight, white); + int unitsX = titleX+5; // set the unit string X position within the title box. titleX + 5 means to the right 5 pixels from the title X position, so nudging it a bit to the right + int unitsY = titleY+30; // set the unit string Y position within the title box. titleY + 30 means 30 pixels below the title Y position. For titleTextHeight of 14, this seems excessive, you would think that this would lead to a really huge gap between the two strings, but it actually results in barely any spacing. Apparently this is because PrintString() adds a crap ton of extra pixels for the letter p, it also adds a few pixels in general when calling DrawArc() for some of the other letters + PrintString(legend, legendUnits.c_str() , unitsX, unitsY, titleTextHeight, white); + + // start to end pixel positions at which to the given cbar color line int x1, x2; - double x; int y1, y2; - double y; + // top left pixel position at which to draw the cbar label text int textX; int textY; - x = 0.05; - y = 0.30; + double cbarBoxXstart_percent = leftMarginPad + 5.0/60.0; // percent of total image width, the 5.0/60.0 = 0.08333333333 represents adding an empty space to the left of the cbar region in addition to the empty space to the left of the title region. For x of 0.05 and the additional padding of 5.0/60.0 = 0.08333333333, this comes out to be 0.05 + 5.0/60.0 = 0.13333333333. Probably could have used 0.084, let int(legendWidth*cbarBoxXstart_percent) round down to 24, but using 0.05 + 5.0/60.0 specifically divides evenly into the legendWidth of 180 to get 24 + double cbarBoxYstart_percent = 0.297; // percent of total image height, care to choose a value for this that gives the title box enough space, with a little bit of padding between the title box and the cbar box region. The original expected value was 0.3 but this gave just a hint too much space, so it was adjusted back a bit + int cbarBoxXstart = legendWidth * cbarBoxXstart_percent; // top left x pixel position for the cbar box region. Note the int rounds it down. For a legendWidth of 180 and a cbarBoxXstart_percent of 5.0/60.0 = 0.08333333333, this comes out to be 24 + int cbarBoxYstart = legendHeight * cbarBoxYstart_percent; // top left y pixel position for the cbar box region. Note the int rounds it down. For a a legendHeight of 240 and a cbarBoxYstart_percent of 0.297, this comes out to be int(71.28) = 71 - titleX = x * legendWidth; - titleY = (y / 3) * legendHeight; + int cbarLabelXpadding = 15; // number of pixels of empty space to pad between the cbar and the cbar labels - PrintString(legend, legendTitle.c_str() , titleX, titleY-10, titleTextHeight, w); - PrintString(legend, legendUnits.c_str() , titleX+5, titleY+15, titleTextHeight, w); + // now the fun part, turns out that while DrawLine() does exact numbers of pixels, DrawArc() in PrintString() seems to randomly pad up to +1 above and below, in addition to always padding +1 above and below to the value. So you would think that a value of 1 above and 1 below would be good, but the gaps came out looking uneven. Turned out that a value of 1 above and 2 below seemed to look the best, with the gaps coming out looking even at least for this case. Note that a value of 2 above and 2 below looked similar to a value of 1 above and 1 below where the gap came out looking uneven, but with a bit more gap. Note that a value of 2 above and 3 below looked good and even just like a value of 1 above and 2 below, but with even more gap. In addition, using gaps 2 above and below or greater starts to add too many color lines overall, running out of space for the cbar + int upperLabelPad = 1; + int lowerLabelPad = 2; - x1 = int(legendWidth * x); - x2 = x1 + arrowLength; - y1 = int(legendHeight * y); - y2 = y1; - int k = -1; + int nColorsToUse = 10; // for 4 cbar colors - for(int i=0;i < 10;i++) + int labelIdx = 0; + int yPos = cbarBoxYstart; // yPos is the current pixel y position + int yStart = yPos; // yStart is the first pixel y position of the given cbar colorbox + for(int colorIdx = 0; colorIdx < nColorsToUse; colorIdx++) { - for(int j=0; j<20; j++) + yStart = yPos; // store the first pixel y position of the given cbar colorbox + for(int colorLineIdx = 0; colorLineIdx < (upperLabelPad+textHeight+lowerLabelPad); colorLineIdx++) { - x1 = int(legendWidth * (x+0.1)); - x2 = x1 + arrowLength; - y1 = int(legendHeight * (y-0.03)); + x1 = cbarBoxXstart; + x2 = x1 + cbarWidth; + y1 = yPos; y2 = y1; - DrawLine(legend, x1, y1, x2, y2, colors[i]); - y+=0.003; + DrawLine(legend, x1, y1, x2, y2, colors[colorIdx]); + yPos+=1; } - if(i==0 || i==3 || i==6 || i==9) + + // only add labels to the main colors + if(colorIdx == 0 || colorIdx == 3 || colorIdx == 6 || colorIdx == 9) { - k+=1; - textX = x2 + 15; - textY = y2 - int(textHeight * 0.5) - 8; - PrintString(legend, legendStrings[k].c_str(), textX, textY, textHeight, w); + textX = x2 + cbarLabelXpadding; + textY = yStart + upperLabelPad; + PrintString(legend, legendStrings[labelIdx].c_str(), textX, textY, textHeight, white); + labelIdx+=1; } } legend.WriteToFile(scalarLegendFilename.c_str()); } - /* -------------------------------------------------------------------- */ - /* close and reopen with GDALOpenShared and set color table */ - /* -------------------------------------------------------------------- */ - - AsciiGrid poDS_grid(poDS, 1); - poDS_grid.write_Grid("poDS_grid", 0); - - GDALDataset *srcDS; - //reopen with GDALOpenShared() for GDALAutoCreateWarpedVRT() - srcDS = (GDALDataset*)GDALOpenShared( "poDS_grid", GA_ReadOnly ); - if( srcDS == NULL ) { - CPLDebug( "ascii_grid::ascii2png()", "cannot open poDS_grid"); - } - - GDALRasterBand *srcBand = srcDS->GetRasterBand(1); - srcBand->SetColorTable(poCT); - /* -------------------------------------------------------------------- */ /* Warp the image */ /* -------------------------------------------------------------------- */ @@ -2013,10 +2013,20 @@ void AsciiGrid::ascii2png(std::string outFilename, GDALDataset *wrpDS; - wrpDS = (GDALDataset*)GDALAutoCreateWarpedVRT(srcDS, pszSRS_WKT, pszDST_WKT, + wrpDS = (GDALDataset*)GDALAutoCreateWarpedVRT(poDS, pszSRS_WKT, pszDST_WKT, GRA_NearestNeighbour, 0.0, NULL); + /* -------------------------------------------------------------------- */ + /* Write the warped tiff */ + /* -------------------------------------------------------------------- */ + + GDALDataset *poDstDS_tiff; + + CPLPushErrorHandler(CPLQuietErrorHandler); //silence TIFF dirver data type error + poDstDS_tiff = tiffDriver->CreateCopy(tiff_latlon_fileout.c_str(), wrpDS, FALSE, NULL, NULL, NULL); + CPLPopErrorHandler(); + /* -------------------------------------------------------------------- */ /* Write the png */ /* -------------------------------------------------------------------- */ @@ -2054,15 +2064,101 @@ void AsciiGrid::ascii2png(std::string outFilename, GDALDestroyColorTable((GDALColorTableH) poCT); GDALClose((GDALDatasetH) poDS); - GDALClose((GDALDatasetH) srcDS); if( poDstDS != NULL){ GDALClose((GDALDatasetH) poDstDS); } + if( poDstDS_tiff != NULL){ + GDALClose((GDALDatasetH) poDstDS_tiff); + } GDALClose((GDALDatasetH) wrpDS); - VSIUnlink("poDS_grid"); - VSIUnlink("poDS_grid.aux.xml"); - VSIUnlink("temp_fileout"); + VSIUnlink(tiff_utm_fileout.c_str()); + VSIUnlink(tiff_latlon_fileout.c_str()); + + + /* -------------------------------------------------------------------- */ + /* temporary, create additional raw tiff file */ + /* -------------------------------------------------------------------- */ + + GDALDataset *raw_poDS; + GDALDriver *raw_tiffDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + char** raw_papszOptions = NULL; + + raw_poDS = tiffDriver->Create(rawTiff_utm_fileout.c_str(), get_nCols(), get_nRows(), 1, + GDT_Float64, papszOptions); + + double raw_adfGeoTransform[6] = {get_xllCorner(), get_cellSize(), 0, + get_yllCorner()+(get_nRows()*get_cellSize()), + 0, -(get_cellSize())}; + raw_poDS->SetGeoTransform(raw_adfGeoTransform); + int raw_nXSize = raw_poDS->GetRasterXSize(); + int raw_nYSize = raw_poDS->GetRasterYSize(); + + GDALRasterBand *raw_poBand = raw_poDS->GetRasterBand(1); + + AsciiGridrawDataGrid(*this); + + double *raw_padfScanline; + raw_padfScanline = new double[raw_nXSize]; + + for(int i=raw_nYSize-1;i>=0;i--) + { + for(int j=0;jRasterIO(GF_Write, 0, i, raw_nXSize, 1, raw_padfScanline, raw_nXSize, 1, GDT_Float64, 0, 0)); + } + } + + /* -------------------------------------------------------------------- */ + /* Warp the image */ + /* -------------------------------------------------------------------- */ + + OGRSpatialReference raw_oSRS; + char *raw_pszSRS_WKT = NULL; + + const char* raw_prj2 = (const char*)prjString.c_str(); + raw_oSRS.importFromWkt((char **)&raw_prj2); + raw_oSRS.exportToWkt(&raw_pszSRS_WKT); + + char *raw_pszDST_WKT = NULL; + raw_oSRS.importFromEPSG(4326); + raw_oSRS.exportToWkt(&raw_pszDST_WKT); + + GDALDataset *raw_wrpDS; + + raw_wrpDS = (GDALDataset*)GDALAutoCreateWarpedVRT(raw_poDS, raw_pszSRS_WKT, raw_pszDST_WKT, + GRA_NearestNeighbour, + 0.0, NULL); + + /* -------------------------------------------------------------------- */ + /* Write the warped tiff */ + /* -------------------------------------------------------------------- */ + + GDALDataset *poDstDS_rawTiff; + + CPLPushErrorHandler(CPLQuietErrorHandler); //silence TIFF driver data type error + poDstDS_rawTiff = raw_tiffDriver->CreateCopy(rawTiff_latlon_fileout.c_str(), raw_wrpDS, FALSE, NULL, NULL, NULL); + CPLPopErrorHandler(); + + /* -------------------------------------------------------------------- */ + /* clean up */ + /* -------------------------------------------------------------------- */ + + CPLFree(raw_pszSRS_WKT); + CPLFree(raw_pszDST_WKT); + + delete [] raw_padfScanline; + + GDALClose((GDALDatasetH) raw_poDS); + if( poDstDS_rawTiff != NULL){ + GDALClose((GDALDatasetH) poDstDS_rawTiff); + } + GDALClose((GDALDatasetH) raw_wrpDS); + + VSIUnlink(rawTiff_utm_fileout.c_str()); + if( keepTiff == false ) + VSIUnlink(rawTiff_latlon_fileout.c_str()); } diff --git a/src/ninja/ascii_grid.h b/src/ninja/ascii_grid.h index f954dd6d..ff64d534 100644 --- a/src/ninja/ascii_grid.h +++ b/src/ninja/ascii_grid.h @@ -207,7 +207,7 @@ class AsciiGrid std::string legendTitle, std::string legendUnits, std::string scalarLegendFilename, - bool writeLegend); + bool writeLegend, bool keepTiff); void exportToTiff( std::string outFilename, tiffType type = tiffGray ); diff --git a/src/ninja/cli.cpp b/src/ninja/cli.cpp index a73b3ded..12c78ed8 100644 --- a/src/ninja/cli.cpp +++ b/src/ninja/cli.cpp @@ -923,6 +923,7 @@ int windNinjaCLI(int argc, char* argv[]) conflicting_options(vm, "momentum_flag", "input_points_file"); conflicting_options(vm, "momentum_flag", "write_vtk_output"); option_dependency(vm, "turbulence_output_flag", "momentum_flag"); + option_dependency(vm, "turbulence_output_flag", "write_goog_output"); #ifdef FRICTION_VELOCITY conflicting_options(vm, "momentum_flag", "compute_friction_velocity"); #endif diff --git a/src/ninja/gdal_util.cpp b/src/ninja/gdal_util.cpp index 7899f5d6..c4160e13 100644 --- a/src/ninja/gdal_util.cpp +++ b/src/ninja/gdal_util.cpp @@ -94,6 +94,7 @@ bool GDALGetCenter( GDALDataset *poDS, double *longitude, double *latitude ) #ifdef GDAL_COMPUTE_VERSION #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) OSRSetAxisMappingStrategy(hTargetSRS, OAMS_TRADITIONAL_GIS_ORDER); + OSRSetAxisMappingStrategy(hSrcSRS, OAMS_TRADITIONAL_GIS_ORDER); #endif /* GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) */ #endif /* GDAL_COMPUTE_VERSION */ @@ -896,10 +897,7 @@ bool GDALWarpToUtm (const char* filename, GDALDatasetH& hSrcDS, GDALDatasetH& hD OGRSpatialReference oDstSRS; - double lat, lon; - GDALGetCenter( (GDALDataset*)hSrcDS, &lon, &lat); - int utmZone = gdalGetUtmZone(lat, lon); - int nUtmZone = GetUTMZoneInEPSG(lat, lon); + int nUtmZone = GDALGetUtmZone( (GDALDataset*)hSrcDS ); oDstSRS.importFromEPSG(nUtmZone); oDstSRS.exportToWkt((char**)&pszDstWKT); diff --git a/src/ninja/ninja.cpp b/src/ninja/ninja.cpp index 2ea31bf3..4c285d22 100644 --- a/src/ninja/ninja.cpp +++ b/src/ninja/ninja.cpp @@ -106,6 +106,7 @@ ninja::ninja(const ninja &rhs) , CloudGrid(rhs.CloudGrid) #ifdef NINJAFOAM , TurbulenceGrid(rhs.TurbulenceGrid) +, colMaxGrid(rhs.colMaxGrid) #endif , outputSpeedArray(rhs.outputSpeedArray) , outputDirectionArray(rhs.outputDirectionArray) @@ -187,6 +188,7 @@ ninja &ninja::operator=(const ninja &rhs) CloudGrid = rhs.CloudGrid; #ifdef NINJAFOAM TurbulenceGrid = rhs.TurbulenceGrid; + colMaxGrid = rhs.colMaxGrid; #endif outputSpeedArray=rhs.outputSpeedArray; outputDirectionArray = rhs.outputDirectionArray; @@ -592,6 +594,7 @@ if(input.frictionVelocityFlag == 1){ CloudGrid.deallocate(); #ifdef NINJAFOAM TurbulenceGrid.deallocate(); + colMaxGrid.deallocate(); #endif #ifdef FRICTION_VELOCITY if(input.frictionVelocityFlag == 1){ @@ -2203,8 +2206,15 @@ void ninja::prepareOutput() #ifdef NINJAFOAM if(input.writeTurbulence) { - TurbulenceGrid.set_headerData(input.dem.get_nCols(),input.dem.get_nRows(), input.dem.get_xllCorner(), - input.dem.get_yllCorner(), input.dem.get_cellSize(), input.dem.get_noDataValue(), 0, input.dem.prjString); + // um, the only way it ever gets to this point is for a diurnal run on a ninjafoam run, represented by if input.writeTurbulence is set for a mass solver run + // BUT, set_headerData() is NOT the appropriate way to deal with the data here, the data already comes in filled. Resampling to the new input resolution + // would be more appropriate, but since the diurnal run defines mesh and dem with the same resolution as the ninjafoam run, they are on the same grid + // this is evidenced by not seeing any calls to resample to a different resolution in the foam initialization classes + // leave this here as a reminder just in case though, probably should revisit when trying to do a final double check ALL grid consistencies in ALL places + //TurbulenceGrid.set_headerData(input.dem.get_nCols(),input.dem.get_nRows(), input.dem.get_xllCorner(), + // input.dem.get_yllCorner(), input.dem.get_cellSize(), input.dem.get_noDataValue(), 0, input.dem.prjString); + //colMaxGrid.set_headerData(input.dem.get_nCols(),input.dem.get_nRows(), input.dem.get_xllCorner(), + // input.dem.get_yllCorner(), input.dem.get_cellSize(), input.dem.get_noDataValue(), 0, input.dem.prjString); } #endif @@ -2232,6 +2242,7 @@ void ninja::prepareOutput() AngleGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); #ifdef NINJAFOAM TurbulenceGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); + colMaxGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); #endif //Clip cloud cover grid if it's a wxModel intitialization (since it's gridded) @@ -2246,6 +2257,7 @@ void ninja::prepareOutput() if(input.writeTurbulence) { velocityUnits::fromBaseUnits(TurbulenceGrid, input.outputSpeedUnits); + velocityUnits::fromBaseUnits(colMaxGrid, input.outputSpeedUnits); } #endif @@ -2702,7 +2714,7 @@ void ninja::computeDustEmissions() dust.ComputePM10(UstarGrid, DustGrid); - //DustGrid.ascii2png("dust_out_shear.png", "pm10", "ug/m3", "legend", false); + //DustGrid.ascii2png("dust_out_shear.png", "pm10", "ug/m3", "legend", false, false); } #endif //EMISISONS @@ -2996,13 +3008,14 @@ void ninja::writeOutputFiles() { AsciiGrid *velTempGrid, *angTempGrid; #ifdef NINJAFOAM - AsciiGrid *turbTempGrid; + AsciiGrid *turbTempGrid, *colMaxTempGrid; #endif velTempGrid=NULL; angTempGrid=NULL; #ifdef NINJAFOAM turbTempGrid=NULL; + colMaxTempGrid=NULL; #endif KmlVector ninjaKmlFiles; @@ -3013,11 +3026,18 @@ void ninja::writeOutputFiles() #ifdef NINJAFOAM if(input.writeTurbulence) { - turbTempGrid = new AsciiGrid (TurbulenceGrid.resample_Grid(input.kmzResolution, + //turbTempGrid = new AsciiGrid (TurbulenceGrid.resample_Grid(input.kmzResolution, + // AsciiGrid::order0)); + // + //ninjaKmlFiles.setTurbulenceFlag("true"); + //ninjaKmlFiles.setTurbulenceGrid(*turbTempGrid, input.outputSpeedUnits); + + + colMaxTempGrid = new AsciiGrid (colMaxGrid.resample_Grid(input.kmzResolution, AsciiGrid::order0)); - ninjaKmlFiles.setTurbulenceFlag("true"); - ninjaKmlFiles.setTurbulenceGrid(*turbTempGrid, input.outputSpeedUnits); + ninjaKmlFiles.setColMaxFlag("true"); + ninjaKmlFiles.setColMaxGrid(*colMaxTempGrid, input.outputSpeedUnits, input.colMax_colHeightAGL, input.colMax_colHeightAGL_units); } #endif //NINJAFOAM @@ -3097,6 +3117,11 @@ void ninja::writeOutputFiles() delete turbTempGrid; turbTempGrid=NULL; } + if(colMaxTempGrid) + { + delete colMaxTempGrid; + colMaxTempGrid=NULL; + } #endif } }catch (exception& e) @@ -3669,6 +3694,12 @@ void ninja::set_writeTurbulenceFlag(bool flag) { input.writeTurbulence = flag; } + +void ninja::set_colMaxSampleHeightAGL( double colMaxSampleHeightAGL, lengthUnits::eLengthUnits units ) +{ + input.colMax_colHeightAGL = colMaxSampleHeightAGL; + input.colMax_colHeightAGL_units = units; +} #endif void ninja::set_speedFile(std::string speedFile, velocityUnits::eVelocityUnits units) diff --git a/src/ninja/ninja.h b/src/ninja/ninja.h index 0a07db93..31d5fbf2 100644 --- a/src/ninja/ninja.h +++ b/src/ninja/ninja.h @@ -138,6 +138,7 @@ class ninja AsciiGridCloudGrid; #ifdef NINJAFOAM AsciiGridTurbulenceGrid; //this needs to be a member of ninja since we need to write this for ninjafoam runs with diurnal + AsciiGrid colMaxGrid; // same as for TurbulenceGrid #endif wn_3dScalarField alphaVfield; //store spatially varying alphaV variable @@ -285,6 +286,7 @@ class ninja void set_foamVelocityGrid(AsciiGrid velocityGrid); void set_foamAngleGrid(AsciiGrid angleGrid); void set_writeTurbulenceFlag(bool flag); + void set_colMaxSampleHeightAGL( double colMaxSampleHeightAGL, lengthUnits::eLengthUnits units ); #endif void set_speedFile(std::string speedFile, velocityUnits::eVelocityUnits units); diff --git a/src/ninja/ninja_init.cpp b/src/ninja/ninja_init.cpp index 501b64f5..41c24ba4 100644 --- a/src/ninja/ninja_init.cpp +++ b/src/ninja/ninja_init.cpp @@ -27,53 +27,116 @@ *****************************************************************************/ #include "ninja_init.h" +#include +#include +#include +#include +#include "cpl_http.h" +#include +#include +#include +#include "ninja_version.h" boost::local_time::tz_database globalTimeZoneDB; /* -** Query the windninja.org server and ask for the most up to date version. The -** return value is a set of key value pairs stored in a GDAL string list. +** Query the ninjastorm.firelab.org/sqlitetest/messages.txt and ask for the most up to date version. ** There are three current values: ** ** VERSION -> a semantic version string, comparable with strcmp() ** MESSAGE -> One or more messages to display to the user ** ABORT -> There is a fundamental problem with windninja, and it should call ** abort() after displaying a message. -** -** The response is currently a semi-colon delimeted list of key value pairs as -** in: -** -** VERSION=3.4.0;MESSAGE=some message for the user, may not have -** semi-colons;ABORT:TRUE -** -** The returned string list must be freed by the caller using CSLDestroy(). */ -char ** NinjaCheckVersion(void) { - CPLHTTPResult *poResult; - char **papszTokens = NULL; - char *pszResp = NULL; - CPLPushErrorHandler(CPLQuietErrorHandler); - poResult = CPLHTTPFetch("http://windninja.org/version/", NULL); - CPLPopErrorHandler(); - if (!poResult || poResult->nStatus != 0 || poResult->nDataLen == 0) { + +bool NinjaCheckVersions(char * mostrecentversion, char * localversion) { + char comparemostrecentversion[256]; + char comparelocalversion[256]; + int mostrecentversionIndex = 0; + int localversionIndex = 0; + while (*mostrecentversion) { + if (*mostrecentversion != '.') { + comparemostrecentversion[mostrecentversionIndex++] = *mostrecentversion; + } + mostrecentversion++; + } + comparemostrecentversion[mostrecentversionIndex] = '\0'; + + while (*localversion) { + if (*localversion != '.') { + comparelocalversion[localversionIndex++] = *localversion; + } + localversion++; + } + + comparelocalversion[localversionIndex] = '\0'; + return strcmp(comparemostrecentversion, comparelocalversion) == 0; + +} + +char * NinjaQueryServerMessages(bool checkAbort) { + const char* url = "https://ninjastorm.firelab.org/sqlitetest/messages.txt"; + CPLHTTPResult *poResult = NULL; + try{ + CPLHTTPResult *poResult = CPLHTTPFetch(url, NULL); + + if (poResult != NULL) { + const char* pszTextContent = reinterpret_cast(poResult->pabyData); + std::vector messages; + std::istringstream iss(pszTextContent); + std::string message; + + // Read all lines into the vector + while (std::getline(iss, message)) { + messages.push_back(message); + } + + // Process all lines except the last two + std::ostringstream oss; + if (checkAbort) { + for (size_t i = 0; i < messages.size(); ++i) { + if (i == messages.size()-1) { // check final line + oss << messages[i] << "\n"; + } + } + } + else { + bool versionisuptodate = NinjaCheckVersions(const_cast(messages[1].c_str()), const_cast(NINJA_VERSION_STRING)); + if (!versionisuptodate) { + oss << messages[0] << "\n"; + oss << "You are using an outdated WindNinja version, please update to version: " << messages[1] << "\n" << "\n"; + } + + if (messages[4].empty() == false) { + for (size_t i = 3; i < messages.size() - 2; ++i) { + if (!messages[i].empty()) { + oss << messages[i] << "\n"; + } + } + } + if (messages[4].empty() && versionisuptodate) { + return NULL; + } + } + + std::string resultingmessage = oss.str(); + char* returnString = new char[resultingmessage.length() + 1]; + std::strcpy(returnString, resultingmessage.c_str()); + + CPLHTTPDestroyResult(poResult); + return returnString; + } + } + catch (std::exception& e) { + std::cout << "can't fetch" << std::endl; + } return NULL; - } - pszResp = (char *)malloc(poResult->nDataLen + 1); - if (pszResp == 0) { - return NULL; - } - /* - ** Copy the message body into a null terminated string for - ** CSLTokenizeString() - */ - memcpy(pszResp, poResult->pabyData, poResult->nDataLen); - pszResp[poResult->nDataLen] = '\0'; - papszTokens = CSLTokenizeString2((const char *)pszResp, ";", 0); - free(pszResp); - CPLHTTPDestroyResult( poResult ); - return papszTokens; + } + + + void NinjaCheckThreddsData( void *rc ) { int *r; @@ -120,7 +183,6 @@ int NinjaInitialize(const char *pszGdalData, const char *pszWindNinjaData) GDALAllRegister(); OGRRegisterAll(); - CPLSetConfigOption( "GDAL_HTTP_UNSAFESSL", "YES"); if(!CPLCheckForFile(CPLStrdup(CPLFormFilename(CPLStrdup(pszGdalData), "gdalicon.png", NULL)), NULL)) { @@ -147,8 +209,10 @@ int NinjaInitialize(const char *pszGdalData, const char *pszWindNinjaData) /* ** Initialize global singletons and environments. */ -int NinjaInitialize() -{ + +int NinjaInitialize(const char* typeofrun) { + + GDALAllRegister(); OGRRegisterAll(); /* @@ -156,8 +220,7 @@ int NinjaInitialize() ** but that doesn't mean we are in trouble. */ CPLPushErrorHandler(CPLQuietErrorHandler); - int rc = 0; - + int rc = 0; /* ** Setting the CURL_CA_BUNDLE variable through GDAL doesn't seem to work, ** but could be investigated in the future. CURL_CA_BUNDLE can only be set in GDAL @@ -166,6 +229,46 @@ int NinjaInitialize() ** For now, just skip the SSL verification with GDAL_HTTP_UNSAFESSL. */ CPLSetConfigOption( "GDAL_HTTP_UNSAFESSL", "YES"); + + if (strcmp(typeofrun, "") != 0) { + + time_t now = time(0); + + // convert now to tm struct for UTC + tm *gmtm = gmtime(&now); + char* dt = asctime(gmtm); + std::string cpp_string(dt); + + + std::string url = "https://ninjastorm.firelab.org/sqlitetest/?time="; + cpp_string.erase(std::remove_if(cpp_string.begin(), cpp_string.end(), ::isspace), + cpp_string.end()); + + + std::string full = url + cpp_string + "&runtype=" + typeofrun; + + + const char *charStr = full.data(); + + CPLHTTPResult *poResult; + CPLSetConfigOption("GDAL_HTTP_UNSAFESSL", "YES"); + char **papszOptions = NULL; + + // Fetch the URL with custom headers + try { + poResult = CPLHTTPFetch(charStr, papszOptions); + } + catch (std::exception& e) { + std::cout << "can't fetch" << std::endl; + } + + if (poResult) { + CPLHTTPDestroyResult(poResult); + + } + } + + #ifdef WIN32 CPLDebug( "WINDNINJA", "Setting GDAL_DATA..." ); @@ -210,19 +313,24 @@ int NinjaInitialize() CPLFree( (void*)pszExecPath ); #endif /* defined(NINJAFOAM) && defined(FIRELAB_PACKAGE)*/ + #endif - /* - ** Set windninja data if it isn't set. - */ - if( !CSLTestBoolean( CPLGetConfigOption( "WINDNINJA_DATA", "FALSE" ) ) ) - { - std::string osDataPath; - osDataPath = FindDataPath( "tz_world.zip" ); - if( osDataPath != "" ) - { - CPLSetConfigOption( "WINDNINJA_DATA", CPLGetPath( osDataPath.c_str() ) ); - } + + + + +/* +** Set windninja data if it isn't set. +*/ +if (!CSLTestBoolean(CPLGetConfigOption("WINDNINJA_DATA", "FALSE"))) { + std::string osDataPath; + osDataPath = FindDataPath("tz_world.zip"); + if (osDataPath != "") { + CPLSetConfigOption("WINDNINJA_DATA", CPLGetPath(osDataPath.c_str())); } +} + + globalTimeZoneDB.load_from_file(FindDataPath("date_time_zonespec.csv")); CPLPopErrorHandler(); return 0; diff --git a/src/ninja/ninja_init.h b/src/ninja/ninja_init.h index f7cd6cac..71185ada 100644 --- a/src/ninja/ninja_init.h +++ b/src/ninja/ninja_init.h @@ -43,9 +43,10 @@ #include "boost/date_time/local_time/local_time.hpp" #endif -int NinjaInitialize(); +int NinjaInitialize(const char* typeofrun = ""); int NinjaInitialize(const char *pszGdalData, const char *pszWindNinjaData); -char ** NinjaCheckVersion(void); +char * NinjaQueryServerMessages(bool checkAbort); +bool NinjaCheckVersions(char * mostrecentversion, char * localversion); #endif /* NINJA_INIT_H_ */ diff --git a/src/ninja/ninjafoam.cpp b/src/ninja/ninjafoam.cpp index 709694fc..df54cd8b 100644 --- a/src/ninja/ninjafoam.cpp +++ b/src/ninja/ninjafoam.cpp @@ -73,6 +73,8 @@ NinjaFoam::NinjaFoam() : ninja() startStlConversion = 0.0; endStlConversion = 0.0; + writeMassMesh = false; + writeMassMeshVtk = false; } @@ -126,7 +128,11 @@ double NinjaFoam::get_meshResolution() bool NinjaFoam::simulate_wind() { + #ifdef WIN32 + foamVersion = "2.2.0"; + #else foamVersion = CPLGetConfigOption("WM_PROJECT_VERSION", ""); + #endif CPLDebug("NINJAFOAM", CPLSPrintf("foamVersion = \"%s\"",foamVersion.c_str())); if(CSLTestBoolean(CPLGetConfigOption("WRITE_TURBULENCE", "FALSE"))) @@ -138,6 +144,31 @@ bool NinjaFoam::simulate_wind() { writeMassMeshVtk = CPLGetConfigOption("WRITE_FOAM_MASSMESH_VTK", "FALSE"); } + + if(input.writeTurbulence == true) + { + std::string found_colHeightAGL_str = CPLGetConfigOption("COLMAX_HEIGHT_AGL", ""); + // if read it, but no value, don't want a 0.0 put in, leave it as default value + // Well at least it didn't break with a 0.0, just grabbed surf values, but still, want it as the default value unless they specify it with an actual value + if( found_colHeightAGL_str != "" ) + { + double found_colHeightAGL = atof(found_colHeightAGL_str.c_str()); + std::string found_colHeightAGL_units = CPLGetConfigOption("COLMAX_HEIGHT_AGL_UNITS", "m"); + if( found_colHeightAGL_units == "" ){ + found_colHeightAGL_units = "m"; // default value if not set + } + std::cout << "found CPL config option COLMAX_HEIGHT_AGL, setting colMax_colHeightAGL to " << found_colHeightAGL << " " << found_colHeightAGL_units << std::endl; + set_colMaxSampleHeightAGL( found_colHeightAGL, lengthUnits::getUnit(found_colHeightAGL_units) ); + } + CPLDebug("NINJAFOAM", "Writing turbulence colMax output, using colHeightAGL %f %s",input.colMax_colHeightAGL,lengthUnits::getString(input.colMax_colHeightAGL_units).c_str()); + } + + + if( writeMassMeshVtk == true || input.writeTurbulence == true ) + { + writeMassMesh = true; + } + #ifdef _OPENMP startTotal = omp_get_wtime(); @@ -355,6 +386,11 @@ bool NinjaFoam::simulate_wind() #ifdef _OPENMP endOutputSampling = omp_get_wtime(); #endif + + /*-------------------------------------------------------------------*/ + /* Generate and Sample mass mesh */ + /*-------------------------------------------------------------------*/ + GenerateAndSampleMassMesh(); /*----------------------------------------*/ /* write output files */ @@ -397,6 +433,12 @@ bool NinjaFoam::simulate_wind() AngleGrid.deallocate(); VelocityGrid.deallocate(); TurbulenceGrid.deallocate(); + colMaxGrid.deallocate(); + + massMesh_u.deallocate(); + massMesh_v.deallocate(); + massMesh_w.deallocate(); + massMesh_k.deallocate(); } if(input.diurnalWinds == true){ @@ -2143,143 +2185,6 @@ const char * NinjaFoam::GetGridFilename() return pszGridFilename; } -void NinjaFoam::SetOutputResolution() -{ - //Set output file resolutions now - if( input.kmzResolution <= 0.0 ) //if negative, use DEM resolution - input.kmzResolution = input.dem.get_cellSize(); - if( input.shpResolution <= 0.0 ) //if negative, use DEM resolution - input.shpResolution = input.dem.get_cellSize(); - if( input.velResolution <= 0.0 ) //if negative, use DEMresolution - input.velResolution = input.dem.get_cellSize(); - if( input.angResolution <= 0.0 ) //if negative, use DEM resolution - input.angResolution = input.dem.get_cellSize(); - if( input.pdfResolution <= 0.0 ) - input.pdfResolution = input.dem.get_cellSize(); -} - -void NinjaFoam::SetOutputFilenames() -{ - //Do file naming string stuff for all output files - std::string rootFile, rootName, timeAppend, wxModelTimeAppend, fileAppend, kmz_fileAppend, \ - shp_fileAppend, ascii_fileAppend, mesh_units, kmz_mesh_units, \ - shp_mesh_units, ascii_mesh_units, pdf_fileAppend, pdf_mesh_units; - - boost::local_time::local_time_facet* timeOutputFacet; - timeOutputFacet = new boost::local_time::local_time_facet(); - //NOTE: WEIRD ISSUE WITH THE ABOVE 2 LINES OF CODE! DO NOT CALL DELETE ON THIS BECAUSE THE LOCALE OBJECT BELOW DOES. - // THIS IS A "PROBLEM" IN THE STANDARD LIBRARY. SEE THESE WEB SITES FOR MORE INFO: - // https://collab.firelab.org/software/projects/windninja/wiki/KnownIssues - // http://rhubbarb.wordpress.com/2009/10/17/boost-datetime-locales-and-facets/#comment-203 - - std::ostringstream timestream; - timestream.imbue(std::locale(std::locale::classic(), timeOutputFacet)); - timeOutputFacet->format("%m-%d-%Y_%H%M_"); - - if( input.initializationMethod == WindNinjaInputs::wxModelInitializationFlag ){ - timestream << input.ninjaTime; - } - - std::string pathName; - std::string baseName(CPLGetBasename(input.dem.fileName.c_str())); - - if(input.customOutputPath == "!set"){ - //prepend directory paths to rootFile for wxModel run - if( input.initializationMethod == WindNinjaInputs::wxModelInitializationFlag ){ - pathName = CPLGetPath(input.forecastFilename.c_str()); - //if it's a .tar, write to directory containing the .tar file - if( strstr(pathName.c_str(), ".tar") ){ - pathName.erase( pathName.rfind("/") ); - } - }else{ - pathName = CPLGetPath(input.dem.fileName.c_str()); - } - } - else{ - pathName = input.customOutputPath; - } - - rootFile = CPLFormFilename(pathName.c_str(), baseName.c_str(), NULL); - - /* set the output path member variable */ - input.outputPath = pathName; - - timeAppend = timestream.str(); - - ostringstream wxModelTimestream; - boost::local_time::local_time_facet* wxModelOutputFacet; - wxModelOutputFacet = new boost::local_time::local_time_facet(); - wxModelTimestream.imbue(std::locale(std::locale::classic(), wxModelOutputFacet)); - wxModelOutputFacet->format("%m-%d-%Y_%H%M"); - if( input.initializationMethod == WindNinjaInputs::wxModelInitializationFlag) - { - wxModelTimestream << input.ninjaTime; - } - wxModelTimeAppend = wxModelTimestream.str(); - mesh_units = "m"; - kmz_mesh_units = lengthUnits::getString( input.kmzUnits ); - shp_mesh_units = lengthUnits::getString( input.shpUnits ); - ascii_mesh_units = lengthUnits::getString( input.velOutputFileDistanceUnits ); - pdf_mesh_units = lengthUnits::getString( input.pdfUnits ); - - ostringstream os, os_kmz, os_shp, os_ascii, os_pdf; - - if( input.initializationMethod == WindNinjaInputs::domainAverageInitializationFlag ){ - double tempSpeed = input.inputSpeed; - velocityUnits::fromBaseUnits(tempSpeed, input.inputSpeedUnits); - os << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); - os_kmz << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); - os_shp << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); - os_ascii << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); - os_pdf << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); - } - - double meshResolutionTemp = input.dem.get_cellSize(); - double kmzResolutionTemp = input.kmzResolution; - double shpResolutionTemp = input.shpResolution; - double velResolutionTemp = input.velResolution; - double pdfResolutionTemp = input.pdfResolution; - - lengthUnits::eLengthUnits meshResolutionUnits = lengthUnits::meters; - - lengthUnits::fromBaseUnits(meshResolutionTemp, meshResolutionUnits); - lengthUnits::fromBaseUnits(kmzResolutionTemp, input.kmzUnits); - lengthUnits::fromBaseUnits(shpResolutionTemp, input.shpUnits); - lengthUnits::fromBaseUnits(velResolutionTemp, input.velOutputFileDistanceUnits); - lengthUnits::fromBaseUnits(pdfResolutionTemp, input.pdfUnits); - - os << "_" << timeAppend << (long) (meshResolutionTemp+0.5) << mesh_units; - os_kmz << "_" << timeAppend << (long) (kmzResolutionTemp+0.5) << kmz_mesh_units; - os_shp << "_" << timeAppend << (long) (shpResolutionTemp+0.5) << shp_mesh_units; - os_ascii << "_" << timeAppend << (long) (velResolutionTemp+0.5) << ascii_mesh_units; - os_pdf << "_" << timeAppend << (long) (pdfResolutionTemp+0.5) << pdf_mesh_units; - - fileAppend = os.str(); - kmz_fileAppend = os_kmz.str(); - shp_fileAppend = os_shp.str(); - ascii_fileAppend = os_ascii.str(); - pdf_fileAppend = os_pdf.str(); - - input.kmlFile = rootFile + kmz_fileAppend + ".kml"; - input.kmzFile = rootFile + kmz_fileAppend + ".kmz"; - - input.shpFile = rootFile + shp_fileAppend + ".shp"; - input.dbfFile = rootFile + shp_fileAppend + ".dbf"; - - input.pdfFile = rootFile + pdf_fileAppend + ".pdf"; - - input.cldFile = rootFile + ascii_fileAppend + "_cld.asc"; - input.velFile = rootFile + ascii_fileAppend + "_vel.asc"; - input.angFile = rootFile + ascii_fileAppend + "_ang.asc"; - input.atmFile = rootFile + ascii_fileAppend + ".atm"; - - input.legFile = rootFile + kmz_fileAppend + ".bmp"; - if( input.ninjaTime.is_not_a_date_time() ) //date and time not set? - input.dateTimeLegFile = ""; - else - input.dateTimeLegFile = rootFile + kmz_fileAppend + ".date_time" + ".bmp"; -} - void NinjaFoam::SampleRawOutput() { /*-------------------------------------------------------------------*/ @@ -2370,359 +2275,66 @@ void NinjaFoam::SampleRawOutput() } } GDALClose( hDS ); + } -void NinjaFoam::WriteOutputFiles() + +void NinjaFoam::generateMassMesh() { - /*-------------------------------------------------------------------*/ - /* prepare output */ - /*-------------------------------------------------------------------*/ + CPLDebug("NINJAFOAM", "generating NINJAFOAM mass mesh grid"); - //Clip off bounding doughnut if desired - VelocityGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); - AngleGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); - if(input.writeTurbulence) - { - TurbulenceGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); - } + + massMesh.buildStandardMesh(input); + + + // no longer need to resize any of the ascii grids, even L and bl_height, as they are already at the mass mesh resolution + // after the dem is resampled to the mesh resolution they are set using the dem resolution, + // and the mesh resolution now is expected to always match the mass solver mesh resolution + CPLDebug("NINJAFOAM", "mass mesh resolution = %f %s", massMesh.meshResolution, lengthUnits::getString(massMesh.meshResolutionUnits).c_str()); + CPLDebug("NINJAFOAM", "mass mesh nrows = %d, ncols = %d, nlayers = %d", input.dem.get_nRows(), input.dem.get_nCols(), massMesh.nlayers); + CPLDebug("NINJAFOAM", "mass mesh minX = %f, maxX = %f, minY = %f, maxY = %f", massMesh.get_minX(), massMesh.get_maxX(), massMesh.get_minY(), massMesh.get_maxY()); + + + writeProbeSampleFile( massMesh.XORD, massMesh.YORD, massMesh.ZORD, input.dem.xllCorner, input.dem.yllCorner, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers ); + + runProbeSample(); + + massMesh_u.allocate(&massMesh); + massMesh_v.allocate(&massMesh); + massMesh_w.allocate(&massMesh); + readInProbeData( massMesh.XORD, massMesh.YORD, massMesh.ZORD, input.dem.xllCorner, input.dem.yllCorner, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers, massMesh_u, massMesh_v, massMesh_w ); + + fillEmptyProbeVals( massMesh.ZORD, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers, massMesh_u, massMesh_v, massMesh_w ); + + + massMesh_k.allocate(&massMesh); + readInProbeData( massMesh.XORD, massMesh.YORD, massMesh.ZORD, input.dem.xllCorner, input.dem.yllCorner, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers, massMesh_k ); + + fillEmptyProbeVals( massMesh.ZORD, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers, massMesh_k ); + + + CPLDebug("NINJAFOAM", "finished generating NINJAFOAM mass mesh grid"); +} - //change windspeed units back to what is specified by speed units switch - velocityUnits::fromBaseUnits(VelocityGrid, input.outputSpeedUnits); - if(input.writeTurbulence) - { - velocityUnits::fromBaseUnits(TurbulenceGrid, input.outputSpeedUnits); +void NinjaFoam::writeProbeSampleFile( const wn_3dArray& x, const wn_3dArray& y, const wn_3dArray& z, + const double dem_xllCorner, const double dem_yllCorner, + const int ncols, const int nrows, const int nlayers) +{ + CPLDebug("NINJAFOAM", "writing probes sample file"); + + const char *probes_filename; + if ( foamVersion == "2.2.0" ) { + probes_filename = CPLFormFilename(pszFoamPath, "system/sampleDict_probes", ""); + } else { + probes_filename = CPLFormFilename(pszFoamPath, "system/probes", ""); } - - //resample to requested output resolutions - SetOutputResolution(); - - //set up filenames - SetOutputFilenames(); - - /*-------------------------------------------------------------------*/ - /* write output files */ - /*-------------------------------------------------------------------*/ - - try{ - if(input.asciiOutFlag==true) - { - AsciiGrid *velTempGrid, *angTempGrid; - velTempGrid=NULL; - angTempGrid=NULL; - - angTempGrid = new AsciiGrid (AngleGrid.resample_Grid(input.angResolution, - AsciiGrid::order0)); - velTempGrid = new AsciiGrid (VelocityGrid.resample_Grid(input.velResolution, - AsciiGrid::order0)); - - //Set cloud grid - int longEdge = input.dem.get_nRows(); - if(input.dem.get_nRows() < input.dem.get_nCols()) - longEdge = input.dem.get_nCols(); - double tempCloudCover; - if(input.cloudCover < 0){ - tempCloudCover = 0.0; - } - else{ - tempCloudCover = input.cloudCover; - } - - CloudGrid.set_headerData(1, 1, input.dem.get_xllCorner(), - input.dem.get_yllCorner(), (longEdge * input.dem.cellSize), - -9999.0, tempCloudCover, input.dem.prjString); - - AsciiGrid tempCloud(CloudGrid); - tempCloud *= 100.0; //Change to percent, which is what FARSITE needs - - // if output clipping was set by the user, don't buffer to overlap the DEM - // but only if writing atm file for farsite grids - if(!input.outputBufferClipping > 0.0 && input.writeAtmFile == true) - { - //ensure grids cover original DEM extents for FARSITE - AsciiGrid demGrid; - GDALDatasetH hDS; - hDS = GDALOpen( input.dem.fileName.c_str(), GA_ReadOnly ); - if( hDS == NULL ) - { - input.Com->ninjaCom(ninjaComClass::ninjaNone, - "Problem reading DEM during output writing." ); - } - GDAL2AsciiGrid( (GDALDataset *)hDS, 1, demGrid ); - tempCloud.BufferToOverlapGrid(demGrid); - angTempGrid->BufferToOverlapGrid(demGrid); - velTempGrid->BufferToOverlapGrid(demGrid); - } - - tempCloud.write_Grid(input.cldFile.c_str(), 1); - angTempGrid->write_Grid(input.angFile.c_str(), 0); - velTempGrid->write_Grid(input.velFile.c_str(), 2); - - if(angTempGrid) - { - delete angTempGrid; - angTempGrid=NULL; - } - if(velTempGrid) - { - delete velTempGrid; - velTempGrid=NULL; - } - - //Write .atm file for this run. Only has one time value in file. - if(input.writeAtmFile) - { - farsiteAtm atmosphere; - atmosphere.push(input.ninjaTime, input.velFile, input.angFile, input.cldFile); - atmosphere.writeAtmFile(input.atmFile, input.outputSpeedUnits, input.outputWindHeight); - } - } - }catch (exception& e) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during ascii file writing: %s", e.what()); - }catch (...) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during ascii file writing: Cannot determine exception type."); - } - - //write text file comparing measured to simulated winds (measured read from file, filename, etc. hard-coded in function) - try{ - if(input.txtOutFlag==true) - write_compare_output(); - }catch (exception& e) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during text file writing: %s", e.what()); - }catch (...) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during text file writing: Cannot determine exception type."); - } - - //write shape files - try{ - if(input.shpOutFlag==true) - { - AsciiGrid *velTempGrid, *angTempGrid; - velTempGrid=NULL; - angTempGrid=NULL; - - ShapeVector ninjaShapeFiles; - - angTempGrid = new AsciiGrid (AngleGrid.resample_Grid(input.shpResolution, AsciiGrid::order0)); - velTempGrid = new AsciiGrid (VelocityGrid.resample_Grid(input.shpResolution, AsciiGrid::order0)); - - ninjaShapeFiles.setDirGrid(*angTempGrid); - ninjaShapeFiles.setSpeedGrid(*velTempGrid); - ninjaShapeFiles.setDataBaseName(input.dbfFile); - ninjaShapeFiles.setShapeFileName(input.shpFile); - ninjaShapeFiles.makeShapeFiles(); - - if(angTempGrid) - { - delete angTempGrid; - angTempGrid=NULL; - } - if(velTempGrid) - { - delete velTempGrid; - velTempGrid=NULL; - } - } - }catch (exception& e) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during shape file writing: %s", e.what()); - }catch (...) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during shape file writing: Cannot determine exception type."); - } - - //write kmz file - try{ - if(input.googOutFlag==true) - - { - AsciiGrid *velTempGrid, *angTempGrid, *turbTempGrid; - velTempGrid=NULL; - angTempGrid=NULL; - turbTempGrid=NULL; - - KmlVector ninjaKmlFiles; - - angTempGrid = new AsciiGrid (AngleGrid.resample_Grid(input.kmzResolution, - AsciiGrid::order0)); - velTempGrid = new AsciiGrid (VelocityGrid.resample_Grid(input.kmzResolution, - AsciiGrid::order0)); - if(input.writeTurbulence) - { - turbTempGrid = new AsciiGrid (TurbulenceGrid.resample_Grid(input.kmzResolution, - AsciiGrid::order0)); - - ninjaKmlFiles.setTurbulenceFlag("true"); - ninjaKmlFiles.setTurbulenceGrid(*turbTempGrid, input.outputSpeedUnits); - } - - ninjaKmlFiles.setKmlFile(input.kmlFile); - ninjaKmlFiles.setKmzFile(input.kmzFile); - ninjaKmlFiles.setDemFile(input.dem.fileName); - - ninjaKmlFiles.setLegendFile(input.legFile); - ninjaKmlFiles.setDateTimeLegendFile(input.dateTimeLegFile, input.ninjaTime); - ninjaKmlFiles.setSpeedGrid(*velTempGrid, input.outputSpeedUnits); - ninjaKmlFiles.setDirGrid(*angTempGrid); - - ninjaKmlFiles.setLineWidth(input.googLineWidth); - ninjaKmlFiles.setTime(input.ninjaTime); - - if(ninjaKmlFiles.writeKml(input.googSpeedScaling,input.googColor,input.googVectorScale)) - { - if(ninjaKmlFiles.makeKmz()) - ninjaKmlFiles.removeKmlFile(); - } - if(angTempGrid) - { - delete angTempGrid; - angTempGrid=NULL; - } - if(velTempGrid) - { - delete velTempGrid; - velTempGrid=NULL; - } - if(turbTempGrid) - { - delete turbTempGrid; - turbTempGrid=NULL; - } - } - }catch (exception& e) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during Google Earth file writing: %s", e.what()); - }catch (...) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during Google Earth file writing: Cannot determine exception type."); - } - - try{ - if(input.pdfOutFlag==true) - { - AsciiGrid *velTempGrid, *angTempGrid; - velTempGrid=NULL; - angTempGrid=NULL; - OutputWriter output; - - angTempGrid = new AsciiGrid (AngleGrid.resample_Grid(input.pdfResolution, AsciiGrid::order0)); - velTempGrid = new AsciiGrid (VelocityGrid.resample_Grid(input.pdfResolution, AsciiGrid::order0)); - - output.setDirGrid(*angTempGrid); - output.setSpeedGrid(*velTempGrid, input.outputSpeedUnits); - output.setDEMfile(input.pdfDEMFileName); - output.setLineWidth(input.pdfLineWidth); - output.setDPI(input.pdfDPI); - output.setSize(input.pdfWidth, input.pdfHeight); - output.write(input.pdfFile, "PDF"); - - - if(angTempGrid) - { - delete angTempGrid; - angTempGrid=NULL; - } - if(velTempGrid) - { - delete velTempGrid; - velTempGrid=NULL; - } - } - }catch (exception& e) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during pdf file writing: %s", e.what()); - }catch (...) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during pdf file writing: Cannot determine exception type."); - } - - - try{ - if ( writeMassMeshVtk == true ) { - CPLDebug("NINJAFOAM", "writing mass mesh vtk output for foam simulation."); - writeMassMeshVtkOutput(); - } - }catch (exception& e) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during NINJAFOAM mass mesh vtk file writing: %s", e.what()); - }catch (...) - { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during NINJAFOAM mass mesh vtk file writing: Cannot determine exception type."); - } - -} - - -void NinjaFoam::writeMassMeshVtkOutput() -{ - - massMesh.buildStandardMesh(input); - - - // no longer need to resize any of the ascii grids, even L and bl_height, as they are already at the mass mesh resolution - // after the dem is resampled to the mesh resolution they are set using the dem resolution, - // and the mesh resolution now is expected to always match the mass solver mesh resolution - CPLDebug("NINJAFOAM", "mass mesh resolution = %f %s", massMesh.meshResolution, lengthUnits::getString(massMesh.meshResolutionUnits).c_str()); - CPLDebug("NINJAFOAM", "mass mesh nrows = %d, ncols = %d, nlayers = %d", input.dem.get_nRows(), input.dem.get_nCols(), massMesh.nlayers); - CPLDebug("NINJAFOAM", "mass mesh minX = %f, maxX = %f, minY = %f, maxY = %f", massMesh.get_minX(), massMesh.get_maxX(), massMesh.get_minY(), massMesh.get_maxY()); - - - writeProbeSampleFile( massMesh.XORD, massMesh.YORD, massMesh.ZORD, input.dem.xllCorner, input.dem.yllCorner, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers ); - - runProbeSample(); - - wn_3dScalarField u, v, w; - u.allocate(&massMesh); - v.allocate(&massMesh); - w.allocate(&massMesh); - readInProbeData( massMesh.XORD, massMesh.YORD, massMesh.ZORD, input.dem.xllCorner, input.dem.yllCorner, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers, u, v, w ); + CPLDebug("NINJAFOAM", "probes sample filename = \"%s\"", probes_filename); - fillEmptyProbeVals( massMesh.ZORD, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers, u, v, w ); + FILE *fout; - std::string massMeshVtkFilename = CPLFormFilename(pszFoamPath, "massMesh", "vtk"); - try { - CPLDebug("NINJAFOAM", "writing vtk file"); - bool vtk_out_as_utm = false; - if(CSLTestBoolean(CPLGetConfigOption("VTK_OUT_AS_UTM", "FALSE"))) - { - vtk_out_as_utm = CPLGetConfigOption("VTK_OUT_AS_UTM", "FALSE"); - } - // can pick between "ascii" and "binary" format for the vtk write format - std::string vtkWriteFormat = "ascii";//"binary";//"ascii"; - volVTK VTK(u, v, w, massMesh.XORD, massMesh.YORD, massMesh.ZORD, input.dem.get_xllCorner(), input.dem.get_yllCorner(), input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers, massMeshVtkFilename, vtkWriteFormat, vtk_out_as_utm); - } catch (exception& e) { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during volume VTK file writing: %s", e.what()); - } catch (...) { - input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during volume VTK file writing: Cannot determine exception type."); - } - - u.deallocate(); - v.deallocate(); - w.deallocate(); -} - -void NinjaFoam::writeProbeSampleFile( const wn_3dArray& x, const wn_3dArray& y, const wn_3dArray& z, - const double dem_xllCorner, const double dem_yllCorner, - const int ncols, const int nrows, const int nlayers) -{ - CPLDebug("NINJAFOAM", "writing probes sample file"); - - const char *probes_filename; - if ( foamVersion == "2.2.0" ) { - probes_filename = CPLFormFilename(pszFoamPath, "system/sampleDict_probes", ""); - } else { - probes_filename = CPLFormFilename(pszFoamPath, "system/probes", ""); - } - CPLDebug("NINJAFOAM", "probes sample filename = \"%s\"", probes_filename); - - FILE *fout; - - fout = fopen(probes_filename, "w"); - if( fout == NULL ) - throw std::runtime_error("probes_filename cannot be opened for writing."); + fout = fopen(probes_filename, "w"); + if( fout == NULL ) + throw std::runtime_error("probes_filename cannot be opened for writing."); // Write header stuff if ( foamVersion == "2.2.0" ) { @@ -2774,8 +2386,8 @@ void NinjaFoam::writeProbeSampleFile( const wn_3dArray& x, const wn_3dArray& y, fprintf(fout, "\n"); fprintf(fout, "\n"); fprintf(fout, "// choice of variables\n"); - //fprintf(fout, "fields (%s %s);\n", "U", "k"); - fprintf(fout, "fields (%s %s %s %s %s);\n", "U", "k", "epsilon", "nut", "p"); + fprintf(fout, "fields (%s %s);\n", "U", "k"); + //fprintf(fout, "fields (%s %s %s %s %s);\n", "U", "k", "epsilon", "nut", "p"); fprintf(fout, "\n"); fprintf(fout, "\n"); fprintf(fout, "// Sampling and I/O settings\n"); @@ -3103,126 +2715,889 @@ void NinjaFoam::readInProbeData( const wn_3dArray& x, const wn_3dArray& y, const } -void NinjaFoam::fillEmptyProbeVals(const wn_3dArray& z, - const int ncols, const int nrows, const int nlayers, - wn_3dScalarField& u, wn_3dScalarField& v, wn_3dScalarField& w) +void NinjaFoam::readInProbeData( const wn_3dArray& x, const wn_3dArray& y, const wn_3dArray& z, + const double dem_xllCorner, const double dem_yllCorner, + const int ncols, const int nrows, const int nlayers, + wn_3dScalarField& k ) +{ + CPLDebug("NINJAFOAM", "reading in probes sample data"); + + const char *probesPostProcessDirname = "probes"; + if ( foamVersion == "2.2.0" ) { + probesPostProcessDirname = "sets"; + } + + // method from surfaces sampling, to find the time directory for the surfaces file + char **papszOutputProbeDataPath; + papszOutputProbeDataPath = VSIReadDir(CPLSPrintf("%s/postProcessing/%s/", pszFoamPath, probesPostProcessDirname)); + + const char *probeSampleData_filename; + const char *timeDir; + for(int i = 0; i < CSLCount( papszOutputProbeDataPath ); i++){ + if(std::string(papszOutputProbeDataPath[i]) != "." && + std::string(papszOutputProbeDataPath[i]) != "..") { + timeDir = papszOutputProbeDataPath[i]; + probeSampleData_filename = CPLSPrintf("%s/postProcessing/%s/%s/points_k.xy", pszFoamPath, probesPostProcessDirname, timeDir); + break; + } + else{ + continue; + } + } + CPLDebug("NINJAFOAM", "probes sample data file path = \"%s\"", probeSampleData_filename); + + + // read the full data file into a string separated by "\n" chars for each line + VSILFILE *fin; + fin = VSIFOpenL( probeSampleData_filename, "r" ); + + char *data; + + vsi_l_offset offset; + VSIFSeekL(fin, 0, SEEK_END); + offset = VSIFTellL(fin); + + VSIRewindL(fin); + data = (char*)CPLMalloc(offset * sizeof(char) + 1); + VSIFReadL(data, offset, 1, fin); + data[offset] = '\0'; + + std::string probeSampleData_allLines(data); + + CPLFree(data); + VSIFCloseL( fin ); + + + // now go through the data string line by line from where the probe data starts, + // comparing the probe data points with the mesh points to detect whether any data points were silently dropped for being outside the mesh + // filling any dropped probe data with NoData vals to be filled at later steps + double noDataVal = -9999; + int startLinePos = 0; // data starts right at the first character of the file, so no need to set it by finding a specific spot/character in the file + int endLinePos; + for(int layerIdx=0; layerIdx noDataCheckVal or if < noDataCheckVal depending on the style of the noDataVal + + // going to first do a log profile fill for all cells with noData vals up to the lowest known z value for each column, + // then do an ascii grid style nan filling to fill in any remaining noData vals + // + // currently, the log profile filling seems to get all the noData vals, but keep in mind the method in case it ever becomes relevant + // + // the lowest known z value is the first non noData val in a given column, or the 1stCellHeight value, whichever is higher up in the column + // so if values below 1stCellHeight have data, they are also replaced with log profile interpolation + + windProfile profile; + profile.profile_switch = windProfile::monin_obukov_similarity; + + for(int rowIdx=0; rowIdx noDataCheckVal ) { // they always get set together for the same idx, so only need to check one of them + //if ( u(rowIdx,colIdx,layerIdx) > noDataCheckVal && v(rowIdx,colIdx,layerIdx) > noDataCheckVal && w(rowIdx,colIdx,layerIdx) > noDataCheckVal ) { + lowestKnown_zIdx = layerIdx; + break; + } + } + + if ( lowestKnown_zIdx == nlayers-1 ) { + // is a column of no data values, skip it for the log profile part of the fill + // also warn because the method for filling no data values past the log profile part hasn't yet been implemented + std::cout << "!!! no lowest known zIdx for column of data !!! for rowIdx = " << rowIdx << ", colIdx = " << colIdx << std::endl; + continue; + } + + + double current_z_ground = z(rowIdx,colIdx,0); + + double firstCellHeight_zIdx = 0; + for(int layerIdx=0; layerIdx= finalFirstCellHeight ) + { + firstCellHeight_zIdx = layerIdx; + break; + } + } + + + // for debugging purposes, want to know how different the indices and resulting values ended up + // turns out to be easier if change back and forth between when/why it is printed rather than printing each instance + //if ( lowestKnown_zIdx == firstCellHeight_zIdx ) { + //if ( lowestKnown_zIdx > firstCellHeight_zIdx ) { + // std::cout << "rowIdx = " << rowIdx << ", colIdx = " << colIdx << std::endl; + // std::cout << "lowestKnown_zIdx = " << lowestKnown_zIdx << ", firstCellHeight_zIdx = " << firstCellHeight_zIdx << std::endl; + // std::cout << "z_ASL[lowestKnown_zIdx] = " << z(rowIdx,colIdx,lowestKnown_zIdx) << ", current_z_ground = " << current_z_ground << ", z_AGL[lowestKnown_zIdx] = " << z(rowIdx,colIdx,lowestKnown_zIdx)-current_z_ground << ", firstCellHeight z = " << finalFirstCellHeight << std::endl; + //} + + + // this one is NOT for debugging purposes, has to be run each and every time to make the code run safely + if ( lowestKnown_zIdx < firstCellHeight_zIdx ) { + //std::cout << "lowestKnown zIdx is less than firstCellHeight zIdx for rowIdx = " << rowIdx << ", colIdx = " << colIdx << std::endl; + //std::cout << "old lowestKnown_zIdx = " << lowestKnown_zIdx << std::endl; + lowestKnown_zIdx = firstCellHeight_zIdx; + //std::cout << "new lowestKnown_zIdx = " << lowestKnown_zIdx << std::endl; + //std::cout << "new z_ASL[lowestKnown_zIdx] = " << z(rowIdx,colIdx,lowestKnown_zIdx) << ", current_z_ground = " << current_z_ground << ", new z_AGL[lowestKnown_zIdx] = " << z(rowIdx,colIdx,lowestKnown_zIdx)-current_z_ground << ", firstCellHeight z = " << finalFirstCellHeight << std::endl; + // looks like this is getting triggered a LOT, if not for each and every occurence, hard to tell without commenting it out and printing for the other occurences. + // looks like the resulting value is rarely equal to firstCellHeight though, seems to always be just a hint above that value + } + + + + double z_ref = z(rowIdx,colIdx,lowestKnown_zIdx); + double u_ref = u(rowIdx,colIdx,lowestKnown_zIdx); + double v_ref = v(rowIdx,colIdx,lowestKnown_zIdx); + double w_ref = w(rowIdx,colIdx,lowestKnown_zIdx); + + + profile.ObukovLength = init->L(rowIdx,colIdx); + profile.ABL_height = init->bl_height(rowIdx,colIdx); + profile.Roughness = input.surface.Roughness(rowIdx,colIdx); + double current_rough_h = input.surface.Rough_h(rowIdx,colIdx); + double current_rough_d = input.surface.Rough_d(rowIdx,colIdx); + profile.Rough_h = current_rough_h; + profile.Rough_d = current_rough_d; + + //this is height above the vegetation + profile.inputWindHeight = z_ref - current_z_ground; // needs to be AGL + + for(int zIdx=1; zIdx noDataCheckVal or if < noDataCheckVal depending on the style of the noDataVal + + // find the first non nan value for a given column, and fill all values below with that value + // so no log profile correction for this dataset + + for(int rowIdx=0; rowIdx noDataCheckVal ) { + lowestKnown_zIdx = layerIdx; + break; + } + } + + if ( lowestKnown_zIdx == nlayers-1 ) { + // is a column of no data values, skip it for the log profile part of the fill + // also warn because the method for filling no data values past the log profile part hasn't yet been implemented + std::cout << "!!! no lowest known zIdx for column of data !!! for rowIdx = " << rowIdx << ", colIdx = " << colIdx << std::endl; + continue; + } + + + double lowestKnown_zVal = k(rowIdx,colIdx,lowestKnown_zIdx); + + for(int zIdx=0; zIdx colHeightAGL ) { + // looked at all values up to the colHeightAGL, move on + break; + } + if ( current_colMaxVal < k(rowIdx,colIdx,layerIdx) ) { + current_colMaxVal = k(rowIdx,colIdx,layerIdx); + } + //// for debugging + //if ( rowIdx == 0 && colIdx == 0 ) { + // std::cout << "zIdx = " << layerIdx << ", z = " << z(rowIdx,colIdx,layerIdx) << ", z_AGL = " << z(rowIdx,colIdx,layerIdx) - current_z_ground << ", val = " << k(rowIdx,colIdx,layerIdx) << std::endl; + //} + } + + colMaxGrid.set_cellValue( rowIdx, colIdx, current_colMaxVal ); + } + } + + //// for debugging + //CPLDebug("NINJAFOAM", "writing ascii file"); + //const char *colMaxOutputFile_ascii; + //colMaxOutputFile_ascii = CPLSPrintf("%s/colMax.asc", pszFoamPath); + //colMaxGrid.write_Grid(colMaxOutputFile_ascii, 5); + + CPLDebug("NINJAFOAM", "finished generating turbulence column max ascii grid from NINJAFOAM mass mesh"); +} + + +void NinjaFoam::GenerateAndSampleMassMesh() +{ + + // generate massMesh if required for other outputs + try{ + if ( writeMassMesh == true ) { + generateMassMesh(); + } + }catch (exception& e) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during NINJAFOAM mass mesh generation: %s", e.what()); + }catch (...) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during NINJAFOAM mass mesh generation: Cannot determine exception type."); + } + + + // prep colMaxGrid + try{ + if ( input.writeTurbulence == true ) { + generateColMaxGrid(massMesh.ZORD, input.dem.xllCorner, input.dem.yllCorner, input.dem.get_nCols(), input.dem.get_nRows(), massMesh.nlayers, massMesh.meshResolution, input.dem.prjString, massMesh_k); + } + }catch (exception& e) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during turbulence column max from NINJAFOAM mass mesh ascii grid generation: %s", e.what()); + }catch (...) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during turbulence column max from NINJAFOAM mass mesh ascii grid generation: Cannot determine exception type."); + } + +} + + +void NinjaFoam::SetOutputResolution() +{ + //Set output file resolutions now + if( input.kmzResolution <= 0.0 ) //if negative, use DEM resolution + input.kmzResolution = input.dem.get_cellSize(); + if( input.shpResolution <= 0.0 ) //if negative, use DEM resolution + input.shpResolution = input.dem.get_cellSize(); + if( input.velResolution <= 0.0 ) //if negative, use DEMresolution + input.velResolution = input.dem.get_cellSize(); + if( input.angResolution <= 0.0 ) //if negative, use DEM resolution + input.angResolution = input.dem.get_cellSize(); + if( input.pdfResolution <= 0.0 ) + input.pdfResolution = input.dem.get_cellSize(); +} + +void NinjaFoam::SetOutputFilenames() +{ + //Do file naming string stuff for all output files + std::string rootFile, rootName, timeAppend, wxModelTimeAppend, fileAppend, kmz_fileAppend, \ + shp_fileAppend, ascii_fileAppend, mesh_units, kmz_mesh_units, \ + shp_mesh_units, ascii_mesh_units, pdf_fileAppend, pdf_mesh_units; + + boost::local_time::local_time_facet* timeOutputFacet; + timeOutputFacet = new boost::local_time::local_time_facet(); + //NOTE: WEIRD ISSUE WITH THE ABOVE 2 LINES OF CODE! DO NOT CALL DELETE ON THIS BECAUSE THE LOCALE OBJECT BELOW DOES. + // THIS IS A "PROBLEM" IN THE STANDARD LIBRARY. SEE THESE WEB SITES FOR MORE INFO: + // https://collab.firelab.org/software/projects/windninja/wiki/KnownIssues + // http://rhubbarb.wordpress.com/2009/10/17/boost-datetime-locales-and-facets/#comment-203 + + std::ostringstream timestream; + timestream.imbue(std::locale(std::locale::classic(), timeOutputFacet)); + timeOutputFacet->format("%m-%d-%Y_%H%M_"); + + if( input.initializationMethod == WindNinjaInputs::wxModelInitializationFlag ){ + timestream << input.ninjaTime; + } + + std::string pathName; + std::string baseName(CPLGetBasename(input.dem.fileName.c_str())); + + if(input.customOutputPath == "!set"){ + //prepend directory paths to rootFile for wxModel run + if( input.initializationMethod == WindNinjaInputs::wxModelInitializationFlag ){ + pathName = CPLGetPath(input.forecastFilename.c_str()); + //if it's a .tar, write to directory containing the .tar file + if( strstr(pathName.c_str(), ".tar") ){ + pathName.erase( pathName.rfind("/") ); + } + }else{ + pathName = CPLGetPath(input.dem.fileName.c_str()); + } + } + else{ + pathName = input.customOutputPath; + } + + rootFile = CPLFormFilename(pathName.c_str(), baseName.c_str(), NULL); + + /* set the output path member variable */ + input.outputPath = pathName; + + timeAppend = timestream.str(); + + ostringstream wxModelTimestream; + boost::local_time::local_time_facet* wxModelOutputFacet; + wxModelOutputFacet = new boost::local_time::local_time_facet(); + wxModelTimestream.imbue(std::locale(std::locale::classic(), wxModelOutputFacet)); + wxModelOutputFacet->format("%m-%d-%Y_%H%M"); + if( input.initializationMethod == WindNinjaInputs::wxModelInitializationFlag) + { + wxModelTimestream << input.ninjaTime; + } + wxModelTimeAppend = wxModelTimestream.str(); + mesh_units = "m"; + kmz_mesh_units = lengthUnits::getString( input.kmzUnits ); + shp_mesh_units = lengthUnits::getString( input.shpUnits ); + ascii_mesh_units = lengthUnits::getString( input.velOutputFileDistanceUnits ); + pdf_mesh_units = lengthUnits::getString( input.pdfUnits ); + + ostringstream os, os_kmz, os_shp, os_ascii, os_pdf; + + if( input.initializationMethod == WindNinjaInputs::domainAverageInitializationFlag ){ + double tempSpeed = input.inputSpeed; + velocityUnits::fromBaseUnits(tempSpeed, input.inputSpeedUnits); + os << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); + os_kmz << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); + os_shp << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); + os_ascii << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); + os_pdf << "_" << (long) (input.inputDirection+0.5) << "_" << (long) (tempSpeed+0.5); + } + + double meshResolutionTemp = input.dem.get_cellSize(); + double kmzResolutionTemp = input.kmzResolution; + double shpResolutionTemp = input.shpResolution; + double velResolutionTemp = input.velResolution; + double pdfResolutionTemp = input.pdfResolution; + + lengthUnits::eLengthUnits meshResolutionUnits = lengthUnits::meters; + + lengthUnits::fromBaseUnits(meshResolutionTemp, meshResolutionUnits); + lengthUnits::fromBaseUnits(kmzResolutionTemp, input.kmzUnits); + lengthUnits::fromBaseUnits(shpResolutionTemp, input.shpUnits); + lengthUnits::fromBaseUnits(velResolutionTemp, input.velOutputFileDistanceUnits); + lengthUnits::fromBaseUnits(pdfResolutionTemp, input.pdfUnits); + + os << "_" << timeAppend << (long) (meshResolutionTemp+0.5) << mesh_units; + os_kmz << "_" << timeAppend << (long) (kmzResolutionTemp+0.5) << kmz_mesh_units; + os_shp << "_" << timeAppend << (long) (shpResolutionTemp+0.5) << shp_mesh_units; + os_ascii << "_" << timeAppend << (long) (velResolutionTemp+0.5) << ascii_mesh_units; + os_pdf << "_" << timeAppend << (long) (pdfResolutionTemp+0.5) << pdf_mesh_units; + + fileAppend = os.str(); + kmz_fileAppend = os_kmz.str(); + shp_fileAppend = os_shp.str(); + ascii_fileAppend = os_ascii.str(); + pdf_fileAppend = os_pdf.str(); + + input.kmlFile = rootFile + kmz_fileAppend + ".kml"; + input.kmzFile = rootFile + kmz_fileAppend + ".kmz"; + + input.shpFile = rootFile + shp_fileAppend + ".shp"; + input.dbfFile = rootFile + shp_fileAppend + ".dbf"; + + input.pdfFile = rootFile + pdf_fileAppend + ".pdf"; + + input.cldFile = rootFile + ascii_fileAppend + "_cld.asc"; + input.velFile = rootFile + ascii_fileAppend + "_vel.asc"; + input.angFile = rootFile + ascii_fileAppend + "_ang.asc"; + input.atmFile = rootFile + ascii_fileAppend + ".atm"; + + input.legFile = rootFile + kmz_fileAppend + ".bmp"; + if( input.ninjaTime.is_not_a_date_time() ) //date and time not set? + input.dateTimeLegFile = ""; + else + input.dateTimeLegFile = rootFile + kmz_fileAppend + ".date_time" + ".bmp"; +} + +void NinjaFoam::WriteOutputFiles() +{ + /*-------------------------------------------------------------------*/ + /* prepare output */ + /*-------------------------------------------------------------------*/ + + //Clip off bounding doughnut if desired + VelocityGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); + AngleGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); + if(input.writeTurbulence) + { + TurbulenceGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); + colMaxGrid.clipGridInPlaceSnapToCells(input.outputBufferClipping); + } + + //change windspeed units back to what is specified by speed units switch + velocityUnits::fromBaseUnits(VelocityGrid, input.outputSpeedUnits); + if(input.writeTurbulence) + { + velocityUnits::fromBaseUnits(TurbulenceGrid, input.outputSpeedUnits); + velocityUnits::fromBaseUnits(colMaxGrid, input.outputSpeedUnits); + } + + //resample to requested output resolutions + SetOutputResolution(); + + //set up filenames + SetOutputFilenames(); + + /*-------------------------------------------------------------------*/ + /* write output files */ + /*-------------------------------------------------------------------*/ + + try{ + if(input.asciiOutFlag==true) + { + AsciiGrid *velTempGrid, *angTempGrid; + velTempGrid=NULL; + angTempGrid=NULL; + + angTempGrid = new AsciiGrid (AngleGrid.resample_Grid(input.angResolution, + AsciiGrid::order0)); + velTempGrid = new AsciiGrid (VelocityGrid.resample_Grid(input.velResolution, + AsciiGrid::order0)); + + //Set cloud grid + int longEdge = input.dem.get_nRows(); + if(input.dem.get_nRows() < input.dem.get_nCols()) + longEdge = input.dem.get_nCols(); + double tempCloudCover; + if(input.cloudCover < 0){ + tempCloudCover = 0.0; + } + else{ + tempCloudCover = input.cloudCover; + } + + CloudGrid.set_headerData(1, 1, input.dem.get_xllCorner(), + input.dem.get_yllCorner(), (longEdge * input.dem.cellSize), + -9999.0, tempCloudCover, input.dem.prjString); + + AsciiGrid tempCloud(CloudGrid); + tempCloud *= 100.0; //Change to percent, which is what FARSITE needs + + // if output clipping was set by the user, don't buffer to overlap the DEM + // but only if writing atm file for farsite grids + if(!input.outputBufferClipping > 0.0 && input.writeAtmFile == true) + { + //ensure grids cover original DEM extents for FARSITE + AsciiGrid demGrid; + GDALDatasetH hDS; + hDS = GDALOpen( input.dem.fileName.c_str(), GA_ReadOnly ); + if( hDS == NULL ) + { + input.Com->ninjaCom(ninjaComClass::ninjaNone, + "Problem reading DEM during output writing." ); + } + GDAL2AsciiGrid( (GDALDataset *)hDS, 1, demGrid ); + tempCloud.BufferToOverlapGrid(demGrid); + angTempGrid->BufferToOverlapGrid(demGrid); + velTempGrid->BufferToOverlapGrid(demGrid); + } + + tempCloud.write_Grid(input.cldFile.c_str(), 1); + angTempGrid->write_Grid(input.angFile.c_str(), 0); + velTempGrid->write_Grid(input.velFile.c_str(), 2); + + if(angTempGrid) + { + delete angTempGrid; + angTempGrid=NULL; + } + if(velTempGrid) + { + delete velTempGrid; + velTempGrid=NULL; + } + + //Write .atm file for this run. Only has one time value in file. + if(input.writeAtmFile) + { + farsiteAtm atmosphere; + atmosphere.push(input.ninjaTime, input.velFile, input.angFile, input.cldFile); + atmosphere.writeAtmFile(input.atmFile, input.outputSpeedUnits, input.outputWindHeight); + } + } + }catch (exception& e) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during ascii file writing: %s", e.what()); + }catch (...) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during ascii file writing: Cannot determine exception type."); + } + + //write text file comparing measured to simulated winds (measured read from file, filename, etc. hard-coded in function) + try{ + if(input.txtOutFlag==true) + write_compare_output(); + }catch (exception& e) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during text file writing: %s", e.what()); + }catch (...) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during text file writing: Cannot determine exception type."); + } + + //write shape files + try{ + if(input.shpOutFlag==true) + { + AsciiGrid *velTempGrid, *angTempGrid; + velTempGrid=NULL; + angTempGrid=NULL; + + ShapeVector ninjaShapeFiles; + + angTempGrid = new AsciiGrid (AngleGrid.resample_Grid(input.shpResolution, AsciiGrid::order0)); + velTempGrid = new AsciiGrid (VelocityGrid.resample_Grid(input.shpResolution, AsciiGrid::order0)); + + ninjaShapeFiles.setDirGrid(*angTempGrid); + ninjaShapeFiles.setSpeedGrid(*velTempGrid); + ninjaShapeFiles.setDataBaseName(input.dbfFile); + ninjaShapeFiles.setShapeFileName(input.shpFile); + ninjaShapeFiles.makeShapeFiles(); + + if(angTempGrid) + { + delete angTempGrid; + angTempGrid=NULL; + } + if(velTempGrid) + { + delete velTempGrid; + velTempGrid=NULL; + } + } + }catch (exception& e) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during shape file writing: %s", e.what()); + }catch (...) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during shape file writing: Cannot determine exception type."); + } + + //write kmz file + try{ + if(input.googOutFlag==true) + + { + AsciiGrid *velTempGrid, *angTempGrid, *turbTempGrid, *colMaxTempGrid; + velTempGrid=NULL; + angTempGrid=NULL; + turbTempGrid=NULL; + colMaxTempGrid=NULL; + + KmlVector ninjaKmlFiles; + + angTempGrid = new AsciiGrid (AngleGrid.resample_Grid(input.kmzResolution, + AsciiGrid::order0)); + velTempGrid = new AsciiGrid (VelocityGrid.resample_Grid(input.kmzResolution, + AsciiGrid::order0)); + if(input.writeTurbulence) + { + //turbTempGrid = new AsciiGrid (TurbulenceGrid.resample_Grid(input.kmzResolution, + // AsciiGrid::order0)); + // + //ninjaKmlFiles.setTurbulenceFlag("true"); + //ninjaKmlFiles.setTurbulenceGrid(*turbTempGrid, input.outputSpeedUnits); + + + colMaxTempGrid = new AsciiGrid (colMaxGrid.resample_Grid(input.kmzResolution, + AsciiGrid::order0)); + + ninjaKmlFiles.setColMaxFlag("true"); + ninjaKmlFiles.setColMaxGrid(*colMaxTempGrid, input.outputSpeedUnits, input.colMax_colHeightAGL, input.colMax_colHeightAGL_units); + + //// for debugging + //colMaxGrid.ascii2png("colMaxGrid.png", "Speed Fluctuation", "(mph)", "colMax_legend", true, false); + //colMaxTempGrid->ascii2png("colMaxTempGrid.png", "Speed Fluctuation", "(mph)", "colMax_legend", true, false); + } + + ninjaKmlFiles.setKmlFile(input.kmlFile); + ninjaKmlFiles.setKmzFile(input.kmzFile); + ninjaKmlFiles.setDemFile(input.dem.fileName); + + ninjaKmlFiles.setLegendFile(input.legFile); + ninjaKmlFiles.setDateTimeLegendFile(input.dateTimeLegFile, input.ninjaTime); + ninjaKmlFiles.setSpeedGrid(*velTempGrid, input.outputSpeedUnits); + ninjaKmlFiles.setDirGrid(*angTempGrid); + + ninjaKmlFiles.setLineWidth(input.googLineWidth); + ninjaKmlFiles.setTime(input.ninjaTime); + + if(ninjaKmlFiles.writeKml(input.googSpeedScaling,input.googColor,input.googVectorScale)) + { + if(ninjaKmlFiles.makeKmz()) + ninjaKmlFiles.removeKmlFile(); + } + if(angTempGrid) + { + delete angTempGrid; + angTempGrid=NULL; + } + if(velTempGrid) + { + delete velTempGrid; + velTempGrid=NULL; + } + if(turbTempGrid) + { + delete turbTempGrid; + turbTempGrid=NULL; + } + if(colMaxTempGrid) + { + delete colMaxTempGrid; + colMaxTempGrid=NULL; + } + } + }catch (exception& e) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during Google Earth file writing: %s", e.what()); + }catch (...) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during Google Earth file writing: Cannot determine exception type."); + } + + try{ + if(input.pdfOutFlag==true) + { + AsciiGrid *velTempGrid, *angTempGrid; + velTempGrid=NULL; + angTempGrid=NULL; + OutputWriter output; + + angTempGrid = new AsciiGrid (AngleGrid.resample_Grid(input.pdfResolution, AsciiGrid::order0)); + velTempGrid = new AsciiGrid (VelocityGrid.resample_Grid(input.pdfResolution, AsciiGrid::order0)); + + output.setDirGrid(*angTempGrid); + output.setSpeedGrid(*velTempGrid, input.outputSpeedUnits); + output.setDEMfile(input.pdfDEMFileName); + output.setLineWidth(input.pdfLineWidth); + output.setDPI(input.pdfDPI); + output.setSize(input.pdfWidth, input.pdfHeight); + output.write(input.pdfFile, "PDF"); + + + if(angTempGrid) + { + delete angTempGrid; + angTempGrid=NULL; + } + if(velTempGrid) + { + delete velTempGrid; + velTempGrid=NULL; + } + } + }catch (exception& e) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during pdf file writing: %s", e.what()); + }catch (...) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during pdf file writing: Cannot determine exception type."); + } + + + try{ + if ( writeMassMeshVtk == true ) { + writeMassMeshVtkOutput(); + } + }catch (exception& e) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during NINJAFOAM mass mesh vtk file writing: %s", e.what()); + }catch (...) + { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during NINJAFOAM mass mesh vtk file writing: Cannot determine exception type."); + } + +} + + +void NinjaFoam::writeMassMeshVtkOutput() { - CPLDebug("NINJAFOAM", "filling in probes sample no data vals"); - - //double noDataVal = -9999; // make sure this one matches the one used in the readInProbeData() function - double noDataCheckVal = -9998; // instead of using the if == noDataVal, use if > noDataCheckVal or if < noDataCheckVal depending on the style of the noDataVal - - // going to first do a log profile fill for all cells with noData vals up to the lowest known z value for each column, - // then do an ascii grid style nan filling to fill in any remaining noData vals - // - // currently, the log profile filling seems to get all the noData vals, but keep in mind the method in case it ever becomes relevant - // - // the lowest known z value is the first non noData val in a given column, or the 1stCellHeight value, whichever is higher up in the column - // so if values below 1stCellHeight have data, they are also replaced with log profile interpolation - - windProfile profile; - profile.profile_switch = windProfile::monin_obukov_similarity; + CPLDebug("NINJAFOAM", "writing mass mesh vtk output for foam simulation."); - for(int rowIdx=0; rowIdx noDataCheckVal ) { // they always get set together for the same idx, so only need to check one of them - //if ( u(rowIdx,colIdx,layerIdx) > noDataCheckVal && v(rowIdx,colIdx,layerIdx) > noDataCheckVal && w(rowIdx,colIdx,layerIdx) > noDataCheckVal ) { - lowestKnown_zIdx = layerIdx; - break; - } - } - - if ( lowestKnown_zIdx == nlayers-1 ) { - // is a column of no data values, skip it for the log profile part of the fill - // also warn because the method for filling no data values past the log profile part hasn't yet been implemented - std::cout << "!!! no lowest known zIdx for column of data !!! for rowIdx = " << rowIdx << ", colIdx = " << colIdx << std::endl; - continue; - } - - - double current_z_ground = z(rowIdx,colIdx,0); - - double firstCellHeight_zIdx = 0; - for(int layerIdx=0; layerIdx= finalFirstCellHeight ) - { - firstCellHeight_zIdx = layerIdx; - break; - } - } - - - // for debugging purposes, want to know how different the indices and resulting values ended up - // turns out to be easier if change back and forth between when/why it is printed rather than printing each instance - //if ( lowestKnown_zIdx == firstCellHeight_zIdx ) { - //if ( lowestKnown_zIdx > firstCellHeight_zIdx ) { - // std::cout << "rowIdx = " << rowIdx << ", colIdx = " << colIdx << std::endl; - // std::cout << "lowestKnown_zIdx = " << lowestKnown_zIdx << ", firstCellHeight_zIdx = " << firstCellHeight_zIdx << std::endl; - // std::cout << "z_ASL[lowestKnown_zIdx] = " << z(rowIdx,colIdx,lowestKnown_zIdx) << ", current_z_ground = " << current_z_ground << ", z_AGL[lowestKnown_zIdx] = " << z(rowIdx,colIdx,lowestKnown_zIdx)-current_z_ground << ", firstCellHeight z = " << finalFirstCellHeight << std::endl; - //} - - - // this one is NOT for debugging purposes, has to be run each and every time to make the code run safely - if ( lowestKnown_zIdx < firstCellHeight_zIdx ) { - //std::cout << "lowestKnown zIdx is less than firstCellHeight zIdx for rowIdx = " << rowIdx << ", colIdx = " << colIdx << std::endl; - //std::cout << "old lowestKnown_zIdx = " << lowestKnown_zIdx << std::endl; - lowestKnown_zIdx = firstCellHeight_zIdx; - //std::cout << "new lowestKnown_zIdx = " << lowestKnown_zIdx << std::endl; - //std::cout << "new z_ASL[lowestKnown_zIdx] = " << z(rowIdx,colIdx,lowestKnown_zIdx) << ", current_z_ground = " << current_z_ground << ", new z_AGL[lowestKnown_zIdx] = " << z(rowIdx,colIdx,lowestKnown_zIdx)-current_z_ground << ", firstCellHeight z = " << finalFirstCellHeight << std::endl; - // looks like this is getting triggered a LOT, if not for each and every occurence, hard to tell without commenting it out and printing for the other occurences. - // looks like the resulting value is rarely equal to firstCellHeight though, seems to always be just a hint above that value - } - - - - double z_ref = z(rowIdx,colIdx,lowestKnown_zIdx); - double u_ref = u(rowIdx,colIdx,lowestKnown_zIdx); - double v_ref = v(rowIdx,colIdx,lowestKnown_zIdx); - double w_ref = w(rowIdx,colIdx,lowestKnown_zIdx); - - - profile.ObukovLength = init->L(rowIdx,colIdx); - profile.ABL_height = init->bl_height(rowIdx,colIdx); - profile.Roughness = input.surface.Roughness(rowIdx,colIdx); - double current_rough_h = input.surface.Rough_h(rowIdx,colIdx); - double current_rough_d = input.surface.Rough_d(rowIdx,colIdx); - profile.Rough_h = current_rough_h; - profile.Rough_d = current_rough_d; - - //this is height above the vegetation - profile.inputWindHeight = z_ref - current_z_ground; // needs to be AGL - - for(int zIdx=1; zIdxninjaCom(ninjaComClass::ninjaWarning, "Exception caught during volume VTK file writing: %s", e.what()); + } catch (...) { + input.Com->ninjaCom(ninjaComClass::ninjaWarning, "Exception caught during volume VTK file writing: Cannot determine exception type."); + } + } + void NinjaFoam::UpdateExistingCase() { if(!CheckForValidDem()) @@ -3742,7 +4117,7 @@ void NinjaFoam::SetMeshResolutionAndResampleDem() } - if ( writeMassMeshVtk == true ) { + if ( writeMassMesh == true ) { // need to setup mesh sizing BEFORE the dem gets resampled, but AFTER the mesh resolution gets set massMesh.set_numVertLayers(20); // done in cli.cpp calling ninja_army calling ninja calling this function, with windsim.setNumVertLayers( i_, 20); where i_ is ninjaIdx CPLDebug("NINJAFOAM", "mass mesh vtk output set by mesh resolution, %f %s", meshResolution, lengthUnits::getString(meshResolutionUnits).c_str()); diff --git a/src/ninja/ninjafoam.h b/src/ninja/ninjafoam.h index 2aba8b1d..92ae4bd9 100644 --- a/src/ninja/ninjafoam.h +++ b/src/ninja/ninjafoam.h @@ -81,7 +81,6 @@ class NinjaFoam : public ninja static int GenerateFoamDirectory(std::string demName); static void SetFoamPath(const char *pszPath); - AsciiGrid TurbulenceGrid; private: std::string foamVersion; @@ -173,10 +172,13 @@ class NinjaFoam : public ninja void SetOutputResolution(); void SetOutputFilenames(); bool CheckIfOutputWindHeightIsResolved(); - - bool writeMassMeshVtk; + + bool writeMassMesh; Mesh massMesh; - void writeMassMeshVtkOutput(); + wn_3dScalarField massMesh_u, massMesh_v, massMesh_w; + wn_3dScalarField massMesh_k; + void GenerateAndSampleMassMesh(); + void generateMassMesh(); void writeProbeSampleFile(const wn_3dArray& x, const wn_3dArray& y, const wn_3dArray& z, const double dem_xllCorner, const double dem_yllCorner, const int ncols, const int nrows, const int nlayers); @@ -185,9 +187,25 @@ class NinjaFoam : public ninja const double dem_xllCorner, const double dem_yllCorner, const int ncols, const int nrows, const int nlayers, wn_3dScalarField& u, wn_3dScalarField& v, wn_3dScalarField& w); + void readInProbeData(const wn_3dArray& x, const wn_3dArray& y, const wn_3dArray& z, + const double dem_xllCorner, const double dem_yllCorner, + const int ncols, const int nrows, const int nlayers, + wn_3dScalarField& k); void fillEmptyProbeVals(const wn_3dArray& z, const int ncols, const int nrows, const int nlayers, wn_3dScalarField& u, wn_3dScalarField& v, wn_3dScalarField& w); + void fillEmptyProbeVals(const wn_3dArray& z, + const int ncols, const int nrows, const int nlayers, + wn_3dScalarField& k); + + void generateColMaxGrid(const wn_3dArray& z, + const double dem_xllCorner, const double dem_yllCorner, + const int ncols, const int nrows, const int nlayers, + const double massMeshResolution, std::string prjString, + wn_3dScalarField& k); + + bool writeMassMeshVtk; + void writeMassMeshVtkOutput(); const char *pszVrtMem; const char *pszVrtMemTurbulence; diff --git a/src/ninja/srtmclient.cpp b/src/ninja/srtmclient.cpp index f6b7fd8d..97b6059b 100644 --- a/src/ninja/srtmclient.cpp +++ b/src/ninja/srtmclient.cpp @@ -153,8 +153,10 @@ SURF_FETCH_E SRTMClient::FetchBoundingBox( double *bbox, double resolution, return SURF_FETCH_E_BAD_INPUT; } + //make a temporary file to hold the result of the HTTP request (in lat/lon coordinates) + CPLDebug( "SRTM_CLIENT", "Writing temporary file to: %s", CPLFormFilename(CPLGetPath(filename), "NINJA_SRTM", ".tif")); VSILFILE *fout; - fout = VSIFOpenL( "NINJA_SRTM.tif", "wb" ); + fout = VSIFOpenL( CPLFormFilename(CPLGetPath(filename), "NINJA_SRTM", ".tif"), "wb" ); if( !fout ) { CPLError( CE_Failure, CPLE_AppDefined, @@ -170,7 +172,7 @@ SURF_FETCH_E SRTMClient::FetchBoundingBox( double *bbox, double resolution, * Warp to UTM coords *-----------------------------------------------------------------------------*/ CPLDebug( "SRTM_CLIENT", "Warping to UTM..." ); - GDALDatasetH hDS = (GDALDatasetH) GDALOpen("NINJA_SRTM.tif", GA_ReadOnly); + GDALDatasetH hDS = (GDALDatasetH) GDALOpen(CPLFormFilename(CPLGetPath(filename), "NINJA_SRTM", ".tif"), GA_ReadOnly); GDALDatasetH hUtmDS; bool rc = GDALWarpToUtm( filename, hDS, hUtmDS); if(rc != true) @@ -180,7 +182,7 @@ SURF_FETCH_E SRTMClient::FetchBoundingBox( double *bbox, double resolution, GDALClose(hDS); GDALClose(hUtmDS); - VSIUnlink("NINJA_SRTM.tif"); + VSIUnlink(CPLFormFilename(CPLGetPath(filename), "NINJA_SRTM", ".tif")); return SURF_FETCH_E_NONE; } diff --git a/src/ninja/wrf3dInitialization.cpp b/src/ninja/wrf3dInitialization.cpp index c02fd39f..93831f82 100644 --- a/src/ninja/wrf3dInitialization.cpp +++ b/src/ninja/wrf3dInitialization.cpp @@ -423,8 +423,9 @@ void wrf3dInitialization::set3dGrids( WindNinjaInputs &input, Mesh const& mesh ) std::string legendTitle = "tempGrid"; std::string legendUnits = "(?)"; bool writeLegend = TRUE; + bool keepTiff = FALSE; - //tempGrid.ascii2png( outFilename, scalarLegendFilename, legendUnits, legendTitle, writeLegend ); + //tempGrid.ascii2png( outFilename, scalarLegendFilename, legendUnits, legendTitle, writeLegend, keepTiff ); //Make final grids with same header as dem @@ -437,8 +438,9 @@ void wrf3dInitialization::set3dGrids( WindNinjaInputs &input, Mesh const& mesh ) legendTitle = "temp2Grid"; legendUnits = "(?)"; writeLegend = TRUE; + keepTiff = FALSE; - //temp2Grid.ascii2png( outFilename, scalarLegendFilename, legendUnits, legendTitle, writeLegend );*/ + //temp2Grid.ascii2png( outFilename, scalarLegendFilename, legendUnits, legendTitle, writeLegend, keepTiff );*/ //=======end testing=================================// @@ -743,7 +745,7 @@ void wrf3dInitialization::set3dGrids( WindNinjaInputs &input, Mesh const& mesh ) } testGrid.write_Grid(filename.c_str(), 2); //if(k = 10){ - //testGrid.ascii2png( outFilename, scalarLegendFilename, legendUnits, legendTitle, writeLegend ); + //testGrid.ascii2png( outFilename, scalarLegendFilename, legendUnits, legendTitle, writeLegend, keepTiff ); //} } testGrid.deallocate();*/ @@ -1118,7 +1120,7 @@ void wrf3dInitialization::buildWxMeshes(WindNinjaInputs &input, Mesh const& mesh /*if(k == 0){ testGrid.write_Grid(filename.c_str(), 2); testGrid.replaceNan( -9999.0 ); - //testGrid.ascii2png( outFilename, scalarLegendFilename, legendUnits, legendTitle, writeLegend ); + //testGrid.ascii2png( outFilename, scalarLegendFilename, legendUnits, legendTitle, writeLegend, keepTiff ); }*/ //} //testGrid.deallocate(); @@ -1136,10 +1138,10 @@ void wrf3dInitialization::buildWxMeshes(WindNinjaInputs &input, Mesh const& mesh z = 0; // ----- x-staggered mesh ------------------------------------------------------------------------- - for(int i = 0; i < zStaggerWxMesh.NUMEL; i++) //Start loop over elements - { - elemElevation.node0 = zStaggerWxMesh.get_node0(i); //get the global node number of local node 0 of element i - for(int j = 0; j < elemElevation.NUMQPTV; j++) //Start loop over quadrature points in the element + for(int i = 0; i < zStaggerWxMesh.NUMEL; i++) //Start loop over elements + { + elemElevation.node0 = zStaggerWxMesh.get_node0(i); //get the global node number of local node 0 of element i + for(int j = 0; j < elemElevation.NUMQPTV; j++) //Start loop over quadrature points in the element { zStaggerWxMesh.get_elemIndex(i, elem_i, elem_j, elem_k); @@ -1194,10 +1196,10 @@ void wrf3dInitialization::buildWxMeshes(WindNinjaInputs &input, Mesh const& mesh } // ----- y-staggered mesh ------------------------------------------------------------------------- - for(int i = 0; i < zStaggerWxMesh.NUMEL; i++) //Start loop over elements - { - elemElevation.node0 = zStaggerWxMesh.get_node0(i); //get the global node number of local node 0 of element i - for(int j = 0; j < elemElevation.NUMQPTV; j++) //Start loop over quadrature points in the element + for(int i = 0; i < zStaggerWxMesh.NUMEL; i++) //Start loop over elements + { + elemElevation.node0 = zStaggerWxMesh.get_node0(i); //get the global node number of local node 0 of element i + for(int j = 0; j < elemElevation.NUMQPTV; j++) //Start loop over quadrature points in the element { zStaggerWxMesh.get_elemIndex(i, elem_i, elem_j, elem_k);