Skip to content

Commit

Permalink
Adds bilinear interpolation as an option in -d gradient
Browse files Browse the repository at this point in the history
  • Loading branch information
@tnipen
tnipen committed Nov 20, 2017
1 parent f3d10cc commit 954950b
Show file tree
Hide file tree
Showing 4 changed files with 118 additions and 21 deletions.
54 changes: 33 additions & 21 deletions src/Downscaler/Gradient.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ DownscalerGradient::DownscalerGradient(Variable::Type iVariable, const Options&
mDefaultGradient(0),
mMinElevDiff(30),
mAverageNeighbourhood(false),
mDownscalerName("nearestNeighbour"),
mHasIssuedWarningUnstable(false) {

iOptions.getValue("minGradient", mMinGradient);
Expand All @@ -23,6 +24,7 @@ DownscalerGradient::DownscalerGradient(Variable::Type iVariable, const Options&
iOptions.getValue("constantGradient", mConstantGradient);
iOptions.getValue("minElevDiff", mMinElevDiff);
iOptions.getValue("averageNeighbourhood", mAverageNeighbourhood);
iOptions.getValue("downscaler", mDownscalerName);
}

void DownscalerGradient::downscaleCore(const File& iInput, File& iOutput) const {
Expand All @@ -49,6 +51,17 @@ void DownscalerGradient::downscaleCore(const File& iInput, File& iOutput) const
Field& ifield = *iInput.getField(mVariable, t);
Field& ofield = *iOutput.getField(mVariable, t);

vec2 elevsInterp;
vec2 lafsInterp;
if(mDownscalerName == "nearestNeighbour") {
elevsInterp = DownscalerNearestNeighbour::downscaleVec(ielevs, ilats, ilons, olats, olons, nearestI, nearestJ);
DownscalerNearestNeighbour::downscaleField(ifield, ofield, ilats, ilons, olats, olons, nearestI, nearestJ);
}
else if (mDownscalerName == "bilinear") {
elevsInterp = DownscalerBilinear::downscaleVec(ielevs, ilats, ilons, olats, olons, nearestI, nearestJ);
DownscalerBilinear::downscaleField(ifield, ofield, ilats, ilons, olats, olons, nearestI, nearestJ);
}

#pragma omp parallel for
for(int i = 0; i < nLat; i++) {
for(int j = 0; j < nLon; j++) {
Expand All @@ -58,35 +71,34 @@ void DownscalerGradient::downscaleCore(const File& iInput, File& iOutput) const
assert(Jcenter < ielevs[Icenter].size());
for(int e = 0; e < nEns; e++) {
float currElev = oelevs[i][j];
float nearestElev = ielevs[Icenter][Jcenter];
if(!Util::isValid(currElev) || !Util::isValid(nearestElev)) {
// Can't adjust if we don't have an elevation, use nearest neighbour
ofield(i,j,e) = ifield(Icenter,Jcenter,e);
float baseElev = elevsInterp[i][j];
float baseValue = ofield(i, j, e);
if(!Util::isValid(currElev) || !Util::isValid(baseElev)) {
// Can't adjust if we don't have an elevation, use uncorrected value
}
else {
float gradient = mDefaultGradient;
float totalValue = 0;
float totalElev = 0;
int counter = 0;
int averagingRadius = 0;
if(mAverageNeighbourhood)
averagingRadius = mSearchRadius;
for(int ii = std::max(0, Icenter-averagingRadius); ii <= std::min(iInput.getNumLat()-1, Icenter+averagingRadius); ii++) {
for(int jj = std::max(0, Jcenter-averagingRadius); jj <= std::min(iInput.getNumLon()-1, Jcenter+averagingRadius); jj++) {
float currValue = ifield(ii,jj,e);
float currElev = ielevs[ii][jj];
if(Util::isValid(currValue) && Util::isValid(currElev)) {
totalValue += currValue;
totalElev += currElev;
counter++;

// Compute the average elevation and value within the neighbourhood
if(mAverageNeighbourhood) {
for(int ii = std::max(0, Icenter-mSearchRadius); ii <= std::min(iInput.getNumLat()-1, Icenter+mSearchRadius); ii++) {
for(int jj = std::max(0, Jcenter-mSearchRadius); jj <= std::min(iInput.getNumLon()-1, Jcenter+mSearchRadius); jj++) {
float currValue = ifield(ii,jj,e);
float currElev = ielevs[ii][jj];
if(Util::isValid(currValue) && Util::isValid(currElev)) {
totalValue += currValue;
totalElev += currElev;
counter++;
}
}
}
}
float baseValue = Util::MV;
float baseElev = Util::MV;
if(counter > 0) {
baseValue = totalValue / counter;
baseElev = totalElev / counter;
if(counter > 0) {
baseValue = totalValue / counter;
baseElev = totalElev / counter;
}
}
float dElev = currElev - baseElev;

Expand Down
1 change: 1 addition & 0 deletions src/Downscaler/Gradient.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,5 +34,6 @@ class DownscalerGradient : public Downscaler {
float mDefaultGradient;
bool mAverageNeighbourhood;
mutable bool mHasIssuedWarningUnstable;
std::string mDownscalerName;
};
#endif
70 changes: 70 additions & 0 deletions src/Downscaler/NearestNeighbour.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,76 @@ void DownscalerNearestNeighbour::downscaleCore(const File& iInput, File& iOutput
}
}
}
void DownscalerNearestNeighbour::downscaleField(const Field& iInput, Field& iOutput,
const vec2& iInputLats, const vec2& iInputLons,
const vec2& iOutputLats, const vec2& iOutputLons,
const vec2Int& nearestI, const vec2Int& nearestJ) {

int nLat = iOutput.getNumLat();
int nLon = iOutput.getNumLon();
int nEns = iOutput.getNumEns();

#pragma omp parallel for
for(int i = 0; i < nLat; i++) {
for(int j = 0; j < nLon; j++) {
int I = nearestI[i][j];
int J = nearestJ[i][j];
for(int e = 0; e < nEns; e++) {
if(Util::isValid(I) && Util::isValid(J)) {
iOutput(i,j,e) = iInput(I, J, e);
}
else
iOutput(i,j,e) = Util::MV;
}
}
}
}

vec2 DownscalerNearestNeighbour::downscaleVec(const vec2& iInput,
const vec2& iInputLats, const vec2& iInputLons,
const vec2& iOutputLats, const vec2& iOutputLons,
const vec2Int& nearestI, const vec2Int& nearestJ) {

int nLat = iOutputLats.size();
int nLon = iOutputLats[0].size();

assert(nearestI.size() == iOutputLats.size());
assert(nearestJ.size() == iOutputLats.size());
assert(nearestI[0].size() == iOutputLats[0].size());
assert(nearestJ[0].size() == iOutputLats[0].size());

vec2 output;
output.resize(nLat);
for(int i = 0; i < nLat; i++)
output[i].resize(nLon);

#pragma omp parallel for
for(int i = 0; i < nLat; i++) {
for(int j = 0; j < nLon; j++) {
assert(i < nearestI.size());
assert(j < nearestI[i].size());
assert(i < nearestJ.size());
assert(j < nearestJ[i].size());
int I = nearestI[i][j];
int J = nearestJ[i][j];
assert(I < iInput.size());
assert(J < iInput[I].size());
assert(I >= 0);
assert(J >= 0);
assert(i < output.size());
assert(j < output[i].size());
assert(i >= 0);
assert(j >= 0);
if(Util::isValid(I) && Util::isValid(J)) {
output[i][j] = iInput[I][J];
}
else
output[i][j] = Util::MV;
}
}
return output;
}


std::string DownscalerNearestNeighbour::description() {
std::stringstream ss;
Expand Down
14 changes: 14 additions & 0 deletions src/Downscaler/NearestNeighbour.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,24 @@
#include "Downscaler.h"
#include "../Variable.h"
#include "../Util.h"
#include "../Field.h"
class File;
class DownscalerNearestNeighbour : public Downscaler {
public:
DownscalerNearestNeighbour(Variable::Type iVariable, const Options& iOptions);

// Interpolate a whole field
static void downscaleField(const Field& iInput, Field& iOutput,
const vec2& iInputLats, const vec2& iInputLons,
const vec2& iOutputLats, const vec2& iOutputLons,
const vec2Int& nearestI, const vec2Int& nearestJ);

// Interpolate a whole vec2
static vec2 downscaleVec(const vec2& iInput,
const vec2& iInputLats, const vec2& iInputLons,
const vec2& iOutputLats, const vec2& iOutputLons,
const vec2Int& nearestI, const vec2Int& nearestJ);

static std::string description();
std::string name() const {return "nearestNeighbour";};
private:
Expand Down

0 comments on commit 954950b

Please sign in to comment.