All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.neo4j.values.storable.CRSCalculator Maven / Gradle / Ivy

There is a newer version: 5.25.1
Show newest version
/*
 * Copyright (c) "Neo4j"
 * Neo4j Sweden AB [http://neo4j.com]
 *
 * This file is part of Neo4j.
 *
 * Neo4j is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see .
 */
package org.neo4j.values.storable;

import java.util.Arrays;
import java.util.Collections;
import java.util.List;

import org.neo4j.internal.helpers.collection.Pair;

import static java.lang.Math.asin;
import static java.lang.Math.atan2;
import static java.lang.Math.cos;
import static java.lang.Math.pow;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
import static java.lang.Math.toDegrees;
import static java.lang.Math.toRadians;

public abstract class CRSCalculator
{
    public abstract double distance( PointValue p1, PointValue p2 );

    public abstract List> boundingBox( PointValue center, double distance );

    public abstract boolean withinBBox( PointValue point, PointValue lowerLeft, PointValue upperRight );

    public abstract List> computeBBoxes( PointValue lowerLeft, PointValue upperRight );

    protected static double pythagoras( double[] a, double[] b )
    {
        double sqrSum = 0.0;
        for ( int i = 0; i < a.length; i++ )
        {
            double diff = a[i] - b[i];
            sqrSum += diff * diff;
        }
        return sqrt( sqrSum );
    }

    public static class CartesianCalculator extends CRSCalculator
    {
        int dimension;

        CartesianCalculator( int dimension )
        {
            this.dimension = dimension;
        }

        @Override
        public double distance( PointValue p1, PointValue p2 )
        {
            assert p1.getCoordinateReferenceSystem().getDimension() == dimension;
            assert p2.getCoordinateReferenceSystem().getDimension() == dimension;
            return pythagoras( p1.coordinate(), p2.coordinate() );
        }

        @Override
        public List> boundingBox( PointValue center, double distance )
        {
            assert center.getCoordinateReferenceSystem().getDimension() == dimension;
            double[] coordinates = center.coordinate();
            double[] min = new double[dimension];
            double[] max = new double[dimension];
            for ( int i = 0; i < dimension; i++ )
            {
                min[i] = coordinates[i] - distance;
                max[i] = coordinates[i] + distance;
            }
            CoordinateReferenceSystem crs = center.getCoordinateReferenceSystem();
            return Collections.singletonList( Pair.of( Values.pointValue( crs, min ), Values.pointValue( crs, max ) ) );
        }

        @Override
        public boolean withinBBox( PointValue point, PointValue lowerLeft, PointValue upperRight )
        {
            assert point.getCoordinateReferenceSystem().equals( lowerLeft.getCoordinateReferenceSystem() ) &&
                   point.getCoordinateReferenceSystem().equals( upperRight.getCoordinateReferenceSystem() );

            var check = point.coordinate();
            var ll = lowerLeft.coordinate();
            var ur = upperRight.coordinate();
            for ( int i = 0; i < check.length; i++ )
            {
                if ( check[i] < ll[i] || check[i] > ur[i] )
                {
                    return false;
                }
            }
            return true;
        }

        @Override
        public List> computeBBoxes( PointValue lowerLeft, PointValue upperRight )
        {
            return List.of( Pair.of( lowerLeft, upperRight ) );
        }
    }

    public static class GeographicCalculator extends CRSCalculator
    {
        public static final double EARTH_RADIUS_METERS = 6378140.0;
        private static final double EXTENSION_FACTOR = 1.0001;
        int dimension;

        GeographicCalculator( int dimension )
        {
            this.dimension = dimension;
        }

