Содержание


Выполнять преобразование координат стало проще

Пособие по преобразованию между различными системами координат

Comments

Службы позиционирования, включая навигацию на основе GPS и картографические сайты, такие как Google Maps и Yahoo! Maps, становятся популярными среди пользователей. Множество организаций уже используют сервисы, связанные с использованием информации о географических координатах, и еще больше компаний пойдут по этому пути, как только осознают преимущества и потенциал подобных приложений. В 2006 году аналитическая компания Gartner отметила, что "приложения, связанные с позиционированием, станут массовыми в течение следующих двух-пяти лет" и что уже "значительное число организаций развернули мобильные бизнес-приложения, использующие позиционирование". (В разделе Ресурсы приведена ссылка на этот отчет.)

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

Две различные системы координат

Прежде чем погрузиться в код, представленный в этой статье, необходимо обсудить координатные системы, для поддержки которых этот код предназначен: известная система широты и долготы и универсальная поперечная проекция Меркатора (Universal Transverse Mercator - UTM). Также нужно коснуться военной системы координат (Military Grid Reference System - MGRS), которая основана на UTM.

Широта и долгота

Системы широты и долготы, вероятно, самый известный способ определения географических координат. В ней местоположение представляется двумя числами. Широта - это угол от центра земли к некоторой параллели на земной поверхности. Долгота - это угол от центра земли к некоторому меридиану на земной поверхности. Широта и долгота могут быть выражены в десятичных градусах (DD) или градусах, минутах и секунда (DMS); в последнем случае получаются числа в таком формате - 49°30'00" S 12°30'00" E. Этот формат обычно используется в GPS-навигаторах.

Земля разделена экватором (0° широты) на северное и южное полушария, и нулевым меридианом (0° долготы), воображаемой линией от северного к южному полюсу, которая проходитчерез город Гринвич в Великобритании и делит планету на восточное и западное полушарие. Диапазон широт в северном полушарии от 0 до 90 градусов, а в южном полушарии - от 0 до -90 градусов. Диапазон восточного полушария от 0 до 180 градусов, а западного полушария от 0 до -180 градусов.

Например, точка с координатами 61.44, 25.40 (в формате DD) или 61°26'24''N, 25°23'60''E (в формате DMS) находится в южной Финляндии. А точка с координатами -47.04, -73.48 (DD) или 47°02'24''S, 73°28'48''W (DMS) находится в южном Чили. На рисунке 1 приведено изображение Земли с перекрывающимися линиями параллелей и меридианов:

Рисунок 1. Земля с показанными линиями параллелей и меридианов
Рисунок 1. Земля с показанными линиями параллелей и меридианов
Рисунок 1. Земля с показанными линиями параллелей и меридианов

Дополнительную информацию можно найти в разделе Ресурсы.

Поперечная проекция Меркатора

Система координат UTM - это метод, использующий сетку для определения координат. Система UTM делит Землю на 60 зон, каждая из которых основана на поперечной проекции Меркатора. Проекция карты в картографии - это способ представить двумерную неровную поверхность на плоскости, как обычную карту. На рисунке 2 приведена поперечная проекция Меркатора:

Рисунок 2. Поперечная проекция Меркатора
Рисунок 2. Поперечная проекция Меркатора
Рисунок 2. Поперечная проекция Меркатора

Зоны долготы в UTM пронумерованы от 1 до 60; все зоны кроме двух, о которых будет рассказано позже, имеют ширину 6° от востока до запада. Зоны долготы полностью покрывают поверхность Земли между широтами 80°S и 84°N.

Также есть 20 зон широты, каждая в 8° высотой; эти зоны пронумерованы от C до X, буквы I и O пропущены. Зоны A, B, Y и Z находятся за пределами этой системы, они покрывают Арктику и Антарктику. На рисунке 3 приведены UTM зоны для Европы. На этом рисунке видны две нестандартные зоны долготы: зона 32V расширена для покрытия всей южной Норвегии, а зона 31V сокращена, чтобы покрывать только водное пространство.

Рисунок 3. Зоны UTM в Европе
Рисунок 3. Зоны UTM в Европе
Рисунок 3. Зоны UTM в Европе

Координаты в UTM представлены в формате зона долготы зона широты восточное склонение северное склонение, где восточное склонение - это проекционное расстояние от центрального меридиана зоны долготы, северное склонение - это проекционное расстояние от экватора. Значения восточного и северного склонений задаются в метрах. Например, координаты широты/долготы 61.44, 25.40 в UTM представлены как 35 V 414668 6812844; координаты широты/долготы -47.04, -73.48 соответствуют координатам 18 G 615471 4789269 в UTM.

