42

Given a coordinate (lat, long), I am trying to calculate a square bounding box that is a given distance (e.g. 50km) away from the coordinate. So as input I have lat, long and distance and as output I would like two coordinates; one being the south-west (bottom-left) corner and one being the north-east (top-right) corner. I have seen a couple of answers on here that try to address this question in Python, but I am looking for a Java implementation in particular.

Just to be clear, I intend on using the algorithm on Earth only and so I don't need to accommodate a variable radius.

It doesn't have to be hugely accurate (+/-20% is fine) and it'll only be used to calculate bounding boxes over small distances (no more than 150km). So I'm happy to sacrifice some accuracy for an efficient algorithm. Any help is much appreciated.

Edit: I should have been clearer, I really am after a square, not a circle. I understand that the distance between the center of a square and various points along the square's perimeter is not a constant value like it is with a circle. I guess what I mean is a square where if you draw a line from the center to any one of the four points on the perimeter that results in a line perpendicular to a side of the perimeter, then those 4 lines have the same length.

Bryce Thomas
  • 10,479
  • 26
  • 77
  • 126
  • This question is just straight trigonometry; I'd suggest removing the java and algorithm tags. (If that's even possible.) – Kevin Bourrillion Nov 07 '09 at 03:12
  • The spherical nature of the earth needs to be considered. What answer do you want if the input position is the North pole? – MarkJ Nov 09 '09 at 11:39
  • If you have a Python implementation I don't see why you couldn't convert it to Java. The algorithm should be the same. – Alexandru Luchian Nov 06 '09 at 18:01
  • The best Python example I found had been marked as untested by the author. I imagine given that I'm after something that is fairly forgiving on accuracy that a different algorithm might be appropriate also. – Bryce Thomas Nov 06 '09 at 18:21
  • It's just a formula you need; I'm sure someone will provide it for you and hopefully give you the reference link. – Kevin Bourrillion Nov 07 '09 at 03:16

6 Answers6

59

I wrote an article about finding the bounding coordinates:

http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates

The article explains the formulae and also provides a Java implementation. (It also shows why IronMan's formula for the min/max longitude is inaccurate.)

  • Thank you for writing this; I especially appreciate the SQL implementation details. – Dave Jarvis Jun 25 '10 at 02:17
  • 2
    I am confused by why these spherical approximations are taken as the standard solution. There is an abominable amount of trigonometric functions involved here. Plus, all GPS devices produce lat/long coordinates based on the WGS84 ellipsoid -- You simply cannot use a spherical approximation if you want complete accuracy... Adapting those equations in the spherical approximations to account for the ellipsoid will result in an unimaginable monstrocity. The method must then be to perform calculations and arc length integrals in 3D space, it seems. – Steven Lu Feb 13 '14 at 07:51
  • @StephanKlein The bounding box is a box that for sure covers all coordinates in distance but not all coordinates it covers are equal or below the distance. The bounding box is like a rectangle around a circle. It covers all points within the circle, but near the edges it also covers some points beyond the circle. But only with a box (a rect) you can do simple `<` and `>` comparison within a SQL statement that's why the full SQL select shown on the page first tests if a point is within the bounding box (fast!) and then it checks if it really is equal or below the desired distance (slow!) – Mecki Dec 11 '15 at 12:40
  • @StevenLu What exactly is your problem with the calculations on that page? They seem very accurate to me and there isn't much more trigonometric than in the answers of IronMan or susheel, except for for one call to `sin()` and one call to `asin()` - the rest is identical for calculating the bounding box. And the exact distance calculation used is the standard formula that everyone uses for that. So I really see no reason why this reply should not be the accepted answer, after all it is the most accurate one. – Mecki Dec 11 '15 at 12:46
  • 1
    @Mecki sorry I wasnt being completely clear and it took me a bit of reflecting to remember the original reason for my complaint. Basically the "issue" is that computationally for a computer the trig operations are not very fast (compared to simpler operations like greater than/less than). This question asks for a basic bounding box check. OP even specifies with multiple sentences that he does not want extreme accuracy, he wants speed, in fact. So using a crap load of trig on each test point to get an exact result for whether or not it is within a given distance to another point is overkill – Steven Lu Dec 12 '15 at 01:00
  • 1
    What OP wants to do (as I did in my last project I worked on in which we worked with lat/longs) is to apply a basic approximation that is super fast: You need two constants, the kilometers in a degree of latitude, and the kilometers in a degree of longitude, at the equator. Then just use this to figure out *within what amount of degrees* your points should pass the check. Yes this approximation degrades near the poles. It shouldn't matter much for most applications, though. Geographical points of interest usually are very far from the poles. (actually i did oversimplify it, its a bit harder) – Steven Lu Dec 12 '15 at 01:03
  • 1
    it's the difference between having a marginal amount of false positives (even less than you'd expect, if your screen area being displayed is rectangular which is almost always is, in fact OP *also* explicitly points this out) while doing an AABB check per point vs having zero false positives doing 3 or like 9 trig function calls per point. the tricky bit is that the nearer you get to the poles the larger the range of longitudes you have to test to safely cover your box – Steven Lu Dec 12 '15 at 01:05
  • 1
    I will say this. As long as someone actually reads the linked article they should be able to appreciate how to implement this the super accurate way and also to think about how to tweak it and relax the condition to make the code run fewer trig functions. Sometimes you can guarantee that there is always enough time to run tons and tons of trig functions for each query point. – Steven Lu Dec 12 '15 at 01:19
  • On testing this with a distance of 1 km, I got the haversine distance between lat_min,lon_min and lat_max_lon_max as 2.82 km; Why is it not 2; Also the distance between lat,lon and these points is 1.41 instead of closer to 1 ; Gist here https://gist.github.com/alexcpn/f95ae83a7ee0293a5225 – Alex Punnen Jan 16 '16 at 02:33
  • This calculates a very large bounding rectangle for coordinates in Australia. That is not acceptable since my app has thousands of users in Melbourne. I found the Drupal Earth API from Rochester IOT work better – Sacky San Aug 30 '17 at 01:08
