¿Cómo calcular el cuadro delimitador para una ubicación lat/lng dada?

Escribí un artículo sobre la búsqueda de las coordenadas delimitadoras:

Http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates

El artículo explica las fórmulas y también proporciona una implementación de Java. (También muestra por qué la fórmula de Federico para la longitud min/max es inexacta.)

Aquí he convertido la respuesta de Federico A. Ramponi a C# para cualquiera interesado:

public class MapPoint
{
    public double Longitude { get; set; } // In Degrees
    public double Latitude { get; set; } // In Degrees
}

public class BoundingBox
{
    public MapPoint MinPoint { get; set; }
    public MapPoint MaxPoint { get; set; }
}        

// Semi-axes of WGS-84 geoidal reference
private const double WGS84_a = 6378137.0; // Major semiaxis [m]
private const double WGS84_b = 6356752.3; // Minor semiaxis [m]

// 'halfSideInKm' is the half length of the bounding box you want in kilometers.
public static BoundingBox GetBoundingBox(MapPoint point, double halfSideInKm)
{            
    // Bounding box surrounding the point at given coordinates,
    // assuming local approximation of Earth surface as a sphere
    // of radius given by WGS84
    var lat = Deg2rad(point.Latitude);
    var lon = Deg2rad(point.Longitude);
    var halfSide = 1000 * halfSideInKm;

    // Radius of Earth at given latitude
    var radius = WGS84EarthRadius(lat);
    // Radius of the parallel at given latitude
    var pradius = radius * Math.Cos(lat);

    var latMin = lat - halfSide / radius;
    var latMax = lat + halfSide / radius;
    var lonMin = lon - halfSide / pradius;
    var lonMax = lon + halfSide / pradius;

    return new BoundingBox { 
        MinPoint = new MapPoint { Latitude = Rad2deg(latMin), Longitude = Rad2deg(lonMin) },
        MaxPoint = new MapPoint { Latitude = Rad2deg(latMax), Longitude = Rad2deg(lonMax) }
    };            
}

// degrees to radians
private static double Deg2rad(double degrees)
{
    return Math.PI * degrees / 180.0;
}

// radians to degrees
private static double Rad2deg(double radians)
{
    return 180.0 * radians / Math.PI;
}

// Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
private static double WGS84EarthRadius(double lat)
{
    // http://en.wikipedia.org/wiki/Earth_radius
    var An = WGS84_a * WGS84_a * Math.Cos(lat);
    var Bn = WGS84_b * WGS84_b * Math.Sin(lat);
    var Ad = WGS84_a * Math.Cos(lat);
    var Bd = WGS84_b * Math.Sin(lat);
    return Math.Sqrt((An*An + Bn*Bn) / (Ad*Ad + Bd*Bd));
}

Escribí una función JavaScript que devuelve las cuatro coordenadas de un cuadro delimitador cuadrado, dada una distancia y un par de coordenadas:

'use strict';

/**
 * @param {number} distance - distance (km) from the point represented by centerPoint
 * @param {array} centerPoint - two-dimensional array containing center coords [latitude, longitude]
 * @description
 *   Computes the bounding coordinates of all points on the surface of a sphere
 *   that has a great circle distance to the point represented by the centerPoint
 *   argument that is less or equal to the distance argument.
 *   Technique from: Jan Matuschek <http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates>
 * @author Alex Salisbury
*/

