/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package usda.weru.gis;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Envelope;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryCollection;
import com.vividsolutions.jts.geom.GeometryCollectionIterator;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.Point;


import java.text.DecimalFormat;
import java.text.NumberFormat;
import java.text.ParseException;
import java.util.List;
import javax.measure.Measurable;
import javax.measure.Measure;
import javax.measure.quantity.Length;
import javax.measure.unit.NonSI;
import javax.measure.unit.SI;
import org.apache.log4j.Logger;
import org.jscience.geography.coordinates.LatLong;
import org.opengis.feature.Feature;
import org.opengis.feature.GeometryAttribute;
import org.opengis.geometry.BoundingBox;
import org.opengis.geometry.DirectPosition;



import org.opengis.geometry.coordinate.Position;


/**
 *
 * @author Joseph Levin <joelevin@weru.ksu.edu>
 */
public class GISUtil {
    private static final Logger LOGGER = Logger.getLogger(GISData.class);

    public static Coordinate toCoordinate(LatLong latlong){
        return new Coordinate(latlong.longitudeValue(NonSI.DEGREE_ANGLE), latlong.latitudeValue(NonSI.DEGREE_ANGLE));
    }

    public static LatLong toLatLong(Position position){
        DirectPosition dp = position.getDirectPosition();
        return LatLong.valueOf(dp.getOrdinate(0), dp.getOrdinate(1), NonSI.DEGREE_ANGLE);
    }

    public static LatLong toLatLong(Coordinate coord){
        return LatLong.valueOf(coord.y, coord.x, NonSI.DEGREE_ANGLE);
    }

    public static LatLong parse(String text){
        String[] parts = text.split(";");
        assert parts.length == 2;
        NumberFormat format = new DecimalFormat("+0.####;-0.####");
        try{
            double lat = format.parse(parts[0]).doubleValue();
            double lon = format.parse(parts[1]).doubleValue();
            return LatLong.valueOf(lat, lon, NonSI.DEGREE_ANGLE);
        }
        catch(ParseException pe){
            return null;
        }
    }

    public static String format(LatLong latlong){
        double lat = latlong.latitudeValue(NonSI.DEGREE_ANGLE);
        double lon = latlong.longitudeValue(NonSI.DEGREE_ANGLE);

        NumberFormat format = new DecimalFormat("0.00#");
        StringBuffer buffer = new StringBuffer();
        buffer.append(lat >= 0 ? "+" : "-");
        buffer.append(format.format(Math.abs(lat)));
        buffer.append(";");
        buffer.append(lon >= 0 ? "+" : "-");
        buffer.append(format.format(Math.abs(lon)));


        return buffer.toString();
    }
    

