Bulletin17B.java [src/java/m/cfa/flood] Revision: Date:
package m.cfa.flood;
import m.cfa.CircleDrawer;
import m.cfa.DoubleArray;
import m.cfa.DoubleMath;
import m.cfa.Graphing;
import java.awt.BasicStroke;
import java.awt.Color;
import java.awt.geom.Ellipse2D;
import java.awt.geom.Rectangle2D;
import java.io.File;
import java.io.IOException;
import java.text.NumberFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.Comparator;
import java.util.Iterator;
import org.apache.commons.math.ArgumentOutsideDomainException;
import org.codehaus.jettison.json.JSONArray;
import org.codehaus.jettison.json.JSONException;
import org.codehaus.jettison.json.JSONObject;
import org.jfree.chart.ChartUtilities;
import org.jfree.chart.JFreeChart;
import org.jfree.chart.annotations.XYAnnotation;
import org.jfree.chart.annotations.XYDrawableAnnotation;
import org.jfree.chart.annotations.XYPointerAnnotation;
import org.jfree.chart.annotations.XYTextAnnotation;
import org.jfree.chart.annotations.XYTitleAnnotation;
import org.jfree.chart.axis.LogarithmicAxis;
import org.jfree.chart.axis.NumberAxis;
import org.jfree.chart.axis.NumberTickUnit;
import org.jfree.chart.axis.TickUnits;
import org.jfree.chart.axis.ValueAxis;
import org.jfree.chart.block.BlockBorder;
import org.jfree.chart.plot.IntervalMarker;
import org.jfree.chart.plot.XYPlot;
import org.jfree.chart.renderer.xy.XYItemRenderer;
import org.jfree.chart.renderer.xy.XYLineAndShapeRenderer;
import org.jfree.chart.title.LegendTitle;
import org.jfree.data.xy.XYDataset;
import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;
import org.jfree.ui.Layer;
import org.jfree.ui.RectangleAnchor;
import org.jfree.ui.RectangleEdge;
import org.jfree.ui.TextAnchor;
class sort1_smallToLarge implements Comparator<double[]>{
//Compares the first entry of sorts smallest to largest
public int compare(final double[] entry1, final double[] entry2) {
double value1 = entry1[0];
double value2 = entry2[0];
//Compare and return the second entries
return Double.compare(value1,value2);
}
}
class sort2_smallToLarge implements Comparator<double[]>{
//Compares the second entry of sorts largest to smallest
public int compare(final double[] entry1, final double[] entry2) {
double value1 = entry1[1];
double value2 = entry2[1];
//Compare and return the second entries
return Double.compare(value1,value2);
}
}
class sort2_largeToSmall implements Comparator<double[]>{
//Compares the second entry of sorts largest to smallest
public int compare(final double[] entry1, final double[] entry2) {
double value1 = entry1[1];
double value2 = entry2[1];
//Compare and return the second entries
if(Double.compare(value1,value2) > 0){
return -1;
}else if(Double.compare(value1,value2) < 0){
return 1;
}else{
return 0;
}
}
}
/**
* Last Updated: 31-May-2017
* @author Tyler Wible (converted from Matlab to Java retaining comments from the Matlab code and the original Matlab code segments as comments before/after their corresponding Java code segments)
* @author Jeff Burkey (originally written in Matlab)
* @since Java version: 16-May-2012
*/
public class Bulletin17B {
//These are global variables for offsets for the custom graph
double[] tickLabelOffset = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
double tickMarkOffset = 0.95;
double yAxisOffset = 1;
double yAxisLimit = -1;
double legendOffset = 0;
String skewErrorMessage = "";
/**
* @param datain Nx2 double
* datain[][0] = year or datenum
* datain[][1] = peak annual flow rate
* @param gg a generalized weighted skew
* @param MSEGbar the error value for the generalized weighted skew
* @param directory the file location of the output file
* @param database the database of the current station
* @param stationId the station ID of the current station
* @param stationName the name of the current station
* @param showLargeFloods if true then the 5 largest floods will be labeled with their corresponding years, if false only the largest flood is labeled.
* @param plotref boolean set to true means reference lines will be plotted, any other value and reference lines will not be plotted. By default, the function will assume set to true if not provided.
* @param plottype boolean, set to true means original plot, other than true creates plot with legend outside of figure and summary table of frequencies outside of figure.
* @return returns a String[][] containing the a summary table of the return periods and flows for the graph and data used
* @throws IOException
*/
public Object[] b17(double[][] datain,
Double gg,
double MSEGbar,
String directory,
String database,
String stationId,
String stationName,
boolean showLargeFloods,
boolean plotref,
boolean plottype) throws IOException, JSONException{
//Note that a number of lines/formulas are commented out after non-commented formulas
//The comment formulas are the original Matlab code lines written by Jeff Burkey
//The non-commented part that comes first is the Java equivalent as written by Tyler Wible
//Check for early return due to data size being too large or too small
//This is for limits on the values in the KNtable provided by USGS (Tyler)
if(datain.length < 10){
ArrayList<String> errorMessage = new ArrayList<>();
errorMessage.add("The flood frequency analysis tool requires a minimum of 10 data points required for the Bulletin 17B statistical analysis. The record of data, based on the specified dates, for this station has " + datain.length + " annual maxima flood observations.");
writeError(errorMessage);
}else if(datain.length > 149){
ArrayList<String> errorMessage = new ArrayList<>();
errorMessage.add("The flood frequency analysis tool requires a maximum of 149 data points for use with the Bulletin 17B statistical analysis. The record of data, based on the specified dates, for this station has " + datain.length + " annual maxima flood observations.");
writeError(errorMessage);
}
//Remove all zero flood years before any calculations are done. (Jeff)
//Matlab code: nonZeroFlood = datain(datain(:,2)>0,:);
ArrayList<Integer> rowIndexList = DoubleArray.findRowsConditional(datain, 1, 0, false);
double[][] nonZeroFlood = DoubleArray.getDataArrayRows(datain, rowIndexList);
//Calculate stats for the B17 station data (Jeff)
//Matlab code: [G N S Xmean] = stationStats(nonZeroFlood);
double[] statsArray = stationStats(nonZeroFlood);
double G = statsArray[0];
double N = statsArray[1];
double S = statsArray[2];
double Xmean = statsArray[3];
double[] skews = {G};
double[][] KNtable = Bulletin17Tables.KNtable;
double[][] ktable = Bulletin17Tables.ktable();
//structure of ktable is:
// Column 1 = Skew
// Column 2 = Probability
// Column 3 = k value
// (Jeff)
//Initialize a number of variables for use outside of the following if statements (Jeff)
int n = datain.length;
int QLcnt = 0;
int QHcnt = 0;
double xl = 0, ql = 0, xh = 0, qh = 0;
double[][] datafilter = new double[0][0];
JSONArray highOutlierList = new JSONArray ();
JSONArray lowOutlierList = new JSONArray ();
if(G >= -.40 && G <= 0.40){
//Estimate outliers thresholds on full record
//Per Appendix 4, one-sided 10-percent KNtable (Jeff)
double KNvalue = findKNvalue(KNtable,N);
xh = Xmean + KNvalue*S; //xh = Xmean + KNtable(KNtable(:,1)==N,2)*S;
qh = Math.pow(10,xh); //qh = 10^xh; //High threshold
xl = Xmean - KNvalue*S; //xl = Xmean - KNtable(KNtable(:,1)==N,2)*S;
ql = Math.pow(10,xl); //ql = 10^xl; //Low threshold
//remove high and low outliers and recompute Station statistics (Jeff)
//Overall Matlab code: datafilter = nonZeroFlood(ql < nonZeroFlood(:,2) & nonZeroFlood(:,2) < qh,[1 2]);
//Broken into two parts for Java (Tyler)
rowIndexList = DoubleArray.findRowsConditional(nonZeroFlood, 1, qh, true);
datafilter = DoubleArray.getDataArrayRows(nonZeroFlood, rowIndexList); //datafilter = nonZeroFlood(nonZeroFlood(:,2) < qh,[1 2]);
rowIndexList = DoubleArray.findRowsConditional(datafilter, 1, ql, false);
datafilter = DoubleArray.getDataArrayRows(datafilter, rowIndexList); //datafilter = datafilter(ql < datafilter(:,2),[1 2]);
//check for low/zero records removed (Jeff)
//Matlab code: QLcnt = sum(datain(:,2) < ql);
for(int i=0; i<datain.length; i++){
if(datain[i][1] < ql){
QLcnt++;
JSONObject low = new JSONObject();
low.put("year", datain[i][0]);
low.put("peak flow", datain[i][1]);
lowOutlierList.put(low);
}
}
//Matlab code: [G N S Xmean] = stationStats(datafilter);
if(QLcnt > 0){
statsArray = stationStats(datafilter);
G = statsArray[0];
N = statsArray[1];
S = statsArray[2];
Xmean = statsArray[3];
}
//check for high records removed (Jeff)
//Matlab code: QHcnt = sum(datain(:,2) > qh);
for(int i=0; i<datain.length; i++){
if(datain[i][1] > qh){
QHcnt++;
JSONObject high = new JSONObject();
high.put("year", datain[i][0]);
high.put("peak flow", datain[i][1]);
highOutlierList.put(high);
}
}
}else if(G > 0.40){
//Test for high outliers, recompute stats, then test for low
//outliers, then recompute stats again (Jeff)
double KNvalue = findKNvalue(KNtable,N);
xh = Xmean + KNvalue*S; //xh = Xmean + KNtable(KNtable(:,1)==N,2)*S;
qh = Math.pow(10,xh); //qh = 10^xh; //High threshold
//Matlab code: datafilter = nonZeroFlood(nonZeroFlood(:,2) < qh,[1 2]);
rowIndexList = DoubleArray.findRowsConditional(nonZeroFlood, 1, qh, true);
datafilter = DoubleArray.getDataArrayRows(nonZeroFlood, rowIndexList);
xl = Xmean - KNvalue*S; //xl = Xmean - KNtable(KNtable(:,1)==N,2)*S;
ql = Math.pow(10,xl); //ql = 10^xl; //Low threshold
//remove outliers and recompute Station statistics (Jeff)
//Matlab code: datafilter = datafilter(ql < datafilter(:,2),[1 2]);
rowIndexList = DoubleArray.findRowsConditional(datafilter, 1, ql, false);
datafilter = DoubleArray.getDataArrayRows(datafilter, rowIndexList);
//check for low/zero records removed (Jeff)
//Matlab code: QLcnt = sum(datain(:,2) < ql);
for(int i=0; i<datain.length; i++){
if(datain[i][1] < ql){
QLcnt++;
JSONObject low = new JSONObject();
low.put("year", datain[i][0]);
low.put("peak flow", datain[i][1]);
lowOutlierList.put(low);
}
}
//Matlab code: [G N S Xmean] = stationStats(datafilter);
if(QLcnt > 0){
statsArray = stationStats(datafilter);
G = statsArray[0];
N = statsArray[1];
S = statsArray[2];
Xmean = statsArray[3];
}
//check for high records removed (Jeff)
//Matlab code: QHcnt = sum(datain(:,2) > qh);
for(int i=0; i<datain.length; i++){
if(datain[i][1] > qh){
QHcnt++;
JSONObject high = new JSONObject();
high.put("year", datain[i][0]);
high.put("peak flow", datain[i][1]);
highOutlierList.put(high);
}
}
}else if(G < -0.40){
//Test for low outliers, recompute stats, then test for high
//outliers, then recompute stats again (Jeff)
double KNvalue = findKNvalue(KNtable,N);
xl = Xmean - KNvalue*S; //xl = Xmean - KNtable(KNtable(:,1)==N,2)*S;
ql = Math.pow(10,xl); //ql = 10^xl; % Low threshold
//remove outliers and recompute Station statistics (Jeff)
//Matlab code: datafilter = nonZeroFlood(ql < nonZeroFlood(:,2),[1 2]);
rowIndexList = DoubleArray.findRowsConditional(nonZeroFlood, 1, ql, false);
datafilter = DoubleArray.getDataArrayRows(nonZeroFlood, rowIndexList);
//check for low/zero records removed (Jeff)
//Matlab code: QLcnt = sum(datain(:,2) < ql);
for(int i=0; i<datain.length; i++){
if(datain[i][1] < ql){
QLcnt++;
JSONObject low = new JSONObject();
low.put("year", datain[i][0]);
low.put("peak flow", datain[i][1]);
lowOutlierList.put(low);
}
}
//Matlab code: [G N S Xmean] = stationStats(datafilter);
if(QLcnt > 0){
statsArray = stationStats(datafilter);
G = statsArray[0];
N = statsArray[1];
S = statsArray[2];
Xmean = statsArray[3];
}
xh = Xmean + KNvalue*S; //xh = Xmean + KNtable(KNtable(:,1)==N,2)*S;
qh = Math.pow(10,xh); //qh = 10^xh; //High threshold
//Matlab code: datafilter = nonZeroFlood(ql < nonZeroFlood(:,2),[1 2]);
rowIndexList = DoubleArray.findRowsConditional(datafilter, 1, qh, true);
datafilter = DoubleArray.getDataArrayRows(datafilter, rowIndexList);
//check for high records removed
//Matlab code: QHcnt = sum(datain(:,2) > qh);
for(int i=0; i<datain.length; i++){
if(datain[i][1] > qh){
QHcnt++;
JSONObject high = new JSONObject();
high.put("year", datain[i][0]);
high.put("peak flow", datain[i][1]);
highOutlierList.put(high);
}
}
}
//Recompute with Historical/high outlier data (Jeff)
//Matlab code: Xz = log10(datain(datain(:,2) > qh,:));
rowIndexList = DoubleArray.findRowsConditional(datain, 1, qh, false);
double[][] highFilter = DoubleArray.getDataArrayRows(datain, rowIndexList);
double[] Xz = DoubleMath.Log10(highFilter, 1);
//Xz = Log10 of historic peaks and high outliers (not low outliers) (Jeff)
//X = Log10 of peaks not including outliers, zero floods, etc.
double[] X = DoubleMath.Log10(datafilter, 1); //X = log10(datafilter);
skews = DoubleArray.appendcolumn(skews, G); //skews = [skews G];
int Z = QHcnt; //Z = Number of historic peaks including high outliers
int H = n; //H = number of years in historic record
int L = QLcnt; //L = Number of low outliers and/or zero flood years
System.out.println("B17: Computing historical correction");
System.out.println("B17: \tH=" + H + "\tL=" + L + "\tZ=" + Z );
//Matlab code: [G M S hp] = historical(H, L, Z, X(:,2), Xz(:,2));
Object[] historicalArray = historical(H, L, Z, X, Xz);
G = (Double) historicalArray[0];
//double M = (Double) historicalArray[1];
S = (Double) historicalArray[2];
//double[][] hp = (double[][]) historicalArray[3]; //hp = historical plot with outliers removed (Jeff)
skews = DoubleArray.appendcolumn(skews, G);// skews = [skews G];
if(!gg.isNaN() && Math.abs(G - gg) > 0.5){
//Notify user more weight should be given to station skew (Jeff)
System.out.println("\nThere is a large discrepency (> 0.5) between the calculated station skew (G = " + G + ") and the generalized skew (GG = " + gg + ").");
System.out.println("More weight should be given to the Station skew.\n");
}
KNtable = null; //clear KNtable
double Pest = 0;
if(QLcnt > 0){
//If low flow outliers or zero flow records occur adjust P with Pest.
//Apply conditional probability adjustment (Appendix 5) (Jeff)
Pest = N/n;
if(Pest < .75){
String[][] errorMessage = {{"Error B170003: Too many years are considered outliers or with zero flow.\nFunction is terminated. Calculated thresholds are low= " + String.valueOf(ql) + " and high = " + String.valueOf(qh)}};
return errorMessage;
}
}else{
Pest = 1;
}
if(N != n){
//notify user of outliers (Jeff)
System.out.println("\nNote: These years were removed as outliers %s");
for(int i=0; i<datain.length; i++){
if(ql > datain[i][1] || datain[i][1] > qh){// outliers = datain(ql > datain(:,2) | datain(:,2) > qh,[1 2]);
System.out.println("Date: " + String.valueOf(datain[i][0]) + "\tFlow: " + String.valueOf(datain[i][1]));
}
}
}
double[][][] freqCurveArray = freqCurve(Xmean, S, G, ktable);//[K floodfreq Gzero] = freqCurve(Xmean, S, G, ktable);
double[][] K = freqCurveArray[0];
double[][] floodfreq = freqCurveArray[1];
double[][] Gzero = freqCurveArray[2];
//Matlab code: floodfreq = [floodfreq floodfreq(:,1)*Pest];
double[] floodfreqPest = new double[floodfreq.length];
for(int i=0; i<floodfreq.length; i++){
floodfreqPest[i] = floodfreq[i][0] * Pest;
}
floodfreq = DoubleArray.appendcolumn_Matrix(floodfreq, floodfreqPest);
//P's outside of adjusted probabilities are extrapolated when using spline method. (Jeff)
double[] floodfreq_column1 = new double[floodfreq.length];
double[] floodfreq_column2 = new double[floodfreq.length];
double[] floodfreq_column3 = new double[floodfreq.length];
for(int i=0; i<floodfreq.length; i++){
floodfreq_column1[i] = floodfreq[i][0];
floodfreq_column2[i] = floodfreq[i][1];
floodfreq_column3[i] = floodfreq[i][2];
}
//Matlab code: adjfreq = [floodfreq(:,1) interp1(floodfreq(:,3), floodfreq(:,2), floodfreq(:,1), 'pchip')];
double[] adjustColum2 = DoubleMath.linearInterpolation(floodfreq_column3, floodfreq_column2, floodfreq_column1);
//double[] adjustColum2 = splineInterpolation(floodfreq_column3, floodfreq_column2, floodfreq_column1);
double[][] adjfreq = DoubleArray.appendcolumn_Matrix(floodfreq_column1, adjustColum2);
double GS = 0, SS = 0, XS = 0;
if(QLcnt > 0){
//Interpolate P back from adjusted P (Jeff)
//Matlab code: Q01 = 10^adjfreq(adjfreq(:,1)==.01,2);
double adjfreqValue = adjfreq[DoubleArray.findRowsEqual(adjfreq, 0, 0.01).get(0)][1];
double Q01 = Math.pow(10, adjfreqValue);
//Matlab code: Q10 = 10^adjfreq(adjfreq(:,1)==.10,2);
adjfreqValue = adjfreq[DoubleArray.findRowsEqual(adjfreq, 0, 0.10).get(0)][1];
double Q10 = Math.pow(10, adjfreqValue);
//Matlab code: Q50 = 10^adjfreq(adjfreq(:,1)==.50,2);
adjfreqValue = adjfreq[DoubleArray.findRowsEqual(adjfreq, 0, 0.50).get(0)][1];
double Q50 = Math.pow(10, adjfreqValue);
//Generate synthetic Log-Pearson statistics the below equation is valid for skews between -2.0 < GS < +2.5 (Jeff)
//Matlab code: GS = -2.5 + 3.12*(log10(Q01/Q10)/log10(Q10/Q50));
GS = -2.5 + 3.12*((Math.log10(Q01/Q10))/(Math.log10(Q10/Q50)));
skews = DoubleArray.appendcolumn(skews, GS); //skews = [skews GS];
//Only using freqCurve to grab K, variables Xmean, S, have no baring for this instance. (Jeff)
//Matlab code: KS = freqCurve(Xmean, S, GS, ktable);
double[][][] freqCurveArray2 = freqCurve(Xmean, S, GS, ktable);
double[][] KS = freqCurveArray2[0];
double K01 = KS[DoubleArray.findRowsEqual(KS, 1, 0.01).get(0)][0]; //K01 = KS(KS(:,2)==.01);
double K50 = KS[DoubleArray.findRowsEqual(KS, 1, 0.50).get(0)][0]; //K50 = KS(KS(:,2)==.50);
SS = (Math.log10(Q01/Q50))/(K01 - K50);
XS = Math.log10(Q50) - K50*SS;
//Per step 4 of Appendix 5
if(GS < -2 || GS > 2.5){
skewErrorMessage = "Skew Error: Synthetic skew exceeds acceptable limits.\nUser should plot data by hand for further evaluation.\n";
}
}else{
GS = G;
XS = Xmean;
SS = S;
}
//Based on 17B Plate I documentation. If generalized skew is not read
//from Plate I, then MSEGbar will need to be calculated. (Jeff)
double A = 0, B = 0;
if(Math.abs(G) <= 0.90){
A = -.33 + 0.08*(Math.abs(G));
}else{
A = -.52 + .30*(Math.abs(G));
}
if(Math.abs(G) <= 1.50){
B = 0.94 - .26*(Math.abs(G));
}else{
B = .55;
}
double MSEG = Math.pow(10, A-B*(Math.log10(H/10))); //10^(A-B*log10(H/10));
//Weighted Skew calculation (Jeff)
double GW = Double.NaN;
if (gg.isNaN()){
// If no generalized skew is provided, assume we are performing these calculations with station skew
GW = G;
} else {
// assume we are performing these calculations with weighted general skew
GW = (MSEGbar*GS + MSEG*gg)/(MSEGbar + MSEG);
}
//Adopted skew (i.e. round skew to nearest tenth (Jeff)
double GD = GW; //GD = round(GW*10)/10;
skews = DoubleArray.appendcolumn(skews, GD);// skews = [skews GD];
//Matlab code: [K finalfreq Gzero] = freqCurve(XS, SS, GD, ktable);
double[][][] freqCurveArray2 = freqCurve(XS, SS, GD, ktable);
K = freqCurveArray2[0];
double[][] finalfreq = freqCurveArray2[1];
Gzero = freqCurveArray2[2];
//Equation 11-1 in Appendix 11 (Jeff)
//Matlab code: pnn = unique(Pntable(:,2));
double[][] PNtable = Bulletin17Tables.PNtable;
double[] pnn = DoubleArray.getUnique(PNtable, 1);
double[] pexp_single = new double[pnn.length];
//lookup and interpolate adjusted probability based on sample size (Jeff)
for(int i=0; i<pnn.length; i++){
//Matlab code: pn2 = Pntable(Pntable(:,2)==pnn(i),:);
ArrayList<Integer> rowIndex2 = DoubleArray.findRowsEqual(PNtable,1,pnn[i]);
double[] pn2_column1 = DoubleArray.getColumn(PNtable, 0, rowIndex2);// pn2(:,1)
double[] pn2_column3 = DoubleArray.getColumn(PNtable, 2, rowIndex2);// pn2(:,3)
//Perform a linear interpolation (Tyler)
//Matlab code: pexp(i) = interp1(pn2(:,1),pn2(:,3),N-1,"pchip");
pexp_single[i] = DoubleMath.linearInterpolation(pn2_column1, pn2_column3, N-1);
//pexp_single[i] = splineInterpolation(pn2_column1, pn2_column3, N-1);
}
double[][] pexp = DoubleArray.appendcolumn_Matrix(pnn, pexp_single); //pexp = [pnn pexp'];
//Matlab code: [c, ia, ib] = intersect(pexp(:,1),finalfreq(:,1),"rows");
ArrayList<Integer> ib = DoubleArray.intersect(DoubleArray.getColumn(pexp, 0), DoubleArray.getColumn(finalfreq, 0));
//Recalculate P using expected probabilities
//'pchip' is used because it better plots linearly in log space w.r.t.
//identified expected probabilities of: .9999, .999, .99, .95, .9, .7,
//.5, .3, .1, .05, .01, .001, and .0001 (Jeff)
//Matlab code: expectedP = [finalfreq(:,1) interp1(pexp(:,2), finalfreq(ib,2), finalfreq(:,1),"pchip")];
double[] column2 = DoubleMath.linearInterpolation(DoubleArray.getColumn(pexp, 1), DoubleArray.getColumn(finalfreq, 1, ib), DoubleArray.getColumn(finalfreq, 0));
//double[] column2 = splinInterpolation(getColumn(pexp, 1), getColumn(finalfreq, 1, ib), getColumn(finalfreq, 0)); old method
double[][] expectedP = DoubleArray.appendcolumn_Matrix(DoubleArray.getColumn(finalfreq, 0), column2);
//Not sure if the below comment is still valid. 1/5/2009.
//USGS PKFQWin seems to not use the adopted, rather the weighted. At
//least with version 5.0 Beta 8 5/6/2005. It also infers a generalized
//skew with more resolution than is available with the published 17B
//documentation. (Jeff)
//In the next line when Find is working again, replace 0.05 with
//cialpha variable. (Jeff)
/////////////////////////////////////// (Jeff)
//Matlab code: Galpha = Gzero(Gzero(:,2)==0.05,:); // (Jeff)
ArrayList<Integer> rowIndex3 = DoubleArray.findRowsEqual(Gzero, 1, 0.05);
double[][] Galpha = DoubleArray.getDataArrayRows(Gzero, rowIndex3);
/////////////////////////////////////// (Jeff)
//Matlab code: a = 1 - (Galpha(:,3).^2)/(2*(H-1));
double a = 1 - Math.pow(Galpha[0][2],2)/(2*(H-1));
double[] b = new double[K.length];
double[] Ku = new double[K.length];
double[] Kl = new double[K.length];
double[] LQu = new double[K.length];
double[] LQl = new double[K.length];
for(int i=0; i<K.length; i++){
b[i] = Math.pow(K[i][0], 2) - Math.pow(Galpha[0][2], 2)/H; //b = K(:,1).^2 - (Galpha(3).^2)/H;
Ku[i] = (K[i][0] + Math.pow(Math.pow(K[i][0], 2) - a*b[i], 0.5))/a; //Ku = (K(:,1) + sqrt(K(:,1).^2 - a.*b))./a;
Kl[i] = (K[i][0] - Math.pow(Math.pow(K[i][0], 2) - a*b[i], 0.5))/a; //Kl = (K(:,1) - sqrt(K(:,1).^2 - a.*b))./a;
LQu[i] = XS + Ku[i]*SS; //LQu = XS + Ku*SS;
LQl[i] = XS + Kl[i]*SS; //LQl = XS + Kl*SS;
}
//expected probabilities breaks down for probabilities less than 0.002,
//replace with NaN's. (Jeff)
//Matlab code: expectedP(expectedP(:,1) < .002,2) = nan;
ArrayList<Integer> rowIndex4 = DoubleArray.findRowsConditional(expectedP, 0, 0.002, true);
expectedP = DoubleArray.replaceRows(expectedP, rowIndex4);
//Create ouput data (Tyler)
//Matlab code: dataout = [1./finalfreq(:,1) finalfreq(:,1) 10.^finalfreq(:,2) 10.^LQu 10.^LQl 10.^expectedP(:,2)];
double[][] dataout = new double[finalfreq.length][6];
for(int i=0; i<dataout.length; i++){
dataout[i][0] = 1/finalfreq[i][0]; //1./finalfreq(:,1)
dataout[i][1] = finalfreq[i][0]; //finalfreq(:,1)
dataout[i][2] = Math.pow(10, finalfreq[i][1]); //10.^finalfreq(:,2)
dataout[i][3] = Math.pow(10,LQu[i]); //10.^LQu
dataout[i][4] = Math.pow(10, LQl[i]); //10.^LQl
dataout[i][5] = Math.pow(10, expectedP[i][1]); //10.^expectedP(:,2)
}
Arrays.sort(dataout, new sort1_smallToLarge()); //dataout = sortrows(dataout,1);
//Matlab code: pp = plotpos(datain(datain(:,2)>0,:),H);
ArrayList<Integer> rowIndex5 = DoubleArray.findRowsConditional(datain, 0, 0, false);
double[][] plotData = DoubleArray.getDataArrayRows(datain, rowIndex5);
double[][] pp = plotpos(plotData, H);
//Graph the data (Tyler)
//Matlab code: pplot(pp(:,2:4),K,dataout(:,2:end),GD,imgfile,gaugeName,plotref,plottype);
pplot(pp, K, dataout, GD, directory, database, stationId, stationName, showLargeFloods, plotref, plottype);
//Assemble a summary of return periods and flow values (Tyler)
String[][] summaryTable = probFreqData(dataout);
Object[] returnArray = {summaryTable, GW, lowOutlierList, highOutlierList};
return returnArray;
}
/**
* stationStats: used to calculate Station Statistics Skew and Standard Deviation
* Computes the statistics for the given data set providing an output of G, N, S, and Xmean
* @param datain data array in (first column = years, 2nd column = flow values for stats)
* @return an Object[] array containing the G (double), N (double), S (double), and Xmean (double) variables in that order
*/
private double[] stationStats(double[][] datain){
//Used to calculate Station Statistics Skew and Standard Deviation
//cialpha = round(1 - cii*100)/100;
double[] data = new double[datain.length]; //data = datain(:,2:end);
for(int i=0; i<datain.length; i++){
data[i] = datain[i][1];
}
//yr = datain(:,1);
double N = data.length;
double Xmean = DoubleMath.Average_Log10(data); //Xmean = mean(log10(data));
double[] X = DoubleMath.Log10(data); //X = log10(data);
double S = calcS(X, N); //S = sqrt((sum(X.^2) - sum(X)^2/N)/(N-1));
//for now check for outliers and notify user
//this function will go away after B17.m is completed.
//outliers = idoutliers(datain, Xmean, S);
double G = calcG(X, N, S); //G = (N^2*sum(X.^3)- 3*N*sum(X)*sum(X.^2) + 2*sum(X)^3)/((N*(N-1)*(N-2))*S^3);
double[] StationStatsArray = {G, N, S, Xmean};
return StationStatsArray;
}
/**
* freqCurve: Retrieve Frequency Curve Coordinates
* Uses ktable.mat to retrieve the skew for a set of probabilities.
* @param Xmean the mean of the dataset (double)
* @param S the S parameter of the data set (double)
* @param G the G parameter of the data set (double)
* @param ktable the ktable containing skew values (double[][])
* @return
* @throws ArgumentOutsideDomainException
*/
private double[][][] freqCurve(double Xmean, double S, double G, double[][] ktable){
//function [K_final floodfreq Gzero] = freqCurve(Xmean, S, G, ktable)
//Retrieve Frequency Curve Coordinates (Jeff)
//Grab list of available probabilites from ktable (Jeff)
double[] P = DoubleArray.getUnique(ktable, 1); //P = unique(ktable(:,2));
//the below loop uses the adopted Skew GD, not the station skew (Jeff)
double[] K = new double[P.length];
for(int i=0; i<P.length; i++){
ArrayList<Integer> p = DoubleArray.findRowsEqual(ktable, 1, P[i]); //p = find(ktable(:,2) == P(i));
double[] g = DoubleArray.getColumn(ktable, 0, p); //g = ktable(p,1);
double[] k = DoubleArray.getColumn(ktable, 2, p); //k = ktable(p,3);
//Replaced lagrange.m
//(http://www.mathworks.com/matlabcentral/fileexchange/14398-a-para
//bolic-lagrangian-interpolating-polynomial-function) for matlab
//spline interpretor. Using interp1 this loop could be reduced to a
//single sequence, but that is for another day. (Jeff)
//
//K(i) = lagrange([g k], G, 4);
//Spline Interpolate using plugin Apache Commons Math Analysis package (Tyler)
//SplineInterpolator splineInterp = new SplineInterpolator();
//PolynomialSplineFunction splineFunction = splineInterp.interpolate(g, k);
//K[i] = splineFunction.value(G);// K[i] = interp1(g,k,G,"spline");
//Linear Interpolation (Tyler)
K[i] = DoubleMath.linearInterpolation(g, k, G);
}
//Matlab code: floodfreq = [P (Xmean + K.*S)];
double[] Xmean_KS = new double[K.length];
for(int i=0; i<K.length; i++){
Xmean_KS[i] = Xmean + K[i]*S;
}
double[][] floodfreq = DoubleArray.appendcolumn_Matrix(P, Xmean_KS);
//add probability to K for looking up values (Jeff)
double[][] K_final = DoubleArray.appendcolumn_Matrix(K, P); //K = [K P];
//Compute confidence limits for estimate
//retrieve zero skew probability (Jeff)
//Matlab code: Gzero = sortrows(ktable(ktable(:,1)==0.0,:),2)
ArrayList<Integer> rowIndex = DoubleArray.findRowsEqual(ktable,0,0);
double[][] Gzero = DoubleArray.getDataArrayRows(ktable,rowIndex);
Arrays.sort(Gzero, new sort2_smallToLarge());
double[][][] freqCurveArray = {K_final, floodfreq, Gzero};
return freqCurveArray;
}
/**
* plotpos: Calculates plotting positions for peak annual flows
* Plotting positions currently calculated are: Gringorten and Weibull.
* Others can be added but the Weibull plotting position is the one used for
* plotting data on the figure.
* @param datain double[][] array of the years and flow values for the input data
* @param n the number of observations in the dataset
* @return returns a double[][] containing a list of:
* the first column the Gringorten plotting position,
* the second column the Weibull plotting position,
* the third column the sorted values of flow occurances provided in datain
* the fourth column the sorted values of the years for the corresponding flow values in the third column
*/
private double[][] plotpos(double[][] datain, int n){
//Matlab code: function [gplot] = plotpos(datain, n)
double N = n;
double[][] dsort = datain.clone();
Arrays.sort(dsort, new sort2_largeToSmall()); //dsort = sortrows(datain,-2);
ArrayList<Integer> rowIndex = DoubleArray.findRowsConditional(dsort, 0, 0, false); //[r c] = find(dsort(:,1));
Iterator<Integer> iterate = rowIndex.iterator();
double[] gp = new double[rowIndex.size()];
double[] gw = new double[rowIndex.size()];
double[][] gplot = new double[rowIndex.size()][4];
int ctr = 0;
while(iterate.hasNext()){
double r = iterate.next() + 1;
gp[ctr] = (r - 0.44)/(N + 0.12); //(r - .44)./(n+.12)
gw[ctr] = r/(N + 1); //gw = (r./(n+1));
//Matlab code: gplot = [gp gw dsort(:,2) dsort(:,1)];
gplot[ctr][0] = gp[ctr];
gplot[ctr][1] = gw[ctr];
gplot[ctr][2] = dsort[ctr][1];
gplot[ctr][3] = dsort[ctr][0];
ctr++;
}
return gplot;
}
/**
* historical: Adjusts the 17B statistics accounting for Historical Data
* High outliers are converted to Historical events-- per Appendix 6 in 17B.
* Also per the 17B method, one can manually supply historical events,
* however, this function does not provide a means for the user to provide
* them. This lack of function generally will not be an issue for the Puget
* Sound region.
* Appendix 6 in 17B in USGS Bulliten 17
* @param H the number of years in historic record
* @param L the number of low outliers and/or zero flood years
* @param Z the number of historic peaks including high outliers
* @param X Log10 of peaks not including outliers, zero floods, etc.
* @param Xz Log10 of historic peaks and high outliers (not low outliers)
* @return summary of historical parameters, returnArray = {Gbar (double), Mbar (double), Sbar (double), hp (double[H][2])};
*/
private Object[] historical(int H, int L, int Z, double[] X, double[] Xz){
//Matlab code: function [Gbar Mbar Sbar hp] = historical(H, L, Z, X, Xz)
//Number of X events
int N = X.length;
//Systematic record weight
double W = (H-Z)/(N+L);
//Historically adjusted mean
double Mbar = (W*DoubleMath.sum(X) + DoubleMath.sum(Xz))/(H-W*L);
//Historically adjusted standard deviation
//Matlab code: Sbar = sqrt((W*sum((X-Mbar).^2) + sum((Xz-Mbar).^2))/(H-W*L-1));
double Sbar = 0;
//Calculate the sum((X-Mbar).^2)
double sumSquareX_Mbar = 0;
for(int i=0; i<X.length; i++){
sumSquareX_Mbar = sumSquareX_Mbar + Math.pow((X[i] - Mbar),2);
}
//Calculate sum((Xz-Mbar).^2))
double sumSquareXz_Mbar = 0;
for(int i=0; i<Xz.length; i++){
sumSquareXz_Mbar = sumSquareXz_Mbar + Math.pow((Xz[i] - Mbar),2);
}
//Calculate Sbar
Sbar = Math.pow(((W*sumSquareX_Mbar + sumSquareXz_Mbar)/(H-W*L-1)),0.5);
//Historically adjusted Skew
//MatLab code: Gbar = ((H-W*L)/((H-W*L-1)*(H-W*L-2)))*((W*sum((X-Mbar).^3) + sum((Xz-Mbar).^3))/Sbar^3);
double Gbar = 0;
//Calculate the sum((X-Mbar).^3)
double sumCubeX_Mbar = 0;
for(int i=0; i<X.length; i++){
sumCubeX_Mbar = sumCubeX_Mbar + Math.pow((X[i] - Mbar),3);
}
//Calculate sum((Xz-Mbar).^3)
double sumCubeXz_Mbar = 0;
for(int i=0; i<Xz.length; i++){
sumCubeXz_Mbar = sumCubeXz_Mbar + Math.pow((Xz[i] - Mbar),3);
}
//Calculate Sbar^3
double SbarCube = Math.pow(Sbar, 3);
//Calculate Gbar
Gbar = ((H-W*L)/((H-W*L-1)*(H-W*L-2)))*((W*sumCubeX_Mbar + sumCubeXz_Mbar)/SbarCube);
//Number of events (Z + N)
//Matlab code: E = 1:(Z+N);
double[] E = new double[Z + N];
for(int i=0; i<E.length; i++){
E[i] = i + 1;
}
//combining historical and systematic peaks
Double[] Xh = DoubleArray.appendcolumn(X, Xz); //Xh = [X; Xz];
//sort from largest to smallest
Arrays.sort(Xh,Collections.reverseOrder()); //Xh = sort(Xh,'descend');
//Define ranks 1 through Z+N
double[] m = new double[E.length];
for(int i=0; i<E.length; i++){
m[i] = E[i];
}
//Historically adjusted ranks starting after historical peaks
double[] mh = new double[E.length];
for(int i=0; i<mh.length; i++){
mh[i] = W*E[i] - ((W-1)*(Z+.50)); //mh = W.*E-((W-1)*(Z+.50));
}
m = DoubleArray.replaceRows(m, mh, Z); //m(Z+1:end) = mh(Z+1:end);
//Matlab code: wp = m./(H+1); //Weibul Plotting Position
double[] wp = new double[m.length];
for(int i=0; i<m.length; i++){
wp[i] = m[i]/(H + 1);
}
//Historical and systematic peaks combined with weibull and gringorten
//plotting positions.
//Matlab code: hp = [wp' 10.^Xh];
double[][] hp = new double[wp.length][2];
for(int i=0; i<wp.length; i++){
hp[i][0] = wp[i];
hp[i][1] = Math.pow(10, Xh[i]);
}
Object[] returnArray = {Gbar, Mbar, Sbar, hp};
return returnArray;
}
/**
* pplot- will create a custom probability plot. The probability spacing
* is based on the weighted skew- it is NOT a normal probability unless
* the weighted skew = 0.0.
*
* On the plot, the last year in the dataset will be plotted with a
* different symbol to highlight. So each year a new plot is created,
* the peak for that last year can be identified on the plot.
*
* //datain = data with plotting position
* //K = K values for final frequency
* //curves = calculated regressions
* //skew = weighted skew used for final frequency curve
* //imgfile = full path and file name to export figure
* //
* //pscale.mat = a vector of user specified tick marks to plot. They do
* // not have to be sorted. One can append any value. If a value
* // specified is beyond the limits of the probabilities calculated in
* // 17B, then the limits will be reset.
* //
* Requirements: none
* Everything used to generate this probability plot can be done
* without toolboxes (e.g. statistics, line, etc.).
*
* Written by
* Jeff Burkey
* King County- Department of Natural Resources and Parks
* Seattle, WA
* email: jeff.burkey@kingcounty.gov
* January 8, 2009
* @param datain data with plotting position (double[][]), ignore the first column
* @param K K values for final frequency (double[][])
* @param curves calculated regressions (double[][]), ignore the first column
* @param skew weighted skew used for final frequency curve (double)
* @param directory the file location of the output file
* @param database the database of the current station
* @param stationId the station ID of the current station
* @param stationName the name of the current station
* @param showLargeFloods if true then the 5 largest floods will be labeled with their corresponding years, if false only the largest flood is labeled.
* @param plotref true if the reference lines are to be plotted, otherwise not plotted (boolean)
* @param plottype boolean, set to true means original plot, other than true creates plot with legend outside of figure and summary table of frequencies outside of figure.
* @param filePath the file location of the pscale file
* @throws IOException
*/
private void pplot( double[][] datain,
double[][] K,
double[][] curves,
double skew,
String directory,
String database,
String stationId,
String stationName,
boolean showLargeFloods,
boolean plotref,
boolean plottype) throws IOException{
//Matlab code: function pplot(datain, K,curves,skew,imgfile,gaugeName,plotref,plottype)
double[][] obs = datain.clone();
Arrays.sort(obs, new sort2_smallToLarge()); //obs = sortrows(datain,2);
//Matlab code: yRnk= [obs interp1(K(:,2),K(:,1), obs(:,1),"pchip")];
double[] column4 = DoubleMath.linearInterpolation(DoubleArray.getColumn(K, 1), DoubleArray.getColumn(K, 0), DoubleArray.getColumn(obs, 1));
//double[] column4 = splineInterpolation(getColumn(K, 1), getColumn(K, 0), getColumn(obs, 1));
double[][] yRnk = new double[column4.length][4];
for(int i=0; i<column4.length; i++){
yRnk[i][0] = obs[i][1];
yRnk[i][1] = obs[i][2];
yRnk[i][2] = obs[i][3];
yRnk[i][3] = column4[i];
}
//build grid for figure (Jeff)
//maybe overkill, but sort and filter out any duplicates. This does let
//the user get lazy defining grid tick marks though. (Jeff)
//Matlab code: sort(unique(pscale_mat),"descend");
double[] pscale = Bulletin17Tables.pscale;
Double[] pscale2 = new Double[pscale.length];
for(int i=0; i<pscale.length; i++){
pscale2[i] = pscale[i];
}
Arrays.sort(pscale2, Collections.reverseOrder());
for(int i=0; i<pscale.length; i++){
pscale[i] = pscale2[i];
}
if(DoubleMath.max(pscale) > DoubleMath.max(DoubleArray.getColumn(K,1))){
double xmin = DoubleMath.max(DoubleArray.getColumn(K,1));
System.out.println("Users specified an exceedance beyond available probabilities.\nLower limit was reset.\n");
ArrayList<Integer> rowIndex = DoubleArray.findRowsConditional(pscale, xmin, true);
pscale = DoubleArray.getRows(pscale, rowIndex); //pscale = pscale(pscale<=xmin);
}
if(DoubleMath.min(pscale) < DoubleMath.min(DoubleArray.getColumn(K,1))){
double xmax = DoubleMath.min(DoubleArray.getColumn(K,1));
System.out.println("Users specified an exceedance beyond available probabilities.\nUpper limit was reset.\n");
ArrayList<Integer> rowIndex = DoubleArray.findRowsConditional(pscale, xmax, false);
pscale = DoubleArray.getRows(pscale, rowIndex); //pscale = pscale(pscale>=xmax);
}
//Matlab code: gridMajor = sort(interp1(K(:,2),K(:,1),pscale,'pchip'));
double[] gridMajor = DoubleMath.linearInterpolation(DoubleArray.getColumn(K, 1), DoubleArray.getColumn(K, 0), pscale);
//double[] gridMajor = splineInterpolation(getColumn(K, 1), getColumn(K, 0), pscale);
Arrays.sort(gridMajor);
double[] xcurve = DoubleArray.getColumn(K, 0);
Arrays.sort(xcurve); //xcurve = sort(K(:,1));
//Start graphing in Java (Tyler)
//Create XYSeries of each line to be graphed (Tyler)
XYSeries B17_xy = new XYSeries("Bulletin 17B");
XYSeries upperCI_95_xy = new XYSeries("95% CI");
XYSeries lowerCI_95_xy = new XYSeries("95% lower CI, Do not show this title");
for(int i=0; i<xcurve.length; i++){
B17_xy.add(xcurve[i], curves[i][2]); //(x,y)
upperCI_95_xy.add(xcurve[i], curves[i][3]); //(x,y)
lowerCI_95_xy.add(xcurve[i], curves[i][4]); //(x,y)
}
XYSeries Weibull_xy = new XYSeries("Data");
for(int i=0; i<yRnk.length; i++){
Weibull_xy.add(yRnk[i][3], yRnk[i][1]); //(x,y)
}
//Determine y axis limits of the graph
//Matlab code: interp1(curves(:,1),curves(:,4),max(pscale),'pchip');
double ymin = DoubleMath.linearInterpolation(DoubleArray.getColumn(curves,1), DoubleArray.getColumn(curves, 4), DoubleMath.max(pscale));
double ymax = DoubleMath.linearInterpolation(DoubleArray.getColumn(curves, 1), DoubleArray.getColumn(curves, 3), DoubleMath.min(pscale));
double[] xLimits = {DoubleMath.min(gridMajor), 1.02*DoubleMath.max(gridMajor)};
//Check if there is any historical data below the calculated ymin, if so reset ymin to the smaller of the two values
//System.out.println("Y max: " + ymax);
//System.out.println("Y min: " + ymin);
if(yRnk[yRnk.length - 1][1] < ymin){
ymin = 0.8*yRnk[yRnk.length - 1][1];
}
//System.out.println("Y min: " + ymin);
double[] yLimits = {ymin, ymax};
//Determine the necessary offsets to create the proper graph
//Note: this resets values of global variables:
yLimits = determineOffsets(yLimits);
//Create XYSeries for each square (=outline) and circle to be plotted on top (Tyler)
//Place a symbol on top of the most recent year (Jeff)
double maxYear = DoubleMath.max(DoubleArray.getColumn(yRnk, 2));
ArrayList<Integer> index = DoubleArray.findRowsEqual(yRnk, 2, maxYear);
XYSeries lastYear_xy = new XYSeries(String.valueOf(Math.round(maxYear)));
lastYear_xy.add(yRnk[index.get(0)][3], yRnk[index.get(0)][1]);
//In addition to plotting the last water year with a square symbol to highlight...
//Draw a drop down reference line from the Final curve representing the return period of most recent year in data (i.e. square symbol data point) (Jeff)
XYSeries lastYear_yLine = new XYSeries("lastYear_yLine (Don't show in Legend)");
double xlast = DoubleMath.linearInterpolation(DoubleArray.getColumn(curves,2), DoubleArray.getColumn(curves, 1), yRnk[index.get(index.size()-1)][1]);
double xLscale = DoubleMath.linearInterpolation(DoubleArray.getColumn(K,1), DoubleArray.getColumn(K,0), xlast);
double xrefMin = DoubleMath.min(gridMajor);
lastYear_yLine.add(xLscale, yLimits[0]);
lastYear_yLine.add(xLscale, yRnk[index.get(0)][1]);
//Higlight the largest flood on record with its year
XYSeries year1_xy = new XYSeries(String.valueOf(Math.round(yRnk[0][2])));
year1_xy.add(yRnk[0][3], yRnk[0][1]);
//Graph each element using Jfreechart
XYPlot plot = new XYPlot();
//Graph line series for plot
plot = graphData(plot, B17_xy, true, false, Color.blue, false, false, true, "", 0, "");
plot = graphData(plot, upperCI_95_xy, true, false, Color.red, true, false, true, "", 1, "");
plot = graphData(plot, lowerCI_95_xy, true, false, Color.red, true, false, false, "", 2, "");
plot = graphData(plot, Weibull_xy, false, true, Color.gray, false, false, true, "circle", 10, "");
//Graph extra years for clarification (Tyler)
plot = graphData(plot, lastYear_xy, false, true, Color.magenta, false, false, false, "outline", 8, String.valueOf(Math.round(maxYear)));
plot = graphData(plot, year1_xy, false, true, Color.blue, false, false, false, "outline", 9, String.valueOf(Math.round(yRnk[0][2])));
//If the user has requested it graph the rest of the 5 largest floods (Tyler)
if(showLargeFloods){
//Add a series for each of the 5 largest floods, except the largest which is already created(Tyler)
XYSeries year2_xy = new XYSeries(String.valueOf(Math.round(yRnk[1][2])));
XYSeries year3_xy = new XYSeries(String.valueOf(Math.round(yRnk[2][2])));
XYSeries year4_xy = new XYSeries(String.valueOf(Math.round(yRnk[3][2])));
XYSeries year5_xy = new XYSeries(String.valueOf(Math.round(yRnk[4][2])));
year2_xy.add(yRnk[1][3], yRnk[1][1]);
year3_xy.add(yRnk[2][3], yRnk[2][1]);
year4_xy.add(yRnk[3][3], yRnk[3][1]);
year5_xy.add(yRnk[4][3], yRnk[4][1]);
plot = graphData(plot, year2_xy, false, true, Color.black, false, false, false, "outline", 4, String.valueOf(Math.round(yRnk[1][2])));
plot = graphData(plot, year3_xy, false, true, Color.black, false, false, false, "outline", 5, String.valueOf(Math.round(yRnk[2][2])));
plot = graphData(plot, year4_xy, false, true, Color.black, false, false, false, "outline", 6, String.valueOf(Math.round(yRnk[3][2])));
plot = graphData(plot, year5_xy, false, true, Color.black, false, false, false, "outline", 7, String.valueOf(Math.round(yRnk[4][2])));
}
//Graph line for clarification of return period of lastYear_xy (Tyler)
plot = graphData(plot, lastYear_yLine, true, false, Color.gray, true, false, false, "", 11, "");
//Add lable to line (Tyler)
String lineLabel1 = DoubleMath.round(1/xlast, 1) + " - year";
double yLocation = yLimits[0]*0.88 + (Math.abs(yLimits[0]*0.88) + Math.abs(yRnk[index.get(0)][1]))/4;
XYTextAnnotation lastYearAnnotation = new XYTextAnnotation(lineLabel1, xLscale, yLocation);
lastYearAnnotation.setFont(Graphing.floodLabelFont);
lastYearAnnotation.setRotationAngle(3*Math.PI/2);
lastYearAnnotation.setBackgroundPaint(Color.white);
plot.addAnnotation(lastYearAnnotation);
//Plot reference lines
// Currently hard coded to create reference lines and display the Flow rate for the
//Upper/Lower 100-yr CI, 100-yr, 25-yr, 10-yr, 5-yr, and 2-yr. If you pick a probability
//not in the output frequency table, you'll have to use an interpolation function to
//extract the value. (Jeff)
if(plotref){
//Create 100year flood line for the 95% confidence interval (Tyler)
XYSeries x100_xy95 = new XYSeries("100 yr CI");
ArrayList<Integer> list = DoubleArray.findRowsEqual(curves, 1, 0.01);
x100_xy95.add(xrefMin, curves[list.get(0)][3]);
x100_xy95.add(xcurve[list.get(0)], curves[list.get(0)][3]);
plot = graphData(plot, x100_xy95, true, false, Color.gray, true, false, false, "", 13, "");
//Add text label for line (Tyler)
String lineLabel = DoubleMath.round(curves[list.get(0)][3], 1) + " - cfs";
double xLocation = xrefMin + (Math.abs(xrefMin) + Math.abs(xcurve[list.get(0)]))/4;
XYTextAnnotation lineAnnotation = new XYTextAnnotation(lineLabel, xLocation, curves[list.get(0)][3]);
lineAnnotation.setFont(Graphing.floodLabelFont);
lineAnnotation.setBackgroundPaint(Color.white);
plot.addAnnotation(lineAnnotation);
//100 year flood line (Tyler)
XYSeries x100_xy = new XYSeries("100 yr");
list = DoubleArray.findRowsEqual(curves, 1, 0.01);
x100_xy.add(xrefMin, curves[list.get(0)][2]);
x100_xy.add(xcurve[list.get(0)], curves[list.get(0)][2]);
plot = graphData(plot, x100_xy, true, false, Color.gray, true, false, false, "", 14, "");
//Add text label for line (Tyler)
lineLabel = DoubleMath.round(curves[list.get(0)][2], 1) + " - cfs";
xLocation = xrefMin + (Math.abs(xrefMin) + Math.abs(xcurve[list.get(0)]))/4;
lineAnnotation = new XYTextAnnotation(lineLabel, xLocation, curves[list.get(0)][2]);
lineAnnotation.setFont(Graphing.floodLabelFont);
lineAnnotation.setBackgroundPaint(Color.white);
plot.addAnnotation(lineAnnotation);
//25 year flood line (Tyler)
XYSeries x25_xy = new XYSeries("25 yr");
list = DoubleArray.findRowsEqual(curves, 1, 0.04);
x25_xy.add(xrefMin, curves[list.get(0)][2]);
x25_xy.add(xcurve[list.get(0)], curves[list.get(0)][2]);
plot = graphData(plot, x25_xy, true, false, Color.gray, true, false, false, "", 15, "");
//Add text label for line
lineLabel = DoubleMath.round(curves[list.get(0)][2], 1) + " - cfs";
lineAnnotation = new XYTextAnnotation(lineLabel, xLocation, curves[list.get(0)][2]);
lineAnnotation.setFont(Graphing.floodLabelFont);
lineAnnotation.setBackgroundPaint(Color.white);
plot.addAnnotation(lineAnnotation);
//10 year flood line (Tyler)
XYSeries x10_xy = new XYSeries("10 yr");
list = DoubleArray.findRowsEqual(curves, 1, 0.1);
x10_xy.add(xrefMin, curves[list.get(0)][2]);
x10_xy.add(xcurve[list.get(0)], curves[list.get(0)][2]);
plot = graphData(plot, x10_xy, true, false, Color.gray, true, false, false, "",16, "");
//Add text label for line (Tyler)
lineLabel = DoubleMath.round(curves[list.get(0)][2], 1) + " - cfs";
lineAnnotation = new XYTextAnnotation(lineLabel, xLocation, curves[list.get(0)][2]);
lineAnnotation.setFont(Graphing.floodLabelFont);
lineAnnotation.setBackgroundPaint(Color.white);
plot.addAnnotation(lineAnnotation);
//5 year flood line (Tyler)
XYSeries x5_xy = new XYSeries("5 yr");
list = DoubleArray.findRowsEqual(curves, 1, 0.2);
x5_xy.add(xrefMin, curves[list.get(0)][2]);
x5_xy.add(xcurve[list.get(0)], curves[list.get(0)][2]);
plot = graphData(plot, x5_xy, true, false, Color.gray, true, false, false, "", 17, "");
//Add text label for line (Tyler)
lineLabel = DoubleMath.round(curves[list.get(0)][2], 1) + " - cfs";
lineAnnotation = new XYTextAnnotation(lineLabel, xLocation, curves[list.get(0)][2]);
lineAnnotation.setFont(Graphing.floodLabelFont);
lineAnnotation.setBackgroundPaint(Color.white);
plot.addAnnotation(lineAnnotation);
//2 year flood line (Tyler)
XYSeries x2_xy = new XYSeries("2 yr");
list = DoubleArray.findRowsEqual(curves, 1, 0.5);
x2_xy.add(xrefMin, curves[list.get(0)][2]);
x2_xy.add(xcurve[list.get(0)], curves[list.get(0)][2]);
plot = graphData(plot, x2_xy, true, false, Color.gray, true, false, false, "", 18, "");
//Add text label for line (Tyler)
lineLabel = DoubleMath.round(curves[list.get(0)][2], 1) + " - cfs";
lineAnnotation = new XYTextAnnotation(lineLabel, xLocation, curves[list.get(0)][2]);
lineAnnotation.setFont(Graphing.floodLabelFont);
lineAnnotation.setBackgroundPaint(Color.white);
plot.addAnnotation(lineAnnotation);
}
//Add custom outline to the graph since the tick marks offset the y axis
plot = addCustomOutline(plot, xLimits, yLimits, 19);
//Add fake tick marks using XYTextAnnotations
plot = addCustomTickMarks(plot, pscale, gridMajor, yLimits, xLimits, 23);
//Create log-scale y axis with scientific notation (Tyler)
LogarithmicAxis yAxis = new LogarithmicAxis("Annual Max. Flow Rate [cfs]");
yAxis.setAxisLinePaint(new Color(233, 233, 233));
yAxis.setTickMarkPaint(Color.black);
yAxis.setRange(yAxisOffset*yLimits[0], yLimits[1]);
yAxis.setAllowNegativesFlag(true);
yAxis.setLog10TickLabelsFlag(true);
plot.setRangeAxis(0,yAxis);
//Create x Axis
ValueAxis domain2 = new NumberAxis("Exceedance Probability (Return Period)");
domain2.setAxisLinePaint(new Color(233, 233, 233));
TickUnits unit1 = new TickUnits();
unit1.add(new NumberTickUnit(-10, NumberFormat.getNumberInstance()));
domain2.setStandardTickUnits(unit1);
domain2.setRange(xLimits[0], xLimits[1]);
plot.setDomainAxis(0, domain2);
//Adjust legend location inside the graph
if(plottype){
//Create the legend inside of the graph
LegendTitle lt = new LegendTitle(plot);
lt.setItemFont(Graphing.masterFont);
lt.setBackgroundPaint(new Color(255, 255, 255, 255));
lt.setFrame(new BlockBorder(Color.black));
lt.setPosition(RectangleEdge.LEFT);
XYTitleAnnotation ta = new XYTitleAnnotation(1, legendOffset, lt, RectangleAnchor.BOTTOM_RIGHT);
plot.addAnnotation(ta);
}
//Set other graph preferences
plot = Graphing.setAxisPreferences(plot);
//Graph plot onto JfreeChart (Tyler)
String graphTitle = "Flood Data for " + database + " Station: " + stationId + "; " + stationName +
"\nWeighted Skew (G=" + DoubleMath.round(skew, 5) + ") Probability Plot";
JFreeChart parentChart = new JFreeChart(graphTitle, Graphing.titleFont, plot, !plottype);
//Adjust the location of the legend outside the graph (Tyler)
if(!plottype){
//Create the legend outside the graph
LegendTitle legend = parentChart.getLegend();
legend.setItemFont(Graphing.masterFont);
legend.setPosition(RectangleEdge.RIGHT);
}
//Save resulting graph (Tyler)
try{
guiFlood_Model model = new guiFlood_Model();
String path = directory + File.separator + model.getGraph();
ChartUtilities.saveChartAsJPEG(new File(path), parentChart, 1280, 800);
System.out.println("JFreeChart created properly at: " + path);
}catch(IOException e) {
System.err.println("A problem occurred while trying to creating the chart.");
}
}
/**
* Creates a summary table of flood flows and their correpsonding return periods
* @param curves the dataout containing the calculated expected and B17 datasets
* to be used for interpolating the flows for each return period
* @return a String[][] containing the data summary table, with headers. The
* first column is the return period, second column is the B17 flood flow, third
* column is the expected flood flow
*/
private String[][] probFreqData(double[][] curves){
//Matlab code: function frqlabel = procFreqData(curves) (Tyler)
//User can edit this to adjust which Return periods are displayed (Jeff)
//Matlab code frqDisp = [200 100 50 40 25 20 10 5 2 1.5 1.25 1.01]'; (Tyler)
double[] frqDisp = {500, 200, 100, 50, 40, 25, 20, 10, 5, 2, 1.5, 1.25, 1.05, 1.01};
String[] frqDispString = {"500", "200", "100", "50", "40", "25", "20", "10", "5.0", "2.0", "1.5", "1.25", "1.05", "1.01"};
//Return periods need to be in terms of probabilities for
//interpolating using the 'curves' variable. (Jeff)
//Matlab code invfrqDisp = 1./frqDisp; (Tyler)
double[] invfrqDisp = new double[frqDisp.length];
for(int i=0; i<invfrqDisp.length; i++){
invfrqDisp[i] = 1/frqDisp[i];
}
//remove any NaN's (Jeff)
//curves(any(isnan(curves),2),:) = [];
double[] yFinal = DoubleMath.linearInterpolation(DoubleArray.getColumn(curves,1), DoubleArray.getColumn(curves,2), invfrqDisp); //yFinal = interp1(curves(:,1),curves(:,2),invfrqDisp,'pchip');
//double[] yExpect = DoubleMath.linearInterpolation(DoubleArray.getColumn(curves,1), DoubleArray.getColumn(curves,5), invfrqDisp); //yExpect = interp1(curves(:,1),curves(:,5),invfrqDisp,'pchip');
double[] yFinal95L = DoubleMath.linearInterpolation(DoubleArray.getColumn(curves,1), DoubleArray.getColumn(curves,4), invfrqDisp); //Added for extra information (Tyler)
double[] yFinal95U = DoubleMath.linearInterpolation(DoubleArray.getColumn(curves,1), DoubleArray.getColumn(curves,3), invfrqDisp); //Added for extra information (Tyler)
//Round all values to one decimal place for summary table
yFinal = DoubleMath.roundColumn(yFinal, 2);
//yExpect = DoubleMath.roundColumn(yExpect, 2);
yFinal95L = DoubleMath.roundColumn(yFinal95L, 2);
yFinal95U = DoubleMath.roundColumn(yFinal95U, 2);
//Matlab code: mm = num2str([frqDisp yFinal yExpect],'%5.2f %7.0f %7.0f');
//below is very specific to mm format above (Jeff)
//Matlab code: frqlabel = [...
//'Return 17B Expectd';...
//'(year) (cfs) (cfs) ';...
//'----------------------';...
//mm;];
String[][] summaryTable_new = new String[invfrqDisp.length + 2][4];
//Row 1 (Tyler)
summaryTable_new[0][0] = "Return Period";
//summaryTable_new[0][1] = "Expected";
summaryTable_new[0][1] = "Lower 95% CI for B17";
summaryTable_new[0][2] = "B17";
summaryTable_new[0][3] = "Upper 95% CI for B17";
//Row 2 (Tyler)
summaryTable_new[1][0] = "(year)";
//summaryTable_new[1][1] = "(cfs)";
summaryTable_new[1][1] = "(cfs)";
summaryTable_new[1][2] = "(cfs)";
summaryTable_new[1][3] = "(cfs)";
for(int i=0; i<invfrqDisp.length; i++){
summaryTable_new[i+2][0] = frqDispString[i]; // return period
//summaryTable_new[i+2][1] = String.valueOf(yExpect[i]); // Expected value
summaryTable_new[i+2][1] = String.valueOf(yFinal95L[i]); // Lower 95% CI from B17 value
summaryTable_new[i+2][2] = String.valueOf(yFinal[i]); // B17 value
summaryTable_new[i+2][3] = String.valueOf(yFinal95U[i]); // Upper 95% CI from B17 value
}
return summaryTable_new;
}
/**
* Calculates the S value given X and N for b17
* @param X a double[] array
* @param N the length of the above array
* @return the S value for the corresponding X and N values
*/
private double calcS(double[] X, double N){
//Matlab code: S = sqrt((sum(X.^2) - sum(X)^2/N)/(N-1));
double S = 0;
double sum = 0, sumSquare = 0;
//Calculate the S value
for(int i=0; i<X.length; i++){//Loop rows of X
sum = sum + X[i];
sumSquare = sumSquare + Math.pow(X[i],2);
}
S = Math.pow((sumSquare - (Math.pow(sum, 2)/N))/(N-1), 0.5);
return S;
}
/**
* Calculates the G parameter based on X, N and S
* @param X double[] array containing X
* @param N double containing the value of X.length
* @param S double containing S
* @return double containing G
*/
private double calcG(double[] X, double N, double S){
//Matlab code: G = (N^2*sum(X.^3)- 3*N*sum(X)*sum(X.^2) + 2*sum(X)^3)/((N*(N-1)*(N-2))*S^3);
double G = 0;
double sumX = 0, sumX2 = 0, sumX3 = 0;
for(int i=0; i<X.length; i++){
sumX = sumX + X[i];//sum the elements of the current row of X
sumX2 = sumX2 + Math.pow(X[i],2);//sum the square of the elements of the current row of X
sumX3 = sumX3 + Math.pow(X[i],3);//sum the cube of the elements of the current row of X
}
G = (((Math.pow(N,2))*sumX3) - (3*N*sumX*sumX2) + ((2*(Math.pow(sumX, 3))) / ((N*(N-1)*(N-2))*(Math.pow(S, 3)))));
return G;
}
/**
* Searches the KNtable's first column for a matching value of N and returns the matching row's second column
* @param KNtable the double[][] array containing the KNtable
* @param N the current data size
* @return the value in the second column of KNtable correspoding to the row with KNtable's first column equaling N
*/
private double findKNvalue(double[][]KNtable, double N){
//Matlab code: KNtable(KNtable(:,1)==N,2)*S;
double KNvalue = 0;
for(int i=0; i<KNtable.length; i++){
if(KNtable[i][0] == N){
KNvalue = KNtable[i][1];
}
}
return KNvalue;
}
/**
* Graphs the XYSeries on to the provided plot with the below properties for rendering, color, series index
* @param plot The XYPlot to graph the XYSeries on
* @param currentLine the XYSeries containing the x,y points of the current line series
* @param lineTrue a boolean, true if the series is to have lines connecting the points, otherwise no line
* @param shapeTrue a boolean, true if the series is to have shapes on each of the points, otherwise no points.
* Note that if lineTrue is "true" then shapeTrue needs to be "false", and vise-versa
* @param seriesColor The desired color of the series (Java.Color)
* @param dashedLine if true then the XYSeries will have a dashed line instead of a solid one, if false it will be the normal solid line
* @param shortDash if true a dashed line like the major grid lines will be drawn, otherwise a long dashed line will be used
* @param visibleInLegend if true then the XYSeries' name will be visible in the legend, if false this XYSeries' name will not be visible in the legend
* @param shapeType the type of shape, if any, to be rendered, either circle or rectangle
* @param graphIndex the series index for the current line, this should be incrementally increased
* each time this function is call so that the new dataset does not replace the old one
* @param extraText an extra text string used to fill the label for outline objects and...
* @return the XYPlot with the added XYSeries
*/
private XYPlot graphData( XYPlot plot,
XYSeries currentLine,
boolean lineTrue,
boolean shapeTrue,
Color seriesColor,
boolean dashedLine,
boolean shortDash,
boolean visibleInLegend,
String shapeType,
int graphIndex,
String extraText){
//Create Line Data and renderer
XYDataset currentDataset = new XYSeriesCollection(currentLine);
//Check if this series should have a dashed line or solid line and change the renderer accordingly
XYItemRenderer currentRenderer = new XYLineAndShapeRenderer(lineTrue, shapeTrue);
currentRenderer.setSeriesPaint(0, seriesColor);
currentRenderer.setSeriesVisibleInLegend(0, visibleInLegend);
//Check if this series should have a dashed line, if so change the renderer
if(dashedLine){
if(shortDash){
//Change the renderer's stroke to a short dashed line
currentRenderer.setSeriesStroke(0,
new BasicStroke(
1.0f, BasicStroke.CAP_ROUND, BasicStroke.JOIN_ROUND,
1.0f, new float[] {1.0f, 2.0f}, 0.0f
));
}else{
//Change the renderer's stroke to a long dashed line
currentRenderer.setSeriesStroke(0,
new BasicStroke(
1.0f, BasicStroke.CAP_ROUND, BasicStroke.JOIN_ROUND,
1.0f, new float[] {6.0f, 6.0f}, 0.0f
));
}
}
//Check if this series should have lines or shapes and change the renderer accordingly
if(!lineTrue && shapeTrue){
if(shapeType.equalsIgnoreCase("circle")){
currentRenderer.setSeriesShape(0, new Ellipse2D.Double(-3.0, 6.0, 6.0, 6.0));
}else if(shapeType.equalsIgnoreCase("outline")){
//This was commented out because this function is used elsewhere and the icon is unessesary
currentRenderer.setSeriesShape(0, new Ellipse2D.Double(-3.0, 3.0, 6.0, 6.0));
double x = currentLine.getX(0).doubleValue();
double y = currentLine.getY(0).doubleValue();
//Add an outline around the point
CircleDrawer cd = new CircleDrawer(seriesColor, new BasicStroke(1.0f), null);
XYAnnotation bestBid = new XYDrawableAnnotation(x, y, 11, 11, cd);
plot.addAnnotation(bestBid);
//Add an arrow pointing to the point with the year of that flood
final XYPointerAnnotation pointer = new XYPointerAnnotation(extraText, x, y, Math.PI / 4.0);
pointer.setBaseRadius(35.0);
pointer.setTipRadius(10.0);
pointer.setFont(Graphing.floodLabelFont);
pointer.setPaint(Color.black);
pointer.setTextAnchor(TextAnchor.HALF_ASCENT_LEFT);
plot.addAnnotation(pointer);
}else{
currentRenderer.setSeriesShape(0, new Rectangle2D.Double(-3.0, 3.0, 6.0, 6.0));
}
}
//Put the LDC line data, renderer, and axis into plot
plot.setDataset(graphIndex, currentDataset);
plot.setRenderer(graphIndex, currentRenderer);
//Put the line on the first Domain and first Range
plot.mapDatasetToDomainAxis(0, 0);
plot.mapDatasetToRangeAxis(0, 0);
return plot;
}
/**
* Subfunction used to set the offset limits for the y axis allowing the fake x-tick marks to be created in the proper place
* @param yLimits a double[] containing the following yLimits[0] = minimum of y data, yLimits[1] = maximum of y data
* @return the input yLimits altered slightly for visual purposes
*/
private double[] determineOffsets(double[] yLimits){
//Determine the log-limits of the current data set
//Note: should expect between 10^0 and 10^7 values for minLimit and maxLimit
double minLimit = Math.log10(yLimits[0]);
double maxLimit = Math.log10(yLimits[1]);
minLimit = Math.floor(minLimit);
maxLimit = Math.ceil(maxLimit);
//System.out.println("Log Limits: " + minLimit + "\t" + maxLimit);
//Determine based on current data range what the offsets should be
if(Double.isInfinite(minLimit)){
double[] tempArray = {1, 4, 4, 4, 3.8, 4.1, 3.8, 4.3, 4, 4, 4.3, 4.3, 4.6, 5.3, 5.3, 5.7};
tickLabelOffset = tempArray;
yAxisOffset = 20;
yAxisLimit = -25;
legendOffset = 0.01;
tickMarkOffset = 1.25;
//System.out.println("Resetting yLimit:" + yAxisLimit);
yLimits[0] = -1;
}else if(minLimit <= 0 && maxLimit == 3){//USGS 09085400
double[] tempArray = {1, 0.56, 0.56, 0.59, 0.59, 0.56, 0.6, 0.56, 0.57, 0.57, 0.56, 0.57, 0.54, 0.5, 0.5, 0.49};
tickLabelOffset = tempArray;
yAxisOffset = 0.25;
yAxisLimit = 0;
legendOffset = 0.02;
}else if(minLimit <= 0 && maxLimit == 4){//USGS 06754500
double[] tempArray = {1, 0.41, 0.41, 0.44, 0.44, 0.41, 0.47, 0.4, 0.43, 0.43, 0.4, 0.4, 0.38, 0.33, 0.33, 0.31};
tickLabelOffset = tempArray;
yAxisOffset = 0.1;
yAxisLimit = 0;
legendOffset = 0.02;
}else if(minLimit <= 0 && maxLimit == 5){//USGS 06755960
double[] tempArray = {1, 0.4, 0.4, 0.43, 0.43, 0.4, 0.46, 0.4, 0.43, 0.43, 0.4, 0.4, 0.37, 0.32, 0.32, 0.3};
tickLabelOffset = tempArray;
yAxisOffset = 0.1;
yAxisLimit = 0;
legendOffset = 0.02;
}else if(minLimit <= 0 && maxLimit >= 6){//USGS 07103703, 07105900
double[] tempArray = {1, 0.13, 0.13, 0.15, 0.15, 0.13, 0.18, 0.13, 0.15, 0.15, 0.13, 0.13, 0.11, 0.075, 0.075, 0.06};
tickLabelOffset = tempArray;
yAxisOffset = 0.005;
yAxisLimit = 0;
legendOffset = 0.00002;
}else if(minLimit == 1 && maxLimit == 3){//USGS 06748200
double[] tempArray = {1, 0.63, 0.63, 0.65, 0.65, 0.63, 0.67, 0.63, 0.65, 0.65, 0.63, 0.63, 0.61, 0.56, 0.56, 0.54};
tickLabelOffset = tempArray;
yAxisOffset = 0.3;
yAxisLimit = 0;
legendOffset = 0.03;
}else if(minLimit == 1 && maxLimit == 4){//USGS 16557000
double[] tempArray = {1, 0.59, 0.59, 0.61, 0.61, 0.59, 0.63, 0.59, 0.61, 0.61, 0.58, 0.58, 0.56, 0.515, 0.515, 0.5};
tickLabelOffset = tempArray;
yAxisOffset = 0.3;
yAxisLimit = 0;
legendOffset = 0.02;
}else if(minLimit == 1 && maxLimit == 5){//USGS 02221000
double[] tempArray = {1, 0.5, 0.5, 0.53, 0.53, 0.5, 0.55, 0.5, 0.52, 0.52, 0.5, 0.5, 0.47, 0.43, 0.43, 0.41};
tickLabelOffset = tempArray;
yAxisOffset = 0.18;
yAxisLimit = 0;
legendOffset = 0.02;
}else if(minLimit == 2 && maxLimit == 4){//USGS 06747500
double[] tempArray = {1, 0.65, 0.65, 0.67, 0.67, 0.65, 0.69, 0.65, 0.67, 0.67, 0.65, 0.65, 0.63, 0.59, 0.59, 0.57};
tickLabelOffset = tempArray;
yAxisOffset = 0.3;
yAxisLimit = 0;
legendOffset = 0.1;
}else if(minLimit == 2 && maxLimit == 5){//USGS 11501000
double[] tempArray = {1, 0.48, 0.48, 0.5, 0.5, 0.48, 0.51, 0.48, 0.49, 0.49, 0.46, 0.46, 0.43, 0.39, 0.39, 0.37};
tickLabelOffset = tempArray;
yAxisOffset = 0.15;
yAxisLimit = 0;
legendOffset = 0.04;
}else if(minLimit == 2 && maxLimit == 6){//USGS 06756995
double[] tempArray = {1, 0.43, 0.43, 0.46, 0.46, 0.43, 0.48, 0.43, 0.46, 0.46, 0.43, 0.43, 0.4, 0.35, 0.35, 0.32};
tickLabelOffset = tempArray;
yAxisOffset = 0.1;
yAxisLimit = 0;
legendOffset = 0.02;
}else if(minLimit == 2 && maxLimit == 7){//USGS 06759500
double[] tempArray = {1, 0.31, 0.31, 0.35, 0.35, 0.31, 0.39, 0.31, 0.35, 0.35, 0.31, 0.31, 0.29, 0.24, 0.24, 0.22};
tickLabelOffset = tempArray;
yAxisOffset = 0.05;
yAxisLimit = 0;
legendOffset = 0.02;
}else if(minLimit == 3 && maxLimit == 5){//USGS 09072500
double[] tempArray = {1, 0.65, 0.65, 0.67, 0.67, 0.65, 0.68, 0.65, 0.67, 0.67, 0.65, 0.65, 0.63, 0.59, 0.59, 0.57};
tickLabelOffset = tempArray;
yAxisOffset = 0.3;
yAxisLimit = 0;
legendOffset = 0.1;
}else if(minLimit == 3 && maxLimit == 6){//USGS 06796000
double[] tempArray = {1, 0.57, 0.57, 0.59, 0.59, 0.57, 0.61, 0.57, 0.59, 0.59, 0.57, 0.57, 0.55, 0.5, 0.5, 0.48};
tickLabelOffset = tempArray;
yAxisOffset = 0.23;
yAxisLimit = 0;
legendOffset = 0.06;
}else if(minLimit == 3 && maxLimit == 7){//USGS 06337500
double[] tempArray = {1, 0.45, 0.45, 0.49, 0.49, 0.45, 0.5, 0.45, 0.49, 0.49, 0.45, 0.45, 0.42, 0.38, 0.38, 0.35};
tickLabelOffset = tempArray;
yAxisOffset = 0.1;
yAxisLimit = 0;
legendOffset = 0.02;
}else if(minLimit == 4 && maxLimit == 6){//USGS 06807000
double[] tempArray = {1, 0.67, 0.67, 0.69, 0.69, 0.67, 0.7, 0.67, 0.69, 0.69, 0.67, 0.67, 0.65, 0.61, 0.61, 0.59};
tickLabelOffset = tempArray;
yAxisOffset = 0.35;
yAxisLimit = 0;
legendOffset = 0.1;
}else if(minLimit < 2 && maxLimit > 6){//USGS 06758200
double[] tempArray = {1, 0.06, 0.06, 0.075, 0.075, 0.06, 0.08, 0.06, 0.075, 0.075, 0.06, 0.06, 0.05, 0.03, 0.03, 0.025};
tickLabelOffset = tempArray;
yAxisOffset = 0.001;
yAxisLimit = 0;
legendOffset = 0.0000001;
}else if(maxLimit > 4){//USGS 07374000
double[] tempArray = {1, 0.82, 0.82, 0.83, 0.83, 0.82, 0.84, 0.82, 0.83, 0.83, 0.82, 0.82, 0.81, 0.785, 0.785, 0.77};
tickLabelOffset = tempArray;
yAxisOffset = 0.6;
yAxisLimit = 0;
legendOffset = 0.2;
tickMarkOffset = 0.99;
}
return yLimits;
}
/**
* Create fau tick labels for the cdf values and their corresponding return periods.
* @param plot the XYPlot to add the tick marks to
* @param pscale double[] the cdf values to be displayed on the bottom and 1/pscale = the return period to be shown at the top
* @param xLocationArray double[], the corresponding array of x locations of the tick marks for for the above pscale labels
* @param yLimits a double[] = {the lower y limit, the upper y limit} of the graph
* @param xLimits a double[] = {the lower x limit, the upper x limit} of the graph
* @param graphIndex the starting graph series index for the major tick markers
* @return the XYPlot with the tick marks added
*/
private XYPlot addCustomTickMarks(XYPlot plot, double[] pscale, double[] xLocationArray, double[] yLimits, double[] xLimits, int graphIndex){
//Loop through and create a label for each element in pscale except the first one at location x = xLocationArray[i]
double xOffset = 1, decimalAmount = 0;
for(int i=1; i<xLocationArray.length; i++){
graphIndex++;
if(i <= 7){
xOffset = 1;
decimalAmount = 2;//100;
}else{
xOffset = 1;
decimalAmount = -1;//0.1;
}
//Create text of return period value and CDF value for x-axis label
String label = DoubleMath.round(pscale[i], 3) + "\t (" + DoubleMath.round(1/pscale[i], decimalAmount) + "- yr)";
//Add grid line to mimic grid lines for probability plots
XYSeries majorTickLine = new XYSeries(label);
majorTickLine.add(xLocationArray[i], yLimits[0]);
majorTickLine.add(xLocationArray[i], yLimits[1]);
plot = graphData(plot, majorTickLine, true, false, Color.lightGray, true, true, false, "", graphIndex, "");
graphIndex++;
//Add text annotation to the grid line as a "tick label"
XYTextAnnotation bottomAnnotation = new XYTextAnnotation(label, xOffset*xLocationArray[i], tickLabelOffset[i]*yLimits[0]);
bottomAnnotation.setFont(Graphing.floodLabelFont);
bottomAnnotation.setRotationAngle(3*Math.PI/2);
plot.addAnnotation(bottomAnnotation);
//Add tick line to the tick label
XYSeries tickMarks = new XYSeries(label + "tick Mark");
tickMarks.add(xLocationArray[i], yLimits[0]);
tickMarks.add(xLocationArray[i], tickMarkOffset*yLimits[0]);
plot = graphData(plot, tickMarks, true, false, Color.black, false, false, false, "", graphIndex, "");
}
return plot;
}
/**
* Adds a custom outline at the specified x-y limits of the graph
* @param plot the XYPlot to which the outline will be added
* @param xLimits a double array containing the x limits of the outline to be created (x-min, x-max)
* @param yLimits a double array containing the y limits of the outline to be created (y-min, y-max)
* @param graphIndex the current series index for the graph, this will be increased by 4 for the lines created to bound the axis
* @return the XYPlot provided with the outline added at the x-min/x-max and y-min/y-max locations
*/
private XYPlot addCustomOutline(XYPlot plot, double[] xLimits, double[] yLimits, int graphIndex){
//To complete custom tick markes add an IntervalMarker as a background the same color as the graph background
final IntervalMarker label = new IntervalMarker(yAxisLimit, yLimits[0]);
label.setLabel(" ");
label.setLabelAnchor(RectangleAnchor.BOTTOM);
label.setLabelTextAnchor(TextAnchor.BASELINE_CENTER);
label.setPaint(new Color(233, 233, 233));
label.setOutlinePaint(Color.black);
label.setOutlineStroke(null);
plot.addRangeMarker(label, Layer.BACKGROUND);
plot.setOutlineStroke(null);
//Create an outline around the desired portion of the graph
//Add left boundary
XYSeries leftBoundary = new XYSeries("LeftBound");
leftBoundary.add(xLimits[0], yLimits[0]);
leftBoundary.add(xLimits[0], yLimits[1]);
plot = graphData(plot, leftBoundary, true, false, Color.black, false, false, false, "", graphIndex, "");
graphIndex++;
//Add right boundary
XYSeries rightBoundary = new XYSeries("RightBound");
rightBoundary.add(0.9999*xLimits[1], yLimits[0]);
rightBoundary.add(0.9999*xLimits[1], yLimits[1]);
plot = graphData(plot, rightBoundary, true, false, Color.black, false, false, false, "", graphIndex, "");
graphIndex++;
//Add top boundary
XYSeries topBoundary = new XYSeries("TopBound");
topBoundary.add(xLimits[0], yLimits[1]);
topBoundary.add(xLimits[1], yLimits[1]);
plot = graphData(plot, topBoundary, true, false, Color.black, false, false, false, "", graphIndex, "");
graphIndex++;
//Add bottom boundary
XYSeries bottomBoundary = new XYSeries("BottomBound");
bottomBoundary.add(xLimits[0], yLimits[0]);
bottomBoundary.add(xLimits[1], yLimits[0]);
plot = graphData(plot, bottomBoundary, true, false, Color.black, false, false, false, "", graphIndex, "");
graphIndex++;
return plot;
}
/**
* Writes out the error message, if any, for finding the file and then exits the program
* @param error string array to be written as each line of an error message
* @throws IOException
*/
private void writeError(ArrayList<String> error) throws IOException{
//Output data to text file
String errorContents = error.get(0);
for(int i=1; i<error.size(); i++){
errorContents = errorContents + "\n" + error.get(i);
}
throw new IOException("Error encountered. Please see the following message for details: \n" + errorContents);
}
// //The following comments are from from Jeff from the Matlab version of this program:
//
// //Flood Flow Frequency Analysis - Bulletin #17B USGS
// //b17.m - This function estimates Flood Frequencies
// //Using U.S. Water Resources
// //Council publication Flood Flow Frequencies, Bulletin #17B (1976, revised
// //1981,1982). This method uses a Log Pearson Type-3 distribution for estimating
// //quantiles. See url: http://water.usgs.gov/osw/bulletin17b/bulletin_17B.html
// //for further documentation.
// //
// //NaN need to be removed in the dataset. If any years have missing data, it
// //will still assume to include that year as part of the sample size-- as
// //stipulated in 17B guidelines. An exmaple MAT file is provided for the
// //user to test. Further down in these comments is some sample script to
// //pre-process the examples.mat file provided.
// //
// //There are only a couple of loops in this function and subfunctions, most
// //of this is vectorized for speed.
// //
// //A nice enhancement to this function is a plot of the analyses. It is
// //*plotted in probability space-- SKEWED probability space!* Because data
// //may not be normally distributed, plotting in skewed space maintains a
// //straight line for the final frequency curve. Again, no need of any
// //special toolboxes to create this figure. All self contained in this
// //function.
// //
// //Modifications:
// //
// //* January 28, 2009 - Tweaked the figure created by cleaning up the
// // legend, reference lines can be turned on/off and are hard coded to
// // display: Upper/Lower 100yr CIs, 100yr, 25yr, 10yr, 5yr, and 2yr. Also
// // changed a few of the default tick placements used to make the grid.
// //
// //* March 13, 2009 - Removed the lower confidence reference line on the
// // figure. Also now will project the last water year in the data set to
// // the final frequency curve and plot a drop down reference line
// // annotating what the interpreted return period would be.
// //
// //* April 20, 2009 - Modified plotting function pplot to either create a
// // plot like the original of this funciton, or to create a plot and
// // located the legend outside the figure and create a table summarizing a
// // subset of flow frequencies likely to be of interest. Added a bit more
// // in the legend, cleaned up the title block. Made the comments in this
// // m-file Report friendly as well.
// //
// //* December 2, 2009 - People were mistakenly downloading a different
// // version of lagrange. My version
// // (http://www.mathworks.com/matlabcentral/fileexchange/14398-a-parabolic
// //-lagrangian-interpolating-polynomial-function) should have been the
// //one to use. I've removed the lagrange interpolator and used the
// //intrinsic function in Matlab (interp1 'spline'). This should have no
// //significant effect of change, but if it does it only will make it more
// //accurate. Thanks to Shan Zou for sleuthing out this problem. Also,
// //if you're not familiar with matlab reporting, you'll notice somewhat
// // extraneous symbols (e.g. #, *, etc). Those are symbols used in the
// // reporting interpretor for organizing the text.
// //
// //Outputs of this function include:
// //
// //# estimates of a final frequency (based on a weighted skew),
// //# confidence intervals (95) for the final frequency,
// //# expected frequency curve based on an adjusted conditional
// // probability,
// //# observed data with computed plotting positions using Gringorten and
// // Weibull techniques (no toolbox required),
// //# various Skews,
// //# mean of log10(Q),
// //# standard deviation of log10(Q),
// //# and the coup de grace,
// //# a probability plot that does not require a toolbox to create, but
// //also plots the probability space using the computed weighted skew
// //and not just the normal probability.
// //
// //*Note:*
// //This added feature yields a straight line plot for the final
// // frequency curve even if the data are not normally distributed.
// //
// //*Important*
// //
// //The one important aspect not included in this funtion is the assumed
// //generalized skew (which is variable throughout the country), which can be
// //obtained from Plate 1 in the bulletin. Using the USGS program PKFQWin,
// //this generalized skew is automatically estimated with given lat/long
// //coordinates. For this function, the user must specify a generalized
// //skew, if no generalized skew is provided, 0.0 is assumed.
// //Even though this function computes probabilities, skews, etc., no
// //toolboxes are required. All necessary tables are provided as additional
// //MAT files supporting this function. These tables are created from the
// //published USGS 17B manual, and not taken from any Matlab toolboxes, so there are no
// //conflicts or copyright violations.
// //
// //Other files required to support this function are:
// //
// //# KNtable.mat - using normal distribution, a table of 10-percent
// // significance level of K values per N sample size.
// //# ktable.mat - Appendix 3 table Pearson distributions
// //# PNtable.mat - Table 11-1 in Appendix 11. Table of probabilities
// // adjusted for sample size.
// //# pscale.mat - table used to define tick/grid intervals when creating a
// // probability plot of the results. Can be modified by user if other than
// // the default values.
// //# examples.mat - dataset presented in the 17B publication.
// //
// //*NOTE: lagrange is no longer needed for this function.
// //Parabolic interpolation of Pearson Distribution is dependant on
// //function: lagrange.m (written by Jeff Burkey 3/23/2007). Can be
// //downloaded from Mathworks user community.
// //
// //Syntax
// // [dataout skews pp XS SS hp] = b17(datain, gg, imgfile, gaugeName, plotref, plottype);
// //
// //Where
// // datain = Nx2 double
// // datain(:,1) = year or datenum
// // datain(:,2) = peak annual flow rate
// // gg = a generalized weighted skew
// // imgfile = full path and file name to export the figure
// // example: imgfile = 'd:\temp\figure.png'
// // gaugeName = string used to populate title of figure
// // example: 'USGS 12108500'
// // plotref = integer set to 1 means reference lines will be plotted,
// // any other value and reference lines will not be plotted. By default,
// // the function will assume set to 1 if not provided.
// // plottype = integer, set to 1 means original plot, other than 1 creates
// // plot with legend outside of figure and summary table of frequencies
// // outside of figure.
// //
// // *note: confidence intervals are hard coded to 0.95 at present
// //
// //Output is in the form of a N x 6 double, and all flows are in Log10
// //
// //* dataout(:,1) = Return Period (in years)
// //* dataout(:,2) = Probability
// //* dataout(:,3) = Final Frequency curve
// //* dataout(:,4) = upper CI frequency curve
// //* dataout(:,5) = lower CI frequency curve
// //* dataout(:,6) = Expected Probability (N-1)
// //
// //* skews(1) = Station Skew
// //* skews(2) = Station Skew outliers, zero removed
// //* skews(3) = Synthetic Skew
// //* skews(4) = Weight adjusted skew
// //
// //* pp(:,1) = gringorten plotting position
// //* pp(:,2) = weibull plotting position
// //* pp(:,3) = observations
// //
// //* XS = mean of Log10(Flow rates)
// //* SS = standard deviation Log10(flow rates)
// //
// //* hp = historically adjusted weibull plotting position where historical
// //peaks our high outliers are ranked in integer values (e.g. 1,2,3,etc.)
// //and systematic observations are adjusted.
// //
// //Subfunctions contained in this function are:
// // stationStats- Used to calculate Station Statistics Skew and Standard
// // Deviation
// // gw = generalized regional skew coefficient
// // for Puget Sound = 0.0
// // freqCurve- compute curve statistics
// // plotpos - calculate plotting positions of data using Gringorten and
// // Weibull.
// // historical- adjust statistics using historical peaks, also recomputes
// // plotting positions when including historical.
// // pplot - creates a customized probability plot.
// // procFreqData - creates text table with summary of frequencies
// //
// //Some example script to process the examples.mat file also provided
// //
// // ex1 = examples(~isnan(examples(:,2)),1:2);
// // ex2 = [examples(~isnan(examples(:,3)),1) examples(~isnan(examples(:,3)),3)];
// // ex3 = [examples(~isnan(examples(:,4)),1) examples(~isnan(examples(:,4)),4)];
// // ex4 = [examples(~isnan(examples(:,5)),1) examples(~isnan(examples(:,5)),5)];
// //
// // written by
// // Jeff Burkey
// // King County- Department of Natural Resources and Parks
// // Seattle, WA
// // email: jeff.burkey@kingcounty.gov
// // January 6, 2009
// //
// // [dataout skews pp XS SS hp] = b17(ex1, .590537, 'c:\temp\test.png', 'Demo Station',1,2);
}