GeoUtils.java

package eu.simstadt.geo;

import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.text.DecimalFormat;
import java.text.DecimalFormatSymbols;
import java.util.Optional;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import java.util.stream.Stream;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Polygon;
import org.osgeo.proj4j.BasicCoordinateTransform;
import org.osgeo.proj4j.CRSFactory;
import org.osgeo.proj4j.CoordinateReferenceSystem;
import org.osgeo.proj4j.ProjCoordinate;


public class GeoUtils
{
	private static final CRSFactory CRS_FACTORY = new CRSFactory();
	public static final CoordinateReferenceSystem WGS84 = CRS_FACTORY.createFromName("EPSG:4326");
	private static final Pattern srsNamePattern = Pattern.compile("(?i)(?<=srsName=[\"'])[^\"']+(?=[\"'])");

	private GeoUtils() {
		// only static use
	}

	/**
	 * Writes down the input coordinates in a string, without "." or "-". E.g. "N48_77__E9_18" for (48.77, 9.18) and
	 * "S12_345__W6_789" for (-12.345, -6.789).
	 * 
	 * Useful for caching geolocated information. The process should be reversible.
	 * 
	 * WARNING: LATITUDE FIRST!
	 * 
	 * @param latitude
	 * @param longitude
	 * @return Filename friendly string for given coordinates
	 */
	public static String coordinatesForFilename(double latitude, double longitude) {
		DecimalFormatSymbols symbols = DecimalFormatSymbols.getInstance();
		symbols.setDecimalSeparator('.');
		DecimalFormat df = new DecimalFormat("+#.######;-#.######", symbols);
		String lat = df.format(latitude).replace('-', 'S').replace('+', 'N');
		String lon = df.format(longitude).replace('-', 'W').replace('+', 'E');
		return (lat + "__" + lon).replace('.', '_');
	}

	/**
	 * The srsName (The name by which this reference system is identified) inside the CityGML file can have multiple
	 * formats. This method tries to parse the string and detect the corresponding reference system. If it is found, it
	 * returns a proj4j.CoordinateReferenceSystem. It throws an IllegalArgumentException otherwise.
	 * 
	 * This method should be able to parse any EPSG id : e.g. "EPSG:1234". German Citygmls might also have "DE_DHDN_3GK3"
	 * or "ETRS89_UTM32" as srsName, so those are also included. It isn't guaranteed that those formats are correctly
	 * parsed, though.
	 * 
	 * The EPSG ids and parameters are defined in resources ('nad/epsg') inside proj4j-0.1.0.jar. Some EPSG ids are
	 * missing though, e.g. 7415
	 * 
	 * @param srsName
	 * @return CoordinateReferenceSystem
	 */
	public static CoordinateReferenceSystem crsFromSrsName(String srsName) {
		// EPSG:31467
		Pattern pEPSG = Pattern.compile("^(EPSG:\\d+)$");
		Matcher mEPSG = pEPSG.matcher(srsName);
		if (mEPSG.find()) {
			return CRS_FACTORY.createFromName(srsName);
		}

		// urn:ogc:def:crs,crs:EPSG:6.12:31467,crs:EPSG:6.12:5783
		//   or
		// urn:ogc:def:crs,crs:EPSG::28992
		Pattern pOGC = Pattern.compile("urn:ogc:def:crs(?:,crs)?:EPSG:[\\d\\.]*:([\\d]+)\\D*");
		Matcher mOGC = pOGC.matcher(srsName);
		if (mOGC.find()) {
			return CRS_FACTORY.createFromName("EPSG:" + mOGC.group(1));
		}
		// urn:adv:crs:DE_DHDN_3GK3*DE_DHHN92_NH
		// urn:adv:crs:ETRS89_UTM32*DE_DHHN92_NH
		Pattern pURN = Pattern.compile("urn:adv:crs:([^\\*]+)");
		Matcher mURN = pURN.matcher(srsName);
		//NOTE: Could use a HashMap if the switch/case becomes too long.
		if (mURN.find()) {
			switch (mURN.group(1)) {
			case "DE_DHDN_3GK2":
				return CRS_FACTORY.createFromName("EPSG:31466");
			case "DE_DHDN_3GK3":
				return CRS_FACTORY.createFromName("EPSG:31467");
			case "DE_DHDN_3GK4":
				return CRS_FACTORY.createFromName("EPSG:31468");
			case "DE_DHDN_3GK5":
				return CRS_FACTORY.createFromName("EPSG:31469");
			case "ETRS89_UTM32":
				return CRS_FACTORY.createFromName("EPSG:25832");
			default:
				// nothing found
			}
		}
		throw new IllegalArgumentException("Unknown srsName format: " + srsName);
	}

