diff --git a/docs/NTG_71.pdf b/docs/NTG_71.pdf new file mode 100644 index 0000000..06b5b95 Binary files /dev/null and b/docs/NTG_71.pdf differ diff --git a/docs/NTG_80.pdf b/docs/NTG_80.pdf new file mode 100644 index 0000000..2d654df Binary files /dev/null and b/docs/NTG_80.pdf differ diff --git a/docs/TransformationsCoordonneesGeodesiques.pdf b/docs/TransformationsCoordonneesGeodesiques.pdf new file mode 100644 index 0000000..5ba8de6 Binary files /dev/null and b/docs/TransformationsCoordonneesGeodesiques.pdf differ diff --git a/src/main/java/net/yageek/lambert/Lambert.java b/src/main/java/net/yageek/lambert/Lambert.java index 8f0739b..289be90 100644 --- a/src/main/java/net/yageek/lambert/Lambert.java +++ b/src/main/java/net/yageek/lambert/Lambert.java @@ -1,38 +1,137 @@ package net.yageek.lambert; +import sun.reflect.generics.reflectiveObjects.NotImplementedException; + import static java.lang.Math.*; import static net.yageek.lambert.LambertZone.*; -public class Lambert { +/* +https://github.com/yageek/lambert-java +https://bintray.com/yageek/maven/lambert-java/view/files/net/yageek/lambert/lambert-java/1.1 + +Online samples : +http://geofree.fr/gf/coordinateConv.asp#listSys + +-------------------------------------------------------------------------------------- +Install cs2cs on Ubuntu : +http://www.sarasafavi.com/installing-gdalogr-on-ubuntu.html + +-------------------------------------------------------------------------------------- +http://cs2cs.mygeodata.eu/ +Conversion From Lambert Zone II to WGS 84 : +$>cs2cs +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs +to +proj=longlat +datum=WGS84 +no_defs -f "%.11f" < 618115 2430676 +> EOF + +2.58331732871 48.87414278182 43.05512374267 + +-------------------------------------------------------------------------------------- +Conversion From WGS 84 To Lambert Zone II: +$>cs2cs +proj=longlat +datum=WGS84 +no_defs +to +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs -f "%.11f" < eps) - { + while (delta > eps) { phi0 = phiI; - phiI = 2*atan(pow((1+e*sin(phi0))/(1-e*sin(phi0)),e/2.0)*exp(latISo)) - M_PI_2; + phiI = 2 * atan(pow((1 + e * sin(phi0)) / (1 - e * sin(phi0)), e / 2d) * exp(latISo)) - M_PI_2; delta = abs(phiI - phi0); } return phiI; + } + + + /* + * ALGO0003 + */ + public static LambertPoint geographicToLambertAlg003(double latitude, double longitude, LambertZone zone, double lonMeridian, double e) { + + double n = zone.n(); + double C = zone.c(); + double xs = zone.xs(); + double ys = zone.ys(); + + double latIso = latitudeISOFromLat(latitude, e); + + double eLatIso = exp(-n * latIso); + + double nLon = n * (longitude - lonMeridian); + + double x = xs + C * eLatIso * sin(nLon); + double y = ys - C * eLatIso * cos(nLon); + + return new LambertPoint(x, y, 0); + } + + /* + * http://geodesie.ign.fr/contenu/fichiers/documentation/pedagogiques/TransformationsCoordonneesGeodesiques.pdf + * 3.4 Coordonnées géographiques Lambert + */ + public static LambertPoint geographicToLambert(double latitude, double longitude, LambertZone zone, double lonMeridian, double e) { + + double n = zone.n(); + double C = zone.c(); + double xs = zone.xs(); + double ys = zone.ys(); + + double sinLat = sin(latitude); + double eSinLat = (e * sinLat); + double elt1 = (1 + sinLat) / (1 - sinLat); + double elt2 = (1 + eSinLat) / (1 - eSinLat); + + double latIso = (1 / 2d) * log(elt1) - (e / 2d) * log(elt2); + double R = C * exp(-(n * latIso)); + + double LAMBDA = n * (longitude - lonMeridian); + + double x = xs + (R * sin(LAMBDA)); + double y = ys - (R * cos(LAMBDA)); + + return new LambertPoint(x, y, 0); } /* * ALGO0004 - Lambert vers geographiques */ - private static LambertPoint lambertToGeographic(LambertPoint org, LambertZone zone, double lonMeridian, double e, double eps) - { + public static LambertPoint lambertToGeographic(LambertPoint org, LambertZone zone, double lonMeridian, double e, double eps) { double n = zone.n(); double C = zone.c(); double xs = zone.xs(); @@ -46,16 +145,15 @@ private static LambertPoint lambertToGeographic(LambertPoint org, LambertZone zo R = sqrt((x - xs) * (x - xs) + (y - ys) * (y - ys)); - gamma = atan((x-xs)/(ys-y)); + gamma = atan((x - xs) / (ys - y)); - lon = lonMeridian + gamma/n; + lon = lonMeridian + gamma / n; - latIso = -1/n*log(abs(R/C)); + latIso = -1 / n * log(abs(R / C)); double lat = latitudeFromLatitudeISO(latIso, e, eps); - LambertPoint dest = new LambertPoint(lon,lat,0); - return dest; + return new LambertPoint(lon, lat, 0); } /* @@ -63,10 +161,9 @@ private static LambertPoint lambertToGeographic(LambertPoint org, LambertZone zo * */ - private static double lambertNormal(double lat, double a, double e) - { + private static double lambertNormal(double lat, double a, double e) { - return a/sqrt(1-e*e*sin(lat)*sin(lat)); + return a / sqrt(1 - e * e * sin(lat) * sin(lat)); } /* @@ -74,15 +171,14 @@ private static double lambertNormal(double lat, double a, double e) * */ - private static LambertPoint geographicToCartesian(double lon, double lat, double he, double a, double e) - { + private static LambertPoint geographicToCartesian(double lon, double lat, double he, double a, double e) { double N = lambertNormal(lat, a, e); - LambertPoint pt = new LambertPoint(0,0,0); + LambertPoint pt = new LambertPoint(0, 0, 0); - pt.setX((N+he)*cos(lat)*cos(lon)); - pt.setY((N+he)*cos(lat)*sin(lon)); - pt.setZ((N*(1-e*e)+he)*sin(lat)); + pt.setX((N + he) * cos(lat) * cos(lon)); + pt.setY((N + he) * cos(lat) * sin(lon)); + pt.setZ((N * (1 - e * e) + he) * sin(lat)); return pt; @@ -92,64 +188,98 @@ private static LambertPoint geographicToCartesian(double lon, double lat, double * ALGO0012 - Passage des coordonnées cartésiennes aux coordonnées géographiques */ - private static LambertPoint cartesianToGeographic(LambertPoint org, double meridien, double a, double e, double eps) - { + private static LambertPoint cartesianToGeographic(LambertPoint org, double meridien, double a, double e, double eps) { double x = org.getX(), y = org.getY(), z = org.getZ(); - double lon = meridien + atan(y/x); + double lon = meridien + atan(y / x); - double module = sqrt(x*x + y*y); + double module = sqrt(x * x + y * y); - double phi0 = atan(z/(module*(1-(a*e*e)/sqrt(x*x+y*y+z*z)))); - double phiI = atan(z/module/(1-a*e*e*cos(phi0)/(module * sqrt(1-e*e*sin(phi0)*sin(phi0))))); - double delta= abs(phiI - phi0); - while(delta > eps) - { + double phi0 = atan(z / (module * (1 - (a * e * e) / sqrt(x * x + y * y + z * z)))); + double phiI = atan(z / module / (1 - a * e * e * cos(phi0) / (module * sqrt(1 - e * e * sin(phi0) * sin(phi0))))); + double delta = abs(phiI - phi0); + while (delta > eps) { phi0 = phiI; - phiI = atan(z/module/(1-a*e*e*cos(phi0)/(module * sqrt(1-e*e*sin(phi0)*sin(phi0))))); - delta= abs(phiI - phi0); + phiI = atan(z / module / (1 - a * e * e * cos(phi0) / (module * sqrt(1 - e * e * sin(phi0) * sin(phi0))))); + delta = abs(phiI - phi0); } - double he = module/cos(phiI) - a/sqrt(1-e*e*sin(phiI)*sin(phiI)); - - LambertPoint pt = new LambertPoint(lon,phiI,he); + double he = module / cos(phiI) - a / sqrt(1 - e * e * sin(phiI) * sin(phiI)); - return pt; + return new LambertPoint(lon, phiI, he); } + /* * Convert Lambert -> WGS84 * http://geodesie.ign.fr/contenu/fichiers/documentation/pedagogiques/transfo.pdf * */ - public static LambertPoint convertToWGS84(LambertPoint org, LambertZone zone){ + public static LambertPoint convertToWGS84(LambertPoint org, LambertZone zone) { - if(zone == Lambert93) - { - return lambertToGeographic(org,Lambert93,LON_MERID_IERS,E_WGS84,DEFAULT_EPS); - } - else { - LambertPoint pt1 = lambertToGeographic(org, zone, LON_MERID_PARIS, E_CLARK_IGN, DEFAULT_EPS); + if (zone == Lambert93) { + return lambertToGeographic(org, Lambert93, LON_MERID_IERS, E_WGS84, DEFAULT_EPS); + } else { + LambertPoint pt1 = lambertToGeographic(org, zone, LON_MERID_PARIS, E_CLARK_IGN, DEFAULT_EPS); LambertPoint pt2 = geographicToCartesian(pt1.getX(), pt1.getY(), pt1.getZ(), A_CLARK_IGN, E_CLARK_IGN); - pt2.translate(-168,-60,320); + pt2.translate(-168, -60, 320); //WGS84 refers to greenwich return cartesianToGeographic(pt2, LON_MERID_GREENWICH, A_WGS84, E_WGS84, DEFAULT_EPS); } } - public static LambertPoint convertToWGS84(double x, double y, LambertZone zone){ + /* + * Convert WGS84 -> Lambert + * http://geodesie.ign.fr/contenu/fichiers/documentation/pedagogiques/transfo.pdf + * + */ + + public static LambertPoint convertToLambert(double latitude, double longitude, LambertZone zone) throws NotImplementedException { + + if (zone == Lambert93) { + throw new NotImplementedException(); + } else { + LambertPoint pt1 = geographicToCartesian(longitude - LON_MERID_GREENWICH, latitude, 0, A_WGS84, E_WGS84); + + pt1.translate(168, 60, -320); + + LambertPoint pt2 = cartesianToGeographic(pt1, LON_MERID_PARIS, A_WGS84, E_WGS84, DEFAULT_EPS); + + return geographicToLambert(pt2.getY(), pt2.getX(), zone, LON_MERID_PARIS, E_WGS84); + } + } + + /* + Method not really usefull, just to have two ways of doing the same conversion. + */ + public static LambertPoint convertToLambertByAlg003(double latitude, double longitude, LambertZone zone) throws NotImplementedException { + + if (zone == Lambert93) { + throw new NotImplementedException(); + } else { + LambertPoint pt1 = geographicToCartesian(longitude - LON_MERID_GREENWICH, latitude, 0, A_WGS84, E_WGS84); + + pt1.translate(168, 60, -320); + + LambertPoint pt2 = cartesianToGeographic(pt1, LON_MERID_PARIS, A_WGS84, E_WGS84, DEFAULT_EPS); + + return geographicToLambertAlg003(pt2.getY(), pt2.getX(), zone, LON_MERID_PARIS, E_WGS84); + } + } + + public static LambertPoint convertToWGS84(double x, double y, LambertZone zone) { - LambertPoint pt = new LambertPoint(x,y,0); + LambertPoint pt = new LambertPoint(x, y, 0); return convertToWGS84(pt, zone); } - public static LambertPoint convertToWGS84Deg(double x, double y, LambertZone zone){ + public static LambertPoint convertToWGS84Deg(double x, double y, LambertZone zone) { - LambertPoint pt = new LambertPoint(x,y,0); + LambertPoint pt = new LambertPoint(x, y, 0); return convertToWGS84(pt, zone).toDegree(); } diff --git a/src/test/java/net/yageek/lambert/LambertTest.java b/src/test/java/net/yageek/lambert/LambertTest.java index c7b2e30..9677a08 100644 --- a/src/test/java/net/yageek/lambert/LambertTest.java +++ b/src/test/java/net/yageek/lambert/LambertTest.java @@ -1,9 +1,10 @@ package net.yageek.lambert; +import static junit.framework.Assert.assertNotNull; +import static net.yageek.lambert.LambertZone.*; import static org.junit.Assert.assertEquals; import org.junit.After; -import org.junit.Assert; import org.junit.Before; import org.junit.Test; import org.junit.runner.RunWith; @@ -29,8 +30,9 @@ public void cleanUpStreams() { System.setOut(null); System.setErr(null); } + @Test - public void ResultTest(){ + public void ResultTest() { LambertPoint pt = Lambert.convertToWGS84Deg(994272.661, 113467.422, LambertZone.LambertI); System.out.println("Point latitude:" + pt.getY() + " longitude:" + pt.getX()); @@ -38,11 +40,137 @@ public void ResultTest(){ } @Test - public void Lamber93BugTest(){ - LambertPoint pt = Lambert.convertToWGS84Deg(668832.5384, 6950138.7285,LambertZone.Lambert93); - assertEquals(2.56865,pt.getX(),0.0001); - assertEquals(49.64961,pt.getY(),0.0001); + public void Lambert93BugTest() { + LambertPoint pt = Lambert.convertToWGS84Deg(668832.5384, 6950138.7285, LambertZone.Lambert93); + assertEquals(2.56865, pt.getX(), 0.0001); + assertEquals(49.64961, pt.getY(), 0.0001); + } + + @Test + public void LambertIIExtendedToWgs84Test() { + LambertPoint pt = Lambert.convertToWGS84Deg(618115, 2430676, LambertZone.LambertIIExtended); + assertEquals(2.58331732871, pt.getX(), 0.0001); // Longitude 2.5832231178521186 + assertEquals(48.8741427818, pt.getY(), 0.0001); // Latitude 48.87412734248018 + } + + + @Test + public void LambertAlg0001Test() { + double lat = Lambert.latitudeISOFromLat(0.87266462600, 0.08199188998); + assertNotNull(lat); + assertEquals(1.00552653649, lat, 0.0001); + + double lat2 = Lambert.latitudeISOFromLat(-0.30000000000, 0.08199188998); + assertNotNull(lat2); + assertEquals(-0.30261690063, lat2, 0.0001); + } + + /* + Avec les données de tests pour valider l'algorythme + http://geodesie.ign.fr/contenu/fichiers/documentation/algorithmes/notice/NTG_71.pdf + ALG0003 + */ + @Test + public void LambertAlg0003Test() { + double latitude = 0.87266462600; + double longitude = 0.14551209900; + + LambertPoint lambertPoint = Lambert.geographicToLambertAlg003(latitude, longitude, LambertZone.LambertI, LambertZone.LON_MERID_GREENWICH, LambertZone.E_CLARK_IGN); + + assertEquals(1029705.0818, lambertPoint.getX(), 0.0001); + assertEquals(272723.84730, lambertPoint.getY(), 0.0001); + } + + /* + Calcul par Algorythme ALG003 + */ + @Test + public void ConvertWGS84ToLambertByAlg0003Test() { + double latitude = 48.87412734248018; // 48.8741427818; + double latitude2 = latitude * 400 / 360; // Grad + double radLat = Math.toRadians(latitude); + double longitude = 2.5832231178521186; //2.58331732871; + double longitude2 = longitude * 400 / 360; // Grad + double radLong = Math.toRadians(longitude); + + LambertPoint lambertPoint = Lambert.convertToLambertByAlg003(radLat, radLong, LambertZone.LambertIIExtended); + + assertEquals(618115, lambertPoint.getX(), 1); + assertEquals(2430676, lambertPoint.getY(), 1); } + /* + Methode provenant de + http://geodesie.ign.fr/contenu/fichiers/documentation/pedagogiques/TransformationsCoordonneesGeodesiques.pdf + 3.4 attaquée en direct. + avec les valeurs calculées précedemment sans la translation +towgs84=-168,-60,320,0,0,0,0 + */ + @Test + public void LambertGeographicToLambertTest() { + double latitude = 48.8741427818; + double radLat = Math.toRadians(latitude); + double longitude = 2.58331732871; + double radLong = Math.toRadians(longitude); + + LambertPoint lambertPoint = Lambert.geographicToLambert(radLat, radLong, LambertZone.LambertIIExtended, LambertZone.LON_MERID_GREENWICH, LambertZone.E_CLARK_IGN); + + assertEquals(618062, lambertPoint.getX(), 1); + assertEquals(2430668, lambertPoint.getY(), 1); + } + + /* + Methode provenant de + http://geodesie.ign.fr/contenu/fichiers/documentation/pedagogiques/TransformationsCoordonneesGeodesiques.pdf + 3.4 attaquée en direct + avec les données de test + */ + @Test + public void LambertConvertNTFToLambertTest() { + double latitude = 51.8072313; // Grad + double radLat = Math.toRadians(latitude * 360d / 400d); // Deg before Rad + double longitude = 0.4721669; //Grad + double radLong = Math.toRadians(longitude * 360d / 400d); // Deg before Rad + + LambertPoint lambertPoint = Lambert.geographicToLambert(radLat, radLong, LambertZone.LambertII, LON_MERID_PARIS, E_CLARK_IGN); + + assertEquals(632542.058, lambertPoint.getX(), 0.001); + assertEquals(180804.145, lambertPoint.getY(), 0.01); + } + + /* + Validation de la méthode + http://geodesie.ign.fr/contenu/fichiers/documentation/pedagogiques/TransformationsCoordonneesGeodesiques.pdf + 3.3 attaquée en direct + */ + @Test + public void LamberConvertLambertToNTFTest() { + double X = 1029705.083; + double Y = 272723.849; + + LambertPoint lambertPoint = Lambert.lambertToGeographic(new LambertPoint(X, Y, 0), LambertZone.LambertI, LON_MERID_PARIS, E_CLARK_IGN, DEFAULT_EPS); + + assertEquals(0.145512099, lambertPoint.getX(), 10); // Longitude en rad + assertEquals(0.872664626, lambertPoint.getY(), 10); // Latitude en trad + } + + /* + Test avec translation avant conversion + .translate(168,60,-320); + */ + @Test + public void LambertConvertToLambertTest() { + //double latitude = 48.8741427818; + double latitude = 48.87412734248018; // 48.8741427818; + double radLat = Math.toRadians(latitude); + //double longitude = 2.58331732871; + double longitude = 2.5832231178521186; //2.58331732871; + double radLong = Math.toRadians(longitude); + + LambertPoint lambertPoint = Lambert.convertToLambert(radLat, radLong, LambertZone.LambertIIExtended); + + assertEquals(618115, lambertPoint.getX(), 1); + assertEquals(2430676, lambertPoint.getY(), 1); + } + }