Java tutorial
/* * Copyright (C) 2014 Ehsan Roshani * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ package whitebox.stats; import Jama.*; import java.awt.Color; import java.io.BufferedWriter; import java.io.File; import java.io.FileNotFoundException; import java.io.FileWriter; import java.io.IOException; import java.io.PrintWriter; import java.util.ArrayList; import java.util.Date; import java.util.List; import java.util.logging.Level; import java.util.logging.Logger; import java.text.DecimalFormat; import javax.swing.JOptionPane; import java.beans.PropertyChangeSupport; import java.beans.PropertyChangeListener; import net.finmath.optimizer.LevenbergMarquardt; import net.finmath.optimizer.SolverException; import org.jfree.chart.ChartFrame; import org.jfree.chart.JFreeChart; import org.jfree.chart.axis.NumberAxis; import org.jfree.data.xy.XYSeries; import org.jfree.data.xy.XYSeriesCollection; import org.jfree.chart.renderer.xy.*; import org.jfree.chart.plot.*; import org.jfree.chart.renderer.LookupPaintScale; import org.jfree.data.xy.DefaultXYZDataset; import org.jfree.data.xy.XYDataset; import org.jfree.ui.RectangleAnchor; import whitebox.geospatialfiles.ShapeFile; import whitebox.geospatialfiles.WhiteboxRaster; import whitebox.geospatialfiles.WhiteboxRasterBase; import whitebox.geospatialfiles.shapefile.MultiPoint; import whitebox.geospatialfiles.shapefile.MultiPointM; import whitebox.geospatialfiles.shapefile.MultiPointZ; import whitebox.geospatialfiles.shapefile.PointM; import whitebox.geospatialfiles.shapefile.PointZ; import whitebox.geospatialfiles.shapefile.ShapeFileRecord; import whitebox.geospatialfiles.shapefile.ShapeType; import static whitebox.geospatialfiles.shapefile.ShapeType.MULTIPOINT; import static whitebox.geospatialfiles.shapefile.ShapeType.MULTIPOINTM; import static whitebox.geospatialfiles.shapefile.ShapeType.MULTIPOINTZ; import static whitebox.geospatialfiles.shapefile.ShapeType.POINT; import static whitebox.geospatialfiles.shapefile.ShapeType.POINTM; import static whitebox.geospatialfiles.shapefile.ShapeType.POINTZ; import whitebox.geospatialfiles.shapefile.attributes.AttributeTable; import whitebox.geospatialfiles.shapefile.attributes.DBFException; import whitebox.geospatialfiles.shapefile.attributes.DBFField; import whitebox.structures.KdTree; import java.util.Random; import jmetal.util.JMException; import whitebox.geospatialfiles.shapefile.attributes.DBFWriter; import whitebox.interfaces.WhiteboxPluginHost; /** * * @author Ehsan Roshani, Ph.D. Department of Geography University of Guelph * Guelph, Ont. N1G 2W1 CANADA Phone: (519) 824-4120 x53527 Email: * eroshani@uoguelph.ca * * modified by John Lindsay */ public class Kriging { public boolean Anisotropic; public double BandWidth; public double Angle; public double Tolerance; public double resolution; //public double public double bMinX; //Bounding box Minimum X public double bMinY; public double bMaxX; public double bMaxY; private PropertyChangeSupport changes = new PropertyChangeSupport(this); public void SetBoundary(double MinimumX, double MaximumX, double MinimumY, double MaximumY) { bMinX = MinimumX; bMinY = MinimumY; bMaxX = MaximumX; bMaxY = MaximumY; } public double MinX; //Minimum X Coordinate in the Points public double MinY; //Minimum Y Coordinate in the Points public double MaxX; //Maximum X Coordinate in the Points public double MaxY; //Maximum Y Coordinate in the Points public int NumberOfLags; public double LagSize; //public double public KdTree<Double> pointsTree; //This is the point tree which will be filled in the calcPair method public Matrix DistanceMatrix; //The matrix that contains the distance of each known point to all other known points public int nKown; //Number of known points //public double[][] Points; //Array of points location x=0, y = 1, z = 2 public double MaximumDistance; public bin[][] BinSurface; //n*3 matrix to store all the bins // public class point // { // public double x; // public double y; // public double z; // } // public point point(double x , double y, double z){ // point p = new point(); // p.x = x; // p.y = y; // p.z = z; // return p; // } public List<KrigingPoint> points = new ArrayList(); public class pair { int FirstP; int SecondP; double Distance; double Direction; double MomentI; //Moment of Inertia double VerDistance; //Vertical Distance (Y Axes) double HorDistance; //Horizontal Distance (X Axes) } //List<bin> bins = new ArrayList(); public bin[][] bins; // = new bin[] public class bin { double GridHorDistance; double GridVerDistance; double HorDistance; double VerDistance; double Distance; double Value; double Weight; int Size; } public List<pair> Pairs = new ArrayList(); private KdTree<Double> PairsTree; public SemivariogramType SemiVariogramModel; public enum SemivariogramType { SPHERICAL, EXPONENTIAL, GAUSSIAN } public double Range; //SemiVariogram a value public double Sill; //SemiVariogram h value public double Nugget; //Semivariogram Nugget value public boolean ConsiderNugget; //If nugget should be considered or not public SemivariogramType SVType; private int nthSVariogram; //this is the nth SV for Anisotropic private double[] x;// This is the x value for fitting the theoritical SV public class Variogram { public double Range; public double Sill; public double Nugget; public SemivariogramType Type; public double mse; } public void DrawShapeFile(String outputFile, List<KrigingPoint> pnts) throws DBFException, IOException { File file = new File(outputFile); if (file.exists()) { file.delete(); } // set up the output files of the shapefile and the dbf ShapeFile output = new ShapeFile(outputFile, ShapeType.POINT); DBFField fields[] = new DBFField[2]; fields[0] = new DBFField(); fields[0].setName("FID"); fields[0].setDataType(DBFField.DBFDataType.NUMERIC); fields[0].setFieldLength(10); fields[0].setDecimalCount(0); fields[1] = new DBFField(); fields[1].setName("Z"); fields[1].setDataType(DBFField.DBFDataType.NUMERIC); fields[1].setFieldLength(10); fields[1].setDecimalCount(3); String DBFName = output.getDatabaseFile(); DBFWriter writer = new DBFWriter(new File(DBFName)); writer.setFields(fields); int numPointsInFile = pnts.size(); double x, y, z; for (int a = 0; a < numPointsInFile; a++) { x = pnts.get(a).x; y = pnts.get(a).y; z = pnts.get(a).z; whitebox.geospatialfiles.shapefile.Point wbGeometry = new whitebox.geospatialfiles.shapefile.Point(x, y); output.addRecord(wbGeometry); Object[] rowData = new Object[2]; rowData[0] = new Double(a + 1); rowData[1] = new Double(z); writer.addRecord(rowData); } output.write(); writer.write(); } /** * Calculates the bin values based on sector classification and for * isotropic model * * @param Range */ void calcBins4Sec(double Range) { int ad = 0; if (Range % this.LagSize == 0) { ad = 0; } if (!this.Anisotropic) { bins = new bin[(int) Math.ceil(Range / this.LagSize) + ad][1]; int r = 0; for (int i = 0; i < Pairs.size(); i++) { if (Pairs.get(i).Distance < Range && Pairs.get(i).HorDistance >= 0) { r = (int) Math.floor(Pairs.get(i).Distance / LagSize); if (bins[r][0] == null) { bin bb = new bin(); bins[r][0] = bb; } bins[r][0].Distance += Pairs.get(i).Distance; bins[r][0].Value += Pairs.get(i).MomentI; bins[r][0].Size++; } } for (int i = 0; i < bins.length; i++) { if (bins[i][0] == null) { bin bb = new bin(); bins[i][0] = bb; } bins[i][0].Distance = bins[i][0].Distance / bins[i][0].Size; bins[i][0].Value = bins[i][0].Value / bins[i][0].Size; } } } /** * Calculates the bin values based on sector classification and for * AnIsotropic model * * @param Range * @param Angle * @param Tolerance * @param BandWidth */ void calcBins4Sec(double Range, double Angle, double Tolerance, double BandWidth) { int ad = 0; if (Range % this.LagSize == 0) { ad = 0; } double width = 0; if (this.Anisotropic) { bins = new bin[(int) Math.ceil(Range / this.LagSize) + ad][1]; int r = 0; for (int i = 0; i < Pairs.size(); i++) { boolean tt = Between(Angle, Tolerance, Pairs.get(i).Direction); width = Pairs.get(i).Distance * Math.cos((Math.PI / 2) - Angle + Pairs.get(i).Direction); if (tt && Pairs.get(i).Distance < Range && Math.abs(width) <= BandWidth) { r = (int) Math.floor(Pairs.get(i).Distance / LagSize); if (bins[r][0] == null) { bin bb = new bin(); bins[r][0] = bb; } bins[r][0].Distance += Pairs.get(i).Distance; bins[r][0].Value += Pairs.get(i).MomentI; bins[r][0].Size++; } } for (int i = 0; i < bins.length; i++) { if (bins[i][0] == null) { bin bb = new bin(); bins[i][0] = bb; } bins[i][0].Distance = bins[i][0].Distance / bins[i][0].Size; bins[i][0].Value = bins[i][0].Value / bins[i][0].Size; } } //========================== } /** * Checks to see if a pair is located in a sector of interest or not * * @param Angle * @param Tolerance * @param Direction * @return */ private boolean Between(double Angle, double Tolerance, double Direction) { boolean flag = false; double la = (Angle - Tolerance); if (la < 0) { la = 2 * Math.PI + la; flag = true; } double ha = (Angle + Tolerance); if (ha >= 2 * Math.PI) { ha = ha - 2 * Math.PI; flag = true; } if (flag) { if (Direction >= ha && Direction <= la) { return false; } else { return true; } } else { if (Direction >= la && Direction <= ha) { return true; } else { return false; } } } /** * Calculates the Bin list for SV Map * * @param Range */ void CalcBins4Map(double Range) { //bins = new bin[this.NumberOfLags][this.NumberOfLags]; //bins Category on the axies //2 . 1 //3 4 //Only 1 and 4 are calculated the rest are mirror int ad = 0; if (Range % this.LagSize == 0) { ad = 0; } bin[][] bins1 = new bin[(int) Math.ceil(Range / this.LagSize) + ad][(int) Math.ceil(Range / this.LagSize + ad)]; bin[][] bins4 = new bin[(int) Math.ceil(Range / this.LagSize) + ad][(int) Math.ceil(Range / this.LagSize + ad)]; bin[][] bins1c = new bin[(int) Math.ceil(Range / this.LagSize) + ad][(int) Math.ceil(Range / this.LagSize + ad)]; bin[][] bins4c = new bin[(int) Math.ceil(Range / this.LagSize) + ad][(int) Math.ceil(Range / this.LagSize + ad)]; BinSurface = new bin[2 * ((int) Math.ceil(Range / this.LagSize) + ad)][2 * ((int) Math.ceil(Range / this.LagSize + ad))]; //double radious =Math.sqrt(2*this.LagSize*this.LagSize); double radious = this.LagSize * 2 / Math.sqrt(2); double halfLagSize = this.LagSize; List<pair> prs = new ArrayList(); double w = 0; for (int r = 0; r < bins1.length; r++) { for (int c = 0; c < bins1[r].length; c++) { if (bins1[r][c] == null) { bin bb = new bin(); bin bbc = new bin(); bins1[r][c] = bb; bins1c[r][c] = bbc; } bins1[r][c].GridHorDistance = 0.5 * this.LagSize + c * this.LagSize; bins1[r][c].GridVerDistance = 0.5 * this.LagSize + r * this.LagSize; bins1c[r][c].GridHorDistance = -0.5 * this.LagSize - c * this.LagSize; bins1c[r][c].GridVerDistance = -0.5 * this.LagSize - r * this.LagSize; double[] center = new double[] { bins1[r][c].GridVerDistance, bins1[r][c].GridHorDistance }; prs = getBinNNPairs4Map(PairsTree, center, halfLagSize, radious); for (int n = 0; n < prs.size(); n++) { bins1[r][c].HorDistance += prs.get(n).HorDistance; bins1[r][c].VerDistance += prs.get(n).VerDistance; w = (1 - (Math.abs(bins1[r][c].GridHorDistance - prs.get(n).HorDistance) / this.LagSize)) * (1 - (Math.abs(bins1[r][c].GridVerDistance - prs.get(n).VerDistance) / this.LagSize)); bins1[r][c].Weight += w; bins1[r][c].Value += prs.get(n).MomentI * w; bins1[r][c].Size += 1; bins1c[r][c].HorDistance += prs.get(n).HorDistance; bins1c[r][c].VerDistance += prs.get(n).VerDistance; bins1c[r][c].Weight += w; bins1c[r][c].Value += prs.get(n).MomentI * w; bins1c[r][c].Size += 1; } } } for (int i = 0; i < bins1.length; i++) { for (int j = 0; j < bins1[i].length; j++) { if (bins1[i][j] == null) { bin bb = new bin(); bins1[i][j] = bb; bins1[i][j].HorDistance = i * this.LagSize; bins1[i][j].VerDistance = j * this.LagSize; bins1[i][j].Value = -1; bin bbc = new bin(); bins1c[i][j] = bbc; bins1c[i][j].HorDistance = -i * this.LagSize; bins1c[i][j].VerDistance = -j * this.LagSize; bins1c[i][j].Value = -1; } else { bins1[i][j].HorDistance = bins1[i][j].HorDistance / bins1[i][j].Size; bins1[i][j].VerDistance = bins1[i][j].VerDistance / bins1[i][j].Size; bins1[i][j].Value = bins1[i][j].Value / bins1[i][j].Weight; bins1c[i][j].HorDistance = bins1c[i][j].HorDistance / bins1c[i][j].Size; bins1c[i][j].VerDistance = bins1c[i][j].VerDistance / bins1c[i][j].Size; bins1c[i][j].Value = bins1c[i][j].Value / bins1c[i][j].Weight; } //System.out.println( (0.5*this.LagSize + j*this.LagSize) + " , " + (0.5*this.LagSize+i*this.LagSize) + " , " + // Bins1[i][j].HorDistance + " , " + Bins1[i][j].VerDistance + " , " + Bins1[i][j].Value); } } //========================== for (int r = 0; r < bins4.length; r++) { for (int c = 0; c < bins4[r].length; c++) { if (bins4[r][c] == null) { bin bb = new bin(); bin bbc = new bin(); bins4[r][c] = bb; bins4c[r][c] = bbc; } bins4[r][c].GridHorDistance = 0.5 * this.LagSize + c * this.LagSize; bins4[r][c].GridVerDistance = -0.5 * this.LagSize - r * this.LagSize; bins4c[r][c].GridHorDistance = -0.5 * this.LagSize - c * this.LagSize; bins4c[r][c].GridVerDistance = 0.5 * this.LagSize + r * this.LagSize; double[] center = new double[] { bins4[r][c].GridVerDistance, bins4[r][c].GridHorDistance }; prs = getBinNNPairs4Map(PairsTree, center, halfLagSize, radious); for (int n = 0; n < prs.size(); n++) { bins4[r][c].HorDistance += prs.get(n).HorDistance; bins4[r][c].VerDistance += prs.get(n).VerDistance; w = (1 - (Math.abs(bins4[r][c].GridHorDistance - prs.get(n).HorDistance) / this.LagSize)) * (1 - (Math.abs(bins4[r][c].GridVerDistance - prs.get(n).VerDistance) / this.LagSize)); bins4[r][c].Weight += w; bins4[r][c].Value += prs.get(n).MomentI * w; bins4[r][c].Size += 1; bins4c[r][c].HorDistance += prs.get(n).HorDistance; bins4c[r][c].VerDistance += prs.get(n).VerDistance; bins4c[r][c].Weight += w; bins4c[r][c].Value += prs.get(n).MomentI * w; bins4c[r][c].Size += 1; } } } for (int i = 0; i < bins4.length; i++) { for (int j = 0; j < bins4[i].length; j++) { if (bins4[i][j] == null) { bin bb = new bin(); bins4[i][j] = bb; bins4[i][j].HorDistance = i * this.LagSize; bins4[i][j].VerDistance = j * this.LagSize; bins4[i][j].Value = -1; bin bbc = new bin(); bins4c[i][j] = bbc; bins4c[i][j].HorDistance = -i * this.LagSize; bins4c[i][j].VerDistance = -j * this.LagSize; bins4c[i][j].Value = -1; } else { bins4[i][j].HorDistance = bins4[i][j].HorDistance / bins4[i][j].Size; bins4[i][j].VerDistance = bins4[i][j].VerDistance / bins4[i][j].Size; bins4[i][j].Value = bins4[i][j].Value / bins4[i][j].Weight; bins4c[i][j].HorDistance = bins4c[i][j].HorDistance / bins4c[i][j].Size; bins4c[i][j].VerDistance = bins4c[i][j].VerDistance / bins4c[i][j].Size; bins4c[i][j].Value = bins4c[i][j].Value / bins4c[i][j].Weight; } //System.out.println( (0.5*this.LagSize + j*this.LagSize) + " , " + (0.5*this.LagSize+i*this.LagSize) + " , " + // bins1[i][j].HorDistance + " , " + bins1[i][j].VerDistance + " , " + bins1[i][j].Value); } } int stI = BinSurface.length / 2; int stJ = BinSurface[0].length / 2; //bins1c = bins1.clone(); for (int i = 0; i < bins1.length; i++) { for (int j = 0; j < bins1[i].length; j++) { BinSurface[stI + i][stJ + j] = bins1[i][j]; BinSurface[stI - 1 - i][stJ - 1 - j] = bins1c[i][j]; } } stI = BinSurface.length / 2; stJ = BinSurface[0].length / 2; for (int i = 0; i < bins4.length; i++) { for (int j = 0; j < bins4[i].length; j++) { BinSurface[stI - 1 - i][stJ + j] = bins4[i][j]; BinSurface[stI + i][stJ - 1 - j] = bins4c[i][j]; } } // for (int i = 0; i < BinSurface.length; i++) { // for (int j = 0; j < BinSurface[i].length; j++) { // System.out.println(BinSurface[i][j].GridHorDistance + " , " + BinSurface[i][j].GridVerDistance // + " , " + BinSurface[i][j].HorDistance+ " , " + BinSurface[i][j].VerDistance // + " , " + BinSurface[i][j].Value); // } // } int resd = 0; } /** * Draw Semivariogram surface map and also draw the search are if * Anisotropic * * @param Radius * @param AnIsotropic */ public void DrawSemivariogramSurface(double Radius, boolean AnIsotropic) { double[][] data = new double[3][BinSurface.length * BinSurface[0].length]; int n = 0; double max = Double.MIN_VALUE; for (int i = 0; i < BinSurface.length; i++) { for (int j = 0; j < BinSurface[i].length; j++) { data[0][n] = BinSurface[i][j].GridHorDistance; data[1][n] = BinSurface[i][j].GridVerDistance; if ((Math.pow(data[0][n], 2) + Math.pow(data[1][n], 2)) <= Radius * Radius && !Double.isNaN(BinSurface[i][j].Value)) { data[2][n] = BinSurface[i][j].Value; if (max < data[2][n]) { max = data[2][n]; } } else { data[2][n] = -1; } n++; } } DefaultXYZDataset dataset = new DefaultXYZDataset(); dataset.addSeries("Value", data); NumberAxis xAxis = new NumberAxis(); xAxis.setStandardTickUnits(NumberAxis.createIntegerTickUnits()); xAxis.setLowerMargin(0.0); xAxis.setUpperMargin(0.0); NumberAxis yAxis = new NumberAxis(); yAxis.setStandardTickUnits(NumberAxis.createIntegerTickUnits()); yAxis.setLowerMargin(0.0); yAxis.setUpperMargin(0.0); XYBlockRenderer renderer = new XYBlockRenderer(); renderer.setBlockWidth(LagSize); renderer.setBlockHeight(LagSize); renderer.setBlockAnchor(RectangleAnchor.CENTER); LookupPaintScale paintScale = new LookupPaintScale(0, max, Color.white); double colorRange = max / 6; //double colorRange = 23013; paintScale.add(0.0, Color.blue); paintScale.add(1 * colorRange, Color.green); paintScale.add(2 * colorRange, Color.cyan); paintScale.add(3 * colorRange, Color.yellow); paintScale.add(4 * colorRange, Color.ORANGE); paintScale.add(5 * colorRange, Color.red); renderer.setPaintScale(paintScale); XYPlot plot = new XYPlot(dataset, xAxis, yAxis, renderer); plot.setBackgroundPaint(Color.lightGray); plot.setDomainGridlinesVisible(false); plot.setRangeGridlinePaint(Color.white); if (AnIsotropic) { CombinedRangeXYPlot combinedrangexyplot = new CombinedRangeXYPlot(); XYSeries seriesT1 = new XYSeries("1"); XYSeriesCollection AngleCollct = new XYSeriesCollection(); double bw = BandWidth; double r = bw / Math.sin(Tolerance); if (r > Radius) { bw = Radius * Math.sin(Tolerance); r = Radius; } seriesT1.add(r * Math.cos(Angle + Tolerance), r * Math.sin(Angle + Tolerance)); if ((double) Math.round(Math.sin(Angle) * 10000) / 10000 != 0) { if ((double) Math.round(Math.cos(Angle) * 10000) / 10000 != 0) { double a = (1 + Math.pow(Math.tan(Angle), 2)); double b = 2 * bw / Math.sin(Angle) * Math.pow(Math.tan(Angle), 2); double c = Math.pow(Math.tan(Angle), 2) * Math.pow(bw / Math.sin(Angle), 2) - Math.pow(Radius, 2); double x1 = (-b + Math.sqrt(Math.pow(b, 2) - 4 * a * c)) / (2 * a); double y1 = Math.tan(Angle) * (x1 + bw / Math.sin(Angle)); double x2 = (-b - Math.sqrt(Math.pow(b, 2) - 4 * a * c)) / (2 * a); double y2 = Math.tan(Angle) * (x2 + bw / Math.sin(Angle)); double d1 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x1), 2)) + (Math.pow((Radius * Math.sin(Angle) - y1), 2))); double d2 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x2), 2)) + (Math.pow((Radius * Math.sin(Angle) - y2), 2))); if (d1 < d2) { seriesT1.add(x1, y1); } else { seriesT1.add(x2, y2); } } else { double x1 = -bw * Math.sin(Angle); double y1 = Math.sqrt(Math.pow(Radius, 2) - Math.pow(x1, 2)); double y2 = -Math.sqrt(Math.pow(Radius, 2) - Math.pow(x1, 2)); double d1 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x1), 2)) + (Math.pow((Radius * Math.sin(Angle) - y1), 2))); double d2 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x1), 2)) + (Math.pow((Radius * Math.sin(Angle) - y2), 2))); if (d1 < d2) { seriesT1.add(x1, y1); } else { seriesT1.add(x1, y2); } } } else { double y1 = bw * Math.cos(Angle); double x1 = Math.sqrt(Math.pow(Radius, 2) - Math.pow(y1, 2)); double x2 = -Math.sqrt(Math.pow(Radius, 2) - Math.pow(y1, 2)); double d1 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x1), 2)) + (Math.pow((Radius * Math.sin(Angle) - y1), 2))); double d2 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x2), 2)) + (Math.pow((Radius * Math.sin(Angle) - y1), 2))); if (d1 < d2) { seriesT1.add(x1, y1); } else { seriesT1.add(x2, y1); } } AngleCollct.addSeries(seriesT1); XYSeries seriesT2 = new XYSeries("2"); seriesT2.add(r * Math.cos(Angle + Tolerance), r * Math.sin(Angle + Tolerance)); seriesT2.add(0.0, 0.0); AngleCollct.addSeries(seriesT2); XYSeries seriesT3 = new XYSeries("3"); seriesT3.add(Radius * Math.cos(Angle), Radius * Math.sin(Angle)); seriesT3.add(0, 0); AngleCollct.addSeries(seriesT3); XYSeries seriesT4 = new XYSeries("4"); seriesT4.add(r * Math.cos(Angle - Tolerance), r * Math.sin(Angle - Tolerance)); seriesT4.add(0, 0); AngleCollct.addSeries(seriesT4); XYSeries seriesT5 = new XYSeries("5"); seriesT5.add(r * Math.cos(Angle - Tolerance), r * Math.sin(Angle - Tolerance)); if ((double) Math.round(Math.sin(Angle) * 10000) / 10000 != 0) { if ((double) Math.round(Math.cos(Angle) * 10000) / 10000 != 0) { double a = (1 + Math.pow(Math.tan(Angle), 2)); double b = -2 * bw / Math.sin(Angle) * Math.pow(Math.tan(Angle), 2); double c = Math.pow(Math.tan(Angle), 2) * Math.pow(bw / Math.sin(Angle), 2) - Math.pow(Radius, 2); double x1 = (-b + Math.sqrt(Math.pow(b, 2) - 4 * a * c)) / (2 * a); double y1 = Math.tan(Angle) * (x1 - bw / Math.sin(Angle)); double x2 = (-b - Math.sqrt(Math.pow(b, 2) - 4 * a * c)) / (2 * a); double y2 = Math.tan(Angle) * (x2 - bw / Math.sin(Angle)); double d1 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x1), 2)) + (Math.pow((Radius * Math.sin(Angle) - y1), 2))); double d2 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x2), 2)) + (Math.pow((Radius * Math.sin(Angle) - y2), 2))); if (d1 < d2) { seriesT5.add(x1, y1); } else { seriesT5.add(x2, y2); } } else { double x1 = bw * Math.sin(Angle); double y1 = Math.sqrt(Math.pow(Radius, 2) - Math.pow(x1, 2)); double y2 = -Math.sqrt(Math.pow(Radius, 2) - Math.pow(x1, 2)); double d1 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x1), 2)) + (Math.pow((Radius * Math.sin(Angle) - y1), 2))); double d2 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x1), 2)) + (Math.pow((Radius * Math.sin(Angle) - y2), 2))); if (d1 < d2) { seriesT5.add(x1, y1); } else { seriesT5.add(x1, y2); } } } else { double y1 = -bw * Math.cos(Angle); double x1 = Math.sqrt(Math.pow(Radius, 2) - Math.pow(y1, 2)); double x2 = -Math.sqrt(Math.pow(Radius, 2) - Math.pow(y1, 2)); double d1 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x1), 2)) + (Math.pow((Radius * Math.sin(Angle) - y1), 2))); double d2 = Math.sqrt((Math.pow((Radius * Math.cos(Angle) - x2), 2)) + (Math.pow((Radius * Math.sin(Angle) - y1), 2))); if (d1 < d2) { seriesT5.add(x1, y1); } else { seriesT5.add(x2, y1); } } AngleCollct.addSeries(seriesT5); plot.setDataset(1, AngleCollct); XYLineAndShapeRenderer lineshapRend = new XYLineAndShapeRenderer(true, false); for (int i = 0; i < AngleCollct.getSeriesCount(); i++) { //plot.getRenderer().setSeriesPaint(i , Color.BLUE); lineshapRend.setSeriesPaint(i, Color.BLACK); } plot.setRenderer(1, lineshapRend); combinedrangexyplot.add(plot); } plot.setDatasetRenderingOrder(DatasetRenderingOrder.FORWARD); JFreeChart chart = new JFreeChart("Semivariogram Surface", plot); chart.removeLegend(); chart.setBackgroundPaint(Color.white); // create and display a frame... ChartFrame frame = new ChartFrame("", chart); frame.pack(); //frame.setSize(100, 50); frame.setVisible(true); } /** * This method uses NSGA algorithm to fit the Semi Variogram * * @param semiType * @param n * @return */ Variogram TheoryVariogramNSGA(SemivariogramType semiType, int n) { // Set solver parameters double[] y = new double[bins.length]; for (int i = 0; i < y.length; i++) { y[i] = bins[i][n].Value; } int nNan = 0; for (int i = 0; i < y.length; i++) { if (!Double.isNaN(y[i])) { nNan++; } } x = new double[nNan]; double[] y2 = new double[nNan]; int ntmp = 0; for (int i = 0; i < y.length; i++) { if (!Double.isNaN(y[i])) { y2[ntmp] = y[i]; x[ntmp] = bins[i][nthSVariogram].Distance; ntmp++; } } y = y2; double[][] pnts = new double[y.length][2]; for (int i = 0; i < y.length; i++) { pnts[i][1] = y[i]; pnts[i][0] = x[i]; } Variogram var = new Variogram(); var.Type = semiType; SemivariogramCurveFitter svcf = new SemivariogramCurveFitter(); try { var = svcf.Run(pnts, semiType, ConsiderNugget); } catch (JMException ex) { //Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } catch (SecurityException ex) { //Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } catch (IOException ex) { //Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } catch (ClassNotFoundException ex) { //Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } return var; } /** * * @param semiType * @param n is the nth sector for anisotropic * @return */ Variogram TheoryVariogram(SemivariogramType semiType, int n) { SVType = semiType; nthSVariogram = n; //x = xValue ; //TheoryVariogramNSGA(SVType,n); LevenbergMarquardt optimizer = new LevenbergMarquardt() { // Override your objective function here @Override public void setValues(double[] parameters, double[] values) { //parameters[0] = sill, parameters[1] Range, parameters[2] nugget // double [] x = new double[values.length]; // for (int i = 0; i < values.length; i++) { // x[i]=bins[i][ nthSVariogram].Distance; // } switch (SVType) { case EXPONENTIAL: for (int i = 0; i < x.length; i++) { if (x[i] != 0) { values[i] = (ConsiderNugget ? parameters[2] : 0) + parameters[0] * (1 - Math.exp(-x[i] / parameters[1])); } else { values[i] = 0; } } break; case GAUSSIAN: for (int i = 0; i < x.length; i++) { if (x[i] != 0) { values[i] = (ConsiderNugget ? parameters[2] : 0) + parameters[0] * (1 - Math.exp(-(Math.pow(x[i], 2)) / (Math.pow(parameters[1], 2)))); } else { values[i] = 0; } } break; case SPHERICAL: for (int i = 0; i < x.length; i++) { if (x[0] > parameters[1]) { values[i] = (ConsiderNugget ? parameters[2] : 0) + parameters[0]; } else if (0 < x[0] && x[0] <= parameters[1]) { values[i] = (ConsiderNugget ? parameters[2] : 0) + parameters[0] * (1.5 * x[i] / parameters[1] - 0.5 * Math.pow((x[i] / parameters[1]), 3)); } else { values[i] = 0; } } break; } } }; // Set solver parameters double[] y = new double[bins.length]; for (int i = 0; i < y.length; i++) { y[i] = bins[i][n].Value; } int nNan = 0; for (int i = 0; i < y.length; i++) { if (!Double.isNaN(y[i])) { nNan++; } } x = new double[nNan]; double[] y2 = new double[nNan]; int ntmp = 0; for (int i = 0; i < y.length; i++) { if (!Double.isNaN(y[i])) { y2[ntmp] = y[i]; x[ntmp] = bins[i][nthSVariogram].Distance; ntmp++; } } y = y2; double[] iniPar = new double[y.length]; double[] w = new double[y.length]; for (int i = 0; i < y.length; i++) { iniPar[i] = 1; w[i] = 1; } double tmp = 0; int tmpN = 0; for (int i = 0; i < y.length; i++) { if (!Double.isNaN(y[i])) { tmp += y[i]; tmpN++; } else { w[i] = 0; } } iniPar[1] = this.LagSize; iniPar[0] = tmp / tmpN; optimizer.setInitialParameters(iniPar); optimizer.setWeights(w); optimizer.setMaxIteration(100); optimizer.setErrorTolerance(0.1); optimizer.setTargetValues(y); try { optimizer.run(); } catch (SolverException ex) { Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } double[] bestParameters = optimizer.getBestFitParameters(); // this.Sill = bestParameters[0]; // this.Range=bestParameters[1]; // this.Nugget = (ConsiderNugget ? bestParameters[2] : 0 ); Variogram var = new Variogram(); var.Sill = bestParameters[0]; var.Range = bestParameters[1]; if (var.Sill < 0) { var.Sill = 0; } var.Nugget = (ConsiderNugget ? bestParameters[2] : 0); var.Type = semiType; return var; } /** * Calculates the Theoretical SV value to be drawn on the SV * * @param Distance * @param vario * @return */ public double getTheoreticalSVValue(double Distance, Variogram vario) { double res = 0.0; switch (vario.Type) { case EXPONENTIAL: if (Distance != 0) { res = vario.Nugget + vario.Sill * (1 - Math.exp(-Distance / vario.Range)); } else { res = 0; } break; case GAUSSIAN: if (Distance != 0) { res = vario.Nugget + vario.Sill * (1 - Math.exp(-3 * (Math.pow(Distance, 2)) / (Math.pow(vario.Range, 2)))); } else { res = 0; } break; case SPHERICAL: if (Distance > vario.Range) { res = vario.Nugget + vario.Sill; } else if (0 < Distance && Distance <= vario.Range) { res = vario.Nugget + vario.Sill * (1.5 * Distance / vario.Range - 0.5 * Math.pow((Distance / vario.Range), 3)); } else { res = 0; } break; } return res; } /** * Calculates the points for drawing the theoretical variogram * * @param SVType * @return */ double[][] CalcTheoreticalSVValues(Variogram vario, double MaximumDisplyDistanst) { double[][] res = new double[2 * NumberOfLags + 1][2]; //0=X, 1= Y for (int i = 0; i < res.length; i++) { res[i][0] = i * MaximumDisplyDistanst / (2 * NumberOfLags); switch (vario.Type) { case EXPONENTIAL: if (res[i][0] != 0) { res[i][1] = vario.Nugget + vario.Sill * (1 - Math.exp(-res[i][0] / vario.Range)); } else { res[i][1] = vario.Nugget; } break; case GAUSSIAN: if (res[i][0] != 0) { res[i][1] = vario.Nugget + vario.Sill * (1 - Math.exp(-3 * (Math.pow(res[i][0], 2)) / (Math.pow(vario.Range, 2)))); } else { res[i][1] = vario.Nugget; } break; case SPHERICAL: if (res[i][0] > vario.Range) { res[i][1] = vario.Nugget + vario.Sill; } else if (res[i][0] > 0 && res[i][0] <= vario.Range) { res[i][1] = vario.Nugget + vario.Sill * (1.5 * res[i][0] / vario.Range - 0.5 * Math.pow((res[i][0] / vario.Range), 3)); } else { res[i][1] = vario.Nugget; } break; } } return res; } /** * Reads points from a shapefile retrieving the attribute data from the * z-values. This method assumes that the shapefile is of z type. * * @param inputFile String containing the file name of the shapefile. */ public void readPointFile(String inputFile) { int fieldNum = 0; WhiteboxRasterBase.DataType dataType = WhiteboxRasterBase.DataType.INTEGER; boolean useRecID = false; // initialize the shapefile input ShapeFile input = null; try { input = new ShapeFile(inputFile); } catch (IOException ex) { System.out.println(ex.getMessage().toString()); Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } ShapeType shapeType = input.getShapeType(); if (shapeType != ShapeType.POINTZ && shapeType != ShapeType.MULTIPOINTZ) { //showFeedback("The input shapefile must be of a 'point' data type."); JOptionPane.showMessageDialog(null, "The input shapefile must be of a 'point' data type."); return; } Object[] data = null; double[][] geometry; points = new ArrayList<>(); PointZ ptz; MultiPointZ mptz; double x, y, z; double[][] point; for (ShapeFileRecord record : input.records) { if (shapeType.getBaseType() == ShapeType.POINT) { ptz = (PointZ) (record.getGeometry()); z = ptz.getZ(); x = ptz.getX(); y = ptz.getY(); points.add(new KrigingPoint(x, y, z)); } else if (shapeType.getBaseType() == ShapeType.MULTIPOINT) { mptz = (MultiPointZ) (record.getGeometry()); point = record.getGeometry().getPoints(); double[] zArray = mptz.getzArray(); for (int p = 0; p < point.length; p++) { x = point[p][0]; y = point[p][1]; z = zArray[p]; points.add(new KrigingPoint(x, y, z)); } } } } /** * Reads points from a shapefile retrieving the attribute data from a field * within the shapefile's database. * * @param inputFile String containing the file name of the shapefile. * @param fieldName String specifying the name of a field in the database * file. * @return */ public void readPointFile(String inputFile, String fieldName) { int fieldNum = 0; WhiteboxRasterBase.DataType dataType = WhiteboxRasterBase.DataType.INTEGER; boolean useRecID = false; // initialize the shapefile input ShapeFile input = null; try { input = new ShapeFile(inputFile); } catch (IOException ex) { System.out.println(ex.getMessage().toString()); Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } if (input.getShapeType() != ShapeType.POINT && input.getShapeType() != ShapeType.POINTZ && input.getShapeType() != ShapeType.POINTM && input.getShapeType() != ShapeType.MULTIPOINT && input.getShapeType() != ShapeType.MULTIPOINTZ && input.getShapeType() != ShapeType.MULTIPOINTM) { //showFeedback("The input shapefile must be of a 'point' data type."); JOptionPane.showMessageDialog(null, "The input shapefile must be of a 'point' data type."); return; } /////////////// // what type of data is contained in fieldName? AttributeTable reader = input.getAttributeTable(); //new DBFReader(input.getDatabaseFile()); int numberOfFields = reader.getFieldCount(); for (int i = 0; i < numberOfFields; i++) { DBFField field = reader.getField(i); if (field.getName().equals(fieldName)) { fieldNum = i; if (field.getDataType() == DBFField.DBFDataType.NUMERIC || field.getDataType() == DBFField.DBFDataType.FLOAT) { if (field.getDecimalCount() == 0) { dataType = WhiteboxRasterBase.DataType.INTEGER; } else { dataType = WhiteboxRasterBase.DataType.FLOAT; } } else { useRecID = true; } } } if (fieldNum < 0) { useRecID = true; } ////////////////////// Object[] data = null; double[][] geometry; points = new ArrayList<>(); for (ShapeFileRecord record : input.records) { try { data = reader.nextRecord(); } catch (DBFException ex) { Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } geometry = getXYFromShapefileRecord(record); for (int i = 0; i < geometry.length; i++) { KrigingPoint p = new KrigingPoint(geometry[i][0], geometry[i][1], Double.valueOf(data[fieldNum].toString())); points.add(p); } } } public void setPoints(List<KrigingPoint> krigingPoints) { points = krigingPoints; } public List<KrigingPoint> getPoints() { return points; } private double[][] getXYFromShapefileRecord(ShapeFileRecord record) { double[][] ret; ShapeType shapeType = record.getShapeType(); switch (shapeType) { case POINT: whitebox.geospatialfiles.shapefile.Point recPoint = (whitebox.geospatialfiles.shapefile.Point) (record .getGeometry()); ret = new double[1][2]; ret[0][0] = recPoint.getX(); ret[0][1] = recPoint.getY(); break; case POINTZ: PointZ recPointZ = (PointZ) (record.getGeometry()); ret = new double[1][2]; ret[0][0] = recPointZ.getX(); ret[0][1] = recPointZ.getY(); break; case POINTM: PointM recPointM = (PointM) (record.getGeometry()); ret = new double[1][2]; ret[0][0] = recPointM.getX(); ret[0][1] = recPointM.getY(); break; case MULTIPOINT: MultiPoint recMultiPoint = (MultiPoint) (record.getGeometry()); return recMultiPoint.getPoints(); case MULTIPOINTZ: MultiPointZ recMultiPointZ = (MultiPointZ) (record.getGeometry()); return recMultiPointZ.getPoints(); case MULTIPOINTM: MultiPointM recMultiPointM = (MultiPointM) (record.getGeometry()); return recMultiPointM.getPoints(); default: ret = new double[1][2]; ret[1][0] = -1; ret[1][1] = -1; break; } return ret; } /** * It calculates the location of each grid cell. the resolution should be * set before calling this method * * @return a point list */ public List<KrigingPoint> calcInterpolationPoints() { double north, south, east, west; int nrows, ncols; double northing, easting; west = bMinX - 0.5 * resolution; north = bMaxY + 0.5 * resolution; nrows = (int) (Math.ceil((north - bMinY) / resolution)); ncols = (int) (Math.ceil((bMaxX - west) / resolution)); south = north - nrows * resolution; east = west + ncols * resolution; int row, col; List<KrigingPoint> pnts = new ArrayList(); // Create the whitebox raster object. double halfResolution = resolution / 2; for (row = 0; row < nrows; row++) { for (col = 0; col < ncols; col++) { easting = (col * resolution) + (west + halfResolution); northing = (north - halfResolution) - (row * resolution); pnts.add(new KrigingPoint(easting, northing, 0.0)); } } return pnts; } public void addPropertyChangeListener(PropertyChangeListener listener) { changes.addPropertyChangeListener(listener); } WhiteboxPluginHost host = null; public void setPluginHost(WhiteboxPluginHost host) { this.host = host; } /** * Produces the output raster * * @param outputRaster * @param pnts * @param drawKrigingVariance */ public void buildRaster(String outputRaster, List<KrigingPoint> pnts, boolean drawKrigingVariance) { double north, south, east, west; int nrows, ncols; double northing, easting; west = bMinX - 0.5 * resolution; north = bMaxY + 0.5 * resolution; nrows = (int) (Math.ceil((north - bMinY) / resolution)); ncols = (int) (Math.ceil((bMaxX - west) / resolution)); south = north - nrows * resolution; east = west + ncols * resolution; String outputHeader = outputRaster; FileWriter fw = null; BufferedWriter bw = null; PrintWriter out = null; String str1; double noData = -32768; // see if the output files already exist, and if so, delete them. if ((new File(outputHeader)).exists()) { (new File(outputHeader)).delete(); (new File(outputHeader.replace(".dep", ".tas"))).delete(); } try { // create the whitebox header file. fw = new FileWriter(outputHeader, false); bw = new BufferedWriter(fw); out = new PrintWriter(bw, true); str1 = "Min:\t" + Double.toString(Integer.MAX_VALUE); out.println(str1); str1 = "Max:\t" + Double.toString(Integer.MIN_VALUE); out.println(str1); str1 = "North:\t" + Double.toString(north); out.println(str1); str1 = "South:\t" + Double.toString(south); out.println(str1); str1 = "East:\t" + Double.toString(east); out.println(str1); str1 = "West:\t" + Double.toString(west); out.println(str1); str1 = "Cols:\t" + Integer.toString(ncols); out.println(str1); str1 = "Rows:\t" + Integer.toString(nrows); out.println(str1); str1 = "Data Type:\t" + "float"; out.println(str1); str1 = "Z Units:\t" + "not specified"; out.println(str1); str1 = "XY Units:\t" + "not specified"; out.println(str1); str1 = "Projection:\t" + "not specified"; out.println(str1); str1 = "Data Scale:\tcontinuous"; out.println(str1); str1 = "Preferred Palette:\t" + "rgb.pal"; out.println(str1); str1 = "NoData:\t" + noData; out.println(str1); if (java.nio.ByteOrder.nativeOrder() == java.nio.ByteOrder.LITTLE_ENDIAN) { str1 = "Byte Order:\t" + "LITTLE_ENDIAN"; } else { str1 = "Byte Order:\t" + "BIG_ENDIAN"; } out.println(str1); out.close(); } catch (Exception e) { return; } int row, col; // Create the whitebox raster object. WhiteboxRaster image = new WhiteboxRaster(outputHeader, "rw"); double halfResolution = resolution / 2; int nn = 0; int progress; int oldProgress = -1; for (row = 0; row < nrows; row++) { for (col = 0; col < ncols; col++) { easting = (col * resolution) + (west + halfResolution); northing = (north - halfResolution) - (row * resolution); if (!drawKrigingVariance) { image.setValue(row, col, pnts.get(nn).z); } else { image.setValue(row, col, pnts.get(nn).v); } nn++; } progress = (int) (100f * row / (nrows - 1)); if (progress > oldProgress) { changes.firePropertyChange("progress", oldProgress, progress); //host.updateProgress("Interpolating Data:", progress); oldProgress = progress; } } image.addMetadataEntry("Created by the Kriging Interpolation Tool."); image.addMetadataEntry("Created on " + new Date()); image.addMetadataEntry("Semivariogram Model = " + SemiVariogramModel); image.addMetadataEntry("Range = " + Range); image.addMetadataEntry("Sill = " + Sill); image.addMetadataEntry("Nugget = " + Nugget); image.close(); } /** * Gets the variogram and unknown point list and returns the interpolated * values for the known points This is to calculate the predicted value for * each known point, the result would be used for cross validation. * * @param variogram * @param pnts * @return */ public List<KrigingPoint> CrossValidationPoints(Variogram variogram, List<KrigingPoint> pnts, int NumberOfNearestPoints) { double[] res = new double[NumberOfNearestPoints]; double[][] D = new double[NumberOfNearestPoints + 1][1]; List<KrigingPoint> NNPoitns = new ArrayList(); List<KrigingPoint> outPnts = new ArrayList(); for (int n = 0; n < pnts.size(); n++) { NNPoitns = getNNpoints(this.pointsTree, pnts.get(n), NumberOfNearestPoints + 1); for (int ni = 0; ni < NumberOfNearestPoints + 1; ni++) { if (pnts.get(n).x == NNPoitns.get(ni).x && pnts.get(n).y == NNPoitns.get(ni).y && pnts.get(n).z == NNPoitns.get(ni).z) { NNPoitns.remove(ni); break; } } double[][] C = CalcConstantCoef(variogram, NNPoitns); double[] tm = CalcVariableCoef(variogram, pnts.get(n), NNPoitns); ///------------ for (int i = 0; i < tm.length; i++) { D[i][0] = tm[i]; } //double[][] d = {{1,2,3},{4,5,6,},{7,8,10}}; Matrix tmp = Matrix.constructWithCopy(C); Matrix VariableCoef = Matrix.constructWithCopy(D); Matrix w = null; boolean flag = false; try { w = tmp.solve(VariableCoef); double[][] Wi = w.getArray(); double s = 0; for (int i = 0; i < Wi.length - 1; i++) { s = s + Wi[i][0] * NNPoitns.get(i).z; } KrigingPoint pnt = new KrigingPoint(pnts.get(n).x, pnts.get(n).y, s); outPnts.add(pnt); //pnts.get(n).z = s; //res[n]=s; s = 0; } catch (Exception ex) { SingularValueDecomposition svd = tmp.svd(); Matrix u = svd.getU(); Matrix s = svd.getS(); Matrix v = svd.getV(); //u.print(u.getRowDimension(), u.getColumnDimension()); //s.print(s.getRowDimension(), s.getColumnDimension()); //v.print(v.getRowDimension(), v.getColumnDimension()); int rrr = svd.rank(); double[][] stemp = s.getArray(); for (int nn = 0; nn < stemp.length; nn++) { if (stemp[nn][nn] > 0.03) { stemp[nn][nn] = 1 / stemp[nn][nn]; } else { stemp[nn][nn] = 0; } } Matrix sp = new Matrix(stemp); w = v.times(sp).times(u.transpose()).times(VariableCoef); //Matrix test = tmp.times(w).minus(VariableCoef); double[][] Wi = w.getArray(); double ss = 0; for (int i = 0; i < Wi.length - 1; i++) { ss = ss + Wi[i][0] * NNPoitns.get(i).z; } KrigingPoint pnt = new KrigingPoint(pnts.get(n).x + 1, pnts.get(n).y, ss); outPnts.add(pnt); //pnts.get(n).z = ss; ss = 0; } } return outPnts; } public void interpolateRaster(Variogram variogram, int numberOfNearestPoints, WhiteboxRaster raster, boolean mapError) { int row, col; int rows = raster.getNumberRows(); int cols = raster.getNumberColumns(); double easting, northing; double e; double[][] D = new double[numberOfNearestPoints + 1][1]; List<KrigingPoint> nnPoints = new ArrayList(); int progress; int oldProgress = -1; for (row = 0; row < rows; row++) { for (col = 0; col < cols; col++) { easting = raster.getXCoordinateFromColumn(col); northing = raster.getYCoordinateFromRow(row); KrigingPoint pnt = new KrigingPoint(easting, northing, 0.0d); nnPoints = getNNpoints(this.pointsTree, pnt, numberOfNearestPoints); double[][] C = CalcConstantCoef(variogram, nnPoints); double[] tm = CalcVariableCoef(variogram, pnt, nnPoints); for (int i = 0; i < tm.length; i++) { D[i][0] = tm[i]; } Matrix tmp = Matrix.constructWithCopy(C); Matrix VariableCoef = Matrix.constructWithCopy(D); Matrix w = null; //boolean flag = false; try { double vs = 0; w = tmp.solve(VariableCoef); double[][] Wi = w.getArray(); double s = 0; for (int i = 0; i < Wi.length - 1; i++) { s = s + Wi[i][0] * nnPoints.get(i).z; vs = vs + Wi[i][0] * D[i][0]; } e = vs + Wi[Wi.length - 1][0]; if (!mapError) { raster.setValue(row, col, s); } else { raster.setValue(row, col, e); } } catch (Exception ex) { SingularValueDecomposition svd = tmp.svd(); Matrix u = svd.getU(); Matrix s = svd.getS(); Matrix v = svd.getV(); int rrr = svd.rank(); double[][] stemp = s.getArray(); for (int nn = 0; nn < stemp.length; nn++) { if (stemp[nn][nn] > 0.003) { stemp[nn][nn] = 1 / stemp[nn][nn]; } else { stemp[nn][nn] = 0; } } Matrix sp = new Matrix(stemp); w = v.times(sp).times(u.transpose()).times(VariableCoef); //Matrix test = tmp.times(w).minus(VariableCoef); double[][] Wi = w.getArray(); double ss = 0; double vs = 0; for (int i = 0; i < Wi.length - 1; i++) { ss = ss + Wi[i][0] * nnPoints.get(i).z; vs = vs + Wi[i][0] * D[i][0]; } e = vs + Wi[Wi.length - 1][0]; if (!mapError) { raster.setValue(row, col, ss); } else { raster.setValue(row, col, e); } } } progress = (int) (100f * row / (rows - 1.0)); if (progress > oldProgress) { changes.firePropertyChange("progress", oldProgress, progress); if (host != null) { host.updateProgress("Interpolating Data:", progress); } oldProgress = progress; } } } public void interpolateRaster(Variogram variogram, int numberOfNearestPoints, WhiteboxRaster raster, WhiteboxRaster errorRaster) { int row, col; int rows = raster.getNumberRows(); int cols = raster.getNumberColumns(); double easting, northing; double e; double[][] D = new double[numberOfNearestPoints + 1][1]; List<KrigingPoint> nnPoints = new ArrayList(); int progress; int oldProgress = -1; for (row = 0; row < rows; row++) { for (col = 0; col < cols; col++) { easting = raster.getXCoordinateFromColumn(col); northing = raster.getYCoordinateFromRow(row); KrigingPoint pnt = new KrigingPoint(easting, northing, 0.0d); nnPoints = getNNpoints(this.pointsTree, pnt, numberOfNearestPoints); double[][] C = CalcConstantCoef(variogram, nnPoints); double[] tm = CalcVariableCoef(variogram, pnt, nnPoints); for (int i = 0; i < tm.length; i++) { D[i][0] = tm[i]; } Matrix tmp = Matrix.constructWithCopy(C); Matrix VariableCoef = Matrix.constructWithCopy(D); Matrix w = null; //boolean flag = false; try { double vs = 0; w = tmp.solve(VariableCoef); double[][] Wi = w.getArray(); double s = 0; for (int i = 0; i < Wi.length - 1; i++) { s = s + Wi[i][0] * nnPoints.get(i).z; vs = vs + Wi[i][0] * D[i][0]; } e = vs + Wi[Wi.length - 1][0]; raster.setValue(row, col, s); errorRaster.setValue(row, col, e); } catch (Exception ex) { SingularValueDecomposition svd = tmp.svd(); Matrix u = svd.getU(); Matrix s = svd.getS(); Matrix v = svd.getV(); int rrr = svd.rank(); double[][] stemp = s.getArray(); for (int nn = 0; nn < stemp.length; nn++) { if (stemp[nn][nn] > 0.003) { stemp[nn][nn] = 1 / stemp[nn][nn]; } else { stemp[nn][nn] = 0; } } Matrix sp = new Matrix(stemp); w = v.times(sp).times(u.transpose()).times(VariableCoef); //Matrix test = tmp.times(w).minus(VariableCoef); double[][] Wi = w.getArray(); double ss = 0; double vs = 0; for (int i = 0; i < Wi.length - 1; i++) { ss = ss + Wi[i][0] * nnPoints.get(i).z; vs = vs + Wi[i][0] * D[i][0]; } e = vs + Wi[Wi.length - 1][0]; raster.setValue(row, col, ss); errorRaster.setValue(row, col, e); } } progress = (int) (100f * row / (rows - 1.0)); if (progress > oldProgress) { changes.firePropertyChange("progress", oldProgress, progress); if (host != null) { host.updateProgress("Interpolating Data:", progress); } oldProgress = progress; } } } /** * Gets the variogram and unknown point list and returns the interpolated * values for the unknown points It also calculates the Kriging Variance and * sets the KrigingPoint.V * * @param variogram * @param pnts * @param numberOfNearestPoints * @return */ public List<KrigingPoint> interpolatePoints(Variogram variogram, List<KrigingPoint> pnts, int numberOfNearestPoints) { //double[] res = new double[numberOfNearestPoints]; double[][] D = new double[numberOfNearestPoints + 1][1]; List<KrigingPoint> nnPoints = new ArrayList(); List<KrigingPoint> outPnts = new ArrayList(); int progress; int oldProgress = -1; int numPoints = pnts.size(); for (int n = 0; n < numPoints; n++) { nnPoints = getNNpoints(this.pointsTree, pnts.get(n), numberOfNearestPoints); double[][] C = CalcConstantCoef(variogram, nnPoints); double[] tm = CalcVariableCoef(variogram, pnts.get(n), nnPoints); for (int i = 0; i < tm.length; i++) { D[i][0] = tm[i]; } Matrix tmp = Matrix.constructWithCopy(C); Matrix VariableCoef = Matrix.constructWithCopy(D); Matrix w = null; //boolean flag = false; try { double vs = 0; w = tmp.solve(VariableCoef); double[][] Wi = w.getArray(); double s = 0; for (int i = 0; i < Wi.length - 1; i++) { s = s + Wi[i][0] * nnPoints.get(i).z; vs = vs + Wi[i][0] * D[i][0]; } KrigingPoint pnt = new KrigingPoint(pnts.get(n).x, pnts.get(n).y, s); pnt.v = vs + Wi[Wi.length - 1][0]; // if (pnt.v <= 0) { // pnt.v = pnt.v; // } outPnts.add(pnt); //pnts.get(n).z = s; //res[n]=s; s = 0; } catch (Exception ex) { SingularValueDecomposition svd = tmp.svd(); Matrix u = svd.getU(); Matrix s = svd.getS(); Matrix v = svd.getV(); int rrr = svd.rank(); double[][] stemp = s.getArray(); for (int nn = 0; nn < stemp.length; nn++) { if (stemp[nn][nn] > 0.003) { stemp[nn][nn] = 1 / stemp[nn][nn]; } else { stemp[nn][nn] = 0; } } Matrix sp = new Matrix(stemp); w = v.times(sp).times(u.transpose()).times(VariableCoef); //Matrix test = tmp.times(w).minus(VariableCoef); double[][] Wi = w.getArray(); double ss = 0; double vs = 0; for (int i = 0; i < Wi.length - 1; i++) { ss = ss + Wi[i][0] * nnPoints.get(i).z; vs = vs + Wi[i][0] * D[i][0]; } KrigingPoint pnt = new KrigingPoint(pnts.get(n).x, pnts.get(n).y, ss); pnt.v = vs + Wi[Wi.length - 1][0]; if (pnt.v <= 0) { pnt.v = pnt.v; for (int i = 0; i < nnPoints.size(); i++) { System.out.println(nnPoints.get(i).x + " " + nnPoints.get(i).y + " " + nnPoints.get(i).z); } } outPnts.add(pnt); //pnts.get(n).z = ss; ss = 0; } progress = (int) (100f * n / (numPoints - 1.0)); if (progress > oldProgress) { changes.firePropertyChange("progress", oldProgress, progress); if (host != null) { host.updateProgress("Interpolating Data:", progress); } oldProgress = progress; } } return outPnts; } /** * Returns the list of Pairs which are in the Nearest Neighborhood of the * bin center point * * @param Tree * @param entry (y,x) * @param HalfBinSize * @param Range is the search radius * @return */ private List<pair> getBinNNPairs4Map(KdTree<Double> Tree, double[] entry, double BinSize, double Range) { List<KdTree.Entry<Double>> results; results = Tree.neighborsWithinRange(entry, Range); List<pair> res = new ArrayList(); double xd = 0; double yd = 0; for (int i = 0; i < results.size(); i++) { xd = Math.sqrt(Math.pow((Pairs.get(results.get(i).value.intValue()).HorDistance - entry[1]), 2)); yd = Math.sqrt(Math.pow((Pairs.get(results.get(i).value.intValue()).VerDistance - entry[0]), 2)); if (xd <= BinSize && yd <= BinSize) { res.add(Pairs.get(results.get(i).value.intValue())); } } return res; } /** * Returns the list of nearest neighbor points * * @param Tree * @param pnt * @param numPointsToUse * @return */ private List<KrigingPoint> getNNpoints(KdTree<Double> Tree, KrigingPoint pnt, int numPointsToUse) { double[] entry; //double[] outentry; entry = new double[] { pnt.y, pnt.x }; List<KdTree.Entry<Double>> results; results = Tree.nearestNeighbor(entry, numPointsToUse, false); List<KrigingPoint> pnts = new ArrayList(); List<KrigingPoint> res = new ArrayList(); for (int i = 0; i < results.size(); i++) { //KrigingPoint tmp = new KrigingPoint(); //int id = results.get(i).value.intValue(); res.add(points.get(results.get(i).value.intValue())); } return res; } /** * calculates the D matrix for Kriging system * * @param variogram * @param p is the unknown point * @param NNPoints is the list of nearest neighbor points * @return */ private double[] CalcVariableCoef(Variogram variogram, KrigingPoint p, List<KrigingPoint> NNPoints) { int n = NNPoints.size(); double[] mat = new double[n + 1]; double dist = 0.0; for (int i = 0; i < n; i++) { dist = Math.sqrt(Math.abs(Math.pow(NNPoints.get(i).x - p.x, 2)) + Math.abs(Math.pow(NNPoints.get(i).y - p.y, 2))); mat[i] = getTheoreticalSVValue(dist, variogram); } mat[n] = 1; return mat; } /** * This prepares the known points matrix for ordinary Kriging * * @param variogarm * @return */ private double[][] CalcConstantCoef(Variogram variogarm, List<KrigingPoint> NNPoints) { int n = NNPoints.size(); double[][] mat = new double[n + 1][n + 1]; double dist = 0.0; for (int i = 0; i < n; i++) { for (int j = i; j < n; j++) { dist = Math.sqrt(Math.abs(Math.pow(NNPoints.get(i).x - NNPoints.get(j).x, 2)) + Math.abs(Math.pow(NNPoints.get(i).y - NNPoints.get(j).y, 2))); mat[i][j] = getTheoreticalSVValue(dist, variogarm); mat[j][i] = mat[i][j]; } } for (int i = 0; i < n; i++) { mat[i][n] = 1; mat[n][i] = 1; } // // String s= new String(); // try { // PrintWriter pr = new PrintWriter("G:\\test.txt"); // for (int i = 0; i < mat.length; i++) { // for (int j = 0; j < mat.length; j++) { // s = s + "," + mat[i][j]; // } // pr.println(s); // s = ""; // } // pr.close(); // } catch (FileNotFoundException ex) { // Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); // } // return mat; } /** * This just to use when the semivariogram model is provided by the user. * (Kriging Optimizer) */ void BuildPointTree() { pointsTree = new KdTree.SqrEuclid<Double>(2, new Integer(this.points.size())); double[] entry; for (int i = 0; i < this.points.size(); i++) { entry = new double[] { this.points.get(i).y, this.points.get(i).x }; pointsTree.addPoint(entry, (double) i); } } /** * Creates the pairs list based on sector classification. calcs the distance * and moment of inertia for each pair It also calculates the min and max * points and boundary It also build the KDTree object to be used with the * Kriging */ void calPairs4Sec() throws FileNotFoundException { MaximumDistance = 0; MinX = Double.POSITIVE_INFINITY; MinY = Double.POSITIVE_INFINITY; MaxX = Double.NEGATIVE_INFINITY; MaxY = Double.NEGATIVE_INFINITY; pointsTree = new KdTree.SqrEuclid<>(2, new Integer(this.points.size())); //PairsTree = new KdTree.SqrEuclid<Double>(2, new Integer(this.Points.size()*(this.Points.size()-1)/2)); PairsTree = new KdTree.SqrEuclid<>(2, new Integer(this.points.size() * (this.points.size()))); double[] entry; double[] pairentry; // String s = new String(); double dx = 0; double dy = 0; for (int i = 0; i < this.points.size(); i++) { if (this.points.get(i).x < MinX) { MinX = this.points.get(i).x; } if (this.points.get(i).y < MinY) { MinY = this.points.get(i).y; } if (this.points.get(i).x > MaxX) { MaxX = this.points.get(i).x; } if (this.points.get(i).y > MaxY) { MaxY = this.points.get(i).y; } entry = new double[] { this.points.get(i).y, this.points.get(i).x }; pointsTree.addPoint(entry, (double) i); for (int j = 0; j < this.points.size(); j++) { pair pr = new pair(); if (i != j) { pr.FirstP = i; pr.SecondP = j; pr.Distance = Math.sqrt(Math.pow((points.get(i).x - points.get(j).x), 2) + Math.pow((points.get(i).y - points.get(j).y), 2)); pr.HorDistance = (points.get(j).x - points.get(i).x); pr.VerDistance = (points.get(j).y - points.get(i).y); if (MaximumDistance < pr.Distance) { MaximumDistance = pr.Distance; } dx = points.get(j).x - points.get(i).x; dy = points.get(j).y - points.get(i).y; if (dx != 0) { if ((dx > 0 && dy >= 0)) { pr.Direction = Math.atan(dy / dx); } if (dx < 0 && dy >= 0) { pr.Direction = Math.atan(dy / dx) + Math.PI; } if (dx > 0 && dy < 0) { pr.Direction = Math.atan(dy / dx) + 2 * Math.PI; } if (dx < 0 && dy < 0) { pr.Direction = Math.atan(dy / dx) + Math.PI; ; } } else { if (dy >= 0) { pr.Direction = Math.PI / 2; } else { pr.Direction = 3 * Math.PI / 2; } } pr.MomentI = Math.pow((points.get(i).z - points.get(j).z), 2) / 2; Pairs.add(pr); pairentry = new double[] { pr.VerDistance, pr.HorDistance }; PairsTree.addPoint(pairentry, (double) Pairs.size() - 1.0); } } } // pw.close(); //LagSize = MaximumDistance/NumberOfLags; bMaxX = MaxX; bMaxY = MaxY; bMinX = MinX; bMinY = MinY; } /** * Creates the pairs list based for Map classification * * @throws FileNotFoundException */ void CalPairs4Map() throws FileNotFoundException { MaximumDistance = 0; MinX = Double.POSITIVE_INFINITY; MinY = Double.POSITIVE_INFINITY; MaxX = Double.NEGATIVE_INFINITY; MaxY = Double.NEGATIVE_INFINITY; pointsTree = new KdTree.SqrEuclid<Double>(2, new Integer(this.points.size())); PairsTree = new KdTree.SqrEuclid<Double>(2, new Integer(this.points.size() * (this.points.size() - 1) / 2)); double[] entry; double[] pairentry; // String s= new String(); // PrintWriter pw ; // pw = new PrintWriter("G:\\test.txt"); double dx = 0; double dy = 0; for (int i = 0; i < this.points.size(); i++) { if (this.points.get(i).x < MinX) { MinX = this.points.get(i).x; } if (this.points.get(i).y < MinY) { MinY = this.points.get(i).y; } if (this.points.get(i).x > MaxX) { MaxX = this.points.get(i).x; } if (this.points.get(i).y > MaxY) { MaxY = this.points.get(i).y; } entry = new double[] { this.points.get(i).y, this.points.get(i).x }; pointsTree.addPoint(entry, (double) i); for (int j = 0; j < this.points.size(); j++) { pair pr = new pair(); if (points.get(i).x <= points.get(j).x && i != j) { pr.FirstP = i; pr.SecondP = j; pr.Distance = Math.sqrt(Math.pow((points.get(i).x - points.get(j).x), 2) + Math.pow((points.get(i).y - points.get(j).y), 2)); pr.HorDistance = (points.get(j).x - points.get(i).x); pr.VerDistance = (points.get(j).y - points.get(i).y); if (MaximumDistance < pr.Distance) { MaximumDistance = pr.Distance; } dx = points.get(j).x - points.get(i).x; dy = points.get(j).y - points.get(i).y; if (dx != 0) { if ((dx > 0 && dy >= 0)) { pr.Direction = Math.atan(dy / dx); } if (dx < 0 && dy >= 0) { pr.Direction = Math.atan(dy / dx) + Math.PI; } if (dx > 0 && dy < 0) { pr.Direction = Math.atan(dy / dx) + 2 * Math.PI; } if (dx < 0 && dy < 0) { pr.Direction = Math.atan(dy / dx) + Math.PI; ; } } else { if (dy >= 0) { pr.Direction = Math.PI / 2; } else { pr.Direction = 3 * Math.PI / 2; } } pr.MomentI = Math.pow((points.get(i).z - points.get(j).z), 2) / 2; Pairs.add(pr); pairentry = new double[] { pr.VerDistance, pr.HorDistance }; PairsTree.addPoint(pairentry, (double) Pairs.size() - 1.0); // s = Double.toString(pr.Distance) + "," + Double.toString(pr.Direction)+ // "," + Double.toString(pr.MomentI)+ // "," + Double.toString(pr.HorDistance)+ // ","+Double.toString(pr.VerDistance)+ // "," + Integer.toString(pr.FirstP)+ // "," + Integer.toString(pr.SecondP); // // pw.println(s); } } } // pw.close(); //LagSize = MaximumDistance/NumberOfLags; bMaxX = MaxX; bMaxY = MaxY; bMinX = MinX; bMinY = MinY; } /** * It gets the semivariogram type and bins list and draw a graph for them * TheoryVariogram should be called first * * @param bins * @param variogram */ public void DrawSemivariogram(bin[][] bins, Variogram variogram) { XYSeriesCollection sampleCollct = new XYSeriesCollection(); XYSeries series = new XYSeries("Sample Variogram"); // for (Iterator<bin> i = bins.iterator(); i.hasNext(); ) // { // series.add(bins.get(j).Distance,bins.get(j).Value); // i.next(); // j++; // } XYLineAndShapeRenderer xylineshapRend = new XYLineAndShapeRenderer(false, true); CombinedRangeXYPlot combinedrangexyplot = new CombinedRangeXYPlot(); for (int i = 0; i < bins[0].length; i++) { for (int k = 0; k < bins.length; k++) { if (!Double.isNaN(bins[k][i].Value)) { series.add(bins[k][i].Distance, bins[k][i].Value); } } sampleCollct.addSeries(series); double[][] res = CalcTheoreticalSVValues(variogram, series.getMaxX()); XYSeries seriesTSV = new XYSeries("Theoretical Variogram"); for (int l = 0; l < res.length; l++) { seriesTSV.add(res[l][0], res[l][1]); } XYSeriesCollection theorCollct = new XYSeriesCollection(); theorCollct.addSeries(seriesTSV); XYDataset xydataset = sampleCollct; XYPlot xyplot1 = new XYPlot(xydataset, new NumberAxis(), null, xylineshapRend); xyplot1.setDataset(1, theorCollct); XYLineAndShapeRenderer lineshapRend = new XYLineAndShapeRenderer(true, false); xyplot1.setRenderer(1, lineshapRend); xyplot1.setDatasetRenderingOrder(DatasetRenderingOrder.FORWARD); combinedrangexyplot.add(xyplot1); } DecimalFormat df = new DecimalFormat("###,##0.000"); String title = "Semivariogram (RMSE = " + df.format(Math.sqrt(variogram.mse)) + ")"; JFreeChart chart = new JFreeChart(title, JFreeChart.DEFAULT_TITLE_FONT, combinedrangexyplot, true); // JFreeChart chart = ChartFactory.createScatterPlot( // "Semivariogram", // chart title // "Distance", // x axis label // "Moment of Inertia", // y axis label // result, // data // PlotOrientation.VERTICAL, // true, // include legend // true, // tooltips // false // urls // ); // create and display a frame... ChartFrame frame = new ChartFrame("Semivariogram", chart); frame.pack(); frame.setVisible(true); } /** * This is the main method to classify the pairs for the map and to calc the * bin average on the map * * @param Type * @param DistanseRatio * @param NumberOfLags * @param Anisotropic */ public void calcBinSurface(SemivariogramType Type, double DistanseRatio, int NumberOfLags, boolean Anisotropic) { this.NumberOfLags = NumberOfLags; try { CalPairs4Map(); } catch (FileNotFoundException ex) { Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } if (this.LagSize == 0) { this.LagSize = (this.MaximumDistance * DistanseRatio) / this.NumberOfLags; } CalcBins4Map(this.LagSize * this.NumberOfLags); } /** * This Calculates the Sill and Range Value for the Theoretical Semi * Variogram Points list should be filled first This function fills the Sill * and Range in the Kriging object * * @param type * @param distanceRatio is the ratio of the maximum distance in point to the * maximum distance of the variogram * @param numberOfLags * @param anisotropic * @param useNSGA * @return */ public Variogram getSemivariogram(SemivariogramType type, double distanceRatio, int numberOfLags, boolean anisotropic, boolean useNSGA) { this.NumberOfLags = numberOfLags; try { calPairs4Sec(); } catch (FileNotFoundException ex) { Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); } if (this.LagSize == 0) { this.LagSize = (this.MaximumDistance * distanceRatio) / this.NumberOfLags; } int n = 0; if (!anisotropic) { n = 0; calcBins4Sec(this.LagSize * this.NumberOfLags); } else { n = 0; calcBins4Sec(this.LagSize * this.NumberOfLags, this.Angle, this.Tolerance, this.BandWidth); } return (useNSGA) ? TheoryVariogramNSGA(type, n) : TheoryVariogram(type, n); } public Variogram getSemivariogram(SemivariogramType type, double range, double sill, double nugget, boolean anisotropic) { Variogram var = new Variogram(); var.Type = type; var.Range = range; var.Sill = sill; var.Nugget = nugget; return var; } /** * Randomly selects the n points from the entered point list This is not a * necessary method to use but with large point list (More than 1000 points) * It is better to apply it * * @param pnts * @param n * @return */ public List<KrigingPoint> RandomizePoints(List<KrigingPoint> pnts, int n) { Random rnd = new Random(); List<KrigingPoint> res = new ArrayList(); double drnd = 0.0; for (int i = 0; i < n; i++) { res.add(pnts.get(rnd.nextInt(pnts.size()))); } return res; } public static void main(String[] args) { try { //ChartPanel(createChart(createDataset())); Kriging k = new Kriging(); k.readPointFile("/Users/johnlindsay/Documents/Data/Krigging Test Data/test.shp", "Z"); k.ConsiderNugget = false; k.LagSize = 50; k.NumberOfLags = 100; k.Anisotropic = true; k.Tolerance = Math.PI / 4; k.BandWidth = 5000; PrintWriter pw = new PrintWriter("/Users/johnlindsay/Documents/Data/Krigging Test Data/test.txt"); for (int i = 0; i < 13; i++) { k.Angle = Math.PI / 12.0 * i; Variogram var = k.getSemivariogram(SemivariogramType.SPHERICAL, 1, k.NumberOfLags, true, true); pw.println(k.Angle + "," + var.Range + "," + var.Sill + "," + var.Nugget); pw.flush(); System.out.println((i + 1) + " of 12"); } System.out.println("Done!"); } catch (Exception e) { } // k.ConsiderNugget = false; // k.Points = k.ReadPointFile( // "G:\\Optimized Sensory Network\\PALS\\20120607\\test.shp","Z"); // k.LagSize = 2000; // k.Anisotropic = false; // Variogram var = k.SemiVariogram(SemiVariogramType.Spherical, 0.27, 50,false, true); // // //var.Range = 4160.672768; // //var.Sill = 1835.571948; // // k.resolution = 914; // k.BMinX = 588450 + k.resolution/2; // k.BMaxX = 601246 - k.resolution/2; // k.BMinY = 5474650 + k.resolution/2; // k.BMaxY = 5545942 - k.resolution/2; // // List<KrigingPoint> outPnts = k.calcInterpolationPoints() ; // outPnts = k.InterpolatePoints(var, outPnts, 5); // k.BuildRaster("G:\\Optimized Sensory Network\\PALS\\20120607\\Pnts60.dep", outPnts,false); // k.BuildRaster("G:\\Optimized Sensory Network\\PALS\\20120607\\PntsVar60.dep", outPnts,true); // // k.DrawSemiVariogram(k.bins, var); //k.calcBinSurface(SemiVariogramType.Spherical, 1, 99,false); //k.DrawSemiVariogramSurface(k.LagSize*(k.NumberOfLags), false); // k.LagSize = 2500; // k.Anisotropic = false; //// k.Angle = Math.PI*3/4; //// k.Tolerance = Math.PI/4; //// k.BandWidth = 3*k.LagSize; // for (int i = 0; i < 50; i++) { // k.Points = k.ReadPointFile("G:\\Optimized Sensory Network\\PALS\\Pals Shapefiles\\PALS_TA_20120607_HiAlt_v100.shp","h"); // k.Points = k.RandomizePoints(k.Points, 200); // Variogram var = k.SemiVariogram(SemiVariogramType.Spherical, 1, 25,false); // k.resolution = 900; // List<KrigingPoint> pnts = k.calcInterpolationPoints(); // pnts = k.InterpolatePoints(var, pnts, 5); // k.BuildRaster("G:\\Optimized Sensory Network\\PALS\\Pals Shapefiles\\PALS_TA_20120607_HiAlt_v100"+ i +".dep", pnts); // } //k.DrawSemiVariogram(k.bins, var); //k.calcBinSurface(SemiVariogramType.Spherical, 0.27, 25,false); //k.DrawSemiVariogramSurface(k.LagSize*(k.NumberOfLags),false); //k.Points = k.ReadPointFile("G:\\Papers\\AGU 2013\\Sample\\Sample.shp","V"); //k.Points = k.ReadPointFile("G:\\Papers\\AGU 2013\\WakerLake\\WakerLake.shp","V"); // String s= new String(); // PrintWriter pw = null ; // try { // pw = new PrintWriter("G:\\test.txt"); // } catch (FileNotFoundException ex) { // Logger.getLogger(Kriging.class.getName()).log(Level.SEVERE, null, ex); // } // // for (int i = 0; i < 500; i++) { // Kriging k = new Kriging(); // k.ConsiderNugget = false; // k.LagSize = 50; // k.Anisotropic = true; // k.Angle = Math.PI*3/4; // k.Tolerance = Math.PI/4; // k.BandWidth = 3*k.LagSize; // k.Points = k.ReadPointFile("G:\\Optimized Sensory Network\\PALS\\AGU\\SV_Test.shp","v"); // k.Points = k.RandomizePoints(k.Points, 500); // Variogram var = k.SemiVariogram(SemiVariogramType.Exponential, 0.27, 99,true); // s = var.Range + " , " + var.Sill; // // // pw.println(s); // pw.flush(); // k = null; // } // pw.close();; //k.resolution = 2.5; //k.DrawSemiVariogram(k.bins, var); //List<point> pnts = k.calcInterpolationPoints(); // var.Range = 50; // var.Sill = 104843.2; // var.Type = SemiVariogramType.Exponential; // // //pnts = k.InterpolatePoints(var, pnts,10); //k.BuildRaster("G:\\Papers\\AGU 2013\\WakerLake\\WakerLakeOut15.dep", pnts); //k.calcBinSurface(SemiVariogramType.Exponential, 0.27, 99,false); //k.DrawSemiVariogramSurface(k.LagSize*(k.NumberOfLags),true); //Kriging.point p = k.point(65, 137, 0); // List<Kriging.point> pnts = new ArrayList(); // pnts.add(p); } }