В разделе Ресурсы приведена дополнительная информация по UTM и поперечной проекции Меркатора.

Военная система координат

Система координат MGRS - стандарт, используемый военными НАТО. MGRS основана на UTM и дальше делит каждую зону на квадраты 100х100 км. Эти квадраты идентифицируются двухбуквенным кодом, первая буква - восточно-западная позиция в зоне долготы, а вторая буква - северо-южная позиция.

Например, координата в UTM 35 V 414668 6812844 эквивалента координате MGRS 35VMJ1466812844. Точность координаты в MGRS задается с точностью в один метр и представлена с помощью 15 символов, где последние 10 символов представляют значения восточного и северного склонений в указанной сетке. В MGRS координаты могут быть представлены 15 символами, как в прошлом примере, или 13, 11, 9 или 7 символами; представленные таким образом значения будут соответственно иметь точность 1, 10, 100, 1000 или 10000 метров.

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

Преобразования координат

Чтобы определить ширину и долготу - координаты местоположения на Земле, как минимум нужно иметь возможность видеть звезды или Солнце, иметь секстан и часы, показывающие время на меридиане GMT. Можно определить широту из угла между небесным телом и горизонтом, а долготу можно вычислить из вращения Земли. Эта статья не погружается в подобные подробности, но в разделе Ресурсы о них можно узнать больше. Вместо этого, предположим, что у нас уже есть координаты в формате DD, DMS или UTM.

Преобразование десятичных градусов в градусы/минуты/секунды и обратно

Крайне просто преобразовать координаты из формата DD в DMS. Ниже приведена формула для подобного преобразования:

DD: dd.ff
DMS: dd mm ss

dd=dd
mm.gg=60*ff
ss=60*gg

В этом примере gg - это дробная часть вычисления. Отрицательная широта означает местоположение в южном полушарии (S), а отрицательная долгота - в западном полушарии (W). Например, предположим, что имеются координаты в формате DD - 61.44, 25.40. Их можно преобразовать следующим образом:

lat dd=61
lat mm.gg=60*0.44=26.4
lat ss=60*0.4=24

Далее:

lon dd=25
lon mm.gg=60*0.40=24.0
lon ss=60*0.0=0

Таким образом, в формате DMS получаем следующие координаты - 61°26'24''N 25°24'00''E.

Ниже приведена формула для перехода от DMS к DD:

DD: dd.ff
DMS: dd mm ss

dd.ff=dd + mm/60 + ss/3600

Напомним, что места, расположенные в южном полушарии (S), имеют отрицательную широту, а места в западном полушарии (W) имеют отрицательную долготу.

Теперь выполним преобразование DMS координат 47°02'24''S, 73°28'48''W в формат DD:

lat dd.ff= - (47 + 2/60 + 24/3600 )=-47.04
lon dd.ff= - (73 + 28/60 + 48/3600)=-73.48

Таким образом, координаты в DD равны -47.04, -73.48.

Преобразование долготы/широты в UTM и обратно

В отличие от десятичных координат, которые можно определить с помощью хронометра и часов, координаты UTM невозможно определить без помощи вычислений. Хотя эти вычисления ничто иное, как простая тригонометрия и алгебра, формулы у них достаточно сложные. Если ознакомиться со статьей "The Universal Grids: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)" (ссылка на неё приведена в разделе Ресурсы), то станет понятно, что я имею ввиду.

Формулы для преобразования UTM здесь не приводятся, но исходный код в следующем разделе немного освещает эту проблему, а в разделе Ресурсы есть ссылки на дополнительную информацию.

Преобразование координат в Java коде

В этом разделе приведен исходный код класса библиотеки, который выполняет преобразования координат между десятичными градусами и UTM. Этот Java-класс называется com.ibm.util.CoordinateConversion и я хотел создать отдельный класс с методами для преобразования. Этот класс включает внутренние классы, которые на самом деле выполняют преобразования, и при необходимости эти классы могут быть вынесены из класса CoordinateConversion с помощью рефакторинга, чтобы создать пакет библиотеки или добавить классы в существующий пакет. Преобразования, выполняемые этим классом, имеют точность в 1 метр.

