Public
Snippet $34902 authored by Ian Turton

Search for a projection

ProjectionGuesser.java
package spike;

import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStreamReader;
import java.net.URL;
import java.net.URLEncoder;
import java.util.Collection;
import java.util.HashSet;
import java.util.Set;

import javax.measure.converter.UnitConverter;
import javax.measure.unit.SI;
import javax.measure.unit.Unit;
import javax.xml.parsers.DocumentBuilder;
import javax.xml.parsers.DocumentBuilderFactory;
import javax.xml.parsers.ParserConfigurationException;

import org.apache.commons.cli.CommandLine;
import org.apache.commons.cli.CommandLineParser;
import org.apache.commons.cli.DefaultParser;
import org.apache.commons.cli.HelpFormatter;
import org.apache.commons.cli.Options;
import org.apache.commons.cli.ParseException;
import org.geotools.factory.Hints;
import org.geotools.geometry.jts.JTS;
import org.geotools.metadata.iso.IdentifierImpl;
import org.geotools.metadata.iso.extent.GeographicBoundingBoxImpl;
import org.geotools.referencing.CRS;
import org.geotools.referencing.ReferencingFactoryFinder;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.opengis.geometry.MismatchedDimensionException;
import org.opengis.metadata.Identifier;
import org.opengis.metadata.extent.Extent;
import org.opengis.metadata.extent.GeographicExtent;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.ReferenceIdentifier;
import org.opengis.referencing.crs.CRSAuthorityFactory;
import org.opengis.referencing.crs.ProjectedCRS;
import org.opengis.referencing.cs.AxisDirection;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import org.w3c.dom.Document;
import org.w3c.dom.Element;
import org.w3c.dom.Node;
import org.w3c.dom.NodeList;
import org.xml.sax.SAXException;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.Point;

/**
 * Answer questions of the type "I have a point(459579,9894182) and it should be
 * in Arizona, what projection is it?"
 * 
 * @author ian
 *
 */
public class ProjectionGuesser {

  public static void main(String[] args) {
    double x = 0, y = 0;
    String state = "";
    String country = "";
    ProjectionGuesser guesser = new ProjectionGuesser();
    if (args.length == 0) {
      BufferedReader reader = new BufferedReader(new InputStreamReader(System.in));

      try {
        System.out.print("Enter point x,y: ");
        String line = reader.readLine();
        String[] p = line.split(",");
        x = Double.parseDouble(p[0]);
        y = Double.parseDouble(p[1]);
        System.out.print("Enter State, Country: ");
        line = reader.readLine();
        String s[] = line.split(",");
        if (s.length > 1) {
          state = s[0];
          country = s[1];
        } else {
          country = s[0];
        }
      } catch (IOException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
      }
    } else {
      CommandLineParser parser = new DefaultParser();

      Options options = new Options();
      options.addOption("x", true, "X point");
      options.addOption("y", true, "Y point");
      options.addOption("s", true, "State/Subdivision");
      options.addOption("c", true, "Country");
      options.addOption("f", false, "Limit search to just EPSG");
      CommandLine commandLine;
      try {
        commandLine = parser.parse(options, args);
        if (commandLine.getOptions().length == 0) {
          HelpFormatter formatter = new HelpFormatter();
          formatter.printHelp("ProjectionGuesser", options, true);

          System.exit(1);
        }
        if(commandLine.hasOption('f')) {
          guesser.quick = true;
        }
        if (commandLine.hasOption("x")) {
          x = Double.parseDouble(commandLine.getOptionValue("x"));
        }
        if (commandLine.hasOption("y")) {
          y = Double.parseDouble(commandLine.getOptionValue("y"));
        }
        if (commandLine.hasOption("s")) {
          state = commandLine.getOptionValue("s");
        }
        if (commandLine.hasOption("c")) {
          country = commandLine.getOptionValue("c");
        }

      } catch (ParseException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
      }

    }
    
    guesser.setCountry(country);
    if (!state.isEmpty()) {
      guesser.setState(state);
    }
    Coordinate c = new Coordinate(x, y);
    guesser.setCoordinate(c);
    ProjectedCRS crs = guesser.guess();
    StringBuilder idB = new StringBuilder();
    idB.append(crs.getName()).append(" ");
    Set<ReferenceIdentifier> identifiers = crs.getIdentifiers();
    for (ReferenceIdentifier i : identifiers) {
      idB.append(i.toString()).append(" ");
    }
    String id = idB.toString().trim();
    System.out.println("best guess is: " + id);
    System.out.println("Error: " + guesser.close + "m");
    System.out.println("\nWKT:\n" + crs.toWKT());
  }

  private double close;
  private boolean quick=false;

