sábado, 11 de maio de 2013

MySQL & GIS & A fórmula Haversine e distância entre dois pontos.

Original post: http://anothermysqldba.blogspot.com/2013/05/mysql-gis-haversine-formula-and.html


MySQL não é o banco de dados que vem à mente em primeiro lugar quando as pessoas pensam de GIS. As bases de dados são listadas abaixo:

  • Oráculo
  • Microsoft SQL Server
  • IBM DB2
  • IBM Informix
  • PostgreSQL

http://webhelp.esri.com/arcgisserver/9.3/java/index.htm # geodatabase / types_of_geodatabases.htm
http://webhelp.esri.com/arcgisserver/9.3/java/index.htm # geodatabases/determ-1479045992.htm


MySQL vem trabalhando com GIS já há algum tempo, ele vai voltar para o MySQL 4.1:
O post de Mark Maunder acima é um grande olhar para GIS. Este blog vai focar mais em fórmulas, mas eu queria salientar bom post de Marcos.

O uso de um índice espacial trabalha com o mecanismo de armazenamento MyISAM, por isso, se você está no MySQL 5.5 ou superior manter isso em mente como InnoDB é o mecanismo de armazenamento padrão.

> CREATE TABLE geom (
-> lat float(10,7) NOT NULL,
-> lon float(10,7) NOT NULL,
-> g GEOMETRY NOT NULL,
-> SPATIAL INDEX(g)
-> ) ENGINE=Innodb;
ERROR 1464 (HY000): The used table type doesn't support SPATIAL indexes
> CREATE TABLE geom (
-> lat float(10,7) NOT NULL,
-> lon float(10,7) NOT NULL,
-> g GEOMETRY NOT NULL,
-> SPATIAL INDEX(g)
-> ) ENGINE=MyISAM; 

Assim, com o projeto do esquema de Mark I povoaram a latitude / longitude de algumas cidades na China.
Estes dados foram recolhidos através de aquihttp://www.infoplease.com/ipa/A0001769.html

CREATE TABLE china (
cityname varchar(50) DEFAULT NULL,
lat float(10,7) NOT NULL,
lon float(10,7) NOT NULL,
g GEOMETRY NOT NULL,
SPATIAL INDEX(g)
) ENGINE=MyISAM; 
INSERT INTO china VALUES ('Beijing', 39.55, 116.25, GeomFromText('POINT(39.55 116.25)'));
INSERT INTO china VALUES ('Canton', 23.7, 113.15, GeomFromText('POINT(23.7 113.15)'));
INSERT INTO china VALUES ('Chongqing', 29.46, 106.34, GeomFromText('POINT(29.46, 106.34)'));
INSERT INTO china VALUES ('Hong Kong', 22.20, 114.11, GeomFromText('POINT( 22.20 114.11)'));
INSERT INTO china VALUES ('Shanghai', 31.10, 121.28, GeomFromText('POINT(31.10 121.28)'));

Só para ter certeza que tudo funcionou ....

select lat, lon from china;
+------------+-------------+
| lat | lon |
+------------+-------------+
| 39.5499992 | 116.2500000 |
| 23.7000008 | 113.1500015 |
| 22.2000008 | 114.1100006 |
| 31.1000004 | 121.2799988 |
+------------+-------------+


Assim, permite-nos verificar a distância de Pequim para Hong Kong.
Uma rápida olhada on-line nos diz que é 1963 km, mas vamos verificar nossos dados.


Agora era hora de cavar o mundo de fórmulas de distância .. Eu não sou um GIS DBA focado, eu admito que não há problema. Eu gostava mais do leque de respostas para as fórmulas matemáticas que eu encontrei, caso contrário, eu teria usado as informações da tabela acima mais. Em vez disso, é apenas uma referência para você. Tomei a olhar para as formas amplamente debatido para calcular a distância entre dois Lat / Long pontos. Você vai encontrar um monte de diferentes funções, procedimentos e etc on-line para calcular isso.

Primeiro eu definir algumas variáveis, porque eu queria testar as fórmulas

SET @lat1 =39.55;
SET @long1 =116.25;
SET @lat2 =22.20;
SET @long2 =114.11; 

Por exemplo:

Este site tem uma função para a Distância. Eu adicionei-a abaixo para garantir que você pode cortar e colá-lo com as opções de delimitador lugar para você. Mas será que isso funciona? BTW eu também atualizou o raio do valor de terra para 3959.
http://www.sqlexamples.info/SPAT/mysql_distance.htm

delimiter //
CREATE FUNCTION fn_distance
(p_x1 FLOAT, p_y1 FLOAT, p_x2 FLOAT, p_y2 FLOAT)
RETURNS FLOAT
DETERMINISTIC
BEGIN
DECLARE v_dist FLOAT;
DECLARE A FLOAT; DECLARE B FLOAT;
DECLARE C FLOAT; DECLARE D FLOAT;
/*
returns distance calculation between two points in LAT-LONG coordinates
*/

SET v_dist = 0;

-- convert to radians
SET A = p_x1 / 57.29577951;
SET B = p_y1 / 57.29577951;
SET C = p_x2 / 57.29577951;
SET D = p_y2 / 57.29577951;

IF (A = C && B = D) THEN
SET v_dist = 0;
ELSEIF ((sin(A)*sin(C)+cos(A)*cos(C)*cos(B - D)) > 1) THEN
SET v_dist = 3959 * acos(1);
ELSE
SET v_dist = 3959 *acos(sin(A)*sin(C) + cos(A)*cos(C)*cos(B - D));
END IF;

SET v_dist = v_dist * 1.609;

