ingo@420: package de.intevation.gnv.utils;
ingo@420: 
ingo@420: import com.vividsolutions.jts.geom.Coordinate;
ingo@420: import com.vividsolutions.jts.geom.Point;
ingo@420: import com.vividsolutions.jts.geom.LineString;
ingo@420: import com.vividsolutions.jts.io.ParseException;
ingo@420: import com.vividsolutions.jts.io.WKTReader;
ingo@420: 
ingo@420: import de.intevation.gnv.geobackend.base.DefaultResult;
ingo@420: import de.intevation.gnv.geobackend.base.DefaultResultDescriptor;
ingo@420: import de.intevation.gnv.geobackend.base.Result;
ingo@420: import de.intevation.gnv.geobackend.base.ResultDescriptor;
ingo@420: import de.intevation.gnv.geobackend.base.query.QueryExecutor;
ingo@420: import de.intevation.gnv.geobackend.base.query.QueryExecutorFactory;
ingo@420: import de.intevation.gnv.geobackend.base.query.exception.QueryException;
ingo@420: import de.intevation.gnv.math.Interpolation2D;
ingo@420: import de.intevation.gnv.math.LinearFunction;
ingo@420: import de.intevation.gnv.math.LinearMetrics;
ingo@420: import de.intevation.gnv.math.Point2d;
ingo@420: import de.intevation.gnv.state.InputData;
ingo@420: 
ingo@420: import java.util.Arrays;
ingo@420: import java.util.ArrayList;
ingo@420: import java.util.Collection;
ingo@420: import java.util.List;
ingo@420: import java.util.Map;
ingo@420: 
ingo@420: import org.apache.commons.math.optimization.OptimizationException;
ingo@420: import org.apache.commons.math.optimization.fitting.CurveFitter;
ingo@420: import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
ingo@420: import org.apache.commons.math.FunctionEvaluationException;
ingo@420: 
ingo@420: import org.apache.log4j.Logger;
ingo@420: 
ingo@420: public abstract class WKTUtils {
ingo@420: 
ingo@420:     private static Logger log = Logger.getLogger(WKTUtils.class);
ingo@420: 
ingo@420:     private static final String [] DIFF_COLUMS = {
ingo@420:         "GROUP1",
ingo@420:         "GROUP2",
ingo@420:         "GROUP3"
ingo@420:     };
ingo@420: 
ingo@420:     private static final String [] COLUMN_BLACKLIST = {
ingo@420:         "MEDIAN.MESHPOINT.JPOSITION",
ingo@420:         "MEDIAN.MESHPOINT.IPOSITION"
ingo@420:     };
ingo@420: 
ingo@420:     public static final double NAUTICAL_MILE = 1852.216d;
ingo@420:     public static final double KILOMETER     = 1000d;
ingo@420: 
ingo@420:     public static final double EPSILON = 1e-5d;
ingo@420:     public static final int    INTERPOLATION_STEPS =
ingo@420:         Integer.getInteger("interpolation.steps", 500).intValue();
ingo@420: 
ingo@420: 
ingo@420:     public static final class SectionHandler
ingo@420:     implements                Interpolation2D.Consumer
ingo@420:     {
ingo@420:         private ArrayList<Point2d> points;
ingo@420:         private List<Coordinate>   path;
ingo@420:         private Collection<Result> output;
ingo@420:         private Result             prototyp;
ingo@420:         private ResultDescriptor   descriptor;
ingo@420:         private boolean            lastWasSuccess;
ingo@420: 
ingo@420:         public SectionHandler() {
ingo@420:         }
ingo@420: 
ingo@420:         public SectionHandler(
ingo@420:             List<Coordinate>   path,
ingo@420:             Collection<Result> output,
ingo@420:             ResultDescriptor   descriptor
ingo@420:         ) {
ingo@420:             this.path       = path;
ingo@420:             this.output     = output;
ingo@420:             this.descriptor = descriptor;
ingo@420:             points          = new ArrayList<Point2d>();
ingo@420:             lastWasSuccess  = true;
ingo@420:         }
ingo@420: 
ingo@420:         public void finish() {
ingo@420:             if (!points.isEmpty()) {
ingo@420:                 double distance = toKM(
ingo@420:                     DistanceCalculator.calculateDistance(path));
ingo@420: 
ingo@420:                 if (distance > EPSILON) {
ingo@420: 
ingo@420:                     Interpolation2D.interpolate(
ingo@420:                         path,
ingo@420:                         points,
ingo@420:                         0d,
ingo@420:                         distance,
ingo@420:                         INTERPOLATION_STEPS,
ingo@420:                         LinearMetrics.INSTANCE,
ingo@420:                         this);
ingo@420:                 }
ingo@420: 
ingo@420:                 points.clear();
ingo@420:             }
ingo@420:             lastWasSuccess = true;
ingo@420:         }
ingo@420: 
ingo@420:         public void setPrototyp(Result prototyp) {
ingo@420:             this.prototyp = prototyp;
ingo@420:         }
ingo@420: 
ingo@420:         public void handle(Result result) {
ingo@420:             Coordinate coordinate = 
ingo@420:                 toCoordinate(result.getString("SHAPE"));
ingo@420:             double value = result.getDouble("YORDINATE");
ingo@420:             int iPos     = result.getInteger("MEDIAN.MESHPOINT.JPOSITION");
ingo@420:             int jPos     = result.getInteger("MEDIAN.MESHPOINT.JPOSITION");
ingo@420:             Point2d p = new Point2d(
ingo@420:                 coordinate.x,
ingo@420:                 coordinate.y,
ingo@420:                 value,
ingo@420:                 iPos, jPos);
ingo@420:             points.add(p);
ingo@420:         }
ingo@420: 
ingo@420:         public void interpolated(Coordinate coordinate, boolean success) {
ingo@420:             DefaultResult result = new DefaultResult(descriptor);
ingo@420:             ResultDescriptor pd = prototyp.getResultDescriptor();
ingo@420: 
ingo@420:             int pcolums = pd.getColumnCount();
ingo@420:             for (int i = 0, j = 0; i < pcolums; ++i) {
ingo@420:                 String colname = pd.getColumnName(i);
ingo@420:                 if (blacklisted(colname)) {
ingo@420:                     continue;
ingo@420:                 }
ingo@420:                 if (colname.equals("SHAPE")) {
ingo@420:                     result.addColumnValue(j, OutputUtils.toWKT(coordinate));
ingo@420:                 }
ingo@420:                 else if (colname.equals("YORDINATE")) {
ingo@420:                     if (success) {
ingo@420:                         result.addColumnValue(j, Double.valueOf(coordinate.z));
ingo@420:                     }
ingo@420:                     else if (lastWasSuccess) {
ingo@420:                         // only insert null if last was valid.
ingo@420:                         // This prevents flooding the result set with nulls
ingo@420:                         // if interpolating over a large gap.
ingo@420:                         result.addColumnValue(j, null);
ingo@420:                     }
ingo@420:                 }
ingo@420:                 else {
ingo@420:                     result.addColumnValue(j, prototyp.getObject(i));
ingo@420:                 }
ingo@420:                 ++j;
ingo@420:             }
ingo@420:             output.add(result);
ingo@420:             lastWasSuccess = success;
ingo@420:         }
ingo@420:     }
ingo@420: 
ingo@420: 
ingo@420:     private static final boolean blacklisted(String column) {
ingo@420:         for (int i = 0; i < COLUMN_BLACKLIST.length; ++i) {
ingo@420:             if (COLUMN_BLACKLIST.equals(column)) {
ingo@420:                 return true;
ingo@420:             }
ingo@420:         }
ingo@420:         return false;
ingo@420:     }
ingo@420: 
ingo@420:     private static boolean different(Result a, Result b, int [] indices) {
ingo@420:         for (int i = 0; i < indices.length; ++i) {
ingo@420:             String oa = a.getString(indices[i]);
ingo@420:             String ob = b.getString(indices[i]);
ingo@420: 
ingo@420:             if (oa == null && ob == null)  {
ingo@420:                 continue;
ingo@420:             }
ingo@420: 
ingo@420:             if (oa == null || ob == null) {
ingo@420:                 return true;
ingo@420:             }
ingo@420: 
ingo@420:             if (!oa.equals(ob)) {
ingo@420:                 if (log.isDebugEnabled()) {
ingo@420:                     log.debug("+++++++++++++++ differs ++++++++++++++");
ingo@420:                     log.debug("   " + oa + " != " + ob);
ingo@420:                 }
ingo@420:                 return true;
ingo@420:             }
ingo@420:         }
ingo@420:         return false;
ingo@420:     }
ingo@420: 
ingo@420:     public static Collection<Result> process(
ingo@420:         List<Coordinate>   path,
ingo@420:         Collection<Result> input
ingo@420:     ) {
ingo@420:         log.debug("------  number of points before processing: " + input.size());
ingo@420:         ArrayList<Result> output = new ArrayList<Result>();
ingo@420: 
ingo@420:         Result last = null;
ingo@420: 
ingo@420:         int [] diffColums = null;
ingo@420: 
ingo@420:         SectionHandler sectionHandler = null;
ingo@420: 
ingo@420:         for (Result result: input) {
ingo@420: 
ingo@420:             if (sectionHandler == null) {
ingo@420: 
ingo@420:                 ResultDescriptor rd = result.getResultDescriptor();
ingo@420:                 diffColums = rd.getColumnIndices(DIFF_COLUMS);
ingo@420:                 int columns = rd.getColumnCount();
ingo@420: 
ingo@420:                 DefaultResultDescriptor resultDescriptor =
ingo@420:                     new DefaultResultDescriptor();
ingo@420: 
ingo@420:                 for (int j = 0; j < columns; ++j) {
ingo@420:                     String columnName = rd.getColumnName(j);
ingo@420:                     if (!blacklisted(columnName)) {
ingo@420:                         resultDescriptor.addColumn(
ingo@420:                             columnName,
ingo@420:                             rd.getColumnClassName(j));
ingo@420:                     }
ingo@420:                 }
ingo@420: 
ingo@420:                 sectionHandler = new SectionHandler(
ingo@420:                     path,
ingo@420:                     output,
ingo@420:                     resultDescriptor);
ingo@420: 
ingo@420:                 sectionHandler.setPrototyp(result);
ingo@420:             }
ingo@420: 
ingo@420:             if (last != null && different(last, result, diffColums)) {
ingo@420:                 sectionHandler.finish();
ingo@420:                 sectionHandler.setPrototyp(result);
ingo@420:             }
ingo@420: 
ingo@420:             sectionHandler.handle(result);
ingo@420: 
ingo@420:             last = result;
ingo@420:         }
ingo@420: 
ingo@420:         if (sectionHandler != null) {
ingo@420:             sectionHandler.finish();
ingo@420:         }
ingo@420: 
ingo@420:         log.debug("------  number of points after processing: " + output.size());
ingo@420: 
ingo@420:         return output;
ingo@420:     }
ingo@420: 
ingo@420: 
ingo@420:     public static Coordinate toCoordinate(String shape) {
ingo@420:         try {
ingo@420:             return ((Point)(new WKTReader().read(shape))).getCoordinate();
ingo@420:         }
ingo@420:         catch (ParseException pe) {
ingo@420:             log.error(pe);
ingo@420:         }
ingo@420:         return null;
ingo@420:     }
ingo@420: 
ingo@420: 
ingo@420:     public static final double toKM(double distance) {
ingo@420:         return (distance * NAUTICAL_MILE) / KILOMETER;
ingo@420:     }
ingo@420: 
ingo@420: 
ingo@420:     public static String toWKT(Coordinate coordinate) {
ingo@420:         StringBuilder sb = new StringBuilder("POINT(");
ingo@420:         sb.append(coordinate.x)
ingo@420:           .append(' ')
ingo@420:           .append(coordinate.y)
ingo@420:           .append(')');
ingo@420:         return sb.toString();
ingo@420:     }
ingo@420: 
ingo@420: 
ingo@420:     public static Collection<Result> worldCoordinatesToIndex(
ingo@420:         Collection<Result>     result,
ingo@420:         Map<String, InputData> inputData,
ingo@420:         String                 ijkQueryID,
ingo@420:         String                 queryID,
ingo@420:         String[]               filterValues
ingo@420:     ) throws ParseException, QueryException
ingo@420:     {
ingo@420:         // 1. IJK Anfragen für Stuetzpunkte im Netz ausführen.
ingo@420:         LineString ls = (LineString)new WKTReader().read(
ingo@420:             inputData.get("mesh_linestring").getValue());
ingo@420: 
ingo@420:         Coordinate[] coords = ls.getCoordinates();
ingo@420: 
ingo@420:         List<java.awt.Point> points = new ArrayList<java.awt.Point>(coords.length);
ingo@420: 
ingo@420:         String meshid = inputData.get("meshid").getValue();
ingo@420:         QueryExecutor queryExecutor = QueryExecutorFactory
ingo@420:                                         .getInstance()
ingo@420:                                         .getQueryExecutor();
ingo@420: 
ingo@420:         ArrayList missingPoints = new ArrayList();
ingo@420: 
ingo@420:         String additionWhere = "FEATUREID=FEATUREID";
ingo@420: 
ingo@420:         for (int i = 0; i < coords.length; i++) {
ingo@420: 
ingo@420:             String wkt = toWKT(coords[i]);
ingo@420: 
ingo@420:             result = queryExecutor.executeQuery(ijkQueryID,
ingo@420:                                                new String[]{meshid,wkt});
ingo@420:             if (!result.isEmpty()){
ingo@420:                 Result resultValue = result.iterator().next();
ingo@420:                 int iPos = resultValue.getInteger(1);
ingo@420:                 int jPos = resultValue.getInteger(0);
ingo@420:                 log.debug("Found Pos "+iPos+"/"+jPos +" for "+wkt);
ingo@420:                 points.add(i, new java.awt.Point(iPos,jPos));
ingo@420:             }else{
ingo@420:                 log.debug("No i/j Pos found for "+wkt);
ingo@420:                 missingPoints.add(new Object [] { Integer.valueOf(i), coords[i] });
ingo@420:                 points.add(i, null);
ingo@420:                 // Special Case no i,j found for Coordinate
ingo@420:             }
ingo@420:         }
ingo@420: 
ingo@420:         if (missingPoints.size() == coords.length) {
ingo@420:             log.debug("cannot create index buffer");
ingo@420:         }
ingo@420:         else { // generate index filter
ingo@420:             boolean remainsMissingPoints = !missingPoints.isEmpty();
ingo@420: 
ingo@420:             if (remainsMissingPoints) {
ingo@420:                 // try to guess the missing (i, j)
ingo@420:                 CurveFitter iFitter = new CurveFitter(new GaussNewtonOptimizer(true));
ingo@420:                 CurveFitter jFitter = new CurveFitter(new GaussNewtonOptimizer(true));
ingo@420: 
ingo@420:                 for (int i = 0, N = points.size(); i < N; ++i) {
ingo@420:                     java.awt.Point p = (java.awt.Point)points.get(i);
ingo@420:                     if (p != null) {
ingo@420:                         Coordinate coord = coords[i];
ingo@420:                         iFitter.addObservedPoint(coord.x, p.x);
ingo@420:                         jFitter.addObservedPoint(coord.y, p.y);
ingo@420:                     }
ingo@420:                 }
ingo@420:                 try {
ingo@420:                     // XXX: Assumption: (i, j) are created by componentwise linear function.
ingo@420:                     // This is surely not correct because (x, y) are in a ellipsoid projection.
ingo@420:                     // TODO: use ellipsoid functions and fit with Levenberg Marquardt.
ingo@420:                     double [] iParams = iFitter.fit(
ingo@420:                         LinearFunction.INSTANCE, new double [] { 1d, 1d });
ingo@420: 
ingo@420:                     double [] jParams = jFitter.fit(
ingo@420:                         LinearFunction.INSTANCE, new double [] { 1d, 1d });
ingo@420: 
ingo@420:                     for (int i = missingPoints.size()-1; i >= 0; --i) {
ingo@420:                         Object [] a = (Object [])missingPoints.get(i);
ingo@420:                         Coordinate coord = (Coordinate)a[1];
ingo@420:                         int pi = (int)Math.round(iParams[0]*coord.x + iParams[1]);
ingo@420:                         int pj = (int)Math.round(jParams[0]*coord.y + jParams[1]);
ingo@420:                         points.set(
ingo@420:                             ((Integer)a[0]).intValue(),
ingo@420:                             new java.awt.Point(pi, pj));
ingo@420:                     }
ingo@420: 
ingo@420:                     remainsMissingPoints = false; // we filled the gaps
ingo@420:                 }
ingo@420:                 catch (FunctionEvaluationException fee) {
ingo@420:                     log.error(fee);
ingo@420:                 }
ingo@420:                 catch (OptimizationException oe) {
ingo@420:                     log.error(oe);
ingo@420:                 }
ingo@420:             }
ingo@420: 
ingo@420:             if (!remainsMissingPoints) {
ingo@420:                 // TODO: Make Tablenames and Columns Configurable
ingo@420:                 IndexBuffer ib = new IndexBuffer(
ingo@420:                     points, 
ingo@420:                     "MEDIAN.MESHPOINT.IPOSITION", 
ingo@420:                     "MEDIAN.MESHPOINT.JPOSITION" );
ingo@420:                 additionWhere = ib.toWhereClause();
ingo@420:                 log.debug("Additional Where Clause = "+additionWhere);
ingo@420:                 // 2. Aus diesen Stuetzpunkten den Resultset generieren.
ingo@420:             }
ingo@420:         } // if generate index filter
ingo@420: 
ingo@420:         String[] addedFilterValues = new String[filterValues.length+1];
ingo@420:         System.arraycopy(filterValues, 0, addedFilterValues, 0, filterValues.length);
ingo@420:         addedFilterValues[filterValues.length] = additionWhere;
ingo@420: 
ingo@420:         result = process(
ingo@420:             Arrays.asList(coords),
ingo@420:             queryExecutor.executeQuery(
ingo@420:                 queryID,
ingo@420:                 addedFilterValues));
ingo@420: 
ingo@420:         return result;
ingo@420:     }
ingo@420: }