[gdal-dev] Perforce issue in coordinate transformations with Gdal 2.3.2

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Perforce issue in coordinate transformations with Gdal 2.3.2

Stefan Moebius-2
Hi,

As we do lots of coordinate transformations in our (Java based) product, we
have the performance test below in our code base. We recently tried to
upgrade from Gdal 1.9 to 2.3.2, and suddenly our coordinate transformation
tests failed, because they took about twice the time as before.

Should we be calling anything differently when using Gdal 2.x?
Anything else we do wrong?
What could be the problem here?

Regards,
Stefan

The code (extracted from our code base, tested with Java 8 and Junit 4):

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.Random;
import java.util.stream.Collectors;
import java.util.stream.Stream;

import org.gdal.osr.SpatialReference;

public class CoordinateTransformationTest
{
        public CoordinateTransformationTest() {}

        static class SpatialCoordinate
        {
            public final double northing; /// latitude  or northing in
meters
            public final double easting;  /// longitude or easting in meters

            public SpatialCoordinate(double northing, double easting)
            {
                this.northing = northing;
                this.easting = easting;
            }

            @Override
            public boolean equals(Object obj)
            {
                if (this == obj)  return  true;
                if (!(obj instanceof SpatialCoordinate))  return  false;
                SpatialCoordinate  other = (SpatialCoordinate) obj;
                return (easting == other.easting) && (northing ==
other.northing);
            }

            @Override
            public int hashCode()
            {
                return Double.valueOf(northing).hashCode() ^
Integer.rotateLeft(Double.valueOf(easting).hashCode(), 16);
            }
        }

    static class CoordinateTransformation
    {
        private /*final*/ org.gdal.osr.CoordinateTransformation
coordinateTransformation;

        private CoordinateTransformation(SpatialReference source,
SpatialReference target)
        {
            this.coordinateTransformation = new
org.gdal.osr.CoordinateTransformation(source, target);
        }

        public synchronized SpatialCoordinate transform(SpatialCoordinate
from)
        {
            // the JNI interface expects a 3d array!
            double[] coordinates3D = { from.easting, from.northing, 0.0 };
            coordinateTransformation.TransformPoint(coordinates3D);
            return new SpatialCoordinate(coordinates3D[1],
coordinates3D[0]);
        }

        public List<SpatialCoordinate> transform(List<SpatialCoordinate>
coordinates)
        {
            if (coordinates.isEmpty())
                return Collections.emptyList();

            double[][] eastNorth = new double[coordinates.size()][2];

            for (int i = 0; i < coordinates.size(); i++)
            {
                final SpatialCoordinate coord = coordinates.get(i);
                eastNorth[i][0] = coord.easting;
                eastNorth[i][1] = coord.northing;
            }

            transform(eastNorth);

            final List<SpatialCoordinate> result = new
ArrayList<>(coordinates.size());
            for (int i = 0; i < coordinates.size(); i++)
                result.add(new SpatialCoordinate(eastNorth[i][1],
eastNorth[i][0]));

            return result;
        }

        private synchronized void transform(double[][] coordinates)
        {
            assert (coordinates != null && coordinates.length > 0 &&
coordinates[0] != null && coordinates[0].length >= 2);
            coordinateTransformation.TransformPoints(coordinates);
        }
    }

    static SpatialReference wgs84UtmNorth(final int utm)
    {
        final int epsg = 32600 + utm;
        final int meridian = -177 + 6 * (utm - 1);

        final String wkt = String
                .format(
                        "PROJCS[\"WGS 84 / UTM zone %dN\",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[\"degre
e\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4
326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\"
,0],PARAMETER[\"central_meridian\",%d],PARAMETER[\"scale_factor\",0.9996],PA
RAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"met
re\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\
",NORTH],AUTHORITY[\"EPSG\",\"%d\"]]",
                        utm, meridian, epsg);

        return new SpatialReference(wkt);
    }

    static SpatialReference wgs84()
    {
    return new SpatialReference("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[\"degre
e\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4
326\"]]");
    }

