Skip to content

Commit ab9aa13

Browse files
latwoodlatwood
authored
Merge pull request #754 from firelab/GUI-Redesign-fixWxModelNoData
Fix wxModel NO_DATA warp collar
2 parents 4572b68 + 1919997 commit ab9aa13

17 files changed

Lines changed: 361 additions & 242 deletions

‎src/hrrr_to_kmz/hrrr_to_kmz.cpp‎

Lines changed: 18 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -321,7 +321,7 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx
321321
std::string bandName( gc );
322322
if( bandName.find( "10-HTGL" ) != bandName.npos ){
323323
bandList.push_back( j ); // 10u
324-
dfNoData = poBand->GetNoDataValue( &pbSuccess );
324+
dfNoData = poBand->GetNoDataValue(&pbSuccess);
325325
break;
326326
}
327327
}
@@ -371,15 +371,19 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx
371371
poLatLong->exportToWkt(&dstWkt);
372372
delete poLatLong;
373373

374+
int nBandCount = bandList.size();
375+
376+
if(pbSuccess == false)
377+
{
378+
dfNoData = -9999.0;
379+
}
374380

375381
GDALDataset *wrpDS;
376382
GDALWarpOptions* psWarpOptions;
377383

384+
// Setup warp options
378385
psWarpOptions = GDALCreateWarpOptions();
379386

380-
381-
int nBandCount = bandList.size();
382-
383387
psWarpOptions->nBandCount = nBandCount;
384388

