Le portail francophone de la géomatique

Classique des classiques: la question qui tue le soir du départ en WE, pas l'irremplaçable Miss T.

Pour répondre vite "desfois" que Miss T. ne veuille partir en week end, une réponse provenant de
http://williams.best.vwh.net/avform.htm#Crs
"Course" veut dire Cap en anglais, et cap c'est l'angle / au nord que fait la route du point 1 au point 2

Course between points

We obtain the initial course, tc1, (at point 1) from point 1 to point 2 by the following. The formula fails if the initial point is a pole. We can special case this with:

IF (cos(lat1) < EPS)   // EPS a small number ~ machine precision
  IF (lat1 > 0)
     tc1= pi        //  starting from N pole
  ELSE
     tc1= 2*pi         //  starting from S pole
  ENDIF
ENDIF

For starting points other than the poles:

IF sin(lon2-lon1)<0       
   tc1=acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))   
ELSE       
   tc1=2*pi-acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))   
ENDIF

An alternative formula, not requiring the pre-computation of d, the distance between the points, is:

   tc1=mod(atan2(sin(lon1-lon2)*cos(lat2),
           cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon1-lon2)), 2*pi)

et la distance d se calcule:
Distance between points

The great circle distance d between two points with coordinates {lat1,lon1} and {lat2,lon2} is given by:

d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))

A mathematically equivalent formula, which is less subject to rounding error for short distances is:

d=2*asin(sqrt((sin((lat1-lat2)/2))^2 +
                 cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))

Dernière modification par xavier29 (Thu 07 May 2009 18:07)

Merci les gars pour ces réponses.
Maurice, le champ MN, c'est la direction? Si oui, crois-tu qu'une ligne orientée NNE-SSW puisse avoir une direction de 70°?!?

Si j'ai bien compris (je résume car je suis une authentique biquette en trigo...):
1/ ajout des colonnes sin et cos, les remplir
2/ faire autant de requêtes que de "quartiers" trigonométriques pour trouver le quartier dans lequel on est:
Sin x > 0 et cos x > 0 ==> on est dans quart NE
Sin x < 0 et cos x > 0 ==> SE
Sin x < 0 et cos x < 0 ==> SW
Sin x > 0 et cos x < 0 ==> NW

Mais euh... (désolée pour la question, ça va être hyper... blond) comment on transpose l'angle du radian en degré en passant d'un cercle trigo à l'autre? (je vous avais prévenu...)
Sin x > 0 et cos x > 0 ==> NE => acos(cos x) * rad_2_deg
Sin x < 0 et cos x > 0 ==> SE => ?
Sin x < 0 et cos x < 0 ==> SW => ?
Sin x > 0 et cos x < 0 ==> NW => ?

Dernière modification par Miss T. (Mon 11 May 2009 18:03)

L'outil de base est bien le cercle trigo...
J'ai effectivement calculé sinus et cosinus pour savoir dans quel quartier on se trouve
J'ai aussi créé les "points de départ" pour orienter les polylignes (fichier joint)
Les valeurs des angles sont alors suivant le quartier, en degrés exprimés dans le sens trigo
Sin x > 0 et cos x > 0 ==> NE => angles de 0 à 90°
Sin x < 0 et cos x > 0 ==> SE => angles de 270 à 360°
Sin x < 0 et cos x < 0 ==> SW => angles de 180 à 270°
Sin x > 0 et cos x < 0 ==> NW => angles de 90 à 180°
J'ai ensuite utilisé, dans le champ MN, atn(tan x) = atn((Ymax-Ymin)/(Xmax-Xmin)) et, comme le résultat retourné par Atn( ) est, en radians, entre -Pi/2 radians et Pi/2 radians, j'ai ajouté ce qu'il faut pour tomber sur les valeurs ci-dessus dans chacune des 4 sélections...
Et une ligne qui part du NE pour aboutir dans le SW doit avoir un angle entre 180 et 270°...si on ne s'est pas trompé sur point de départ et point d'arrivée...ce qui a pu m'arriver

Pas facile d'être explicite...
Aller en direction de = de l'origine vers l'extrémité
Une ligne qui va en direction du NE => Sin x > 0 et cos x > 0 => angles de 0 à 90°
Une ligne qui va en direction du NW => Sin x > 0 et cos x < 0 => angles de 90 à 180°
Une ligne qui va en direction du SW => Sin x < 0 et cos x < 0 => angles de 180 à 270° (ou de -90 à -180°)
Une ligne qui va en direction du SE => Sin x < 0 et cos x > 0 => angles de 270 à 360° (ou de 0 à -90°)

Encore une soirée télé de fichue

Pour commencer je ne vois pas bien d'où sort ce 260 dans la "remise au cercle trigonométrique"...270 j'aurais mieux compris
Ensuite la dénomination des directions de la "rose des vents" est curieuse...Pour mémoire, voici une belle rose des vents: [img]http://kramvoussanos.unblog.fr/files/2009/02/rose20des20vents.png[/img]

Le film était bien...
J'ai juste écrit ça, avec ArcTangente (désolé !), en faisant au début une sélection pour créer la table provisoire "tata" dont l'objet est juste de rendre le code un peu plus lisible

include "mapbasic.def"
dim table_name as string
dim deltaY, deltaX as float 

