diff --git a/algorithms/java/src/org/openda/algorithms/kalmanFilter/AbstractSequentialAlgorithm.java b/algorithms/java/src/org/openda/algorithms/kalmanFilter/AbstractSequentialAlgorithm.java index d8e3bb860..bd27ce285 100644 --- a/algorithms/java/src/org/openda/algorithms/kalmanFilter/AbstractSequentialAlgorithm.java +++ b/algorithms/java/src/org/openda/algorithms/kalmanFilter/AbstractSequentialAlgorithm.java @@ -34,6 +34,7 @@ import java.io.FileFilter; import java.text.ParseException; import java.util.HashMap; +import java.util.LinkedHashMap; import static org.openda.utils.performance.OdaGlobSettings.getTimePrecision; @@ -77,9 +78,9 @@ public abstract class AbstractSequentialAlgorithm extends Instance implements IA // Kalman preparations // - the actual analysis is performed by the filter - protected HashMap gainVectors = new HashMap(); - protected HashMap obsId = new HashMap(); - protected HashMap obsTimeOffset = new HashMap(); + protected HashMap gainVectors = new LinkedHashMap<>(); + protected HashMap obsId = new LinkedHashMap<>(); + protected HashMap obsTimeOffset = new LinkedHashMap<>(); //collecting data for writing output diff --git a/algorithms/java/src/org/openda/algorithms/kalmanFilter/DEnKF.java b/algorithms/java/src/org/openda/algorithms/kalmanFilter/DEnKF.java index e788cd7d1..8853f70e0 100644 --- a/algorithms/java/src/org/openda/algorithms/kalmanFilter/DEnKF.java +++ b/algorithms/java/src/org/openda/algorithms/kalmanFilter/DEnKF.java @@ -19,6 +19,7 @@ */ package org.openda.algorithms.kalmanFilter; import org.openda.interfaces.*; +import org.openda.utils.Matrix; import org.openda.utils.Results; /** @@ -101,8 +102,9 @@ public void analysis(IStochObserver obs, IVector obsValues, IVector predictions, // Apply localization to the gain matrix (shur-product) applyLocalizationToGain(obs, Kvecs, ensemblePredictionsForecast, ensembleVectorsForecast); + double[][] hk = getHK(obs, ensemblePredictionsForecast); // Store kalman gain for future use in this object - storeGainMatrix(obs, analysisTime, Kvecs); + storeGainMatrix(obs, analysisTime, Kvecs, hk); // Multiply Kalman gain with innovations and update model states updateModelWithGain(obs, ensemblePredictionsForecast, ensembleVectorsForecast, Kvecs); diff --git a/algorithms/java/src/org/openda/algorithms/kalmanFilter/EnKF.java b/algorithms/java/src/org/openda/algorithms/kalmanFilter/EnKF.java index 9381d25d6..9fe1a93f1 100644 --- a/algorithms/java/src/org/openda/algorithms/kalmanFilter/EnKF.java +++ b/algorithms/java/src/org/openda/algorithms/kalmanFilter/EnKF.java @@ -23,13 +23,16 @@ import org.openda.observers.ObserverUtils; import org.openda.utils.Matrix; import org.openda.utils.Results; +import org.openda.utils.Vector; import org.openda.utils.io.FileBasedModelState; import org.openda.utils.io.KalmanGainStorage; import org.openda.utils.performance.OdaTiming; import java.io.File; import java.io.FileFilter; +import java.util.Arrays; import java.util.HashMap; +import java.util.List; /** * @author Nils van Velzen based in the initial OpenDA implementation by Martin Verlaan @@ -423,6 +426,7 @@ protected IVector[] computeGainMatrix(IStochObserver obs, EnsembleVectors ensemb // are not upadted until after the next forecast // H*K = PRED*PRED'*inv(D) timerLinalg.start(); + // This is H*K and should also be written to Kalman Gain Storage Matrix K_pred = new Matrix(m,m); K_pred.multiply(1.0, predMat, E, 0.0, false, false); // pred_a_linear = predAvg + K_pred*(obsVal-predAvg) @@ -444,9 +448,52 @@ protected IVector[] computeGainMatrix(IStochObserver obs, EnsembleVectors ensemb return Kvecs; } + protected Matrix computeHK(IStochObserver obs, EnsembleVectors ensemblePredictions){ + int m = obs.getCount(); // number of observations + int q = this.ensembleSize; // number of ensemble members + + // compute Kalman gain + // D = HPH+R = (1/(q-1))PRED*PRED'+sqrtR*sqrtR' : covariance of + // innovations + timerLinalg.start(); + //H*A^f_k = predMat = prediction of the observed model values, after removing the mean. + Matrix predMat = new Matrix(ensemblePredictions.ensemble); + predMat.scale(Math.sqrt(1.0 / (q - 1.0))); + // System.out.println("predMat="+predMat); + Matrix D = new Matrix(m, m); + D.multiply(1.0, predMat, predMat, 0.0, false, true); + IMatrix sqrtR = obs.getSqrtCovariance().asMatrix(); + D.multiply(1.0, sqrtR, sqrtR, 1.0, false, true); + + // System.out.println("PHT="+K); + Matrix inverseD = D.inverse(); + // System.out.println("inverseD="+inverseD); + + // version without large matrices + // K = XI * E with E=PRED'*inv(D) + Matrix E = new Matrix(q, m); + E.multiply(1.0, predMat, inverseD, 0.0, true, false); + + // Compute H*K for linear update of predictions, since for blackbox models the predictions + // are not upadted until after the next forecast + // H*K = PRED*PRED'*inv(D) + // This is H*K and should also be written to Kalman Gain Storage + Matrix K_pred = new Matrix(m, m); + K_pred.multiply(1.0, predMat, E, 0.0, false, false); + // pred_a_linear = predAvg + K_pred*(obsVal-predAvg) + IVector innovAvg = obs.getExpectations(); + innovAvg.axpy(-1, ensemblePredictions.mean); + IVector pred_a_linear = ensemblePredictions.mean.clone(); + K_pred.rightMultiply(1.0, innovAvg, 1.0, pred_a_linear); + innovAvg.free(); + timerLinalg.stop(); + + pred_a_linear.free(); + return K_pred; + } - protected void storeGainMatrix(IStochObserver obs, ITime analysisTime, IVector[] Kvecs){ + protected void storeGainMatrix(IStochObserver obs, ITime analysisTime, IVector[] Kvecs, double[][] hkArray){ // store kalman gain for future use in this object if(this.saveGainTimes!=null){ @@ -488,6 +535,7 @@ protected void storeGainMatrix(IStochObserver obs, ITime analysisTime, IVector[] +"model = "+this.mainModel.getClass().getSimpleName()+"\n" +"observer = "+obs.getClass().getSimpleName()+"\n"); + gainStorage.setHk(hkArray); gainStorage.writeKalmanGain(this.gainVectors, this.obsId, this.obsTimeOffset); this.gainVectors.clear(); this.obsId.clear(); @@ -576,8 +624,9 @@ public void analysis(IStochObserver obs, IVector obsValues, IVector predictions, this.smoothedGainMatrix.SmoothGain(obs,Kvecs, this.timeRegularisationPerDay, analysisTime); } + double[][] hk = getHK(obs, ensemblePredictionsForecast); // Store kalman gain for future use in this object - storeGainMatrix(obs, analysisTime, Kvecs); + storeGainMatrix(obs, analysisTime, Kvecs, hk); // Multiply Kalman gain with innovations and update model states updateModelWithGain(obs, ensemblePredictionsForecast, ensembleVectorsForecast, Kvecs); @@ -605,6 +654,27 @@ public void analysis(IStochObserver obs, IVector obsValues, IVector predictions, System.gc(); } + protected double[][] getHK(IStochObserver obs, EnsembleVectors ensemblePredictionsForecast) { + IVector[] rhoForLocalization = getRhoForLocalization(obs); + if (rhoForLocalization == null) return null; + Matrix hk = computeHK(obs, ensemblePredictionsForecast); + double[][] array = new double[hk.getNumberOfRows()][hk.getNumberOfRows()]; + for (int i = 0; i < array.length; i++) { + for (int j = 0; j < array[0].length; j++) { + array[i][j] = rhoForLocalization[i].getValue(j) * hk.getValue(i, j); + } + } + return array; + } + + private IVector[] getRhoForLocalization(IStochObserver obs) { + IObservationDescriptions observationDescriptions = obs.getObservationDescriptions(); + if (localizationMethod == LocalizationMethodType.hamill) return this.ensemble[0].getRhoForLocalization(observationDescriptions, this.distance); + if (localizationMethod == LocalizationMethodType.none) return IModelInstance.getDefaultRhoForLocalization(observationDescriptions); + System.out.println("Not computing HK for Kalman Gain Storage because it is not compatible with localization method " + localizationMethod); + return null; + } + @Override protected void restoreGain(File[] tempFiles, File restartTempDir) { // restore the last Kalman gain, in case of EnKF-GS diff --git a/algorithms/java/src/org/openda/algorithms/kalmanFilter/LocEnKF.java b/algorithms/java/src/org/openda/algorithms/kalmanFilter/LocEnKF.java index d908d23b7..12fd546d7 100644 --- a/algorithms/java/src/org/openda/algorithms/kalmanFilter/LocEnKF.java +++ b/algorithms/java/src/org/openda/algorithms/kalmanFilter/LocEnKF.java @@ -19,11 +19,8 @@ */ package org.openda.algorithms.kalmanFilter; import org.openda.interfaces.*; -import org.openda.observers.ObserverUtils; import org.openda.utils.Matrix; import org.openda.utils.Results; -import org.openda.utils.io.KalmanGainStorage; -import org.openda.utils.performance.OdaTiming; /** * @author Nils van Velzen based in the initial OpenDA implementation by Martin Verlaan @@ -173,8 +170,9 @@ public void analysis(IStochObserver obs, IVector obsValues, IVector predictions, this.smoothedGainMatrix.SmoothGain(obs,Kvecs, this.timeRegularisationPerDay, analysisTime); } + double[][] hk = getHK(obs, ensemblePredictionsForecast); // Store kalman gain for future use in this object - storeGainMatrix(lobs, analysisTime, Kvecs); + storeGainMatrix(lobs, analysisTime, Kvecs, hk); // Multiply Kalman gain with innovations and update model states updateModelWithGain(lobs, ensemblePredictionsForecast, ensembleVectorsForecast, Kvecs); diff --git a/algorithms/java/src/org/openda/algorithms/kalmanFilter/SteadyStateFilter.java b/algorithms/java/src/org/openda/algorithms/kalmanFilter/SteadyStateFilter.java index 4df95e23c..bfb9abe1c 100644 --- a/algorithms/java/src/org/openda/algorithms/kalmanFilter/SteadyStateFilter.java +++ b/algorithms/java/src/org/openda/algorithms/kalmanFilter/SteadyStateFilter.java @@ -22,12 +22,15 @@ import org.openda.interfaces.*; import org.openda.observers.ObserverUtils; import org.openda.utils.ConfigTree; +import org.openda.utils.Matrix; import org.openda.utils.Results; +import org.openda.utils.SortUtils; import org.openda.utils.io.KalmanGainStorage; import org.openda.utils.performance.OdaGlobSettings; import java.io.File; import java.text.ParseException; +import java.util.*; /** * Steady state Kalman filter implementation as introduced already in the original papers of Kalman and Bucy. @@ -54,6 +57,10 @@ public class SteadyStateFilter extends AbstractSequentialAlgorithm { private double[] readGainTime; private int steadyStateTimeCounter = 0; private double skipAssimilationStandardDeviationFactor = Double.POSITIVE_INFINITY; + private boolean estimateMissingObservations = false; + private double[][] hk; + private String[] gainStorageObservationIdsArray; + private String[] gainVectorIdArray; public void initialize(File workingDir, String[] arguments) { super.initialize(workingDir, arguments); @@ -151,7 +158,8 @@ public void initialize(File workingDir, String[] arguments) { } } this.skipAssimilationStandardDeviationFactor = this.configurationAsTree.getAsDouble("skipAssimilationStandardDeviationFactor", Double.POSITIVE_INFINITY); - + this.estimateMissingObservations = this.configurationAsTree.getAsBoolean("estimateMissingObservations", false); + //now read the first gain file into memory steadyStateTimeCounter++; KalmanGainStorage gainStorage = new KalmanGainStorage(this.workingDir, this.gainTimeMjd[0]); @@ -159,14 +167,15 @@ public void initialize(File workingDir, String[] arguments) { gainStorage.setKalmanGainStorageFileName(this.gainFileName); if (gainFileName.endsWith(".nc")) gainStorage.setColumnFileType(KalmanGainStorage.StorageType.netcdf_cf); gainStorage.readKalmanGain( this.getCurrentState()); - String[] obsIds = gainStorage.getObservationIds(); - double[] obsTimeOffsets = gainStorage.getObservationOffsetInDays(); + hk = gainStorage.getHk(); + gainStorageObservationIdsArray = gainStorage.getObservationIds(); + double[] gainStorageObsTimeOffsets = gainStorage.getObservationOffsetInDays(); IVector[] KVecs = gainStorage.getKalmanGainColumns(); - for(int i=0;i skipAssimilationStandardDeviationFactor * standardDeviations.getValue(i)) { - Results.putProgression("Info: Skipping assimilation for " + gainVectorId + " because innovation > (skipAssimilationStandardDeviationFactor * obs standard deviation). Observed value = " + observationValues.getValue(i) + ", model prediction value " + pred_a.getValue(i) + ", skipAssimilationStandardDeviationFactor = "+ skipAssimilationStandardDeviationFactor + ", obs stdv = " + standardDeviations.getValue(i) +"\n"); + Results.putProgression("Info: Skipping assimilation for " + gainVectorId + " because innovation > (skipAssimilationStandardDeviationFactor * obs standard deviation). Observed value = " + observationValues.getValue(i) + ", model prediction value " + pred_a.getValue(i) + ", skipAssimilationStandardDeviationFactor = " + skipAssimilationStandardDeviationFactor + ", obs stdv = " + standardDeviations.getValue(i) + "\n"); continue; } + System.out.println("Innovation value for " + obsIds[i] + ": " + innovationValue); delta.axpy(innovationValue, gainVector); - }else{ - throw new RuntimeException("No matching column found for observation with id=" - +obsIds[i]+"and offset="+obsTimeOffsets[i]+"\n"); + if (hkCalculator != null && hkCalculator.hasObservationsAvailable()) hkCalculator.setDAvailableValue(gainVectorId, innovationValue); + } else { + + if (hk == null) throw new RuntimeException("No matching column found for observation with id=" + obsIds[i] + "and offset=" + obsTimeOffsets[i] + "\n"); } } + System.out.println(); + + if (hkCalculator != null && hkCalculator.hasMissingObservations()) hkCalculator.compensateForMissingObservationsWithHK(obsTimeOffsets, delta); this.mainModel.axpyOnState(1.0, delta); Results.putValue("pred_a", pred_a, pred_a.getSize(), "analysis step", IResultWriter.OutputLevel.Essential, IResultWriter.MessageType.Step); } - private boolean isReplaceGainAtThisTime(ITime time){ + private static String[] getAvailableObsGainVectorIdArray(String[] obsIds, double[] obsTimeOffsets) { + String[] availableObsGainVectorIdArray = new String[obsIds.length]; + for(int i = 0; i< obsIds.length; i++){ + String gainVectorId = obsIds[i]+":"+Math.round(obsTimeOffsets[i]*24.0*3600.0); //conversion to seconds + availableObsGainVectorIdArray[i] = gainVectorId; + } + return availableObsGainVectorIdArray; + } + + private static Matrix createIMinusM2(int numberOfMissingObservations, Matrix m2) { + Matrix iMinusM2 = new Matrix(numberOfMissingObservations, numberOfMissingObservations); + for (int row = 0; row < numberOfMissingObservations; row++) { + for (int column = 0; column < numberOfMissingObservations; column++) { + double m2Value = m2.getValue(row, column); + double iMinusM2Value = row == column ? 1 - m2Value : -m2Value; + iMinusM2.setValue(row, column, iMinusM2Value); + } + } + return iMinusM2; + } + + private void fillM1M2Matrices(ArrayList missingObservationIndices, ArrayList availableObservationIndicesGainStorage, Matrix m1, Matrix m2) { + for (int rowIndex = 0; rowIndex < missingObservationIndices.size(); rowIndex++) { + Integer missingObservationRowIndex = missingObservationIndices.get(rowIndex); + for (int columnIndex = 0; columnIndex < availableObservationIndicesGainStorage.size(); columnIndex++) { + Integer availableObservationColumnIndex = availableObservationIndicesGainStorage.get(columnIndex); + m1.setValue(rowIndex, columnIndex, hk[missingObservationRowIndex][availableObservationColumnIndex]); + } + for (int columnIndex = 0; columnIndex < missingObservationIndices.size(); columnIndex++) { + Integer missingObservationColumnIndex = missingObservationIndices.get(columnIndex); + m2.setValue(rowIndex, columnIndex, hk[missingObservationRowIndex][missingObservationColumnIndex]); + } + } + } + + private boolean isReplaceGainAtThisTime(ITime time){ boolean result=false; if(this.readGainTime!=null){ for(int i=0;i gainVectors; + ArrayList missingObservationIndices = new ArrayList<>(); + ArrayList availableObservationIndicesGainStorage = new ArrayList<>(); + ArrayList availableObservationIndicesObserver = new ArrayList<>(); + private int numberOfMissingObservations; + private int numberOfAvailableObservations; + private Matrix m1; + private Matrix m2; + private Matrix dAvailable; + + public HKCalculator(String[] availableObsGainVectorIdArray, HashMap gainVectors) { + this.availableObsGainVectorIdArray = availableObsGainVectorIdArray; + this.gainVectors = gainVectors; + createMatrices(); + } + + private void fillIndicesArrays() { + for (int i = 0; i < gainVectorIdArray.length; i++) { + String gainVectorId = gainVectorIdArray[i]; + int indexOf = SortUtils.indexOfString(this.availableObsGainVectorIdArray, gainVectorId); + if (indexOf >= 0) { + this.availableObservationIndicesGainStorage.add(i); + this.availableObservationIndicesObserver.add(indexOf); + continue; + } + System.out.printf("Observation for %s missing%n", gainVectorId); + this.missingObservationIndices.add(i); + } + } + + public void createMatrices() { + fillIndicesArrays(); + numberOfMissingObservations = missingObservationIndices.size(); + numberOfAvailableObservations = gainStorageObservationIdsArray.length - numberOfMissingObservations; + m1 = new Matrix(numberOfMissingObservations, numberOfAvailableObservations); + m2 = new Matrix(numberOfMissingObservations, numberOfMissingObservations); + fillM1M2Matrices(missingObservationIndices, availableObservationIndicesGainStorage, m1, m2); + dAvailable = new Matrix(this.numberOfAvailableObservations, 1); + } + + public boolean hasObservationsAvailable() { + return this.numberOfAvailableObservations > 0; + } + + public void setDAvailableValue(String gainVectorId, double innovationValue) { + int indexOf = SortUtils.indexOfString(availableObsGainVectorIdArray, gainVectorId); + int indexOfDAvailable = availableObservationIndicesObserver.indexOf(indexOf); + System.out.println("Index in dAvailable " + indexOfDAvailable); + System.out.println(); + dAvailable.setValue(indexOfDAvailable, 0, innovationValue); + } + + public boolean hasMissingObservations() { + return this.numberOfMissingObservations > 0; + } + + private void compensateForMissingObservationsWithHK(double[] obsTimeOffsets, IVector delta) { + Matrix iMinusM2 = createIMinusM2(this.numberOfMissingObservations, this.m2); + System.out.println("Trying to compensate for missing observations using HK from kalman gain storage"); + Matrix inverseIMinusM2 = iMinusM2.inverse(); + Matrix m1DAvailable = this.m1.mult(this.dAvailable); + Matrix dMissing = inverseIMinusM2.mult(m1DAvailable); + /*System.out.println("M1: " + m1); + System.out.println("M2: " + m2); + System.out.println("inverseIMinusM2: " + inverseIMinusM2); + System.out.println("dAvailable: " + dAvailable); + System.out.println("m1DAvailable: " + m1DAvailable); + System.out.println("dMissing: " + dMissing);*/ + for (int i = 0; i < this.missingObservationIndices.size(); i++) { + Integer missingObservationIndex = this.missingObservationIndices.get(i); + String missingObservationId = gainStorageObservationIdsArray[missingObservationIndex]; + String gainVectorId = missingObservationId + ":" + Math.round(obsTimeOffsets[i] * 24.0 * 3600.0); //conversion to seconds + Results.putProgression("processing obs " + gainVectorId + "\n"); + // add to analysis increment for this obs + if (this.gainVectors.containsKey(gainVectorId)) { + IVector gainVector = this.gainVectors.get(gainVectorId); + double calculatedInnovation = dMissing.getValue(i, 0); + System.out.println("Calculated innovation for " + missingObservationId + ": " + calculatedInnovation + " index in dMissing: " + i); + delta.axpy(calculatedInnovation, gainVector); + } + } + System.out.println(); + } + } } diff --git a/core/java/resources/openDA/core_castor.jar b/core/java/resources/openDA/core_castor.jar index 41b7d6c87..d147c41f3 100644 Binary files a/core/java/resources/openDA/core_castor.jar and b/core/java/resources/openDA/core_castor.jar differ diff --git a/core/java/resources/openDA/core_castor_schemas.zip b/core/java/resources/openDA/core_castor_schemas.zip index fb36b4c4d..23158657d 100644 Binary files a/core/java/resources/openDA/core_castor_schemas.zip and b/core/java/resources/openDA/core_castor_schemas.zip differ diff --git a/core/java/resources/openDA/core_castor_src.zip b/core/java/resources/openDA/core_castor_src.zip index 18911cc0a..ee2019f00 100644 Binary files a/core/java/resources/openDA/core_castor_src.zip and b/core/java/resources/openDA/core_castor_src.zip differ diff --git a/core/java/src/org/openda/blackbox/wrapper/BBStochModelInstance.java b/core/java/src/org/openda/blackbox/wrapper/BBStochModelInstance.java index af3c6849b..0bdbbfc3c 100644 --- a/core/java/src/org/openda/blackbox/wrapper/BBStochModelInstance.java +++ b/core/java/src/org/openda/blackbox/wrapper/BBStochModelInstance.java @@ -1290,9 +1290,6 @@ public IVector[] getObservedLocalization(IObservationDescriptions observationDes ITreeVector modelWeightsTreeVector = getLocalizedCohnWeights(obsId[i], distance, xObs.getValue(i), yObs.getValue(i), zObs.getValue(i)); - // TreeVector to vector -// IVector statesWeightsArray = new Vector(modelWeightsTreeVector.getValues()); -// obsVectorArray[i] = statesWeightsArray; TreeVector weightsForFullState = new TreeVector("State-Weight"); weightsForFullState.addChild(noiseModelWeightsTreeVector); weightsForFullState.addChild(modelWeightsTreeVector); @@ -1301,15 +1298,27 @@ public IVector[] getObservedLocalization(IObservationDescriptions observationDes return obsVectorArray; } - //private IVector getLocalizedCohnWeights(String obsId, double distanceCohnMeters, double xObs, double yObs, double zObs){ + public IVector[] getRhoForLocalization(IObservationDescriptions observationDescriptions, double distance) { + // Distance is delegated to the ExchangeItem and then apply the Gaspari-Cohn function. + IVector xObs = observationDescriptions.getValueProperties("xposition"); + IVector yObs = observationDescriptions.getValueProperties("yposition"); + IVector zObs = observationDescriptions.getValueProperties("height"); + int obsCount = observationDescriptions.getObservationCount(); + + List exchangeItems = observationDescriptions.getExchangeItems(); + IGeometryInfo[] geometryInfos = new IGeometryInfo[obsCount]; + for (int i = 0; i < obsCount; i++) { + geometryInfos[i] = exchangeItems.get(i).getGeometryInfo(); + } + return GeometryUtils.calculateCohnWeightsBetweenPointGeometries(distance, geometryInfos, xObs, yObs, zObs); + } + private ITreeVector getLocalizedCohnWeights(String obsId, double distanceCohnMeters, double xObs, double yObs, double zObs){ - //IVector iWeightsVector = new IVector; TreeVector treeVector = new TreeVector("weights-for " + obsId); Collection vectorCollection = this.bbStochModelVectorsConfig.getStateConfig().getVectorCollection(); - //int k=0; for (BBStochModelVectorConfig vectorConfig : vectorCollection) { String id = vectorConfig.getId(); IExchangeItem prevExchangeItem = getExchangeItem(id); @@ -1322,14 +1331,11 @@ private ITreeVector getLocalizedCohnWeights(String obsId, double distanceCohnMet for (int xy = 0; xy < distancesForExchangeItem.length; xy++) { weightsForExchangeItem[xy] = GeometryUtils.calculateCohnWeight(distancesForExchangeItem[xy], distanceCohnMeters); - //k++; } treeVector.addChild(prevExchangeItem.getId(), weightsForExchangeItem); - //iWeightsVector[k] = weightsForExchangeItem; } } return treeVector; - //return iWeightsVector; } diff --git a/core/java/src/org/openda/exchange/ArrayGeometryInfo.java b/core/java/src/org/openda/exchange/ArrayGeometryInfo.java index 0eb4eda4e..141a40f9a 100644 --- a/core/java/src/org/openda/exchange/ArrayGeometryInfo.java +++ b/core/java/src/org/openda/exchange/ArrayGeometryInfo.java @@ -52,7 +52,7 @@ public class ArrayGeometryInfo implements IArrayGeometryInfo { boolean is3D; // Added to check if distance needs to be converted to meters // Before it always assumed it was WGS84 and calculated the "distances" in lat lon which gives unrealistic results since it depends on the exact location on earth - private boolean isWGS84; + boolean isWGS84; /** * latitudeValueIndices contains for each dimension in the latitudeCoordinateValues array diff --git a/core/java/src/org/openda/exchange/PointGeometryInfo.java b/core/java/src/org/openda/exchange/PointGeometryInfo.java index 790cd127d..8fb06a688 100644 --- a/core/java/src/org/openda/exchange/PointGeometryInfo.java +++ b/core/java/src/org/openda/exchange/PointGeometryInfo.java @@ -54,6 +54,7 @@ public PointGeometryInfo(double longitude, double latitude, double height){ this.setInfo(latitudeArray, latitudeValueIndices, latitudeQuantityInfo, longitudeArray, longitudeValueIndices, longitudeQuantityInfo, heightArray, heightValueIndices, heightQuantityInfo); + isWGS84 = true; } public double getLatitude(){ diff --git a/core/java/src/org/openda/interfaces/IModelInstance.java b/core/java/src/org/openda/interfaces/IModelInstance.java index 5d6305029..42820661a 100644 --- a/core/java/src/org/openda/interfaces/IModelInstance.java +++ b/core/java/src/org/openda/interfaces/IModelInstance.java @@ -20,7 +20,10 @@ package org.openda.interfaces; +import org.openda.utils.Vector; + import java.io.File; +import java.util.Arrays; /** * Model Instance interface. @@ -100,6 +103,20 @@ public interface IModelInstance extends IDataObject, IInstance { */ IVector[] getObservedLocalization(IObservationDescriptions observationDescriptions, double distance, int iDomain); + default IVector[] getRhoForLocalization(IObservationDescriptions observationDescriptions, double distance) { + return getDefaultRhoForLocalization(observationDescriptions); + } + + static IVector[] getDefaultRhoForLocalization(IObservationDescriptions observationDescriptions) { + int obsCount = observationDescriptions.getObservationCount(); + IVector[] rhoForLocalization = new IVector[obsCount]; + for (int i = 0; i < obsCount; i++) { + double[] weightsObservation = new double[obsCount]; + Arrays.fill(weightsObservation, 1); + rhoForLocalization[i] = new Vector(weightsObservation); + } + return rhoForLocalization; + } /************************************** *** Save/restore full internal state **************************************/ diff --git a/core/java/src/org/openda/utils/SortUtils.java b/core/java/src/org/openda/utils/SortUtils.java index 5dcb2fe53..736b9ec78 100644 --- a/core/java/src/org/openda/utils/SortUtils.java +++ b/core/java/src/org/openda/utils/SortUtils.java @@ -20,6 +20,7 @@ package org.openda.utils; import java.util.ArrayList; import java.util.Arrays; +import java.util.Objects; import org.openda.interfaces.IVector; @@ -30,6 +31,32 @@ */ public class SortUtils { + public static int indexOfString(String[] array, String string) { + return indexOfString(array, 0, array.length, string); + } + + private static int indexOfString(String[] array, int pos, int length, String string) { + if (string == null) { + for (int i = pos, n = pos + length; i < n; i++) { + if (array[i] == null) return i; + } + } else { + for (int i = pos, n = pos + length; i < n; i++) { + if (stringEquals(array[i], string)) return i; + } + } + return -1; + } + + private static boolean stringEquals(String s1, String s2) { + if (s1 == null && s2 == null) return true; + if (s1 == null || s2 == null) return false; + int length = s1.length(); + if (length != s2.length()) return false; + if (s1.hashCode() != s2.hashCode()) return false; + return s1.regionMatches(0, s2, 0, length); // 3 times faster than equals, no cast required + } + //merge type public enum MergeType {left, union, intersection}; diff --git a/core/java/src/org/openda/utils/geometry/GeometryUtils.java b/core/java/src/org/openda/utils/geometry/GeometryUtils.java index 70ba94429..e5837cec0 100644 --- a/core/java/src/org/openda/utils/geometry/GeometryUtils.java +++ b/core/java/src/org/openda/utils/geometry/GeometryUtils.java @@ -557,4 +557,21 @@ public static double calculateCohnWeight(double r, double cohnDistance) { return p1_12*z5 - p1_2*z4 + p5_8*z3 + p5_3*z2 - 5*z1 + 4 - p2_3/z1; } } + + public static IVector[] calculateCohnWeightsBetweenPointGeometries(double distance, IGeometryInfo[] geometryInfos, IVector xObs, IVector yObs, IVector zObs) { + int obsCount = xObs.getSize(); + IVector[] obsVectorArray = new IVector[xObs.getSize()]; + for (int i = 0; i < geometryInfos.length; i++) { + double[] weightsObservation = new double[obsCount]; + IGeometryInfo geometryInfo = geometryInfos[i]; + for (int j = 0; j < obsCount; j++) { + IArray distanceToPoint = geometryInfo.distanceToPoint(xObs.getValue(j), yObs.getValue(j), zObs.getValue(j)); + if (distanceToPoint.length() != 1) throw new RuntimeException("Calculating Rho for localization only supported for single point geometries."); + double obsDistance = distanceToPoint.getValueAsDouble(0); + weightsObservation[j] = calculateCohnWeight(obsDistance, distance); + } + obsVectorArray[i] = new Vector(weightsObservation); + } + return obsVectorArray; + } } diff --git a/core/java/src/org/openda/utils/io/KalmanGainStorage.java b/core/java/src/org/openda/utils/io/KalmanGainStorage.java index 8f64cb78a..922d23718 100644 --- a/core/java/src/org/openda/utils/io/KalmanGainStorage.java +++ b/core/java/src/org/openda/utils/io/KalmanGainStorage.java @@ -57,6 +57,7 @@ public class KalmanGainStorage { public static final Attribute STATION_IDENTIFICATION_CODE_ATT = new Attribute("long_name", "station identification code"); public static final Attribute TIME_SERIES_ID_ATT = new Attribute("cf_role", "timeseries_id"); public static final Attribute OBSERVATION_OFFSET_LONG_NAME_ATT = new Attribute("long_name", "offset in days for observations"); + public static final Attribute HK_LONG_NAME_ATT = new Attribute("long_name", "HK"); public static final Attribute KALMAN_GAIN_LONG_NAME_ATT = new Attribute("long_name", "kalman gain"); public static final Attribute FRACTIONS_UNIT_ATT = new Attribute("units", "1"); public static int DefaultMaxKeepVectorInXMLSize = 40; @@ -71,6 +72,7 @@ public class KalmanGainStorage { private String kalmanGainStorageFileName = "kalmanGainStorage.xml"; private int maxKeepVectorInXMLSize = DefaultMaxKeepVectorInXMLSize; private boolean useTimeStampInDirectoryName = true; + private double[][] hk = null; public enum StorageType { xml, netcdf, automatic, netcdf_cf @@ -200,25 +202,29 @@ public void writeKalmanGain(HashMap gainVectors, IVector[] kalmanGainColumns = new IVector[n]; Set gainIds = gainVectors.keySet(); int index=0; + //System.out.println("Observation order for writing kalman gain"); for(String gainId : gainIds){ observationIds[index] = obsId.get(gainId); observationOffsetsInDays[index] = obsTimeOffset.get(gainId); + //System.out.printf("ObsId:offset %s:%f%n", observationIds[index], observationOffsetsInDays[index]); kalmanGainColumns[index] = gainVectors.get(gainId); index++; } - writeKalmanGain(observationIds,observationOffsetsInDays,kalmanGainColumns); + writeKalmanGain(observationIds,observationOffsetsInDays,kalmanGainColumns, hk); } /** * Write the kalman gain matrix to xml file and, in case of large column vectors, to netdcf files. - * @param observationIds The observation identifiers + * + * @param observationIds The observation identifiers * @param observationOffsetsInDays The time offset of the observations, expressed in days, * relative to the time stamp for this kalman gain * (0 means: same time stamp as the gain, - * negative means before the kalman gain time stamp) - * @param kalmanGainColumns The vectors to be written to the kalman gain column files. + * negative means before the kalman gain time stamp) + * @param kalmanGainColumns The vectors to be written to the kalman gain column files. + * @param hk */ - public void writeKalmanGain(String[] observationIds, double[] observationOffsetsInDays, IVector[] kalmanGainColumns) { + public void writeKalmanGain(String[] observationIds, double[] observationOffsetsInDays, IVector[] kalmanGainColumns, double[][] hk) { if (kalmanGainColumns == null || kalmanGainColumns.length == 0 ) { throw new IllegalArgumentException(this.getClass().getName() + ": at least one kalman gain column must be provided"); @@ -236,10 +242,14 @@ public void writeKalmanGain(String[] observationIds, double[] observationOffsets File directoryForStorage = determineStorageDirectory(false); if (this.gainFileType == StorageType.netcdf_cf) { - writeKalmanGainToNetcdfCF(observationIds, observationOffsetsInDays, kalmanGainColumns, directoryForStorage); + writeKalmanGainToNetcdfCF(observationIds, observationOffsetsInDays, kalmanGainColumns, directoryForStorage, hk); return; } + writeKalmanGainToXML(observationIds, observationOffsetsInDays, kalmanGainColumns, directoryForStorage); + } + + private void writeKalmanGainToXML(String[] observationIds, double[] observationOffsetsInDays, IVector[] kalmanGainColumns, File directoryForStorage) { long start = System.currentTimeMillis(); // store general info KalmanGainStorageXML kgStorageXML = new KalmanGainStorageXML(); @@ -249,6 +259,14 @@ public void writeKalmanGain(String[] observationIds, double[] observationOffsets KalmanGainObservationsXML observationsXML = new KalmanGainObservationsXML(); kgStorageXML.setObservations(observationsXML); + fillObservationsXML(observationIds, observationOffsetsInDays, kalmanGainColumns, directoryForStorage, observationsXML, kgStorageXML); + File kgStorageXmlFile = new File(directoryForStorage, kalmanGainStorageFileName); + CastorUtils.write(kgStorageXML, kgStorageXmlFile, "opendaKalmanGainStorage", null, null); + long timePassed = System.currentTimeMillis() - start; + Results.putMessage("Writing kalman gain to xml took: " + timePassed); + } + + private void fillObservationsXML(String[] observationIds, double[] observationOffsetsInDays, IVector[] kalmanGainColumns, File directoryForStorage, KalmanGainObservationsXML observationsXML, KalmanGainStorageXML kgStorageXML) { for (int i = 0; i < observationIds.length; i++) { KalmanGainObservationXML observationXML = new KalmanGainObservationXML(); observationsXML.addObservation(observationXML); @@ -270,24 +288,20 @@ public void writeKalmanGain(String[] observationIds, double[] observationOffsets String vectorString = TreeVectorWriter.createVectorXML(kalmanGainColumns[i]); observationXMLChoice.setVector(vectorString); } - } else { - // netCdf-file - String netcdfFileName = columnFilePrefix + String.valueOf(i+1) + ".nc"; - observationXMLChoice.setFileName(netcdfFileName); - if (!kgStorageXML.hasStateSize()) { - kgStorageXML.setStateSize(kalmanGainColumns[i].getSize()); - } - netcdfFileName = new File(directoryForStorage, netcdfFileName).getAbsolutePath(); - writeKalmanGainColumnToNetCdfFile(netcdfFileName, kalmanGainColumns[i]); + continue; + } + // netCdf-file + String netcdfFileName = columnFilePrefix + String.valueOf(i + 1) + ".nc"; + observationXMLChoice.setFileName(netcdfFileName); + if (!kgStorageXML.hasStateSize()) { + kgStorageXML.setStateSize(kalmanGainColumns[i].getSize()); } + netcdfFileName = new File(directoryForStorage, netcdfFileName).getAbsolutePath(); + writeKalmanGainColumnToNetCdfFile(netcdfFileName, kalmanGainColumns[i]); } - File kgStorageXmlFile = new File(directoryForStorage, kalmanGainStorageFileName); - CastorUtils.write(kgStorageXML, kgStorageXmlFile, "opendaKalmanGainStorage", null, null); - long timePassed = System.currentTimeMillis() - start; - Results.putMessage("Writing kalman gain to xml took: " + timePassed); } - private void writeKalmanGainToNetcdfCF(String[] observationIds, double[] observationOffsetsInDays, IVector[] kalmanGainColumns, File directoryForStorage) { + private void writeKalmanGainToNetcdfCF(String[] observationIds, double[] observationOffsetsInDays, IVector[] kalmanGainColumns, File directoryForStorage, double[][] hk) { NetcdfFileWriter netcdfFileWriter = null; try { File file = new File(directoryForStorage, kalmanGainStorageFileName); @@ -297,6 +311,7 @@ private void writeKalmanGainToNetcdfCF(String[] observationIds, double[] observa Dimension stationDimension = netcdfFileWriter.addDimension(null, STATION_DIMENSION, observationIds.length); netcdfFileWriter.addDimension(null, CHAR_LENGTH_ID, STATION_ID_CHAR_LENGTH); Variable observationOffsetVariable = createObservationOffsetVariable(netcdfFileWriter); + Variable hkVariable = hk != null ? createHKVariable(netcdfFileWriter) : null; Variable stationVariable = createStationVariable(netcdfFileWriter); KalmanGainVariableData[] kalmanGainVariableData = getKalmanGainVariableData(kalmanGainColumns, netcdfFileWriter, stationDimension); @@ -305,6 +320,7 @@ private void writeKalmanGainToNetcdfCF(String[] observationIds, double[] observa writeStationData(observationIds, netcdfFileWriter, stationVariable); writeTimeStampData(netcdfFileWriter, timeStampVariable); writeObservationOffsetData(observationOffsetsInDays, netcdfFileWriter, observationOffsetVariable); + if (hkVariable != null) writeHK(hk, netcdfFileWriter, hkVariable); for (int observationIndex = 0; observationIndex < kalmanGainColumns.length; observationIndex++) { @@ -367,6 +383,16 @@ private void writeObservationOffsetData(double[] observationOffsetsInDays, Netcd netcdfFileWriter.write(observationOffsetVariable, observationOffsetArray); } + private void writeHK(double[][] hk, NetcdfFileWriter netcdfFileWriter, Variable hkVariable) throws IOException, InvalidRangeException { + ArrayDouble.D2 hkArray = new ArrayDouble.D2(hk.length, hk.length); + for (int i = 0; i < hk.length; i++) { + for (int j = 0; j < hk.length; j++) { + hkArray.set(i, j, hk[i][j]); + } + } + netcdfFileWriter.write(hkVariable, hkArray); + } + private void writeStationData(String[] observationIds, NetcdfFileWriter netcdfFileWriter, Variable stationVariable) throws IOException, InvalidRangeException { ArrayChar.D2 stationArray = new ArrayChar.D2(observationIds.length, STATION_ID_CHAR_LENGTH); for (int i = 0; i < observationIds.length; i++) { @@ -463,6 +489,13 @@ private Variable createObservationOffsetVariable(NetcdfFileWriter netcdfFileWrit return observationOffsetVariable; } + private Variable createHKVariable(NetcdfFileWriter netcdfFileWriter) { + Variable observationOffsetVariable = netcdfFileWriter.addVariable(null, "HK", DataType.DOUBLE, STATION_DIMENSION + ' ' + STATION_DIMENSION); + observationOffsetVariable.addAttribute(HK_LONG_NAME_ATT); + observationOffsetVariable.addAttribute(FRACTIONS_UNIT_ATT); + return observationOffsetVariable; + } + private Variable createTimeStampVariable(NetcdfFileWriter netcdfFileWriter) { Dimension timeStampDimension = netcdfFileWriter.addDimension(null, "time_stamp_dimension", 1); ArrayList dimensionList = new ArrayList<>(); @@ -526,6 +559,14 @@ public IVector[] getKalmanGainColumns() { return this.kalmanGainColumns; } + public double[][] getHk() { + return hk; + } + + public void setHk(double[][] hk) { + this.hk = hk; + } + /** * Read the kalman gain matrix from xml file and, in case of large column vectors, to netdcf files. */ @@ -623,6 +664,10 @@ private void readVariables(List variables, int stationLength, TreeVect this.timeStampAsMJD = ((double[]) variable.read().get1DJavaArray(double.class))[0]; continue; } + if ("HK".equals(shortName)) { + this.hk = (double[][]) variable.read().copyToNDJavaArray(); + continue; + } int[] shape = variable.getShape(); int[] shapeForRead = variable.getShape(); shapeForRead[0] = 1; diff --git a/core/java/test/org/openda/utils/geometry/GeometryUtilsTest.java b/core/java/test/org/openda/utils/geometry/GeometryUtilsTest.java index a8310ea47..25549d0c8 100644 --- a/core/java/test/org/openda/utils/geometry/GeometryUtilsTest.java +++ b/core/java/test/org/openda/utils/geometry/GeometryUtilsTest.java @@ -23,6 +23,7 @@ import junit.framework.TestCase; import org.openda.exchange.ArrayGeometryInfo; import org.openda.exchange.IrregularGridGeometryInfo; +import org.openda.exchange.PointGeometryInfo; import org.openda.interfaces.IArray; import org.openda.interfaces.IGeometryInfo; import org.openda.interfaces.IVector; @@ -256,4 +257,32 @@ public void testToCurvilinearGeometryInfo() { //if geometryInfo is already curvilinear. assertEquals(geometryInfo, geometryInfo.toCurvilinearGeometryInfo()); } + + public void testCalculateCohnWeightsForPointGeometries() { + // Mock inputs + double distance = 300000.0; + IGeometryInfo[] geometryInfos = new IGeometryInfo[3]; + geometryInfos[0] = new PointGeometryInfo(1.0, 4.0, 0.0); + geometryInfos[1] = new PointGeometryInfo(2.0, 5.0, 0.0); + geometryInfos[2] = new PointGeometryInfo(3.0, 6.0, 0.0); + + int obsCount = 3; + IVector xObs = new Vector(new double[]{1.0, 2.0, 3.0}); + IVector yObs = new Vector(new double[]{4.0, 5.0, 6.0}); + IVector zObs = new Vector(new double[]{0.0, 0.0, 0.0}); + + IVector[] weights = GeometryUtils.calculateCohnWeightsBetweenPointGeometries(distance, geometryInfos, xObs, yObs, zObs); + + assertNotNull(weights[0]); + assertEquals(obsCount, weights[0].getSize()); + assertEquals(1.0, weights[0].getValue(0), 1e-6); + assertEquals(0.661056, weights[0].getValue(1), 1e-6); + assertEquals(0.177380, weights[0].getValue(2), 1e-6); + assertEquals(0.661056, weights[1].getValue(0), 1e-6); + assertEquals(1.0, weights[1].getValue(1), 1e-6); + assertEquals(0.661469, weights[1].getValue(2), 1e-6); + assertEquals(0.177380, weights[2].getValue(0), 1e-6); + assertEquals(0.661469, weights[2].getValue(1), 1e-6); + assertEquals(1.0, weights[2].getValue(2), 1e-6); + } } diff --git a/core/java/test/org/openda/utils/io/KalmanGainStorageTest.java b/core/java/test/org/openda/utils/io/KalmanGainStorageTest.java index e4255f5d3..190f6f790 100644 --- a/core/java/test/org/openda/utils/io/KalmanGainStorageTest.java +++ b/core/java/test/org/openda/utils/io/KalmanGainStorageTest.java @@ -22,6 +22,7 @@ import org.openda.blackbox.config.BBUtils; import org.openda.costa.CtaTreeVector; import org.openda.costa.CtaVector; +import org.openda.exchange.timeseries.TimeUtils; import org.openda.interfaces.ITreeVector; import org.openda.interfaces.IVector; import org.openda.utils.OpenDaTestSupport; @@ -30,6 +31,7 @@ import java.io.File; import java.io.IOException; +import java.text.ParseException; import java.util.ArrayList; /** @@ -59,7 +61,7 @@ public void testReadWriteKalmanGainXML() throws Exception { KalmanGainStorage kgStorageOut = new KalmanGainStorage(testRunDataDir, timeAsMJD); kgStorageOut.setKalmanGainStorageFileName("kalmanGainStorage_out.xml"); kgStorageOut.writeKalmanGain(kgStorageIn.getObservationIds(), - kgStorageIn.getObservationOffsetInDays(), kgStorageIn.getKalmanGainColumns()); + kgStorageIn.getObservationOffsetInDays(), kgStorageIn.getKalmanGainColumns(), null); testData.FilesAreIdentical( new File(kgStorageDir, "kalmanGainStorage_out.xml"), @@ -85,12 +87,67 @@ public void testReadWriteKalmanGainNetcdfCF() { KalmanGainStorage kgStorageOut = new KalmanGainStorage(testRunDataDir, timeAsMJD); kgStorageOut.setColumnFileType(KalmanGainStorage.StorageType.netcdf_cf); kgStorageOut.setKalmanGainStorageFileName("KalmanGainStorage.nc"); - kgStorageOut.writeKalmanGain(observationIds, kgStorageIn.getObservationOffsetInDays(), kalmanGainColumnsIn); + kgStorageOut.writeKalmanGain(observationIds, kgStorageIn.getObservationOffsetInDays(), kalmanGainColumnsIn, kgStorageIn.getHk()); kgStorageOut.readKalmanGain(); checkKalmanGainContents(timeAsMJD, kgStorageOut, kgStorageOut.getKalmanGainColumns(), kgStorageOut.getObservationIds()); } + public void testReadWriteKalmanGainNetcdfCFHK() throws ParseException { + + double timeAsMJD = TimeUtils.date2Mjd("202209011200"); + + KalmanGainStorage kgStorageIn = new KalmanGainStorage(testRunDataDir, timeAsMJD); + kgStorageIn.setColumnFileType(KalmanGainStorage.StorageType.netcdf_cf); + kgStorageIn.setKalmanGainStorageFileName("kalmanGainStorageHK.nc"); + kgStorageIn.readKalmanGain(); + + IVector[] kalmanGainColumnsIn = kgStorageIn.getKalmanGainColumns(); + assertEquals(4, kalmanGainColumnsIn.length); + + String[] observationIds = kgStorageIn.getObservationIds(); + String[] expectedObsIds = {"WESTKPLE.waterlevel", "VLISSGN.waterlevel", "SCHEVNGN.waterlevel", "ROOMPBTN.waterlevel"}; + assertEquals(4, observationIds.length); + for (int i = 0; i < observationIds.length; i++) { + assertEquals(expectedObsIds[i], observationIds[i]); + } + double[][] hk = kgStorageIn.getHk(); + assertNotNull(hk); + assertEquals(4, hk.length); + assertEquals(4, hk[0].length); + double[][] expectedHK = {{0.19371548763628316, 0.21185091852367516, 0.23657124885992412, 0.21748081376956335} , {0.21185091852367666, 0.23168416850382453, 0.2587187889765871, 0.23784112835038507},{0.2365712488599269, 0.2587187889765884, 0.28890800869895084, 0.26559418838590676}, {0.2174808137695654, 0.23784112835038568, 0.2655941883859062, 0.24416170815767535}}; + for (int i = 0; i < hk.length; i++) { + double[] hkRow = hk[i]; + for (int j = 0; j < hkRow.length; j++) { + assertEquals(expectedHK[i][j], hkRow[j]); + } + } + } + + public void testReadWriteKalmanGainNetcdfCFDoubleState() { + + double timeAsMJD = 60310; + + KalmanGainStorage kgStorageIn = new KalmanGainStorage(testRunDataDir, timeAsMJD); + kgStorageIn.setColumnFileType(KalmanGainStorage.StorageType.netcdf_cf); + kgStorageIn.setKalmanGainStorageFileName("KalmanGainStorage.nc"); + kgStorageIn.readKalmanGain(); + + IVector[] kalmanGainColumnsIn = kgStorageIn.getKalmanGainColumns(); + assertEquals(8, kalmanGainColumnsIn.length); + + String[] observationIds = kgStorageIn.getObservationIds(); + checkKalmanGainContentsDoubleState(timeAsMJD, kgStorageIn, kalmanGainColumnsIn, observationIds); + + KalmanGainStorage kgStorageOut = new KalmanGainStorage(testRunDataDir, timeAsMJD); + kgStorageOut.setColumnFileType(KalmanGainStorage.StorageType.netcdf_cf); + kgStorageOut.setKalmanGainStorageFileName("KalmanGainStorage.nc"); + kgStorageOut.writeKalmanGain(observationIds, kgStorageIn.getObservationOffsetInDays(), kalmanGainColumnsIn, null); + kgStorageOut.readKalmanGain(); + + checkKalmanGainContentsDoubleState(timeAsMJD, kgStorageOut, kgStorageOut.getKalmanGainColumns(), kgStorageOut.getObservationIds()); + } + private void checkKalmanGainContents(double timeAsMJD, KalmanGainStorage kgStorageIn, IVector[] kalmanGainColumns, String[] observationIds) { String[] expectedObservationIds = {"WICK.waterlevel", "VLISSGN.waterlevel", "DENHDR.waterlevel", "NORTHSS.waterlevel", "DOVR.waterlevel", "SHEERNS.waterlevel", "HOEKVHLD.waterlevel", "LOWST.waterlevel"}; for (int i = 0; i < expectedObservationIds.length; i++) { @@ -136,30 +193,16 @@ private void checkKalmanGainContents(double timeAsMJD, KalmanGainStorage kgStora assertEquals(-0.052612, values[434], delta); assertEquals(0.015112, values[435], delta); assertEquals(-0.038791, values[436], delta); - } - - public void testReadWriteKalmanGainNetcdfCFDoubleState() { - - double timeAsMJD = 60310; - - KalmanGainStorage kgStorageIn = new KalmanGainStorage(testRunDataDir, timeAsMJD); - kgStorageIn.setColumnFileType(KalmanGainStorage.StorageType.netcdf_cf); - kgStorageIn.setKalmanGainStorageFileName("KalmanGainStorage.nc"); - kgStorageIn.readKalmanGain(); - - IVector[] kalmanGainColumnsIn = kgStorageIn.getKalmanGainColumns(); - assertEquals(8, kalmanGainColumnsIn.length); - String[] observationIds = kgStorageIn.getObservationIds(); - checkKalmanGainContentsDoubleState(timeAsMJD, kgStorageIn, kalmanGainColumnsIn, observationIds); + double[][] hk = kgStorageIn.getHk(); - KalmanGainStorage kgStorageOut = new KalmanGainStorage(testRunDataDir, timeAsMJD); - kgStorageOut.setColumnFileType(KalmanGainStorage.StorageType.netcdf_cf); - kgStorageOut.setKalmanGainStorageFileName("KalmanGainStorage.nc"); - kgStorageOut.writeKalmanGain(observationIds, kgStorageIn.getObservationOffsetInDays(), kalmanGainColumnsIn); - kgStorageOut.readKalmanGain(); - - checkKalmanGainContentsDoubleState(timeAsMJD, kgStorageOut, kgStorageOut.getKalmanGainColumns(), kgStorageOut.getObservationIds()); + assertEquals(8, hk.length); + assertEquals(8, hk[0].length); + for (int i = 0; i < hk.length; i++) { + for (int j = 0; j < hk[0].length; j++) { + assertEquals(i * 10 + j, hk[i][j], 0.000001); + } + } } private void checkKalmanGainContentsDoubleState(double timeAsMJD, KalmanGainStorage kgStorageIn, IVector[] kalmanGainColumns, String[] observationIds) { @@ -244,7 +287,7 @@ private void writeAndReadKalmanGain(double timeAsMJD, int columnSize) { KalmanGainStorage kgStorageOut = new KalmanGainStorage(testRunDataDir, timeAsMJD); kgStorageOut.setComment(commentOut); kgStorageOut.setKalmanGainStorageFileName("fourColumnStorage.xml"); - kgStorageOut.writeKalmanGain(observationIdsOut, observationOffsetsInDaysOut, kalmanGainColumnsOut); + kgStorageOut.writeKalmanGain(observationIdsOut, observationOffsetsInDaysOut, kalmanGainColumnsOut, null); // read the kalman gain KalmanGainStorage kgStorageIn = new KalmanGainStorage(testRunDataDir, timeAsMJD); diff --git a/core/java/test/org/openda/utils/io/testData/kgStorage_200701010000/KalmanGainStorage.nc b/core/java/test/org/openda/utils/io/testData/kgStorage_200701010000/KalmanGainStorage.nc index 4c462075a..d94d1c9b0 100644 Binary files a/core/java/test/org/openda/utils/io/testData/kgStorage_200701010000/KalmanGainStorage.nc and b/core/java/test/org/openda/utils/io/testData/kgStorage_200701010000/KalmanGainStorage.nc differ diff --git a/core/java/test/org/openda/utils/io/testData/kgStorage_202209011200/kalmanGainStorageHK.nc b/core/java/test/org/openda/utils/io/testData/kgStorage_202209011200/kalmanGainStorageHK.nc new file mode 100644 index 000000000..f2660ab20 Binary files /dev/null and b/core/java/test/org/openda/utils/io/testData/kgStorage_202209011200/kalmanGainStorageHK.nc differ diff --git a/core/xmlSchemas/algorithm/steadyStateFilter.xsd b/core/xmlSchemas/algorithm/steadyStateFilter.xsd index 58079e30c..b085065a6 100644 --- a/core/xmlSchemas/algorithm/steadyStateFilter.xsd +++ b/core/xmlSchemas/algorithm/steadyStateFilter.xsd @@ -77,6 +77,11 @@ Skip gain when the model results differ from the observations more than this factor times the standard deviation of the observations + + + When this option is true AND in the kalman gain storage the term HK has been stored, for locations with missing observations the innovation value + + diff --git a/model_metaswap/java/src/org/openda/model_metaswap/SwapResultFile.java b/model_metaswap/java/src/org/openda/model_metaswap/SwapResultFile.java index 09568d999..7f0ac8b4b 100644 --- a/model_metaswap/java/src/org/openda/model_metaswap/SwapResultFile.java +++ b/model_metaswap/java/src/org/openda/model_metaswap/SwapResultFile.java @@ -3,6 +3,7 @@ import org.openda.exchange.timeseries.TimeUtils; import org.openda.interfaces.IDataObject; import org.openda.interfaces.IExchangeItem; +import org.openda.utils.SortUtils; import org.openda.utils.io.CsvReader; import java.io.File; @@ -47,10 +48,10 @@ public void initialize(File workingDir, String[] arguments) { try { CsvReader csvReader = new CsvReader(sourceFile); String[] headerSplit = csvReader.readCSVLineTrimElements(); - int s01Index = indexOfString(headerSplit, S01MM_COLUMN_NAME); - int dprztbIndex = indexOfString(headerSplit, DPRZTB_COLUMN_NAME); - int yearIndex = indexOfString(headerSplit, "yr"); - int dayIndex = indexOfString(headerSplit, "td(d)"); + int s01Index = SortUtils.indexOfString(headerSplit, S01MM_COLUMN_NAME); + int dprztbIndex = SortUtils.indexOfString(headerSplit, DPRZTB_COLUMN_NAME); + int yearIndex = SortUtils.indexOfString(headerSplit, "yr"); + int dayIndex = SortUtils.indexOfString(headerSplit, "td(d)"); List doubleValuesList = new ArrayList<>(); List doubleTimesList = new ArrayList<>(); String[] lineElements = csvReader.readCSVLineTrimElements(); @@ -77,32 +78,6 @@ public void initialize(File workingDir, String[] arguments) { } } - private static int indexOfString(String[] array, String string) { - return indexOfString(array, 0, array.length, string); - } - - private static int indexOfString(String[] array, int pos, int length, String string) { - if (string == null) { - for (int i = pos, n = pos + length; i < n; i++) { - if (array[i] == null) return i; - } - } else { - for (int i = pos, n = pos + length; i < n; i++) { - if (stringEquals(array[i], string)) return i; - } - } - return -1; - } - - private static boolean stringEquals(String s1, String s2) { - if (s1 == s2) return true; - if (s1 == null || s2 == null) return false; - int length = s1.length(); - if (length != s2.length()) return false; - if (s1.hashCode() != s2.hashCode()) return false; - return s1.regionMatches(0, s2, 0, length); // 3 times faster than equals, no cast required - } - private static double[] unboxDoubleArray(Double[] array) { if (array == null)