        @Override
        public double distance( PointValue p1, PointValue p2 )
        {
            assert p1.getCoordinateReferenceSystem().getDimension() == dimension;
            assert p2.getCoordinateReferenceSystem().getDimension() == dimension;
            double[] c1Coord = p1.coordinate();
            double[] c2Coord = p2.coordinate();
            double[] c1 = new double[]{toRadians( c1Coord[0] ), toRadians( c1Coord[1] )};
            double[] c2 = new double[]{toRadians( c2Coord[0] ), toRadians( c2Coord[1] )};
            double dx = c2[0] - c1[0];
            double dy = c2[1] - c1[1];
            double alpha = pow( sin( dy / 2 ), 2.0 ) + cos( c1[1] ) * cos( c2[1] ) * pow( sin( dx / 2.0 ), 2.0 );
            double greatCircleDistance = 2.0 * atan2( sqrt( alpha ), sqrt( 1 - alpha ) );
            if ( dimension == 2 )
            {
                return EARTH_RADIUS_METERS * greatCircleDistance;
            }
            else if ( dimension == 3 )
            {
                // get average height
                double avgHeight = (p1.coordinate()[2] + p2.coordinate()[2]) / 2;
                double distance2D = (EARTH_RADIUS_METERS + avgHeight) * greatCircleDistance;

                double[] a = new double[dimension - 1];
                double[] b = new double[dimension - 1];
                a[0] = distance2D;
                b[0] = 0.0;
                for ( int i = 1; i < dimension - 1; i++ )
                {
                    a[i] = 0.0;
                    b[i] = c1Coord[i + 1] - c2Coord[i + 1];
                }
                return pythagoras( a, b );
            }
            else
            {
                // The above calculation works for more than 3D if all higher dimensions are orthogonal to the 3rd dimension.
                // This might not be true in the general case, and so until we genuinely support higher dimensions fullstack
                // we will explicitly disabled them here for now.
                throw new UnsupportedOperationException( "More than 3 dimensions are not supported for distance calculations." );
            }
        }

        @Override
        // http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates
        // But calculating in degrees instead of radians to avoid rounding errors
        public List> boundingBox( PointValue center, double distance )
        {
            if ( distance == 0.0 )
            {
                return Collections.singletonList( Pair.of( center, center ) );
            }

            // Extend the distance slightly to assure that all relevant points lies inside the bounding box,
            // with rounding errors taken into account
            double extendedDistance = distance * EXTENSION_FACTOR;

            double lat = center.coordinate()[1];
            double lon = center.coordinate()[0];

            double r = extendedDistance / EARTH_RADIUS_METERS;

            double latMin = lat - toDegrees( r );
            double latMax = lat + toDegrees( r );

            // If your query circle includes one of the poles
            if ( latMax >= 90 )
            {
                return Collections.singletonList( boundingBoxOf( -180, 180, latMin, 90, center, distance ) );
            }
            else if ( latMin <= -90 )
            {
                return Collections.singletonList( boundingBoxOf( -180, 180, -90, latMax, center, distance ) );
            }
            else
            {
                double deltaLon = toDegrees( asin( sin( r ) / cos( toRadians( lat ) ) ) );
                double lonMin = lon - deltaLon;
                double lonMax = lon + deltaLon;

                // If you query circle wraps around the dateline
                if ( lonMin < -180 && lonMax > 180 )
                {
                    // Large rectangle covering all longitudes
                    return Collections.singletonList( boundingBoxOf( -180, 180, latMin, latMax, center, distance ) );
                }
                else if ( lonMin < -180 )
                {
                    // two small rectangles east and west of dateline
                    Pair box1 = boundingBoxOf( lonMin + 360, 180, latMin, latMax, center, distance );
                    Pair box2 = boundingBoxOf( -180, lonMax, latMin, latMax, center, distance );
                    return Arrays.asList( box1, box2 );
                }
                else if ( lonMax > 180 )
                {
                    // two small rectangles east and west of dateline
                    Pair box1 = boundingBoxOf( lonMin, 180, latMin, latMax, center, distance );
                    Pair box2 = boundingBoxOf( -180, lonMax - 360, latMin, latMax, center, distance );
                    return Arrays.asList( box1, box2 );
                }
                else
                {
                    return Collections.singletonList( boundingBoxOf( lonMin, lonMax, latMin, latMax, center, distance ) );
                }
            }
        }

