Convert Longitude and Latitude to image X, Y pixel coordinates on Alber's projection

I get many coordinates from a radar, these coordinates formatted pretty much the same as Google maps coordinates, I am not really sure but I asked the company and they told me in order to get the coordinates on Alber's projection, I need to do the following:

You will then be able to create 2 CoordinateReferenceSystem (i.e. coordinate systems), set one as default (which is lon/lat), and set the other with the WKT string (which will be projected x/y). Then you can easily create 2 MathTransform to convert in both directions.

This is the OGC WKT for the Alber's projection:

PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",31.996308],
    PARAMETER["standard_parallel_2",33.996308],
    PARAMETER["latitude_of_center",32.996308],
    PARAMETER["longitude_of_center",35.415901],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,AUTHORITY["EPSG","9001"]]]

So from what I get is, I need to transform from long/lat to that WKT projection to show on Alber's projection map image.

So in GeoTools I used the following code:

    CoordinateReferenceSystem source = CRS.decode("EPSG:4326");
    CoordinateReferenceSystem target = CRS.parseWKT("PROJCS[\"unnamed\", GEOGCS[\"WGS 84\", DATUM[\"WGS_1984\", SPHEROID[\"WGS 84\",6378137,298.2572235629972, AUTHORITY[\"EPSG\",\"7030\"]], AUTHORITY[\"EPSG\",\"6326\"]], PRIMEM[\"Greenwich\",0], UNIT[\"degree\",0.0174532925199433], AUTHORITY[\"EPSG\",\"4326\"]], PROJECTION[\"Albers_Conic_Equal_Area\"], PARAMETER[\"standard_parallel_1\",31.996308], PARAMETER[\"standard_parallel_2\",33.996308], PARAMETER[\"latitude_of_center\",32.996308], PARAMETER[\"longitude_of_center\",35.415901], PARAMETER[\"false_easting\",0], PARAMETER[\"false_northing\",0], UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]");
    MathTransform transform = CRS.findMathTransform(source, target, true);

    Coordinate c = JTS.transform(new Coordinate(34, 35), new Coordinate(), transform);
    System.out.println(c.toString());

And that's the output I get:

(-38422.86847540497, 111410.0483012808, NaN)

Now, it could be because of wrong source coordinate system, but what did he mean by default long/lat system?

Even if I solve this up, how can I make it show these points on my map image? I mean it has to know the width/height of the image doesn't it?

Solution

This albers projection WKT projects from EPSG 4326 to distance in meters from the zero-point (see the parameters). How do we know if it's meters? because the WKT there is a set UNIT label showing it as metre unit.

So how did I use it?

My company gave me two files, the map JPEG file and an XML file for map's constants.

That map constants XML file contains the distance from the zero-point of that map to the corners of the map. So if you have atleast one point on the map, you can find everything.

Things that you need to know to convert it to map X/Y:

How many Kilometers is each pixel (in my case its 1.6 exactly 1 mile)
atleast one known point on the map in distance from zero point

This is how I did it:

        MathTransform transform = CRS.findMathTransform(epsg4326, targetWKT, true);

        DirectPosition2D srcDirectPosition2D
                = new DirectPosition2D(epsg4326, latitude, longitude);
        DirectPosition2D destDirectPosition2D
                = new DirectPosition2D();
        transform.transform(srcDirectPosition2D, destDirectPosition2D);

        double transX = destDirectPosition2D.x;
        double transY = destDirectPosition2D.y;

        int kmPerPixel = mapImage.getWidth / 1024; // It is known to me that my map is 1024x1024km ...

        double x = zeroPointX + ( (transX * 0.001) * kmPerPixel);
        double y = zeroPointY + ( ( (transX * -1) * 0.001) * kmPerPixel);

The zeroPointX and Y you can get from doing the same calculation without adding, on the point you have on the map in distance, like I did with my corner.

Might help some people so I posted what I figured out.