Here is the ouput of my test case:
INFO [main] (TestFeatureTransforming.java:28): Original WGS84 coordinates:
(7.0, 43.0, NaN)
INFO [main] (TestFeatureTransforming.java:33): (999870.608360379,
1789915.8576742434, NaN) <-- GeoTools transform
INFO [main] (TestFeatureTransforming.java:37): (980808.5279030021,
1788755.1365062557, NaN) <-- SIGForum.org transform
INFO [main] (TestFeatureTransforming.java:40): (980808.52797554,
1788755.13499408, NaN) <-- Convers sofware transform
And here is the source:
import org.apache.log4j.Logger;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import com.vividsolutions.jts.geom.Coordinate;
public class TestFeatureTransforming {
private static final Logger logger =
Logger.getLogger(TestFeatureTransforming.class);
public static void main(String[] args) throws Exception {
CoordinateReferenceSystem crsWGS84 = CRS.decode("EPSG:4326");
CoordinateReferenceSystem crsLambert2e = CRS.decode("EPSG:27572");
logger.info("WGS84: " + crsWGS84.toWKT());
logger.info("Lambert2: " + crsLambert2e.toWKT());
MathTransform transformFromWSGS4ToLambert2e = CRS.findMathTransform(crsWGS84,
crsLambert2e);
MathTransform transformFromLambert2eToWGS84 =
CRS.findMathTransform(crsLambert2e, crsWGS84);
Coordinate c1 = new Coordinate(7.0, 43.0);
Coordinate c2 = null;
// transforming from WGS84 to LAMBERT2E with GeoTools tranform
c2 = new Coordinate();
JTS.transform(c1, c2, transformFromWSGS4ToLambert2e);
logger.info(c2 + " <-- GeoTools transform");
// transforming from WGS84 to LAMBERT2E with an algo fetched from
www.forumsig.org
double[] coords = WGS842ToLambert2e(c1.x, c1.y);
logger.info("(" + coords[0] + ", " + coords[1] + ", NaN) <-- SIGForum.org
transform");
// tranforming result for (43.0, 7.0) in WGS84 returned by the Convers software
logger.info("(980808.52797554, 1788755.13499408, NaN) <-- Convers sofware
transform");
}
/**
* Tranforming algo fetched from http://www.forumsig.org/showthread.php?t=7418.
*
* @param d_long WGS84 longitude
* @param d_lat WGS84 latitude
* @return Lamber2E cartesian coordinates
*/
public final static double[] WGS842ToLambert2e(double d_long, double d_lat) {
double lambda_w, phi_w;
/**
***********************************************************************************************************
*/
/** 0) degres-minutes-secondes + orientation (d,m,s,o) -> radian * */
/**
***********************************************************************************************************
*/
// Pour la longitude
lambda_w = Math.toRadians(d_long);
// Pour la latitude
phi_w = Math.toRadians(d_lat);
/**
***********************************************************************************************************
*/
/** 1) coordonnées géographiques WGS84 (phi_w,lambda_w) -> coordonnées
cartésiennes WGS84 (X_w,Y_w,Z_w) * */
/**
***********************************************************************************************************
*/
// J'ai utilisé 2 formules données par l'IGN dans un de leur document ainsi que
deux constantes de
// l'ellipsoide de référence du système WGS84 (les deux demi-axes) :
double a_w = 6378137.0;
double b_w = 6356752.314;
// d'où
double e2_w = (a_w * a_w - b_w * b_w) / (a_w * a_w);
// et on a la grande normale de l'ellipsoide WGS84
double N = a_w / Math.sqrt(1 - e2_w * Math.pow(Math.sin(phi_w), 2));
// ainsi on obtient
double X_w = N * Math.cos(phi_w) * Math.cos(lambda_w);
double Y_w = N * Math.cos(phi_w) * Math.sin(lambda_w);
double Z_w = N * (1 - e2_w) * Math.sin(phi_w);
/**
***********************************************************************************************************
*/
/** 2) coordonnées cartésiennes WGS84 (X_w,Y_w,Z_w) -> coordonnées cartésiennes
NTF (X_n,Y_n,Z_n) * */
/**
***********************************************************************************************************
*/
// Là il n'y a qu'un translation à effectuer :
// on a donc
double dX = 168.0;
double dY = 60.0;
double dZ = -320.0;
// et...
double X_n = X_w + dX;
double Y_n = Y_w + dY;
double Z_n = Z_w + dZ;
/**
***********************************************************************************************************
*/
/** 3) coordonnées cartésiennes NTF (X_n,Y_n,Z_n) -> coordonnées géographiques
NTF (phi_n,lambda_n) * */
/**
***********************************************************************************************************
*/
// J'ai utilisé 1 formule donnée par l'IGN toujours dans le même document ainsi
que deux constantes de
// l'ellipsoide
// de référence du système NTF, Clarke 1880 :
double a_n = 6378249.2;
double b_n = 6356515.0;
// d'où
double e2_n = (a_n * a_n - b_n * b_n) / (a_n * a_n);
// on définit une tolérance de convergence
double epsilon = Math.pow(10, -10);
// puis on amorce une boucle de calcul
double p0 = Math.atan(Z_n / Math.sqrt(X_n * X_n + Y_n * Y_n)
* (1 - (a_n * e2_n) / (Math.sqrt(X_n * X_n + Y_n * Y_n + Z_n * Z_n))));
double p1 = Math.atan((Z_n / Math.sqrt(X_n * X_n + Y_n * Y_n))
/ (1 - (a_n * e2_n * Math.cos(p0))
/ (Math.sqrt((X_n * X_n + Y_n * Y_n) * (1 - e2_n * Math.pow(Math.sin(p0),
2))))));
while (!(Math.abs(p1 - p0) < epsilon)) {
p0 = p1;
p1 = Math.atan((Z_n / Math.sqrt(X_n * X_n + Y_n * Y_n))
/ (1 - (a_n * e2_n * Math.cos(p0))
/ (Math.sqrt((X_n * X_n + Y_n * Y_n) * (1 - e2_n * Math.pow(Math.sin(p0),
2))))));
}
double phi_n = p1;
double lambda_n = Math.atan(Y_n / X_n);
/**
*******************************************************************************************************************
*/
/** 4) coordonnées géographiques NTF (phi_n,lambda_n) coordonnées projetées en
Lambert II étendu (X_l2e, Y_l2e) * */
/**
*******************************************************************************************************************
*/
// J'utilise les formules de projection et les constantes fournies par l'IGN
dans un autre document :
// avant tout on définit quelques constantes
double n = 0.7289686274;
double c = 11745793.39;
double Xs = 600000.0;
double Ys = 8199695.768;
double e_n = Math.sqrt(e2_n);
double lambda0 = 0.04079234433198; // correspond à la longitude en radian de
Paris (2°20'14.025" E) par
// rapport à Greenwich
// puis on calcule la latitude isométrique
double L = Math.log(Math.tan(Math.PI / 4 + phi_n / 2)
* Math.pow(((1 - e_n * Math.sin(phi_n)) / (1 + e_n * Math.sin(phi_n))), (e_n /
2)));
// enfin on projette
double X_l2e = Xs + c * Math.exp((-n * L)) * Math.sin(n * (lambda_n - lambda0));
double Y_l2e = Ys - c * Math.exp((-n * L)) * Math.cos(n * (lambda_n - lambda0));
double[] tabXY = new double[2];
tabXY[0] = X_l2e;
tabXY[1] = Y_l2e;
return tabXY;
}
}
BTW I have this interesting error message at the begining of the output:
18 avr. 2007 14:38:27 FactoryRegistry scanForPlugins
ATTENTION: Échec lors de l'initialisation d'un service de catégorie
"MathTransformProvider". La cause est "NoClassDefFoundError:
javax/media/jai/WarpAffine".
Could it explain the imprecision?
Thanks, Vincent.
-----Message d'origine-----
De : Martin Desruisseaux [mailto:[EMAIL PROTECTED]
Envoyé : mercredi 18 avril 2007 13:06
À : zze-M2S FRISON V ext RD-BIZZ-SOP
Cc : [email protected]; Vincent Frison
Objet : Re: [Geotools-gt2-users] WGS84 to Lambert2e transforming
zze-M2S FRISON V ext RD-BIZZ-SOP a écrit :
> Do you mean that sometimes EPSG database could be wrong (even for very usual
> SRS like WGS84)???
No, EPSG database is usually correct. For more details on this "axis order
issue", see:
http://docs.codehaus.org/display/GEOTOOLS/The+axis+order+issue
Could you provide a test case please? With full code source, what are the
expected values, which value you get and how did you computed the expected
values?
Martin
-------------------------------------------------------------------------
This SF.net email is sponsored by DB2 Express
Download DB2 Express C - the FREE version of DB2 express and take
control of your XML. No limits. Just data. Click to get it now.
http://sourceforge.net/powerbar/db2/
_______________________________________________
Geotools-gt2-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/geotools-gt2-users
