From b3ea15845ec41caf8fcec4e5414a3e087fc49dde Mon Sep 17 00:00:00 2001 From: latwood Date: Wed, 7 Aug 2024 13:16:26 -0600 Subject: [PATCH 1/5] added standalone wrf_to_kmz script, for reading in wrf netcdf files and outputting the corresponding kmz files, without having to run WindNinja itself each and every time --- CMakeLists.txt | 2 + src/CMakeLists.txt | 3 + src/wrf_to_kmz/CMakeLists.txt | 50 ++ src/wrf_to_kmz/wrf_to_kmz.cpp | 1101 +++++++++++++++++++++++++++++++++ 4 files changed, 1156 insertions(+) create mode 100644 src/wrf_to_kmz/CMakeLists.txt create mode 100644 src/wrf_to_kmz/wrf_to_kmz.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 7c76fa87..75823c0a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -317,6 +317,8 @@ endif(NOT WIN32) option(BUILD_FETCH_DEM "Build a standalone command line interface DEM utility" OFF) option(BUILD_STL_CONVERTER "Build a standalone command line interface for STL file conversions" OFF ) option(BUILD_CONVERT_OUTPUT "Build a standalone command line interface for xyz file conversions" OFF ) +option(BUILD_WRF_TO_KMZ "Build a standalone command line interface for converting wrf output runs to kmz, without running full WindNinja" OFF ) +mark_as_advanced(BUILD_WRF_TO_KMZ) option(BUILD_SLOPE_ASPECT_GRID "Build an application for building slope and aspect grids from a dem" OFF) mark_as_advanced(BUILD_SLOPE_ASPECT_GRID) option(BUILD_SOLAR_GRID "Build a application for building solar grids" OFF) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 80696d83..74e6a609 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -40,3 +40,6 @@ endif(BUILD_STL_CONVERTER) if(BUILD_CONVERT_OUTPUT) add_subdirectory(output_converter) endif(BUILD_CONVERT_OUTPUT) +if(BUILD_WRF_TO_KMZ) + add_subdirectory(wrf_to_kmz) +endif(BUILD_WRF_TO_KMZ) diff --git a/src/wrf_to_kmz/CMakeLists.txt b/src/wrf_to_kmz/CMakeLists.txt new file mode 100644 index 00000000..da0687f4 --- /dev/null +++ b/src/wrf_to_kmz/CMakeLists.txt @@ -0,0 +1,50 @@ +# THIS SOFTWARE WAS DEVELOPED AT THE ROCKY MOUNTAIN RESEARCH STATION (RMRS) +# MISSOULA FIRE SCIENCES LABORATORY BY EMPLOYEES OF THE FEDERAL GOVERNMENT +# IN THE COURSE OF THEIR OFFICIAL DUTIES. PURSUANT TO TITLE 17 SECTION 105 +# OF THE UNITED STATES CODE, THIS SOFTWARE IS NOT SUBJECT TO COPYRIGHT +# PROTECTION AND IS IN THE PUBLIC DOMAIN. RMRS MISSOULA FIRE SCIENCES +# LABORATORY ASSUMES NO RESPONSIBILITY WHATSOEVER FOR ITS USE BY OTHER +# PARTIES, AND MAKES NO GUARANTEES, EXPRESSED OR IMPLIED, ABOUT ITS QUALITY, +# RELIABILITY, OR ANY OTHER CHARACTERISTIC. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. + +cmake_minimum_required(VERSION 2.6) + +include_directories(${PROJECT_SOURCE_DIR}/src + ${PROJECT_SOURCE_DIR}/src/ninja + ${Boost_INCLUDE_DIRS} + ${NETCDF_INCLUDES} + ${GDAL_SYSTEM_INCLUDE} ${GDAL_INCLUDE_DIR}) + +set(LINK_LIBS ${NETCDF_LIBRARIES_C} + ${GDAL_LIBRARY} + ${Boost_LIBRARIES}) + +set(WRF_TO_KMZ_SRC wrf_to_kmz.cpp + + ${PROJECT_SOURCE_DIR}/src/ninja/Array2D.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ascii_grid.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/EasyBMP.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/EasyBMP_Geometry.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/EasyBMP_Font.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ninja_conv.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ninja_init.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/gdal_util.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/KmlVector.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/Style.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/LineStyle.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ninjaMathUtility.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ninjaUnits.cpp) + +add_executable(wrf_to_kmz ${WRF_TO_KMZ_SRC}) +target_link_libraries(wrf_to_kmz ${LINK_LIBS}) + +install(TARGETS wrf_to_kmz DESTINATION bin COMPONENT apps) + diff --git a/src/wrf_to_kmz/wrf_to_kmz.cpp b/src/wrf_to_kmz/wrf_to_kmz.cpp new file mode 100644 index 00000000..5a5137e6 --- /dev/null +++ b/src/wrf_to_kmz/wrf_to_kmz.cpp @@ -0,0 +1,1101 @@ +/****************************************************************************** + * + * $Id$ + * + * Project: WindNinja + * Purpose: Executable for converting xyz output from OpenFOAM + * Author: Natalie Wagenbrenner + * + ****************************************************************************** + * + * THIS SOFTWARE WAS DEVELOPED AT THE ROCKY MOUNTAIN RESEARCH STATION (RMRS) + * MISSOULA FIRE SCIENCES LABORATORY BY EMPLOYEES OF THE FEDERAL GOVERNMENT + * IN THE COURSE OF THEIR OFFICIAL DUTIES. PURSUANT TO TITLE 17 SECTION 105 + * OF THE UNITED STATES CODE, THIS SOFTWARE IS NOT SUBJECT TO COPYRIGHT + * PROTECTION AND IS IN THE PUBLIC DOMAIN. RMRS MISSOULA FIRE SCIENCES + * LABORATORY ASSUMES NO RESPONSIBILITY WHATSOEVER FOR ITS USE BY OTHER + * PARTIES, AND MAKES NO GUARANTEES, EXPRESSED OR IMPLIED, ABOUT ITS QUALITY, + * RELIABILITY, OR ANY OTHER CHARACTERISTIC. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + * + *****************************************************************************/ + + +#include "ninja_init.h" +#include "ninja_conv.h" + +#include "ascii_grid.h" + +#include "netcdf.h" +#include "gdal_util.h" + +#include "ninjaUnits.h" +#include "KmlVector.h" + + + +/** +* Static identifier to determine if the netcdf file is a WRF forecast. +* Uses netcdf c api +* @param fileName netcdf filename +* @return true if the forecast is a WRF forecast +*/ +bool identify( std::string fileName ) +{ + bool identified = true; + + //Acquire a lock to protect the non-thread safe netCDF library +//#ifdef _OPENMP +// omp_guard netCDF_guard(netCDF_lock); +//#endif + + /* + * Open the dataset + */ + + int status, ncid, ndims, nvars, ngatts, unlimdimid; + status = nc_open( fileName.c_str(), 0, &ncid ); + if ( status != NC_NOERR ) { + identified = false; + } + + /* + * Check the global attributes for the following tag: + * :TITLE = "...WRF..."" + */ + size_t len; + nc_type type; + char* model; + status = nc_inq_att( ncid, NC_GLOBAL, "TITLE", &type, &len ); + if( status != NC_NOERR ) + identified = false; + else { + model = new char[len + 1]; + status = nc_get_att_text( ncid, NC_GLOBAL, "TITLE", model ); + model[len] = '\0'; + std::string s( model ); + if( s.find("WRF") == s.npos ) { + identified = false; + } + + delete[] model; + } + + status = nc_close( ncid ); + + return identified; +} + + +/** +* Fetch the variable names +* @return a vector of variable names +*/ +std::vector getVariableList() +{ + std::vector varList; + varList.push_back( "T2" ); // T at 2 m + varList.push_back( "V10" ); // V at 10 m + varList.push_back( "U10" ); // U at 10 m + varList.push_back( "QCLOUD" ); // cloud water mixing ratio + return varList; +} + + +/** +* Checks the downloaded data to see if it is all valid. +* May not be functional yet for this class... +*/ +void checkForValidData( std::string wxModelFileName ) +{ + // open ds variable by variable + GDALDataset *srcDS; + std::string temp; + std::string srcWkt; + int nBands = 0; + bool noDataValueExists; + bool noDataIsNan; + + std::vector varList = getVariableList(); + + //Acquire a lock to protect the non-thread safe netCDF library +//#ifdef _OPENMP +// omp_guard netCDF_guard(netCDF_lock); +//#endif + + for( unsigned int i = 0;i < varList.size();i++ ) { + + temp = "NETCDF:\"" + wxModelFileName + "\":" + varList[i]; + //std::cout << "temp = \"" << temp.c_str() << "\"" << std::endl; + + CPLPushErrorHandler(&CPLQuietErrorHandler); + srcDS = (GDALDataset*)GDALOpen( temp.c_str(), GA_ReadOnly ); + CPLPopErrorHandler(); + if( srcDS == NULL ) + throw badForecastFile("Cannot open forecast file."); + + //Get total bands (time steps) + nBands = srcDS->GetRasterCount(); + //std::cout << "nBands = " << nBands << std::endl; + nBands = srcDS->GetRasterCount(); + int nXSize, nYSize; + GDALRasterBand *poBand; + int pbSuccess; + double dfNoData; + double *padfScanline; + + nXSize = srcDS->GetRasterXSize(); + nYSize = srcDS->GetRasterYSize(); + + //std::cout << "nXsize = " << nXSize << std::endl; + //std::cout << "nYsize = " << nYSize << std::endl; + + //loop over all bands for this variable (bands are time steps) + for(int j = 1; j <= nBands; j++) + { + poBand = srcDS->GetRasterBand( j ); + + pbSuccess = 0; + dfNoData = poBand->GetNoDataValue( &pbSuccess ); + if( pbSuccess == false ) + noDataValueExists = false; + else + { + noDataValueExists = true; + noDataIsNan = CPLIsNan(dfNoData); + } + + //set the data + padfScanline = new double[nXSize*nYSize]; + poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, padfScanline, nXSize, nYSize, + GDT_Float64, 0, 0); + for(int k = 0;k < nXSize*nYSize; k++) + { + //Check if value is no data (if no data value was defined in file) + if(noDataValueExists) + { + if(noDataIsNan) + { + if(CPLIsNan(padfScanline[k])) + throw badForecastFile("Forecast file contains no_data values."); + }else + { + if(padfScanline[k] == dfNoData) + throw badForecastFile("Forecast file contains no_data values."); + } + } + if( varList[i] == "T2" ) //units are Kelvin + { + if(padfScanline[k] < 180.0 || padfScanline[k] > 340.0) //these are near the most extreme temperatures ever recored on earth + throw badForecastFile("Temperature is out of range in forecast file."); + } + else if( varList[i] == "V10" ) //units are m/s + { + if(std::abs(padfScanline[k]) > 220.0) + throw badForecastFile("V-velocity is out of range in forecast file."); + } + else if( varList[i] == "U10" ) //units are m/s + { + if(std::abs(padfScanline[k]) > 220.0) + throw badForecastFile("U-velocity is out of range in forecast file."); + } + else if( varList[i] == "QCLOUD" ) //units are kg/kg + { + if(padfScanline[k] < -0.0001 || padfScanline[k] > 100.0) + throw badForecastFile("Total cloud cover is out of range in forecast file."); + } + } + + delete [] padfScanline; + } + + GDALClose((GDALDatasetH) srcDS ); + } +} + + +//void getNcGlobalAttributes( const std::string &wxModelFileName, int &mapProj, float &dx, float &dy, float &cenLat, float &cenLon, float &moadCenLat, float &standLon, float &trueLat1, float &trueLat2, int &wxModel_nLayers, std::string &projString ) +void getNcGlobalAttributes( const std::string &wxModelFileName, float &dx, float &dy, float &cenLat, float &cenLon, std::string &projString ) +{ + + //==========get global attributes to set projection=========================== + //Acquire a lock to protect the non-thread safe netCDF library +//#ifdef _OPENMP +// omp_guard netCDF_guard(netCDF_lock); +//#endif + + /* + * Open the dataset + */ + int status, ncid; + status = nc_open( wxModelFileName.c_str(), 0, &ncid ); + if ( status != NC_NOERR ) { + ostringstream os; + os << "The netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + + /* + * Get global attribute MAP_PROJ + * 1 = Lambert Conformal Conic + * 2 = Polar Stereographic + * 3 = Mercator + * 6 = Lat/Long + */ + + int mapProj; + nc_type type; + size_t len; + status = nc_inq_att( ncid, NC_GLOBAL, "MAP_PROJ", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute 'MAP_PROJ' in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_int( ncid, NC_GLOBAL, "MAP_PROJ", &mapProj ); + } + + //cout<<"MAP_PROJ = "<(standLon)+"],\ + PARAMETER[\"latitude_of_origin\","+boost::lexical_cast(moadCenLat)+"],\ + PARAMETER[\"standard_parallel_1\","+boost::lexical_cast(trueLat1)+"],\ + PARAMETER[\"false_easting\",0.0],PARAMETER[\"false_northing\",0.0],\ + PARAMETER[\"standard_parallel_2\","+boost::lexical_cast(trueLat2)+"],\ + UNIT[\"m\",1.0],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"; + } + else if(mapProj == 2){ //polar stereographic + projString = "PROJCS[\"WGS 84 / Antarctic Polar Stereographic\",GEOGCS[\"WGS 84\",\ + DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],\ + AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],\ + UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],\ + AUTHORITY[\"EPSG\",\"4326\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],\ + PROJECTION[\"Polar_Stereographic\"],\ + PARAMETER[\"latitude_of_origin\","+boost::lexical_cast(moadCenLat)+"],\ + PARAMETER[\"central_meridian\","+boost::lexical_cast(standLon)+"],\ + PARAMETER[\"scale_factor\",1],\ + PARAMETER[\"false_easting\",0],\ + PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"3031\"],\ + AXIS[\"Easting\",UNKNOWN],AXIS[\"Northing\",UNKNOWN]]"; + } + else if(mapProj == 3){ //mercator + projString = "PROJCS[\"World_Mercator\",GEOGCS[\"GCS_WGS_1984\",DATUM[\"WGS_1984\",\ + SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],\ + UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Mercator_1SP\"],\ + PARAMETER[\"False_Easting\",0],\ + PARAMETER[\"False_Northing\",0],\ + PARAMETER[\"Central_Meridian\","+boost::lexical_cast(standLon)+"],\ + PARAMETER[\"latitude_of_origin\","+boost::lexical_cast(moadCenLat)+"],\ + UNIT[\"Meter\",1]]"; + } + else if(mapProj == 6){ //lat/long + throw badForecastFile("Cannot initialize with a forecast file in lat/long spacing. \ + Forecast file must be in a projected coordinate system."); + } + else throw badForecastFile("Cannot determine projection from the forecast file information."); + +} + + +//std::string getTimeZoneString( const std::string &wxModelFileName ) +//std::string getTimeZoneString( const double &cenLat, const double &cenLon, const std::string &projString ) +std::string getTimeZoneString( const double &cenLat, const double &cenLon ) +{ + std::string timeZoneString; + + /* + //// old way, from cli.cpp, works on a dem, but does not work on the wrf file + double longitude = 0; + double latitude = 0; + GDALDataset *poDS = (GDALDataset*)GDALOpen(wxModelFileName.c_str(), GA_ReadOnly); + if(poDS == NULL) + { + GDALClose((GDALDatasetH)poDS); + fprintf(stderr, "Unable to open input wrf file to get timezone string!!!\n"); + std::exit(1); + } + GDALGetCenter(poDS, &longitude, &latitude); + GDALClose((GDALDatasetH)poDS); + */ + + double longitude = cenLon; + double latitude = cenLat; + + std::string tz = FetchTimeZone(longitude, latitude, NULL); + if(tz == "") + { + fprintf(stderr, "Could not detect timezone from input wrf file!!!\n"); + std::exit(1); + } + else{ + timeZoneString = tz; + } + + return timeZoneString; +} + + +/** + * Fetch the list of times that the wrf file holds. It is assumed that + * the time variable is "Times" and the units string is "units". If this + * is not the case, this function needs to be overridden. + * Uses the netcdf api, so we need omp critical sections. + * + * @param timeZoneString Time zone name from the file "date_time_zonespec.csv". + * @param wxModelFileName the input wrf file to open and read times from. + * @throw runtime_error for bad file i/o + * @return a vector of boost::posix_time::ptime objects for the forecast + */ +std::vector +getTimeList( const std::string &timeZoneString, const std::string &wxModelFileName ) +{ + + const char *pszVariable = NULL; + + boost::local_time::time_zone_ptr timeZonePtr; + timeZonePtr = globalTimeZoneDB.time_zone_from_region(timeZoneString.c_str()); + if( timeZonePtr == NULL ) { + ostringstream os; + os << "The time zone string: " << timeZoneString.c_str() + << " does not match any in " + << "the time zone database file: date_time_zonespec.csv."; + throw std::runtime_error( os.str() ); + } + + std::vectortimeList; + +// //Acquire a lock to protect the non-thread safe netCDF library +//#ifdef _OPENMP +// omp_guard netCDF_guard(netCDF_lock); +//#endif + + int status, ncid, ndims, nvars, ngatts, unlimdimid; + nc_type vartype; + int varndims, varnatts; + double *varvals; + int vardimids[NC_MAX_VAR_DIMS]; /* dimension IDs */ + + /* + * Open the dataset + */ + status = nc_open( wxModelFileName.c_str(), 0, &ncid ); + if ( status != NC_NOERR ) { + ostringstream os; + os << "The netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + + + // note that "Times" is usually set per dataset type, is usually "time" for many of the other files + std::string timename = "Times"; + + + /* + * If we can't get simple data from the file, return false + */ + status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid); + if ( status != NC_NOERR ) { + ostringstream os; + os << "The variables in the netcdf file: " << wxModelFileName + << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + /* + * Check the variable names, return false if any aren't found + */ + int varid; + status = nc_inq_varid( ncid, timename.c_str(), &varid ); + if( status != NC_NOERR ) { + ostringstream os; + //os << "The variable \"time\" in the netcdf file: " + os << "The variable \"Times\" in the netcdf file: " + << wxModelFileName << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + //==============If the forecast is WRF, parse Times============================ + + /* + * Check to see if we can read Times variable + */ + status = nc_inq_var( ncid, varid, 0, &vartype, &varndims, vardimids, + &varnatts ); + if( status != NC_NOERR ) { + ostringstream os; + os << "The values for variable \"Times\" " + << " in the netcdf file: " << wxModelFileName + << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + + /* + * Get varid for U10 --> use this to get length of time dimension + */ + status = nc_inq_varid( ncid, "U10", &varid ); + if( status != NC_NOERR ) { + ostringstream os; + os << "The variable \"U10\" in the netcdf file: " + << wxModelFileName << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + + /* + * Get dimid for Time in U10 + */ + int dimid; + status = nc_inq_dimid(ncid, "Time", &dimid ); + if( status != NC_NOERR ) { + ostringstream os; + os << "The dimension \"Time\" " + << " in the netcdf file: " << wxModelFileName + << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + + /* + * Get length of the time dimension in U10 + */ + size_t time_len; + status = nc_inq_dimlen( ncid, dimid, &time_len ); + if( status != NC_NOERR ) { + ostringstream os; + os << "The time dimension for variable \"Times\" in the netcdf file: " + << wxModelFileName << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + /* + * Reset varid to 'Times' + */ + //int varid; + status = nc_inq_varid( ncid, timename.c_str(), &varid ); + if( status != NC_NOERR ) { + ostringstream os; + os << "The variable \"U10\" in the netcdf file: " + << wxModelFileName << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + /* + * Get dimid for DateStrLen + */ + //int dimid; + status = nc_inq_dimid(ncid, "DateStrLen", &dimid ); + //cout<<"dimid =" <(i)}; /* where to get value from */ + status = nc_get_var1_text( ncid, varid, varindex, tp+i ); + } + if( status != NC_NOERR ) { + ostringstream os; + os << "The values for variable \"time\" " + << " in the netcdf file: " << wxModelFileName + << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + tp[t_len] = '\0'; + std::string refString( tp ); + delete[] tp; + + /* + * Clean the string for the boost constructor. Remove the prefix, + * suffix, and delimiters + */ + int pos = -1; + if( refString.find( "_" ) != refString.npos ) { + pos = refString.find( "_" ); + refString.replace(pos, 1, 1, 'T'); + //refString.erase( pos, 1 ); + } + while( refString.find( "-" ) != refString.npos ) { + pos = refString.find( "-" ); + refString.erase( pos, 1 ); + } + while( refString.find( ":" ) != refString.npos ) { + pos = refString.find( ":" ); + refString.erase( pos, 1 ); + } + pos = refString.find("Hour since "); + if(pos != refString.npos) + { + refString.erase(pos, strlen("Hour since ")); + } + + /* + * Make a posix time in UTC/GMT + */ + bpt::ptime reference_pt; + reference_pt = bpt::from_iso_string( refString ); + bpt::ptime first_pt( reference_pt ); + blt::local_date_time first_pt_local(first_pt, timeZonePtr); + timeList.push_back( first_pt_local ); + + } // end for loop + + + return timeList; +} + + +/** +* Sets the surface grids based on a WRF (surface only!) forecast. +* @param wxModelFileName The wrf input filename from which data are read from. +* @param timeBandIdx The band index from which data are read from, corresponding to a specific time from the getTimeList() +* @param airGrid The air temperature grid to be filled. +* @param cloudGrid The cloud cover grid to be filled. +* @param uGrid The u velocity grid to be filled. +* @param vGrid The v velocity grid to be filled. +* @param wGrid The w velocity grid to be filled (filled with zeros here?). +*/ +void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx, const float &dx, const float &dy, const float &cenLat, const float &cenLon, const std::string &projString, AsciiGrid &airGrid, AsciiGrid &cloudGrid, AsciiGrid &uGrid, AsciiGrid &vGrid, AsciiGrid &wGrid ) +{ + + // looks like the bands go in the same order as the timeList + int bandNum = timeBandIdx+1; // timeIdx is 0 to N-1, bandNum is 1 to N. + + + GDALDataset* poDS; + CPLPushErrorHandler(&CPLQuietErrorHandler); + poDS = (GDALDataset*)GDALOpenShared( wxModelFileName.c_str(), GA_ReadOnly ); + CPLPopErrorHandler(); + if( poDS == NULL ) { + CPLDebug( "wrfInitialization::setSurfaceGrids()", + "Bad forecast file"); + throw badForecastFile("Cannot open forecast file in wrfInitialization::setSurfaceGrids()"); + } + else { + GDALClose((GDALDatasetH) poDS ); // close original wxModel file + } + + + // open ds one by one, set geotranform, set projection, then write to grid + GDALDataset *srcDS; + std::string temp; + std::vector varList = getVariableList(); + + for( unsigned int i = 0;i < varList.size();i++ ) { + + temp = "NETCDF:\"" + wxModelFileName + "\":" + varList[i]; + + CPLPushErrorHandler(&CPLQuietErrorHandler); + srcDS = (GDALDataset*)GDALOpenShared( temp.c_str(), GA_ReadOnly ); + CPLPopErrorHandler(); + if( srcDS == NULL ) { + CPLDebug( "wrfInitialization::setSurfaceGrids()", + "Bad forecast file" ); + } + + CPLDebug("WX_MODEL_INITIALIZATION", "varList[i] = %s", varList[i].c_str()); + + /* + * Set up spatial reference stuff for setting projections + * and geotransformations + */ + + OGRSpatialReference oSRS, *poLatLong; + char *srcWKT = NULL; + char* prj2 = (char*)projString.c_str(); + oSRS.importFromWkt(&prj2); + oSRS.exportToWkt(&srcWKT); + + CPLDebug("WX_MODEL_INITIALIZATION", "srcWKT= %s", srcWKT); + + poLatLong = oSRS.CloneGeogCS(); + char *dstWkt = NULL; + poLatLong->exportToWkt(&dstWkt); + + /* + * Transform domain center from lat/long to WRF space + */ + double zCenter; + zCenter = 0; + double xCenter, yCenter; + xCenter = (double)cenLon; + yCenter = (double)cenLat; + +//#ifdef GDAL_COMPUTE_VERSION +//#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) + oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); +//#endif /* GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) */ +//#endif /* GDAL_COMPUTE_VERSION */ + + OGRCoordinateTransformation *poCT; + poCT = OGRCreateCoordinateTransformation(poLatLong, &oSRS); + delete poLatLong; + + if(poCT==NULL || !poCT->Transform(1, &xCenter, &yCenter)) + printf("Transformation failed.\n"); + + CPLDebug("WX_MODEL_INITIALIZATION", "xCenter, yCenter= %f, %f", xCenter, yCenter); + + /* + * Set the geostransform for the WRF file + * upper corner is calculated from transformed x, y + * (in WRF space) + */ + + double ncols, nrows; + int nXSize = srcDS->GetRasterXSize(); + int nYSize = srcDS->GetRasterYSize(); + ncols = (double)(nXSize)/2; + nrows = (double)(nYSize)/2; + + double adfGeoTransform[6] = {(xCenter-(ncols*dx)), dx, 0, + (yCenter+(nrows*dy)), + 0, -(dy)}; + + CPLDebug("WX_MODEL_INITIALIZATION", "ulcornerX, ulcornerY= %f, %f", (xCenter-(ncols*dx)), (yCenter+(nrows*dy))); + CPLDebug("WX_MODEL_INITIALIZATION", "nXSize, nYsize= %d, %d", nXSize, nYSize); + CPLDebug("WX_MODEL_INITIALIZATION", "dx= %f", dx); + + srcDS->SetGeoTransform(adfGeoTransform); + + /* + * get the noDataValue from the current band + */ + + GDALRasterBand *poBand = srcDS->GetRasterBand( bandNum ); + int pbSuccess; + double dfNoData = poBand->GetNoDataValue( &pbSuccess ); + + int nBandCount = srcDS->GetRasterCount(); + + CPLDebug("WX_MODEL_INITIALIZATION", "band count = %d", nBandCount); + + if( pbSuccess == false ) + dfNoData = -9999.0; + + /* + * set the dataset projection + */ + + int rc = srcDS->SetProjection( projString.c_str() ); + + /* + * final setting of the datasets to ascii grids, in the past usually done using a wrp dataset + */ + + CPLDebug("WX_MODEL_INITIALIZATION", "band number to write = %d", bandNum); + + if( varList[i] == "T2" ) { + GDAL2AsciiGrid( srcDS, bandNum, airGrid ); + if( CPLIsNan( dfNoData ) ) { + airGrid.set_noDataValue(-9999.0); + airGrid.replaceNan( -9999.0 ); + } + } + else if( varList[i] == "V10" ) { + GDAL2AsciiGrid( srcDS, bandNum, vGrid ); + if( CPLIsNan( dfNoData ) ) { + vGrid.set_noDataValue(-9999.0); + vGrid.replaceNan( -9999.0 ); + } + } + else if( varList[i] == "U10" ) { + GDAL2AsciiGrid( srcDS, bandNum, uGrid ); + if( CPLIsNan( dfNoData ) ) { + uGrid.set_noDataValue(-9999.0); + uGrid.replaceNan( -9999.0 ); + } + } + else if( varList[i] == "QCLOUD" ) { + GDAL2AsciiGrid( srcDS, bandNum, cloudGrid ); + if( CPLIsNan( dfNoData ) ) { + cloudGrid.set_noDataValue(-9999.0); + cloudGrid.replaceNan( -9999.0 ); + } + } + CPLFree(srcWKT); + CPLFree(dstWkt); + delete poCT; + GDALClose((GDALDatasetH) srcDS ); + } + //don't allow small negative values in cloud cover + for(int i=0; i &uGrid_wxModel, AsciiGrid &vGrid_wxModel ) +{ + +/* + //// wxModelInitialization grid setting style + //// causes the input ascii grids to need to drop the "const" in front, claims set_headerData() discards qualifiers + AsciiGrid speedInitializationGrid_wxModel; + AsciiGrid dirInitializationGrid_wxModel; + + speedInitializationGrid_wxModel.set_headerData(uGrid_wxModel); + dirInitializationGrid_wxModel.set_headerData(uGrid_wxModel); + + //now make speed and direction from u,v components + for(int i=0; i speedInitializationGrid_wxModel( uGrid_wxModel ); + AsciiGrid dirInitializationGrid_wxModel( uGrid_wxModel ); + + for(int i=0; iformat("%m-%d-%Y_%H%M"); + wxModelTimestream << forecastTime; + std::string forecastIdentifier = "WRF-SURFACE"; + std::string rootname = forecastIdentifier + "-" + wxModelTimestream.str(); + + + KmlVector ninjaKmlFiles; + + ninjaKmlFiles.setKmlFile( CPLFormFilename(path.c_str(), rootname.c_str(), "kml") ); + ninjaKmlFiles.setKmzFile( CPLFormFilename(path.c_str(), rootname.c_str(), "kmz") ); + ////ninjaKmlFiles.setDemFile(dem_filename); // turns out to be redundant and doesn't do anything, which is good because don't want this dependency + + ninjaKmlFiles.setLegendFile( CPLFormFilename(path.c_str(), rootname.c_str(), "bmp") ); + ninjaKmlFiles.setSpeedGrid(speedInitializationGrid_wxModel, outputSpeedUnits); + ninjaKmlFiles.setDirGrid(dirInitializationGrid_wxModel); + + //ninjaKmlFiles.setLineWidth(1.0); // input.googLineWidth value + ninjaKmlFiles.setLineWidth(3.0); // input.wxModelGoogLineWidth value + std::string dateTimewxModelLegFileTemp = CPLFormFilename(path.c_str(), (rootname+"_date_time").c_str(), "bmp"); + ninjaKmlFiles.setTime(forecastTime); + ninjaKmlFiles.setDateTimeLegendFile(dateTimewxModelLegFileTemp, forecastTime); + ninjaKmlFiles.setWxModel(forecastIdentifier, forecastTime); + + // default values for input.wxModelGoogSpeedScaling/input.googSpeedScaling,input.googColor,input.googVectorScale are KmlVector::equal_interval,"default",false respectively + if(ninjaKmlFiles.writeKml(KmlVector::equal_interval,"default",false)) + { + if(ninjaKmlFiles.makeKmz()) + ninjaKmlFiles.removeKmlFile(); + } +} + + +int main( int argc, char* argv[] ) +{ + /* parse input arguments */ + if( argc != 3 ) + { + std::cout << "Invalid arguments!" << std::endl; + std::cout << "wrf_to_kmz [input_wrf_filename] [output_speed_units]" << std::endl; + return 1; + } + std::string input_wrf_filename = std::string( argv[1] ); + std::string outputSpeedUnits_str = std::string( argv[2] ); + + std::cout << "input_wrf_filename = \"" << input_wrf_filename.c_str() << "\"" << std::endl; + std::cout << "output_speed_units = \"" << outputSpeedUnits_str.c_str() << "\"" << std::endl; + + + NinjaInitialize(); // needed for GDALAllRegister() + + + velocityUnits::eVelocityUnits outputSpeedUnits = velocityUnits::getUnit(outputSpeedUnits_str); + + + if ( identify( input_wrf_filename ) == false ) + { + throw badForecastFile("input input_wrf_filename is not a valid WRF file!!!"); + } + checkForValidData( input_wrf_filename ); + + + float dx; + float dy; + float cenLat; + float cenLon; + std::string projString; + getNcGlobalAttributes( input_wrf_filename, dx, dy, cenLat, cenLon, projString ); + + //std::string timeZoneString = getTimeZoneString( input_wrf_filename ); + //std::string timeZoneString = getTimeZoneString( cenLat, cenLon, projString ); + std::string timeZoneString = getTimeZoneString( cenLat, cenLon ); + + std::vector timeList = getTimeList( timeZoneString, input_wrf_filename ); + + for(unsigned int timeIdx = 0; timeIdx < timeList.size(); timeIdx++) + { + boost::local_time::local_date_time forecastTime = timeList[timeIdx]; + + AsciiGrid airGrid; + AsciiGrid cloudGrid; + AsciiGrid uGrid; + AsciiGrid vGrid; + AsciiGrid wGrid; + + setSurfaceGrids( input_wrf_filename, timeIdx, dx, dy, cenLat, cenLon, projString, airGrid, cloudGrid, uGrid, vGrid, wGrid ); + + writeWxModelGrids( input_wrf_filename, forecastTime, outputSpeedUnits, uGrid, vGrid ); + } + + + return 0; +} + + + + + + + + From 7fabc84824a23ffc602dbc42fbb8be6f11a29257 Mon Sep 17 00:00:00 2001 From: latwood Date: Fri, 16 Aug 2024 13:09:04 -0600 Subject: [PATCH 2/5] some more cleanup of th wrf_to_kmz script. Mostly tidying up some of the comments and aesthetics, cout and debug statements, but also a fix to the usage display of the possible units and an update to the comment header. Also, caught a bug where I forgot to add the line to convert the kmz output grids to the output units. --- src/wrf_to_kmz/wrf_to_kmz.cpp | 476 ++++++++++++---------------------- 1 file changed, 170 insertions(+), 306 deletions(-) diff --git a/src/wrf_to_kmz/wrf_to_kmz.cpp b/src/wrf_to_kmz/wrf_to_kmz.cpp index 5a5137e6..e85f9ff6 100644 --- a/src/wrf_to_kmz/wrf_to_kmz.cpp +++ b/src/wrf_to_kmz/wrf_to_kmz.cpp @@ -3,8 +3,8 @@ * $Id$ * * Project: WindNinja - * Purpose: Executable for converting xyz output from OpenFOAM - * Author: Natalie Wagenbrenner + * Purpose: Executable for converting wrf netcdf files to kmz without running WindNinja itself in a simulation + * Author: Loren Atwood * ****************************************************************************** * @@ -27,6 +27,8 @@ * *****************************************************************************/ +// got many of these functions from wxModelInitialization.cpp and wrfSurfInitialization.cpp + #include "ninja_init.h" #include "ninja_conv.h" @@ -43,33 +45,24 @@ /** * Static identifier to determine if the netcdf file is a WRF forecast. -* Uses netcdf c api -* @param fileName netcdf filename -* @return true if the forecast is a WRF forecast +* Uses netcdf c api. +* @param fileName netcdf filename. +* @return true if the forecast is a WRF forecast. */ bool identify( std::string fileName ) { bool identified = true; - //Acquire a lock to protect the non-thread safe netCDF library -//#ifdef _OPENMP -// omp_guard netCDF_guard(netCDF_lock); -//#endif - - /* - * Open the dataset - */ - + // Open the dataset int status, ncid, ndims, nvars, ngatts, unlimdimid; status = nc_open( fileName.c_str(), 0, &ncid ); if ( status != NC_NOERR ) { identified = false; } - /* - * Check the global attributes for the following tag: - * :TITLE = "...WRF..."" - */ + // Check the global attributes for the following tag: + // :TITLE = "...WRF..."" + size_t len; nc_type type; char* model; @@ -95,8 +88,8 @@ bool identify( std::string fileName ) /** -* Fetch the variable names -* @return a vector of variable names +* Fetch the variable names. +* @return a vector of variable names. */ std::vector getVariableList() { @@ -111,7 +104,8 @@ std::vector getVariableList() /** * Checks the downloaded data to see if it is all valid. -* May not be functional yet for this class... +* @param wxModelFileName wrf netcdf weather model filename. +* Comment From WindNinja code: May not be functional yet for this class... */ void checkForValidData( std::string wxModelFileName ) { @@ -125,11 +119,6 @@ void checkForValidData( std::string wxModelFileName ) std::vector varList = getVariableList(); - //Acquire a lock to protect the non-thread safe netCDF library -//#ifdef _OPENMP -// omp_guard netCDF_guard(netCDF_lock); -//#endif - for( unsigned int i = 0;i < varList.size();i++ ) { temp = "NETCDF:\"" + wxModelFileName + "\":" + varList[i]; @@ -141,7 +130,7 @@ void checkForValidData( std::string wxModelFileName ) if( srcDS == NULL ) throw badForecastFile("Cannot open forecast file."); - //Get total bands (time steps) + // Get total bands (time steps) nBands = srcDS->GetRasterCount(); //std::cout << "nBands = " << nBands << std::endl; nBands = srcDS->GetRasterCount(); @@ -157,7 +146,7 @@ void checkForValidData( std::string wxModelFileName ) //std::cout << "nXsize = " << nXSize << std::endl; //std::cout << "nYsize = " << nYSize << std::endl; - //loop over all bands for this variable (bands are time steps) + // loop over all bands for this variable (bands are time steps) for(int j = 1; j <= nBands; j++) { poBand = srcDS->GetRasterBand( j ); @@ -172,13 +161,13 @@ void checkForValidData( std::string wxModelFileName ) noDataIsNan = CPLIsNan(dfNoData); } - //set the data + // set the data padfScanline = new double[nXSize*nYSize]; poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, padfScanline, nXSize, nYSize, GDT_Float64, 0, 0); for(int k = 0;k < nXSize*nYSize; k++) { - //Check if value is no data (if no data value was defined in file) + // Check if value is no data (if no data value was defined in file) if(noDataValueExists) { if(noDataIsNan) @@ -191,22 +180,22 @@ void checkForValidData( std::string wxModelFileName ) throw badForecastFile("Forecast file contains no_data values."); } } - if( varList[i] == "T2" ) //units are Kelvin + if( varList[i] == "T2" ) // units are Kelvin { - if(padfScanline[k] < 180.0 || padfScanline[k] > 340.0) //these are near the most extreme temperatures ever recored on earth + if(padfScanline[k] < 180.0 || padfScanline[k] > 340.0) // these are near the most extreme temperatures ever recored on earth throw badForecastFile("Temperature is out of range in forecast file."); } - else if( varList[i] == "V10" ) //units are m/s + else if( varList[i] == "V10" ) // units are m/s { if(std::abs(padfScanline[k]) > 220.0) throw badForecastFile("V-velocity is out of range in forecast file."); } - else if( varList[i] == "U10" ) //units are m/s + else if( varList[i] == "U10" ) // units are m/s { if(std::abs(padfScanline[k]) > 220.0) throw badForecastFile("U-velocity is out of range in forecast file."); } - else if( varList[i] == "QCLOUD" ) //units are kg/kg + else if( varList[i] == "QCLOUD" ) // units are kg/kg { if(padfScanline[k] < -0.0001 || padfScanline[k] > 100.0) throw badForecastFile("Total cloud cover is out of range in forecast file."); @@ -221,19 +210,25 @@ void checkForValidData( std::string wxModelFileName ) } -//void getNcGlobalAttributes( const std::string &wxModelFileName, int &mapProj, float &dx, float &dy, float &cenLat, float &cenLon, float &moadCenLat, float &standLon, float &trueLat1, float &trueLat2, int &wxModel_nLayers, std::string &projString ) +/** +* Uses netcfd c api commands to get wrf netcdf file global attributes, to later use +* to set the gdal dataset projection and geotransform information, which are required +* to allow the wrf netcdf file data to be accessible by standard gdal commands. +* @param wxModelFileName wrf netcdf weather model filename, from which the global attributes data are read in and processed. +* @param dx The cell size x value of the dataset to be filled. +* @param dy The cell size y value of the dataset to be filled. +* @param cenLat The center latitude value of the dataset to be filled. +* @param cenLon The center longitude value of the dataset to be filled. +* @param projString The projection string of the dataset to be filled. +* note that there are additional wrf netcdf file global attributes that are accessed and used to get and process the projection string, +* but are not used outside this function, so they are not returned. These include mapProj, moadCenLat, standLon, trueLat1, trueLat2. +*/ void getNcGlobalAttributes( const std::string &wxModelFileName, float &dx, float &dy, float &cenLat, float &cenLon, std::string &projString ) { //==========get global attributes to set projection=========================== - //Acquire a lock to protect the non-thread safe netCDF library -//#ifdef _OPENMP -// omp_guard netCDF_guard(netCDF_lock); -//#endif - - /* - * Open the dataset - */ + + // Open the dataset int status, ncid; status = nc_open( wxModelFileName.c_str(), 0, &ncid ); if ( status != NC_NOERR ) { @@ -243,13 +238,11 @@ void getNcGlobalAttributes( const std::string &wxModelFileName, float &dx, float throw std::runtime_error( os.str() ); } - /* - * Get global attribute MAP_PROJ - * 1 = Lambert Conformal Conic - * 2 = Polar Stereographic - * 3 = Mercator - * 6 = Lat/Long - */ + // Get global attribute MAP_PROJ + // 1 = Lambert Conformal Conic + // 2 = Polar Stereographic + // 3 = Mercator + // 6 = Lat/Long int mapProj; nc_type type; @@ -264,14 +257,10 @@ void getNcGlobalAttributes( const std::string &wxModelFileName, float &dx, float else { status = nc_get_att_int( ncid, NC_GLOBAL, "MAP_PROJ", &mapProj ); } + //std::cout << "MAP_PROJ = " << mapProj << std::endl; - //cout<<"MAP_PROJ = "< getTimeList( const std::string &timeZoneString, const std::string &wxModelFileName ) @@ -551,20 +460,14 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa std::vectortimeList; -// //Acquire a lock to protect the non-thread safe netCDF library -//#ifdef _OPENMP -// omp_guard netCDF_guard(netCDF_lock); -//#endif int status, ncid, ndims, nvars, ngatts, unlimdimid; nc_type vartype; int varndims, varnatts; double *varvals; - int vardimids[NC_MAX_VAR_DIMS]; /* dimension IDs */ + int vardimids[NC_MAX_VAR_DIMS]; // dimension IDs - /* - * Open the dataset - */ + // Open the dataset status = nc_open( wxModelFileName.c_str(), 0, &ncid ); if ( status != NC_NOERR ) { ostringstream os; @@ -574,13 +477,11 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa } - // note that "Times" is usually set per dataset type, is usually "time" for many of the other files + // note that "Times" is the time variable name found in wrf files, this variable is usually set to "time" for many of the other weather model file types std::string timename = "Times"; - /* - * If we can't get simple data from the file, return false - */ + // If we can't get simple data from the file, return false status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid); if ( status != NC_NOERR ) { ostringstream os; @@ -589,9 +490,7 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - /* - * Check the variable names, return false if any aren't found - */ + // Check the variable names, return false if any aren't found int varid; status = nc_inq_varid( ncid, timename.c_str(), &varid ); if( status != NC_NOERR ) { @@ -602,11 +501,10 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - //==============If the forecast is WRF, parse Times============================ - /* - * Check to see if we can read Times variable - */ + //==============The forecast is WRF, parse Times============================ + + // Check to see if we can read Times variable status = nc_inq_var( ncid, varid, 0, &vartype, &varndims, vardimids, &varnatts ); if( status != NC_NOERR ) { @@ -618,9 +516,7 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa } - /* - * Get varid for U10 --> use this to get length of time dimension - */ + // Get varid for U10 --> use this to get length of time dimension status = nc_inq_varid( ncid, "U10", &varid ); if( status != NC_NOERR ) { ostringstream os; @@ -630,9 +526,7 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa } - /* - * Get dimid for Time in U10 - */ + // Get dimid for Time in U10 int dimid; status = nc_inq_dimid(ncid, "Time", &dimid ); if( status != NC_NOERR ) { @@ -644,9 +538,7 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa } - /* - * Get length of the time dimension in U10 - */ + // Get length of the time dimension in U10 size_t time_len; status = nc_inq_dimlen( ncid, dimid, &time_len ); if( status != NC_NOERR ) { @@ -656,10 +548,8 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - /* - * Reset varid to 'Times' - */ - //int varid; + // Reset varid to 'Times' + //int varid; // already set into the scope above status = nc_inq_varid( ncid, timename.c_str(), &varid ); if( status != NC_NOERR ) { ostringstream os; @@ -668,12 +558,10 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - /* - * Get dimid for DateStrLen - */ - //int dimid; + // Get dimid for DateStrLen + //int dimid; // already set into the scope above status = nc_inq_dimid(ncid, "DateStrLen", &dimid ); - //cout<<"dimid =" <(i)}; /* where to get value from */ + const size_t varindex[] = {t,static_cast(i)}; // where to get value from status = nc_get_var1_text( ncid, varid, varindex, tp+i ); } if( status != NC_NOERR ) { @@ -717,10 +601,7 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa std::string refString( tp ); delete[] tp; - /* - * Clean the string for the boost constructor. Remove the prefix, - * suffix, and delimiters - */ + // Clean the string for the boost constructor. Remove the prefix, suffix, and delimiters int pos = -1; if( refString.find( "_" ) != refString.npos ) { pos = refString.find( "_" ); @@ -741,9 +622,7 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa refString.erase(pos, strlen("Hour since ")); } - /* - * Make a posix time in UTC/GMT - */ + // Make a posix time in UTC/GMT, to call the boost::local_time::local_date_time input UTC time constructor bpt::ptime reference_pt; reference_pt = bpt::from_iso_string( refString ); bpt::ptime first_pt( reference_pt ); @@ -760,17 +639,21 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa /** * Sets the surface grids based on a WRF (surface only!) forecast. * @param wxModelFileName The wrf input filename from which data are read from. -* @param timeBandIdx The band index from which data are read from, corresponding to a specific time from the getTimeList() +* @param timeBandIdx The band index from which data are read from, corresponding to a specific expected time from getTimeList(). +* @param dx The cell size x value of the dataset, used to set the gdal dataset geotransform. +* @param dy The cell size y value of the dataset, used to set the gdal dataset geotransform. +* @param cenLat The center latitude value of the dataset, used to set the gdal dataset geotransform. +* @param cenLon The center longitude value of the dataset, used to set the gdal dataset geotransform. * @param airGrid The air temperature grid to be filled. * @param cloudGrid The cloud cover grid to be filled. * @param uGrid The u velocity grid to be filled. * @param vGrid The v velocity grid to be filled. -* @param wGrid The w velocity grid to be filled (filled with zeros here?). +* @param wGrid The w velocity grid to be filled (filled with zeros). */ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx, const float &dx, const float &dy, const float &cenLat, const float &cenLon, const std::string &projString, AsciiGrid &airGrid, AsciiGrid &cloudGrid, AsciiGrid &uGrid, AsciiGrid &vGrid, AsciiGrid &wGrid ) { - // looks like the bands go in the same order as the timeList + // looks like the bands go in the same order as the timeList, the past function looped through to find the idx that corresponded to an input time int bandNum = timeBandIdx+1; // timeIdx is 0 to N-1, bandNum is 1 to N. @@ -779,9 +662,7 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx poDS = (GDALDataset*)GDALOpenShared( wxModelFileName.c_str(), GA_ReadOnly ); CPLPopErrorHandler(); if( poDS == NULL ) { - CPLDebug( "wrfInitialization::setSurfaceGrids()", - "Bad forecast file"); - throw badForecastFile("Cannot open forecast file in wrfInitialization::setSurfaceGrids()"); + throw badForecastFile("Cannot open forecast file in setSurfaceGrids()"); } else { GDALClose((GDALDatasetH) poDS ); // close original wxModel file @@ -801,16 +682,13 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx srcDS = (GDALDataset*)GDALOpenShared( temp.c_str(), GA_ReadOnly ); CPLPopErrorHandler(); if( srcDS == NULL ) { - CPLDebug( "wrfInitialization::setSurfaceGrids()", - "Bad forecast file" ); + throw badForecastFile("Cannot open forecast file in setSurfaceGrids()"); } - CPLDebug("WX_MODEL_INITIALIZATION", "varList[i] = %s", varList[i].c_str()); + //std::cout << "varList[" << i << "] = " << varList[i].c_str() << std::endl; - /* - * Set up spatial reference stuff for setting projections - * and geotransformations - */ + // Set up spatial reference stuff for setting projections + // and geotransformations OGRSpatialReference oSRS, *poLatLong; char *srcWKT = NULL; @@ -818,41 +696,37 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx oSRS.importFromWkt(&prj2); oSRS.exportToWkt(&srcWKT); - CPLDebug("WX_MODEL_INITIALIZATION", "srcWKT= %s", srcWKT); + //std::cout << "srcWKT = " << srcWKT << std::endl; poLatLong = oSRS.CloneGeogCS(); char *dstWkt = NULL; poLatLong->exportToWkt(&dstWkt); - /* - * Transform domain center from lat/long to WRF space - */ + // Transform domain center from lat/long to WRF space double zCenter; zCenter = 0; double xCenter, yCenter; xCenter = (double)cenLon; yCenter = (double)cenLat; -//#ifdef GDAL_COMPUTE_VERSION -//#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) +#ifdef GDAL_COMPUTE_VERSION +#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); -//#endif /* GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) */ -//#endif /* GDAL_COMPUTE_VERSION */ +#endif /* GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) */ +#endif /* GDAL_COMPUTE_VERSION */ OGRCoordinateTransformation *poCT; poCT = OGRCreateCoordinateTransformation(poLatLong, &oSRS); delete poLatLong; if(poCT==NULL || !poCT->Transform(1, &xCenter, &yCenter)) - printf("Transformation failed.\n"); + throw std::runtime_error("Transformation of lat/lon center to dataset coordinate system failed!"); - CPLDebug("WX_MODEL_INITIALIZATION", "xCenter, yCenter= %f, %f", xCenter, yCenter); + //std::cout << "xCenter = " << xCenter << ", yCenter = " << yCenter << std::endl; - /* - * Set the geostransform for the WRF file - * upper corner is calculated from transformed x, y - * (in WRF space) - */ + // Set the geostransform for the WRF file + // upper corner is calculated from transformed x, y + // (in WRF space) double ncols, nrows; int nXSize = srcDS->GetRasterXSize(); @@ -864,15 +738,13 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx (yCenter+(nrows*dy)), 0, -(dy)}; - CPLDebug("WX_MODEL_INITIALIZATION", "ulcornerX, ulcornerY= %f, %f", (xCenter-(ncols*dx)), (yCenter+(nrows*dy))); - CPLDebug("WX_MODEL_INITIALIZATION", "nXSize, nYsize= %d, %d", nXSize, nYSize); - CPLDebug("WX_MODEL_INITIALIZATION", "dx= %f", dx); + //std::cout << "ulcornerX, ulcornerY = " << xCenter-(ncols*dx) << ", " << yCenter+(nrows*dy) << std::endl; + //std::cout << "nXSize, nYsize = " << nXSize << ", " << nYSize << std::endl; + //std::cout << "dx = " << dx << std::endl; srcDS->SetGeoTransform(adfGeoTransform); - /* - * get the noDataValue from the current band - */ + // get the noDataValue from the current band GDALRasterBand *poBand = srcDS->GetRasterBand( bandNum ); int pbSuccess; @@ -880,57 +752,55 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx int nBandCount = srcDS->GetRasterCount(); - CPLDebug("WX_MODEL_INITIALIZATION", "band count = %d", nBandCount); + //std::cout << "band count = " << nBandCount << std::endl; if( pbSuccess == false ) dfNoData = -9999.0; - /* - * set the dataset projection - */ + // set the dataset projection int rc = srcDS->SetProjection( projString.c_str() ); - /* - * final setting of the datasets to ascii grids, in the past usually done using a wrp dataset - */ + // final setting of the datasets to ascii grids, in the past usually done using a wrp dataset - CPLDebug("WX_MODEL_INITIALIZATION", "band number to write = %d", bandNum); + //std::cout << "band number to write = " << bandNum << std::endl; + // data should already in programming units, can just feed them in without a unit conversion if( varList[i] == "T2" ) { GDAL2AsciiGrid( srcDS, bandNum, airGrid ); - if( CPLIsNan( dfNoData ) ) { - airGrid.set_noDataValue(-9999.0); - airGrid.replaceNan( -9999.0 ); + if( CPLIsNan( dfNoData ) ) { + airGrid.set_noDataValue(-9999.0); + airGrid.replaceNan( -9999.0 ); + } } - } else if( varList[i] == "V10" ) { GDAL2AsciiGrid( srcDS, bandNum, vGrid ); - if( CPLIsNan( dfNoData ) ) { - vGrid.set_noDataValue(-9999.0); - vGrid.replaceNan( -9999.0 ); + if( CPLIsNan( dfNoData ) ) { + vGrid.set_noDataValue(-9999.0); + vGrid.replaceNan( -9999.0 ); + } } - } else if( varList[i] == "U10" ) { GDAL2AsciiGrid( srcDS, bandNum, uGrid ); - if( CPLIsNan( dfNoData ) ) { - uGrid.set_noDataValue(-9999.0); - uGrid.replaceNan( -9999.0 ); + if( CPLIsNan( dfNoData ) ) { + uGrid.set_noDataValue(-9999.0); + uGrid.replaceNan( -9999.0 ); + } } - } else if( varList[i] == "QCLOUD" ) { GDAL2AsciiGrid( srcDS, bandNum, cloudGrid ); - if( CPLIsNan( dfNoData ) ) { - cloudGrid.set_noDataValue(-9999.0); - cloudGrid.replaceNan( -9999.0 ); + if( CPLIsNan( dfNoData ) ) { + cloudGrid.set_noDataValue(-9999.0); + cloudGrid.replaceNan( -9999.0 ); + } } - } CPLFree(srcWKT); CPLFree(dstWkt); delete poCT; GDALClose((GDALDatasetH) srcDS ); - } - //don't allow small negative values in cloud cover + } // end for loop + + // don't allow small negative values in cloud cover for(int i=0; i &uGrid_wxModel, AsciiGrid &vGrid_wxModel ) +void writeWxModelGrids( const std::string &forecastFilename, const boost::local_time::local_date_time &forecastTime, const velocityUnits::eVelocityUnits &outputSpeedUnits, const AsciiGrid &uGrid_wxModel, const AsciiGrid &vGrid_wxModel ) { -/* - //// wxModelInitialization grid setting style - //// causes the input ascii grids to need to drop the "const" in front, claims set_headerData() discards qualifiers - AsciiGrid speedInitializationGrid_wxModel; - AsciiGrid dirInitializationGrid_wxModel; - - speedInitializationGrid_wxModel.set_headerData(uGrid_wxModel); - dirInitializationGrid_wxModel.set_headerData(uGrid_wxModel); - - //now make speed and direction from u,v components - for(int i=0; i speedInitializationGrid_wxModel( uGrid_wxModel ); AsciiGrid dirInitializationGrid_wxModel( uGrid_wxModel ); @@ -988,13 +835,22 @@ void writeWxModelGrids( const std::string &forecastFilename, const boost::local_ { for(int j=0; j timeList = getTimeList( timeZoneString, input_wrf_filename ); for(unsigned int timeIdx = 0; timeIdx < timeList.size(); timeIdx++) + //for(unsigned int timeIdx = 0; timeIdx < 1; timeIdx++) // if just want first time { boost::local_time::local_date_time forecastTime = timeList[timeIdx]; From 5e30cf1b00c5ef245bd54d0dde808ff0bfee08e8 Mon Sep 17 00:00:00 2001 From: latwood Date: Fri, 16 Aug 2024 18:16:10 -0600 Subject: [PATCH 3/5] added standalone hrrr_to_kmz script, for reading in hrrr grib2 files and outputting the corresponding kmz files, without having to run WindNinja itself each and every time this script differs from the wrf_to_kmz script, in that the dataset starts out as the entire conus region instead of already being downsized to a region of interest. So this script has an additional step to clip the data to a smaller region of interest there are still a few things to work out. Currently the script uses hard coded lat/lon points for the bounding box, the script still needs the inputs setup for various clipping box methods. In addition, the clipping is currently done using ascii grid clipping, after warping the dataset to lat/lon, but the clipping should probably be done using gdal itself in earlier steps --- CMakeLists.txt | 2 + src/CMakeLists.txt | 3 + src/hrrr_to_kmz/CMakeLists.txt | 48 ++ src/hrrr_to_kmz/hrrr_to_kmz.cpp | 630 ++++++++++++++++++++ src/ninja/ascii_grid.cpp | 65 ++ src/ninja/ascii_grid.h | 1 + src/slope_aspect_grid/slope_aspect_grid.cpp | 2 +- 7 files changed, 750 insertions(+), 1 deletion(-) create mode 100644 src/hrrr_to_kmz/CMakeLists.txt create mode 100644 src/hrrr_to_kmz/hrrr_to_kmz.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 75823c0a..360de3d4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -319,6 +319,8 @@ option(BUILD_STL_CONVERTER "Build a standalone command line interface for STL fi option(BUILD_CONVERT_OUTPUT "Build a standalone command line interface for xyz file conversions" OFF ) option(BUILD_WRF_TO_KMZ "Build a standalone command line interface for converting wrf output runs to kmz, without running full WindNinja" OFF ) mark_as_advanced(BUILD_WRF_TO_KMZ) +option(BUILD_HRRR_TO_KMZ "Build a standalone command line interface for converting hrrr output runs to kmz, without running full WindNinja" OFF ) +mark_as_advanced(BUILD_HRRR_TO_KMZ) option(BUILD_SLOPE_ASPECT_GRID "Build an application for building slope and aspect grids from a dem" OFF) mark_as_advanced(BUILD_SLOPE_ASPECT_GRID) option(BUILD_SOLAR_GRID "Build a application for building solar grids" OFF) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 74e6a609..1b51847b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -43,3 +43,6 @@ endif(BUILD_CONVERT_OUTPUT) if(BUILD_WRF_TO_KMZ) add_subdirectory(wrf_to_kmz) endif(BUILD_WRF_TO_KMZ) +if(BUILD_HRRR_TO_KMZ) + add_subdirectory(hrrr_to_kmz) +endif(BUILD_HRRR_TO_KMZ) diff --git a/src/hrrr_to_kmz/CMakeLists.txt b/src/hrrr_to_kmz/CMakeLists.txt new file mode 100644 index 00000000..495ebe95 --- /dev/null +++ b/src/hrrr_to_kmz/CMakeLists.txt @@ -0,0 +1,48 @@ +# THIS SOFTWARE WAS DEVELOPED AT THE ROCKY MOUNTAIN RESEARCH STATION (RMRS) +# MISSOULA FIRE SCIENCES LABORATORY BY EMPLOYEES OF THE FEDERAL GOVERNMENT +# IN THE COURSE OF THEIR OFFICIAL DUTIES. PURSUANT TO TITLE 17 SECTION 105 +# OF THE UNITED STATES CODE, THIS SOFTWARE IS NOT SUBJECT TO COPYRIGHT +# PROTECTION AND IS IN THE PUBLIC DOMAIN. RMRS MISSOULA FIRE SCIENCES +# LABORATORY ASSUMES NO RESPONSIBILITY WHATSOEVER FOR ITS USE BY OTHER +# PARTIES, AND MAKES NO GUARANTEES, EXPRESSED OR IMPLIED, ABOUT ITS QUALITY, +# RELIABILITY, OR ANY OTHER CHARACTERISTIC. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. + +cmake_minimum_required(VERSION 2.6) + +include_directories(${PROJECT_SOURCE_DIR}/src + ${PROJECT_SOURCE_DIR}/src/ninja + ${Boost_INCLUDE_DIRS} + ${GDAL_SYSTEM_INCLUDE} ${GDAL_INCLUDE_DIR}) + +set(LINK_LIBS ${GDAL_LIBRARY} + ${Boost_LIBRARIES}) + +set(HRRR_TO_KMZ_SRC hrrr_to_kmz.cpp + + ${PROJECT_SOURCE_DIR}/src/ninja/Array2D.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ascii_grid.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/EasyBMP.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/EasyBMP_Geometry.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/EasyBMP_Font.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ninja_conv.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ninja_init.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/gdal_util.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/KmlVector.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/Style.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/LineStyle.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ninjaMathUtility.cpp + ${PROJECT_SOURCE_DIR}/src/ninja/ninjaUnits.cpp) + +add_executable(hrrr_to_kmz ${HRRR_TO_KMZ_SRC}) +target_link_libraries(hrrr_to_kmz ${LINK_LIBS}) + +install(TARGETS hrrr_to_kmz DESTINATION bin COMPONENT apps) + diff --git a/src/hrrr_to_kmz/hrrr_to_kmz.cpp b/src/hrrr_to_kmz/hrrr_to_kmz.cpp new file mode 100644 index 00000000..8de0dc41 --- /dev/null +++ b/src/hrrr_to_kmz/hrrr_to_kmz.cpp @@ -0,0 +1,630 @@ +/****************************************************************************** + * + * $Id$ + * + * Project: WindNinja + * Purpose: Executable for converting hrrr grib2 files to kmz without running WindNinja itself in a simulation + * Author: Loren Atwood + * + ****************************************************************************** + * + * THIS SOFTWARE WAS DEVELOPED AT THE ROCKY MOUNTAIN RESEARCH STATION (RMRS) + * MISSOULA FIRE SCIENCES LABORATORY BY EMPLOYEES OF THE FEDERAL GOVERNMENT + * IN THE COURSE OF THEIR OFFICIAL DUTIES. PURSUANT TO TITLE 17 SECTION 105 + * OF THE UNITED STATES CODE, THIS SOFTWARE IS NOT SUBJECT TO COPYRIGHT + * PROTECTION AND IS IN THE PUBLIC DOMAIN. RMRS MISSOULA FIRE SCIENCES + * LABORATORY ASSUMES NO RESPONSIBILITY WHATSOEVER FOR ITS USE BY OTHER + * PARTIES, AND MAKES NO GUARANTEES, EXPRESSED OR IMPLIED, ABOUT ITS QUALITY, + * RELIABILITY, OR ANY OTHER CHARACTERISTIC. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + * + *****************************************************************************/ + +// got many of these functions from wxModelInitialization.cpp and ncepHrrrSurfInitialization.cpp + + +#include "ninja_init.h" +#include "ninja_conv.h" + +#include "ascii_grid.h" + +#include "ninjaUnits.h" +#include "KmlVector.h" + + + +/** +* Static identifier to determine if the file is a HRRR forecast. +* If don't have access to grib api, identificaion is done based +* on filename and id of 10u band. +* @param fileName grib2 filename +* +* PCM 01/04/23: account for different variable sets/versions of HRRR. We cannot rely on specific band numbers as the sources and/or filters +* for the HRRR datasets might differ +* TODO - this should check for all bands we need and keep respective band numbers so that we don't have to re-check later +* +* @return true if the forecast is a NCEP HRRR forecast +*/ +bool identify( std::string fileName ) +{ + GDALDataset *srcDS = (GDALDataset*)GDALOpenShared( fileName.c_str(), GA_ReadOnly ); + + if( srcDS == NULL ) { + throw badForecastFile("during identify(), Cannot open forecast file."); + return false; + + } else { + bool identified = false; + const int nRasterSets = srcDS->GetRasterCount(); + + for (int i=1; i<=nRasterSets && !identified; i++) { + GDALRasterBand *poBand = srcDS->GetRasterBand(i); + const char* comment = poBand->GetMetadataItem("GRIB_COMMENT"); + if (comment && strcmp(comment, "u-component of wind [m/s]") == 0){ + const char* description = poBand->GetDescription(); + if (strncmp(description, "10[m] ", 6) == 0){ + identified = true; + } + } + } + + GDALClose( (GDALDatasetH)srcDS ); + return identified; + } +} + + +/** +* Fetch the variable names +* @return a vector of variable names +*/ +std::vector getVariableList() +{ + std::vector varList; + varList.push_back( "2t" ); + varList.push_back( "10v" ); + varList.push_back( "10u" ); + varList.push_back( "tcc" ); // Total cloud cover + return varList; +} + + +/** +* Checks the downloaded data to see if it is all valid. +*/ +void checkForValidData( std::string wxModelFileName ) +{ + +} + + +/** +* Gets a time zone name string corresponding to a latitude and longitude point +* usually the input lat/lon point is the center latitude and longitude of the dataset, +* or the center latitude and longitude of a bounding box subset of the dataset. +* @param lat The latitude used to locate the timezone. +* @param lon The longitude used to locate the timezone. +* @return timeZoneString The time zone name for the lat/lon point, found from the file "date_time_zonespec.csv". +*/ +std::string getTimeZoneString( const double &lat, const double &lon ) +{ + std::string timeZoneString = FetchTimeZone(lon, lat, NULL); + if( timeZoneString == "" ) + { + fprintf(stderr, "Could not get timezone for lat,lon %f,%f location!!!\n", lat, lon); + std::exit(1); + } + + return timeZoneString; +} + + +// PCM - we can't assume fixed band numbers since HRRR formats change and might be filtered to reduce size +GDALRasterBand* get10UBand (GDALDataset *srcDS) { + const int nRasterSets = srcDS->GetRasterCount(); + + for (int i=1; i<=nRasterSets; i++) { + GDALRasterBand *poBand = srcDS->GetRasterBand(i); + const char* comment = poBand->GetMetadataItem("GRIB_COMMENT"); + if (comment && strcmp(comment, "u-component of wind [m/s]") == 0){ + const char* description = poBand->GetDescription(); + if (strncmp(description, "10[m] ", 6) == 0){ + return poBand; + } + } + } + + return NULL; +} + + +/** + * Fetch the list of times that the hrrr file holds. It is assumed that + * the time variable is "time" and the units string is "units". If this + * is not the case, this function needs to be rewritten. + * + * @param timeZoneString Time zone name from the file "date_time_zonespec.csv". + * @param wxModelFileName the input hrrr file to open and read times from. + * @throw runtime_error for bad file i/o + * @return a vector of boost::local_time::local_date_time objects for the forecast. + */ +std::vector +getTimeList( const std::string &timeZoneString, const std::string &wxModelFileName ) +{ + + const char *pszVariable = NULL; + + boost::local_time::time_zone_ptr timeZonePtr; + timeZonePtr = globalTimeZoneDB.time_zone_from_region(timeZoneString.c_str()); + if( timeZonePtr == NULL ) { + ostringstream os; + os << "The time zone string: " << timeZoneString.c_str() + << " does not match any in " + << "the time zone database file: date_time_zonespec.csv."; + throw std::runtime_error( os.str() ); + } + + std::vectortimeList; + + + // Just 1 time step per forecast file for now. + GDALDataset *srcDS; + + if( strstr( wxModelFileName.c_str(), ".tar" ) ){ + srcDS = (GDALDataset*)GDALOpenShared(CPLSPrintf( "/vsitar/%s", wxModelFileName.c_str() ), GA_ReadOnly); + } + else{ + srcDS = (GDALDataset*)GDALOpenShared( wxModelFileName.c_str(), GA_ReadOnly ); + } + + + if( srcDS == NULL ) { + throw badForecastFile("during getTimeList(), Cannot open forecast file."); + } + + // get time info from 10u band + GDALRasterBand *poBand = get10UBand(srcDS); + if (poBand == NULL) { + throw badForecastFile(" getTimeList() failed - no 10m UGRD band"); + } + + const char *vt; + vt = poBand->GetMetadataItem( "GRIB_VALID_TIME" ); + + std::string validTimeString( vt ); + + // Clean the string for the boost constructor. Remove the prefix, suffix, and delimiters + int pos = -1; + while( validTimeString.find( " " ) != validTimeString.npos ) { + pos = validTimeString.find( " " ); + validTimeString.erase( pos, 1 ); + } + pos = validTimeString.find( "secUTC" ); + if( pos != validTimeString.npos ) { + validTimeString.erase( pos, strlen( "secUTC" ) ); + } + + // Make a posix time in UTC/GMT, to call the boost::local_time::local_date_time input UTC time constructor + bpt::ptime time_t_epoch( boost::gregorian::date( 1970,1,1 ) ); + long validTime = atoi( validTimeString.c_str() ); + + bpt::time_duration duration( 0,0,validTime ); + + bpt::ptime first_pt( time_t_epoch + duration ); // UTC time + blt::local_date_time first_pt_local( first_pt, timeZonePtr ); // local time + timeList.push_back( first_pt_local ); + + + return timeList; +} + + +/** +* Sets the surface grids based on a hrrr forecast. +* @param wxModelFileName The hrrr input filename from which data are read from. +* @param timeBandIdx The band index from which data are read from, corresponding to a specific expected time from getTimeList(). +* @param airGrid The air temperature grid to be filled. +* @param cloudGrid The cloud cover grid to be filled. +* @param uGrid The u velocity grid to be filled. +* @param vGrid The v velocity grid to be filled. +* @param wGrid The w velocity grid to be filled (filled with zeros). +* @param west The west side of the clipping box +* @param east The east side of the clipping box +* @param south The south side of the clipping box +* @param north The north side of the clipping box +*/ +void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx, const std::vector &timeList, AsciiGrid &airGrid, AsciiGrid &cloudGrid, AsciiGrid &uGrid, AsciiGrid &vGrid, AsciiGrid &wGrid, const double &west, const double &east, const double &south, const double &north ) +{ + + // looks like the bands go in the same order as the timeList + int bandNum = timeBandIdx+1; // timeIdx is 0 to N-1, bandNum is 1 to N. + + + GDALDataset *srcDS; + srcDS = (GDALDataset*)GDALOpenShared( wxModelFileName.c_str(), GA_ReadOnly ); + + if( srcDS == NULL ) { + throw badForecastFile("Cannot open forecast file in setSurfaceGrids()"); + } + + GDALRasterBand *poBand; + const char *gc; + + double dfNoData; + int pbSuccess = false; + bool convertToKelvin = false; + + // Search time list for our time to identify our band number for cloud/speed/dir + // Right now, just one time step per file + std::vector bandList; + for(unsigned int i = 0; i < timeList.size(); i++) + { + if(i+1 == bandNum) + { + //std::cout << "bandNum = " << bandNum << std::endl; + for(unsigned int j = 1; j <= srcDS->GetRasterCount(); j++) + { + poBand = srcDS->GetRasterBand( j ); + gc = poBand->GetMetadataItem( "GRIB_COMMENT" ); + std::string bandName( gc ); + + if ( bandName.find("Temperature") == 0) { + CPLDebug("HRRR", "2-m T found..."); + gc = poBand->GetMetadataItem( "GRIB_SHORT_NAME" ); + std::string shortName( gc); + if (shortName == "2-HTGL") { + if (bandName == "Temperature [C]") { + convertToKelvin = true; + bandList.push_back( j); // 2t + + } else if (bandName != "Temperature [K]") { + bandList.push_back( j); // 2t + } else { + cout << "skipping unsupported forecast temperature unit: " << bandName << endl; + } + } + } + } + for(unsigned int j = 1; j <= srcDS->GetRasterCount(); j++) + { + poBand = srcDS->GetRasterBand( j ); + gc = poBand->GetMetadataItem( "GRIB_COMMENT" ); + std::string bandName( gc ); + + if( bandName.find( "v-component of wind [m/s]" ) != bandName.npos ){ + CPLDebug("HRRR", "v-component of wind found..."); + gc = poBand->GetMetadataItem( "GRIB_SHORT_NAME" ); + std::string bandName( gc ); + if( bandName.find( "10-HTGL" ) != bandName.npos ){ + bandList.push_back( j ); // 10v + break; + } + } + } + for(unsigned int j = 1; j <= srcDS->GetRasterCount(); j++) + { + poBand = srcDS->GetRasterBand( j ); + gc = poBand->GetMetadataItem( "GRIB_COMMENT" ); + std::string bandName( gc ); + + if( bandName.find( "u-component of wind [m/s]" ) != bandName.npos ){ + CPLDebug("HRRR", "u-component of wind found..."); + gc = poBand->GetMetadataItem( "GRIB_SHORT_NAME" ); + std::string bandName( gc ); + if( bandName.find( "10-HTGL" ) != bandName.npos ){ + bandList.push_back( j ); // 10u + dfNoData = poBand->GetNoDataValue( &pbSuccess ); + break; + } + } + } + for(unsigned int j = 1; j <= srcDS->GetRasterCount(); j++) + { + poBand = srcDS->GetRasterBand( j ); + gc = poBand->GetMetadataItem( "GRIB_COMMENT" ); + std::string bandName( gc ); + + if( bandName.find( "Total cloud cover [%]" ) != bandName.npos ){ + CPLDebug("HRRR", "cloud cover found..."); + gc = poBand->GetMetadataItem( "GRIB_SHORT_NAME" ); + std::string bandName( gc ); + if( bandName.find( "0-RESERVED" ) != bandName.npos || + bandName.find( "0-EATM" ) != bandName.npos){ + bandList.push_back( j ); // Total cloud cover in % + break; + } + } + } + } + } + + if(bandList.size() < 4) { + GDALClose((GDALDatasetH) srcDS ); + throw std::runtime_error("Not enough bands detected in HRRR forecast file."); + } + + // print out the bandList + // useful for knowing which specific bands to use when replicating the gdalwarp with command line utilities + //std::cout << "bandList ="; + //for ( unsigned int idx = 0; idx < bandList.size(); idx++ ) + //{ + // std::cout << " " << bandList[idx]; + //} + //std::cout << std::endl; + + + std::string srcWkt; + srcWkt = srcDS->GetProjectionRef(); + + OGRSpatialReference oSRS, *poLatLong; + oSRS.importFromWkt(srcWkt.c_str()); + poLatLong = oSRS.CloneGeogCS(); + char *dstWkt = NULL; + poLatLong->exportToWkt(&dstWkt); + delete poLatLong; + + + GDALDataset *wrpDS; + GDALWarpOptions* psWarpOptions; + + psWarpOptions = GDALCreateWarpOptions(); + + + int nBandCount = bandList.size(); + + psWarpOptions->nBandCount = nBandCount; + + psWarpOptions->panSrcBands = + (int *) CPLMalloc(sizeof(int) * nBandCount ); + psWarpOptions->panSrcBands[0] = bandList[0]; + psWarpOptions->panSrcBands[1] = bandList[1]; + psWarpOptions->panSrcBands[2] = bandList[2]; + psWarpOptions->panSrcBands[3] = bandList[3]; + + psWarpOptions->panDstBands = + (int *) CPLMalloc(sizeof(int) * nBandCount ); + psWarpOptions->panDstBands[0] = 1; + psWarpOptions->panDstBands[1] = 2; + psWarpOptions->panDstBands[2] = 3; + psWarpOptions->panDstBands[3] = 4; + + psWarpOptions->padfDstNoDataReal = + (double*) CPLMalloc( sizeof( double ) * nBandCount ); + psWarpOptions->padfDstNoDataImag = + (double*) CPLMalloc( sizeof( double ) * nBandCount ); + + + if( pbSuccess == false ) + dfNoData = -9999.0; + + + wrpDS = (GDALDataset*) GDALAutoCreateWarpedVRT( srcDS, srcWkt.c_str(), + dstWkt, + GRA_NearestNeighbour, + 1.0, psWarpOptions ); + + if (wrpDS == NULL) { + throw std::runtime_error("Warp operation failed!"); + } + + std::vector varList = getVariableList(); + + for( unsigned int i = 0; i < varList.size(); i++ ) { + if( varList[i] == "2t" ) { + GDAL2AsciiGrid( wrpDS, i+1, airGrid ); + if( CPLIsNan( dfNoData ) ) { + airGrid.set_noDataValue( -9999.0 ); + airGrid.replaceNan( -9999.0 ); + } + + if (convertToKelvin) { + airGrid += 273.15; + } + } + else if( varList[i] == "10v" ) { + GDAL2AsciiGrid( wrpDS, i+1, vGrid ); + if( CPLIsNan( dfNoData ) ) { + vGrid.set_noDataValue( -9999.0 ); + vGrid.replaceNan( -9999.0 ); + } + } + else if( varList[i] == "10u" ) { + GDAL2AsciiGrid( wrpDS, i+1, uGrid ); + if( CPLIsNan( dfNoData ) ) { + uGrid.set_noDataValue( -9999.0 ); + uGrid.replaceNan( -9999.0 ); + } + } + else if( varList[i] == "tcc" ) { + GDAL2AsciiGrid( wrpDS, i+1, cloudGrid ); + if( CPLIsNan( dfNoData ) ) { + cloudGrid.set_noDataValue( -9999.0 ); + cloudGrid.replaceNan( -9999.0 ); + } + } + } // end for loop + + // if there are any clouds set cloud fraction to 1, otherwise set to 0. + for(int i = 0; i < cloudGrid.get_nRows(); i++){ + for(int j = 0; j < cloudGrid.get_nCols(); j++){ + if(cloudGrid(i,j) < 0.0){ + cloudGrid(i,j) = 0.0; + } + else{ + cloudGrid(i,j) = 1.0; + } + } + } + + wGrid.set_headerData( uGrid ); + wGrid = 0.0; + + GDALDestroyWarpOptions( psWarpOptions ); + GDALClose((GDALDatasetH) srcDS ); + GDALClose((GDALDatasetH) wrpDS ); + + + // clip the ascii grids to the input bounding box + airGrid.clipGridInPlaceSnapToCells( west, east, south, north ); + cloudGrid.clipGridInPlaceSnapToCells( west, east, south, north ); + uGrid.clipGridInPlaceSnapToCells( west, east, south, north ); + vGrid.clipGridInPlaceSnapToCells( west, east, south, north ); + wGrid.clipGridInPlaceSnapToCells( west, east, south, north ); + +} + + +/** +* write the ascii grids to kmz for a wrf surface forecast +* @param forecastFilename The input forecast filename from which data were read from, to get the output path from. +* @param forecastTime The boost::local_time::local_date_time corresponding to the data to be plotted, corresponding to a specific time from getTimeList(). +* @param outputSpeedUnits The speed units to write the output kmz speed data to. +* @param uGrid The u velocity grid to be plotted. +* @param vGrid The v velocity grid to be plotted. +*/ +void writeWxModelGrids( const std::string &forecastFilename, const boost::local_time::local_date_time &forecastTime, const velocityUnits::eVelocityUnits &outputSpeedUnits, const AsciiGrid &uGrid_wxModel, const AsciiGrid &vGrid_wxModel ) +{ + + // make speed and direction grids from u,v components + + AsciiGrid speedInitializationGrid_wxModel( uGrid_wxModel ); + AsciiGrid dirInitializationGrid_wxModel( uGrid_wxModel ); + + for(int i=0; iformat("%m-%d-%Y_%H%M"); + wxModelTimestream << forecastTime; + std::string forecastIdentifier = "NCEP-HRRR-3km-SURFACE"; + std::string rootname = forecastIdentifier + "-" + wxModelTimestream.str(); + + + // don't forget the to output units unit conversion + velocityUnits::fromBaseUnits(speedInitializationGrid_wxModel, outputSpeedUnits); + + + // now do the kmz preparation and writing stuff + + KmlVector ninjaKmlFiles; + + ninjaKmlFiles.setKmlFile( CPLFormFilename(path.c_str(), rootname.c_str(), "kml") ); + ninjaKmlFiles.setKmzFile( CPLFormFilename(path.c_str(), rootname.c_str(), "kmz") ); + ////ninjaKmlFiles.setDemFile(dem_filename); // turns out to be redundant and doesn't do anything, which is good because don't want this dependency + + ninjaKmlFiles.setLegendFile( CPLFormFilename(path.c_str(), rootname.c_str(), "bmp") ); + ninjaKmlFiles.setSpeedGrid(speedInitializationGrid_wxModel, outputSpeedUnits); + ninjaKmlFiles.setDirGrid(dirInitializationGrid_wxModel); + + //ninjaKmlFiles.setLineWidth(1.0); // input.googLineWidth value + ninjaKmlFiles.setLineWidth(3.0); // input.wxModelGoogLineWidth value + std::string dateTimewxModelLegFileTemp = CPLFormFilename(path.c_str(), (rootname+"_date_time").c_str(), "bmp"); + ninjaKmlFiles.setTime(forecastTime); + ninjaKmlFiles.setDateTimeLegendFile(dateTimewxModelLegFileTemp, forecastTime); + ninjaKmlFiles.setWxModel(forecastIdentifier, forecastTime); + + // default values for input.wxModelGoogSpeedScaling/input.googSpeedScaling,input.googColor,input.googVectorScale are KmlVector::equal_interval,"default",false respectively + if(ninjaKmlFiles.writeKml(KmlVector::equal_interval,"default",false)) + { + if(ninjaKmlFiles.makeKmz()) + ninjaKmlFiles.removeKmlFile(); + } +} + + +int main( int argc, char* argv[] ) +{ + /* parse input arguments */ + if( argc != 3 ) + { + std::cout << "Invalid arguments!" << std::endl; + std::cout << "hrrr_to_kmz [input_hrrr_filename] [output_speed_units mph/mps/kph/kts]" << std::endl; + //std::cout << "hrrr_to_kmz [input_hrrr_filename] [output_speed_units mph/mps/kph/kts] [--bbox north south east west] [--point x y x_buf y_buf] [--buf_units mi/km/ft/m]" << std::endl; + return 1; + } + std::string input_hrrr_filename = std::string( argv[1] ); + std::string outputSpeedUnits_str = std::string( argv[2] ); + + std::cout << "input_hrrr_filename = \"" << input_hrrr_filename.c_str() << "\"" << std::endl; + std::cout << "output_speed_units = \"" << outputSpeedUnits_str.c_str() << "\"" << std::endl; + + double north = 47.18597932702905; + double south = 46.54752767224308; + double east = -113.45031738281251; + double west = -114.49401855468751; + + double cenLon = (west+east)/2; + double cenLat = (south+north)/2; + + + NinjaInitialize(); // needed for GDALAllRegister() + + + velocityUnits::eVelocityUnits outputSpeedUnits = velocityUnits::getUnit(outputSpeedUnits_str); + + + if ( identify( input_hrrr_filename ) == false ) + { + throw badForecastFile("input input_hrrr_filename is not a valid hrrr file!!!"); + } + checkForValidData( input_hrrr_filename ); + + + std::string timeZoneString = getTimeZoneString( cenLat, cenLon ); + + std::vector timeList = getTimeList( timeZoneString, input_hrrr_filename ); + + for(unsigned int timeIdx = 0; timeIdx < timeList.size(); timeIdx++) + //for(unsigned int timeIdx = 0; timeIdx < 1; timeIdx++) + { + boost::local_time::local_date_time forecastTime = timeList[timeIdx]; + + AsciiGrid airGrid; + AsciiGrid cloudGrid; + AsciiGrid uGrid; + AsciiGrid vGrid; + AsciiGrid wGrid; + + setSurfaceGrids( input_hrrr_filename, timeIdx, timeList, airGrid, cloudGrid, uGrid, vGrid, wGrid, west, east, south, north ); + + writeWxModelGrids( input_hrrr_filename, forecastTime, outputSpeedUnits, uGrid, vGrid ); + } + + + return 0; +} + + + + + + + + diff --git a/src/ninja/ascii_grid.cpp b/src/ninja/ascii_grid.cpp index f2043b9d..cdf5bef8 100644 --- a/src/ninja/ascii_grid.cpp +++ b/src/ninja/ascii_grid.cpp @@ -1113,6 +1113,71 @@ void AsciiGrid::clipGridInPlaceSnapToCells(double percentClip) } } +template +void AsciiGrid::clipGridInPlaceSnapToCells(double west, double east, double south, double north) +{ + if(west > east || south > north) + { + std::ostringstream buff_str; + buff_str << "The AsciiGrid clipping grid is improperly set to west,east,south,north " << west << "," << east << ", " << south << "," << north << ". bbox should be east > west and north > south."; + throw std::out_of_range(buff_str.str().c_str()); + } + + if( check_inBounds( west, south ) == false || check_inBounds( east, north ) == false ) + { + std::ostringstream buff_str; + buff_str << "The AsciiGrid clipping grid extends outside the bounds of the data!!"; + throw std::out_of_range(buff_str.str().c_str()); + } + + int westIdx = 0; + int eastIdx = 0; + int southIdx = 0; + int northIdx = 0; + // the get_cellIndex function, if in between cells, tends to round down rather than to round up + // this means slightly bigger than the bounding box for south and west, but slightly smaller than the bounding box for north and east + get_cellIndex(west, south, &southIdx, &westIdx); + get_cellIndex(east, north, &northIdx, &eastIdx); + + //std::cout << "westIdx,eastIdx,southIdx,northIdx = " << westIdx << "," << eastIdx << "," << southIdx << "," << northIdx << std::endl; + + if ( westIdx == eastIdx && southIdx == northIdx ) + { + // if there is zero clipping to be done, just return without doing anything + } else + { + int newNumCols, newNumRows; + int xClipCells, yClipCells; + double newXllCorner, newYllCorner; + + xClipCells = westIdx; + yClipCells = southIdx; + + newNumCols = eastIdx - westIdx; + newNumRows = northIdx - southIdx; + //std::cout << "xClipCells,yClipCells = " << xClipCells << "," << yClipCells << std::endl; + //std::cout << "newNumCols,newNumRows = " << newNumCols << "," << newNumRows << std::endl; + //std::cout << "xClipCells+newNumCols,yClipCells+newNumRows = " << xClipCells+newNumCols << "," << yClipCells+newNumRows << std::endl; + + newXllCorner = xllCorner + xClipCells*cellSize; + newYllCorner = yllCorner + yClipCells*cellSize; + //std::cout << "newXllCorner,newXllCorner+newNumCols*cellSize = " << newXllCorner << "," << newXllCorner+newNumCols*cellSize << std::endl; + //std::cout << "newYllCorner,newYllCorner+newNumRows*cellSize = " << newYllCorner << "," << newYllCorner+newNumRows*cellSize << std::endl; + + AsciiGridA(newNumCols, newNumRows, newXllCorner, newYllCorner, cellSize, data.getNoDataValue(), data.getNoDataValue(), prjString); + + for(int i=0; i AsciiGrid AsciiGrid::normalize_Grid(T lowBound, T highBound) { diff --git a/src/ninja/ascii_grid.h b/src/ninja/ascii_grid.h index ff64d534..28de11d0 100644 --- a/src/ninja/ascii_grid.h +++ b/src/ninja/ascii_grid.h @@ -173,6 +173,7 @@ class AsciiGrid double interpDistPower); void clipGridInPlaceSnapToCells(double percentClip); + void clipGridInPlaceSnapToCells(double north, double east, double south, double west); AsciiGrid normalize_Grid(T lowBound, T highBound); diff --git a/src/slope_aspect_grid/slope_aspect_grid.cpp b/src/slope_aspect_grid/slope_aspect_grid.cpp index 49b94e5f..926bf2cf 100644 --- a/src/slope_aspect_grid/slope_aspect_grid.cpp +++ b/src/slope_aspect_grid/slope_aspect_grid.cpp @@ -4,7 +4,7 @@ * * Project: WindNinja * Purpose: Application for creating a slope and aspect grid - * Author: Loren Atwoodr + * Author: Loren Atwood * ****************************************************************************** * From d8ead7850f3db33f58fa4679ec9de223182b597a Mon Sep 17 00:00:00 2001 From: nwagenbrenner Date: Mon, 19 Aug 2024 08:11:24 -0600 Subject: [PATCH 4/5] cleanup --- src/wrf_to_kmz/CMakeLists.txt | 1 - src/wrf_to_kmz/wrf_to_kmz.cpp | 80 +---------------------------------- 2 files changed, 2 insertions(+), 79 deletions(-) diff --git a/src/wrf_to_kmz/CMakeLists.txt b/src/wrf_to_kmz/CMakeLists.txt index da0687f4..c68d6942 100644 --- a/src/wrf_to_kmz/CMakeLists.txt +++ b/src/wrf_to_kmz/CMakeLists.txt @@ -28,7 +28,6 @@ set(LINK_LIBS ${NETCDF_LIBRARIES_C} ${Boost_LIBRARIES}) set(WRF_TO_KMZ_SRC wrf_to_kmz.cpp - ${PROJECT_SOURCE_DIR}/src/ninja/Array2D.cpp ${PROJECT_SOURCE_DIR}/src/ninja/ascii_grid.cpp ${PROJECT_SOURCE_DIR}/src/ninja/EasyBMP.cpp diff --git a/src/wrf_to_kmz/wrf_to_kmz.cpp b/src/wrf_to_kmz/wrf_to_kmz.cpp index e85f9ff6..d6aa426e 100644 --- a/src/wrf_to_kmz/wrf_to_kmz.cpp +++ b/src/wrf_to_kmz/wrf_to_kmz.cpp @@ -27,9 +27,6 @@ * *****************************************************************************/ -// got many of these functions from wxModelInitialization.cpp and wrfSurfInitialization.cpp - - #include "ninja_init.h" #include "ninja_conv.h" @@ -41,8 +38,6 @@ #include "ninjaUnits.h" #include "KmlVector.h" - - /** * Static identifier to determine if the netcdf file is a WRF forecast. * Uses netcdf c api. @@ -122,7 +117,6 @@ void checkForValidData( std::string wxModelFileName ) for( unsigned int i = 0;i < varList.size();i++ ) { temp = "NETCDF:\"" + wxModelFileName + "\":" + varList[i]; - //std::cout << "temp = \"" << temp.c_str() << "\"" << std::endl; CPLPushErrorHandler(&CPLQuietErrorHandler); srcDS = (GDALDataset*)GDALOpen( temp.c_str(), GA_ReadOnly ); @@ -132,7 +126,6 @@ void checkForValidData( std::string wxModelFileName ) // Get total bands (time steps) nBands = srcDS->GetRasterCount(); - //std::cout << "nBands = " << nBands << std::endl; nBands = srcDS->GetRasterCount(); int nXSize, nYSize; GDALRasterBand *poBand; @@ -143,9 +136,6 @@ void checkForValidData( std::string wxModelFileName ) nXSize = srcDS->GetRasterXSize(); nYSize = srcDS->GetRasterYSize(); - //std::cout << "nXsize = " << nXSize << std::endl; - //std::cout << "nYsize = " << nYSize << std::endl; - // loop over all bands for this variable (bands are time steps) for(int j = 1; j <= nBands; j++) { @@ -220,7 +210,7 @@ void checkForValidData( std::string wxModelFileName ) * @param cenLat The center latitude value of the dataset to be filled. * @param cenLon The center longitude value of the dataset to be filled. * @param projString The projection string of the dataset to be filled. -* note that there are additional wrf netcdf file global attributes that are accessed and used to get and process the projection string, +* Note: there are additional wrf netcdf file global attributes that are accessed and used to get and process the projection string, * but are not used outside this function, so they are not returned. These include mapProj, moadCenLat, standLon, trueLat1, trueLat2. */ void getNcGlobalAttributes( const std::string &wxModelFileName, float &dx, float &dy, float &cenLat, float &cenLon, std::string &projString ) @@ -257,10 +247,8 @@ void getNcGlobalAttributes( const std::string &wxModelFileName, float &dx, float else { status = nc_get_att_int( ncid, NC_GLOBAL, "MAP_PROJ", &mapProj ); } - //std::cout << "MAP_PROJ = " << mapProj << std::endl; // Get global attributes DX, DY - ////float dx, dy; // these are function inputs, to be filled status = nc_inq_att( ncid, NC_GLOBAL, "DX", &type, &len ); if( status != NC_NOERR ){ ostringstream os; @@ -283,7 +271,6 @@ void getNcGlobalAttributes( const std::string &wxModelFileName, float &dx, float } // Get global attributes CEN_LAT, CEN_LON - //float cenLat, cenLon; // these are function inputs, to be filled status = nc_inq_att( ncid, NC_GLOBAL, "CEN_LAT", &type, &len ); if( status != NC_NOERR ){ ostringstream os; @@ -406,10 +393,8 @@ void getNcGlobalAttributes( const std::string &wxModelFileName, float &dx, float Forecast file must be in a projected coordinate system."); } else throw badForecastFile("Cannot determine projection from the forecast file information."); - } - /** * Gets a time zone name string corresponding to a latitude and longitude point * usually the input lat/lon point is the center latitude and longitude of the dataset, @@ -430,7 +415,6 @@ std::string getTimeZoneString( const double &lat, const double &lon ) return timeZoneString; } - /** * Fetch the list of times that the wrf file holds. It is assumed that * the time variable is "Times" and the units string is "units". If this @@ -445,7 +429,6 @@ std::string getTimeZoneString( const double &lat, const double &lon ) std::vector getTimeList( const std::string &timeZoneString, const std::string &wxModelFileName ) { - const char *pszVariable = NULL; boost::local_time::time_zone_ptr timeZonePtr; @@ -460,7 +443,6 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa std::vectortimeList; - int status, ncid, ndims, nvars, ngatts, unlimdimid; nc_type vartype; int varndims, varnatts; @@ -476,11 +458,9 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - // note that "Times" is the time variable name found in wrf files, this variable is usually set to "time" for many of the other weather model file types std::string timename = "Times"; - // If we can't get simple data from the file, return false status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid); if ( status != NC_NOERR ) { @@ -501,7 +481,6 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - //==============The forecast is WRF, parse Times============================ // Check to see if we can read Times variable @@ -515,7 +494,6 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - // Get varid for U10 --> use this to get length of time dimension status = nc_inq_varid( ncid, "U10", &varid ); if( status != NC_NOERR ) { @@ -525,7 +503,6 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - // Get dimid for Time in U10 int dimid; status = nc_inq_dimid(ncid, "Time", &dimid ); @@ -537,7 +514,6 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa throw std::runtime_error( os.str() ); } - // Get length of the time dimension in U10 size_t time_len; status = nc_inq_dimlen( ncid, dimid, &time_len ); @@ -549,7 +525,6 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa } // Reset varid to 'Times' - //int varid; // already set into the scope above status = nc_inq_varid( ncid, timename.c_str(), &varid ); if( status != NC_NOERR ) { ostringstream os; @@ -559,9 +534,7 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa } // Get dimid for DateStrLen - //int dimid; // already set into the scope above status = nc_inq_dimid(ncid, "DateStrLen", &dimid ); - //std::cout << "dimid =" << dimid << std::endl; if( status != NC_NOERR ) { ostringstream os; os << "The dimension \"DateStrLen\" " @@ -584,7 +557,6 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa // Get value for one Times variable char* tp = new char[t_len + 1]; - //char tp[t_len + 1]; for (int i=0; i(i)}; // where to get value from status = nc_get_var1_text( ncid, varid, varindex, tp+i ); @@ -606,7 +578,6 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa if( refString.find( "_" ) != refString.npos ) { pos = refString.find( "_" ); refString.replace(pos, 1, 1, 'T'); - //refString.erase( pos, 1 ); } while( refString.find( "-" ) != refString.npos ) { pos = refString.find( "-" ); @@ -631,11 +602,9 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa } // end for loop - return timeList; } - /** * Sets the surface grids based on a WRF (surface only!) forecast. * @param wxModelFileName The wrf input filename from which data are read from. @@ -653,10 +622,8 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx, const float &dx, const float &dy, const float &cenLat, const float &cenLon, const std::string &projString, AsciiGrid &airGrid, AsciiGrid &cloudGrid, AsciiGrid &uGrid, AsciiGrid &vGrid, AsciiGrid &wGrid ) { - // looks like the bands go in the same order as the timeList, the past function looped through to find the idx that corresponded to an input time int bandNum = timeBandIdx+1; // timeIdx is 0 to N-1, bandNum is 1 to N. - GDALDataset* poDS; CPLPushErrorHandler(&CPLQuietErrorHandler); poDS = (GDALDataset*)GDALOpenShared( wxModelFileName.c_str(), GA_ReadOnly ); @@ -668,7 +635,6 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx GDALClose((GDALDatasetH) poDS ); // close original wxModel file } - // open ds one by one, set geotranform, set projection, then write to grid GDALDataset *srcDS; std::string temp; @@ -685,19 +651,14 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx throw badForecastFile("Cannot open forecast file in setSurfaceGrids()"); } - //std::cout << "varList[" << i << "] = " << varList[i].c_str() << std::endl; - // Set up spatial reference stuff for setting projections // and geotransformations - OGRSpatialReference oSRS, *poLatLong; char *srcWKT = NULL; char* prj2 = (char*)projString.c_str(); oSRS.importFromWkt(&prj2); oSRS.exportToWkt(&srcWKT); - //std::cout << "srcWKT = " << srcWKT << std::endl; - poLatLong = oSRS.CloneGeogCS(); char *dstWkt = NULL; poLatLong->exportToWkt(&dstWkt); @@ -722,12 +683,9 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx if(poCT==NULL || !poCT->Transform(1, &xCenter, &yCenter)) throw std::runtime_error("Transformation of lat/lon center to dataset coordinate system failed!"); - //std::cout << "xCenter = " << xCenter << ", yCenter = " << yCenter << std::endl; - // Set the geostransform for the WRF file // upper corner is calculated from transformed x, y // (in WRF space) - double ncols, nrows; int nXSize = srcDS->GetRasterXSize(); int nYSize = srcDS->GetRasterYSize(); @@ -738,34 +696,23 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx (yCenter+(nrows*dy)), 0, -(dy)}; - //std::cout << "ulcornerX, ulcornerY = " << xCenter-(ncols*dx) << ", " << yCenter+(nrows*dy) << std::endl; - //std::cout << "nXSize, nYsize = " << nXSize << ", " << nYSize << std::endl; - //std::cout << "dx = " << dx << std::endl; - srcDS->SetGeoTransform(adfGeoTransform); // get the noDataValue from the current band - GDALRasterBand *poBand = srcDS->GetRasterBand( bandNum ); int pbSuccess; double dfNoData = poBand->GetNoDataValue( &pbSuccess ); int nBandCount = srcDS->GetRasterCount(); - //std::cout << "band count = " << nBandCount << std::endl; - if( pbSuccess == false ) dfNoData = -9999.0; // set the dataset projection - int rc = srcDS->SetProjection( projString.c_str() ); // final setting of the datasets to ascii grids, in the past usually done using a wrp dataset - - //std::cout << "band number to write = " << bandNum << std::endl; - - // data should already in programming units, can just feed them in without a unit conversion + // TODO: data must be in SI units, need to check units here and convert if necessary if( varList[i] == "T2" ) { GDAL2AsciiGrid( srcDS, bandNum, airGrid ); if( CPLIsNan( dfNoData ) ) { @@ -814,7 +761,6 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx wGrid = 0.0; } - /** * write the ascii grids to kmz for a wrf surface forecast * @param forecastFilename The input forecast filename from which data were read from, to get the output path from. @@ -827,7 +773,6 @@ void writeWxModelGrids( const std::string &forecastFilename, const boost::local_ { // make speed and direction grids from u,v components - AsciiGrid speedInitializationGrid_wxModel( uGrid_wxModel ); AsciiGrid dirInitializationGrid_wxModel( uGrid_wxModel ); @@ -846,7 +791,6 @@ void writeWxModelGrids( const std::string &forecastFilename, const boost::local_ } } - // get path from input forecast filename std::string path = CPLGetPath(forecastFilename.c_str()); @@ -863,20 +807,16 @@ void writeWxModelGrids( const std::string &forecastFilename, const boost::local_ // don't forget the to output units unit conversion velocityUnits::fromBaseUnits(speedInitializationGrid_wxModel, outputSpeedUnits); - // now do the kmz preparation and writing stuff - KmlVector ninjaKmlFiles; ninjaKmlFiles.setKmlFile( CPLFormFilename(path.c_str(), rootname.c_str(), "kml") ); ninjaKmlFiles.setKmzFile( CPLFormFilename(path.c_str(), rootname.c_str(), "kmz") ); - ////ninjaKmlFiles.setDemFile(dem_filename); // turns out to be redundant and doesn't do anything, which is good because don't want this dependency ninjaKmlFiles.setLegendFile( CPLFormFilename(path.c_str(), rootname.c_str(), "bmp") ); ninjaKmlFiles.setSpeedGrid(speedInitializationGrid_wxModel, outputSpeedUnits); ninjaKmlFiles.setDirGrid(dirInitializationGrid_wxModel); - //ninjaKmlFiles.setLineWidth(1.0); // input.googLineWidth value ninjaKmlFiles.setLineWidth(3.0); // input.wxModelGoogLineWidth value std::string dateTimewxModelLegFileTemp = CPLFormFilename(path.c_str(), (rootname+"_date_time").c_str(), "bmp"); ninjaKmlFiles.setTime(forecastTime); @@ -891,7 +831,6 @@ void writeWxModelGrids( const std::string &forecastFilename, const boost::local_ } } - int main( int argc, char* argv[] ) { // parse input arguments @@ -907,14 +846,11 @@ int main( int argc, char* argv[] ) std::cout << "input_wrf_filename = \"" << input_wrf_filename.c_str() << "\"" << std::endl; std::cout << "output_speed_units = \"" << outputSpeedUnits_str.c_str() << "\"" << std::endl; - NinjaInitialize(); // needed for GDALAllRegister() - // test and set units velocityUnits::eVelocityUnits outputSpeedUnits = velocityUnits::getUnit(outputSpeedUnits_str); - // check dataset to verify it is a wrf dataset, and to verify it is readable if ( identify( input_wrf_filename ) == false ) { @@ -922,9 +858,7 @@ int main( int argc, char* argv[] ) } checkForValidData( input_wrf_filename ); - // now start processing the data - float dx; float dy; float cenLat; @@ -937,7 +871,6 @@ int main( int argc, char* argv[] ) std::vector timeList = getTimeList( timeZoneString, input_wrf_filename ); for(unsigned int timeIdx = 0; timeIdx < timeList.size(); timeIdx++) - //for(unsigned int timeIdx = 0; timeIdx < 1; timeIdx++) // if just want first time { boost::local_time::local_date_time forecastTime = timeList[timeIdx]; @@ -952,14 +885,5 @@ int main( int argc, char* argv[] ) writeWxModelGrids( input_wrf_filename, forecastTime, outputSpeedUnits, uGrid, vGrid ); } - return 0; } - - - - - - - - From 02633a46d8a62e9c3a55b395b2de58a23e537eec Mon Sep 17 00:00:00 2001 From: nwagenbrenner Date: Mon, 19 Aug 2024 08:15:58 -0600 Subject: [PATCH 5/5] cleanup --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 360de3d4..51dbd63a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -317,7 +317,7 @@ endif(NOT WIN32) option(BUILD_FETCH_DEM "Build a standalone command line interface DEM utility" OFF) option(BUILD_STL_CONVERTER "Build a standalone command line interface for STL file conversions" OFF ) option(BUILD_CONVERT_OUTPUT "Build a standalone command line interface for xyz file conversions" OFF ) -option(BUILD_WRF_TO_KMZ "Build a standalone command line interface for converting wrf output runs to kmz, without running full WindNinja" OFF ) +option(BUILD_WRF_TO_KMZ "Build a standalone command line interface for converting WRF output to kmz" OFF ) mark_as_advanced(BUILD_WRF_TO_KMZ) option(BUILD_HRRR_TO_KMZ "Build a standalone command line interface for converting hrrr output runs to kmz, without running full WindNinja" OFF ) mark_as_advanced(BUILD_HRRR_TO_KMZ)