385389
psWarpOptions->panSrcBands =
@@ -400,10 +404,17 @@ void setSurfaceGrids( const std::string &wxModelFileName, const int &timeBandIdx
400404
(double*) CPLMalloc( sizeof( double ) * nBandCount );
401405
psWarpOptions->padfDstNoDataImag =
402406
(double*) CPLMalloc( sizeof( double ) * nBandCount );
407+
for(int b = 0; b < nBandCount; b++)
408+
{
409+
psWarpOptions->padfDstNoDataReal[b] = dfNoData;
410+
psWarpOptions->padfDstNoDataImag[b] = dfNoData;
411+
}
403412

404-
405-
if( pbSuccess == false )
406-
dfNoData = -9999.0;
413+
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA");
414+
if(pbSuccess == false) // if GDALGetRasterNoDataValue() fails to return that a NO_DATA value is in the source dataset
415+
{
416+
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", boost::lexical_cast<std::string>(dfNoData).c_str());
417+
}
407418

408419
// compute the coordinateTransformationAngle, the angle between the y coordinate grid lines of the pre-warped and warped datasets,
409420
// going FROM the y coordinate grid line of the pre-warped dataset TO the y coordinate grid line of the warped dataset

‎src/ninja/KmlVector.cpp‎

Lines changed: 87 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -1836,94 +1836,97 @@ bool KmlVector::writeVectors(VSILFILE *fileOut)
18361836
{
18371837
for(int j = 0;j < nC;j++)
18381838
{
1839-
yScale = 0.5;
1840-
s = spd(i,j);
1841-
// geTheta is the printed value (geographic coordinates), theta is the drawn value (projected coordinates) which is then reprojected (from projected to geographic coordinates)
1842-
// the formula for going from one projection to another is always prj2 = prj1 - coordinateTransformAngle_from_prj1_to_prj2
1843-
// but in this case, prj1 = dem, prj2 = kmz, and coordinateTransformAngle_from_dem_to_kmz = -coordinateTransformAngle_from_kmz_to_dem = -angleFromNorth
1844-
// this is because angleFromNorth is stored as a value going FROM N TO dem, but here we are going FROM dem TO N,
1845-
// so we need to use a negative value for angleFromNorth rather than a positive value
1846-
// so for this case, prj2 = prj1 - (-angleFromNorth) = prj1 + angleFromNorth, the two negative signs cancel
1847-
// But, if using coordinateTransformAngle_from_dem_to_kmz instead of the angleFromNorth value, make sure to go back to only a single "-" sign in the formula
1848-
geTheta = wrap0to360( dir(i,j) + angleFromNorth ); //convert FROM projected TO geographic coordinates
1849-
theta = dir(i,j) + 180.0;
1850-
1851-
if(s <= splitValue[1])
1852-
yScale *= 0.40;
1853-
else if(s <= splitValue[2])
1854-
yScale *= 0.60;
1855-
else if(s <= splitValue[3])
1856-
yScale *= 0.80;
1857-
else if(s <= splitValue[4])
1858-
yScale *= 1.0;
1859-
xScale = yScale * 0.40;
1860-
1861-
spd.get_cellPosition(i, j, &xCenter, &yCenter);
1862-
1863-
//xCenter = (cSize / 2.0) + (j * cSize + spd.get_xllCorner());
1864-
//yCenter = (cSize / 2.0) + ((nR - i - 1) * cSize) + spd.get_yllCorner();
1865-
1866-
if(theta > 360)
1839+
if(spd(i,j) != spd.get_noDataValue() && dir(i,j) != dir.get_noDataValue())
18671840
{
1868-
theta -= 360;
1869-
}
1870-
theta = 360 - theta;
1841+
yScale = 0.5;
1842+
s = spd(i,j);
1843+
1844+
// geTheta is the printed value (geographic coordinates), theta is the drawn value (projected coordinates) which is then reprojected (from projected to geographic coordinates)
1845+
// the formula for going from one projection to another is always prj2 = prj1 - coordinateTransformAngle_from_prj1_to_prj2
1846+
// but in this case, prj1 = dem, prj2 = kmz, and coordinateTransformAngle_from_dem_to_kmz = -coordinateTransformAngle_from_kmz_to_dem = -angleFromNorth
1847+
// this is because angleFromNorth is stored as a value going FROM N TO dem, but here we are going FROM dem TO N,
1848+
// so we need to use a negative value for angleFromNorth rather than a positive value
1849+
// so for this case, prj2 = prj1 - (-angleFromNorth) = prj1 + angleFromNorth, the two negative signs cancel
1850+
// But, if using coordinateTransformAngle_from_dem_to_kmz instead of the angleFromNorth value, make sure to go back to only a single "-" sign in the formula
1851+
geTheta = wrap0to360( dir(i,j) + angleFromNorth ); //convert FROM projected TO geographic coordinates
1852+
theta = dir(i,j) + 180.0;
18711853

1872-
theta = theta * (PI / 180);
1854+
if(s <= splitValue[1])
1855+
yScale *= 0.40;
1856+
else if(s <= splitValue[2])
1857+
yScale *= 0.60;
1858+
else if(s <= splitValue[3])
1859+
yScale *= 0.80;
1860+
else if(s <= splitValue[4])
1861+
yScale *= 1.0;
1862+
xScale = yScale * 0.40;
18731863

1874-
if( areEqual( s, 0.0 ) )
1875-
{
1876-
double square_size = 16;
1877-
xTip = xCenter - cSize / square_size;
1878-
yTip = yCenter + cSize / square_size;
1879-
xTail = xCenter + cSize / square_size;
1880-
yTail = yCenter + cSize / square_size;
1881-
xHeadLeft = xCenter - cSize / square_size;
1882-
yHeadLeft = yCenter - cSize / square_size;
1883-
xHeadRight = xCenter + cSize / square_size;
1884-
yHeadRight = yCenter - cSize / square_size;
1885-
}
1886-
else
1887-
{
1888-
xPoint = 0;
1889-
yPoint = (cSize * yScale);
1890-
1891-
//compute tip coordinates
1892-
xTip = (xPoint * cos(theta)) - (yPoint * sin(theta));
1893-
yTip = (xPoint * sin(theta)) + (yPoint * cos(theta));
1894-
//compute tail coordinates
1895-
xTail = -xTip;
1896-
yTail = -yTip;
1897-
1898-
//compute right and left coordinates for head
1899-
xPoint = (cSize * xScale);
1900-
yPoint = (cSize * yScale)-(cSize * xScale);
1901-
1902-
xHeadRight = (xPoint * cos(theta)) - (yPoint * sin(theta));
1903-
yHeadRight = (xPoint * sin(theta)) + (yPoint * cos(theta));
1904-
1905-
xPoint = -(cSize * xScale);
1906-
yPoint = (cSize * yScale) - (cSize * xScale);
1907-
xHeadLeft = (xPoint * cos(theta)) - (yPoint * sin(theta));
1908-
yHeadLeft = (xPoint * sin(theta)) + (yPoint * cos(theta));
1909-
1910-
//shift to global coordinates
1911-
xTip+=xCenter;
1912-
yTip+=yCenter;
1913-
xTail+=xCenter;
1914-
yTail+=yCenter;
1915-
xHeadRight += xCenter;
1916-
yHeadRight += yCenter;
1917-
xHeadLeft += xCenter;
1918-
yHeadLeft += yCenter;
1919-
}
1920-
coordTransform->Transform(1, &xTip, &yTip);
1921-
coordTransform->Transform(1, &xTail, &yTail);
1922-
coordTransform->Transform(1, &xHeadRight, &yHeadRight);
1923-
coordTransform->Transform(1, &xHeadLeft, &yHeadLeft);
1864+
spd.get_cellPosition(i, j, &xCenter, &yCenter);
1865+
1866+
//xCenter = (cSize / 2.0) + (j * cSize + spd.get_xllCorner());
1867+
//yCenter = (cSize / 2.0) + ((nR - i - 1) * cSize) + spd.get_yllCorner();
1868+
1869+
if(theta > 360)
1870+
{
1871+
theta -= 360;
1872+
}
1873+
theta = 360 - theta;
1874+
1875+
theta = theta * (PI / 180);
1876+
1877+
if( areEqual( s, 0.0 ) )
1878+
{
1879+
double square_size = 16;
1880+
xTip = xCenter - cSize / square_size;
1881+
yTip = yCenter + cSize / square_size;
1882+
xTail = xCenter + cSize / square_size;
1883+
yTail = yCenter + cSize / square_size;
1884+
xHeadLeft = xCenter - cSize / square_size;
1885+
yHeadLeft = yCenter - cSize / square_size;
1886+
xHeadRight = xCenter + cSize / square_size;
1887+
yHeadRight = yCenter - cSize / square_size;
1888+
}
1889+
else
1890+
{
1891+
xPoint = 0;
1892+
yPoint = (cSize * yScale);
1893+
1894+
//compute tip coordinates
1895+
xTip = (xPoint * cos(theta)) - (yPoint * sin(theta));
1896+
yTip = (xPoint * sin(theta)) + (yPoint * cos(theta));
1897+
//compute tail coordinates
1898+
xTail = -xTip;
1899+
yTail = -yTip;
1900+
1901+
//compute right and left coordinates for head
1902+
xPoint = (cSize * xScale);
1903+
yPoint = (cSize * yScale)-(cSize * xScale);
1904+
1905+
xHeadRight = (xPoint * cos(theta)) - (yPoint * sin(theta));
1906+
yHeadRight = (xPoint * sin(theta)) + (yPoint * cos(theta));
1907+
1908+
xPoint = -(cSize * xScale);
1909+
yPoint = (cSize * yScale) - (cSize * xScale);
1910+
xHeadLeft = (xPoint * cos(theta)) - (yPoint * sin(theta));
1911+
yHeadLeft = (xPoint * sin(theta)) + (yPoint * cos(theta));
1912+
1913+
//shift to global coordinates
1914+
xTip+=xCenter;
1915+
yTip+=yCenter;
1916+
xTail+=xCenter;
1917+
yTail+=yCenter;
1918+
xHeadRight += xCenter;
1919+
yHeadRight += yCenter;
1920+
xHeadLeft += xCenter;
1921+
yHeadLeft += yCenter;
1922+
}
1923+
coordTransform->Transform(1, &xTip, &yTip);
1924+
coordTransform->Transform(1, &xTail, &yTail);
1925+
coordTransform->Transform(1, &xHeadRight, &yHeadRight);
1926+
coordTransform->Transform(1, &xHeadLeft, &yHeadLeft);
1927+
1928+
// now start writing to the kml file
19241929

1925-
if(s != spd.get_noDataValue() && theta != dir.get_noDataValue())
1926-
{
19271930
VSIFPrintfL(fileOut, "<Placemark>");
19281931
//fprintf(fileOut, "\n<Icon><href>ffs_icon.ico</href></Icon>");
19291932
VSIFPrintfL(fileOut, "\n\t<name>Cell %d,%d</name>", i, j);

‎src/ninja/gcp_wx_init.cpp‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -668,7 +668,10 @@ void GCPWxModel::setSurfaceGrids(WindNinjaInputs& input,
668668
GDALRasterBandH hBand = GDALGetRasterBand(hSrcDS, 1);
669669
int bSuccess;
670670
dfNoData = GDALGetRasterNoDataValue(hBand, &bSuccess);
671-
if (!bSuccess) dfNoData = -9999.0;
671+
if(!bSuccess)
672+
{
673+
dfNoData = -9999.0;
674+
}
672675

673676
// Setup warp options
674677
psWarpOptions = GDALCreateWarpOptions();
@@ -679,6 +682,12 @@ void GCPWxModel::setSurfaceGrids(WindNinjaInputs& input,
679682
psWarpOptions->padfDstNoDataImag[i] = dfNoData;
680683
}
681684

685+
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA");
686+
if(bSuccess == false) // if GDALGetRasterNoDataValue() fails to return that a NO_DATA value is in the source dataset
687+
{
688+
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", boost::lexical_cast<std::string>(dfNoData).c_str());
689+
}
690+
682691
pszSrcWkt = GDALGetProjectionRef(hSrcDS);
683692

684693
// compute the coordinateTransformationAngle, the angle between the y coordinate grid lines of the pre-warped and warped datasets,

‎src/ninja/gdal_util.cpp‎

Lines changed: 18 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1507,9 +1507,13 @@ bool GDALWarpToUtm(const char* filename, GDALDatasetH& hSrcDS, GDALDatasetH& hDs
15071507
int nBandCount = GDALGetRasterCount( hDstDS );
15081508

15091509
double dfNoData = GDALGetRasterNoDataValue(hSrcBand, NULL);
1510+
// TODO: check to see if this is needed, when checking "INIT_DEST"="NO_DATA" methods
1511+
//int bSuccess;
1512+
//double dfNoData = GDALGetRasterNoDataValue(hSrcBand, &bSuccess);
15101513

15111514
GDALSetRasterNoDataValue(hDstBand, dfNoData);
15121515

1516+
// Setup warp options
15131517
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
15141518

15151519
psWarpOptions->hSrcDS = hSrcDS;
@@ -1522,13 +1526,21 @@ bool GDALWarpToUtm(const char* filename, GDALDatasetH& hSrcDS, GDALDatasetH& hDs
15221526
(double*) CPLMalloc( sizeof( double ) * nBandCount );
15231527
psWarpOptions->padfDstNoDataReal[0] = dfNoData;
15241528
psWarpOptions->padfDstNoDataImag[0] = dfNoData;
1529+
15251530
psWarpOptions->panSrcBands =
15261531
(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
15271532
psWarpOptions->panSrcBands[0] = 1;
15281533
psWarpOptions->panDstBands =
15291534
(int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
15301535
psWarpOptions->panDstBands[0] = 1;
15311536

1537+
// TODO: not sure how this affects surface_fetch products yet. Need to also check if this is needed for the other gdalWarpToUtm() function as well.
1538+
//psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA");
1539+
//if(bSuccess == false) // if GDALGetRasterNoDataValue() fails to return that a NO_DATA value is in the source dataset
1540+
//{
1541+
// psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", boost::lexical_cast<std::string>(dfNoData).c_str());
1542+
//}
1543+
15321544
psWarpOptions->pTransformerArg =
15331545
GDALCreateGenImgProjTransformer( hSrcDS,
15341546
GDALGetProjectionRef(hSrcDS),
@@ -1684,15 +1696,16 @@ bool GDALWarpToWKT_GDALAutoCreateWarpedVRT( GDALDatasetH& hSrcDS, int band, GDAL
16841696

16851697
GDALDatasetH hBand = GDALGetRasterBand( hSrcDS, band );
16861698
int bSuccess = false;
1687-
double dfNoData = GDALGetRasterNoDataValue( hBand, &bSuccess );
1688-
if( bSuccess == false )
1699+
double dfNoData = GDALGetRasterNoDataValue(hBand, &bSuccess);
1700+
if(bSuccess == false)
16891701
{
16901702
dfNoData = -9999.0;
16911703
}
16921704

16931705
const char *pszSrcWkt;
16941706
pszSrcWkt = GDALGetProjectionRef( hSrcDS );
16951707

1708+
// Setup warp options
16961709
GDALWarpOptions *psWarpOptions;
16971710
psWarpOptions = GDALCreateWarpOptions();
16981711

@@ -1716,10 +1729,10 @@ bool GDALWarpToWKT_GDALAutoCreateWarpedVRT( GDALDatasetH& hSrcDS, int band, GDAL
17161729
psWarpOptions->padfDstNoDataImag[i] = dfNoData;
17171730
}
17181731

1719-
psWarpOptions->papszWarpOptions = CSLSetNameValue( psWarpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA" );
1720-
if( bSuccess == false ) // if GDALGetRasterNoDataValue( hBand, &bSuccess ) fails to return that a NO_DATA value is in the source dataset
1732+
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA");
1733+
if(bSuccess == false) // if GDALGetRasterNoDataValue() fails to return that a NO_DATA value is in the source dataset
17211734
{
1722-
psWarpOptions->papszWarpOptions = CSLSetNameValue( psWarpOptions->papszWarpOptions, "INIT_DEST", boost::lexical_cast<std::string>(dfNoData).c_str() );
1735+
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", boost::lexical_cast<std::string>(dfNoData).c_str());
17231736
}
17241737

17251738
hDstDS = GDALAutoCreateWarpedVRT( hSrcDS, pszSrcWkt, pszDstWkt,

‎src/ninja/genericSurfInitialization.cpp‎

Lines changed: 15 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -414,18 +414,21 @@ void genericSurfInitialization::setSurfaceGrids( WindNinjaInputs &input,
414414
throw std::runtime_error("Could not get projection from forecast file, bad forecast file.");
415415
}
416416

417-
/*
418-
* Grab the first band to get the nodata value for the variable,
419-
* assume all bands have the same ndv
420-
*/
417+
int nBandCount = srcDS->GetRasterCount();
418+
419+
// Grab the first band to get the NO_DATA value for the variable,
420+
// assume all bands have the same NO_DATA value
421421
GDALRasterBand *poBand = srcDS->GetRasterBand( 1 );
422422
int pbSuccess;
423-
double dfNoData = poBand->GetNoDataValue( &pbSuccess );
423+
double dfNoData = poBand->GetNoDataValue(&pbSuccess);
424+
if(pbSuccess == false)
425+
{
426+
dfNoData = -9999.0;
427+
}
424428

429+
// Setup warp options
425430
psWarpOptions = GDALCreateWarpOptions();
426431

427-
int nBandCount = srcDS->GetRasterCount();
428-
429432
psWarpOptions->nBandCount = nBandCount;
430433

431434
psWarpOptions->padfDstNoDataReal =
@@ -438,12 +441,11 @@ void genericSurfInitialization::setSurfaceGrids( WindNinjaInputs &input,
438441
psWarpOptions->padfDstNoDataImag[b] = dfNoData;
439442
}
440443

441-
if( pbSuccess == false )
442-
dfNoData = -9999.0;
443-
444-
psWarpOptions->papszWarpOptions =
445-
CSLSetNameValue( psWarpOptions->papszWarpOptions,
446-
"INIT_DEST", "NO_DATA" );
444+
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", "NO_DATA");
445+
if(pbSuccess == false) // if GDALGetRasterNoDataValue() fails to return that a NO_DATA value is in the source dataset
446+
{
447+
psWarpOptions->papszWarpOptions = CSLSetNameValue(psWarpOptions->papszWarpOptions, "INIT_DEST", boost::lexical_cast<std::string>(dfNoData).c_str());
448+
}
447449

448450
// compute the coordinateTransformationAngle, the angle between the y coordinate grid lines of the pre-warped and warped datasets,
449451
// going FROM the y coordinate grid line of the pre-warped dataset TO the y coordinate grid line of the warped dataset

0 commit comments

Comments
 (0)