diff --git a/src/main/java/pl/wat/ms4ds/terenfunkcje/Coord.java b/src/main/java/pl/wat/ms4ds/terenfunkcje/Coord.java index 88ef67e..b8a712a 100644 --- a/src/main/java/pl/wat/ms4ds/terenfunkcje/Coord.java +++ b/src/main/java/pl/wat/ms4ds/terenfunkcje/Coord.java @@ -1,7 +1,13 @@ package pl.wat.ms4ds.terenfunkcje; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; +import pl.wat.ms4ds.terenfunkcje.konwersja.PUWGCoord; + public class Coord { + private static final Logger logger = LoggerFactory.getLogger(Coord.class); + public static class Geo { public double lat; public double lon; @@ -18,6 +24,14 @@ public class Coord { lat = other.lat; lon = other.lon; } + + @Override + public String toString() { + return "Geo{" + + "lat=" + lat + + ", lon=" + lon + + '}'; + } } /** @@ -45,6 +59,14 @@ public class Coord { this.easting = easting; this.northing = northing; } + + @Override + public String toString() { + return "Puwg{" + + "easting=" + easting + + ", northing=" + northing + + '}'; + } } public static class Grid { @@ -74,15 +96,12 @@ public class Coord { y = zamienWspYmsNaIdKwadratuY(yms); } - + @Override public String toString() { - StringBuilder sb = new StringBuilder("(x= "); - sb.append(x); - sb.append(", y= "); - sb.append(y); - sb.append(')'); -// String s = "(x= " + Integer.toString(x) + ", y= " + Integer.toString(y) + ")"; - return sb.toString(); + return "Grid{" + + "x=" + x + + ", y=" + y + + '}'; } private static final double ODWROT_SS_DX_MS = 1.0 / MapConsts.SS_DX_MS; @@ -311,43 +330,198 @@ public class Coord { return (int) yy; } - /** - * Wyznacza wspolrzedne lewego dolnego wierzcholka oraz prawego gornego wierzcholka - * prostokata opisanego na obszarze: rejon. - * - * @param rejon - * @param min - * @param max - */ - public static void minMaxIdKwadratu(Grid[] rejon, Grid min, Grid max) { - if (rejon == null) { - return; - } - if (min == null) { - min = new Grid(rejon[0].x, rejon[0].y); - } - if (max == null) { - max = new Grid(rejon[0].x, rejon[0].y); - } - for (int i = 1; i < rejon.length; i++) { - if (min.x > rejon[i].x) { - min.x = rejon[i].x; - } - if (min.y > rejon[i].y) { - min.y = rejon[i].y; - } - if (max.x < rejon[i].x) { - max.x = rejon[i].x; - } - if (max.y < rejon[i].y) { - max.y = rejon[i].y; - } - } - } - - private static final double ODWROT_SS_DX_MS = 1.0 / MapConsts.SS_DX_MS; private static final double ODWROT_SS_DY_MS = 1.0 / MapConsts.SS_DY_MS; + /** + * Funkcja służy do konwersji współrzednych elipsoidalnych WGS84 (lat/lon) na płaskie X-northing, Y-easting + * odwzorowania kartograficznego 1992. + * + * @param puwgCoord współrzędne PUWG1992 (easting, northing) [metry] + * @param lat szerokość geograficzna WSG-84 [stopnie dziesiętnie] + * @param lon długość geograficzna WSG-84 [stopnie dziesiętnie] + */ + public static void convertWGS84ToPUWG1992(double lat, double lon, PUWGCoord puwgCoord) { + if (lon < 13.5 || lon > 25.5) { + //Błędna wartość długości geograficznej (zwracana wartość 999999999999999) + puwgCoord.easting = 999999999999999.0; + puwgCoord.northing = 999999999999999.0; + return; + } + double latRad = lat * DEG_2_RAD; + double dlam = (lon - 19.0) * DEG_2_RAD; + double dlam_pow_2 = dlam * dlam; + double dlam_pow_3 = dlam_pow_2 * dlam; + double dlam_pow_4 = dlam_pow_3 * dlam; + double s = Math.sin(latRad); + double c = Math.cos(latRad); + double c_pow_2 = c * c; + double c_pow_3 = c_pow_2 * c; + double c_pow_4 = c_pow_3 * c; + double t = s / c; + double t_pow_2 = t * t; + double t_pow_3 = t_pow_2 * t; + double t_pow_4 = t_pow_3 * t; + double t_pow_5 = t_pow_4 * t; + double eta = E2_SQUARED * c_pow_2; + double eta_pow_2 = eta * eta; + double eta_pow_3 = eta_pow_2 * eta; + double eta_pow_4 = eta_pow_3 * eta; + double sn = sphsn(latRad); + double tmd = sphtmd(latRad); + double t1, t2, t3, t4, t5; + t1 = tmd * OK; + double sns = sn * s; + t2 = sns * c * OK / 2.0; + t3 = sns * c_pow_3 * OK * (5.0 - t_pow_2 + 9.0 * eta + 4.0 * eta_pow_2) / 24.0; + t4 = sns * c_pow_4 * c * OK * (61.0 - 58.0 * t_pow_2 + t_pow_4 + + 270.0 * eta - 330.0 * t_pow_2 * eta + 445.0 * eta_pow_2 + + 324.0 * eta_pow_3 - 680.0 * t_pow_2 * eta_pow_2 + + 88.0 * eta_pow_4 - 600.0 * t_pow_2 * eta_pow_3 - 192.0 * t_pow_2 * eta_pow_4) / 720.0; + t5 = sns * c_pow_4 * c_pow_3 * OK * (1385.0 - 3111.0 * t_pow_2 + + 543.0 * t_pow_4 - t_pow_5 * t) / 40320.0; + puwgCoord.northing = -5300000.0 + t1 + dlam_pow_2 * t2 + dlam_pow_4 * t3 + + dlam_pow_4 * dlam_pow_2 * t4 + dlam_pow_4 * dlam_pow_4 * t5; + t1 = sn * c * OK; + t2 = sn * c_pow_3 * OK * (1.0 - t_pow_2 + eta) / 6.0; + t3 = sn * c_pow_4 * c * OK * (5.0 - 18.0 * t_pow_2 + t_pow_4 + 14.0 * eta + - 58.0 * t_pow_2 * eta + 13.0 * eta_pow_2 + 4.0 * eta_pow_3 + - 64.0 * t_pow_2 * eta_pow_2 - 24.0 * t_pow_2 * eta_pow_3) / 120.0; + t4 = sn * c_pow_4 * c_pow_3 * OK * (61.0 - 479.0 * t_pow_2 + 179.0 * t_pow_4 - t_pow_5 * t) / 5040.0; + puwgCoord.easting = 500000.0 + dlam * t1 + dlam_pow_3 * t2 + + dlam_pow_4 * dlam * t3 + dlam_pow_4 * dlam_pow_3 * t4;// + 0.5; + } + + /** + * Funkcja do konwersji współrzędnych płaskich X/Y odwzorowania kartograficznego 1992 na elipsoidalne lat/lon elipsoide WGS84. + *

+ * PUWGCoord.proj: odwzorowanie kartograficzne (proj = 1 odpowiada odwzorowaniu 1992, natomiast każda inna odwzorowaniu 2000) + * + * @param northing współrzędne na osi OY odwzorowania kartograficznego PUWG-1992 do konwersji [metry] + * @param easting współrzędne na osi OX odwzorowania kartograficznego PUWG-1992 do konwersji [metry] + * @param geoCoord współrzędne geograficzne odwzorowania WGS-84 po konwersji [stopnie] + */ + public static void convertPUWG1992ToWGS84(double northing, double easting, Coord.Geo geoCoord) { + double tmd = (northing + 5300000.0) / OK; + double sr = sphsr(0.0); + double ftphi = tmd / sr; + for (int i = 0; i < 5; i++) { + ftphi += (tmd - sphtmd(ftphi)) / sphsr(ftphi); + } + sr = sphsr(ftphi); + double sn = sphsn(ftphi); + double sn_pow_2 = sn * sn; + double sn_pow_3 = sn_pow_2 * sn; + double sn_pow_4 = sn_pow_3 * sn; + double sn_pow_5 = sn_pow_4 * sn; + double sn_pow_7 = sn_pow_5 * sn_pow_2; + double s = Math.sin(ftphi); + double c = Math.cos(ftphi); + double t = s / c; + double t_pow_2 = t * t; + double t_pow_4 = t_pow_2 * t_pow_2; + double t_pow_6 = t_pow_4 * t_pow_2; + double eta = E2_SQUARED * (c * c); + double eta_pow_2 = eta * eta; + double eta_pow_3 = eta_pow_2 * eta; + double eta_pow_4 = eta_pow_2 * eta_pow_2; + double de = easting - 500000.0; + double de_pow_2 = de * de; + double de_pow_3 = de_pow_2 * de; + double de_pow_4 = de_pow_3 * de; + double t0, t1, t2, t3; + t0 = t / (2.0 * sr * sn * OK_POW_2); + t1 = t * (5.0 + 3.0 * t_pow_2 + eta - 4.0 * eta_pow_2 - 9.0 * t_pow_2 * eta) / (24.0 * sr * sn_pow_3 * OK_POW_4); + t2 = t * (61.0 + 90.0 * t_pow_2 + 46.0 * eta + 45.0 * t_pow_4 - 252.0 * t_pow_2 * eta - 3.0 * eta_pow_2 + + 100.0 * eta_pow_3 - 66.0 * t_pow_2 * eta_pow_2 - 90.0 * t_pow_4 * eta + 88.0 * eta_pow_4 + + 225.0 * t_pow_4 * eta_pow_2 + 84.0 * t_pow_2 * eta_pow_3 - 192.0 * t_pow_2 * eta_pow_4) / (720.0 * sr * sn_pow_5 * OK_POW_6); + t3 = t * (1385.0 + 3633 * t_pow_2 + 4095.0 * t_pow_4 + 1575.0 * t_pow_6) / (40320 * sr * sn_pow_7 * (OK_POW_8)); + geoCoord.lat = ftphi - de_pow_2 * t0 + de_pow_4 * t1 - de_pow_3 * de_pow_3 * t2 + de_pow_4 * de_pow_3 * t3; + t0 = 1.0 / (sn * c * OK); + t1 = (1.0 + 2.0 * t_pow_2 + eta) / (6.0 * sn_pow_3 * c * (OK_POW_3)); + t2 = (5.0 + 6.0 * eta + 28.0 * t_pow_2 - 3.0 * eta_pow_2 + 8.0 * t_pow_2 * eta + + 24.0 * t_pow_4 - 4.0 * eta_pow_3 + 4.0 * t_pow_2 * eta_pow_2 + + 24.0 * t_pow_2 * eta_pow_3) / (120.0 * sn_pow_5 * c * (OK_POW_5)); + t3 = (61.0 + 662.0 * t_pow_2 + 1320.0 * t_pow_4 + 720.0 * t_pow_6) / (5040.0 * sn_pow_7 * c * OK_POW_7); + double dlam = de * t0 - de_pow_3 * t1 + de_pow_3 * de_pow_2 * t2 - de_pow_3 * de_pow_4 * t3; + // 19.0 * DEG_2_RAD == 0.33161255787892263; +// geoCoord.lon = 0.33161255787892263 + dlam; +// geoCoord.lon *= RAD_2_DEG; + geoCoord.lon = dlam * RAD_2_DEG; + geoCoord.lon += 19.0; + geoCoord.lat *= RAD_2_DEG; + } + + //////////////////////////////////////////////////////////////////////////////// + // Funkcje pomocnicze i stałe + + /// ///////////////////////////////////////////////////////////////////////////// + + static double calculateESquared(double a, double b) { + a *= a; + b *= b; + return (a - b) / a; + } + + static double calculateE2Squared(double a, double b) { + a *= a; + b *= b; + return (a - b) / b; + } + + static double denom(double sphi) { + double sinSphi = Math.sin(sphi); + return Math.sqrt(1.0 - E_SQUARED * (sinSphi * sinSphi)); + } + + static double sphsr(double sphi) { + double dn = denom(sphi); + return A * (1.0 - E_SQUARED) / (dn * dn * dn); + } + + static double sphsn(double sphi) { + double sinSphi = Math.sin(sphi); + return A / Math.sqrt(1.0 - E_SQUARED * (sinSphi * sinSphi)); + } + + static double sphtmd(double sphi) { + return (AP * sphi) - (BP * Math.sin(2.0 * sphi)) + (CP * Math.sin(4.0 * sphi)) + - (DP * Math.sin(6.0 * sphi)) + (EP * Math.sin(8.0 * sphi)); + } + + private static final double DEG_2_RAD = Math.PI / 180.0; + private static final double RAD_2_DEG = 180.0 / Math.PI; + + /** + * Dlługość dużej półsi, w metrach dla elipsoidy WGS-84. + */ + private static final double A = 6378137.0; + + /** + * double f: spłaszczenie elipsoidalne dla elipsoidy WGS-84, 1 / 298.257223563 + */ +// private static final double F = 1.0 / 298.257223563; + private static final double B = A * (1.0 - 1.0 / 298.257223563); + private static final double E_SQUARED = calculateESquared(A, B); + private static final double E2_SQUARED = calculateE2Squared(A, B); + private static final double TN = (A - B) / (A + B); + private static final double AP = A * (1.0 - TN + 5.0 * (TN * TN - TN * TN * TN) / 4.0 + 81.0 * TN * TN * TN * TN - TN * TN * TN * TN * TN / 64.0); + private static final double BP = 3.0 * A * (TN - TN * TN + 7.0 * (TN * TN * TN - TN * TN * TN * TN) / 8.0 + 55.0 * TN * TN * TN * TN * TN / 64.0) / 2.0; + private static final double CP = 15.0 * A * (TN * TN - TN * TN * TN + 3.0 * (TN * TN * TN * TN - TN * TN * TN * TN * TN) / 4.0) / 16.0; + private static final double DP = 35.0 * A * (TN * TN * TN - TN * TN * TN * TN + 11.0 * TN * TN * TN * TN * TN / 16.0) / 48.0; + private static final double EP = 315.0 * A * (TN * TN * TN * TN - TN * TN * TN * TN * TN) / 512.0; + + /** + * Współczynnik zniekształcenia skali mapy w południku osiowym dla odwzorowania kartograficznego PUWG-1992. + */ + private static final double OK = 0.9993; + private static final double OK_POW_2 = OK * OK; + private static final double OK_POW_3 = OK_POW_2 * OK; + private static final double OK_POW_4 = OK_POW_3 * OK; + private static final double OK_POW_5 = OK_POW_4 * OK; + private static final double OK_POW_6 = OK_POW_5 * OK; + private static final double OK_POW_7 = OK_POW_6 * OK; + private static final double OK_POW_8 = OK_POW_7 * OK; + } diff --git a/src/main/java/pl/wat/ms4ds/terenfunkcje/Square.java b/src/main/java/pl/wat/ms4ds/terenfunkcje/Square.java index 5c5a704..efaaf1b 100644 --- a/src/main/java/pl/wat/ms4ds/terenfunkcje/Square.java +++ b/src/main/java/pl/wat/ms4ds/terenfunkcje/Square.java @@ -213,8 +213,8 @@ public class Square { } } - final int x; - final int y; + public final int x; + public final int y; public double sumaWysokosci; public int count = 1; public double ell; @@ -238,7 +238,6 @@ public class Square { return result; } - public String toString() { StringBuilder linia = new StringBuilder(100); linia.append("["); diff --git a/src/main/java/pl/wat/ms4ds/terenfunkcje/Teren.java b/src/main/java/pl/wat/ms4ds/terenfunkcje/Teren.java index e3d2ee0..43b2bb2 100644 --- a/src/main/java/pl/wat/ms4ds/terenfunkcje/Teren.java +++ b/src/main/java/pl/wat/ms4ds/terenfunkcje/Teren.java @@ -259,7 +259,7 @@ public class Teren { public static Square getKwadratPUWG(double northing, double easting) { Coord.Geo geoCoord = new Coord.Geo(); - CoordUtils.convertPUWG1992ToWGS84(northing, easting, geoCoord); + Coord.convertPUWG1992ToWGS84(northing, easting, geoCoord); return getKwadrat(geoCoord.lat, geoCoord.lon); } diff --git a/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordTest.java b/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordTest.java index 41d8347..19ee620 100644 --- a/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordTest.java +++ b/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordTest.java @@ -18,20 +18,20 @@ public class CoordTest { Coord.Geo geoCoord = new Coord.Geo(); geoCoord.lon = 19; geoCoord.lat = 50; - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); + Coord.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing); double e = puwgCoord.easting; double n = puwgCoord.northing; geoCoord.lon = 20; geoCoord.lat = 51; - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); + Coord.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); double dx = puwgCoord.easting - e; double dy = puwgCoord.northing - n; logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing); - CoordUtils.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, geoCoord); + Coord.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, geoCoord); logger.debug("PUWG (e,n)=({}, {}) => Lat={}, Lon={}", puwgCoord.easting, puwgCoord.northing, geoCoord.lat, geoCoord.lon); // coord.proj = 1; @@ -52,24 +52,24 @@ public class CoordTest { puwgCoord.easting = 500000; puwgCoord.northing = 500000; logger.debug("----------------------------------"); - CoordUtils.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, geoCoord); + Coord.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, geoCoord); logger.debug("PUWG (e,n)=({}, {}) => Lat={}, Lon={}", puwgCoord.easting, puwgCoord.northing, geoCoord.lat, geoCoord.lon); - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); + Coord.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing); puwgCoord.easting = 100000; puwgCoord.northing = 470642; - CoordUtils.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, geoCoord); + Coord.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, geoCoord); logger.debug("PUWG (e,n)=({}, {}) => Lat={}, Lon={}", puwgCoord.easting, puwgCoord.northing, geoCoord.lat, geoCoord.lon); - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); + Coord.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing); puwgCoord.easting = 821310; puwgCoord.northing = 369750; - CoordUtils.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, geoCoord); + Coord.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, geoCoord); logger.debug("PUWG (e,n)=({}, {}) => Lat={}, Lon={}", puwgCoord.easting, puwgCoord.northing, geoCoord.lat, geoCoord.lon); - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); + Coord.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing); diff --git a/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordUtils.java b/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordUtils.java index b1e1b28..9935099 100644 --- a/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordUtils.java +++ b/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordUtils.java @@ -4,6 +4,7 @@ import pl.wat.ms4ds.terenfunkcje.Coord; import pl.wat.ms4ds.terenfunkcje.Teren; import org.slf4j.Logger; import org.slf4j.LoggerFactory; +import pl.wat.ms4ds.terenfunkcje.nmt.NMTData; import java.io.BufferedReader; import java.io.FileReader; @@ -103,15 +104,15 @@ public class CoordUtils { } catch (NumberFormatException e) { logger.warn("Bledne dane w pliku: " + fileName); } - convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, latLon); + Coord.convertPUWG1992ToWGS84(puwgCoord.northing, puwgCoord.easting, latLon); Coord.Grid idKw = new Coord.Grid(latLon.lon, latLon.lat); NMTData daneWysok = daneWysokHashMap.get(idKw); if (daneWysok == null) { - daneWysok = new NMTData(idKw, wysokosc, 1); + daneWysok = new NMTData(idKw.x, idKw.y, wysokosc, 1); daneWysokHashMap.put(idKw, daneWysok); } else { - daneWysok.suma += wysokosc; - daneWysok.licz++; + daneWysok.sum += wysokosc; + daneWysok.count++; } line = br.readLine(); if (m++ % 100000 == 0) { @@ -127,194 +128,5 @@ public class CoordUtils { } } - /** - * Funkcja służy do konwersji współrzednych elipsoidalnych WGS84 (lat/lon) na płaskie X-northing, Y-easting - * odwzorowania kartograficznego 1992. - * - * @param puwgCoord współrzędne PUWG1992 (easting, northing) [metry] - * @param lat szerokość geograficzna WSG-84 [stopnie dziesiętnie] - * @param lon długość geograficzna WSG-84 [stopnie dziesiętnie] - */ - public static void convertWGS84ToPUWG1992(double lat, double lon, PUWGCoord puwgCoord) { - if (lon < 13.5 || lon > 25.5) { - //Błędna wartość długości geograficznej (zwracana wartość 999999999999999) - puwgCoord.easting = 999999999999999.0; - puwgCoord.northing = 999999999999999.0; - return; - } - double latRad = lat * DEG_2_RAD; - double dlam = (lon - 19.0) * DEG_2_RAD; - double dlam_pow_2 = dlam * dlam; - double dlam_pow_3 = dlam_pow_2 * dlam; - double dlam_pow_4 = dlam_pow_3 * dlam; - double s = Math.sin(latRad); - double c = Math.cos(latRad); - double c_pow_2 = c * c; - double c_pow_3 = c_pow_2 * c; - double c_pow_4 = c_pow_3 * c; - double t = s / c; - double t_pow_2 = t * t; - double t_pow_3 = t_pow_2 * t; - double t_pow_4 = t_pow_3 * t; - double t_pow_5 = t_pow_4 * t; - double eta = E2_SQUARED * c_pow_2; - double eta_pow_2 = eta * eta; - double eta_pow_3 = eta_pow_2 * eta; - double eta_pow_4 = eta_pow_3 * eta; - double sn = sphsn(latRad); - double tmd = sphtmd(latRad); - double t1, t2, t3, t4, t5; - t1 = tmd * OK; - double sns = sn * s; - t2 = sns * c * OK / 2.0; - t3 = sns * c_pow_3 * OK * (5.0 - t_pow_2 + 9.0 * eta + 4.0 * eta_pow_2) / 24.0; - t4 = sns * c_pow_4 * c * OK * (61.0 - 58.0 * t_pow_2 + t_pow_4 - + 270.0 * eta - 330.0 * t_pow_2 * eta + 445.0 * eta_pow_2 - + 324.0 * eta_pow_3 - 680.0 * t_pow_2 * eta_pow_2 - + 88.0 * eta_pow_4 - 600.0 * t_pow_2 * eta_pow_3 - 192.0 * t_pow_2 * eta_pow_4) / 720.0; - t5 = sns * c_pow_4 * c_pow_3 * OK * (1385.0 - 3111.0 * t_pow_2 - + 543.0 * t_pow_4 - t_pow_5 * t) / 40320.0; - puwgCoord.northing = -5300000.0 + t1 + dlam_pow_2 * t2 + dlam_pow_4 * t3 - + dlam_pow_4 * dlam_pow_2 * t4 + dlam_pow_4 * dlam_pow_4 * t5; - t1 = sn * c * OK; - t2 = sn * c_pow_3 * OK * (1.0 - t_pow_2 + eta) / 6.0; - t3 = sn * c_pow_4 * c * OK * (5.0 - 18.0 * t_pow_2 + t_pow_4 + 14.0 * eta - - 58.0 * t_pow_2 * eta + 13.0 * eta_pow_2 + 4.0 * eta_pow_3 - - 64.0 * t_pow_2 * eta_pow_2 - 24.0 * t_pow_2 * eta_pow_3) / 120.0; - t4 = sn * c_pow_4 * c_pow_3 * OK * (61.0 - 479.0 * t_pow_2 + 179.0 * t_pow_4 - t_pow_5 * t) / 5040.0; - puwgCoord.easting = 500000.0 + dlam * t1 + dlam_pow_3 * t2 - + dlam_pow_4 * dlam * t3 + dlam_pow_4 * dlam_pow_3 * t4;// + 0.5; - } - - /** - * Funkcja do konwersji współrzędnych płaskich X/Y odwzorowania kartograficznego 1992 na elipsoidalne lat/lon elipsoide WGS84. - *

- * PUWGCoord.proj: odwzorowanie kartograficzne (proj = 1 odpowiada odwzorowaniu 1992, natomiast każda inna odwzorowaniu 2000) - * - * @param northing współrzędne na osi OY odwzorowania kartograficznego PUWG-1992 do konwersji [metry] - * @param easting współrzędne na osi OX odwzorowania kartograficznego PUWG-1992 do konwersji [metry] - * @param geoCoord współrzędne geograficzne odwzorowania WGS-84 po konwersji [stopnie] - */ - public static void convertPUWG1992ToWGS84(double northing, double easting, Coord.Geo geoCoord) { - double tmd = (northing + 5300000.0) / OK; - double sr = sphsr(0.0); - double ftphi = tmd / sr; - for (int i = 0; i < 5; i++) { - ftphi += (tmd - sphtmd(ftphi)) / sphsr(ftphi); - } - sr = sphsr(ftphi); - double sn = sphsn(ftphi); - double sn_pow_2 = sn * sn; - double sn_pow_3 = sn_pow_2 * sn; - double sn_pow_4 = sn_pow_3 * sn; - double sn_pow_5 = sn_pow_4 * sn; - double sn_pow_7 = sn_pow_5 * sn_pow_2; - double s = Math.sin(ftphi); - double c = Math.cos(ftphi); - double t = s / c; - double t_pow_2 = t * t; - double t_pow_4 = t_pow_2 * t_pow_2; - double t_pow_6 = t_pow_4 * t_pow_2; - double eta = E2_SQUARED * (c * c); - double eta_pow_2 = eta * eta; - double eta_pow_3 = eta_pow_2 * eta; - double eta_pow_4 = eta_pow_2 * eta_pow_2; - double de = easting - 500000.0; - double de_pow_2 = de * de; - double de_pow_3 = de_pow_2 * de; - double de_pow_4 = de_pow_3 * de; - double t0, t1, t2, t3; - t0 = t / (2.0 * sr * sn * OK_POW_2); - t1 = t * (5.0 + 3.0 * t_pow_2 + eta - 4.0 * eta_pow_2 - 9.0 * t_pow_2 * eta) / (24.0 * sr * sn_pow_3 * OK_POW_4); - t2 = t * (61.0 + 90.0 * t_pow_2 + 46.0 * eta + 45.0 * t_pow_4 - 252.0 * t_pow_2 * eta - 3.0 * eta_pow_2 - + 100.0 * eta_pow_3 - 66.0 * t_pow_2 * eta_pow_2 - 90.0 * t_pow_4 * eta + 88.0 * eta_pow_4 - + 225.0 * t_pow_4 * eta_pow_2 + 84.0 * t_pow_2 * eta_pow_3 - 192.0 * t_pow_2 * eta_pow_4) / (720.0 * sr * sn_pow_5 * OK_POW_6); - t3 = t * (1385.0 + 3633 * t_pow_2 + 4095.0 * t_pow_4 + 1575.0 * t_pow_6) / (40320 * sr * sn_pow_7 * (OK_POW_8)); - geoCoord.lat = ftphi - de_pow_2 * t0 + de_pow_4 * t1 - de_pow_3 * de_pow_3 * t2 + de_pow_4 * de_pow_3 * t3; - t0 = 1.0 / (sn * c * OK); - t1 = (1.0 + 2.0 * t_pow_2 + eta) / (6.0 * sn_pow_3 * c * (OK_POW_3)); - t2 = (5.0 + 6.0 * eta + 28.0 * t_pow_2 - 3.0 * eta_pow_2 + 8.0 * t_pow_2 * eta - + 24.0 * t_pow_4 - 4.0 * eta_pow_3 + 4.0 * t_pow_2 * eta_pow_2 - + 24.0 * t_pow_2 * eta_pow_3) / (120.0 * sn_pow_5 * c * (OK_POW_5)); - t3 = (61.0 + 662.0 * t_pow_2 + 1320.0 * t_pow_4 + 720.0 * t_pow_6) / (5040.0 * sn_pow_7 * c * OK_POW_7); - double dlam = de * t0 - de_pow_3 * t1 + de_pow_3 * de_pow_2 * t2 - de_pow_3 * de_pow_4 * t3; - // 19.0 * DEG_2_RAD == 0.33161255787892263; -// geoCoord.lon = 0.33161255787892263 + dlam; -// geoCoord.lon *= RAD_2_DEG; - geoCoord.lon = dlam * RAD_2_DEG; - geoCoord.lon += 19.0; - geoCoord.lat *= RAD_2_DEG; - } - - //////////////////////////////////////////////////////////////////////////////// - // Funkcje pomocnicze i stałe - //////////////////////////////////////////////////////////////////////////////// - - static double calculateESquared(double a, double b) { - a *= a; - b *= b; - return (a - b) / a; - } - - static double calculateE2Squared(double a, double b) { - a *= a; - b *= b; - return (a - b) / b; - } - - static double denom(double sphi) { - double sinSphi = Math.sin(sphi); - return Math.sqrt(1.0 - E_SQUARED * (sinSphi * sinSphi)); - } - - static double sphsr(double sphi) { - double dn = denom(sphi); - return A * (1.0 - E_SQUARED) / (dn * dn * dn); - } - - static double sphsn(double sphi) { - double sinSphi = Math.sin(sphi); - return A / Math.sqrt(1.0 - E_SQUARED * (sinSphi * sinSphi)); - } - - static double sphtmd(double sphi) { - return (AP * sphi) - (BP * Math.sin(2.0 * sphi)) + (CP * Math.sin(4.0 * sphi)) - - (DP * Math.sin(6.0 * sphi)) + (EP * Math.sin(8.0 * sphi)); - } - - private static final double DEG_2_RAD = Math.PI / 180.0; - private static final double RAD_2_DEG = 180.0 / Math.PI; - - /** - * Dlługość dużej półsi, w metrach dla elipsoidy WGS-84. - */ - private static final double A = 6378137.0; - - /** - * double f: spłaszczenie elipsoidalne dla elipsoidy WGS-84, 1 / 298.257223563 - */ -// private static final double F = 1.0 / 298.257223563; - private static final double B = A * (1.0 - 1.0 / 298.257223563); - private static final double E_SQUARED = calculateESquared(A, B); - private static final double E2_SQUARED = calculateE2Squared(A, B); - private static final double TN = (A - B) / (A + B); - private static final double AP = A * (1.0 - TN + 5.0 * (TN * TN - TN * TN * TN) / 4.0 + 81.0 * TN * TN * TN * TN - TN * TN * TN * TN * TN / 64.0); - private static final double BP = 3.0 * A * (TN - TN * TN + 7.0 * (TN * TN * TN - TN * TN * TN * TN) / 8.0 + 55.0 * TN * TN * TN * TN * TN / 64.0) / 2.0; - private static final double CP = 15.0 * A * (TN * TN - TN * TN * TN + 3.0 * (TN * TN * TN * TN - TN * TN * TN * TN * TN) / 4.0) / 16.0; - private static final double DP = 35.0 * A * (TN * TN * TN - TN * TN * TN * TN + 11.0 * TN * TN * TN * TN * TN / 16.0) / 48.0; - private static final double EP = 315.0 * A * (TN * TN * TN * TN - TN * TN * TN * TN * TN) / 512.0; - - /** - * Współczynnik zniekształcenia skali mapy w południku osiowym dla odwzorowania kartograficznego PUWG-1992. - */ - private static final double OK = 0.9993; - private static final double OK_POW_2 = OK * OK; - private static final double OK_POW_3 = OK_POW_2 * OK; - private static final double OK_POW_4 = OK_POW_3 * OK; - private static final double OK_POW_5 = OK_POW_4 * OK; - private static final double OK_POW_6 = OK_POW_5 * OK; - private static final double OK_POW_7 = OK_POW_6 * OK; - private static final double OK_POW_8 = OK_POW_7 * OK; - } diff --git a/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/NMTData.java b/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/NMTData.java deleted file mode 100644 index df3224f..0000000 --- a/src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/NMTData.java +++ /dev/null @@ -1,21 +0,0 @@ -package pl.wat.ms4ds.terenfunkcje.konwersja; - -import pl.wat.ms4ds.terenfunkcje.Coord; - -/** - * - */ -public class NMTData { - - public Coord.Grid gridCoord; - - public double suma; - - public int licz; - - public NMTData(Coord.Grid gridCoord, double suma, int licz) { - this.gridCoord = gridCoord; - this.suma = suma; - this.licz = licz; - } -} diff --git a/src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTData.java b/src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTData.java new file mode 100644 index 0000000..3a931d1 --- /dev/null +++ b/src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTData.java @@ -0,0 +1,41 @@ +package pl.wat.ms4ds.terenfunkcje.nmt; + +/** + * + */ +public class NMTData { + + public int x; + public int y; + + public double sum; + + public int count; + + public double ell; + public double nll; + public double eur; + public double nur; + + + public NMTData(int x, int y, double sum, int count) { + this.x = x; + this.y = y; + this.sum = sum; + this.count = count; + } + + @Override + public String toString() { + return "NMTData{" + + "x=" + x + + ", y=" + y + + ", sum=" + sum + + ", count=" + count + + ", ell=" + ell + + ", nll=" + nll + + ", eur=" + eur + + ", nur=" + nur + + '}'; + } +} diff --git a/src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTReader.java b/src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTDataReader.java similarity index 52% rename from src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTReader.java rename to src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTDataReader.java index 4bd7c25..3efeebd 100644 --- a/src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTReader.java +++ b/src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTDataReader.java @@ -13,35 +13,38 @@ import java.util.Set; import java.util.zip.ZipEntry; import java.util.zip.ZipInputStream; -public class NMTReader { +public class NMTDataReader { - private static final Logger logger = LoggerFactory.getLogger(NMTReader.class); + private static final Logger logger = LoggerFactory.getLogger(NMTDataReader.class); static void main(String[] args) { // File dir = new File(System.getProperty("user.home") + "/nmt/gugik_SkorowidzNMT2018.gml"); - int i = 0; + HashMap nmtDataHashMap = new HashMap<>(); + String inDir = "C:/Workspace/nmt/gugik_1m/asc/"; String outDir = "C:/Workspace/nmt/unzipped/"; - String testFn = "D:\\Work\\73771_1025306_NMT-M348Dc41.xyz\\"; -// String testFn = "C:\\Workspace\\nmt\\M-33-7-A-c-3-2.asc"; +// String testFn = "D:\\Work\\73771_1025306_NMT-M348Dc41.xyz\\"; + String testFn = "D:\\Work\\M-33-7-A-c-3-2.asc"; + try { -// readFromFileASC(testFn); - readFromFileXYZ(testFn); + readFromFileASC(testFn, nmtDataHashMap); +// readFromFileXYZ(testFn, nmtDataHashMap); } catch (IOException e) { - throw new RuntimeException(e); + return; } // renameFiles(inDir, inDir); + Set files = NMTDataProvider.listFiles(inDir); for (String file : files) { try { String unzipfn = unzipFile(inDir + file, outDir); if (unzipfn.endsWith(".asc")) { - readFromFileASC(outDir + unzipfn); + readFromFileASC(outDir + unzipfn, nmtDataHashMap); } else if (unzipfn.endsWith(".xyz")) { - readFromFileXYZ(outDir + unzipfn); + readFromFileXYZ(outDir + unzipfn, nmtDataHashMap); } } catch (IOException e) { @@ -85,7 +88,7 @@ public class NMTReader { } } - private static void readFromFileASC(String fn) throws IOException { + private static void readFromFileASC(String fn, HashMap nmtDataHashMap) throws IOException { long start = System.currentTimeMillis(); File file = new File(fn); InputStream inputStream = new FileInputStream(file); @@ -98,10 +101,10 @@ public class NMTReader { int nrows = Integer.parseInt(split[1]); line = br.readLine(); split = line.split(" "); - double xll = Double.parseDouble(split[1]); + double xll_puwg = Double.parseDouble(split[1]); line = br.readLine(); split = line.split(" "); - double yll = Double.parseDouble(split[1]); + double yll_puwg = Double.parseDouble(split[1]); line = br.readLine(); split = line.split(" "); double cellsize = Double.parseDouble(split[1]); @@ -117,103 +120,81 @@ public class NMTReader { } } Coord.Geo geo_ll = new Coord.Geo(); - CoordUtils.convertPUWG1992ToWGS84(yll, xll, geo_ll); + Coord.convertPUWG1992ToWGS84(yll_puwg, xll_puwg, geo_ll); Coord.Geo geo_ur = new Coord.Geo(); - CoordUtils.convertPUWG1992ToWGS84(yll + nrows + cellsize, xll + ncols * cellsize, geo_ur); + Coord.convertPUWG1992ToWGS84(yll_puwg + nrows + cellsize, xll_puwg + ncols * cellsize, geo_ur); int d_x = (int) ((geo_ur.lon - geo_ll.lon) / MapConsts.DELTA_X) + 3; int d_y = (int) ((geo_ur.lat - geo_ll.lat) / MapConsts.DELTA_Y) + 3; - Square[][] kwadraty = new Square[d_x][d_y]; - final int x = Coord.zamienDlugoscGeoNaIdKwadratuX(geo_ll.lon); - final int y = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo_ll.lat); - kwadraty[0][0] = new Square(); - -// Kwadrat kw = Teren.getKwadrat(x, y); + final int x0 = Coord.zamienDlugoscGeoNaIdKwadratuX(geo_ll.lon); + final int y0 = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo_ll.lat); + NMTData nmtData = nmtDataHashMap.computeIfAbsent(new Coord.Grid(x0, y0), k -> new NMTData(x0, y0, 0, 0)); // Wyznacz współrzędne geo środka kwadratu. Coord.Geo geoCoord = new Coord.Geo(); - geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x); - geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y); + geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x0); + geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y0); PUWGCoord puwgCoord = new PUWGCoord(); // Wyznacz współrzędne PUWG lewego dolnego rogu kwadratu. - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat - MapConsts.DELTA_Y / 2, geoCoord.lon - MapConsts.DELTA_X / 2, puwgCoord); - kwadraty[0][0].ell = (int) puwgCoord.easting; - kwadraty[0][0].nll = (int) puwgCoord.northing; + Coord.convertWGS84ToPUWG1992(geoCoord.lat - MapConsts.DELTA_Y / 2, geoCoord.lon - MapConsts.DELTA_X / 2, puwgCoord); + nmtData.ell = (int) puwgCoord.easting; + nmtData.nll = (int) puwgCoord.northing; // Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord); - kwadraty[0][0].eur = (int) puwgCoord.easting; - kwadraty[0][0].nur = (int) puwgCoord.northing; - double dyy = 0; + Coord.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord); + nmtData.eur = (int) puwgCoord.easting; + nmtData.nur = (int) puwgCoord.northing; + double dy_puwg = 0; double h; - int idX = 0; - int idY = 0; -// HashMap kws = new HashMap(); + int x = x0; + int y = y0; for (int i = 0; i < nrows; i++) { - double yy = yll + dyy; - dyy += cellsize; + double y_puwg = yll_puwg + dy_puwg; + dy_puwg += cellsize; // Reset współrzędnej X na gridzie (siatce). - idX = 0; - if (yy >= kwadraty[idX][idY].nur) { + x = x0; + final int x3 = x; + final int y3 = y; + nmtData = nmtDataHashMap.computeIfAbsent(new Coord.Grid(x, y), k -> new NMTData(x3, y3, 0, 0)); + if (y_puwg >= nmtData.nur) { // Przekracza zakres współrzędnych pionowych, zatem kolejny/sąsiedni kwadrat po osi OY. - idY++; - if (kwadraty[idX][idY] == null) { - kwadraty[idX][idY] = new Square(); - } -// kw = Teren.getKwadrat(idX, idY); - if (kwadraty[idX][idY].nur == 0) { + y++; + final int x1 = x; + final int y1 = y; + nmtData = nmtDataHashMap.computeIfAbsent(new Coord.Grid(x, y), k -> new NMTData(x1, y1, 0, 0)); + if (nmtData.nur == 0) { // Świeży kwadrat. // Wyznacz współrzędne geo środka kwadratu. - geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(idX + x); - geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(idY + y); + geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x); + geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y); // Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord); - kwadraty[idX][idY].eur = (int) puwgCoord.easting; - kwadraty[idX][idY].nur = (int) puwgCoord.northing; + Coord.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord); + nmtData.eur = (int) puwgCoord.easting; + nmtData.nur = (int) puwgCoord.northing; } } - double dxx = 0; + double dx_puwg = 0; for (int j = 0; j < ncols; j++) { - double xx = xll + dxx; - dxx += cellsize; - if (xx >= kwadraty[idX][idY].eur) { + double x_puwg = xll_puwg + dx_puwg; + dx_puwg += cellsize; + h = data[i][j]; + if (x_puwg >= nmtData.eur) { // Przekracza zakres współrzędnych poziomych, zatem kolejny/sąsiedni kwadrat po osi OX. - idX++; - if (kwadraty[idX][idY] == null) { - kwadraty[idX][idY] = new Square(); - } -// kw = Teren.getKwadrat(idX, idY); - if (kwadraty[idX][idY].eur == 0) { + x++; + final int x2 = x; + final int y2 = y; + nmtData = nmtDataHashMap.computeIfAbsent(new Coord.Grid(x, y), k -> new NMTData(x2, y2, 0, 0)); + if (nmtData.eur == 0) { // Świeży kwadrat. // Wyznacz współrzędne geo środka kwadratu. - geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(idX + x); - geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(idY + y); + geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x); + geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y); // Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. - CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord); - kwadraty[idX][idY].eur = (int) puwgCoord.easting; - kwadraty[idX][idY].nur = (int) puwgCoord.northing; + Coord.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord); + nmtData.eur = (int) puwgCoord.easting; + nmtData.nur = (int) puwgCoord.northing; } - } else if (xx < kwadraty[idX][idY].ell) { - idX--; } - -// kws.put(kw, kw); - h = data[i][j]; - if (kwadraty[idX][idY] != Square.EMPTY && h > nodata) { - kwadraty[idX][idY].sumaWysokosci += h; - kwadraty[idX][idY].count++; - } - } - } - for (int i = 0; i < d_x; i++) { - for (int j = 0; j < d_y; j++) { - Square kw = kwadraty[i][j]; - if (kw != null) { - kw.wysokoscSrednia = (int) (kw.sumaWysokosci / kw.count); - Square kwOryg = Teren.getKwadrat(i + x, j + y); - // Aktualizacja tylko w przypadku braku danych o wysokości. - if (kwOryg != Square.EMPTY && kwOryg.wysokoscSrednia <= -1000) { -// kwOryg.wysokoscSrednia = kw.wysokoscSrednia; - kwOryg.sumaWysokosci = kw.sumaWysokosci; - kwOryg.count = kw.count; - } + if (h > nodata) { + nmtData.sum += h; + nmtData.count++; } } } @@ -221,59 +202,82 @@ public class NMTReader { } } - private static void readFromFileXYZ(String fn) throws IOException { + private static void readFromFileXYZ(String fn, HashMap nmtDataHashMap) throws IOException { File file = new File(fn); InputStream inputStream = new FileInputStream(file); - Square kw = Square.EMPTY; - PUWGCoord puwgCoord = new PUWGCoord(); - Coord.Geo geo = new Coord.Geo(); - HashMap map = new HashMap<>(); try (BufferedReader br = new BufferedReader(new InputStreamReader(inputStream))) { - String line; - while ((line = br.readLine()) != null) { - String[] split = line.split(" "); - double x = Double.parseDouble(split[0]); - double y = Double.parseDouble(split[1]); - double h = Double.parseDouble(split[2]); - if (kw.ell > x || kw.eur < x || kw.nll > y || kw.nur < y) { - // Punkt poza granicą bieżącego kwadratu. - CoordUtils.convertPUWG1992ToWGS84(y, x, geo); - kw = Teren.getKwadrat(geo.lat, geo.lon); - if (kw == Square.EMPTY) { - continue; - } - if (kw.stopienZalesienia > 0 || kw.stopienZawodnienia > 0 || kw.stopienZabudowy > 0) { - System.out.println(kw); - } - map.put(kw, kw); - if (kw.nur == 0) { - // Kwadrat jeszcze nie był odczytany (czysty). - int idX = Coord.zamienDlugoscGeoNaIdKwadratuX(geo.lon); - int idY = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo.lat); - // Współrzędne geo środka kwadratu. - geo.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(idX); - geo.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(idY); - // Wyznacz współrzędne PUWG lewego dolnego rogu kwadratu. - CoordUtils.convertWGS84ToPUWG1992(geo.lat - MapConsts.DELTA_Y / 2, geo.lon - MapConsts.DELTA_X / 2, puwgCoord); - kw.ell = (int) puwgCoord.easting; - kw.nll = (int) puwgCoord.northing; - // Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. - CoordUtils.convertWGS84ToPUWG1992(geo.lat + MapConsts.DELTA_Y / 2, geo.lon + MapConsts.DELTA_X / 2, puwgCoord); - kw.eur = (int) puwgCoord.easting; - kw.nur = (int) puwgCoord.northing; - } - } -// if (kw != Kwadrat.EMPTY_SQUARE && kw.wysokoscSrednia <= -1000) { - if (kw != Square.EMPTY) { - kw.sumaWysokosci += h; - kw.count++; -// kw.wysokoscSrednia = kw.wysokoscSrednia; - } + PUWGCoord puwgCoord = new PUWGCoord(); + Coord.Geo geo = new Coord.Geo(); + String line = br.readLine(); + if (line == null) { + return; + } + String[] split = line.split(" "); + if (split.length != 3) { + return; + } + double x_puwg = Double.parseDouble(split[0]); + double y_puwg = Double.parseDouble(split[1]); + double h = Double.parseDouble(split[2]); + Coord.convertPUWG1992ToWGS84(y_puwg, x_puwg, geo); + int x = Coord.zamienDlugoscGeoNaIdKwadratuX(geo.lon); + int y = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo.lat); + final int xx1 = x; + final int yy1 = y; + Coord.Grid coordGrid = new Coord.Grid(x, y); + NMTData nmtData = nmtDataHashMap.computeIfAbsent(coordGrid, k -> new NMTData(xx1, yy1, 0, 0)); + if (nmtData.nur == 0) { + // Kwadrat jeszcze nie był odczytany (czysty). + // Współrzędne geo środka kwadratu. + geo.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x); + geo.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y); + // Wyznacz współrzędne PUWG lewego dolnego rogu kwadratu. + Coord.convertWGS84ToPUWG1992(geo.lat - MapConsts.DELTA_Y / 2, geo.lon - MapConsts.DELTA_X / 2, puwgCoord); + nmtData.ell = (int) puwgCoord.easting; + nmtData.nll = (int) puwgCoord.northing; + // Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. + Coord.convertWGS84ToPUWG1992(geo.lat + MapConsts.DELTA_Y / 2, geo.lon + MapConsts.DELTA_X / 2, puwgCoord); + nmtData.eur = (int) puwgCoord.easting; + nmtData.nur = (int) puwgCoord.northing; + } + nmtData.sum += h; + nmtData.count++; + while ((line = br.readLine()) != null) { + split = line.split(" "); + if (split.length != 3) { + continue; + } + x_puwg = Double.parseDouble(split[0]); + y_puwg = Double.parseDouble(split[1]); + h = Double.parseDouble(split[2]); + if (nmtData.ell > x_puwg || nmtData.eur < x_puwg || nmtData.nll > y_puwg || nmtData.nur < y_puwg) { + // Punkt poza granicą bieżącego kwadratu. + Coord.convertPUWG1992ToWGS84(y_puwg, x_puwg, geo); + x = Coord.zamienDlugoscGeoNaIdKwadratuX(geo.lon); + y = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo.lat); + coordGrid.set(x, y); + final int xx2 = x; + final int yy2 = y; + nmtData = nmtDataHashMap.computeIfAbsent(coordGrid, k -> new NMTData(xx2, yy2, 0, 0)); + if (nmtData.nur == 0) { + // Kwadrat jeszcze nie był odczytany (czysty). + // Współrzędne geo środka kwadratu. + geo.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x); + geo.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y); + // Wyznacz współrzędne PUWG lewego dolnego rogu kwadratu. + Coord.convertWGS84ToPUWG1992(geo.lat - MapConsts.DELTA_Y / 2, geo.lon - MapConsts.DELTA_X / 2, puwgCoord); + nmtData.ell = (int) puwgCoord.easting; + nmtData.nll = (int) puwgCoord.northing; + // Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. + Coord.convertWGS84ToPUWG1992(geo.lat + MapConsts.DELTA_Y / 2, geo.lon + MapConsts.DELTA_X / 2, puwgCoord); + nmtData.eur = (int) puwgCoord.easting; + nmtData.nur = (int) puwgCoord.northing; + } + } + nmtData.sum += h; + nmtData.count++; } - - } - } public static String unzipFile(String zipFileName, String destDir) throws IOException {