/* return distance in km. */
RETURN v_dist;

END;
//
delimiter ;

> SELECT fn_distance (@lat1, @long1, @lat2 , @long2) AS dist_km;
+-------------------+
| dist_km |
+-------------------+
| 1939.5457763671875 |
+-------------------+

Outro teste consulta mostra 

> SELECT ( GLength( LineString(( PointFromWKB( POINT( @lat1, @long1 ))), ( PointFromWKB( POINT( @lat2, @long2 ) ))))) * 100 AS distance;
+--------------------+
| distance |
+--------------------+
| 1748.1478770401545 |
+--------------------+ 

Contudo uma outra encontrada aqui http://www.posteet.com/view/1555 : 

set log_bin_trust_function_creators=TRUE;
DELIMITER |
CREATE FUNCTION GeoDistKM( lat1 FLOAT, lon1 FLOAT, lat2 FLOAT, lon2 FLOAT ) RETURNS float
BEGIN
DECLARE pi, q1, q2, q3 FLOAT;
DECLARE rads FLOAT DEFAULT 0;
SET pi = PI();
SET lat1 = lat1 * pi / 180;
SET lon1 = lon1 * pi / 180;
SET lat2 = lat2 * pi / 180;
SET lon2 = lon2 * pi / 180;
SET q1 = COS(lon1-lon2);
SET q2 = COS(lat1-lat2);
SET q3 = COS(lat1+lat2);
SET rads = ACOS( 0.5*((1.0+q1)*q2 - (1.0-q1)*q3) );
RETURN 6378.388 * rads;
END;
|
DELIMITER ;
select geodistkm( set log_bin_trust_function_creators=TRUE;
DELIMITER |
CREATE FUNCTION GeoDistKM( lat1 FLOAT, lon1 FLOAT, lat2 FLOAT, lon2 FLOAT ) RETURNS float
BEGIN
DECLARE pi, q1, q2, q3 FLOAT;
DECLARE rads FLOAT DEFAULT 0;
SET pi = PI();
SET lat1 = lat1 * pi / 180;
SET lon1 = lon1 * pi / 180;
SET lat2 = lat2 * pi / 180;
SET lon2 = lon2 * pi / 180;
SET q1 = COS(lon1-lon2);
SET q2 = COS(lat1-lat2);
SET q3 = COS(lat1+lat2);
SET rads = ACOS( 0.5*((1.0+q1)*q2 - (1.0-q1)*q3) );
RETURN 6378.388 * rads;
END;
|
DELIMITER ;
select geodistkm( @ lat1, @ long1, @ lat2, @ long2 ) as distance;
+----------------------------------------+
| distance |
+----------------------------------------+
| 1942.0909423828125 |
+----------------------------------------+ ) as distance;
+----------------------------------------+
| distance |
+----------------------------------------+
| 1942.0909423828125 |
+----------------------------------------+ 


No entanto, isso ainda é errado. Porque a distância de Pequim para Hong King é 1224,9 milhas com seus 1.963 km, de acordo com este site:http://www.timeanddate.com/worldclock/distances.html?n=102 . Assim, os primeiros e últimos resultados são muito próximos. 

Então, depois de perder muito tempo, sim, eu admito isso, eu ainda queria um resultado usando fórmulas matemáticas que eu poderia facilmente comparar. 

O grande círculo distância d entre dois pontos com coordenadas {lat1, lon1} e {lat2, lon2} é dada por: 
d = acos (sin (lat1) * sin (lat2) + cos (lat1) * cos (lat2) * cos (lon1-lon2)) 
> SELECT ACOS(
-> SIN(@lat1) * SIN(@lat2) + COS(@lat1) * COS(@lat2) * COS( @long1 - @long2 )
-> ) *1000 as Aviation_forumula_DISTANCE;
+----------------------------+
| Aviation_forumula_DISTANCE |
+----------------------------+
| 1923.0473470093848 |
+----------------------------+


Outra fórmula é a fórmula Haversine 





> SELECT 3956* 2 * ASIN ( SQRT (POWER(SIN((@lat1 - @lat2)*pi()/180 / 2),2) + COS(@lat1 * pi()/180) * COS(@lat2 *pi()/180) * POWER(SIN((@long1 - @long2) *pi()/180 / 2), 2) ) ) as Haversine_Formula_distance;
+----------------------------+
| Haversine_Formula_distance |
+----------------------------+
| 1204.5222518763514 |
+----------------------------+ 




Mais uma vez, você verá resultados muito diferentes. Eu esperava melhores resultados com afórmula Haversine. 
Então, depois de brincar com estas fórmulas Fui com a função GeoDistKM porque parece-me ser a mais próxima da fórmula Haversine que eu acredito ser a fórmula correta para usar quando implementado corretamente. Isto é seguido de um segundo para fechar a fórmula I Aviação escreveu com base nas informações de Williams. 

Enquanto 1942 não é a 1963 resultado Juntei através da pesquisa on-line, que é dizer que eles calcularam corretamente também. A curvatura da terra e do lat, lon aproximando juntos nos pólos vai permitir alguns erros em fórmulas diferentes. Então eu vou ficar com esse momento: 

select geodistkm( @ lat1, @ long1, @ lat2, @ long2 ) as distance;
+----------------------------------------+
| distance |
+----------------------------------------+
| 1942.0909423828125 |
+----------------------------------------+ ) as distance;
+----------------------------------------+
| distance |
+----------------------------------------+
| 1942.0909423828125 |
+----------------------------------------+ 











Postado por



às






















Marcadores:
,
,
,
,