Transform cartesian pixel-data-array to lat/lon pixel-data-array

Question

I have an image (basically, I get raw image data as 1024x1024 pixels) and the position in lat/lon of the center pixel of the image.

Each pixel represents the same fixed pixel scale in meters (e.g. 30m per pixel).

Now, I would like to draw the image onto a map which uses the coordinate reference system "EPSG:4326" (WGS84).

When I draw it by defining just corners in lat/lon of the image, depending on a "image size in pixel * pixel scale" calculation and converting the distances from the center point to lat/lon coordinates of each corner, I suppose, the image is not correctly drawn onto the map. By the term "not correctly drawn" I mean, that the image seems to be shifted and also the contents of the image are not at the map location, where I expected them to be.

I suppose this is the case because I "mix" a pixel scaled image and a "EPSG:4326" coordinate reference system.

Now, with the information I have given, can I transform the whole pixel matrix from fixed pixel scale base to a new pixel matrix in the "EPSG:4326" coordinate reference system, using Geotools? Of course, the transformation must be dependant on the center position in lat/lon, that I have been given, and on the pixel scale.

I wonder if using something like this would point me into the correct direction:

MathTransform transform = CRS.findMathTransform(DefaultGeocentricCRS.CARTESIAN, DefaultGeographicCRS.WGS84, true);

  DirectPosition2D srcDirectPosition2D = new DirectPosition2D(DefaultGeocentricCRS.CARTESIAN, degreeLat.getDegree(), degreeLon.getDegree());
  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);

(got this code from another SO thread and already modified it a little bit, but still wonder if this is the correct starting point for my problem.)

I only suppose that my original image coordinate reference system is of the type DefaultGeocentricCRS.CARTESIAN. Can someone confirm this?

And from here on, is this the correct start to use Geotools for this kind of problem solving, or am I on the complete wrong path?

Additionally, I would like to add that this would be used in a quiet dynamic system. So my image update would be about 10Hz and the transormations have to be performed accordingly often. Again, is this initial thought of mine leading to a solution, or do you have other solutions for solving my problem?

Thank you very much, Kiamur

Answer 1

This is not as simple as it might sound. You are essentially trying to define an area on a sphere (ellipsoid technically) using a flat square. As such there is no "correct" way to do it, so you will always end up with some distortion. Without knowing exactly where your image came from there is no way to answer this exactly but the following code provides you with 3 different possible answers:

The first two make use of GeoTools' GeodeticCalculator to calculate the corner points using bearings and distances. These are the blue "square" and the green "square" above. The blue is calculating the corners directly while the green calculates the edges and infers the corners from the intersections (that's why it is squarer).

final int width = 1024, height = 1024;
GeometryFactory gf = new GeometryFactory();
Point centre = gf.createPoint(new Coordinate(0,51));
WKTWriter writer = new WKTWriter();
//direct method
GeodeticCalculator calc = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
calc.setStartingGeographicPoint(centre.getX(), centre.getY());
double height2 = height/2.0;
double width2 = width/2.0;
double dist = Math.sqrt(height2*height2+width2 *width2);
double bearing = 45.0;
Coordinate[] corners = new Coordinate[5];
for (int i=0;i<4;i++) {
  calc.setDirection(bearing, dist*1000.0 );
  Point2D corner = calc.getDestinationGeographicPoint();
  corners[i] = new Coordinate(corner.getX(),corner.getY());
  bearing+=90.0;
}
corners[4] = corners[0];
Polygon bbox = gf.createPolygon(corners);
System.out.println(writer.write(bbox));

double[] edges = new double[4];
bearing = 0;
for(int i=0;i<4;i++) {
  calc.setDirection(bearing, height2*1000.0 );
  Point2D corner = calc.getDestinationGeographicPoint();

  if(i%2 ==0) {
    edges[i] = corner.getY();
  }else {
    edges[i] = corner.getX();
  }
  bearing+=90.0;
}

corners[0] = new Coordinate( edges[1],edges[0]);
corners[1] = new Coordinate( edges[1],edges[2]);
corners[2] = new Coordinate( edges[3],edges[2]);
corners[3] = new Coordinate( edges[3],edges[0]);

corners[4] = corners[0];
bbox = gf.createPolygon(corners);
System.out.println(writer.write(bbox));

Another way to do this is to transform the centre point into a projection that is "flatter" and use simple addition to calculate the corners and then reverse the transformation. To do this we can use the AUTO projection defined by the OGC WMS Specification to generate an Orthographic projection centred on our point, this gives the red "square" which is very similar to the blue one.

String code = "AUTO:42003," + centre.getX() + "," + centre.getY();
// System.out.println(code);
CoordinateReferenceSystem auto = CRS.decode(code);
// System.out.println(auto);
MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84,
    auto);
MathTransform rtransform = CRS.findMathTransform(auto,DefaultGeographicCRS.WGS84);
Point g = (Point)JTS.transform(centre, transform);

width2 *=1000.0;
height2 *= 1000.0;
corners[0] = new Coordinate(g.getX()-width2,g.getY()-height2);
corners[1] = new Coordinate(g.getX()+width2,g.getY()-height2);
corners[2] = new Coordinate(g.getX()+width2,g.getY()+height2);
corners[3] = new Coordinate(g.getX()-width2,g.getY()+height2);
corners[4] = corners[0];
bbox = gf.createPolygon(corners);
bbox = (Polygon)JTS.transform(bbox, rtransform);
System.out.println(writer.write(bbox));

Which solution to use is a matter of taste, and depends on where your image came from but I suspect that either the red or the blue will be best. If you need to do this at 10Hz then you will need to test them for speed, but I suspect that transforming the images will be the bottle neck.

Once you have your bounding box setup to your satisfaction you can convert you (unreferenced) image to a georeferenced coverage using:

GridCoverageFactory factory = CoverageFactoryFinder.getGridCoverageFactory(null);
GridCoverage2D gc = factory.create("name", image, new ReferencedEnvelope(bbox.getEnvelopeInternal(),DefaultGeographicCRS.WGS84));
String fileName = "myImage.tif";
AbstractGridFormat format = GridFormatFinder.findFormat(fileName);
File out = new File(fileName);
GridCoverageWriter writer = format.getWriter(out);
try {
  writer.write(gc, null);
  writer.dispose();
} catch (IllegalArgumentException | IOException e) {
  // TODO Auto-generated catch block
  e.printStackTrace();
}