Algo

Calcul d'un orthodromie et recouvrement d'une sphère.

Recherche :

Mot : Pseudo : Filtrer
Bas de page
Auteur	Sujet : Calcul d'un orthodromie et recouvrement d'une sphère.

Profil supprimé

Bonjour,
Je cherche à savoir si une position sur la terre se trouve à une distance <= à une distance D d'un autre point sur la terre.
J'ai donc vu la fonction orthodromique, mais j'ai un problème a priori, je trouve des objet positionné à une distance supérieur à D.

Voici ma fonction, c'est avec Ada, mais c'est pas bien compliqué, à part peut-être le sens des multiplication et des division pas 100000.0 qui ne sont pas moins obscure pour moi.

D'abord, je convertie des coordonnée GPS en degré minute seconde en degrés décimaux.

Pour la longitude :

Code :

function Longitude (X : in Longitude_Coordonate_Type) return Decimal_Longitude_Type is                        
   begin
      
      if X.Degr_Lon < 0.0 then
         return -Decimal_Longitude_Type(abs(X.Degr_Lon) + (abs(X.Minu_Lon) / 60.0) + (abs(X.Seco_Lon) / 3600.0));
      else
            return Decimal_Longitude_Type(abs(X.Degr_Lon) + (abs(X.Minu_Lon) / 60.0) + (abs(X.Seco_Lon) / 3600.0));
      end if;
   end Longitude;

Pour la latitude.

Code :

function latitude (Y : in Latitude_Coordonate_Type) return Decimal_Latitude_Type is
   begin
      
      if Y.Degr_Lat < 0.0 then
         return -Decimal_Latitude_Type(abs(Y.Degr_Lat) + (abs(Y.Minu_Lat) / 60.0) + (abs(Y.Seco_Lat) / 3600.0));
      else
         return Decimal_Latitude_Type(abs(Y.Degr_Lat) + (abs(Y.Minu_Lat) / 60.0) + (abs(Y.Seco_Lat) / 3600.0));
      end if;         
   end Latitude;

En suite j'applique la fonction d'orthodromie et je compare les distances.

Code :

