diff --git a/src/hrrr_to_kmz/hrrr_to_kmz.cpp b/src/hrrr_to_kmz/hrrr_to_kmz.cpp index 45d12a165..5e96c488f 100644 --- a/src/hrrr_to_kmz/hrrr_to_kmz.cpp +++ b/src/hrrr_to_kmz/hrrr_to_kmz.cpp @@ -361,6 +361,10 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx std::string srcWkt; srcWkt = srcDS->GetProjectionRef(); + if(srcWkt.empty()) + { + throw std::runtime_error("Could not get projection from forecast file, bad forecast file."); + } OGRSpatialReference oSRS, *poLatLong; oSRS.importFromWkt(srcWkt.c_str()); @@ -432,8 +436,9 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx GRA_NearestNeighbour, 1.0, psWarpOptions ); - if (wrpDS == NULL) { - throw std::runtime_error("Warp operation failed!"); + if(wrpDS == NULL) + { + throw std::runtime_error("Could not warp the forecast file, possibly non-uniform grid."); } std::vector varList = getVariableList(); @@ -614,15 +619,15 @@ void writeWxModelGrids( const std::string &outputPath, const boost::local_time:: } ninjaKmlFiles.setLegendFile( CPLFormFilename(outputPath.c_str(), rootname.c_str(), "bmp") ); - ninjaKmlFiles.setSpeedGrid(speedInitializationGrid_wxModel, outputSpeedUnits); + std::string dateTimewxModelLegFileTemp = CPLFormFilename(outputPath.c_str(), (rootname+"_date_time").c_str(), "bmp"); + ninjaKmlFiles.setDateTimeLegendFile(dateTimewxModelLegFileTemp, forecastTime); + ninjaKmlFiles.setSpeedGrid(speedInitializationGrid_wxModel, outputSpeedUnits); ninjaKmlFiles.setAngleFromNorth(angleFromNorth); ninjaKmlFiles.setDirGrid(dirInitializationGrid_wxModel); //ninjaKmlFiles.setLineWidth(1.0); // input.googLineWidth value ninjaKmlFiles.setLineWidth(3.0); // input.wxModelGoogLineWidth value - std::string dateTimewxModelLegFileTemp = CPLFormFilename(outputPath.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 diff --git a/src/ninja/wrf3dInitialization.cpp b/src/ninja/wrf3dInitialization.cpp index a8918afc4..1353992b3 100644 --- a/src/ninja/wrf3dInitialization.cpp +++ b/src/ninja/wrf3dInitialization.cpp @@ -421,6 +421,9 @@ void wrf3dInitialization::set3dGrids( WindNinjaInputs &input, Mesh const& mesh ) psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", boost::lexical_cast(dfNoData).c_str()); } + // set the dataset projection + int rc = srcDS->SetProjection(projString.c_str()); + wrpDS = (GDALDataset*) GDALAutoCreateWarpedVRT( srcDS, srcWKT, dstWkt.c_str(), GRA_NearestNeighbour, diff --git a/src/ninja/wrfSurfInitialization.cpp b/src/ninja/wrfSurfInitialization.cpp index d4e67bbcd..50189ba3a 100644 --- a/src/ninja/wrfSurfInitialization.cpp +++ b/src/ninja/wrfSurfInitialization.cpp @@ -122,6 +122,92 @@ int wrfSurfInitialization::getEndHour() return 84; } +/** +* function for converting the read in netcdf units to WindNinja units +* specifically for WindNinja velocityUnits +* throws an error if the input unit_string does not yet have a conversion written in the code +* if the input string is "", implies the units were not available for read in, the default units will be returned with a warning +* @param unit_string The read in netcdf unit string to be converted to a WindNinja velocityUnit +* @return spd_units The corresponding WindNinja velocityUnit enum to the input unit_string +*/ +velocityUnits::eVelocityUnits wrfSurfInitialization::processVelocityUnits(std::string unit_string) +{ + velocityUnits::eVelocityUnits spd_units; + spd_units = velocityUnits::metersPerSecond; + + if ( unit_string == "" ) + { + std::cout << "processVelocityUnits warning, input unit_string is empty implying the last read in netcdf band had no units to read in" + << "\n returning default velocityUnits \"" << velocityUnits::getString(spd_units).c_str() << "\"" << std::endl; + return spd_units; + } + + if ( unit_string == "m s-1" || unit_string == "mps" || unit_string == "m/s" ) + { + spd_units = velocityUnits::metersPerSecond; + } else if ( unit_string == "mi hr-1" || unit_string == "mph" || unit_string == "mi/hr" ) + { + spd_units = velocityUnits::milesPerHour; + } else if ( unit_string == "km hr-1" || unit_string == "kph" || unit_string == "km/hr" ) + { + spd_units = velocityUnits::kilometersPerHour; + } else if ( unit_string == "kts" || unit_string == "knots" ) + { + spd_units = velocityUnits::knots; + } else + { + ostringstream os; + os << "unknown unit_string \"" << unit_string << "\" input to processVelocityUnits\n" + << " need to update units processing code for input velocityUnit\n"; + throw std::runtime_error( os.str() ); + } + + return spd_units; +} + +/** +* function for converting the read in netcdf units to WindNinja units +* specifically for WindNinja temperatureUnits +* throws an error if the input unit_string does not yet have a conversion written in the code +* if the input string is "", implies the units were not available for read in, the default units will be returned with a warning +* @param unit_string The read in netcdf unit string to be converted to a WindNinja temperatureUnits +* @return T_units The corresponding WindNinja temperatureUnits enum to the input unit_string +*/ +temperatureUnits::eTempUnits wrfSurfInitialization::processTemperatureUnits(std::string unit_string) +{ + temperatureUnits::eTempUnits T_units; + T_units = temperatureUnits::K; + + if ( unit_string == "" ) + { + std::cout << "processTemperatureUnits warning, input unit_string is empty implying the last read in netcdf band had no units to read in" + << "\n returning default temperatureUnits \"" << temperatureUnits::getString(T_units).c_str() << "\"" << std::endl; + return T_units; + } + + if ( unit_string == "K" || unit_string == "degK" || unit_string == "deg K" ) + { + T_units = temperatureUnits::K; + } else if ( unit_string == "C" || unit_string == "degC" || unit_string == "deg C" ) + { + T_units = temperatureUnits::C; + } else if ( unit_string == "R" || unit_string == "degR" || unit_string == "degR" ) + { + T_units = temperatureUnits::R; + } else if ( unit_string == "F" || unit_string == "degF" || unit_string == "deg F" ) + { + T_units = temperatureUnits::F; + } else + { + ostringstream os; + os << "unknown unit_string \"" << unit_string << "\" input to processTemperatureUnits\n" + << " need to update units processing code for input temperatureUnits\n"; + throw std::runtime_error( os.str() ); + } + + return T_units; +} + /** * Checks the downloaded data to see if it is all valid. * May not be functional yet for this class... @@ -143,6 +229,9 @@ void wrfSurfInitialization::checkForValidData() bool noDataValueExists; bool noDataIsNan; + velocityUnits::eVelocityUnits spd_units = velocityUnits::metersPerSecond; // initialize to default units + temperatureUnits::eTempUnits T_units = temperatureUnits::K; + std::vector varList = getVariableList(); //Acquire a lock to protect the non-thread safe netCDF library @@ -158,8 +247,10 @@ void wrfSurfInitialization::checkForValidData() CPLPushErrorHandler(&CPLQuietErrorHandler); srcDS = (GDALDataset*)GDALOpen( temp.c_str(), GA_ReadOnly ); CPLPopErrorHandler(); - if( srcDS == NULL ) + if(srcDS == NULL) + { throw badForecastFile("Could not get NETCDF variable '"+varList[i]+"' from forecast file, bad forecast file."); + } //Get total bands (time steps) nBands = srcDS->GetRasterCount(); @@ -182,52 +273,67 @@ void wrfSurfInitialization::checkForValidData() poBand = srcDS->GetRasterBand( j ); pbSuccess = 0; - dfNoData = poBand->GetNoDataValue( &pbSuccess ); - if( pbSuccess == false ) + dfNoData = poBand->GetNoDataValue(&pbSuccess); + if(pbSuccess == false) + { noDataValueExists = false; + } else { noDataValueExists = true; noDataIsNan = std::isnan(dfNoData); } - //set the data + const char * poBand_units = poBand->GetUnitType(); + if(varList[i] == "T2") + { + T_units = processTemperatureUnits(poBand_units); + } + if(varList[i] == "V10" || varList[i] == "U10") + { + spd_units = processVelocityUnits(poBand_units); + } + + // set the data padfScanline = new double[nXSize*nYSize]; - poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, padfScanline, nXSize, nYSize, - GDT_Float64, 0, 0); + 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) + double current_val = padfScanline[k]; + // Check if value is no data (if no data value was defined in file) if(noDataValueExists) { if(noDataIsNan) { - if(std::isnan(padfScanline[k])) + if(std::isnan(current_val)) throw badForecastFile("Forecast file contains no_data values."); }else { - if(padfScanline[k] == dfNoData) + if(current_val == dfNoData) throw badForecastFile("Forecast file contains no_data values."); } } - if( varList[i] == "T2" ) //units are Kelvin + if( varList[i] == "T2" ) // usual units are Kelvin { - if(padfScanline[k] < 180.0 || padfScanline[k] > 340.0) //these are near the most extreme temperatures ever recored on earth + temperatureUnits::toBaseUnits( current_val, T_units ); // convert to Kelvin if not in Kelvin + if(current_val < 180.0 || current_val > 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" ) // usual units are m/s { - if(std::abs(padfScanline[k]) > 220.0) + velocityUnits::toBaseUnits( current_val, spd_units ); // convert to m/s if not in m/s + if(std::abs(current_val) > 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" ) // usual units are m/s { - if(std::abs(padfScanline[k]) > 220.0) + velocityUnits::toBaseUnits( current_val, spd_units ); // convert to m/s if not in m/s + if(std::abs(current_val) > 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" ) // usual units are kg/kg { - if(padfScanline[k] < -0.0001 || padfScanline[k] > 100.0) + if(current_val < -0.0001 || current_val > 100.0) throw badForecastFile("Total cloud cover is out of range in forecast file."); } } @@ -244,7 +350,6 @@ void wrfSurfInitialization::checkForValidData() * Uses netcdf c api * * @param fileName netcdf filename -* * @return true if the forecast is a WRF forecast */ bool wrfSurfInitialization::identify( std::string fileName ) @@ -293,69 +398,23 @@ bool wrfSurfInitialization::identify( std::string fileName ) return identified; } - /** -* Sets the surface grids based on a WRF (surface only!) forecast. -* @param input The WindNinjaInputs for misc. info. -* @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?). +* Uses netcfd c api commands to get wrf netcdf file global attributes, to use later +* 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 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 wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, - AsciiGrid &airGrid, - AsciiGrid &cloudGrid, - AsciiGrid &uGrid, - AsciiGrid &vGrid, - AsciiGrid &wGrid ) +void wrfSurfInitialization::getNcGlobalAttributes(float &dx, float &dy, float &cenLat, float &cenLon, std::string &projString) { - GDALWarpOptions* psWarpOptions = nullptr; - - int bandNum = -1; - - //get time list - std::vector timeList( getTimeList(input.ninjaTimeZone) ); - //Search time list for our time to identify our band number for cloud/speed/dir - for(unsigned int i = 0; i < timeList.size(); i++) - { - if(input.ninjaTime == timeList[i]) - { - bandNum = i + 1; - CPLDebug("WRF", "Grabbing band number: %d", bandNum); - break; - } - } - if(bandNum < 0) - throw std::runtime_error("Could not match ninjaTime with a band number in the forecast file."); + //==========get global attributes to set projection=========================== - GDALDataset* poDS; - //attempt to grab the projection from the dem? - //check for member prjString first - std::string dstWkt; - dstWkt = input.dem.prjString; - if ( dstWkt.empty() ) { - //try to open original - poDS = (GDALDataset*)GDALOpen( input.dem.fileName.c_str(), GA_ReadOnly ); - if( poDS == NULL ) { - throw std::runtime_error("Could not open input dem file."); - } - dstWkt = poDS->GetProjectionRef(); - if( dstWkt.empty() ) { - throw std::runtime_error("Could not get projection reference from input dem file."); - } - GDALClose((GDALDatasetH) poDS ); - } - -//==========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 ) { @@ -365,13 +424,11 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, 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; @@ -389,11 +446,7 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, CPLDebug("WRF", "MAP_PROJ = %d", mapProj); - /* - * Get global attributes DX, DY - * - */ - float dx, dy; + // Get global attributes DX, DY status = nc_inq_att( ncid, NC_GLOBAL, "DX", &type, &len ); if( status != NC_NOERR ){ ostringstream os; @@ -415,11 +468,7 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, status = nc_get_att_float( ncid, NC_GLOBAL, "DY", &dy ); } - /* - * Get global attributes CEN_LAT, CEN_LON - * - */ - float cenLat, cenLon; + // Get global attributes CEN_LAT, CEN_LON status = nc_inq_att( ncid, NC_GLOBAL, "CEN_LAT", &type, &len ); if( status != NC_NOERR ){ ostringstream os; @@ -441,10 +490,7 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, status = nc_get_att_float( ncid, NC_GLOBAL, "CEN_LON", &cenLon ); } - /* - * Get global attributes MOAD_CEN_LAT, STAND_LON - * - */ + // 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 ){ @@ -467,10 +513,7 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, status = nc_get_att_float( ncid, NC_GLOBAL, "STAND_LON", &standLon ); } - /* - * Get global attributes TRUELAT1, TRUELAT2 - * - */ + // Get global attributes TRUELAT1, TRUELAT2 float trueLat1, trueLat2; status = nc_inq_att( ncid, NC_GLOBAL, "TRUELAT1", &type, &len ); if( status != NC_NOERR ){ @@ -493,10 +536,7 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, status = nc_get_att_float( ncid, NC_GLOBAL, "TRUELAT2", &trueLat2 ); } - /* - * Get global attribute BOTTOM-TOP_GRID_DIMENSION - * - */ + // Get global attribute BOTTOM-TOP_GRID_DIMENSION status = nc_inq_att( ncid, NC_GLOBAL, "BOTTOM-TOP_GRID_DIMENSION", &type, &len ); if( status != NC_NOERR ){ ostringstream os; @@ -508,12 +548,10 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, else { status = nc_get_att_int( ncid, NC_GLOBAL, "BOTTOM-TOP_GRID_DIMENSION", &wxModel_nLayers ); } - - /* - * Get global attribute WEST-EAST_GRID_DIMENSION - * Not currently used. WX model x/y dims set based on - * reprojected image (in DEM space). - */ + + // Get global attribute WEST-EAST_GRID_DIMENSION + // Not currently used. WX model x/y dims set based on + // reprojected image (in DEM space). /*status = nc_inq_att( ncid, NC_GLOBAL, "WEST-EAST_GRID_DIMENSION", &type, &len ); if( status != NC_NOERR ){ ostringstream os; @@ -525,12 +563,10 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, else { status = nc_get_att_int( ncid, NC_GLOBAL, "WEST-EAST_GRID_DIMENSION", &wxModel_nCols ); }*/ - - /* - * Get global attribute SOUTH-NORTH_GRID_DIMENSION - * Not currently used. WX model x/y dims set based on - * reprojected image (in DEM space). - */ + + // Get global attribute SOUTH-NORTH_GRID_DIMENSION + // Not currently used. WX model x/y dims set based on + // reprojected image (in DEM space). /*status = nc_inq_att( ncid, NC_GLOBAL, "SOUTH-NORTH_GRID_DIMENSION", &type, &len ); if( status != NC_NOERR ){ ostringstream os; @@ -543,10 +579,7 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, status = nc_get_att_int( ncid, NC_GLOBAL, "SOUTH-NORTH_GRID_DIMENSION", &wxModel_nRows ); }*/ - /* - * Close the dataset - * - */ + // Close the dataset status = nc_close( ncid ); if( status != NC_NOERR ) { ostringstream os; @@ -555,6 +588,121 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, 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."); + } +} + +/** +* Sets the surface grids based on a WRF (surface only!) forecast. +* @param input The WindNinjaInputs for misc. info. +* @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 wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, + AsciiGrid &airGrid, + AsciiGrid &cloudGrid, + AsciiGrid &uGrid, + AsciiGrid &vGrid, + AsciiGrid &wGrid ) +{ + GDALWarpOptions* psWarpOptions = nullptr; + + int bandNum = -1; + + //get time list + std::vector timeList( getTimeList(input.ninjaTimeZone) ); + //Search time list for our time to identify our band number for cloud/speed/dir + for(unsigned int i = 0; i < timeList.size(); i++) + { + if(input.ninjaTime == timeList[i]) + { + bandNum = i + 1; + CPLDebug("WRF", "Grabbing band number: %d", bandNum); + break; + } + } + if(bandNum < 0) + throw std::runtime_error("Could not match ninjaTime with a band number in the forecast file."); + + GDALDataset* poDS; + //attempt to grab the projection from the dem? + //check for member prjString first + std::string dstWkt; + dstWkt = input.dem.prjString; + if ( dstWkt.empty() ) { + //try to open original + poDS = (GDALDataset*)GDALOpen( input.dem.fileName.c_str(), GA_ReadOnly ); + if( poDS == NULL ) { + throw std::runtime_error("Could not open input dem file."); + } + dstWkt = poDS->GetProjectionRef(); + if( dstWkt.empty() ) { + throw std::runtime_error("Could not get projection reference from input dem file."); + } + GDALClose((GDALDatasetH) poDS ); + } + +//==========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 + float dx; + float dy; + float cenLat; + float cenLon; + std::string projString; + getNcGlobalAttributes(dx, dy, cenLat, cenLon, projString); //======end get global attributes======================================== CPLPushErrorHandler(&CPLQuietErrorHandler); @@ -565,11 +713,14 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, } GDALClose((GDALDatasetH) poDS); // close original wxModel file - // open ds one by one, set projection, warp, then write to grid + // open ds one by one, set geotranform, set projection, warp, then write to grid GDALDataset *srcDS, *wrpDS; std::string temp; std::vector varList = getVariableList(); + velocityUnits::eVelocityUnits spd_units = velocityUnits::metersPerSecond; // initialize to default units + temperatureUnits::eTempUnits T_units = temperatureUnits::K; + double coordinateTransformationAngle = 0.0; /* @@ -592,54 +743,6 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, * Set up spatial reference stuff for setting projections * and geotransformations */ - - std::string projString; - 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."); - } - OGRSpatialReference oSRS, *poLatLong; char *srcWKT = NULL; char* prj2 = (char*)projString.c_str(); @@ -648,7 +751,6 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, CPLDebug("WRF", "srcWKT= %s", srcWKT); - OGRCoordinateTransformation *poCT; poLatLong = oSRS.CloneGeogCS(); /* @@ -666,6 +768,7 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, #endif /* GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,0,0) */ #endif /* GDAL_COMPUTE_VERSION */ + OGRCoordinateTransformation *poCT; poCT = OGRCreateCoordinateTransformation(poLatLong, &oSRS); delete poLatLong; @@ -708,9 +811,8 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, int nBandCount = srcDS->GetRasterCount(); CPLDebug("WRF", "band count = %d", nBandCount); - // Grab the first band to get the NO_DATA value for the variable, - // assume all bands have the same NO_DATA value - GDALRasterBand *poBand = srcDS->GetRasterBand( 1 ); + // get the noDataValue from the current band + GDALRasterBand *poBand = srcDS->GetRasterBand( bandNum ); int pbSuccess; double dfNoData = poBand->GetNoDataValue(&pbSuccess); if(pbSuccess == false) @@ -718,6 +820,16 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, dfNoData = -9999.0; } + const char * poBand_units = poBand->GetUnitType(); + if(varList[i] == "T2") + { + T_units = processTemperatureUnits(poBand_units); + } + if(varList[i] == "V10" || varList[i] == "U10") + { + spd_units = processVelocityUnits(poBand_units); + } + // Setup warp options psWarpOptions = GDALCreateWarpOptions(); psWarpOptions->padfDstNoDataReal = (double*)CPLMalloc(sizeof(double) * nBandCount); @@ -734,6 +846,9 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", boost::lexical_cast(dfNoData).c_str()); } + // set the dataset projection + int rc = srcDS->SetProjection(projString.c_str()); + // compute the coordinateTransformationAngle, the angle between the y coordinate grid lines of the pre-warped and warped datasets, // going FROM the y coordinate grid line of the pre-warped dataset TO the y coordinate grid line of the warped dataset // in this case, going FROM weather model projection coordinates TO dem projection coordinates @@ -783,32 +898,36 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, if( varList[i] == "T2" ) { GDAL2AsciiGrid( wrpDS, bandNum, airGrid ); - if( std::isnan( dfNoData ) ) { - airGrid.set_noDataValue(-9999.0); - airGrid.replaceNan( -9999.0 ); + temperatureUnits::toBaseUnits( airGrid, T_units ); + if( std::isnan( dfNoData ) ) { + airGrid.set_noDataValue(-9999.0); + airGrid.replaceNan( -9999.0 ); + } } - } else if( varList[i] == "V10" ) { GDAL2AsciiGrid( wrpDS, bandNum, vGrid ); - if( std::isnan( dfNoData ) ) { - vGrid.set_noDataValue(-9999.0); - vGrid.replaceNan( -9999.0 ); + velocityUnits::toBaseUnits( vGrid, spd_units ); + if( std::isnan( dfNoData ) ) { + vGrid.set_noDataValue(-9999.0); + vGrid.replaceNan( -9999.0 ); + } } - } else if( varList[i] == "U10" ) { GDAL2AsciiGrid( wrpDS, bandNum, uGrid ); - if( std::isnan( dfNoData ) ) { - uGrid.set_noDataValue(-9999.0); - uGrid.replaceNan( -9999.0 ); + velocityUnits::toBaseUnits( uGrid, spd_units ); + if( std::isnan( dfNoData ) ) { + uGrid.set_noDataValue(-9999.0); + uGrid.replaceNan( -9999.0 ); + } } - } else if( varList[i] == "QCLOUD" ) { GDAL2AsciiGrid( wrpDS, bandNum, cloudGrid ); - if( std::isnan( dfNoData ) ) { - cloudGrid.set_noDataValue(-9999.0); - cloudGrid.replaceNan( -9999.0 ); + if( std::isnan( dfNoData ) ) { + cloudGrid.set_noDataValue(-9999.0); + cloudGrid.replaceNan( -9999.0 ); + } } - } + CPLFree(srcWKT); delete poCT; GDALClose((GDALDatasetH) srcDS ); @@ -826,6 +945,7 @@ void wrfSurfInitialization::setSurfaceGrids( WindNinjaInputs &input, } } cloudGrid /= 100.0; + wGrid.set_headerData( uGrid ); wGrid = 0.0; diff --git a/src/ninja/wrfSurfInitialization.h b/src/ninja/wrfSurfInitialization.h index bfc70e686..81d344d5d 100644 --- a/src/ninja/wrfSurfInitialization.h +++ b/src/ninja/wrfSurfInitialization.h @@ -58,6 +58,11 @@ class wrfSurfInitialization : public wxModelInitialization virtual double Get_Wind_Height(); protected: + + velocityUnits::eVelocityUnits processVelocityUnits(std::string unit_string); + temperatureUnits::eTempUnits processTemperatureUnits(std::string unit_string); + void getNcGlobalAttributes(float &dx, float &dy, float &cenLat, float &cenLon, std::string &projString); + virtual void setSurfaceGrids( WindNinjaInputs &input, AsciiGrid &airGrid, AsciiGrid &cloudGrid, diff --git a/src/ninja/wxModelInitialization.cpp b/src/ninja/wxModelInitialization.cpp index ba324ec7d..c9993ca73 100644 --- a/src/ninja/wxModelInitialization.cpp +++ b/src/ninja/wxModelInitialization.cpp @@ -776,7 +776,7 @@ wxModelInitialization::getTimeList(const char *pszVariable, 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 \"" << timename << "\" in the netcdf file: " << wxModelFileName << " cannot be read\n"; throw std::runtime_error( os.str() ); } diff --git a/src/wrf_to_kmz/wrf_to_kmz.cpp b/src/wrf_to_kmz/wrf_to_kmz.cpp index 045a3b585..0c624ba25 100644 --- a/src/wrf_to_kmz/wrf_to_kmz.cpp +++ b/src/wrf_to_kmz/wrf_to_kmz.cpp @@ -128,6 +128,7 @@ temperatureUnits::eTempUnits processTemperatureUnits( std::string unit_string ) /** * 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. */ @@ -191,8 +192,15 @@ std::vector getVariableList() */ void checkForValidData( std::string wxModelFileName ) { - // open ds variable by variable GDALDataset *srcDS; + srcDS = (GDALDataset*)GDALOpen(wxModelFileName.c_str(), GA_ReadOnly); + if(srcDS == NULL) + { + throw std::runtime_error("Could not open forecast file, bad forecast file."); + } + GDALClose((GDALDatasetH)srcDS); + + // open ds variable by variable std::string temp; std::string srcWkt; int nBands = 0; @@ -211,12 +219,13 @@ void checkForValidData( std::string wxModelFileName ) CPLPushErrorHandler(&CPLQuietErrorHandler); srcDS = (GDALDataset*)GDALOpen( temp.c_str(), GA_ReadOnly ); CPLPopErrorHandler(); - if( srcDS == NULL ) - throw badForecastFile("Cannot open forecast file."); + if(srcDS == NULL) + { + throw badForecastFile("Could not get NETCDF variable '"+varList[i]+"' from forecast file, bad forecast file."); + } // Get total bands (time steps) nBands = srcDS->GetRasterCount(); - nBands = srcDS->GetRasterCount(); int nXSize, nYSize; GDALRasterBand *poBand; int pbSuccess; @@ -232,9 +241,11 @@ void checkForValidData( std::string wxModelFileName ) poBand = srcDS->GetRasterBand( j ); pbSuccess = 0; - dfNoData = poBand->GetNoDataValue( &pbSuccess ); - if( pbSuccess == false ) + dfNoData = poBand->GetNoDataValue(&pbSuccess); + if(pbSuccess == false) + { noDataValueExists = false; + } else { noDataValueExists = true; @@ -242,19 +253,18 @@ void checkForValidData( std::string wxModelFileName ) } const char * poBand_units = poBand->GetUnitType(); - if( varList[i] == "T2" ) + if(varList[i] == "T2") { - T_units = processTemperatureUnits( poBand_units ); + T_units = processTemperatureUnits(poBand_units); } - if( varList[i] == "V10" || varList[i] == "U10" ) + if(varList[i] == "V10" || varList[i] == "U10") { - spd_units = processVelocityUnits( poBand_units ); + spd_units = processVelocityUnits(poBand_units); } // set the data padfScanline = new double[nXSize*nYSize]; - poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, padfScanline, nXSize, nYSize, - GDT_Float64, 0, 0); + poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, padfScanline, nXSize, nYSize, GDT_Float64, 0, 0); for(int k = 0;k < nXSize*nYSize; k++) { double current_val = padfScanline[k]; @@ -305,7 +315,7 @@ void checkForValidData( std::string wxModelFileName ) /** -* Uses netcfd c api commands to get wrf netcdf file global attributes, to later use +* Uses netcfd c api commands to get wrf netcdf file global attributes, to use later * 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. @@ -314,7 +324,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: 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 ) @@ -580,8 +590,7 @@ getTimeList( const std::string &timeZoneString, const std::string &wxModelFileNa 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: " + os << "The variable \"" << timename << "\" in the netcdf file: " << wxModelFileName << " cannot be read\n"; throw std::runtime_error( os.str() ); } @@ -743,11 +752,9 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx 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 + throw std::runtime_error("Could not open forecast file, bad forecast file."); } + GDALClose((GDALDatasetH) poDS); // close original wxModel file // open ds one by one, set geotranform, set projection, then write to grid GDALDataset *srcDS; @@ -768,7 +775,7 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx srcDS = (GDALDataset*)GDALOpenShared( temp.c_str(), GA_ReadOnly ); CPLPopErrorHandler(); if( srcDS == NULL ) { - throw badForecastFile("Cannot open forecast file in setSurfaceGrids()"); + throw std::runtime_error("Could not get NETCDF variable '"+varList[i]+"' from forecast file, bad forecast file."); } // Set up spatial reference stuff for setting projections @@ -801,7 +808,9 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx delete poLatLong; if(poCT==NULL || !poCT->Transform(1, &xCenter, &yCenter)) - throw std::runtime_error("Transformation of lat/lon center to dataset coordinate system failed!"); + { + throw std::runtime_error("Transformation of dataset center from dataset coordinate system to lat/lon failed."); + } // Set the geostransform for the WRF file // upper corner is calculated from transformed x, y @@ -818,30 +827,30 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx srcDS->SetGeoTransform(adfGeoTransform); + int nBandCount = srcDS->GetRasterCount(); + // 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 ) + double dfNoData = poBand->GetNoDataValue(&pbSuccess); + if(pbSuccess == false) + { dfNoData = -9999.0; + } const char * poBand_units = poBand->GetUnitType(); - if( varList[i] == "T2" ) + if(varList[i] == "T2") { - T_units = processTemperatureUnits( poBand_units ); + T_units = processTemperatureUnits(poBand_units); } - if( varList[i] == "V10" || varList[i] == "U10" ) + if(varList[i] == "V10" || varList[i] == "U10") { - spd_units = processVelocityUnits( poBand_units ); + spd_units = processVelocityUnits(poBand_units); } // set the dataset projection int rc = srcDS->SetProjection( projString.c_str() ); - //// compute the coordinateTransformationAngle, the angle between the y coordinate grid lines of the pre-warped and warped datasets, //// going FROM the y coordinate grid line of the pre-warped dataset TO the y coordinate grid line of the warped dataset //// in this case, going FROM weather model projection coordinates TO geographic lat/lon coordinates @@ -860,7 +869,6 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx // 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 ); temperatureUnits::toBaseUnits( airGrid, T_units ); @@ -1043,14 +1051,14 @@ void writeWxModelGrids( const std::string &outputPath, const boost::local_time:: } ninjaKmlFiles.setLegendFile( CPLFormFilename(outputPath.c_str(), rootname.c_str(), "bmp") ); + std::string dateTimewxModelLegFileTemp = CPLFormFilename(outputPath.c_str(), (rootname+"_date_time").c_str(), "bmp"); + ninjaKmlFiles.setDateTimeLegendFile(dateTimewxModelLegFileTemp, forecastTime); ninjaKmlFiles.setSpeedGrid(speedInitializationGrid_wxModel, outputSpeedUnits); ninjaKmlFiles.setAngleFromNorth(angleFromNorth); ninjaKmlFiles.setDirGrid(dirInitializationGrid_wxModel); ninjaKmlFiles.setLineWidth(3.0); // input.wxModelGoogLineWidth value - std::string dateTimewxModelLegFileTemp = CPLFormFilename(outputPath.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 @@ -1260,12 +1268,12 @@ int main( int argc, char* argv[] ) 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 ) + if(identify(input_wrf_filename) == false) { NinjaFinalize(); throw badForecastFile("input input_wrf_filename is not a valid WRF file!!!"); } - checkForValidData( input_wrf_filename ); + checkForValidData(input_wrf_filename); // now start processing the data float dx;