Use high precision to convert image from one EPSG defined CRS to another

Question

This question is based on another question: I'm converting multiple images (.jpg or .png) from one EPSG defined CRS system to another, e.g. from EPSG:31255 to EPSG:25833 (on a server within a Java app).

Currently I'm using this code (GeoTools 22.2, Java 8):

//codeIn file (already includes TOWGS84): http://epsg.io/31255.wkt
String codeOut = "epsg:25833";

public void convert(String pathIn, String pathOut, String codeOut) {
    File fileIn = new File(pathIn);

    if(fileIn.exists() && fileIn.isFile()) {
        try {
            File fileOut = new File(pathOut);
            CoordinateReferenceSystem crsOut = CRS.decode(codeOut);
            AbstractGridFormat format = GridFormatFinder.findFormat(fileIn);
            Hints hints = null;

            AbstractGridCoverage2DReader reader = format.getReader(fileIn, hints);
            GridCoverage2D coverage = reader.read(null);
            reader.dispose();

            GridCoverage2D coverageTransf = (GridCoverage2D) Operations.DEFAULT.resample(coverage, crsOut);
            GeoTiffFormat outFormat = new GeoTiffFormat();
            GridCoverageWriter writer = outFormat.getWriter(fileOut, hints);
            writer.write(coverageTransf, null);
            writer.dispose();
        } catch(Exception e) {
            e.printStackTrace();
        }
    }
}

Now I want to add high precision to it. In Proj4 I'd use this command:

+ellps=bessel +nadgrids=GIS_GRID.gsb +proj=tmerc +lat_0=0 +lon_0=13d20 +k=1 +x_0=0 +y_0=-5000000 +units=m +no_defs +geoidgrids=foo.gtx +to +datum=WGS84 +proj=utm +zone=33 +no_defs

In this similar question there's a link to a website that tells you to copy the two files (.gsb and .gtx) to the user_projections folder, which of course doesn't exist in my app.

My guess is that I'll have to give GridCoverage2D coverage = reader.read(null) some type of GeneralParameterValue[]. Is this the right way to do it and if so, how exactly do I give it the path of the two files and/or the command?

Answer 1

You should be able to solve this by creating a src/main/resources/org/geotools/referencing/factory/gridshift/ folder (it will get rolled up into your jar when you do a mvn install) to drop the GIS_GRID.gsb file into.

You can use some code like this to check the transform is being picked up:

CoordinateReferenceSystem epsg31255 = CRS.decode("epsg:31255");
System.out.println(epsg31255);
System.out.println();
transform = CRS.findMathTransform(CRS.decode("epsg:25833"), epsg31255);
System.out.println(transform);

Currently, this gives me this transform:

CONCAT_MT[INVERSE_MT[PARAM_MT["Transverse_Mercator", 
      PARAMETER["semi_major", 6378137.0], 
      PARAMETER["semi_minor", 6356752.314140356], 
      PARAMETER["central_meridian", 15.0], 
      PARAMETER["latitude_of_origin", 0.0], 
      PARAMETER["scale_factor", 0.9996], 
      PARAMETER["false_easting", 500000.0], 
      PARAMETER["false_northing", 0.0]]], 
  PARAM_MT["Ellipsoid_To_Geocentric", 
    PARAMETER["dim", 2], 
    PARAMETER["semi_major", 6378137.0], 
    PARAMETER["semi_minor", 6356752.314140356]], 
  PARAM_MT["Coordinate Frame Rotation (geog2D domain)", 
    PARAMETER["dx", -601.7055480109889], 
    PARAMETER["dy", -84.2586089323459], 
    PARAMETER["dz", -485.2300592613099], 
    PARAMETER["ex", 4.73539999802488], 
    PARAMETER["ey", 1.3144999994517264], 
    PARAMETER["ez", 5.392999997750596], 
    PARAMETER["ppm", 2.3878715100789094]], 
  PARAM_MT["Geocentric_To_Ellipsoid", 
    PARAMETER["dim", 2], 
    PARAMETER["semi_major", 6377397.155], 
    PARAMETER["semi_minor", 6356078.962818189]], 
  PARAM_MT["Transverse_Mercator", 
    PARAMETER["semi_major", 6377397.155], 
    PARAMETER["semi_minor", 6356078.962818189], 
    PARAMETER["central_meridian", 13.333333333333336], 
    PARAMETER["latitude_of_origin", 0.0], 
    PARAMETER["scale_factor", 1.0], 
    PARAMETER["false_easting", 0.0], 
    PARAMETER["false_northing", -5000000.0]], 
  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]]]

Update

I can confirm that dropping the AT_GIS_GRID.gsb I downloaded from http://www.bev.gv.at/portal/page?_pageid=713,2157075&_dad=portal&_schema=PORTAL into a directory called src/main/resources/org/geotools/referencing/factory/gridshift/ causes the above code to give me this transform:

CONCAT_MT[INVERSE_MT[PARAM_MT["Transverse_Mercator", 
      PARAMETER["semi_major", 6378137.0], 
      PARAMETER["semi_minor", 6356752.314140356], 
      PARAMETER["central_meridian", 15.0], 
      PARAMETER["latitude_of_origin", 0.0], 
      PARAMETER["scale_factor", 0.9996], 
      PARAMETER["false_easting", 500000.0], 
      PARAMETER["false_northing", 0.0]]], 
  INVERSE_MT[PARAM_MT["NTv2", 
      PARAMETER["Latitude and longitude difference file", "AT_GIS_GRID.gsb"]]], 
  PARAM_MT["Transverse_Mercator", 
    PARAMETER["semi_major", 6377397.155], 
    PARAMETER["semi_minor", 6356078.962818189], 
    PARAMETER["central_meridian", 13.333333333333336], 
    PARAMETER["latitude_of_origin", 0.0], 
    PARAMETER["scale_factor", 1.0], 
    PARAMETER["false_easting", 0.0], 
    PARAMETER["false_northing", -5000000.0]], 
  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 found the link to the transform page from the Proj documentation.