locationtech / proj4j

Java port of the Proj.4 library for coordinate reprojection
Other
184 stars 73 forks source link

*2 in stereographic projection near the equator #58

Closed ruzsaz closed 3 years ago

ruzsaz commented 4 years ago

When I set the base point for stereographic projection on the equator, an unexpected *2 appears in the results compared to the case when the base point is only near the equator.

Example code:

CRSFactory csFactory = new CRSFactory();
CoordinateReferenceSystem STERE0 = csFactory.createFromParameters("STERE", "+proj=stere +lat_0=0.0 +lon_0=0.0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs");
CoordinateReferenceSystem STERE1 = csFactory.createFromParameters("STERE", "+proj=stere +lat_0=0.000001 +lon_0=0.0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs");
CoordinateReferenceSystem WGS84 = csFactory.createFromParameters("WGS84", "+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs");
CoordinateTransform transformer0 = new CoordinateTransformFactory().createTransform(WGS84, STERE0);
CoordinateTransform transformer1 = new CoordinateTransformFactory().createTransform(WGS84, STERE1);

ProjCoordinate pc = new ProjCoordinate(1, 1);
ProjCoordinate result0 = new ProjCoordinate();
ProjCoordinate result1 = new ProjCoordinate();

transformer0.transform(pc, result0);
transformer1.transform(pc, result1);
System.out.println(result0.x + ", " + result0.y);
System.out.println(result1.x + ", " + result1.y);

And the result is:

222627.90213684924, 221171.23231657242 111313.95105169504, 110585.50558411982

I'm not good at projections, but it seems weird for me. In the file https://github.com/locationtech/proj4j/blob/master/src/main/java/org/locationtech/proj4j/proj/StereographicAzimuthalProjection.java I can see from line 146:

            case OBLIQUE:
                A = akm1 / (cosphi0 * (1. + sinphi0 * sinX + cosphi0 * cosX * coslam));
                xy.y = A * (cosphi0 * sinX - sinphi0 * cosX * coslam);
                xy.x = A * cosX;
                break;
            case EQUATOR:
                A = 2. * akm1 / (1. + cosX * coslam);
                xy.y = A * sinX;
                xy.x = A * cosX;
                break;

What is that 2* there, and why is it necessary?

(I use version 1.1.0 from maven.)

pomadchin commented 3 years ago

I think it is just a typo, https://github.com/OSGeo/PROJ/blob/8.0.0/src/projections/stere.cpp#L77