    static SpatialReference dhdn3()
    {
    return new SpatialReference("PROJCS[\"DHDN / 3-degree Gauss zone 2
(deprecated)\",GEOGCS[\"DHDN\",DATUM[\"Deutsches_Hauptdreiecksnetz\",SPHEROI
D[\"Bessel
1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],AUTHORITY[\"EPS
G\",\"6314\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"d
egree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\"
,\"4314\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_orig
in\",0],PARAMETER[\"central_meridian\",6],PARAMETER[\"scale_factor\",1],PARA
METER[\"false_easting\",2500000],PARAMETER[\"false_northing\",0],UNIT[\"metr
e\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHO
RITY[\"EPSG\",\"31462\"]]");
    }

    static void testBulkSpeed(
            SpatialReference source,
            SpatialReference target,
            SpatialCoordinate center,
            double epsilon)
    {
        final CoordinateTransformation ct = new
CoordinateTransformation(source, target);
        final Random rand = new Random();
        final int size = 1024;

        for (int i = 0; i < size; ++i)
        {
            final List<SpatialCoordinate> coords = Stream //
                    .generate(() -> new SpatialCoordinate(center.northing +
(rand.nextDouble() * 2 - 1) * epsilon, //
                                                          center.easting +
(rand.nextDouble() * 2 - 1) * epsilon)) //
                    .limit(size) //
                    .collect(Collectors.toList());

            @SuppressWarnings("unused")
            List<SpatialCoordinate> results = ct.transform(coords);
        }
    }

    @org.junit.Test(timeout=1000)
    public void testGeogToProjBulkSpeed()
    {
    SpatialReference utm30 = wgs84UtmNorth(30);
    SpatialReference wgs84 = wgs84();

        testBulkSpeed(wgs84, utm30, new SpatialCoordinate(40.0, -3.0), 360.0
/ 40000000.0);
    }

    @org.junit.Test(timeout=1500)
    public void testProjToGeogBulkSpeed()
    {
    SpatialReference utm30 = wgs84UtmNorth(30);
    SpatialReference wgs84 = wgs84();

        testBulkSpeed(utm30, wgs84, new SpatialCoordinate(0.0, 0.0), 1.0);
    }

    @org.junit.Test(timeout=2000)
    public void testProjToProjBulkSpeed()
    {
    SpatialReference utm30 = wgs84UtmNorth(30);
    SpatialReference utm31 = wgs84UtmNorth(31);

        testBulkSpeed(utm30, utm31, new SpatialCoordinate(0.0, 0.0), 1.0);
    }

    @org.junit.Test(timeout=3000)
    public void testProjToProjViaToWgs84BulkSpeed()
    {
    SpatialReference utm32 = wgs84UtmNorth(32);
    SpatialReference dhdn3 = dhdn3();

        testBulkSpeed(utm32, dhdn3, new SpatialCoordinate(0.0, 0.0), 1.0);
    }
}

This email and the information contained herein is proprietary and confidential and subject to the Amdocs Email Terms of Service, which you may review at https://www.amdocs.com/about/email-terms-of-service


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

smime.p7s (3K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Perforce issue in coordinate transformations with Gdal 2.3.2

Even Rouault-3
Stefan,

I bet you also upgraded your PROJ version in the process ? I can't think of a reason in GDAL itself for a performance change. But on the PROJ side, yes, there might be some explanation
The main one I can see is that in PROJ 4.9.3, the Extended Transverse Mercator is used by default for UTM, which has probably a higher computation time. If you define the GDAL configuration option / environmenet variable OSR_USE_ETMERC to NO, then the older Transverse Mercator method will be used.
PROJ 5.0 has also undergone major changes that might have increased a bit the computation time, but I wouldn't expect them to cause a 2x slowdown.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Perforce issue in coordinate transformations with Gdal 2.3.2

Thomas Knudsen
Prepare/finalize functions were added to pj_inv & pj_fwd  in PROJ PR #731 (https://github.com/OSGeo/proj.4/pull/731), in order to resolve ticket #700 (https://github.com/OSGeo/proj.4/issues/700).
These add some overhead as well but IIRC, more like +15% (much depending on the actual projection) than +100%.

Den man. 28. jan. 2019 kl. 17.38 skrev Even Rouault <[hidden email]>:
Stefan,

I bet you also upgraded your PROJ version in the process ? I can't think of a reason in GDAL itself for a performance change. But on the PROJ side, yes, there might be some explanation
The main one I can see is that in PROJ 4.9.3, the Extended Transverse Mercator is used by default for UTM, which has probably a higher computation time. If you define the GDAL configuration option / environmenet variable OSR_USE_ETMERC to NO, then the older Transverse Mercator method will be used.
PROJ 5.0 has also undergone major changes that might have increased a bit the computation time, but I wouldn't expect them to cause a 2x slowdown.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev