0

I am transforming a coordinate between Lambert Conformal and WGS84. The values returned by GeoTools do not seem correct.

From:

PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unnamed",
            SPHEROID["Sphere",6371229,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["latitude_of_origin",50],
    PARAMETER["central_meridian",-105],
    PARAMETER["standard_parallel_1",50],
    PARAMETER["standard_parallel_2",20],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Metre",1],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

To:

EPSG:4326

When I transform the point -2133754.9447882758 1556902.2763524954 I get differing results between proj4 and GeoTools.

Proj4: -141.538155 60.069556
Geotools: -141.538155 60.235687

As I move toward the equator, the difference in latitude of the transformed point approaches zero. Looks like a datum shift error. I have reason to believe that the proj4 transform is correct as the raster data I am working with is defined with the lambert conformal CRS shown above. When I warp to epsg:4326 with GDAL, the result is correct. When I warp using geotools, the resulting image is shifted north progressively from the equator northward. Any thoughts as to why geotools is incorrect?

echo -2133754.9447882758 1556902.2763524954 | cs2cs -f "%f" +proj=lcc +lat_0=50 +lon_0=-105 +lat_1=50 +lat_2=20 +x_0=0 +y_0=0 +R=6371229 +units=m +no_defs +to +init=epsg:4326
-141.538155 60.069556 0.000000

Geotools code to transform a point v26.3:

        String wktLC = "PROJCS[\"unnamed\",\n" +
                "    GEOGCS[\"Coordinate System imported from GRIB file\",\n" +
                "        DATUM[\"unnamed\",\n" +
                "            SPHEROID[\"Sphere\",6371229,0]],\n" +
                "        PRIMEM[\"Greenwich\",0],\n" +
                "        UNIT[\"degree\",0.0174532925199433,\n" +
                "            AUTHORITY[\"EPSG\",\"9122\"]]],\n" +
                "    PROJECTION[\"Lambert_Conformal_Conic_2SP\"],\n" +
                "    PARAMETER[\"latitude_of_origin\",50],\n" +
                "    PARAMETER[\"central_meridian\",-105],\n" +
                "    PARAMETER[\"standard_parallel_1\",50],\n" +
                "    PARAMETER[\"standard_parallel_2\",20],\n" +
                "    PARAMETER[\"false_easting\",0],\n" +
                "    PARAMETER[\"false_northing\",0],\n" +
                "    UNIT[\"Metre\",1],\n" +
                "    AXIS[\"Easting\",EAST],\n" +
                "    AXIS[\"Northing\",NORTH]]";

        crsLCGraf = CRS.parseWKT(wktLC);
        crsWGS84 = CRS.decode("EPSG:4326");

        CoordinateReferenceSystem from = CRS.parseWKT(wktLC);
        CoordinateReferenceSystem to = CRS.decode("EPSG:4326");

        WKTReader reader = new WKTReader(new GeometryFactory());
        Geometry point = reader.read("POINT (-2133754.9447882758 1556902.2763524954)");

        MathTransform transform = CRS.findMathTransform(from, to, true);

        Point p2 = (Point) JTS.transform(point, transform);

1 Answer 1

0

The fact you have to set the lenient flag in the transform is probably the problem, this is giving you the transform:

CONCAT_MT[INVERSE_MT[PARAM_MT["Lambert_Conformal_Conic_2SP", 
      PARAMETER["semi_major", 6371229.0], 
      PARAMETER["semi_minor", 6371229.0], 
      PARAMETER["central_meridian", -105.0], 
      PARAMETER["latitude_of_origin", 50.0], 
      PARAMETER["standard_parallel_1", 50.0], 
      PARAMETER["false_easting", 0.0], 
      PARAMETER["false_northing", 0.0], 
      PARAMETER["scale_factor", 1.0], 
      PARAMETER["standard_parallel_2", 20.0]]], 
  PARAM_MT["Molodenski", 
    PARAMETER["dim", 2], 
    PARAMETER["dx", 0.0], 
    PARAMETER["dy", 0.0], 
    PARAMETER["dz", 0.0], 
    PARAMETER["src_semi_major", 6371229.0], 
    PARAMETER["src_semi_minor", 6371229.0], 
    PARAMETER["tgt_semi_major", 6378137.0], 
    PARAMETER["tgt_semi_minor", 6356752.314245179]], 
  PARAM_MT["Affine", 
    PARAMETER["num_row", 3], 
    PARAMETER["num_col", 3], 
    PARAMETER["elt_0_0", 0.0], 
    PARAMETER["elt_0_1", 1.0], 
    PARAMETER["elt_1_0", 1.0], 
    PARAMETER["elt_1_1", 0.0]]]

I can't work out how to convince cs2cs to print out it's transform but I suspect that it is picking up a grid transform from somewhere which gives a different (more accurate?) answer.

2
  • Forgive me, I don't understand the details of how geotools performs the transforms from one CRS to another. If I set the lenient flag to false, I get an error that Bursa wolf parameters are required. Adding TOWGS84[0,0,0,0,0,0,0] removes the error, but the transform result is the same (as expected I suppose). Shouldn't the datums described in the WKT above and the datum described in EPSG:4326 be enough for the transform...including the datum transform to be as precise as possible. Apr 8 at 21:28
  • I suspect that cs2cs (and proj in general) is picking some other value for TOWGS84 but I don't know.
    – Ian Turton
    Apr 9 at 9:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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