¿Cómo se calcula el cuadro delimitador para una ubicación lat / lng determinada?

He dado una ubicación definida por latitud y longitud. Ahora quiero calcular un cuadro delimitador dentro de, por ejemplo, 10 kilómetros de ese punto.

El cuadro delimitador debe definirse como latmin, lngmin y latmax, lngmax.

Necesito esto para usar la API de panoramio .

¿Alguien sabe la fórmula de cómo obtener esos puntos?

Edit: Chicos estoy buscando una fórmula / función que tome lat & lng como entrada y devuelva un cuadro delimitador como latmin & lngmin y latmax & latmin. Mysql, php, c #, javascript está bien pero también el pseudocódigo debería estar bien.

Editar: No estoy buscando una solución que me muestre la distancia de 2 puntos

Sugiero aproximar localmente la superficie de la Tierra como una esfera con radio dado por el elipsoide WGS84 en la latitud dada. Sospecho que el cálculo exacto de latMin y latMax requeriría funciones elípticas y no produciría un aumento apreciable en la precisión (WGS84 es en sí mismo una aproximación).

Mi implementación sigue (está escrita en Python, no la he probado):

# degrees to radians def deg2rad(degrees): return math.pi*degrees/180.0 # radians to degrees def rad2deg(radians): return 180.0*radians/math.pi # 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] def WGS84EarthRadius(lat): # http://en.wikipedia.org/wiki/Earth_radius An = WGS84_a*WGS84_a * math.cos(lat) Bn = WGS84_b*WGS84_b * math.sin(lat) Ad = WGS84_a * math.cos(lat) Bd = WGS84_b * math.sin(lat) return math.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 def 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*math.cos(lat) latMin = lat - halfSide/radius latMax = lat + halfSide/radius lonMin = lon - halfSide/pradius lonMax = lon + halfSide/pradius return (rad2deg(latMin), rad2deg(lonMin), rad2deg(latMax), rad2deg(lonMax)) 

EDITAR: El siguiente código convierte (grados, primos, segundos) a grados + fracciones de un grado, y viceversa (no probado):

 def dps2deg(degrees, primes, seconds): return degrees + primes/60.0 + seconds/3600.0 def deg2dps(degrees): intdeg = math.floor(degrees) primes = (degrees - intdeg)*60.0 intpri = math.floor(primes) seconds = (primes - intpri)*60.0 intsec = round(seconds) return (int(intdeg), int(intpri), int(intsec)) 

Escribí un artículo sobre cómo encontrar las coordenadas de límite:

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 mínima / máxima es inexacta).

Aquí he convertido la respuesta de Federico A. Ramponi a C # para cualquiera que esté 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 de 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  * @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á buscando una fórmula elipsoide.

El mejor lugar que he encontrado para comenzar a codificar se basa en la biblioteca Geo :: Ellipsoid de CPAN. Le brinda 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

Eche un vistazo al método de location . Llámalo dos veces y obtendrás tu bbox.

No publicaste el idioma que estabas usando. Es posible que ya haya una biblioteca de geoencoding disponible para usted.

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

Adapte un script PHP que encontré para hacer solo esto. Puede usarlo para encontrar las esquinas de una caja alrededor de un punto (por ejemplo, a 20 km de distancia). 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

Como necesitaba una estimación aproximada, para filtrar algunos documentos innecesarios en una consulta de búsqueda elástica, 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 desde la ubicación dada. Para su caso N = 10

No es preciso pero práctico.

Aquí hay 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 desde una latitud y longitud dada 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 la excelente explicación de @Jan Philip Matuschek. (Por favor, vote-up su respuesta, no esto, estoy agregando esto ya que me tomé un poco de tiempo para entender la respuesta original)

La técnica del cuadro delimitador para optimizar el hallazgo de los vecinos más cercanos necesitaría derivar los pares mínimo y máximo de latitud y longitud para un punto P a la distancia d. Todos los puntos que quedan fuera de estos están definitivamente a una distancia mayor que d desde el punto. Una cosa a tener en cuenta aquí es el cálculo de la latitud de intersección como se destaca en la explicación de Jan Philip Matuschek. La latitud de la intersección no está en la latitud del punto P, sino que está ligeramente desviada de la misma. Esta es una parte que a menudo se pierde, pero es importante para determinar 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.

La esencia de Python aquí https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

enter image description here

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

 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é que esto no daba 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. Al parecer, al no asegurarme de que estaba haciendo (SRCrad / RadEarth) Primero y luego dividir por la parte Cos, estaba dejando fuera algunos puntos de ubicación.

Después de obtener todos los puntos del cuadro delimitador, si tiene una función que calcula la distancia de punto a punto dada latitud y longitud, es fácil obtener solo los puntos que tienen un radio de distancia determinado desde el punto fijo. Aquí esta lo que hice. Sé que tomó algunos pasos extra 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 ir a la página web de panoramio y luego abrir el mapa mundial desde el sitio web de panoramio. Luego, vaya a la ubicación especificada, lo que requiere latitud y longitud.

Luego, encontró la latitud y la 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 API JavaScript de Panoramio crea un cuadro delimitador alrededor de un par latitud / longitud y luego regresa todas las fotos con esos límites.

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

No se muestra en el humor de composición. Se muestra después de la publicación.

 

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

  

Gracias @Fedrico A. por la implementación de Phyton, lo 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 debe ser reemplazado por una clase, lo he usado solo para compartirlo aquí.

Todas las respuestas anteriores son parcialmente correctas . Especialmente en regiones como Australia, siempre incluyen polos y calculan un rectángulo muy grande incluso para 10kms.

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

Descubrí que el algoritmo de la Tierra Drupal API del Rochester Institute of Technology funciona mejor en los polos y 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

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

Y use 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 verificació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), } } 
Intereting Posts