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 RegionChooserUtils { 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 RegionChooserUtils() { // 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 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 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")); } }