getBoundingBox = function (centerPoint, distance) {
  var MIN_LAT, MAX_LAT, MIN_LON, MAX_LON, R, radDist, degLat, degLon, radLat, radLon, minLat, maxLat, minLon, maxLon, deltaLon;
  if (distance < 0) {
    return 'Illegal arguments';
  }
  // helper functions (degrees<–>radians)
  Number.prototype.degToRad = function () {
    return this * (Math.PI / 180);
  };
  Number.prototype.radToDeg = function () {
    return (180 * this) / Math.PI;
  };
  // coordinate limits
  MIN_LAT = (-90).degToRad();
  MAX_LAT = (90).degToRad();
  MIN_LON = (-180).degToRad();
  MAX_LON = (180).degToRad();
  // Earth's radius (km)
  R = 6378.1;
  // angular distance in radians on a great circle
  radDist = distance / R;
  // center point coordinates (deg)
  degLat = centerPoint[0];
  degLon = centerPoint[1];
  // center point coordinates (rad)
  radLat = degLat.degToRad();
  radLon = degLon.degToRad();
  // minimum and maximum latitudes for given distance
  minLat = radLat - radDist;
  maxLat = radLat + radDist;
  // minimum and maximum longitudes for given distance
  minLon = void 0;
  maxLon = void 0;
  // define deltaLon to help determine min and max longitudes
  deltaLon = Math.asin(Math.sin(radDist) / Math.cos(radLat));
  if (minLat > MIN_LAT && maxLat < MAX_LAT) {
    minLon = radLon - deltaLon;
    maxLon = radLon + deltaLon;
    if (minLon < MIN_LON) {
      minLon = minLon + 2 * Math.PI;
    }
    if (maxLon > MAX_LON) {
      maxLon = maxLon - 2 * Math.PI;
    }
  }
  // a pole is within the given distance
  else {
    minLat = Math.max(minLat, MIN_LAT);
    maxLat = Math.min(maxLat, MAX_LAT);
    minLon = MIN_LON;
    maxLon = MAX_LON;
  }
  return [
    minLon.radToDeg(),
    minLat.radToDeg(),
    maxLon.radToDeg(),
    maxLat.radToDeg()
  ];
};

Estás buscando una fórmula elipsoidal.

El mejor lugar que he encontrado para comenzar a codificar está basado en la biblioteca Geo::Ellipsoid de CPAN. Le da una línea de base para crear sus pruebas y comparar sus resultados con sus resultados. Lo utilicé como base para una biblioteca similar para PHP en mi empleador anterior.

Geo:: Elipsoide

Echa un vistazo al método location. Llámalo dos veces y tendrás tu bbox.

No publicaste qué idioma estábamos consumiendo. Es posible que ya haya una biblioteca de geocodificación disponible para usted.

Ah, y si aún no lo has descubierto, Google Maps utiliza el elipsoide WGS84.

Adapté un script PHP que encontré para hacer precisamente esto. Puede usarlo para encontrar las esquinas de una caja alrededor de un punto (digamos, a 20 km). Mi ejemplo específico es para la API de Google Maps:

Http://www.richardpeacock.com/blog/2011/11/draw-box-around-coordinate-google-maps-based-miles-or-kilometers

Dado que necesitaba una estimación muy aproximada, por lo que para filtrar algunos documentos innecesarios en una consulta de elasticsearch, empleé la siguiente fórmula:

Min.lat = Given.Lat - (0.009 x N)
Max.lat = Given.Lat + (0.009 x N)
Min.lon = Given.lon - (0.009 x N)
Max.lon = Given.lon + (0.009 x N)

N = kms requerido de la ubicación dada. Para su caso N = 10

No es preciso pero práctico.

Aquí está una implementación simple usando javascript que se basa en la conversión del grado de latitud a kms donde 1 degree latitude ~ 111.2 km.

Estoy calculando los límites del mapa a partir de una latitud y longitud dadas con un ancho de 10 km.

function getBoundsFromLatLng(lat, lng){
     var lat_change = 10/111.2;
     var lon_change = Math.abs(Math.cos(lat*(Math.PI/180)));
     var bounds = { 
         lat_min : lat - lat_change,
         lon_min : lng - lon_change,
         lat_max : lat + lat_change,
         lon_max : lng + lon_change
     };
     return bounds;
}

Ilustración de @Jan Philip Matuschek excelente explicación.(Por favor vote su respuesta, no esta; estoy agregando esto ya que me tomé un poco de tiempo para entender la respuesta original)

La técnica de caja delimitadora de optimización de encontrar vecinos más cercanos necesitaría derivar la latitud mínima y máxima,pares de longitud, para un punto P a distancia d . Todos los puntos que caen fuera de estos están definitivamente a una distancia mayor que d del punto. Una cosa a tener en cuenta aquí es la cálculo de la latitud de intersección como se destaca en la explicación de Jan Philip Matuschek. La latitud de intersección no está en la latitud del punto P sino ligeramente desplazada de él. Esta es una parte a menudo perdida pero importante en la determinación de la longitud límite mínima y máxima correcta para el punto P para la distancia d.Esto también es útil en la verificación.

