view gnv-artifacts/src/main/java/de/intevation/gnv/utils/WKTUtils.java @ 433:828df3ddb758

Added interpolation capabilities along z axis to XYColumns. gnv-artifacts/trunk@481 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 23 Dec 2009 12:20:25 +0000
parents 2402173a1490
children eb2ac62e853a
line wrap: on
line source
package de.intevation.gnv.utils;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Point;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.io.ParseException;
import com.vividsolutions.jts.io.WKTReader;

import de.intevation.gnv.geobackend.base.Result;
import de.intevation.gnv.geobackend.base.query.QueryExecutor;
import de.intevation.gnv.geobackend.base.query.QueryExecutorFactory;
import de.intevation.gnv.geobackend.base.query.exception.QueryException;
import de.intevation.gnv.math.LinearFunction;
import de.intevation.gnv.state.InputData;

import java.util.ArrayList;
import java.util.Collection;
import java.util.List;
import java.util.Map;

import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.fitting.CurveFitter;
import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
import org.apache.commons.math.FunctionEvaluationException;

import org.apache.log4j.Logger;

public abstract class WKTUtils {

    private static Logger log = Logger.getLogger(WKTUtils.class);

    public static final double NAUTICAL_MILE = 1852.216d;
    public static final double KILOMETER     = 1000d;

    public static boolean different(Result a, Result b, int [] indices) {
        for (int i = 0; i < indices.length; ++i) {
            String oa = a.getString(indices[i]);
            String ob = b.getString(indices[i]);

            if (oa == null && ob == null)  {
                continue;
            }

            if (oa == null || ob == null) {
                return true;
            }

            if (!oa.equals(ob)) {
                if (log.isDebugEnabled()) {
                    log.debug("+++++++++++++++ differs ++++++++++++++");
                    log.debug("   " + oa + " != " + ob);
                }
                return true;
            }
        }
        return false;
    }


    public static Coordinate toCoordinate(String shape) {
        try {
            return ((Point)(new WKTReader().read(shape))).getCoordinate();
        }
        catch (ParseException pe) {
            log.error(pe);
        }
        return null;
    }


    public static final double toKM(double distance) {
        return (distance * NAUTICAL_MILE) / KILOMETER;
    }


    public static String toWKT(Coordinate coordinate) {
        StringBuilder sb = new StringBuilder("POINT(");
        sb.append(coordinate.x)
          .append(' ')
          .append(coordinate.y)
          .append(')');
        return sb.toString();
    }


    public static String worldCoordinatesToIndex(
        Collection<Result>     result,
        Map<String, InputData> inputData,
        String                 ijkQueryID
    ) throws ParseException, QueryException
    {
        // 1. IJK Anfragen für Stuetzpunkte im Netz ausführen.
        LineString ls = (LineString)new WKTReader().read(
            inputData.get("mesh_linestring").getValue());

        Coordinate[] coords = ls.getCoordinates();

        List<java.awt.Point> points = new ArrayList<java.awt.Point>(coords.length);

        String meshid = inputData.get("meshid").getValue();
        QueryExecutor queryExecutor = QueryExecutorFactory
                                        .getInstance()
                                        .getQueryExecutor();

        ArrayList missingPoints = new ArrayList();

        String additionWhere = "FEATUREID=FEATUREID";

        for (int i = 0; i < coords.length; i++) {

            String wkt = toWKT(coords[i]);

            result = queryExecutor.executeQuery(ijkQueryID,
                                               new String[]{meshid,wkt});
            if (!result.isEmpty()){
                Result resultValue = result.iterator().next();
                int iPos = resultValue.getInteger(1);
                int jPos = resultValue.getInteger(0);
                log.debug("Found Pos "+iPos+"/"+jPos +" for "+wkt);
                points.add(i, new java.awt.Point(iPos,jPos));
            }else{
                log.debug("No i/j Pos found for "+wkt);
                missingPoints.add(new Object [] { Integer.valueOf(i), coords[i] });
                points.add(i, null);
                // Special Case no i,j found for Coordinate
            }
        }

        if (missingPoints.size() == coords.length) {
            log.debug("cannot create index buffer");
        }
        else { // generate index filter
            boolean remainsMissingPoints = !missingPoints.isEmpty();

            if (remainsMissingPoints) {
                // try to guess the missing (i, j)
                CurveFitter iFitter = new CurveFitter(new GaussNewtonOptimizer(true));
                CurveFitter jFitter = new CurveFitter(new GaussNewtonOptimizer(true));

                for (int i = 0, N = points.size(); i < N; ++i) {
                    java.awt.Point p = (java.awt.Point)points.get(i);
                    if (p != null) {
                        Coordinate coord = coords[i];
                        iFitter.addObservedPoint(coord.x, p.x);
                        jFitter.addObservedPoint(coord.y, p.y);
                    }
                }
                try {
                    // XXX: Assumption: (i, j) are created by componentwise linear function.
                    // This is surely not correct because (x, y) are in a ellipsoid projection.
                    // TODO: use ellipsoid functions and fit with Levenberg Marquardt.
                    double [] iParams = iFitter.fit(
                        LinearFunction.INSTANCE, new double [] { 1d, 1d });

                    double [] jParams = jFitter.fit(
                        LinearFunction.INSTANCE, new double [] { 1d, 1d });

                    for (int i = missingPoints.size()-1; i >= 0; --i) {
                        Object [] a = (Object [])missingPoints.get(i);
                        Coordinate coord = (Coordinate)a[1];
                        int pi = (int)Math.round(iParams[0]*coord.x + iParams[1]);
                        int pj = (int)Math.round(jParams[0]*coord.y + jParams[1]);
                        points.set(
                            ((Integer)a[0]).intValue(),
                            new java.awt.Point(pi, pj));
                    }

                    remainsMissingPoints = false; // we filled the gaps
                }
                catch (FunctionEvaluationException fee) {
                    log.error(fee);
                }
                catch (OptimizationException oe) {
                    log.error(oe);
                }
            }

            if (!remainsMissingPoints) {
                // TODO: Make Tablenames and Columns Configurable
                IndexBuffer ib = new IndexBuffer(
                    points, 
                    "MEDIAN.MESHPOINT.IPOSITION", 
                    "MEDIAN.MESHPOINT.JPOSITION" );
                additionWhere = ib.toWhereClause();
                log.debug("Additional Where Clause = "+additionWhere);
                // 2. Aus diesen Stuetzpunkten den Resultset generieren.
            }
        } // if generate index filter

        return additionWhere;
    }
}

http://dive4elements.wald.intevation.org