    /*
     * Calculates distance between pairs of coordinations.  Uses the heversine
     * formula.
     */
    public static Measurable<Length> distanceBetweenCoordinates(LatLong latlong1, LatLong latlong2) {
        
        if(latlong1 == null || latlong2 == null){
            return null;
        }

        //The mean radius http://wn.wikipedia.org/wiki/Earth_Radius in KILOMETERS
        double meanEarthRadius = 6371;

        double lat1 = latlong1.latitudeValue(SI.RADIAN);
        double long1 = latlong1.longitudeValue(SI.RADIAN);
        
        double lat2 = latlong2.latitudeValue(SI.RADIAN);
        double long2 = latlong2.longitudeValue(SI.RADIAN);


        double deltaLat = lat2 - lat1;
        double deltaLon = long2 - long1;

        double a = Math.pow(Math.sin(deltaLat / 2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(deltaLon / 2), 2);
        double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
        double d = meanEarthRadius * c;
        return Measure.valueOf(d, SI.KILOMETER);      
    }

    public static <T> LatLong representativeLatLong(T t){
        if(t == null){
            //silly to ask for the centroid of nothing.  Just plain silly.
            return null;
        }
        //use a cast to find the correct lookup
        GISLookup<T> lookup = GISData.getInstance().getLookup((Class<T>) t.getClass());

        if(lookup == null){
            LOGGER.warn("Unable to calculate centroid. No registered GISLookups for " + t.getClass().getCanonicalName());
            return null;
        }

        List<Feature> features = lookup.lookup(t);
        if(features == null || features.isEmpty()){
            LOGGER.warn("Unable to calculate centroid. No gis features found for " + t.toString());
            return null;
        }

        if(features.size() > 1){
            LOGGER.debug("Found " + features.size() + " features for object.  Determining dominant feature.");
        }

        Feature dominantFeature = null;        
        double dominantArea = 0;

        for(Feature feature : features){
            
            Geometry featureGeometry = geometry(feature);            
                        
            if(featureGeometry == null){
                continue;
            }
                        
            double area = areaOrLength(featureGeometry);

            if(dominantFeature == null){
                dominantFeature = feature;
                dominantArea = area;
            }
            else if (area > dominantArea){
                dominantFeature = feature;
                dominantArea = area;
            }
        }
       
       
        if(dominantFeature != null){

            Point point = representativePoint(geometry(dominantFeature));
            
            if(point != null){
                return toLatLong(new Coordinate(point.getX(), point.getY()));
            }

        }

        //couldn't calculate a centroid
        return null;
    }
    
    public static Geometry geometry(Feature feature){
        if(feature == null){
            return null;
        }
        
        GeometryAttribute attribute = feature.getDefaultGeometryProperty();
        if(attribute == null){
            return null;
        }
        
        Object value = attribute.getValue();
        if(value instanceof Geometry){
            return (Geometry) value;
        }
        
        return null;
                
    }

    /**
     * 
     * @param geometry
     * @return
     * @throws NullPointerException if the geometry is null
     */
    public static Point representativePoint(Geometry geometry) throws NullPointerException{
        if(geometry == null){
            throw new NullPointerException("Can not find a representative point for a null geometry.");
        }
              
        Point centroid = geometry.getCentroid();
        
        /*
         * If the centroid is contained by the geometry, we're good to go
         */
        if(geometry.contains(centroid)){            
            return centroid;
        }
        
        
        /*
         * If the geometry is a collection we use the geometry with the greatest area.
         */
        if(geometry.getNumGeometries() > 1 && geometry instanceof GeometryCollection){
            GeometryCollection collection = (GeometryCollection) geometry;
            Geometry max = null;
            for(int i = 0; i < collection.getNumGeometries(); i++){
                Geometry other = collection.getGeometryN(i);
                if(max == null){
                    max = other;
                }
                else if(areaOrLength(other) > areaOrLength(max)){
                    max = other;
                }
            }
            
            //if we found a max area sub geometry use it.
            if(max != null){
                return representativePoint(max);
            }            
        }
        
        
        /*
         * Split the geometry into 4 halves by bisecting around the centroid N/S and E/W
         */
        GeometryFactory gf = new GeometryFactory();
        
        //bounds the entire geometry
        Envelope envelope = geometry.getEnvelopeInternal();
        
        //North
        Envelope northEnvelope = new Envelope(envelope.getMinX(), envelope.getMaxX(), envelope.getMinY(), centroid.getY());
        Geometry north = geometry.intersection(gf.toGeometry(northEnvelope));
        Geometry max = north;
        
        //South
        Envelope southEnvelope = new Envelope(envelope.getMinX(), envelope.getMaxX(), envelope.getMaxY(), centroid.getY());
        Geometry south = geometry.intersection(gf.toGeometry(southEnvelope));
        if(areaOrLength(south) > areaOrLength(max)){
            max = south;
        }
        
        //East
        Envelope eastEnvelope = new Envelope(envelope.getMaxX(), centroid.getX(), envelope.getMinY(), envelope.getMaxY());
        Geometry east = geometry.intersection(gf.toGeometry(eastEnvelope));
        if(areaOrLength(east) > areaOrLength(max)){
            max = east;
        }
        
        //West
        Envelope westEnvelope = new Envelope(envelope.getMinX(), centroid.getX(), envelope.getMinY(), envelope.getMaxY());
        Geometry west = geometry.intersection(gf.toGeometry(westEnvelope));
        if(areaOrLength(west) > areaOrLength(max)){
            max = west;
        }
        
        //use the half with the max area
                        
        return representativePoint(max);
        
    }
    
    

    private static double areaOrLength(Geometry geometry){
        double result = geometry.getArea();
        if(result == 0){
            result = geometry.getLength();
        }
        return result;
    }
}