using System; using System.Collections.Generic; using System.Linq; using System.Text; using MikeSheInOpenDA.Spatial; using OpenDA.DotNet.OpenMI.Bridge; using OpenMI.Standard2; using OpenMI.Standard2.TimeSpace; using DHI.OpenMI2.Sdk.Spatial; using DHI.OpenMI2.Sdk.Backbone; namespace MikeSheInOpenDA { public class MikeSheOpenMITimespaceComponentExtensions : DHI.OpenMI2.MikeShe.WMEngineAccess,ITimeSpaceComponentExtensions { public ITime currentTime() { return this.CurrentTime; } /// /// For easy observation handling a model can provide some additional exhange items especially for model predictions at observation locations /// /// public IList getAdditionalExchangeItem() { IList additional = new List(); // We will hack something here! return additional; } public IList CreateModelIndicesHashTable(OpenDA.DotNet.Interfaces.IObservationDescriptions observationDescriptions) { IList modelIndices = new List(); String[] keys = observationDescriptions.PropertyKeys; String[] quantity = observationDescriptions.GetStringProperties("quantity"); double[] xpos = observationDescriptions.GetValueProperties("xposition").Values; double[] ypos = observationDescriptions.GetValueProperties("yposition").Values; double[] height = observationDescriptions.GetValueProperties("height").Values; int observationCount = observationDescriptions.ObservationCount; int nObs = quantity.Length; // same as observationCount? // The heights should be specified in an array of integers representing the layer. Check if the values are indeed integers or close to integers before converting // the array of doubles to an array of integers. const double tolerance = 1e-5; int[] layer = new int[observationCount]; for (int i = 0; i < observationCount; i++) { layer[i] = Convert.ToInt32(height[i]); if (Math.Abs(layer[i] - height[i]) > tolerance) { throw new Exception("The height specified in the observation was not an integer. Observation \n"); } } // An array of model values corresponding to the observation points. double[] Hx = new double[nObs]; List listExchangeIds = ObsIDtoExchangeId(observationDescriptions); for (int obsC = 0; obsC < nObs; obsC++) { IDictionary modelCoord = GetModelCoordinates(listExchangeIds[obsC]); IXYLayerPoint obsPoint = new XYLayerPoint(xpos[obsC], ypos[obsC], layer[obsC]); int modelVariableIndex = XYZGeometryTools.ModelIndexForPoint(obsPoint, modelCoord); modelIndices.Add(modelVariableIndex); if (modelVariableIndex < 0) { throw new Exception("The observation point was NOT in the model grid! For Point: (" + xpos[obsC].ToString() + "," + ypos[obsC].ToString() + "," + layer[obsC].ToString() + ") \n"); } } return modelIndices; } private List ObsIDtoExchangeId(OpenDA.DotNet.Interfaces.IObservationDescriptions observationDescriptions) { String[] keys = observationDescriptions.PropertyKeys; String[] quantity = observationDescriptions.GetStringProperties("quantity"); double[] xpos = observationDescriptions.GetValueProperties("xposition").Values; double[] ypos = observationDescriptions.GetValueProperties("yposition").Values; double[] height = observationDescriptions.GetValueProperties("height").Values; int observationCount = observationDescriptions.ObservationCount; int nObs = quantity.Length; // same as observationCount? List exchangeItemId = new List(); for (int obsC = 0; obsC < nObs; obsC++) { if (quantity[obsC].Equals("Head", StringComparison.OrdinalIgnoreCase)) { exchangeItemId.Add("head elevation in saturated zone,SZ3DGrid"); } else if (quantity[obsC].Equals("SoilMoisture", StringComparison.OrdinalIgnoreCase)) { exchangeItemId.Add("water content in unsaturated zone,WMUZ3DGrid"); } else if (quantity[obsC].Equals("SurfaceTemperature", StringComparison.OrdinalIgnoreCase)) { exchangeItemId.Add("Surface temperature (effective),BaseGrid"); } else { throw new Exception("Cannot (yet) handle obversvations of quantity (" + quantity[obsC] + ")"); } } return exchangeItemId; } /// /// Hx = Gets the model values at the indices that are passed in the list. /// /// The description of the observations. /// List of indices to be retrieved from the model. /// Array of model values at the given indices. public double[] ModelValuesAtProvidedIndices(OpenDA.DotNet.Interfaces.IObservationDescriptions observationDescriptions, IList indices ) { String[] keys = observationDescriptions.PropertyKeys; String[] quantity = observationDescriptions.GetStringProperties("quantity"); double[] xpos = observationDescriptions.GetValueProperties("xposition").Values; double[] ypos = observationDescriptions.GetValueProperties("yposition").Values; double[] height = observationDescriptions.GetValueProperties("height").Values; int observationCount = observationDescriptions.ObservationCount; int nObs = quantity.Length; // same as observationCount? // The heights should be specified in an array of integers representing the layer. Check if the values are indeed integers or close to integers before converting // the array of doubles to an array of integers. const double tolerance = 1e-5; int[] layer = new int[observationCount]; for (int i = 0; i < observationCount; i++) { layer[i] = Convert.ToInt32(height[i]); if (Math.Abs(layer[i] - height[i]) > tolerance) { throw new Exception("The height specified in the observation was not an integer. Observation \n"); } } List obsIds = ObsIDtoExchangeId(observationDescriptions); // An array of model values corresponding to the observation points. double[] Hx = GetModelValuesDifferentVariables(obsIds, indices); return Hx; } public double[] getObservedValues(OpenDA.DotNet.Interfaces.IObservationDescriptions observationDescriptions) { String[] keys = observationDescriptions.PropertyKeys; String[] quantity = observationDescriptions.GetStringProperties("quantity"); double[] xpos = observationDescriptions.GetValueProperties("xposition").Values; double[] ypos = observationDescriptions.GetValueProperties("yposition").Values; double[] height = observationDescriptions.GetValueProperties("height").Values; int observationCount = observationDescriptions.ObservationCount; int nObs = quantity.Length; // same as observationCount? // The heights should be specified in an array of integers representing the layer. Check if the values are indeed integers or close to integers before converting // the array of doubles to an array of integers. for (int i = 0; i < height.Length; i++) { if (double.IsNaN(height[i])) height[i] = 0.0; } const double tolerance = 1e-5; int[] layer = new int[observationCount]; for (int i = 0; i < observationCount; i++) { layer[i] = Convert.ToInt32(height[i]); if (Math.Abs(layer[i] - height[i]) > tolerance) { throw new Exception("The height specified in the observation was not an integer. Observation \n"); } } // An array of model values corresponding to the observation points. double[] Hx = new double[nObs]; for (int obsC = 0; obsC < nObs; obsC++) { // Set exchangeItem that corresponds to EntityID (no conversion yet) String exchangeItemId; if (quantity[obsC].Equals("Head", StringComparison.OrdinalIgnoreCase)) { exchangeItemId = "head elevation in saturated zone,SZ3DGrid"; } else if (quantity[obsC].Equals("SoilMoisture", StringComparison.OrdinalIgnoreCase)) { exchangeItemId = "water content in unsaturated zone,WMUZ3DGrid"; } else if (quantity[obsC].Equals("SurfaceTemperature", StringComparison.OrdinalIgnoreCase)) { exchangeItemId = "Surface temperature (effective),BaseGrid"; } else { throw new Exception("Cannot (yet) handle obversvations of quantity (" + quantity[obsC] + ")"); } IDictionary modelCoord = GetModelCoordinates(exchangeItemId); IXYLayerPoint obsPoint = new XYLayerPoint(xpos[obsC], ypos[obsC], layer[obsC]); int modelVariableIndex = XYZGeometryTools.ModelIndexForPoint(obsPoint, modelCoord); if(modelVariableIndex >= 0) { Hx[obsC] = GetModelValue(exchangeItemId, modelVariableIndex); } else { throw new Exception("The observation point was NOT in the model grid! For Point: (" + xpos[obsC].ToString() + "," + ypos[obsC].ToString()+ "," + layer[obsC].ToString() + ") \n" ); } } return Hx; } /// /// OpenMI does not know about localization. In this method, the user can implement localization for OpenMI models in this method /// /// the exchange ID string /// OpenDA type of observation description /// the distance for gaussian localization /// public double[][] getLocalization(string exchangeItemId, OpenDA.DotNet.Interfaces.IObservationDescriptions observationDescriptions, double locDistance) { //Get the Keys from the observer String[] keys = observationDescriptions.PropertyKeys; String[] quantity = observationDescriptions.GetStringProperties("quantity"); double[] xpos = observationDescriptions.GetValueProperties("xposition").Values; double[] ypos = observationDescriptions.GetValueProperties("yposition").Values; double[] height = observationDescriptions.GetValueProperties("height").Values; int observationCount = observationDescriptions.ObservationCount; // The heights should be specified in an array of integers representing the layer. Check if the values are indeed integers or close to integers before converting // the array of doubles to an array of integers. for (int i = 0; i < height.Length; i++) { if (double.IsNaN(height[i])) height[i] = 0.0;} const double tolerance = 1e-5; int[] layer = new int[observationCount]; for (int i = 0; i < observationCount; i++) { layer[i] = Convert.ToInt32(height[i]); if ( Math.Abs(layer[i] - height[i] ) > tolerance ) { throw new Exception("The height specified in the observation was not an integer. Observation \n"); } } // Gets the Grid type for the model. Can be a number of possibilities depending on the variable. GeometryTypes geometrytype = GetGridType(exchangeItemId); //BaseGrid if (geometrytype == GeometryTypes.Geometry2D) { return GetLocalized2D(exchangeItemId, observationCount, locDistance, xpos, ypos); } if (geometrytype == GeometryTypes.Geometry3D) { return GetLocalized3D(exchangeItemId, observationCount, locDistance, xpos, ypos, layer); } else { throw new NotImplementedException("Only 3D SZ and 2D BaseGrid supported so far."); } } #region PrivateMethods private double[] GetModelValuesDifferentVariables(List exchangeItemId, IList modelVariableIndices) { double[] modelValues = new double[modelVariableIndices.Count]; for (int i = 0; i < modelVariableIndices.Count; i++) { char[] delimiterChars = { ',' }; string[] words = exchangeItemId[i].Split(delimiterChars); string openMIVariableID = words[0].Trim(); ITimeSpaceOutput output = FindOutputItem(openMIVariableID); var modVars = output.Values.Values2D[0].Cast().ToArray(); modelValues[i] = modVars[modelVariableIndices[i]]; } return modelValues; } /// /// Get the values from the model at the given indices /// /// The exchange id of the variable to retrived /// List of indices in the model to retrive values from. /// Corresponding model values. private double[] GetModelValuesSameVariable(string exchangeItemId, IList modelVariableIndices) { char[] delimiterChars = { ',' }; string[] words = exchangeItemId.Split(delimiterChars); string openMIVariableID = words[0].Trim(); double[] modelValues = new double[modelVariableIndices.Count]; //Get values from model ITimeSpaceOutput output = FindOutputItem(openMIVariableID); var modVars = output.Values.Values2D[0].Cast().ToArray(); for (int i = 0; i < modelVariableIndices.Count; i++) { modelValues[i] = modVars[modelVariableIndices[i]]; } return modelValues; } private double GetModelValue(string exchangeItemId, int modelVariableIndex) { char[] delimiterChars = { ',' }; string[] words = exchangeItemId.Split(delimiterChars); string openMIVariableID = words[0].Trim(); ITimeSpaceOutput output = FindOutputItem(openMIVariableID); return output.Values.Values2D[0].Cast().ToArray()[modelVariableIndex]; } private ITimeSpaceOutput FindOutputItem(string outputItemId) { foreach (ITimeSpaceOutput outputItem in base.Outputs) { if (outputItem.Id.Equals(outputItemId)) { return outputItem; } } throw new Exception("Output item \"" + outputItemId + "not found."); } /// /// BASEGRID /// Returns a localization hash table. /// For each observation point, calculate the gaussian distance to all other points in the model grid. /// /// The exchange ID string /// The number of observation points /// The localization distance (for Gaussian distance) /// array of observation X coordinates /// array of observation Y coordinates /// array of observation Z heights (layer in the model) /// private double[][] GetLocalized3D(string exchangeItemId, int observationCount, double locDistance, double[] xpos, double[] ypos, int[] layer) { double[][] localized2D = new double[observationCount][]; var mshe = base.WMEngine; // SZGrid int n = mshe.Grid.ElementCount; IDictionary modelCoord = GetModelCoordinates(exchangeItemId); for (int obsC = 0; obsC < observationCount; obsC++) { localized2D[obsC] = new double[n]; IXYLayerPoint obsPoint = new XYLayerPoint(xpos[obsC], ypos[obsC], 0); if (XYZGeometryTools.IsPointInModelPlain(obsPoint, modelCoord)) { for (int i = 0; i < modelCoord.Count; i++) { if (Convert.ToInt32(obsPoint.Layer) == modelCoord[i].Layer) { double distanceC = XYZGeometryTools.CalculatePointToPointDistance2D(modelCoord[i].MidPoint, obsPoint); localized2D[obsC][i] = normalCooefs(distanceC, locDistance); } } } } return localized2D; } /// /// BASEGRID /// Returns a localization hash table. /// For each observation point, calculate the gaussian distance to all other points in the model grid. /// /// The exchange ID string /// The number of observation points /// The localization distance (for Gaussian distance) /// array of observation X coordinates /// array of observation Y coordinates /// private double[][] GetLocalized2D(string exchangeItemId, int observationCount, double locDistance, double[] xpos, double[] ypos) { double[][] localized2D = new double[observationCount][]; var mshe = base.WMEngine; // BASEGRID int n = mshe.Grid.ElementCount; IDictionary modelCoord = GetModelCoordinates(exchangeItemId); for (int obsC = 0; obsC < observationCount; obsC++) { localized2D[obsC] = new double[n]; IXYLayerPoint obsPoint = new XYLayerPoint(xpos[obsC], ypos[obsC], 0); if (XYZGeometryTools.IsPointInModelPlain(obsPoint, modelCoord)) { for (int i = 0; i < modelCoord.Count; i++) { if (Convert.ToInt32(obsPoint.Layer) == modelCoord[i].Layer) { double distanceC = XYZGeometryTools.CalculatePointToPointDistance2D(modelCoord[i].MidPoint, obsPoint); localized2D[obsC][i] = normalCooefs(distanceC, locDistance); } } } } return localized2D; } /// /// Returns the grid type based on the exchangeID. The exchange ID is a string delimiated by a comma. /// The second string after the comma contains the information of the Grid. /// Options include: /// 1) SZ3DGrid used for hydraulic head /// 2) BaseGrid used for land surface variabels (surface temperature, evapotranspiration, ...) /// 3) WMUZ3DGrid used for soil moisture content (in the UZ) /// /// string deliminated by a comma where the string after the comma contains the grid information. /// The geometry type of the variable. This is MikeSHE specific GeometryTypes GetGridType(string exchangeItemId) { IBaseOutput baseOut = base._outputExchangeItems.First(vID => string.Compare(vID.Id, exchangeItemId) == 0); char[] delimiterChars = { ',' }; string[] words = baseOut.Description.Split(delimiterChars); string gridTypewords = words[1].Trim(); if (string.Compare(gridTypewords, "SZ3DGrid", 0) == 0) { return GeometryTypes.Geometry3D; } else if (string.Compare(gridTypewords, "BaseGrid", 0) == 0) { return GeometryTypes.Geometry2D; } else if (string.Compare(gridTypewords, "WMUZ3DGrid", 0) == 0) { return GeometryTypes.Geometry3DUZ; } else { throw new Exception("Other types do exisit (UZ...)"); } } /// /// 2D BaseGrid !! /// Creates a dictionary with key equal to the model state index and the value the spatial information of that state index. /// /// The geometric type of the exchange itme (2d or 3d) /// The exchange item base output /// the string id of the exchange item. /// private IDictionary GetModelCoordinates2D(GeometryTypes gType, IBaseOutput baseOut, string elementID) { IDictionary modelEntities = new Dictionary(); int elementIDNumber; int n; try { elementIDNumber = WMEngine.GetElementCount(elementID); n = baseOut.ElementSet().ElementCount; } catch { Console.WriteLine("\nElement {0} does not found in the model\n", elementID); throw new Exception("\nProblem in Model Instance - unable to find exchange item\n"); } for (int i = 0; i < n; i++) { XYPolygon modelpolygon = ElementMapper.CreateXYPolygon(baseOut.ElementSet(), i); //int zLayer = Convert.ToInt32(i % base.WMEngine.NumberOfSZLayers); // Points in Polygon are defined as LL, LR, UR, UL (l/l = lower/left, u = upper, r = right ) // Finds the mid x and mid y point in the polygon (assuming rectangular grid) IXYLayerPoint min = new XYLayerPoint(modelpolygon.GetX(0), modelpolygon.GetY(0), 0); IXYLayerPoint max = new XYLayerPoint(modelpolygon.GetX(1), modelpolygon.GetY(3), 0); modelEntities.Add(i, new SpatialDefine(min, max, GeometryTypes.Geometry2D)); } return modelEntities; } /// /// 3D UZ !!!! /// Creates a dictionary with key equal to the model state index and the value the spatial information of that state index. /// /// The geometric type of the exchange itme (2d or 3d) /// The exchange item base output /// the string id of the exchange item. /// private IDictionary GetModelCoordinates3DUZ(GeometryTypes gType, IBaseOutput baseOut, string elementID) { IDictionary modelEntities = new Dictionary(); int n; try { WMEngine.GetElementCount(elementID); n = baseOut.ElementSet().ElementCount; } catch { Console.WriteLine("\nElement {0} does not found in the model\n", elementID); throw new Exception("\nProblem in Model Instance - unable to find exchange item\n"); } // Determines the number of layers in the UZ Grid. int numLayersInGrid = Convert.ToInt32( Math.Round( (double)base.WMEngine.UzGrid.ElementCount/(double)base.WMEngine.UzGrid.BaseGrid.ElementCount ) ); for (int i = 0; i < n; i++) { XYPolygon modelpolygon = ElementMapper.CreateXYPolygon(baseOut.ElementSet(), i); int zLayer = Convert.ToInt32(i % numLayersInGrid); // Points in Polygon are defined as LL, LR, UR, UL (l/l = lower/left, u = upper, r = right ) // Finds the mid x and mid y point in the polygon (assuming rectangular grid) IXYLayerPoint min = new XYLayerPoint(modelpolygon.GetX(0), modelpolygon.GetY(0), zLayer); IXYLayerPoint max = new XYLayerPoint(modelpolygon.GetX(1), modelpolygon.GetY(3), zLayer); modelEntities.Add(i, new SpatialDefine(min, max, GeometryTypes.Geometry3D)); } return modelEntities; } private IDictionary _modelEntities; /// /// 3D SZ !!!! /// Creates a dictionary with key equal to the model state index and the value the spatial information of that state index. /// /// The geometric type of the exchange itme (2d or 3d) /// The exchange item base output /// the string id of the exchange item. /// private IDictionary GetModelCoordinates3DSZ(GeometryTypes gType, IBaseOutput baseOut, string elementID) { //Run Only once - because it's slow if(_modelEntities == null) { _modelEntities = new Dictionary(); int n; try { WMEngine.GetElementCount(elementID); n = baseOut.ElementSet().ElementCount; } catch { Console.WriteLine("\nElement {0} does not found in the model\n", elementID); throw new Exception("\nProblem in Model Instance - unable to find exchange item\n"); } //int numBaseGrid = Convert.ToInt32(Math.Floor((double)n / (double)_mshe.WMEngine.NumberOfSZLayers)); for (int i = 0; i < n; i++) { XYPolygon modelpolygon = ElementMapper.CreateXYPolygon(baseOut.ElementSet(), i); int zLayer = Convert.ToInt32(i % base.WMEngine.NumberOfSZLayers); // Points in Polygon are defined as LL, LR, UR, UL (l/l = lower/left, u = upper, r = right ) // Finds the mid x and mid y point in the polygon (assuming rectangular grid) IXYLayerPoint min = new XYLayerPoint(modelpolygon.GetX(0), modelpolygon.GetY(0), zLayer); IXYLayerPoint max = new XYLayerPoint(modelpolygon.GetX(1), modelpolygon.GetY(3), zLayer); _modelEntities.Add(i, new SpatialDefine(min, max, GeometryTypes.Geometry3D)); } } return _modelEntities; } /// /// Returns a List of ISpatialDefine with assumed 90deg angles. /// The ISpatialDefine is defined by two coordinates and a layer integer. /// A ISpatialDefine can represent a cuboid, rectangle or a point. /// /// /// private IDictionary GetModelCoordinates(string elementIDAndGrid) { char[] delimiterChars = { ',' }; string[] words = elementIDAndGrid.Split(delimiterChars); string elementID = words[0].Trim(); IBaseOutput baseOut = base._outputExchangeItems.First(vID => string.Compare(vID.Id, elementID) == 0); // Get the Grid Type from the elementID (delimiated by a ','); GeometryTypes gType = GetGridType(elementID); if (gType == GeometryTypes.Geometry3D) { return GetModelCoordinates3DSZ(gType, baseOut, elementID); } if (gType == GeometryTypes.Geometry3DUZ) { return GetModelCoordinates3DUZ(gType, baseOut, elementID); } if (gType == GeometryTypes.Geometry2D) { return GetModelCoordinates2D(gType, baseOut, elementID); } else { throw new Exception("The Rest not Implemented"); } } /// /// Distance to Normal calculator. /// Returns the Gaussian normalized distances. /// /// array of doubles of distances to each other /// radius factor /// private double normalCooefs(double dist, double radius) { // Calculated result saved into iteself return Math.Exp(-0.5 * Math.Pow((dist / radius), 2)); } #endregion PrivateMethods } } /* /// /// Calculates Gaussian Localization Mask. /// Given a point on a grid, return an array of doubles (of the same size as the grid) with localization values between 0 and 1. /// /// A model instance from WMengine /// Exchange item string id /// point around which to calcualte the licalization mask /// the distance radius of the Gaussian mask /// private double[] GaussianLocalization(DHI.OpenMI2.MikeShe.WMEngineAccess modelInstance, string elementSetID, IXYLayerPoint point, double locDistance) { IDictionary modelCoord = GetModelCoordinates(elementSetID); double[] localized2D = new double[modelCoord.Count]; if (XYZGeometryTools.IsPointInModelPlain(point, modelCoord)) { for (int i = 0; i < modelCoord.Count; i++) { if (Convert.ToInt32(point.Layer) == modelCoord[i].Layer) { double distance = XYZGeometryTools.CalculatePointToPointDistance2D(modelCoord[i].MidPoint, point); localized2D[i] = normalCooefs(distance, locDistance); } } } return localized2D; } */