В исходном коде CoordinateConversion около 750 строк, так что в этой статье он представлен не полностью. Существенные методы описаны в следующих разделах, а полностью код доступен в разделе Материалы для скачивания.

Класс CoordinateConversion

CoordinateConversion - это главный класс, объекты которого создаются при необходимости выполнения преобразования координат. В листинге 1 приведены существенные public методы вместе с внутренними private классами, входящими в класс CoordinateConversion:

Листинг 1. CoordinateConversion
public class CoordinateConversion
{

  public CoordinateConversion()
  {

  }
  
  public double[] utm2LatLon(String UTM)
  {
    UTM2LatLon c = new UTM2LatLon();
    return c.convertUTMToLatLong(UTM);
  }

  public String latLon2UTM(double latitude, double longitude)
  {
    LatLon2UTM c = new LatLon2UTM();
    return c.convertLatLonToUTM(latitude, longitude);

  }

  //..реализация пропущена

  private class LatLon2UTM
  {
    public String convertLatLonToUTM(double latitude, double longitude)
    {
      //..реализация пропущена
    }
    //..реализация пропущена
  }
  
  private class LatLon2MGRUTM extends LatLon2UTM
  {
    public String convertLatLonToMGRUTM(double latitude, double longitude)
    {
      //..реализация пропущена
    }
    //..реализация пропущена
  }


  private class MGRUTM2LatLon extends UTM2LatLon
  {
    public double[] convertMGRUTMToLatLong(String mgrutm)
    {
      //..реализация пропущена
    }
    //..реализация пропущена
  }
 
 
  private class UTM2LatLon
  {
    public double[] convertUTMToLatLong(String UTM)
    {
      //..реализация пропущена
    }
    //..реализация пропущена
  }

  private class Digraphs
  {
    //используется для получения двухбуквенных кодов
    //при преобразовании от долготы/широты к MGRS
    //..реализация пропущена
  }