16
double R = 6371;  // earth radius in km

double radius = 50; // km

double x1 = lon - Math.toDegrees(radius/R/Math.cos(Math.toRadians(lat)));

double x2 = lon + Math.toDegrees(radius/R/Math.cos(Math.toRadians(lat)));

double y1 = lat + Math.toDegrees(radius/R);

double y2 = lat - Math.toDegrees(radius/R);

Although I would also recommend JTS.

takrl
  • 6,356
  • 3
  • 60
  • 69
IronMan
  • 161
  • 3
5
import com.vividsolutions.jts.geom.Envelope;

...
Envelope env = new Envelope(centerPoint.getCoordinate());
env.expandBy(distance_in_degrees); 
...

Now env contains your envelope. It's not actually a "square" (whatever that means on the surface of a sphere), but it should do.

You should note that the distance in degrees will depend on the latitude of the center point. At the equator, 1 degree of latitude is about 111km, but in New York, it's only about 75km.

The really cool thing is that you can toss all your points into a com.vividsolutions.jts.index.strtree.STRtree and then use it to quickly calculate points inside that Envelope.

Bart Kiers
  • 166,582
  • 36
  • 299
  • 288
novalis
  • 1,112
  • 6
  • 12
3

All of the previous answers are only partially correct. Specially in region like Australia, they always include pole and calculate a very large rectangle even for 10kms.

Specially the algorithm by Jan Philip Matuschek at http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates#UsingIndex included a very large rectangle from (-37, -90, -180, 180) for almost every point in Australia. This hits a large users in database and distance have to be calculated for all of the users in almost half the country.

I found that the Drupal API Earth Algorithm by Rochester Institute of Technology works better around pole as well as elsewhere and is much easier to implement.

https://www.rit.edu/drupal/api/drupal/sites%21all%21modules%21location%21earth.inc/7.54

Use earth_latitude_range and earth_longitude_range from the above algorithm for calculating bounding rectangle

Here is the implementation is Java

    /**
 * Get bouding rectangle using Drupal Earth Algorithm
 * @see https://www.rit.edu/drupal/api/drupal/sites%21all%21modules%21location%21earth.inc/7.54
 * @param lat
 * @param lng
 * @param distance
 * @return
 */
default BoundingRectangle getBoundingRectangleDrupalEarthAlgo(double lat, double lng, int distance) {
    lng = Math.toRadians(lng);
    lat = Math.toRadians(lat);
    double radius = earth_radius(lat);
    List<Double> retLats = earth_latitude_range(lat, radius, distance);
    List<Double> retLngs = earth_longitude_range(lat, lng, radius, distance);
    return new BoundingRectangle(retLats.get(0), retLats.get(1), retLngs.get(0), retLngs.get(1));
}


/**
 * Calculate latitude range based on earths radius at a given point
 * @param latitude
 * @param longitude
 * @param distance
 * @return
 */
default List<Double> earth_latitude_range(double lat, double radius, double distance) {
      // Estimate the min and max latitudes within distance of a given location.

      double angle = distance / radius;
      double minlat = lat - angle;
      double maxlat = lat + angle;
      double rightangle = Math.PI / 2;
      // Wrapped around the south pole.
      if (minlat < -rightangle) {
        double overshoot = -minlat - rightangle;
        minlat = -rightangle + overshoot;
        if (minlat > maxlat) {
          maxlat = minlat;
        }
        minlat = -rightangle;
      }
      // Wrapped around the north pole.
      if (maxlat > rightangle) {
        double overshoot = maxlat - rightangle;
        maxlat = rightangle - overshoot;
        if (maxlat < minlat) {
          minlat = maxlat;
        }
        maxlat = rightangle;
      }
      List<Double> ret = new ArrayList<>();
      ret.add((minlat));
      ret.add((maxlat));
      return ret;
    }

