[GIS] Converting bbox when importing GeoTIFF with Geotools

coordinate systemextentsgeotiff-tiffgeotools

Here is my code:

String EPSG4326 = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";


File geoTiffFile = new File( "/home/user/files/1geotiff.tif" );

CoordinateReferenceSystem crs = CRS.parseWKT(EPSG4326);
GeoTiffFormat format = new GeoTiffFormat();
Hints hint = new Hints();
hint.put(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, crs );    
hint.put(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);

GeoTiffReader tiffReader = format.getReader( geoTiffFile, hint );
GridCoverage2D coverage = tiffReader.read(null);
Envelope env = coverage.getEnvelope();  

double minX = env.getLowerCorner().getOrdinate(0);
double minY = env.getLowerCorner().getOrdinate(1);

double maxX = env.getUpperCorner().getOrdinate(0);
double maxY = env.getUpperCorner().getOrdinate(1);

System.out.println(coverage.getCoordinateReferenceSystem2D().getName().toString() + " " + coverage.getName() + ": " + minX + "," + minY + "  " + maxX + "," + maxY );

When I run this code, I get the following output:

EPSG:WGS 84 1geotiff: -6345422.834303701,-4399127.000735723  -2688622.8343037013,965672.9992642764

But when creating a Store in my GeoServer with this same GeoTIFF file, I can't preview the Layer until I do a "Compute from data" and "Compute from native bounds", when Geoserver changes the coordinates to:

-58.999479314465475,-38.00123605669721 || -24.998704016910057,9.001542880512202

How can I get the same values using Geotools in my code?

This is the WKT of the TIFF file ( I get this by using geotools toWKT() method):

GEOGCS["WGS 84", 
  DATUM["WGS_1984", 
    SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], 
    AUTHORITY["EPSG","6326"]], 
  PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], 
  UNIT["degree", 0.017453292519943295], 
  AXIS["Longitude", EAST], 
  AXIS["Latitude", NORTH], 
  AUTHORITY["EPSG","4326"]]

Added GDALINFO output:

Driver: GTiff/GeoTIFF
Files: 1geotiff.tif
Size is 6530, 9580
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["latitude_of_origin",-15],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-6345422.834303701296449,965672.999264276353642)
Pixel Size = (560.000000000000000,-560.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2013:04:25 10:34:09
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Produced by CARIS++ 3.0
  TIFFTAG_XRESOLUTION=254
  TIFFTAG_YRESOLUTION=254
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-6345422.834,  965672.999) ( 58d59'58.13"W,  9d 0' 5.55"N)
Lower Left  (-6345422.834,-4399127.001) ( 58d59'58.13"W, 38d 0' 4.45"S)
Upper Right (-2688622.834,  965672.999) ( 24d59'55.33"W,  9d 0' 5.55"N)
Lower Right (-2688622.834,-4399127.001) ( 24d59'55.33"W, 38d 0' 4.45"S)
Center      (-4517022.834,-1716727.001) ( 41d59'56.73"W, 15d51'37.35"S)

Best Answer

I'm still a little unclear on what you are trying to do but if all you need is to simply reproject the native bounds to WGS84 (as GeoServer does) then the following code will work for you:

String wkt = "PROJCS[\"unnamed\"," + 
        "    GEOGCS[\"WGS 84\"," + 
        "        DATUM[\"WGS_1984\"," + 
        "            SPHEROID[\"WGS 84\",6378137,298.257223563," + 
        "                AUTHORITY[\"EPSG\",\"7030\"]]," + 
        "            AUTHORITY[\"EPSG\",\"6326\"]]," + 
        "        PRIMEM[\"Greenwich\",0]," + 
        "        UNIT[\"degree\",0.0174532925199433]," + 
        "        AUTHORITY[\"EPSG\",\"4326\"]]," + 
        "    PROJECTION[\"Mercator_1SP\"]," + 
        "    PARAMETER[\"latitude_of_origin\",-15]," + 
        "    PARAMETER[\"central_meridian\",0]," + 
        "    PARAMETER[\"scale_factor\",1]," + 
        "    PARAMETER[\"false_easting\",0]," + 
        "    PARAMETER[\"false_northing\",0]," + 
        "    UNIT[\"metre\",1," + 
        "        AUTHORITY[\"EPSG\",\"9001\"]]]";

    CoordinateReferenceSystem crs= CRS.parseWKT(wkt);

    ReferencedEnvelope bbox = new ReferencedEnvelope(-6345422.834,-4399127.001,
    -2688622.834,  965672.999, crs);
    System.out.println(bbox);
    boolean lenient = false;
    MathTransform transform = CRS.findMathTransform(crs, DefaultGeographicCRS.WGS84, lenient );
    Envelope res = JTS.transform(bbox, transform);
    System.out.println(res);

Which gives me a slightly different answer (I suspect a copy and paste error somewhere as I've compared it to your solution above and the answers are the same for my inputs):

ReferencedEnvelope[-6345422.834303701 : -4399127.000735723, -2688622.8343037013 : 965672.9992642764]
Env[-58.999479314465475 : -40.90290107674037, -24.38570281128449 : 9.001542880512202]

Update There are two issues with your code (in the question):

Hints hint = new Hints();
hint.put(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, crs );    
hint.put(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);

The first line tells GeoTools to use the WGS84 crs you defined as the default CRS this overrides the one the coverage reads in so it looked like your tif was in WGS84 (and hence reprojecting it to WGS84 did nothing). The second line tells GeoTools to consider projections as longitude, latitude pairs (rather than asking the crs), I think it is unnecessary here. But it just adjusts the order of the resulting bbox axis either without it:

ReferencedEnvelope[-38.00123605669721 : 9.001542880512217, -58.999479314465475 : -24.998704016910057]

or with it on:

ReferencedEnvelope[-58.999479314465475 : -24.998704016910057, -38.00123605669721 : 9.001542880512217]

As you see the answers now match with the GeoTiff providing the CRS.

    GeoTiffFormat format = new GeoTiffFormat();
    Hints hint = new Hints();
   // hint.put(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, DefaultGeographicCRS.WGS84 );
    hint.put(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);

    File geoTiffFile = new File("1geotiff.tif");
    GeoTiffReader tiffReader = format.getReader( geoTiffFile, hint );
    GridCoverage2D coverage = tiffReader.read(null);
    CoordinateReferenceSystem crs = coverage.getCoordinateReferenceSystem();
    System.out.println(crs);
    ReferencedEnvelope bbox = new ReferencedEnvelope(crs);
    bbox.setBounds(coverage.getEnvelope2D());
    boolean lenient = false;
    CoordinateReferenceSystem target = DefaultGeographicCRS.WGS84;

    MathTransform transform = CRS.findMathTransform(crs, target, lenient );
    ReferencedEnvelope res = new ReferencedEnvelope(JTS.transform(bbox, transform),target);
    System.out.println(bbox);
    System.out.println(res);
Related Question