La distancia de haversine entre (latitud de intersección,longitud alta) a (latitud,longitud) de P es igual a la distancia d.

Python gist here https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

Estaba trabajando en el problema del cuadro delimitador como un problema secundario para encontrar todos los puntos dentro del radio SrcRad de un LAT estático, punto LARGO. Ha habido bastantes cálculos que utilizan

maxLon = $lon + rad2deg($rad/$R/cos(deg2rad($lat)));
minLon = $lon - rad2deg($rad/$R/cos(deg2rad($lat)));

Para calcular los límites de longitud, pero encontré esto para no dar todas las respuestas que se necesitaban. Porque lo que realmente quieres hacer es

(SrcRad/RadEarth)/cos(deg2rad(lat))

Lo sé, sé que la respuesta debería ser la misma, pero descubrí que no lo era. (SRCrad/RadEarth) Primero y luego dividiendo por la parte Cos estaba dejando fuera algunos puntos de ubicación.

Después de obtener todos sus puntos de caja delimitadora, si tiene una función que calcula la distancia de Punto a Punto dada lat, long es fácil obtener solo aquellos puntos que están a un cierto radio de distancia desde el punto fijo. Esto es lo que hice. Sé que tomó algunos pasos adicionales, pero me ayudó

-- GLOBAL Constants
gc_pi CONSTANT REAL := 3.14159265359;  -- Pi

-- Conversion Factor Constants
gc_rad_to_degs          CONSTANT NUMBER := 180/gc_pi; -- Conversion for Radians to Degrees 180/pi
gc_deg_to_rads          CONSTANT NUMBER := gc_pi/180; --Conversion of Degrees to Radians

lv_stat_lat    -- The static latitude point that I am searching from 
lv_stat_long   -- The static longitude point that I am searching from 

-- Angular radius ratio in radians
lv_ang_radius := lv_search_radius / lv_earth_radius;
lv_bb_maxlat := lv_stat_lat + (gc_rad_to_deg * lv_ang_radius);
lv_bb_minlat := lv_stat_lat - (gc_rad_to_deg * lv_ang_radius);

--Here's the tricky part, accounting for the Longitude getting smaller as we move up the latitiude scale
-- I seperated the parts of the equation to make it easier to debug and understand
-- I may not be a smart man but I know what the right answer is... :-)

lv_int_calc := gc_deg_to_rads * lv_stat_lat;
lv_int_calc := COS(lv_int_calc);
lv_int_calc := lv_ang_radius/lv_int_calc;
lv_int_calc := gc_rad_to_degs*lv_int_calc;

lv_bb_maxlong := lv_stat_long + lv_int_calc;
lv_bb_minlong := lv_stat_long - lv_int_calc;

-- Now select the values from your location datatable 
SELECT *  FROM (
SELECT cityaliasname, city, state, zipcode, latitude, longitude, 
-- The actual distance in miles
spherecos_pnttopntdist(lv_stat_lat, lv_stat_long, latitude, longitude, 'M') as miles_dist    
FROM Location_Table 
WHERE latitude between lv_bb_minlat AND lv_bb_maxlat
AND   longitude between lv_bb_minlong and lv_bb_maxlong)
WHERE miles_dist <= lv_limit_distance_miles
order by miles_dist
;

Es muy simple solo tienes que ir al sitio web de panoramio y luego abrir el Mapa del mundo desde el sitio web de panoramio.Luego vaya a la ubicación especificada que se requiere latitud y longitud.

Luego encontró latitud y longitud en la barra de direcciones, por ejemplo, en esta dirección.

Http://www.panoramio.com/map#lt=32.739485&ln=70.491211&z=9&k=1&a=1&tab=1&pl=all

Lt = 32.739485 = > latitud ln = 70.491211 = > longitud

Este widget de la API de JavaScript de Panoramio crea un cuadro delimitador alrededor de un par lat/long y luego devolver todas las fotos con en esos límites.

Otro tipo de widget API JAVASCRIPT de Panoramio en el que también puede cambiar el color de fondo con ejemplo y código está aquí.

No se muestra en la composición mood.It mostrar después de publicar.