/**
 * Calculate longitude range based on earths radius at a given point
 * @param lat
 * @param lng
 * @param earth_radius
 * @param distance
 * @return
 */
default List<Double> earth_longitude_range(double lat, double lng, double earth_radius, int distance) {
      // Estimate the min and max longitudes within distance of a given location.
      double radius = earth_radius * Math.cos(lat);

      double angle;
      if (radius > 0) {
        angle = Math.abs(distance / radius);
        angle = Math.min(angle, Math.PI);
      }
      else {
        angle = Math.PI;
      }
      double minlong = lng - angle;
      double maxlong = lng + angle;
      if (minlong < -Math.PI) {
        minlong = minlong + Math.PI * 2;
      }
      if (maxlong > Math.PI) {
        maxlong = maxlong - Math.PI * 2;
      }

      List<Double> ret = new ArrayList<>();
      ret.add((minlong));
      ret.add((maxlong));
      return ret;
    }

/**
 * Calculate earth radius at given latitude
 * @param latitude
 * @return
 */
default Double earth_radius(double latitude) {
      // Estimate the Earth's radius at a given latitude.
      // Default to an approximate average radius for the United States.
      double lat = Math.toRadians(latitude);

      double x = Math.cos(lat) / 6378137.0;
      double y = Math.sin(lat) / (6378137.0 * (1 - (1 / 298.257223563)));

      //Make sure earth's radius is in km , not meters
      return (1 / (Math.sqrt(x * x + y * y)))/1000;
    }

And use the distance calculation formula documented by google maps to calculate distance

https://developers.google.com/maps/solutions/store-locator/clothing-store-locator#outputting-data-as-xml-using-php

To search by kilometers instead of miles, replace 3959 with 6371. For (Lat, Lng) = (37, -122) and a Markers table with columns lat and lng, the formula is:

SELECT id, ( 3959 * acos( cos( radians(37) ) * cos( radians( lat ) ) * cos( radians( lng ) - radians(-122) ) + sin( radians(37) ) * sin( radians( lat ) ) ) ) AS distance FROM markers HAVING distance < 25 ORDER BY distance LIMIT 0 , 20;
Sacky San
  • 1,535
  • 21
  • 26
  • Observe the earth_radius code . Make sure earth's radius is in km , not meters: return (1 / (Math.sqrt(x * x + y * y)))/1000; – Sacky San Aug 30 '17 at 02:47
1

Based on IronMan response:

/**
 * Calculate the lat and len of a square around a point.
 * @return latMin, latMax, lngMin, lngMax
 */
public static double[] calculateSquareRadius(double lat, double lng, double radius) {
    double R = 6371;  // earth radius in km
    double latMin = lat - Math.toDegrees(radius/R);
    double latMax = lat + Math.toDegrees(radius/R);
    double lngMin = lng - Math.toDegrees(radius/R/Math.cos(Math.toRadians(lat)));
    double lngMax = lng + Math.toDegrees(radius/R/Math.cos(Math.toRadians(lat)));

    return new double[] {latMin, latMax, lngMin, lngMax};
}
Laurent
  • 469
  • 3
  • 7
0

Here is a simple solution that I used to generate bounding box coordinates that I use with GeoNames citieJSON API to get nearby big cities from a gps decimal coordinate.

This is a Java method from my GitHub repository: FusionTableModifyJava

I had a decimal GPS location and I needed to find the biggest city/state "near" that location. I needed a relatively accurate bounding box to pass to the citiesJSON GeoNames webservice to get back the biggest city in that bounding box. I pass the location and the "radius" I am interested in (in km) and it gives back the north, south, east, west decimal coordinates needed to pass to citiesJSON.

(I found these resources useful in doing my research:

Calculate distance, bearing and more between Latitude/Longitude points.

Longitude - Wikipedia)

It is not super accurate but accurate enough for what I was using it for:

    // Compute bounding Box coordinates for use with Geonames API.
    class BoundingBox
    {
        public double north, south, east, west;
        public BoundingBox(String location, float km)
        {
             //System.out.println(location + " : "+ km);
            String[] parts = location.replaceAll("\\s","").split(","); //remove spaces and split on ,

            double lat = Double.parseDouble(parts[0]);
            double lng = Double.parseDouble(parts[1]);

            double adjust = .008983112; // 1km in degrees at equator.
            //adjust = 0.008983152770714983; // 1km in degrees at equator.

            //System.out.println("deg: "+(1.0/40075.017)*360.0);


            north = lat + ( km * adjust);
            south = lat - ( km * adjust);

            double lngRatio = 1/Math.cos(Math.toRadians(lat)); //ratio for lng size
            //System.out.println("lngRatio: "+lngRatio);

            east = lng + (km * adjust) * lngRatio;
            west = lng - (km * adjust) * lngRatio;
        }

    }