Commit 7bdbc0cb authored by Marc R.'s avatar Marc R.

structured grid: added methods to extract a vertical profile at a given geographical location

parent 7609140e
......@@ -1015,6 +1015,31 @@ void MSkewTActor::printDebugOutputOnUserRequest()
.arg(tp.x()).arg(tp.y()).arg(xp.x()).arg(xp.y());
}
}
// Print (T, p) profile at current (lon, lat) position.
float lon = diagramConfiguration.geoPosition.x();
float lat = diagramConfiguration.geoPosition.y();
str += "\n\n\nDEBUG output, vertical (T, p) profile at current location: ";
str += QString("(%1, %2).\n").arg(lon).arg(lat);
for (int vi = 0; vi < variables.size(); vi++)
{
MNWPSkewTActorVariable* var =
static_cast<MNWPSkewTActorVariable*> (variables.at(vi));
if ( !var->hasData() ) continue;
QVector<QVector2D> profile = var->grid->extractVerticalProfile(lon, lat);
str += QString("\nGRID: %1\n").arg(var->grid->getGeneratingRequest());
for (QVector2D tpCoordinate : profile)
{
str += QString("T=%1, p=%2\n")
.arg(tpCoordinate.x()).arg(tpCoordinate.y());
}
}
str += "DEBUG output end.\n";
LOG4CPLUS_DEBUG(mlog, str.toStdString());
}
......
......@@ -4,10 +4,13 @@
** three-dimensional visual exploration of numerical ensemble weather
** prediction data.
**
** Copyright 2015-2018 Marc Rautenhaus
** Copyright 2017-2018 Bianca Tost
** Copyright 2015-2019 Marc Rautenhaus [*, previously +]
** Copyright 2017-2018 Bianca Tost [+]
**
** Computer Graphics and Visualization Group
** * Regional Computing Center, Visualization
** Universitaet Hamburg, Hamburg, Germany
**
** + Computer Graphics and Visualization Group
** Technische Universitaet Muenchen, Garching, Germany
**
** Met.3D is free software: you can redistribute it and/or modify
......@@ -315,17 +318,27 @@ void MStructuredGrid::setToValue(float val)
}
float MStructuredGrid::interpolateValue(float lon, float lat, float p_hPa)
void MStructuredGrid::findEnclosingHorizontalIndices(
float lon, float lat, int *i, int *j, int *i1, int *j1,
float *mixI, float *mixJ)
{
// Find horizontal indices i,i+1 and j,j+1 that enclose lon, lat.
float mixI = MMOD(lon - lons[0], 360.) / abs(lons[1]-lons[0]);
float mixJ = (lats[0] - lat) / abs(lats[1]-lats[0]);
int i = int(mixI);
int j = int(mixJ);
*mixI = MMOD(lon - lons[0], 360.) / abs(lons[1]-lons[0]);
*mixJ = (lats[0] - lat) / abs(lats[1]-lats[0]);
*i = int(*mixI);
*j = int(*mixJ);
*i1 = (*i)+1;
if (gridIsCyclicInLongitude()) *i1 %= nlons;
*j1 = (*j)+1;
}
int i1 = i+1;
if (gridIsCyclicInLongitude()) i1 %= nlons;
int j1 = j+1;
float MStructuredGrid::interpolateValue(float lon, float lat, float p_hPa)
{
int i, j, i1, j1;
float mixI, mixJ;
findEnclosingHorizontalIndices(lon, lat, &i, &j, &i1, &j1, &mixI, &mixJ);
if ((i < 0) || (j < 0) || (i1 >= int(nlons)) || (j1 >= int(nlats)))
return M_MISSING_VALUE;
......@@ -355,19 +368,56 @@ float MStructuredGrid::interpolateValue(QVector3D vec3_lonLatP)
}
float MStructuredGrid::interpolateValueOnLevel(
float lon, float lat, unsigned int k)
{
int i, j, i1, j1;
float mixI, mixJ;
findEnclosingHorizontalIndices(lon, lat, &i, &j, &i1, &j1, &mixI, &mixJ);
if ((i < 0) || (j < 0) || (i1 >= int(nlons)) || (j1 >= int(nlats)))
return M_MISSING_VALUE;
// Get scalar values at the four surrounding grid points of the level.
float scalar_i0j0 = getValue(k, j , i);
float scalar_i1j0 = getValue(k, j , i1);
float scalar_i0j1 = getValue(k, j1, i);
float scalar_i1j1 = getValue(k, j1, i1);
// Interpolate horizontally.
mixJ = MFRACT(mixJ); // fract(mixJ) in GLSL
float scalar_i0 = MMIX(scalar_i0j0, scalar_i0j1, mixJ);
float scalar_i1 = MMIX(scalar_i1j0, scalar_i1j1, mixJ);
mixI = MFRACT(mixI);
float scalar = MMIX(scalar_i0, scalar_i1, mixI);
return scalar;
}
QVector<QVector2D> MStructuredGrid::extractVerticalProfile(float lon, float lat)
{
QVector<QVector2D> profile;
for (unsigned int k = 0; k < nlevs; k++)
{
float scalar = interpolateValueOnLevel(lon, lat, k);
float pressure_hPa = levelPressureAtLonLat_hPa(lon, lat, k);
profile.append(QVector2D(scalar, pressure_hPa));
}
return profile;
}
bool MStructuredGrid::findTopGridIndices(
float lon, float lat, float p_hPa,
MIndex3D *nw, MIndex3D *ne, MIndex3D *sw, MIndex3D *se)
{
// Find horizontal indices i,i+1 and j,j+1 that enclose lon, lat.
float mixI = MMOD(lon - lons[0], 360.) / abs(lons[1]-lons[0]);
float mixJ = (lats[0] - lat) / abs(lats[1]-lats[0]);
int i = int(mixI);
int j = int(mixJ);
int i1 = i+1;
if (gridIsCyclicInLongitude()) i1 %= nlons;
int j1 = j+1;
int i, j, i1, j1;
float mixI, mixJ;
findEnclosingHorizontalIndices(lon, lat, &i, &j, &i1, &j1, &mixI, &mixJ);
nw->i = i; nw->j = j; nw->k = findLevel(nw->j, nw->i, p_hPa);
ne->i = i1; ne->j = j; ne->k = findLevel(ne->j, ne->i, p_hPa);
......@@ -1100,6 +1150,15 @@ float MRegularLonLatLnPGrid::interpolateGridColumnToPressure(
}
float MRegularLonLatLnPGrid::levelPressureAtLonLat_hPa(
float lon, float lat, unsigned int k)
{
Q_UNUSED(lon);
Q_UNUSED(lat);
return exp(levels[k]);
}
int MRegularLonLatLnPGrid::findLevel(
unsigned int j, unsigned int i, float p_hPa)
{
......@@ -1305,6 +1364,15 @@ float MRegularLonLatStructuredPressureGrid::interpolateGridColumnToPressure(
}
float MRegularLonLatStructuredPressureGrid::levelPressureAtLonLat_hPa(
float lon, float lat, unsigned int k)
{
Q_UNUSED(lon);
Q_UNUSED(lat);
return levels[k];
}
int MRegularLonLatStructuredPressureGrid::findLevel(
unsigned int j, unsigned int i, float p_hPa)
{
......@@ -1835,6 +1903,18 @@ float MLonLatHybridSigmaPressureGrid::interpolateGridColumnToPressure(
}
float MLonLatHybridSigmaPressureGrid::levelPressureAtLonLat_hPa(
float lon, float lat, unsigned int k)
{
// Interpolate surface pressure to lon/lat position (pressure value is
// ignored by interpolateValue() since surface pressure is a 2D field), ..
float psfc_hPa = surfacePressure->interpolateValue(lon, lat, 0.) / 100.;
// .. then compute pressure of level.
float p_hPa = ak_hPa[k] + bk[k] * psfc_hPa;
return p_hPa;
}
int MLonLatHybridSigmaPressureGrid::findLevel(
unsigned int j, unsigned int i, float p_hPa)
{
......@@ -2120,6 +2200,13 @@ float MLonLatAuxiliaryPressureGrid::interpolateGridColumnToPressure(
}
float MLonLatAuxiliaryPressureGrid::levelPressureAtLonLat_hPa(
float lon, float lat, unsigned int k)
{
return auxPressureField_hPa->interpolateValueOnLevel(lon, lat, k);
}
int MLonLatAuxiliaryPressureGrid::findLevel(
unsigned int j, unsigned int i, float p_hPa)
{
......
......@@ -4,10 +4,13 @@
** three-dimensional visual exploration of numerical ensemble weather
** prediction data.
**
** Copyright 2015-2018 Marc Rautenhaus
** Copyright 2017-2018 Bianca Tost
** Copyright 2015-2019 Marc Rautenhaus [*, previously +]
** Copyright 2017-2018 Bianca Tost [+]
**
** Computer Graphics and Visualization Group
** * Regional Computing Center, Visualization
** Universitaet Hamburg, Hamburg, Germany
**
** + Computer Graphics and Visualization Group
** Technische Universitaet Muenchen, Garching, Germany
**
** Met.3D is free software: you can redistribute it and/or modify
......@@ -275,8 +278,17 @@ public:
inline float getSouthInterfaceLat(unsigned int j)
{ return lats[j] - getDeltaLat()/2.; }
/**
Determine the horizontal grid indices @p i, @p j, @p i1, @p j1 that
enclose the position given by @p lon, @p lat.
*/
void findEnclosingHorizontalIndices(float lon, float lat, int *i, int *j,
int *i1, int *j1, float *mixI, float *mixJ);
/**
Sample the data grid at lon, lat and p, using trilinear interpolation.
Uses @ref interpolateGridColumnToPressure(). For derived grid classes
that are only two-dimensional, the @p p_hPa parameter is ignored.
*/
float interpolateValue(float lon, float lat, float p_hPa);
......@@ -284,12 +296,34 @@ public:
/**
Implement this method in derived classes that know about their vertical
coordinate. It is used by @ref interpolateValue().
coordinate. It is used by @ref interpolateValue(). If the derived class
is two-dimensional, the @p p_hPa parameter can be ignored.
*/
virtual float interpolateGridColumnToPressure(
unsigned int j, unsigned int i, float p_hPa)
{ Q_UNUSED(j); Q_UNUSED(i); Q_UNUSED(p_hPa); return M_MISSING_VALUE; }
/**
Computes the pressure on grid level @p k at position (@p lon, @p lat).
Implement this method in derived classes.
*/
virtual float levelPressureAtLonLat_hPa(float lon, float lat, unsigned int k)
{ Q_UNUSED(lon); Q_UNUSED(lat); Q_UNUSED(k); return M_MISSING_VALUE; }
/**
Samples the data grid on vertical level k and at position (@p lon, @p lat)
using bi-linear interpolation.
*/
float interpolateValueOnLevel(float lon, float lat, unsigned int k);
/**
Extracts a vertical profile of (scalar, p_hPa) tuples from the data
field at position (@p lon, @p lat). Uses @ref interpolateValueOnLevel()
and @ref levelPressureAtLonLat_hPa().
*/
QVector<QVector2D> extractVerticalProfile(float lon, float lat);
/**
Determine the four grid indices that horizontally bound the grid cell
that contains the position specified by @p lon, @p lat, @p p_hPa. In the
......@@ -598,6 +632,8 @@ public:
float interpolateGridColumnToPressure(unsigned int j, unsigned int i,
float p_hPa);
float levelPressureAtLonLat_hPa(float lon, float lat, unsigned int k) override;
int findLevel(unsigned int j, unsigned int i, float p_hPa);
float getPressure(unsigned int k, unsigned int j, unsigned int i);
......@@ -636,6 +672,8 @@ public:
float interpolateGridColumnToPressure(unsigned int j, unsigned int i,
float p_hPa);
float levelPressureAtLonLat_hPa(float lon, float lat, unsigned int k) override;
int findLevel(unsigned int j, unsigned int i, float p_hPa);
float getPressure(unsigned int k, unsigned int j, unsigned int i);
......@@ -669,6 +707,14 @@ public:
inline float getValue(unsigned int j, unsigned int i) const
{ return data[INDEX2yx(j, i, nlons)]; }
/**
2D special case: ignore @p p_hPa parameter and simply map to getValue().
Implementation required for @ref MStructuredGrid::interpolateValue().
*/
float interpolateGridColumnToPressure(unsigned int j, unsigned int i,
float p_hPa)
{ Q_UNUSED(p_hPa); return getValue(j, i); }
GL::MTexture* getTexture(QGLWidget *currentGLContext = nullptr,
bool nullTexture = false);
......@@ -715,6 +761,8 @@ public:
float interpolateGridColumnToPressure(unsigned int j, unsigned int i,
float p_hPa);
float levelPressureAtLonLat_hPa(float lon, float lat, unsigned int k) override;
int findLevel(unsigned int j, unsigned int i, float p_hPa);
float getPressure(unsigned int k, unsigned int j, unsigned int i);
......@@ -781,6 +829,8 @@ public:
float interpolateGridColumnToPressure(unsigned int j, unsigned int i,
float p_hPa);
float levelPressureAtLonLat_hPa(float lon, float lat, unsigned int k) override;
int findLevel(unsigned int j, unsigned int i, float p_hPa);
float getPressure(unsigned int k, unsigned int j, unsigned int i);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment