domingo, 12 de mayo de 2013

MySQL y GIS y La fórmula Haversine y la distancia entre dos puntos.

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


MySQL no es la base de datos que viene a la mente primero cuando la gente piensa de los SIG. Las bases de datos se enumeran a continuación son:

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

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


MySQL ha estado trabajando con GIS por algún tiempo ahora, que se remonta a MySQL 4.1:
El post de Mark Maunder arriba es un gran aspecto en GIS. Este blog se centrará más en las fórmulas, pero quería señalar buen puesto de Mark.

El uso de un índice espacial funciona con el motor de almacenamiento MyISAM, así que si usted está en MySQL 5.5 o superior que tenlo en cuenta como InnoDB es el motor de almacenamiento por defecto.

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

Así que con el diseño de Mark esquema que poblaron la latitud / longitud de algunas ciudades de China.
Estos datos fueron recogidos a través de estahttp://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ólo para asegurarse de que todo funcionaba ....

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


Así que nos permite comprobar la distancia entre Pekín y Hong Kong.
Un vistazo en línea nos dice que es 1.963 kilometros, pero nos dejaron registrarnos nuestros datos.


Ahora era el momento de profundizar en el mundo de las fórmulas distancia .. Yo no soy un DBA centrado GIS, admito que no hay problema. Me gustó más de la gama de respuestas a las fórmulas matemáticas que encontré, de lo contrario habría utilizado la información de la tabla anterior más. En cambio, es sólo una referencia para usted. Que tomé en vistazo a las formas ampliamente debatidos para calcular la distancia entre dos puntos / latitud y longitud. Usted encontrará una gran cantidad de funciones, procedimientos y etc en línea para calcular esto.

Primero me puse algunas variables porque quería probar las fórmulas

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

Por ejemplo:

Este sitio web cuenta con una función de la distancia. He añadido lo siguiente para asegurarse de que usted puede cortar y pegar con las opciones delimitadores en el lugar para usted. Pero, ¿funciona? Por cierto que también se actualiza el radio del valor de la tierra 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 |
+-------------------+

Otra prueba de consulta muestra 

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

Otro encontrar aquí 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 |
+----------------------------------------+ 


Sin embargo, esto sigue siendo incorrecto. Debido a la distancia de Beijing a Hong King es 1224,9 millas con sus 1.963 kilometros, de acuerdo al sitio web:http://www.timeanddate.com/worldclock/distances.html?n=102 . Así que la primera y la última los resultados están muy cerca. 

Así que después de perder mucho tiempo, sí admito que yo todavía quería un resultado utilizando fórmulas matemáticas que podía comparar fácilmente en contra. 

La distancia ortodrómica d entre dos puntos con coordenadas {lat1, lon1} y {lat2, lon2} viene dado 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 |
+----------------------------+


Otra fórmula es la 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 |
+----------------------------+ 




Una vez que vea resultados muy diferentes. Me esperaba mejores resultados con la fórmulaHaversine. 
Así que después de jugar con estas fórmulas Fui con la función GeoDistKM porque me parece ser el más cercano a la fórmula Haversine que yo creo que es la fórmula correcta para utilizar cuando se aplica correctamente. Esto es seguido por un segundo cercano a la fórmula de Aviación escribí sobre la base de información de Williams. 

Mientras que 1942 no es el resultado de 1963 reuní a través de la búsqueda en línea, que quiere decir que ellos calcularon correctamente también. La curvatura de la tierra y el lat, lon acercarse juntos en los polos permitirá algunos errores en diferentes fórmulas. Así que voy a seguir con esta actualmente: 

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











Publicado por



en






















Etiquetas:
,
,
,
,