  private class LatZones
  {
    //включает методы для определения зон широты
	//..реализация пропущена
  }

В следующем разделе подробно рассматриваются преобразования между долготой/широтой и UTM.

Преобразование от широты/долготы к UTM

Координаты преобразуются от широты/долготы в UTM с помощью метода String latLon2UTM(double latitude, double longitude). Реализация этого метода создает экземпляр внутреннего класса LatLon2UTM c = new LatLon2UTM(); и возвращает координаты UTM в виде 15-символьной строки с точностью 1 метр. Реализация методов класса LatLon2UTM приведена в листинге 2:

Листинг 2. public String convertLatLonToUTM(double latitude, double longitude)
public String convertLatLonToUTM(double latitude, double longitude)
{
  validate(latitude, longitude);
  String UTM = "";

  setVariables(latitude, longitude);

  String longZone = getLongZone(longitude);
  LatZones latZones = new LatZones();
  String latZone = latZones.getLatZone(latitude);

  double _easting = getEasting();
  double _northing = getNorthing(latitude);

  UTM = longZone + " " + latZone + " " + ((int) _easting) + " "+ ((int) _northing);

  return UTM;
}

Этот метод выполняет преобразование, вызывая различные методы для получения зоны широты и зоны долготы, вычисления восточного и северного склонения и т.д. Входные данные проверяются в методе validate(), если выражение (latitude < -90.0 || latitude > 90.0 || longitude < -180.0 || longitude >= 180.0) принимает значение true, то сбрасывается исключительная ситуация IllegalArgumentException.

Метод setVariables() в листинге 3 устанавливает различные переменные, требующиеся для вычисления преобразований (дополнительная информация представлена в статье "The Universal Grids", ссылка на которую приведена в разделе Ресурсы:

Листинг 3. protected void setVariables(double latitude, double longitude)
protected void setVariables(double latitude, double longitude)
{
  latitude = degreeToRadian(latitude);
  rho = equatorialRadius * (1 - e * e) / POW(1 - POW(e * SIN(latitude), 2), 3 / 2.0);

  nu = equatorialRadius / POW(1 - POW(e * SIN(latitude), 2), (1 / 2.0));

  double var1;
  if (longitude < 0.0)
  {
    var1 = ((int) ((180 + longitude) / 6.0)) + 1;
  }
  else
  {
    var1 = ((int) (longitude / 6)) + 31;
  }
  double var2 = (6 * var1) - 183;
  double var3 = longitude - var2;
  p = var3 * 3600 / 10000;

  S = A0 * latitude - B0 * SIN(2 * latitude) + C0 * SIN(4 * latitude) - D0
      * SIN(6 * latitude) + E0 * SIN(8 * latitude);

  K1 = S * k0;
  K2 = nu * SIN(latitude) * COS(latitude) * POW(sin1, 2) * k0 * (100000000) / 2;
  K3 = ((POW(sin1, 4) * nu * SIN(latitude) * Math.pow(COS(latitude), 3)) / 24)
      * (5 - POW(TAN(latitude), 2) + 9 * e1sq * POW(COS(latitude), 2) + 4
      * POW(e1sq, 2) * POW(COS(latitude), 4))
      * k0
      * (10000000000000000L);

  K4 = nu * COS(latitude) * sin1 * k0 * 10000;

  K5 = POW(sin1 * COS(latitude), 3) * (nu / 6)
      * (1 - POW(TAN(latitude), 2) + e1sq * POW(COS(latitude), 2)) * k0
      * 1000000000000L;

  A6 = (POW(p * sin1, 6) * nu * SIN(latitude) * POW(COS(latitude), 5) / 720)
      * (61 - 58 * POW(TAN(latitude), 2) + POW(TAN(latitude), 4) + 270
      * e1sq * POW(COS(latitude), 2) - 330 * e1sq
      * POW(SIN(latitude), 2)) * k0 * (1E+24);

}

Метод getLongZone() в листинге 4 и класс LatZones, доступный в исходном коде, используются, чтобы узнать зону долготы и зону широты. Зона долготы вычисляется по параметру longitude, а зоны широты обычно представлены как константы, с помощью массива в классе LatZones.

Листинг 4. protected String getLongZone(double longitude)
protected String getLongZone(double longitude)
{
  double longZone = 0;
  if (longitude < 0.0)
  {
    longZone = ((180.0 + longitude) / 6) + 1;
  }
  else
  {
    longZone = (longitude / 6) + 31;
  }
  String val = String.valueOf((int) longZone);
  if (val.length() == 1)
  {
    val = "0" + val;
  }
  return val;
}

Метод getNorthing() в листинге 5 и метод getEasting() в листинге 6 вычисляют значения северного и восточного склонения. Оба метода используют переменные, установленные в методе setVariables() из листинга 3.

Листинг 5. protected double getNorthing(double latitude)
protected double getNorthing(double latitude)
{
  double northing = K1 + K2 * p * p + K3 * POW(p, 4);
  if (latitude < 0.0)
  {
    northing = 10000000 + northing;
  }
  return northing;
}
Листинг 6. protected double getEasting()
protected double getEasting()
{
  return 500000 + (K4 * p + K5 * POW(p, 3));
}

В листинге 7 приведены несколько примеров результатов работы программы, включая координаты в формате широта/долгота и соответствующие им UTM координаты:

Листинг 7. Тестовые преобразования от широты/долготы к значениям в UTM
( 0.0000    0.0000  )     "31 N 166021 0"
( 0.1300   -0.2324  )     "30 N 808084 14385"
(-45.6456   23.3545 )     "34 G 683473 4942631"
(-12.7650  -33.8765 )     "25 L 404859 8588690"
(-80.5434  -170.6540)     "02 C 506346 1057742"
( 90.0000   177.0000)     "60 Z 500000 9997964"
(-90.0000  -177.0000)     "01 A 500000 2035"
( 90.0000    3.0000 )     "31 Z 500000 9997964"
( 23.4578  -135.4545)     "08 Q 453580 2594272"
( 77.3450   156.9876)     "57 X 450793 8586116"
(-89.3454  -48.9306 )     "22 A 502639 75072"

Преобразование от UTM к широте/долготе

Преобразование от координат в формате UTM к широте/долготе выполняется несколько проще, чем обратный процесс. В статье "The Universal Grids" в разделе Ресурсы) приведены формулы преобразований. В листинге 8 приведен код метода convertUTMToLatLong(). Этот метод возвращает массив double значений, где первый элемент (с индексом массива [0]) - это широта, а второй элемент (с индексом массива [1]) - это долгота. Так как строковый параметр содержит координаты UTM с точностью до 1 метра, то и координаты в широте/долготе будут иметь такую же точность.

Листинг 8. public double[] convertUTMToLatLong(String UTM)
public double[] convertUTMToLatLong(String UTM)
{
  double[] latlon = { 0.0, 0.0 };
  String[] utm = UTM.split(" ");
  zone = Integer.parseInt(utm[0]);
  String latZone = utm[1];
  easting = Double.parseDouble(utm[2]);
  northing = Double.parseDouble(utm[3]);
  String hemisphere = getHemisphere(latZone);
  double latitude = 0.0;
  double longitude = 0.0;

  if (hemisphere.equals("S"))
  {
    northing = 10000000 - northing;
  }
  setVariables();
  latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / Math.PI;

  if (zone > 0)
  {
    zoneCM = 6 * zone - 183.0;
  }
  else
  {
    zoneCM = 3.0;
  }

  longitude = zoneCM - _a3;
  if (hemisphere.equals("S"))
  {
    latitude = -latitude;
  }

  latlon[0] = latitude;
  latlon[1] = longitude;
  return latlon;

}

Метод convertUTMToLatLong() разбивает координаты UTM во входном строковом параметре, которые имеют формат 34 G 683473 4942631, и использует метод getHemisphere() для определения полушария, где находится место с указанными координатами. Определить полушарие просто: зоны широты A, C, D, E, F, G, H, J, K, L и M находятся в южном полушарии, а остальные зоны находятся в северном полушарии.

Метод setVariables(), показанный в листинге 9, устанавливает значения переменных, требуемых для вычисления, и затем немедленно вычисляют широту. Долгота вычисляется с использованием зоны долготы.

Листинг 9. protected void setVariables()
protected void setVariables()
{
  arc = northing / k0;
  mu = arc
      / (a * (1 - POW(e, 2) / 4.0 - 3 * POW(e, 4) / 64.0 - 5 * POW(e, 6) / 256.0));

  ei = (1 - POW((1 - e * e), (1 / 2.0)))
      / (1 + POW((1 - e * e), (1 / 2.0)));

  ca = 3 * ei / 2 - 27 * POW(ei, 3) / 32.0;

  cb = 21 * POW(ei, 2) / 16 - 55 * POW(ei, 4) / 32;
  cc = 151 * POW(ei, 3) / 96;
  cd = 1097 * POW(ei, 4) / 512;
  phi1 = mu + ca * SIN(2 * mu) + cb * SIN(4 * mu) + cc * SIN(6 * mu) + cd
      * SIN(8 * mu);

  n0 = a / POW((1 - POW((e * SIN(phi1)), 2)), (1 / 2.0));

  r0 = a * (1 - e * e) / POW((1 - POW((e * SIN(phi1)), 2)), (3 / 2.0));
  fact1 = n0 * TAN(phi1) / r0;

  _a1 = 500000 - easting;
  dd0 = _a1 / (n0 * k0);
  fact2 = dd0 * dd0 / 2;

  t0 = POW(TAN(phi1), 2);
  Q0 = e1sq * POW(COS(phi1), 2);
  fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * POW(dd0, 4) / 24;

  fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0
          * Q0)
          * POW(dd0, 6) / 720;