<div dir="ltr" style="text-align: center;" trbidi="on">
<script src="https://ssl.panoramio.com/wapi/wapi.js?v=1&amp;hl=en"></script>
<div id="wapiblock" style="float: right; margin: 10px 15px"></div>
<script type="text/javascript">
var myRequest = {
  'tag': 'kahna',
  'rect': {'sw': {'lat': -30, 'lng': 10.5}, 'ne': {'lat': 50.5, 'lng': 30}}
};
  var myOptions = {
  'width': 300,
  'height': 200
};
var wapiblock = document.getElementById('wapiblock');
var photo_widget = new panoramio.PhotoWidget('wapiblock', myRequest, myOptions);
photo_widget.setPosition(0);
</script>
</div>

Aquí he convertido la respuesta de Federico A. Ramponi a PHP si alguien está interesado:

<?php
# deg2rad and rad2deg are already within PHP

# Semi-axes of WGS-84 geoidal reference
$WGS84_a = 6378137.0;  # Major semiaxis [m]
$WGS84_b = 6356752.3;  # Minor semiaxis [m]

# Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
function WGS84EarthRadius($lat)
{
    global $WGS84_a, $WGS84_b;

    $an = $WGS84_a * $WGS84_a * cos($lat);
    $bn = $WGS84_b * $WGS84_b * sin($lat);
    $ad = $WGS84_a * cos($lat);
    $bd = $WGS84_b * sin($lat);

    return sqrt(($an*$an + $bn*$bn)/($ad*$ad + $bd*$bd));
}

# Bounding box surrounding the point at given coordinates,
# assuming local approximation of Earth surface as a sphere
# of radius given by WGS84
function boundingBox($latitudeInDegrees, $longitudeInDegrees, $halfSideInKm)
{
    $lat = deg2rad($latitudeInDegrees);
    $lon = deg2rad($longitudeInDegrees);
    $halfSide = 1000 * $halfSideInKm;

    # Radius of Earth at given latitude
    $radius = WGS84EarthRadius($lat);
    # Radius of the parallel at given latitude
    $pradius = $radius*cos($lat);

    $latMin = $lat - $halfSide / $radius;
    $latMax = $lat + $halfSide / $radius;
    $lonMin = $lon - $halfSide / $pradius;
    $lonMax = $lon + $halfSide / $pradius;

    return array(rad2deg($latMin), rad2deg($lonMin), rad2deg($latMax), rad2deg($lonMax));
}
?>

Gracias @Fedrico A. por la implementación de Phyton, la he portado a una clase de categoría Objective C. Aquí está:

#import "LocationService+Bounds.h"

//Semi-axes of WGS-84 geoidal reference
const double WGS84_a = 6378137.0; //Major semiaxis [m]
const double WGS84_b = 6356752.3; //Minor semiaxis [m]

@implementation LocationService (Bounds)

struct BoundsLocation {
    double maxLatitude;
    double minLatitude;
    double maxLongitude;
    double minLongitude;
};

+ (struct BoundsLocation)locationBoundsWithLatitude:(double)aLatitude longitude:(double)aLongitude maxDistanceKm:(NSInteger)aMaxKmDistance {
    return [self boundingBoxWithLatitude:aLatitude longitude:aLongitude halfDistanceKm:aMaxKmDistance/2];
}

#pragma mark - Algorithm 

+ (struct BoundsLocation)boundingBoxWithLatitude:(double)aLatitude longitude:(double)aLongitude halfDistanceKm:(double)aDistanceKm {
    double radianLatitude = [self degreesToRadians:aLatitude];
    double radianLongitude = [self degreesToRadians:aLongitude];
    double halfDistanceMeters = aDistanceKm*1000;


    double earthRadius = [self earthRadiusAtLatitude:radianLatitude];
    double parallelRadius = earthRadius*cosl(radianLatitude);

    double radianMinLatitude = radianLatitude - halfDistanceMeters/earthRadius;
    double radianMaxLatitude = radianLatitude + halfDistanceMeters/earthRadius;
    double radianMinLongitude = radianLongitude - halfDistanceMeters/parallelRadius;
    double radianMaxLongitude = radianLongitude + halfDistanceMeters/parallelRadius;

    struct BoundsLocation bounds;
    bounds.minLatitude = [self radiansToDegrees:radianMinLatitude];
    bounds.maxLatitude = [self radiansToDegrees:radianMaxLatitude];
    bounds.minLongitude = [self radiansToDegrees:radianMinLongitude];
    bounds.maxLongitude = [self radiansToDegrees:radianMaxLongitude];

    return bounds;
}

