1

My goal is to transform a bounding box from EPSG:4326 to EPSG:32632. For the bounding box a ReferencedEnvelope is used. After the transformation I noticed that the resulting ReferencedEnvelope does not have the expected coordinates.

In the following code I transformed the min and max point of the bounding box separately and I get the expected result. Why is there a difference between transforming the individual points and transforming the ReferencedEnvelope?

package test;

import org.geotools.geometry.*;
import org.geotools.geometry.jts.*;
import org.geotools.referencing.*;
import org.opengis.geometry.*;
import org.opengis.referencing.crs.*;
import org.opengis.referencing.operation.*;

public class GtTransform
{
  public static void main(String[] args)
  {
    try
    {
      CoordinateReferenceSystem crs4326 = CRS.decode("EPSG:4326");
      CoordinateReferenceSystem crs32632 = CRS.decode("EPSG:32632");
      MathTransform transform = CRS.findMathTransform(crs4326, crs32632);
      
      double lon1 = 15.46875;
      double lat1 = 40.979898069620134;
      DirectPosition position4326_1 = new DirectPosition2D(crs4326, lat1, lon1);
      DirectPosition position32632_1 = transform.transform(position4326_1, null);
      System.out.println("Point (EPSG:4326): " + position4326_1);
      System.out.println("Transformed Point (EPSG:32632): " + position32632_1);

      double lon2 = 16.171875;
      double lat2 = 41.50857729743934;
      DirectPosition position4326_2 = new DirectPosition2D(crs4326, lat2, lon2);
      DirectPosition position32632_2 = transform.transform(position4326_2, null);
      System.out.println("Point (EPSG:4326): " + position4326_2);
      System.out.println("Transformed Point (EPSG:32632): " + position32632_2);

      ReferencedEnvelope env4326 = new ReferencedEnvelope(lat1, lat2, lon1, lon2, crs4326);
      System.out.println("Envelope 4326: " + env4326);
      ReferencedEnvelope env32632 = env4326.transform(crs32632, false, 20);
      System.out.println("Envelope 32632: " + env32632);
    }
    catch (Exception e)
    {
      e.printStackTrace();
    }
  }
}

Output:

Point (EPSG:4326): DirectPosition2D[40.979898069620134, 15.46875]
Transformed Point (EPSG:32632): GeneralDirectPosition[1044359.3127961868, 4556723.937322532]
Point (EPSG:4326): DirectPosition2D[41.50857729743934, 16.171875]
Transformed Point (EPSG:32632): GeneralDirectPosition[1098695.5582549274, 4620118.653484692]
Envelope 4326: ReferencedEnvelope[40.979898069620134 : 41.50857729743934, 15.46875 : 16.171875] DefaultGeographicCRS[EPSG:WGS 84] AXIS["Geodetic latitude", NORTH] AXIS["Geodetic longitude", EAST]
Envelope 32632: ReferencedEnvelope[1039968.4188604902 : 1103569.4286610063, 4556723.937322532 : 4620118.653484692] DefaultProjectedCRS[EPSG:WGS 84 / UTM zone 32N] AXIS["Easting", EAST] AXIS["Northing", NORTH]

Output for UTM zone 33N

Same code as above, just replaced 32632 with 32633 in the whole file. There is still a difference.

Point (EPSG:4326): DirectPosition2D[40.979898069620134, 15.46875]
Transformed Point (EPSG:32633): GeneralDirectPosition[539434.6300387884, 4536631.349536699]
Point (EPSG:4326): DirectPosition2D[41.50857729743934, 16.171875]
Transformed Point (EPSG:32633): GeneralDirectPosition[597795.8703118529, 4595879.323949931]
Envelope 4326: ReferencedEnvelope[40.979898069620134 : 41.50857729743934, 15.46875 : 16.171875] DefaultGeographicCRS[EPSG:WGS 84] AXIS["Geodetic latitude", NORTH] AXIS["Geodetic longitude", EAST]
Envelope 32633: ReferencedEnvelope[539118.0648996164 : 598587.3949375926, 4536631.349536699 : 4595879.323949931] DefaultProjectedCRS[EPSG:WGS 84 / UTM zone 33N] AXIS["Easting", EAST] AXIS["Northing", NORTH]
2
  • The first problem is probably that your points are in UTM zone 33 not 32, that brings the offset down from 4km to 300m. Not sure why that remains I'll try to look later
    – Ian Turton
    Jun 2, 2023 at 10:21
  • @IanTurton Sorry about that, I used the CRS provided by a WMTS layer showing whole Italy. For my example code I just picked a point more or less randomly. Unfortunately this does not resolve the problem. Please see my edit.
    – gillesB
    Jun 2, 2023 at 10:37

1 Answer 1

3

Actually both results are correct :-)

enter image description here

What is wrong is your assumption that the corners in WGS84 are still the corners in EPSG:326632. As you can see in the map I've drawn the bounding box becomes slightly twisted when you project it and so you get back slightly different corners which give you the different answers than if you transform the points and then make an envelope.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.