Как рассчитать ограничивающий прямоугольник для данного местоположения широты / долготы?

Question

Я дал местоположение, определенное широтой и долготой. Теперь я хочу вычислить ограничивающую рамку, например, в пределах. 10 километров от этой точки.

Ограничительная рамка должна быть определена как latmin, lngmin и latmax, lngmax.

Мне нужен этот материал, чтобы использовать Panoramio API .

Кто-нибудь знает формулу, как получить эти очки?

Редактировать: Ребята, я ищу формулу / функцию, которая принимает lat & lng в качестве входных данных и возвращает ограничивающий прямоугольник как latmin & lngmin и latmax & latmin. Mysql, php, c #, javascript - это нормально, но с псевдокодом все должно быть в порядке.

Редактировать: Я не ищу решение, которое показывает мне расстояние 2 точки

Answer 1

Я предлагаю локально аппроксимировать поверхность Земли как сферу с радиусом, заданным эллипсоидом WGS84 на данной широте. Я подозреваю, что точное вычисление latMin и latMax потребовало бы эллиптических функций и не дало бы заметного увеличения точности (WGS84 сам по себе является приблизительным).

Моя реализация следующая (написана на Python; я ее не тестировал):

# 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))

РЕДАКТИРОВАТЬ: следующий код преобразует (градусы, простые числа, секунды) в градусы + доли градуса и наоборот (не проверено):

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))

Answer 2

Я написал статью о поиске ограничивающих координат:

http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates

Статья объясняет формулы, а также предоставляет реализацию Java. (Это также показывает, почему формула Федерико для минимальной / максимальной долготы является неточной.)

Answer 3

Здесь я преобразовал ответ Федерико А. Рампони в C # для всех, кто заинтересован:

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));
}

Answer 4

Я написал функцию JavaScript, которая возвращает четыре координаты квадратного ограничивающего прямоугольника с учетом расстояния и пары координат:

'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()
  ];
};

Answer 5

Вы ищете формулу эллипсоида.

Лучшее место, которое я нашел для начала кодирования, основано на библиотеке Geo :: Ellipsoid из CPAN. Это дает вам основу для создания ваших тестов и сравнения ваших результатов с его результатами. Я использовал его в качестве основы для подобной библиотеки для PHP у моего предыдущего работодателя.

Geo :: Эллипсоид

Взгляните на метод location. Звоните дважды, и вы получите свой bbox.

Вы не опубликовали, какой язык вы использовали. Возможно, вам уже доступна библиотека геокодирования.

О, и если вы еще не поняли это, Google maps использует эллипсоид WGS84.

Answer 6

Поскольку мне нужна была очень приблизительная оценка, поэтому для фильтрации некоторых ненужных документов в запросе эластичного поиска я использовал следующую формулу:

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 = км требуется от данного места. Для вашего случая N = 10

Не точно, но удобно.

Answer 7

Вот простая реализация, использующая JavaScript, которая основана на преобразовании градуса широты в км, где 1 degree latitude ~ 111.2 km.

Я рассчитываю границы карты по заданной широте и долготе с шириной 10 км.

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;
}

Answer 8

Я адаптировал скрипт PHP, который нашел именно для этого. Вы можете использовать его, чтобы найти углы коробки вокруг точки (скажем, 20 км). Мой конкретный пример для API Карт Google:

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

Answer 9

Иллюстрация @Jan Philip Matuschek превосходное объяснение. (Пожалуйста, проголосуйте за его ответ, а не за это; я добавляю это, поскольку я потратил немного времени на понимание исходного ответа)

Метод ограничивающего прямоугольника оптимизации поиска ближайших соседей должен был бы получить минимальную и максимальную пары широты и долготы для точки P на расстоянии d. Все точки, которые выходят за их пределы, определенно находятся на расстоянии, большем чем d от точки. Здесь нужно отметить одну широту пересечения, что подчеркивается в объяснении Яна Филиппа Матушека. Широта пересечения не находится на широте точки P, но немного смещена от нее. Это часто пропускаемая, но важная часть в определении правильной минимальной и максимальной ограничивающей долготы для точки P для расстояния d. Это также полезно при проверке.

Расстояние haversine от (широта пересечения, долгота высокая) до (широта, долгота) P равно расстоянию d.

Суть Python здесь https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

Answer 10

Я работал над проблемой ограничивающего прямоугольника как побочной проблемой при поиске всех точек в радиусе SrcRad статической точки LAT, LONG. Было довольно много вычислений, которые используют

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

для расчета границ долготы, но я обнаружил, что это не дает ответы на все вопросы, которые были необходимы. Потому что то, что вы действительно хотите сделать, это

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

Я знаю, я знаю, что ответ должен быть таким же, но я обнаружил, что это не так. Оказалось, что, не убедившись, что сначала я выполняю (SRCrad / RadEarth), а затем делю на часть Cos, я оставлял некоторые точки местоположения.

После того, как вы получите все свои ограничивающие точки, если у вас есть функция, которая вычисляет расстояние от точки до точки по заданному лату, можно легко получить только те точки, которые находятся на определенном расстоянии от фиксированной точки. Вот что я сделал. Я знаю, что потребовалось несколько дополнительных шагов, но это помогло мне

-- 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
;