diff --git a/CMakeLists.txt b/CMakeLists.txt index f59f240d..f8ce0b63 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -327,6 +327,10 @@ 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 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) 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..1b51847b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -40,3 +40,9 @@ 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) +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 * ****************************************************************************** * diff --git a/src/wrf_to_kmz/CMakeLists.txt b/src/wrf_to_kmz/CMakeLists.txt new file mode 100644 index 00000000..c68d6942 --- /dev/null +++ b/src/wrf_to_kmz/CMakeLists.txt @@ -0,0 +1,49 @@ +# 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..d6aa426e --- /dev/null +++ b/src/wrf_to_kmz/wrf_to_kmz.cpp @@ -0,0 +1,889 @@ +/****************************************************************************** + * + * $Id$ + * + * Project: WindNinja + * Purpose: Executable for converting wrf netcdf 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. + * + *****************************************************************************/ + +#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; + + // 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. +* @param wxModelFileName wrf netcdf weather model filename. +* Comment From WindNinja code: 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(); + + for( unsigned int i = 0;i < varList.size();i++ ) { + + temp = "NETCDF:\"" + wxModelFileName + "\":" + varList[i]; + + 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(); + nBands = srcDS->GetRasterCount(); + int nXSize, nYSize; + GDALRasterBand *poBand; + int pbSuccess; + double dfNoData; + double *padfScanline; + + nXSize = srcDS->GetRasterXSize(); + nYSize = srcDS->GetRasterYSize(); + + // 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 ); + } +} + + +/** +* 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: 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=========================== + + // 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 ); + } + + // Get global attributes DX, DY + status = nc_inq_att( ncid, NC_GLOBAL, "DX", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute DX in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_float( ncid, NC_GLOBAL, "DX", &dx ); + } + status = nc_inq_att( ncid, NC_GLOBAL, "DY", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute DY in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_float( ncid, NC_GLOBAL, "DY", &dy ); + } + + // Get global attributes CEN_LAT, CEN_LON + status = nc_inq_att( ncid, NC_GLOBAL, "CEN_LAT", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute CEN_LAT in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_float( ncid, NC_GLOBAL, "CEN_LAT", &cenLat ); + } + status = nc_inq_att( ncid, NC_GLOBAL, "CEN_LON", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute CEN_LON in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_float( ncid, NC_GLOBAL, "CEN_LON", &cenLon ); + } + + // Get global attributes MOAD_CEN_LAT, STAND_LON + float moadCenLat, standLon; + status = nc_inq_att( ncid, NC_GLOBAL, "MOAD_CEN_LAT", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute MOAD_CEN_LAT in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_float( ncid, NC_GLOBAL, "MOAD_CEN_LAT", &moadCenLat ); + } + status = nc_inq_att( ncid, NC_GLOBAL, "STAND_LON", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute STAND_LON in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_float( ncid, NC_GLOBAL, "STAND_LON", &standLon ); + } + + // Get global attributes TRUELAT1, TRUELAT2 + float trueLat1, trueLat2; + status = nc_inq_att( ncid, NC_GLOBAL, "TRUELAT1", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute TRUELAT1 in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_float( ncid, NC_GLOBAL, "TRUELAT1", &trueLat1 ); + } + status = nc_inq_att( ncid, NC_GLOBAL, "TRUELAT2", &type, &len ); + if( status != NC_NOERR ){ + ostringstream os; + os << "Global attribute TRUELAT2 in the netcdf file: " << wxModelFileName + << " cannot be opened\n"; + throw std::runtime_error( os.str() ); + } + else { + status = nc_get_att_float( ncid, NC_GLOBAL, "TRUELAT2", &trueLat2 ); + } + + // Close the dataset + status = nc_close( ncid ); + if( status != NC_NOERR ) { + ostringstream os; + os << "The netcdf file: " << wxModelFileName + << " cannot be closed\n"; + throw std::runtime_error( os.str() ); + } + + //======end get global attributes======================================== + + + //std::string projString; // this is a function input, to be filled + if(mapProj == 1){ //lambert conformal conic + projString = "PROJCS[\"WGC 84 / WRF Lambert\",GEOGCS[\"WGS 84\",DATUM[\"World Geodetic System 1984\",\ + SPHEROID[\"WGS 84\",6378137.0,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],\ + PRIMEM[\"Greenwich\",0.0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.017453292519943295],\ + AXIS[\"Geodetic longitude\",EAST],AXIS[\"Geodetic latitude\",NORTH],AUTHORITY[\"EPSG\",\"4326\"]],\ + PROJECTION[\"Lambert_Conformal_Conic_2SP\"],\ + PARAMETER[\"central_meridian\","+boost::lexical_cast(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."); +} + +/** +* 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; +} + +/** + * 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 rewritten. + * Uses netcdf api commands. + * + * @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::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; + + 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 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 ) { + 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() ); + } + + //==============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' + 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 + status = nc_inq_dimid(ncid, "DateStrLen", &dimid ); + if( status != NC_NOERR ) { + ostringstream os; + os << "The dimension \"DateStrLen\" " + << " in the netcdf file: " << wxModelFileName + << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + // Get the length of the time string + size_t t_len; + status = nc_inq_dimlen( ncid, dimid, &t_len ); + if( status != NC_NOERR ) { + ostringstream os; + os << "The dimension for variable \"Times\" in the netcdf file: " + << wxModelFileName << " cannot be read\n"; + throw std::runtime_error( os.str() ); + } + + for( unsigned int t = 0;t < time_len;t++ ) { + + // Get value for one Times variable + char* tp = new char[t_len + 1]; + for (int i=0; i(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'); + } + 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, 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 ); + 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 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). +*/ +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 ) +{ + + 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 ) { + throw badForecastFile("Cannot open forecast file in 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 ) { + throw badForecastFile("Cannot open forecast file in setSurfaceGrids()"); + } + + // 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); + + 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)) + throw std::runtime_error("Transformation of lat/lon center to dataset coordinate system failed!"); + + // 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)}; + + 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(); + + 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 + // 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 ) ) { + 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 ); + } // end for loop + + // don't allow small negative values in cloud cover + for(int i=0; i &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 = "WRF-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.setLegendFile( CPLFormFilename(path.c_str(), rootname.c_str(), "bmp") ); + ninjaKmlFiles.setSpeedGrid(speedInitializationGrid_wxModel, outputSpeedUnits); + ninjaKmlFiles.setDirGrid(dirInitializationGrid_wxModel); + + 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 mph/mps/kph/kts]" << 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() + + // 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 ) + { + throw badForecastFile("input input_wrf_filename is not a valid WRF file!!!"); + } + checkForValidData( input_wrf_filename ); + + // now start processing the data + float dx; + float dy; + float cenLat; + float cenLon; + std::string projString; + getNcGlobalAttributes( input_wrf_filename, dx, dy, 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; +}