  public ProjectedCRS guess() {
    //
    Coordinate wgsPoint = lookupLocation();
    System.out.println(wgsPoint);
    if (wgsPoint == null) {
      // lookup failed
      System.err.println("couldn't generate a point");
      System.exit(3);
    }
    // Look for a CRS
    //
    Set<ProjectedCRS> possibles = new HashSet<>();
    if (quick) {
      String authority = "EPSG";
      possibles = fetchPossibles(wgsPoint, authority);
    } else {
      Set<String> auths = ReferencingFactoryFinder.getAuthorityNames();
      for (String auth : auths) {
        possibles.addAll(fetchPossibles(wgsPoint, auth));
      }
    }
    GeometryFactory geometryFactory = new GeometryFactory();
    Point testPoint = geometryFactory.createPoint(wgsPoint);
    Point testPoint2 = geometryFactory.createPoint(new Coordinate(wgsPoint.y, wgsPoint.x));
    Point targetPoint = geometryFactory.createPoint(coordinate);
    ProjectedCRS bestCRS = null;
    close = Double.POSITIVE_INFINITY;
    for (ProjectedCRS proj : possibles) {
      try {
        Unit<?> unit = proj.getCoordinateSystem().getAxis(0).getUnit();
        MathTransform testTransform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, proj);
        Point p;
        if (proj.getCoordinateSystem().getAxis(0).getDirection().equals(AxisDirection.NORTH)) {
          p = (Point) JTS.transform(testPoint, testTransform);
        } else {
          p = (Point) JTS.transform(testPoint2, testTransform);
        }
        double dist = targetPoint.distance(p);
        if (!unit.equals(SI.METER)) {// normalise the distance otherwise feet
                                     // never win
          UnitConverter conv = unit.getConverterTo(SI.METER);
          dist = conv.convert(dist);
        }

        if (dist < close) {
          close = dist;
          bestCRS = proj;
        }
      } catch (FactoryException | MismatchedDimensionException | TransformException e) {
        // lots of things can go wrong with random tests so ignore them
      }
    }

    return bestCRS;
  }

  /**
   * @param wgsPoint
   * @return
   */
  private Set<ProjectedCRS> fetchPossibles(Coordinate wgsPoint, String authority) {
    Set<ProjectedCRS> possibles = new HashSet<>();
    Hints hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.FALSE);
    Set<CRSAuthorityFactory> facs = ReferencingFactoryFinder.getCRSAuthorityFactories(hints);
    Set<String> codes = new HashSet<>();
    CRSAuthorityFactory factory = null;
    int bestCodes = 0;
    Set<String> c;
    IdentifierImpl identifier = new IdentifierImpl(authority);
    for (CRSAuthorityFactory fac : facs) {
      Collection<? extends Identifier> identifiers = fac.getAuthority().getIdentifiers();

      if (identifiers.contains(identifier)) {

        try {
          c = fac.getAuthorityCodes(org.opengis.referencing.crs.ProjectedCRS.class);
          if (c.size() > bestCodes) {
            bestCodes = c.size();
            factory = fac;
            codes = c;
          }
        } catch (FactoryException e) {
          // TODO Auto-generated catch block
          e.printStackTrace();
        }
      }
    }
    
    for (String code : codes) {
      try {
        ProjectedCRS crs = factory.createProjectedCRS(code);
        Extent valid = crs.getDomainOfValidity();

        Collection<? extends GeographicExtent> geographicElements = valid.getGeographicElements();

        for (GeographicExtent bbox : geographicElements) {

          if (bbox instanceof GeographicBoundingBoxImpl) {
            GeographicBoundingBoxImpl gbbox = ((GeographicBoundingBoxImpl) bbox);

            if (gbbox.getNorthBoundLatitude() > wgsPoint.x && gbbox.getSouthBoundLatitude() < wgsPoint.x
                && gbbox.getWestBoundLongitude() < wgsPoint.y && gbbox.getEastBoundLongitude() > wgsPoint.y) {
              possibles.add(crs);
            }
          }
        }

      } catch (FactoryException e) {
        // TODO Auto-generated catch block
        // e.printStackTrace();
      }

    }
    return possibles;
  }

  private String country = "";
  private String state = "";
  private Coordinate coordinate;

  public void setCoordinate(Coordinate c) {
    this.coordinate = c;
  }

  public void setState(String state) {
    this.state = state;
  }

  public void setCountry(String country) {
    this.country = country;
  }

  private Coordinate lookupLocation() {
    String base = "http://api.geonames.org/search?style=SHORT&isNameRequired=true&featureClass=A&username=ianturton";
    String name = "";
    if (state.isEmpty()) {
      name = country;
    } else {
      name = state + "," + country;
    }
    String urlEncodedName;
    try {
      urlEncodedName = URLEncoder.encode(name, "UTF-8");

      String request = base + "&q=" + urlEncodedName;
      // System.out.println(request);
      URL url = new URL(request);
      DocumentBuilder builder = DocumentBuilderFactory.newInstance().newDocumentBuilder();
      Document doc = builder.parse(url.openStream());
      NodeList nodes = doc.getElementsByTagName("geoname");
      if (nodes.getLength() == 0) {
        System.out.println("No matches for '" + name + "' did you spell it right?");
        return null;
      }
      Element top = (Element) nodes.item(0);
      NodeList kids = top.getChildNodes();
      double lat = 0, lng = 0;
      for (int i = 0; i < kids.getLength(); i++) {
        Node node = kids.item(i);
        if (!(node instanceof Element))
          continue;
        Element kid = (Element) node;
        if (kid.getTagName().equalsIgnoreCase("name")) {
          System.out.println("Working with a guess of " + kid.getTextContent());
        }
        if (kid.getTagName().equalsIgnoreCase("lat")) {
          lat = Double.parseDouble(kid.getTextContent());
        }
        if (kid.getTagName().equalsIgnoreCase("lng")) {
          lng = Double.parseDouble(kid.getTextContent());
        }
      }
      Coordinate ret = new Coordinate(lat, lng);
      return ret;
    } catch (IOException | ParserConfigurationException | SAXException e) {
      // TODO Auto-generated catch block
      e.printStackTrace();
    }
    return null;
  }

}