        @Override
        public List> computeBBoxes( PointValue lowerLeft, PointValue upperRight )
        {
            assert lowerLeft.getCoordinateReferenceSystem().equals( upperRight.getCoordinateReferenceSystem() );

            double[] ll = lowerLeft.coordinate();
            double[] ur = upperRight.coordinate();
            if ( ll[0] <= ur[0] )
            {   //ll is west of ur or we cross the dateline with ll being east of ur,
                //either way we get a simple continuous range from lower to upper
                return List.of( Pair.of( lowerLeft, upperRight ) );
            }
            else
            {
                //ll > ur either means that we crossed the dateline, and we must
                //scan from lower to the dateline, and then from the dateline to upper
                if ( ll[0] >= 0 && ur[0] < 0 )
                {
                    return List.of( Pair.of( lowerLeft, withLongitude( 180, upperRight ) ),
                                    Pair.of( withLongitude( -180, lowerLeft ), upperRight ) );
                }
                else
                {
                    // or that ll is "east" of ur, in which case we must scan the regions not included between the
                    //two points.
                    return List.of( Pair.of( withLongitude( -180, lowerLeft ), upperRight ),
                                    Pair.of( lowerLeft, withLongitude( 180, upperRight ) ) );
                }
            }
        }

        @Override
        public boolean withinBBox( PointValue point, PointValue lowerLeft, PointValue upperRight )
        {
            assert point.getCoordinateReferenceSystem().equals( lowerLeft.getCoordinateReferenceSystem() ) &&
                   point.getCoordinateReferenceSystem().equals( upperRight.getCoordinateReferenceSystem() );
            double[] check = point.coordinate();

            double[] ll = lowerLeft.coordinate();
            double[] ur = upperRight.coordinate();
            //If lowerLeft is north of upperRight the bbox will always be empty
            if ( isNorthOf( ll, ur ) )
            {
                return false;
            }

            int i = 0;
            //first check latitude and longitude
            for ( ; i < 2; i++ )
            {
                if ( ll[i] <= ur[i] )
                {
                    if ( check[i] < ll[i] || check[i] > ur[i] )
                    {
                        return false;
                    }
                }
                else
                {
                    if ( check[i] > ur[i] && check[i] < ll[i] )
                    {
                        return false;
                    }
                }
            }

            //check potential extra dimensions
            for ( ; i < check.length; i++ )
            {
                if ( check[i] < ll[i] || check[i] > ur[i] )
                {
                    return false;
                }
            }
            return true;
        }

        private static Pair boundingBoxOf( double minLon, double maxLon, double minLat, double maxLat, PointValue center,
                double distance )
        {
            CoordinateReferenceSystem crs = center.getCoordinateReferenceSystem();
            int dimension = center.getCoordinateReferenceSystem().getDimension();
            double[] min = new double[dimension];
            double[] max = new double[dimension];
            min[0] = minLon;
            min[1] = minLat;
            max[0] = maxLon;
            max[1] = maxLat;
            if ( dimension > 2 )
            {
                double[] coordinates = center.coordinate();
                for ( int i = 2; i < dimension; i++ )
                {
                    min[i] = coordinates[i] - distance;
                    max[i] = coordinates[i] + distance;
                }
            }
            return Pair.of( Values.pointValue( crs, min ), Values.pointValue( crs, max ) );
        }

        private PointValue withLongitude( double longitude, PointValue point )
        {
            double[] coordinates = Arrays.copyOf( point.coordinate(), point.coordinate().length );
            coordinates[0] = longitude;
            return Values.pointValue( point.getCoordinateReferenceSystem(), coordinates );
        }

        private boolean isNorthOf( double[] p1, double[] p2 )
        {
            return p1[1] > p2[1];
        }
    }
}




© 2015 - 2024 Weber Informatics LLC | Privacy Policy