function Is_in(Target : in Gps_Coordonates_Type; item : in Gps_Coordonates_Type; Distance : in Float) return Boolean is
            
                              
      --  Formule À Appliquer Pour Calculer Une Orthodromie : Cos P = (Sin LA . Sin LB) + (Cos LA . Cos LB . Cos G)      
      LA, LB, G: Float := 0.0;
      P : Float := 0.0;
      --    LA Correspond À La Latitude Du Point A.
      --    LB Correspond À La Latitude Du Point B.
      --    G Correspond À La Différence De Longitude Entre Le Point A Et Le Point B.       
      Cos_Produit : Float := 0.0;
      sin_Produit : Float := 0.0;
      Sum_Produit : Float := 0.0;
   begin   
      
      La := Float'Rounding(Latitude(Target.Latitude) * 100000.0) / 100000.0;
      
      Lb := Float'Rounding(Latitude(Item.Latitude) * 100000.0) / 100000.0;
      
      G  :=  Float'Rounding(Longitude(item.Longitude) * 100000.0) / 100000.0 - Float'Rounding(Longitude(Target.Longitude) * 100000.0) / 100000.0;
      
      Cos_Produit := (Float'Rounding((Cos(La) * Cos(Lb) * Cos(G)) * 100000.0) /100000.0);
      --Text_Io.Put_Line("cos_produit"  & Float'Image(cos_produit));
      sin_Produit := (Float'Rounding((Sin(La) * Sin(Lb))  * 100000.0) /100000.0);      
      --Text_Io.Put_Line("sin_produit"  & Float'Image(sin_produit));
      Sum_Produit := (Cos_Produit + Sin_produit);      
      --Text_Io.Put_Line("sum_produit"  & Float'Image(sum_produit));
      P := Arccos(Sum_Produit);
      
      P := P * 60.0;
      
      if  P - Distance > 0.0 then
         --Text_Io.Put_Line("False" );
         return False;
      else
         Text_Io.Put_Line("TRUE" );
         Text_Io.Put_Line("P = " & Float'Image(P));      
         return True;
      end if;
   end Is_in;

Qui renvoie vrai pour beaucoup et suutout n'importe quoi.

Si vous pouviez jeter un oeil si vous avez le temps. S'il vous plaît ?
Merci.

Message édité par Profil supprimé le 30-09-2013 à 20:57:06

Publicité

Profil supprimé

Par exemple,

Pour les coordonnée suivante

Wopr coordonates
Long 1.0° 0.0' 0.0", Lat 85.0° 0.0' 0.0"
Target coordonates
Long -118.0° 0.0' 0.0", Lat -90.0° 0.0' 0.0"

La distance P est P = 5.80813E+01

A mon avis j'ai soit pas compris un truc soit mal fais un autre.

verdoux

And I'm still waiting

http://fr.wikipedia.org/wiki/Radian

rufo

Pas me confondre avec Lycos!

http://www.developpez.net/forums/d [...] es-resolu/
http://koti.mbnet.fi/ojalesa/googlepages/circle.htm
http://forum.webrankinfo.com/calcu [...] 74024.html
Ton pb s'apparente à un point qui se trouverait ou pas à l'intérieur d'un cercle connaissant les coordonnées de son centre et son rayon.

---------------
Astres, outil de help-desk GPL : http://sourceforge.net/projects/astres, ICARE, gestion de conf : http://sourceforge.net/projects/icare, Outil Planeta Calandreta : https://framalibre.org/content/planeta-calandreta

Profil supprimé

Pourquoi, l'orthodromie c'est pas bien ?

Totoche17

Toi t'as jamais utilisé de fonctions trigo de ta vie dans un programme informatique, je me trompes ?

Profil supprimé

Presque, c'est la première fois.

Totoche17

Les paramètres de toutes les fonctions trigo doivent être exprimées en rad.

Message cité 1 fois

Profil supprimé

L'orthodromie retourne des mile nautique. :jap:

Profil supprimé

En réalité, j'essaie de raccorder deux fonction lié par leurs exploitation.

D'un côté je fais du recouvrement de surface dont voici le code, et de l'autre je regarde avec l'orthodromie si une position se trouve dans la distance de recouvrement.

Code :

--
      function Solid_List (Side : in Side_Type;
                           Ray : in Float;
                           Distance : in Float) return GPS_Coordonates_Manager.Vector is
 
         Terra : GPS_Coordonates_Manager.Vector;
         Diametre : constant Float := Ray * 2.0;
         Perimetre : constant Float := (Pi * Diametre) / 1.852;
         D : constant Float := (1.0/60.0) * Distance;
 
         Count_Of_Equator : constant Long_Long_integer := Long_Long_integer((Perimetre/Distance));
 
         Y : Decimal_Latitude_Type := -90.0;
         X : Decimal_Longitude_Type := -180.0;
      begin
         Text_Io.Put_Line("D =" & Float'Image(D));
         Text_Io.Put_Line("count_of_equator =" & Long_Long_integer'Image(Count_Of_equator));
         for I in 1..Count_Of_Equator/2 loop
            for J in I..Count_Of_Equator loop
 
               declare
                  Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
               begin
                  if Side_Of(GPS) /= Side then
                     Terra := Terra & GPS;
                     Print(Gps);
 
                  end if;
               exception
                  when Constraint_Error =>
                     null;
               end;
                Text_Io.Put_Line("X := " & Float'Image(X));
               exit when X + D > 180.0;
               X := X + D;
 
            end loop;
            Text_Io.Put_Line("Y := " & Float'Image(Y));
 
            X := -180.0;
            exit when Y + D > 90.0;
            Y := Y + D;
 
         end loop;
 
         return Terra;
      end Solid_List;

Mon problème est qu'avec une même distance je ne trouve pas les même points.

Message édité par Profil supprimé le 30-09-2013 à 19:53:05

Publicité

Profil supprimé

Oups, j'ai dit une bêtise.

Je donne distance à solid_list et (distance * (1.0/60.0)) à Is_in.

rufo

Pas me confondre avec Lycos!

Et comme indiquer, faut travailler en radians Or, je vois des 180 et 90 qui me laissent à penser que tu bosses en degrés

Message cité 1 fois

Profil supprimé

Non non, c'est pas le problème.

L'orthodromie exactement ça prend des coordonnées GPS. et ça retourne des mile nautique.|

Message édité par Profil supprimé le 01-10-2013 à 11:21:16

Totoche17

rufo a écrit :

Et comme indiquer, faut travailler en radians Or, je vois des 180 et 90 qui me laissent à penser que tu bosses en degrés

inutile d'insister

Profil supprimé

Surtout que c'est exactement la même chose. des degrés ou radians, c'est que l'unité qui change.

Bon après j'ai une explication à l'effet de divergence des proximité, peut-être que mon orthodromie est bonne, mais l'effet sur un plan 2D en tout cas, ça fait bizarre.
Faut voir que sur un plan 2D les bord sont en réalité des proximité.
Puis ça doit être tout défomé.

Enfin, ça reste obscure.

Alors, si vous savez ce que j'ai fait. Et que vous me comprenez, parce que c'est pas évident. Vous pouvez peut-être me venir en aide.
Je vous en remercie d'avance, merci pour les conseil.

Message cité 1 fois

rufo

Pas me confondre avec Lycos!

Question bête : tu utilises quel système de projection 2D ? Lambert II, WGS84, autre ?

Profil supprimé

J'utilise, la représentation jovasile.

Voici ma fonction To_Map coordonates qui applique une converstion GPS dms en coordonnées 2D

Code :

function To_Map_Coordonates(GPS_Coordonates : in Gps_Coordonates_Type) return Map_Coordonates_Type is
 
      Map_Coordonates : Map_Coordonates_Type;
   begin
      if Gps_Coordonates.Longitude.Degr_Lon >= 0.0 then
 
         Map_Coordonates.X :=
           Max_X_Coordonate/2.0 +
           ((GPS_Coordonates.Longitude.Degr_Lon) * (3600.0 * Lon_unit)) +
           ((GPS_Coordonates.Longitude.Minu_Lon) * (60.0 * Lon_Unit)) +
           ((GPS_Coordonates.Longitude.Seco_Lon) * Lon_Unit);
      else
 
         Map_Coordonates.X :=
           Max_X_Coordonate/2.0 -
           ((GPS_Coordonates.Longitude.Degr_Lon) * (3600.0 * Lon_unit)  * (-1.0)) -
           ((GPS_Coordonates.Longitude.Minu_Lon) * (60.0 * Lon_Unit)) -
           ((GPS_Coordonates.Longitude.Seco_Lon) * Lon_Unit);
     end if;
 
     if Gps_Coordonates.Latitude.Degr_Lat >= 0.0 then
 
         Map_Coordonates.Y :=
           Max_Y_Coordonate/2.0 -
           (((Gps_Coordonates.Latitude.Degr_Lat) * (3600.0 * Lat_unit)) +
           ((Gps_Coordonates.Latitude.Minu_Lat) * (60.0 * Lat_Unit)) +
           ((Gps_Coordonates.Latitude.Seco_Lat) * Lat_Unit));
      else
         Map_Coordonates.Y :=
           Max_Y_Coordonate/2.0 +
           ((Gps_Coordonates.Latitude.Degr_Lat) * (3600.0 * Lat_unit) * (-1.0)) +
           ((Gps_Coordonates.Latitude.Minu_Lat) * (60.0 * Lat_Unit)) +
           ((Gps_Coordonates.Latitude.Seco_Lat) * Lat_Unit);
 
      end if;
 
 
      return Map_Coordonates;
   end To_Map_Coordonates;

Je ne lui reproche que son manque d'arithmétique. En effet, j'aimerais bien assumer les négatif pour les minutes et les seconde également.

Profil supprimé

Je pourrait mieux faire en représentant en 3D ma map sur un plan 2D.
C'est à dire qu'au delà 90° ou -90° longitudinal, je retourne au centre.

Totoche17

Tout à fait, c'est juste une mesure d'angle

Profil supprimé

Si non, je cherche à améliorer ma procédure solid_list.

Dans l'état actuel, je me trouve avec un horizon de 360 dégrés aux pôles.

J'ai entrepris de faire une troisième boucle pour itérer du nord au sud l'index longitudinal, mais je me suis embrouillé plus qu'autre chose.

Des idées ?

Profil supprimé

Bon, j'ai trouvé quoi faire, :

Code :

function Solid_List (Side : in Side_Type;
               Ray : in Float;
               Distance : in Float) return GPS_Coordonates_Manager.Vector is
     
     Terra : GPS_Coordonates_Manager.Vector;
     Diametre : constant Float := Ray * 2.0;   
     Perimetre : constant Float := ((Pi * Diametre) / 1.852);
     D : constant Float := (1.0/60.0) * Distance;
     
     Count_Of_Equator : constant Long_Long_integer := Long_Long_integer((Perimetre/Distance));
          
     Y : Decimal_Latitude_Type := -90.0 + D/2.0;
     X : Decimal_Longitude_Type := 0.0+ D/ 2.0;            
      begin 
     Text_Io.Put_Line("D =" & Float'Image(D));
     Text_Io.Put_Line("count_of_equator =" & Long_Long_integer'Image(Count_Of_equator));                                 
     
     case Side is     
        when Ouest =>
       North_est :
           for I in 1..Count_Of_Equator/2 loop
          for J in 1..I loop
             for K in 1..Count_Of_Equator loop
            for L in  J..K loop                
               
               declare
                  Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
               begin                      
                  if Side_Of(GPS) /= Side then
                 Terra := Terra & GPS;
                 Print(Gps);
                 
                  end if;
               exception
                  when Constraint_Error =>
                 null;             
               end;        
               --Text_Io.Put_Line("X := " & Float'Image(X));                   
               exit when X + D > 180.0;
               X := X + D;
               
            end loop;          
            --Text_Io.Put_Line("Y := " & Float'Image(Y));
            
            X := 0.0 + D/2.0;            
            exit North_est when Y + D > 0.0;
            Y := Y + D;                        
             end loop;
          end loop;      
           end loop North_est;
           
           Y := 90.0;
       South_Est :
           for I in 1..Count_Of_Equator/2 loop
          for J in 1..I loop
             for K in 1..Count_Of_Equator loop
            for L in  J..K loop                
               
               declare
                  Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
               begin                      
                  if Side_Of(GPS) /= Side then
                 Terra := Terra & GPS;
                 Print(Gps);
                 
                  end if;
               exception
                  when Constraint_Error =>
                 null;             
               end;        
               --Text_Io.Put_Line("X := " & Float'Image(X));                   
               exit when X + D > 180.0;
               X := X + D;
               
            end loop;          
            --Text_Io.Put_Line("Y := " & Float'Image(Y));
            
            X := 0.0 + D/2.0;            
            exit South_est when Y - D < 0.0;
            Y := Y - D;                        
             end loop;
          end loop;      
           end loop South_Est;
           
           
           
           
           
        when Est =>                
           X := -180.0 + D/ 2.0;            
       North_West :
           for I in 1..Count_Of_Equator/2 loop
          for J in 1..I loop
             for K in 1..Count_Of_Equator loop
            for L in  J..K loop                
               
               declare
                  Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
               begin                      
                  if Side_Of(GPS) /= Side then
                 Terra := Terra & GPS;
                 Print(Gps);
                 
                  end if;
               exception
                  when Constraint_Error =>
                 null;             
               end;        
               --Text_Io.Put_Line("X := " & Float'Image(X));                   
               exit when X + D > 180.0;
               X := X + D;
               
            end loop;          
            --Text_Io.Put_Line("Y := " & Float'Image(Y));
            
            X := -180.0 + D/2.0;            
            exit North_west when Y + D > 0.0;
            Y := Y + D;                        
             end loop;
          end loop;      
           end loop North_west;           
           
           Y := 90.0;
       South_west:
           for I in 1..Count_Of_Equator/2 loop
          for J in 1..I loop
             for K in 1..Count_Of_Equator loop
            for L in  J..K loop                
               
               declare
                  Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
               begin                      
                  if Side_Of(GPS) /= Side then
                 Terra := Terra & GPS;
                 Print(Gps);
                 
                  end if;
               exception
                  when Constraint_Error =>
                 null;             
               end;        
               --Text_Io.Put_Line("X := " & Float'Image(X));                   
               exit when X + D > 180.0;
               X := X + D;
               
            end loop;          
            --Text_Io.Put_Line("Y := " & Float'Image(Y));
            
            X := -180.0 + D/2.0;            
            exit South_west when Y - D < 0.0;
            Y := Y - D;                        
             end loop;
          end loop;      
           end loop South_west;           
 
           
     end case;
           
           
     return Terra;     
   end Solid_List;

Ce qui fait 8281 surface de 1° de côté. Les coordonnées obtenue son les centres de ces surface.

Profil supprimé

La revoici corrigée.

Code :

--
    function Solid_List (Side : in Side_Type;
            Ray : in Float;
            Distance : in Float) return GPS_Coordonates_Manager.Vector is
      
      Terra : GPS_Coordonates_Manager.Vector;
      Diametre : constant Float := Ray * 2.0;   
      Perimetre : constant Float := ((Pi * Diametre) / 1.852);
      D : constant Float := (1.0/60.0) * Distance * 2.0;
      
      Count_Of_Equator : constant Long_Long_integer := Long_Long_integer((Perimetre/Distance));
      
      Y : Decimal_Latitude_Type := -90.0 + D /2.0;
      X : Decimal_Longitude_Type := 0.0+ D / 2.0;            
   begin 
      Text_Io.Put_Line("D =" & Float'Image(D));
      Text_Io.Put_Line("count_of_equator =" & Long_Long_integer'Image(Count_Of_equator));                                 
      
      case Side is     
     when Ouest =>
    North_est :
        for I in 1..Count_Of_Equator/2 loop
           for J in 1..Count_Of_Equator/2 loop
          for K in 1..I loop
             for L in  1..J*2 loop                
            
            declare
               Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
            begin                      
               if Side_Of(GPS) /= Side then
                  Terra := Terra & GPS;
                  Print(Gps);
                  
               end if;
            exception
               when Constraint_Error =>
                  null;             
            end;        
            --Text_Io.Put_Line("X := " & Float'Image(X));                   
            exit North_Est when X + D > 180.0;
            X := X + D;
            
             end loop;          
 
             --Text_Io.Put_Line("Y := " & Float'Image(Y));
             X := 0.0 + D/2.0;            
             exit when Y + D > 0.0;
             Y := Y + D;                        
 
             
          end loop;
           end loop;      
        end loop North_est;
        X := 0.0 + D/2.0;            
        Y := 90.0 - D/2.0;
    South_Est :
        for I in 1..Count_Of_Equator/2 loop
           for J in 1..Count_Of_Equator/2 loop
          for K in 1..I loop
             for L in  1..J*2 loop                
            
            declare
               Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
            begin                      
               if Side_Of(GPS) /= Side then
                  Terra := Terra & GPS;
                  Print(Gps);
                  
               end if;
            exception
               when Constraint_Error =>
                  null;             
            end;        
            --Text_Io.Put_Line("X := " & Float'Image(X));                   
            exit South_Est when X + D > 180.0;
            X := X + D;
            
             end loop;          
             --Text_Io.Put_Line("Y := " & Float'Image(Y));
          X := 0.0 + D/2.0;            
          exit when Y - D < 0.0;
          Y := Y - D;                        
          
 
             
          end loop;
                    
           end loop;      
        end loop South_Est;
        
        
        
        
        
     when Est =>                
        X := 0.0 - D/ 2.0;            
    North_West :
        for I in 1..Count_Of_Equator/2 loop
           for J in 1..Count_Of_Equator/2 loop
          for K in 1..I loop
             for L in  1..J*2 loop                
            
            declare
               Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
            begin                      
               if Side_Of(GPS) /= Side then
                  Terra := Terra & GPS;
                  Print(Gps);
                  
               end if;
            exception
               when Constraint_Error =>
                  null;             
            end;        
            --Text_Io.Put_Line("X := " & Float'Image(X));                   
            exit when X - D < -180.0;
            X := X - D;
            
             end loop;          
             
             --Text_Io.Put_Line("Y := " & Float'Image(Y));
             
             X := 0.0 - D/2.0;            
             exit North_west when Y + D > 0.0;
             Y := Y + D;        
             
          end loop;
          
           end loop;      
        end loop North_west;           
        X := -0.0 + D/ 2.0;            
        Y := 90.0 - D/2.0;
    South_west:
        for I in 1..Count_Of_Equator/2 loop
           for J in 1..Count_Of_Equator/2 loop
          for K in 1..I loop
             for L in  1..J*2 loop                
            
            
            declare
               Gps : constant Gps_Coordonates_Type := To_Gps(X, Y);
            begin                      
               if Side_Of(GPS) /= Side then
                  Terra := Terra & GPS;
                  Print(Gps);
                  
               end if;
            exception
               when Constraint_Error =>
                  null;             
            end;        
            --Text_Io.Put_Line("X := " & Float'Image(X));                   
            exit South_West when X - D < -180.0;
            X := X - D;
            
             end loop;          
             --Text_Io.Put_Line("Y := " & Float'Image(Y));
             X := 0.0 - D/2.0;            
             exit when Y - D < 0.0;
             Y := Y - D;                        
 
          end loop;
 
           end loop;      
        end loop South_west;           
 
        
      end case;
      
      
      return Terra;     
   end Solid_List;

Publicité

FORUM HardWare.fr

Programmation

Algo

Calcul d'un orthodromie et recouvrement d'une sphère.

Sujets relatifs
Calcul de d'une orthodromie avec Ada	ecrire un programme (pascal) qui calcule X puissance n ?!!!
[opencv] calcule du jacobien	fonction strlen calcule de la longueur de chaine
aidez moi : que calcule cette procédure	calcule de la longeur de la variable
[Access 2003] champs calculé à partir de 2 tables différentes	[SQL server] Recherche sur champ calculé
Calcule en base 13??	Classement suivant un total calculé
Plus de sujets relatifs à : Calcul d'un orthodromie et recouvrement d'une sphère.

Page générée en 0.189 secondes