  lof1 = _a1 / (n0 * k0);
  lof2 = (1 + 2 * t0 + Q0) * POW(dd0, 3) / 6.0;
  lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * POW(Q0, 2) + 8 * e1sq + 24 * POW(t0, 2))
          * POW(dd0, 5) / 120;
  _a2 = (lof1 - lof2 + lof3) / COS(phi1);
  _a3 = _a2 * 180 / Math.PI;

}

Метод setVariables() использует значения восточного и северного склонений для установки требуемых переменных. Эти переменные принадлежат обоим классам и устанавливаются в методе convertUTMToLatLong(String UTM) из листинга 8.

Другие методы

Исходный код также содержит другие public и private методы и классы. Например, туда включены методы и классы для преобразования координат между широтой/долготой и MGRS вместе с вспомогательными методами, выполняющими преобразование градусов в радианы и наоборот, и различными математическими операциями, такими как POW, SIN, COS, and TAN.

Заключение

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

Мне потребовалось выполнять преобразования между широтой и долготой, UTM и MGRS, так что я выполнил базовые исследования и реализовал эти преобразования в Java-классе. Для меня разработка заняла несколько часов, и я надеюсь, что и другие также смогут сэкономить несколько часов для других задач и сочтут полезным использование класса CoordinateConversion в собственной работе.


Ресурсы для скачивания


Похожие темы


Комментарии

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

static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=40
Zone=Технология Java
ArticleID=302110
ArticleTitle=Выполнять преобразование координат стало проще
publish-date=04212008