3

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?

1 Answer 1

3

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.

12
  • Is it possible to read the gsb file at runtime, so I can change it without exporting the jar again? What about the .gtx file and how do I use the remaining parameters in the command (e.g. the -5 mio. one) on that image?
    – Neph
    Jan 30, 2020 at 15:04
  • It may be, GeoServer manages it - there is probably Spring magic involved. You can add a new projection string (again see the GeoServer docs for details)
    – Ian Turton
    Jan 30, 2020 at 15:09
  • So is this more of a GeoServer problem than a GeoTools one and the latter doesn't really support adding these files at runtime, right? Is there anything I have to add to my code - there isn't too much more to it than what I posted in my question - to make it use the .gsb and .gtx file? Something like a certain call to a function to tell it to use it or does everything happen automatically (as with the .prj and wld. files) as long as the files are in the folder you mentioned?
    – Neph
    Jan 30, 2020 at 15:32
  • GeoTools is the underlying lib of GeoServer (which has more users and docs) - all you should need to do is drop the files into the folder I said - then the referencing engine will pick them up
    – Ian Turton
    Jan 30, 2020 at 16:16
  • Do I have to use a specific name? I copied the files into that folder, renamed them to "foo.gsb" and "foo.gtx" and ran your test code but the resulting transform didn't mention any of the files, so they probably weren't used. Even if they were used, how does GeoServer know what part of EPSG:25833 they should be used on without giving it the additional information of the command I posted in the question?
    – Neph
    Jan 31, 2020 at 12:11

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.