diff --git a/src/Downscaler/Gradient.cpp b/src/Downscaler/Gradient.cpp index 287d121b..14567cbc 100644 --- a/src/Downscaler/Gradient.cpp +++ b/src/Downscaler/Gradient.cpp @@ -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); @@ -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 { @@ -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++) { @@ -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; diff --git a/src/Downscaler/Gradient.h b/src/Downscaler/Gradient.h index df318bc9..587a9649 100644 --- a/src/Downscaler/Gradient.h +++ b/src/Downscaler/Gradient.h @@ -34,5 +34,6 @@ class DownscalerGradient : public Downscaler { float mDefaultGradient; bool mAverageNeighbourhood; mutable bool mHasIssuedWarningUnstable; + std::string mDownscalerName; }; #endif diff --git a/src/Downscaler/NearestNeighbour.cpp b/src/Downscaler/NearestNeighbour.cpp index f84cee01..23e7c8d8 100644 --- a/src/Downscaler/NearestNeighbour.cpp +++ b/src/Downscaler/NearestNeighbour.cpp @@ -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; diff --git a/src/Downscaler/NearestNeighbour.h b/src/Downscaler/NearestNeighbour.h index 4a6850e9..f8a98a90 100644 --- a/src/Downscaler/NearestNeighbour.h +++ b/src/Downscaler/NearestNeighbour.h @@ -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: