jrulka/terrain-utilities
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Czytanie danych nmt cd.

Browse Source
This commit is contained in:
Jarosław Rulka
2026-01-14 20:18:28 +01:00
parent f4d2d60286
commit 1f4c5b52f2
8 changed files with 413 additions and 404 deletions
Download Patch File Download Diff File Expand all files Collapse all files

260
src/main/java/pl/wat/ms4ds/terenfunkcje/Coord.java
View File

@@ -1,7 +1,13 @@
package pl.wat.ms4ds.terenfunkcje; package pl.wat.ms4ds.terenfunkcje;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import pl.wat.ms4ds.terenfunkcje.konwersja.PUWGCoord;
public class Coord { public class Coord {
private static final Logger logger = LoggerFactory.getLogger(Coord.class);
public static class Geo { public static class Geo {
public double lat; public double lat;
public double lon; public double lon;
@@ -18,6 +24,14 @@ public class Coord {
lat = other.lat; lat = other.lat;
lon = other.lon; lon = other.lon;
} }
@Override
public String toString() {
return "Geo{" +
"lat=" + lat +
", lon=" + lon +
'}';
}
} }
/** /**
@@ -45,6 +59,14 @@ public class Coord {
this.easting = easting; this.easting = easting;
this.northing = northing; this.northing = northing;
} }
@Override
public String toString() {
return "Puwg{" +
"easting=" + easting +
", northing=" + northing +
'}';
}
} }
public static class Grid { public static class Grid {
@@ -74,15 +96,12 @@ public class Coord {
y = zamienWspYmsNaIdKwadratuY(yms); y = zamienWspYmsNaIdKwadratuY(yms);
} }
@Override
public String toString() { public String toString() {
StringBuilder sb = new StringBuilder("(x= "); return "Grid{" +
sb.append(x); "x=" + x +
sb.append(", y= "); ", y=" + y +
sb.append(y); '}';
sb.append(')');
// String s = "(x= " + Integer.toString(x) + ", y= " + Integer.toString(y) + ")";
return sb.toString();
} }
private static final double ODWROT_SS_DX_MS = 1.0 / MapConsts.SS_DX_MS; private static final double ODWROT_SS_DX_MS = 1.0 / MapConsts.SS_DX_MS;
@@ -311,43 +330,198 @@ public class Coord {
return (int) yy; 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_DX_MS = 1.0 / MapConsts.SS_DX_MS;
private static final double ODWROT_SS_DY_MS = 1.0 / MapConsts.SS_DY_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.
* <p>
* 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;
} }

5
src/main/java/pl/wat/ms4ds/terenfunkcje/Square.java
View File

@@ -213,8 +213,8 @@ public class Square {
} }
} }
final int x; public final int x;
final int y; public final int y;
public double sumaWysokosci; public double sumaWysokosci;
public int count = 1; public int count = 1;
public double ell; public double ell;
@@ -238,7 +238,6 @@ public class Square {
return result; return result;
} }
public String toString() { public String toString() {
StringBuilder linia = new StringBuilder(100); StringBuilder linia = new StringBuilder(100);
linia.append("["); linia.append("[");

2
src/main/java/pl/wat/ms4ds/terenfunkcje/Teren.java
View File

@@ -259,7 +259,7 @@ public class Teren {
public static Square getKwadratPUWG(double northing, double easting) { public static Square getKwadratPUWG(double northing, double easting) {
Coord.Geo geoCoord = new Coord.Geo(); Coord.Geo geoCoord = new Coord.Geo();
CoordUtils.convertPUWG1992ToWGS84(northing, easting, geoCoord); Coord.convertPUWG1992ToWGS84(northing, easting, geoCoord);
return getKwadrat(geoCoord.lat, geoCoord.lon); return getKwadrat(geoCoord.lat, geoCoord.lon);
} }

18
src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordTest.java
View File

@@ -18,20 +18,20 @@ public class CoordTest {
Coord.Geo geoCoord = new Coord.Geo(); Coord.Geo geoCoord = new Coord.Geo();
geoCoord.lon = 19; geoCoord.lon = 19;
geoCoord.lat = 50; 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); logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing);
double e = puwgCoord.easting; double e = puwgCoord.easting;
double n = puwgCoord.northing; double n = puwgCoord.northing;
geoCoord.lon = 20; geoCoord.lon = 20;
geoCoord.lat = 51; geoCoord.lat = 51;
CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord); Coord.convertWGS84ToPUWG1992(geoCoord.lat, geoCoord.lon, puwgCoord);
double dx = puwgCoord.easting - e; double dx = puwgCoord.easting - e;
double dy = puwgCoord.northing - n; double dy = puwgCoord.northing - n;
logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing); 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); logger.debug("PUWG (e,n)=({}, {}) => Lat={}, Lon={}", puwgCoord.easting, puwgCoord.northing, geoCoord.lat, geoCoord.lon);
// coord.proj = 1; // coord.proj = 1;
@@ -52,24 +52,24 @@ public class CoordTest {
puwgCoord.easting = 500000; puwgCoord.easting = 500000;
puwgCoord.northing = 500000; puwgCoord.northing = 500000;
logger.debug("----------------------------------"); 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); 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); logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing);
puwgCoord.easting = 100000; puwgCoord.easting = 100000;
puwgCoord.northing = 470642; 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); 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); logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing);
puwgCoord.easting = 821310; puwgCoord.easting = 821310;
puwgCoord.northing = 369750; 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); 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); logger.debug("Lat={}, Lon={} => PUWG (e,n)=({}, {})", geoCoord.lat, geoCoord.lon, puwgCoord.easting, puwgCoord.northing);

198
src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/CoordUtils.java
View File

@@ -4,6 +4,7 @@ import pl.wat.ms4ds.terenfunkcje.Coord;
import pl.wat.ms4ds.terenfunkcje.Teren; import pl.wat.ms4ds.terenfunkcje.Teren;
import org.slf4j.Logger; import org.slf4j.Logger;
import org.slf4j.LoggerFactory; import org.slf4j.LoggerFactory;
import pl.wat.ms4ds.terenfunkcje.nmt.NMTData;
import java.io.BufferedReader; import java.io.BufferedReader;
import java.io.FileReader; import java.io.FileReader;
@@ -103,15 +104,15 @@ public class CoordUtils {
} catch (NumberFormatException e) { } catch (NumberFormatException e) {
logger.warn("Bledne dane w pliku: " + fileName); 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); Coord.Grid idKw = new Coord.Grid(latLon.lon, latLon.lat);
NMTData daneWysok = daneWysokHashMap.get(idKw); NMTData daneWysok = daneWysokHashMap.get(idKw);
if (daneWysok == null) { if (daneWysok == null) {
daneWysok = new NMTData(idKw, wysokosc, 1); daneWysok = new NMTData(idKw.x, idKw.y, wysokosc, 1);
daneWysokHashMap.put(idKw, daneWysok); daneWysokHashMap.put(idKw, daneWysok);
} else { } else {
daneWysok.suma += wysokosc; daneWysok.sum += wysokosc;
daneWysok.licz++; daneWysok.count++;
} }
line = br.readLine(); line = br.readLine();
if (m++ % 100000 == 0) { 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.
* <p>
* 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;
} }

21
src/main/java/pl/wat/ms4ds/terenfunkcje/konwersja/NMTData.java
View File

@@ -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;
}
}

41
src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTData.java Normal file
View File

@@ -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 +
'}';
}
}

272
src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTReader.java → src/main/java/pl/wat/ms4ds/terenfunkcje/nmt/NMTDataReader.java
View File

@@ -13,35 +13,38 @@ import java.util.Set;
import java.util.zip.ZipEntry; import java.util.zip.ZipEntry;
import java.util.zip.ZipInputStream; 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) { static void main(String[] args) {
// File dir = new File(System.getProperty("user.home") + "/nmt/gugik_SkorowidzNMT2018.gml"); // File dir = new File(System.getProperty("user.home") + "/nmt/gugik_SkorowidzNMT2018.gml");
int i = 0; HashMap<Coord.Grid, NMTData> nmtDataHashMap = new HashMap<>();
String inDir = "C:/Workspace/nmt/gugik_1m/asc/"; String inDir = "C:/Workspace/nmt/gugik_1m/asc/";
String outDir = "C:/Workspace/nmt/unzipped/"; String outDir = "C:/Workspace/nmt/unzipped/";
String testFn = "D:\\Work\\73771_1025306_NMT-M348Dc41.xyz\\"; // 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\\M-33-7-A-c-3-2.asc";
try { try {
// readFromFileASC(testFn); readFromFileASC(testFn, nmtDataHashMap);
readFromFileXYZ(testFn); // readFromFileXYZ(testFn, nmtDataHashMap);
} catch (IOException e) { } catch (IOException e) {
throw new RuntimeException(e); return;
} }
// renameFiles(inDir, inDir); // renameFiles(inDir, inDir);
Set<String> files = NMTDataProvider.listFiles(inDir); Set<String> files = NMTDataProvider.listFiles(inDir);
for (String file : files) { for (String file : files) {
try { try {
String unzipfn = unzipFile(inDir + file, outDir); String unzipfn = unzipFile(inDir + file, outDir);
if (unzipfn.endsWith(".asc")) { if (unzipfn.endsWith(".asc")) {
readFromFileASC(outDir + unzipfn); readFromFileASC(outDir + unzipfn, nmtDataHashMap);
} else if (unzipfn.endsWith(".xyz")) { } else if (unzipfn.endsWith(".xyz")) {
readFromFileXYZ(outDir + unzipfn); readFromFileXYZ(outDir + unzipfn, nmtDataHashMap);
} }
} catch (IOException e) { } 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<Coord.Grid, NMTData> nmtDataHashMap) throws IOException {
long start = System.currentTimeMillis(); long start = System.currentTimeMillis();
File file = new File(fn); File file = new File(fn);
InputStream inputStream = new FileInputStream(file); InputStream inputStream = new FileInputStream(file);
@@ -98,10 +101,10 @@ public class NMTReader {
int nrows = Integer.parseInt(split[1]); int nrows = Integer.parseInt(split[1]);
line = br.readLine(); line = br.readLine();
split = line.split(" "); split = line.split(" ");
double xll = Double.parseDouble(split[1]); double xll_puwg = Double.parseDouble(split[1]);
line = br.readLine(); line = br.readLine();
split = line.split(" "); split = line.split(" ");
double yll = Double.parseDouble(split[1]); double yll_puwg = Double.parseDouble(split[1]);
line = br.readLine(); line = br.readLine();
split = line.split(" "); split = line.split(" ");
double cellsize = Double.parseDouble(split[1]); double cellsize = Double.parseDouble(split[1]);
@@ -117,103 +120,81 @@ public class NMTReader {
} }
} }
Coord.Geo geo_ll = new Coord.Geo(); 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(); 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_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; int d_y = (int) ((geo_ur.lat - geo_ll.lat) / MapConsts.DELTA_Y) + 3;
Square[][] kwadraty = new Square[d_x][d_y]; final int x0 = Coord.zamienDlugoscGeoNaIdKwadratuX(geo_ll.lon);
final int x = Coord.zamienDlugoscGeoNaIdKwadratuX(geo_ll.lon); final int y0 = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo_ll.lat);
final int y = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo_ll.lat); NMTData nmtData = nmtDataHashMap.computeIfAbsent(new Coord.Grid(x0, y0), k -> new NMTData(x0, y0, 0, 0));
kwadraty[0][0] = new Square();
// Kwadrat kw = Teren.getKwadrat(x, y);
// Wyznacz współrzędne geo środka kwadratu. // Wyznacz współrzędne geo środka kwadratu.
Coord.Geo geoCoord = new Coord.Geo(); Coord.Geo geoCoord = new Coord.Geo();
geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x); geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x0);
geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y); geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y0);
PUWGCoord puwgCoord = new PUWGCoord(); PUWGCoord puwgCoord = new PUWGCoord();
// Wyznacz współrzędne PUWG lewego dolnego rogu kwadratu. // Wyznacz współrzędne PUWG lewego dolnego rogu kwadratu.
CoordUtils.convertWGS84ToPUWG1992(geoCoord.lat - MapConsts.DELTA_Y / 2, geoCoord.lon - MapConsts.DELTA_X / 2, puwgCoord); Coord.convertWGS84ToPUWG1992(geoCoord.lat - MapConsts.DELTA_Y / 2, geoCoord.lon - MapConsts.DELTA_X / 2, puwgCoord);
kwadraty[0][0].ell = (int) puwgCoord.easting; nmtData.ell = (int) puwgCoord.easting;
kwadraty[0][0].nll = (int) puwgCoord.northing; nmtData.nll = (int) puwgCoord.northing;
// Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. // 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); Coord.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord);
kwadraty[0][0].eur = (int) puwgCoord.easting; nmtData.eur = (int) puwgCoord.easting;
kwadraty[0][0].nur = (int) puwgCoord.northing; nmtData.nur = (int) puwgCoord.northing;
double dyy = 0; double dy_puwg = 0;
double h; double h;
int idX = 0; int x = x0;
int idY = 0; int y = y0;
// HashMap<Kwadrat, Kwadrat> kws = new HashMap();
for (int i = 0; i < nrows; i++) { for (int i = 0; i < nrows; i++) {
double yy = yll + dyy; double y_puwg = yll_puwg + dy_puwg;
dyy += cellsize; dy_puwg += cellsize;
// Reset współrzędnej X na gridzie (siatce). // Reset współrzędnej X na gridzie (siatce).
idX = 0; x = x0;
if (yy >= kwadraty[idX][idY].nur) { 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. // Przekracza zakres współrzędnych pionowych, zatem kolejny/sąsiedni kwadrat po osi OY.
idY++; y++;
if (kwadraty[idX][idY] == null) { final int x1 = x;
kwadraty[idX][idY] = new Square(); final int y1 = y;
} nmtData = nmtDataHashMap.computeIfAbsent(new Coord.Grid(x, y), k -> new NMTData(x1, y1, 0, 0));
// kw = Teren.getKwadrat(idX, idY); if (nmtData.nur == 0) {
if (kwadraty[idX][idY].nur == 0) {
// Świeży kwadrat. // Świeży kwadrat.
// Wyznacz współrzędne geo środka kwadratu. // Wyznacz współrzędne geo środka kwadratu.
geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(idX + x); geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x);
geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(idY + y); geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y);
// Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. // 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); Coord.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord);
kwadraty[idX][idY].eur = (int) puwgCoord.easting; nmtData.eur = (int) puwgCoord.easting;
kwadraty[idX][idY].nur = (int) puwgCoord.northing; nmtData.nur = (int) puwgCoord.northing;
} }
} }
double dxx = 0; double dx_puwg = 0;
for (int j = 0; j < ncols; j++) { for (int j = 0; j < ncols; j++) {
double xx = xll + dxx; double x_puwg = xll_puwg + dx_puwg;
dxx += cellsize; dx_puwg += cellsize;
if (xx >= kwadraty[idX][idY].eur) { h = data[i][j];
if (x_puwg >= nmtData.eur) {
// Przekracza zakres współrzędnych poziomych, zatem kolejny/sąsiedni kwadrat po osi OX. // Przekracza zakres współrzędnych poziomych, zatem kolejny/sąsiedni kwadrat po osi OX.
idX++; x++;
if (kwadraty[idX][idY] == null) { final int x2 = x;
kwadraty[idX][idY] = new Square(); final int y2 = y;
} nmtData = nmtDataHashMap.computeIfAbsent(new Coord.Grid(x, y), k -> new NMTData(x2, y2, 0, 0));
// kw = Teren.getKwadrat(idX, idY); if (nmtData.eur == 0) {
if (kwadraty[idX][idY].eur == 0) {
// Świeży kwadrat. // Świeży kwadrat.
// Wyznacz współrzędne geo środka kwadratu. // Wyznacz współrzędne geo środka kwadratu.
geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(idX + x); geoCoord.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x);
geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(idY + y); geoCoord.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y);
// Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. // 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); Coord.convertWGS84ToPUWG1992(geoCoord.lat + MapConsts.DELTA_Y / 2, geoCoord.lon + MapConsts.DELTA_X / 2, puwgCoord);
kwadraty[idX][idY].eur = (int) puwgCoord.easting; nmtData.eur = (int) puwgCoord.easting;
kwadraty[idX][idY].nur = (int) puwgCoord.northing; nmtData.nur = (int) puwgCoord.northing;
} }
} else if (xx < kwadraty[idX][idY].ell) {
idX--;
} }
if (h > nodata) {
// kws.put(kw, kw); nmtData.sum += h;
h = data[i][j]; nmtData.count++;
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;
}
} }
} }
} }
@@ -221,59 +202,82 @@ public class NMTReader {
} }
} }
private static void readFromFileXYZ(String fn) throws IOException { private static void readFromFileXYZ(String fn, HashMap<Coord.Grid, NMTData> nmtDataHashMap) throws IOException {
File file = new File(fn); File file = new File(fn);
InputStream inputStream = new FileInputStream(file); InputStream inputStream = new FileInputStream(file);
Square kw = Square.EMPTY;
PUWGCoord puwgCoord = new PUWGCoord();
Coord.Geo geo = new Coord.Geo();
HashMap<Square, Square> map = new HashMap<>();
try (BufferedReader br = new BufferedReader(new InputStreamReader(inputStream))) { try (BufferedReader br = new BufferedReader(new InputStreamReader(inputStream))) {
String line; PUWGCoord puwgCoord = new PUWGCoord();
while ((line = br.readLine()) != null) { Coord.Geo geo = new Coord.Geo();
String[] split = line.split(" "); String line = br.readLine();
double x = Double.parseDouble(split[0]); if (line == null) {
double y = Double.parseDouble(split[1]); return;
double h = Double.parseDouble(split[2]); }
if (kw.ell > x || kw.eur < x || kw.nll > y || kw.nur < y) { String[] split = line.split(" ");
// Punkt poza granicą bieżącego kwadratu. if (split.length != 3) {
CoordUtils.convertPUWG1992ToWGS84(y, x, geo); return;
kw = Teren.getKwadrat(geo.lat, geo.lon); }
if (kw == Square.EMPTY) { double x_puwg = Double.parseDouble(split[0]);
continue; double y_puwg = Double.parseDouble(split[1]);
} double h = Double.parseDouble(split[2]);
if (kw.stopienZalesienia > 0 || kw.stopienZawodnienia > 0 || kw.stopienZabudowy > 0) { Coord.convertPUWG1992ToWGS84(y_puwg, x_puwg, geo);
System.out.println(kw); int x = Coord.zamienDlugoscGeoNaIdKwadratuX(geo.lon);
} int y = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo.lat);
map.put(kw, kw); final int xx1 = x;
if (kw.nur == 0) { final int yy1 = y;
// Kwadrat jeszcze nie był odczytany (czysty). Coord.Grid coordGrid = new Coord.Grid(x, y);
int idX = Coord.zamienDlugoscGeoNaIdKwadratuX(geo.lon); NMTData nmtData = nmtDataHashMap.computeIfAbsent(coordGrid, k -> new NMTData(xx1, yy1, 0, 0));
int idY = Coord.zamienSzerokoscGeoNaIdKwadratuY(geo.lat); if (nmtData.nur == 0) {
// Współrzędne geo środka kwadratu. // Kwadrat jeszcze nie był odczytany (czysty).
geo.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(idX); // Współrzędne geo środka kwadratu.
geo.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(idY); geo.lon = Coord.zamienIdKwadratuXNaDlugoscGeo(x);
// Wyznacz współrzędne PUWG lewego dolnego rogu kwadratu. geo.lat = Coord.zamienIdKwadratuYNaSzerokoscGeo(y);
CoordUtils.convertWGS84ToPUWG1992(geo.lat - MapConsts.DELTA_Y / 2, geo.lon - MapConsts.DELTA_X / 2, puwgCoord); // Wyznacz współrzędne PUWG lewego dolnego rogu kwadratu.
kw.ell = (int) puwgCoord.easting; Coord.convertWGS84ToPUWG1992(geo.lat - MapConsts.DELTA_Y / 2, geo.lon - MapConsts.DELTA_X / 2, puwgCoord);
kw.nll = (int) puwgCoord.northing; nmtData.ell = (int) puwgCoord.easting;
// Wyznacz współrzędne PUWG prawego górnego rogu kwadratu. nmtData.nll = (int) puwgCoord.northing;
CoordUtils.convertWGS84ToPUWG1992(geo.lat + MapConsts.DELTA_Y / 2, geo.lon + MapConsts.DELTA_X / 2, puwgCoord); // Wyznacz współrzędne PUWG prawego górnego rogu kwadratu.
kw.eur = (int) puwgCoord.easting; Coord.convertWGS84ToPUWG1992(geo.lat + MapConsts.DELTA_Y / 2, geo.lon + MapConsts.DELTA_X / 2, puwgCoord);
kw.nur = (int) puwgCoord.northing; nmtData.eur = (int) puwgCoord.easting;
} nmtData.nur = (int) puwgCoord.northing;
} }
// if (kw != Kwadrat.EMPTY_SQUARE && kw.wysokoscSrednia <= -1000) { nmtData.sum += h;
if (kw != Square.EMPTY) { nmtData.count++;
kw.sumaWysokosci += h; while ((line = br.readLine()) != null) {
kw.count++; split = line.split(" ");
// kw.wysokoscSrednia = kw.wysokoscSrednia; 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 { public static String unzipFile(String zipFileName, String destDir) throws IOException {