	/**
	 * Converts a jts.geom.Polygon from one CoordinateReferenceSystem to another.
	 * 
	 * NOTE: It would be easier with org.geotools.referencing.CRS instead of Proj4J
	 * 
	 * @param polygonInOriginalCRS
	 * @param originalCRS
	 * @param newCRS
	 * @return
	 */
	public static Polygon changePolygonCRS(Polygon polygonInOriginalCRS, CoordinateReferenceSystem originalCRS,
			CoordinateReferenceSystem newCRS) {
		GeometryFactory geometryFactory = new GeometryFactory();
		Coordinate[] convexHullcoordinates = polygonInOriginalCRS.getCoordinates();
		BasicCoordinateTransform transformToWgs84 = new BasicCoordinateTransform(originalCRS, newCRS);
		for (int i = 0; i < convexHullcoordinates.length; i++) {
			ProjCoordinate wgs84Coordinate = transformToWgs84.transform(
					new ProjCoordinate(convexHullcoordinates[i].x, convexHullcoordinates[i].y), new ProjCoordinate());
			convexHullcoordinates[i] = new Coordinate(wgs84Coordinate.x, wgs84Coordinate.y);
		}

		return geometryFactory.createPolygon(convexHullcoordinates);
	}

	/**
	 * 
	 * Fast scan of the 50 first lines of a Citygml file to look for srsName. It might not be as reliable as parsing the
	 * whole CityGML, but it should be much faster and use much less memory. For a more reliable solution, use
	 * GeoCoordinatesAccessor. This solution can be convenient for Web services, RegionChooser or HullExtractor.
	 * 
	 * @param citygmlPath
	 * @return
	 * @throws IOException
	 */
	public static CoordinateReferenceSystem crsFromCityGMLHeader(Path citygmlPath) throws IOException {
		Optional<String> line = Files.lines(citygmlPath).limit(50).filter(srsNamePattern.asPredicate()).findFirst();
		if (line.isPresent()) {
			Matcher matcher = srsNamePattern.matcher(line.get());
			matcher.find();
			return crsFromSrsName(matcher.group());
		} else {
			throw new IllegalArgumentException("No srsName found in the header of " + citygmlPath);
		}
	}

	/**
	 * The haversine formula determines the great-circle distance between two points on a sphere given their longitudes
	 * and latitudes. https://en.wikipedia.org/wiki/Haversine_formula
	 * 
	 * @param latitude1
	 * @param longitude1
	 * @param latitude2
	 * @param longitude2
	 * @return distance in kilometer
	 */
	public static Double greatCircleDistance(Double latitude1, Double longitude1, Double latitude2, Double longitude2) {
		final int earthRadius = 6371; // [km]. Radius of the earth

		double latDistance = Math.toRadians(latitude2 - latitude1);
		double lonDistance = Math.toRadians(longitude2 - longitude1);
		double a = Math.sin(latDistance / 2) * Math.sin(latDistance / 2) + Math.cos(Math.toRadians(latitude1))
				* Math.cos(Math.toRadians(latitude2)) * Math.sin(lonDistance / 2) * Math.sin(lonDistance / 2);
		double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));
		return earthRadius * c;
	}

	/**
	 * Finds every CityGML in every .proj folder in a repository.
	 * 
	 * @param repository
	 * @return a stream of CityGML Paths
	 * @throws IOException
	 */
	public static Stream<Path> everyCityGML(Path repository) throws IOException {
		return Files.walk(repository, 2)
				.sorted()
				.filter(
						p -> Files.isRegularFile(p) &&
								p.getParent().toString().endsWith(".proj") &&
								p.toString().toLowerCase().endsWith(".gml"));
	}

}