Convertion Longitude/Latitude en Pixel [résolu]

[MetalOS]

Ok, merci kernadec pour ton aide.

[MetalOS]

Ta fait comment kernadec pour générer ta carte google maps ?

[kernadec]

Bonjour MetalOs
j'utilise ce programme très pratique : http://www.clubic.com/telecharger-fiche ... wsnap.html

je viens de mettre a jours le code..
car j'avais un ptite erreur dans la procédure de rotation un and un lieu d'un or

Cordialement

[MetalOS]

Je vient de trouver ca pour exporter des carte. C'est Top OpenStreet Map.

http://www.openstreetmap.org/?lat=46.97 ... rs=B000FTF

[Scrat]

Salut

Si tu utilises une projection utm/mercator (gmaps/osm etc..) voici une procedure qui te retourne les coordonnées X/Y d'une Lat/lon

en entrée
la latitude/longitude en degré decimaux
Le zoom un entier

en sortie
Tu recuperes dans la structure TileXY le XY depuis le zero mercator et depuis le coin sup-gauche de la tile contenant le point

Code : Tout sélectionner

;#####################################################

;Recup du X/Y d'une Lat/lon par rapport au 0 et par rapprot an coin NW d'une tile

;#####################################################

Procedure.b GetXYOnTile(LatitudeCalc.d,LongitudeCalc.d,ZoomCarte.b)

  WTiles.d = 256*(Pow(2, ZoomCarte))

  CorrectedLong.d = 180 + LongitudeCalc

  LongTileSize.d = 360/WTiles

  Xori = Int(CorrectedLong/LongTileSize)

  XPos = Xori % 256

  

  CorrectedLat.d = 0

  If (LatitudeCalc < -90) Or (LatitudeCalc > 90)

    CorrectedLat = 90 - LatitudeCalc

  Else

    CorrectedLat = LatitudeCalc

  EndIf

  

  phi.d = #PI * CorrectedLat / 180

  Mercator.d = 0.5 * Log((1+Sin(phi))/(1-Sin(phi)))

  Yori=Int(((1 - (Mercator / #PI)) / 2) * WTiles)

  Ypos = Yori % 256      

  

  TileXY\Xori=Xori

  TileXY\Yori=Yori

  TileXY\Xpos=Xpos

  TileXY\Ypos=Ypos 

  

EndProcedure

A+

[kernadec]

bonjour MetalOS et Scart

je me suis permis de d'adapter en Purebasic le code TCL/TK pour la conversion des coordonnées UTM et Vice Versa
Après plusieurs comparaisons avec TCL/TK, ca semble fonctionner normalement.

correction de la fonction degrés/décimales qui ne prenait pas en compte les décimales pour les secondes.
pour utiliser ses coordonnées sur une carte il faut définir un lieu UTM et projeter les UTM sur la carte à partir de ce lieu.

Cordialement

Code : Tout sélectionner

;##############################################################################
;## adaptation du code tcl.tk de conversion Mercator par kernadec en PB441   ##
;## version tcl.tk : http://wiki.tcl.tk/11109                                ##
;## http://www.dmap.co.uk/utmworld.htm                                       ##
;## http://en.wikipedia.org/wiki/File:Utm-zones.jpg                          ##
;##############################################################################

Global latitude.d,longitude.d,easting.d,northing.d,zone.d,letter.s,xm.d,ym.d 
Macro modulo(xm,ym)                              ; enleve entier garde le reste 
  xm-ym*Int(xm/ym)
EndMacro
Procedure.d Frac(Xfrac.d)
  ProcedureReturn Xfrac-Int(Xfrac)
EndProcedure
Procedure.d dms(Xdms.d)
  ;conversion decimales en degres minutes secondes 
  ProcedureReturn Int(Xdms)+Frac(Int(Frac(Xdms)*60)/100)+Frac(Frac(Frac(Xdms)*60)*0.006)
EndProcedure
Procedure.d dec(xdec.d)
  ;conversion  degres minutes secondes  en decimales
  ProcedureReturn Int(Xdec)+(Int(frac(Xdec)*100)/60)+((((frac(Xdec)*100)-(Int(frac(Xdec)*100)))*100)/3600)
EndProcedure
;##+##########################################################################
;#
; # ll2utm -- Converts latitude And longitude into Universal Transverse
; # Mercator (UTM) coordinates. Lots of fun math which I got off the Web
;#
Procedure ll2utm(latitude.d,longitude.d)
  Protected K0.d,es2.d,es4.d,es6.d,es2x.d,er.l,l.d,N.d,T.d,C.d,A.d,M.d
  Protected long_origin.d,long_origin_rad.d,eccPrimeSquared.d,lat_rad.d,long_rad.d
  K0=0.9996
  ;# WGS-84
  er=6378137                              ;# EquatorialRadius
  es2=0.00669438                          ;# EccentricitySquared
  es4=Pow(0.00669438,2)                   
  es6=Pow(0.00669438,3)                   
  While longitude < 0-180
    longitude=longitude+360
  Wend 
  While longitude >= 180
    longitude=longitude-360
  Wend     
  ;# Now convert
  latitude=dec(latitude)
  longitude=dec(longitude)
  Debug "> Lat/long en decimales"
  Debug latitude
  Debug longitude
  lat_rad=latitude * #PI/180.0
  long_rad=longitude * #PI/180.0
  zone=Int((longitude + 180) / 6) + 1 
  If latitude >= 56.0 And latitude < 64.0 And longitude >= 3.0 And longitude < 12.0
    zone = 32
  EndIf
  If latitude >= 72.0 And latitude < 84.0
    If longitude >= 0.0  And  longitude <  9.0:zone = 31:EndIf
    If longitude >= 9.0  And  longitude < 21.0:zone = 33:EndIf
    If longitude >= 21.0 And  longitude < 33.0:zone = 35:EndIf
    If longitude >= 33.0 And  longitude < 42.0:zone = 37:EndIf
  EndIf
  ;# +3 puts origin in middle of zone
  long_origin=(zone-1)*6-180+3
  long_origin_rad=long_origin*#PI/180.0
  eccPrimeSquared=es2/(1.0-es2)
  N=er/Sqr(1.0-es2*Sin(lat_rad)*Sin(lat_rad))
  T=Tan(lat_rad)*Tan(lat_rad)
  C=eccPrimeSquared*Cos(lat_rad)*Cos(lat_rad)
  A=Cos(lat_rad)*(long_rad-long_origin_rad)
  M=er*((1.0-es2/4-3*es4/64-5*es6/256)*lat_rad-(3*es2/8+3*es4/32+45*es6/1024)*Sin(2*lat_rad)+(15*es4/256+45*es6/1024)*Sin(4*lat_rad)-(35*es6/3072)*Sin(6*lat_rad))
  easting=K0*N*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)+500000.0
  northing=K0*(M+N*Tan(lat_rad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24+(61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720))
  If latitude < 0      ;# 1e7 meter offset for southern hemisphere
    northing=northing+10000000.0
  EndIf
  ; northing=Int(northing)
  ; easting=Int(easting)
  If latitude > 84.0 Or latitude < -80.0
    letter="Z"
  Else 
    l=Int((latitude + 80) / 8.0)
    letter=Mid("CDEFGHJKLMNPQRSTUVWXX",l+1,1) 
  EndIf
  
EndProcedure

; ##+##########################################################################
; #
; # utm2ll -- Converts Universal Transverse Mercator (UTM) into
; # latitude And longitude coordinates. More fun math which I got off the Web
; #
Procedure  utm2ll(northing.d,easting.d,zone.d,letter.s)
  Protected K0.d,es2.d,es4.d,es6.d,es2x.d,er.l,l.d
  Protected x.d,y.d,northernHemisphere.d,long_origin.d
  Protected ep2.d,e1.d,M.d,mu.d,phi.d,N1.d,T1.d,C1.d,R1.d,D.d 
  K0=0.9996
  ;# WGS-84
  er=6378137                              ;# EquatorialRadius
  es2=0.00669438                          ;# EccentricitySquared
  es4=Pow(0.00669438,2)                   
  es6=Pow(0.00669438,3)                  
  es2x=1.0-es2
  x=easting-500000.0
  If Asc(letter)>=78
    northernHemisphere=1
    y=northing
  Else
    northernHemisphere=0
    y=northing-10000000.0
  EndIf
  long_origin=(zone-1)*6-180+3                  ;# +3 puts in middle
  ep2=es2/es2x
  e1=(1.0-Sqr(es2x))/(1.0+Sqr(es2x))
  M=y/K0
  mu=M/(er*(1.0-es2/4.0-3*es2*es2/64.0-5*es2*es2*es2/256.0))
  phi=mu+(3*e1/2-27*e1*e1*e1/32)*Sin(2*mu)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*Sin(4*mu)+(151*e1*e1*e1/96)*Sin(6*mu)
  N1=er/Sqr(1.0-es2*Sin(phi)*Sin(phi))
  T1=Tan(phi)*Tan(phi)
  C1=ep2*Cos(phi)*Cos(phi)
  R1=er*es2x/Pow(1.0-es2*Sin(phi)*Sin(phi),1.5)
  D=x/(N1*K0)
  latitude=phi-(N1*Tan(phi)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*ep2)*D*D*D*D/24+(61+90*T1+298*C1+45*T1*T1-252*ep2-3*C1*C1)*D*D*D*D*D*D/720)
  latitude=latitude*180.0/#PI
  longitude=(D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*ep2+24*T1*T1)*D*D*D*D*D/120)/Cos(phi)
  longitude=long_origin+longitude*180.0/#PI
  
EndProcedure
Debug "UTM"
;43°38'33.24?N 79°23'13.7?W? / ?43.6425667°N 79.387139°W? / 43.6425667; -79.387139? (CN Tower).  This is in zone 17, and the grid position is 630084m east, 4833438m north. 

ll2utm(43.383324,-79.23137)
  Debug "> Coordonnées décimales"
  Debug northing
  Debug easting
  Debug "> Coordonnées entieres"
  Debug Int(northing)
  Debug Int(easting)
  Debug zone
  Debug letter
utm2ll(northing.d,easting.d,zone.d,letter.s)
  Debug "> Sorties Lat/long en decimales"
  Debug latitude
  Debug longitude
  Debug "> Lat/long en degres"
  Debug dms(latitude)
  Debug dms(longitude)
  Debug "> comparaison avec les données de départ"
  Debug "43°38'33.24?N 79°23'13.7?W"
; #  1. UTM Nord, fuseau 30 : entre 6 degrés ouest et 0 degré Greenwich ;
; #  2. UTM Nord, fuseau 31 : entre 0 degré et 6 degrés est Greenwich ;
; #  3. UTM Nord, fuseau 32 : entre 6 degrés est et 12 degrés est Greenwich.


; ################################################################
; ## Other datums which could be used here instead:
; ##
; ## Name                         EquatorialRadius        EccentricitySquared
; ## ----                         ----------------        -------------------
; ## Airy                         6377563                 0.00667054
; ## Australian National          6378160                 0.006694542
; ## Bessel 1841                  6377397                 0.006674372
; ## Bessel 1841 (Nambia)         6377484                 0.006674372
; ## Clarke 1866                  6378206                 0.006768658
; ## Clarke 1880                  6378249                 0.006803511
; ## Everest                      6377276                 0.006637847
; ## Fischer 1960 (Mercury)       6378166                 0.006693422
; ## Fischer 1968                 6378150                 0.006693422
; ## GRS 1967                     6378160                 0.006694605
; ## GRS 1980                     6378137                 0.00669438
; ## Helmert 1906                 6378200                 0.006693422
; ## Hough                        6378270                 0.00672267
; ## International                6378388                 0.00672267
; ## Krassovsky                   6378245                 0.006693422
; ## Modified Airy                6377340                 0.00667054
; ## Modified Everest             6377304                 0.006637847
; ## Modified Fischer 1960        6378155                 0.006693422
; ## South American 1969          6378160                 0.006694542
; ## WGS 60                       6378165                 0.006693422
; ## WGS 66                       6378145                 0.006694542
; ## WGS-72                       6378135                 0.006694318
; ## WGS-84                       6378137                 0.00669438

; ################################################################

[MetalOS]

Salut Kernadec, je revient vert toi car tu semble très bien maitrisé le sujet. Pense tu qu'il est possible d'obtenir une distance en Km entre deux points en complément de ton code.

Code : Tout sélectionner

;#############################################################################
;## kernadec Conversion latitude longitude en pixels et vice versa  05/2010 ##
;#############################################################################
UseJPEGImageDecoder()
Global echelle.l,latitudepix.d,longitudepix.d,latitude.d,longitude.d,xdessin.l,ydessin.l,imgx.l,imgy.l,posx.d,posy.d,angle.d,rayon.d
xdessin=221
ydessin=272
echelle.l=3800
img.l=LoadImage(1,"c:\purebasic\Francegooglemap.jpg")
If img<>0
  imgx=ImageWidth(1):imgy=ImageHeight(1) 
Else
  imgx=650:imgy=650
EndIf
Procedure rotation(x.d,y.d,decalage.d) 
  ;correction de l'angle pour l'ensemble des points d'un demi degrés pour la zone France
  posx=xdessin-x:posy=ydessin-y
  rayon=Sqr(Pow(posx,2)+Pow(posy,2))
  If posx<>0 or posy<>0
    If x>xdessin And y>ydessin:angle=Abs(ATan(posx/posy)/#PI*180):EndIf
    If x>xdessin And y<=ydessin:angle=Abs(ATan(posy/posx)/#PI*180)+90:EndIf
    If x<=xdessin And y<=ydessin:angle=Abs(ATan(posx/posy)/#PI*180)+180:EndIf
    If x<=xdessin And y>ydessin:angle=Abs(ATan(posy/posx)/#PI*180)+270:EndIf
  Else:angle=0:EndIf
  posx=xdessin+rayon*Sin((angle+decalage)*#PI/180)
  posy=ydessin+rayon*Cos((angle+decalage)*#PI/180)
EndProcedure   
Procedure convert_lat_long(y.d,x.d)
  ;longitudepix=xdessin+((Sin(x*#PI/180)*Cos(y*#PI/180))*echelle)
  longitudepix=xdessin+((Sin(x*#PI/180)*Cos(47*#PI/180))*echelle)  ; latitude moyenne de 47 pour calcule de la longitude
  latitudepix=ydessin-((Sin((y-47)*#PI/180))*echelle)
  ;Circle(longitudepix,latitudepix,3,RGB(255,255,255))
  rotation(longitudepix,latitudepix,0.5)
  Circle(posx,posy,2,RGB(255,0,0))
EndProcedure
Procedure convert_pixel(my.l,mx.l)
  rotation(mx,my,-0.5)
  latitude=47+(ASin(((ydessin-posy)/echelle))/#PI*180) 
  ;longitude=ASin(((mx-xdessin)/echelle)/Cos(latitude*#PI/180))/#PI*180
  longitude=ASin(((posx-xdessin)/echelle)/Cos(47*#PI/180))/#PI*180  ; latitude moyenne = 47 pour le retour
EndProcedure

OpenWindow(0, 0,0,imgx,imgy, "Fenetre Topo ", #PB_Window_SystemMenu|#PB_Window_ScreenCentered)
StartDrawing(WindowOutput(0)) 
  If img<>0:DrawImage(img,0,0):EndIf
  LineXY(0,ydessin,imgx,ydessin,RGB(255,0,0))  ; latitude de 47
  LineXY(xdessin,0,xdessin,imgy,RGB(255,0,0))  ; méridien 0
  ; ville en degrés décimaux
  convert_lat_long(49.633,-1.617)     ; Cherbourg
  convert_lat_long(48.8139,-3.4429)   ; Peros Guirec
  convert_lat_long(48.40,-4.5166)     ; Brest
  convert_lat_long(47.7500,-3.3636)   ; Lorient
  convert_lat_long(46.167,-1.167)     ; La Rochelle
  convert_lat_long(43.4833,-1.55)     ; Biarritz
  convert_lat_long(42.7000,2.900)     ; Perpignan
  convert_lat_long(43.3000,5.3670)    ; Marseille
  convert_lat_long(41.917,8.717)      ; Ajaccio
  convert_lat_long(42.666,9.50)       ; Bastia
  convert_lat_long(43.700,7.25)       ; Nice   
  convert_lat_long(45.7830,3.0830)    ; Clermont_ferrand
  convert_lat_long(45.8325,6.8643)    ; Mont Blanc   
  convert_lat_long(48.5833,7.75)      ; Strasbourg
  convert_lat_long(48.9655,8.2305)    ; point nord est france allemagne            
  convert_lat_long(49.1196,6.1761)    ; Metz
  convert_lat_long(48.8333,2.3333)    ; Paris   
  convert_lat_long(50.9581,1.8520)    ; Calais
  convert_lat_long(49.4970,0.1276)    ; Havre (le)
  convert_lat_long(47,0)              ; centre xy
  convert_lat_long(51,0)              ; meridien zero à 51
  convert_lat_long(42,0)              ; meridien zero à 42
StopDrawing()  
Repeat: 
  event = WaitWindowEvent()   
  Select event
    Case #PB_Event_CloseWindow
      End
  EndSelect
  y=WindowMouseY(0):x=WindowMouseX(0)
  convert_pixel(y,x)
  StartDrawing(WindowOutput(0))
    DrawText(10,10,StrD(latitude,4),RGB(0,0,0),RGB(255,255,255))
    DrawText(10,30,StrD(longitude,4)+" ",RGB(0,0,0),RGB(255,255,255))
  StopDrawing()
ForEver

Le but serai de prendre une position de départ que l'utilisateur choisirais (Paris par exemple) et de calculer la distance sur la crate entre la position de départ et la position du pointeur de la souris sur la carte. Je sais que je te demande de me masher un peut le travail mais je t'avoue que je suis complètement largé à ce sujet. Merci d'avance pour ton aide.

[kernadec]

bonsoir Metalos
pour cela, il y a les coordonnées Lambert pour le calcul des distances utilisées dans les cartes IGN,
ou la méthode Mercator de projection de la sphère a plat.

mais je te conseil un bout de code api java dans ton code tu auras la précision du gps..
que l'on peut intégrer avec Purebasic, il y des exemples sur le forum anglais et allemand
il sont différent mais il sont très instructif sur la façon de procéder.
pour infos "hjbremer" du forum allemand aime bricoler les apis google maps

sinon voila l'adresse des api pour le calcul des distances
https://developers.google.com/maps/docu ... rixService
https://developers.google.com/maps/docu ... ncematrix/ en francais
Cordialement

[MetalOS]

Je vais regarder tous ca. Merci Kernadec

[kernadec]

bonjour MetalOS
pour le calcul de l' orthodromie ou distance entre deux point sur la sphère..

Info-Wiki Le coefficient 60 devant l'arc cosinus provient du fait que le mile marin (1852 m) correspond à une minute d'arc
sur un grand cercle de la surface terrestre, et donc 60 milles

((60*Degree(ACos(Sin(Radian(Lat1))*Sin(Radian(Lat2))+Cos(Radian(Lat1))*Cos(Radian(Lat2))*Cos(Radian(Lg1-Lg2))))))*1.852

calcul de la distance Paris à Marseille en orthodromie * (1.852km) = Mile nautique

Code : Tout sélectionner

Debug  ((60*Degree(ACos(Sin(Radian(48.80))*Sin(Radian(43.39))+Cos(Radian(48.80))*Cos(Radian(43.39))*Cos(Radian(2.30-5.39))))))*1.852

calcul de la distance Paris à new York en orthodromie * (1.852km) = Mile nautique l
attention la longitude de NewYork étant négative et celle de Paris positive j'ai utilise une addition

Code : Tout sélectionner

Debug  ((60*Degree(ACos(Sin(Radian(48.80))*Sin(Radian(40.72))+Cos(Radian(48.80))*Cos(Radian(40.72))*Cos(Radian(74+2.30))))))*1.852

Cordialement

[MetalOS]

Merci de ton aide une fois de plus Kernadec. Peut on appliquer ces calcul au code que j'ai posté 2 message plus haut ?

[kernadec]

yes, c'est bon
le programme retourne des coordonnées géographiques, c'est avec ces données que la formule fonctionne.
mais c'est la distance à vol d'oiseau basée sur une circonférence de la terre de 40 000km

Attention les coordonnées de longitude ouest (W) de Greenwich sont toujours négative.

Cordialement

[GallyHC]

Bonjour,

Ce sujet est très instructif et j'avous qu'il me sera très utile . Pour la distance entre deux point j'avais fait une petite routine ici http://koakdesign.info/php/form_sources.php?num=20.

Cordialement,
GallyHC

[kernadec]

bonjour GallyHC
mile = 1.6093470878864446 qui ne correspond pas à une minute d' arc
attention.. ta routine n'est pas avec le mile nautique = 1852m environ
pour comprendrez que c'est le mile nautique voilà la raison.

Code : Tout sélectionner

;rayon moyen de la terre = 6366 km
;perimetre = 6366*2*#PI = 39998.757665505247
perimetre.d = 6366*2*#PI 
; 1 degré = 39998.7576 / 360 = 111.10766018195902km
degre.d = perimetre / 360
; 1 minute d'arc = 111.10766018195902 / 60 = 1.8517943363659837 km
minute.d = degre / 60
Debug perimetre
Debug degre
Debug minute 
;après arrondi de tout ça à 40 000 / 360 / 60 = 1.8518518518518519 ou 1852 m 

donc avec ce calcul pour 1000km, on est a plus ou moins 150m de la distance à vol d'oiseau

Cordialement

[MetalOS]

Tu ne serais pas prof ou employé chez Google Map Kernadec ? En étant sérieux, tu maîtrise quand même bien le sujet.