table_name = TableInfo(1, 1)   
Select Nb_Noeuds, Sinuosite, GeoAngle, Axe, (ObjectNodeY(obj,1,ObjectInfo(obj,20)) - ObjectNodeY(obj,1,1)) "deltaY", (ObjectNodeX(obj,1,ObjectInfo(obj,20)) - ObjectNodeX(obj,1,1)) "deltaX" from Pour_GeoRezo into tata
 
    'remise au cercle trigonométrique        
 Select * from tata where deltaY > 0 and deltaX > 0  Into Orientation  '==> NE => angles de 0 à 90°
  Update Orientation Set GeoAngle = atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg
 Select * from tata where deltaY > 0 and deltaX < 0  Into Orientation  '==> NW => angles de 90 à 180°
  Update Orientation Set GeoAngle = atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg+90
 Select * from tata where deltaY < 0 and deltaX < 0  Into Orientation  '==> SW => angles de 180 à 270°
  Update Orientation Set GeoAngle = atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg+180
 Select * from tata where deltaY < 0 and deltaX > 0  Into Orientation  '==> SE => angles de 270 à 360°
  Update Orientation Set GeoAngle = atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg
 Close Table tata
 Commit Table table_name

    'direction dans rose des vents
 Select * from  table_name where GeoAngle between 0 AND 22.5 Into Orientation
  Update Orientation set Axe = "E"  
 Select * from  table_name where GeoAngle between 22.5 AND 67.5 Into Orientation
  Update Orientation set Axe = "NE"          
 Select * from  table_name where GeoAngle between 67.5 AND 112.5 Into Orientation
  Update Orientation set Axe = "N" 
 Select * from  table_name where GeoAngle between 112.5 AND 157.5 Into Orientation
  Update Orientation set Axe = "NW"
 Select * from  table_name where GeoAngle between 157.5 AND 202.5 Into Orientation
  Update Orientation set Axe = "W"
 Select * from  table_name where GeoAngle between 202.5 AND 247.5 Into Orientation
  Update Orientation set Axe = "SW"         
 Select * from  table_name where GeoAngle between 247.5 AND 292.5 Into Orientation
  Update Orientation set Axe = "S"    
 Select * from  table_name where GeoAngle between 292.5  AND 337.5 Into Orientation
  Update Orientation set Axe = "SE"    
 Select * from  table_name where GeoAngle between 337.5  AND 360 Into Orientation
  Update Orientation set Axe = "E" 
 Close Table Orientation
 Commit Table table_name
 Run Menu Command 304

Reste à vérifier que ce code (intitulé MissTrigo.mb ) donne satisfaction..

Effectivement non !!
Les directions sont "orientées", soit un 0° pointant vers l'Est (l'orient)
Le code est donc écrit pour un Nord à 90°

Désolé mais non,

L'axe des X pointe toujours vers l'est, ensuite, pour l'écriture il faut parcourir la rose des vents dans le sens trigonométrique
donc le sens opposé au sens des aiguilles d'une montre... et non pas dans le sens des aiguilles d'une montre comme tu l'as écrit.

Ce n'est pas grave d'utiliser arctangente à la place d'arcsin, sauf pour les cas ou cos X = 0...
d'ou la modif pour éviter le 0 dans les calculs...
J'ai utilisé le sinus pour rester dans le même type de formule pour la détermination de l'angle réel...
pour la tangente ,il faut faire attention...

J'ai fait une petite erreur hier dans mon code ... Mon résultat variais de -90 à 270 et non pas de 0 à 360

poursuivons...
de 90 à 180° le sinus est positif et le cosinus est négatif. Donc arctangente varie de -90 à 0 donc + 180
Si + 90 alors géoangle dans ce cas variera entre 0 et 90 désolé maurice...

déterminer l'angle avec le cosinus est plus facile...compris entre -1 et 1 le cosinus varie de 0 à 180°
si le sinus est négatif l'angle est négatif donc -180 à 0... reste plus qu'a lui ajouter 360.. pour une variation de 180 à 360

Select * from tata where deltaY > 0 and deltaX < 0  Into Orientation  '==> NW => angles de 90 à 180°
  Update Orientation Set GeoAngle = atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg+180
'tangente (X) <0  attention -90 <= arctangente (X) < 0'
Select * from tata where deltaY < 0 and deltaX < 0  Into Orientation  '==> SW => angles de 180 à 270°
  Update Orientation Set GeoAngle = atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg+180
' tangente (X)>=0 donc 0 <= arctangente (X) <= 90
Select * from tata where deltaY < 0 and deltaX > 0  Into Orientation  '==> SE => angles de 270 à 360°
  Update Orientation Set GeoAngle = atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg+360
'dans ce cas tangente négative. comme on veut une valeur positive, on rajoute un tour de roue

Dernière modification par Damien BEAUSEIGNEUR (Wed 13 May 2009 12:04)

méa culpa...
il faut enlever (abs()) dans la formule....
tel que c'est écrit la arctangente donnera toujours une valeur comprise entre 0 et pi/2
ou changer le calcul de l'angle.
compte tenu du réglage le + 0.001 pour éviter le 0 je garde la valeur absolue...

Select * from tata where deltaY > 0 and deltaX < 0  Into Orientation  '==> NW => angles de 90 à 180°
  Update Orientation Set GeoAngle = 180 - atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg
'cela revient à retirer l'angle de 180 puisque l'on retourne en arrière
Select * from tata where deltaY < 0 and deltaX < 0  Into Orientation  '==> SW => angles de 180 à 270°
  Update Orientation Set GeoAngle = 180 + atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg
' cas inverse
Select * from tata where deltaY < 0 and deltaX > 0  Into Orientation  '==> SE => angles de 270 à 360°
  Update Orientation Set GeoAngle = 360 - atn(abs(0.001+deltaY)/abs(0.001+deltaX))* rad_2_deg
' pour retirer l'angle au tour complet...

On y arrive, les formules deviennent compliquée... je trouve tout ça pour éviter le cos(X) = 0... et la distance nulle...
enfin là ça devrait marcher...