+ (double)earthRadiusAtLatitude:(double)aRadianLatitude {
    double An = WGS84_a * WGS84_a * cosl(aRadianLatitude);
    double Bn = WGS84_b * WGS84_b * sinl(aRadianLatitude);
    double Ad = WGS84_a * cosl(aRadianLatitude);
    double Bd = WGS84_b * sinl(aRadianLatitude);
    return sqrtl( ((An * An) + (Bn * Bn))/((Ad * Ad) + (Bd * Bd)) );
}

+ (double)degreesToRadians:(double)aDegrees {
    return M_PI*aDegrees/180.0;
}

+ (double)radiansToDegrees:(double)aRadians {
    return 180.0*aRadians/M_PI;
}



@end

Lo he probado y parece estar funcionando bien. Struct BoundsLocation debería ser reemplazado por una clase, la he usado solo para compartirla aquí.

Todas las respuestas anteriores son solo parcialmente correctas. Especialmente en la región como Australia, siempre incluyen polo y calcular un rectángulo muy grande, incluso para 10 km.

Especialmente el algoritmo de Jan Philip Matuschek en http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates#UsingIndex incluye un rectángulo muy grande de (-37, -90, -180, 180) para casi todos los puntos de Australia. Esto golpea a usuarios grandes en base de datos y la distancia tiene que ser calculada para todos los usuarios en casi la mitad del país.

Encontré que el Algoritmo de Drupal API Earth del Rochester Institute of Technology funciona mejor alrededor de pole, así como en otros lugares y es mucho más fácil de implementar.

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

Use earth_latitude_range y earth_longitude_range del algoritmo anterior para calcular el rectángulo delimitador

Y utilizar la fórmula de cálculo de distancia documentada por google maps para calcular la distancia

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

Para buscar por kilómetros en lugar de millas, reemplace 3959 con 6371. Para (Lat, Lng) = (37, -122) y una tabla de marcadores con columnas lat y lng , la fórmula es:

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;

Lea mi respuesta detallada en https://stackoverflow.com/a/45950426/5076414

Aquí está la respuesta de Federico Ramponi en Go. Nota: sin comprobación de errores: (

import (
    "math"
)

// Semi-axes of WGS-84 geoidal reference
const (
    // Major semiaxis (meters)
    WGS84A = 6378137.0
    // Minor semiaxis (meters)
    WGS84B = 6356752.3
)

// BoundingBox represents the geo-polygon that encompasses the given point and radius
type BoundingBox struct {
    LatMin float64
    LatMax float64
    LonMin float64
    LonMax float64
}

// Convert a degree value to radians
func deg2Rad(deg float64) float64 {
    return math.Pi * deg / 180.0
}

// Convert a radian value to degrees
func rad2Deg(rad float64) float64 {
    return 180.0 * rad / math.Pi
}

// Get the Earth's radius in meters at a given latitude based on the WGS84 ellipsoid
func getWgs84EarthRadius(lat float64) float64 {
    an := WGS84A * WGS84A * math.Cos(lat)
    bn := WGS84B * WGS84B * math.Sin(lat)

    ad := WGS84A * math.Cos(lat)
    bd := WGS84B * math.Sin(lat)

    return math.Sqrt((an*an + bn*bn) / (ad*ad + bd*bd))
}

// GetBoundingBox returns a BoundingBox encompassing the given lat/long point and radius
func GetBoundingBox(latDeg float64, longDeg float64, radiusKm float64) BoundingBox {
    lat := deg2Rad(latDeg)
    lon := deg2Rad(longDeg)
    halfSide := 1000 * radiusKm

    // Radius of Earth at given latitude
    radius := getWgs84EarthRadius(lat)

    pradius := radius * math.Cos(lat)

    latMin := lat - halfSide/radius
    latMax := lat + halfSide/radius
    lonMin := lon - halfSide/pradius
    lonMax := lon + halfSide/pradius

    return BoundingBox{
        LatMin: rad2Deg(latMin),
        LatMax: rad2Deg(latMax),
        LonMin: rad2Deg(lonMin),
        LonMax: rad2Deg(lonMax),
    }
}