Bonjour à tous !

C'est la première fois que je poste sur le forum, j'espère ne pas me tromper en postant dans cette section.

Voici mon problème. Je cherche à réaliser un programme qui serait en mesure de prendre une carte en projection Mercator en entré pour en extraire automatiquement les différentes zones en fonction du code couleur appliqué (localisation et estimation de la surface).

J'ai par exemple récupérer une carte d'un fichier KML (KMG compressé) destiné à être exploité dans Google Earth comme un calque. Je placerai l'image à la fin de mon post.
À partir de là, je me suis appuyé sur les formules de projection Mercator disponibles un peu partout sur le web et notamment sur la page wiki : https://fr.wikipedia.org/wiki/Projection_de_Mercator

J'ai compris qu'il existait plusieurs projections normées se basant sur le procédé Mercator donnant lieu à différentes représentation. Les positions X et Y doivent donc être pondérées selon les dimensions de la carte source.

J'ai donc appliqué les mêmes méthodes de calcul que l'on retrouve à la fin de l'article wiki, mais je retrouve des décalage important, surtout en Y.

Malgré mes efforts, je ne parviens jamais à un niveau de précision suffisant pour pointé l'emplacement XY à partir de géolocalisation latitude et longitude.

Par exemple :

New York : long=-74° ; lat=40.4° ; h=3240 ; w=6480
x = w * ((long+180)/360) = 6480 * ((180-74)/360) = 1908 px
yrad = lat * (Pi/180) = 40.4 * (Pi/180) = 0.705113
y = h/2 - w/(2*Pi) * LN(TAN(Pi/4 + yrad/2)) = 824 px
Normalement (environ) : x=1913 et y=880

Djibouti : long=43.08° ; lat=11.35° ; h=3240 ; w=6480
Résultat du calcul : x=4015 et y=1414
Normalement (environ) : x=4011 et y=1396

George Island (Malouine) : long=-59.44° ; lat=52.22° ; h=3240 ; w=6480
Résultat du calcul : x=2170 et y=2133
Normalement (environ) : x=2165 et y=2563

Wellington : long=174.46° ; lat=41.17° ; h=3240 ; w=6480
Résultat du calcul : x=6380 et y=2425
Normalement (environ) : x=6398 et y=2361

D'après ce que j'ai compris, il s'agit d'un problème au niveau du facteur avant LN(TAN(...)), malgré tout je ne parviens pas à définir précisément la cause.

Sauriez-vous me dire quelle est la source de l'imprécision que je rencontre ?

Merci !

---

Voici le fichier extrait du KML que j'utilise :

-
Edité par Athanor68 14 janvier 2021 à 2:21:14

Bonjour ! La carte que tu montres n'est pas une projection de Mercator.

Est-ce que tu as un moyen de trouver de quelle projection il s'agit ?

-
Edité par robun 9 janvier 2021 à 16:42:59

Bonjour Robun

Merci pour ta réponse !

Au vu de la représentation générale, j'en ai déduit qu'il s'agissait d'un représentation Mercator adaptée pour être appliquée comme calque dans Google Earth.

Je suis entrain de regarder sans parvenir à trouver de réponse quant au fait qu'il s'agisse bien d'une projection normée. Je vais continuer à chercher car puisque l'image s'adapte parfaitement à la représentation sphérique de la terre, je suppose que le format est effectivement construit sur une norme quelconque.

----------------

J'ai fais mes petites recherches et je pense que la carte sur laquelle je suis partie est en réalité une projection cylindrique équidistante dites quadratique. Article wiki : https://fr.wikipedia.org/wiki/Projection_cylindrique_%C3%A9quidistante

Voici les formules :

En revanche je ne comprend pas comment exploiter les équations car la longitude 0 et la latitude 0, ne devrait pas justement valoir 0 ?

-
Edité par Athanor68 10 janvier 2021 à 0:48:43

Oui, je crois que tu as la bonne projection. Elle a un nom compliqué, mais c'est la plus simple ! (J'aurais dû la reconnaître immédiatement : c'est une de celles que j'utilise... Il est vrai que je fais dans la cartographie céleste, avec des constellations à la place des continents.) L'article de Wikipédia explique bien l'intérêt de \(\lambda_0\) et \(\varphi_0\) :

• Pour afficher la carte d'une région particulière et non du globe au complet (dans l'article ils prennent l'exemple de la France), il faut utiliser une latitude de cette région : \(\varphi_0\).
• Pour afficher le globe au complet, on prend en général \(\varphi_0 = 0\) . Par contre on n'est pas obligé de centrer la carte sur le méridien de Greenwich, d'où l'utilité de \(\lambda_0\).
• La carte que tu utilises est centrée sur le méridien de Greenwich, donc \(\lambda_0 = 0\). Du coup c'est la projection la plus simple qui soit : x = longitude, y = latitude (à un facteur d'échelle près).

-
Edité par robun 10 janvier 2021 à 13:31:18

Merci pour ton aide Robun. En exploitant ce que tu m'as dis j'ai défini les équations suivantes :

Lg = longitude en ° ; Lt = latitude en ° ; w = largeur de l'image ; h = hauteur de l'image

x = Lg / 180 * (w/2)

y = Lt / 89.59 * (h/2)

J'ai défini 89.59° comme latitude maximale pour coïncider avec Google Earth qui n'affiche pas jusqu'à 90°. Cette petit adaptation me permet d'obtenir des résultats plus précis. J'ai donc pu comparer mes positions calculées avec celle dans Google Earth et malheureusement je continue à obtenir des décalages parfois très limités, mais parfois suffisamment important pour placer une ville en pleine mer

New York : lat=40,4° ; long=-74° ; h=3240 ; w=6480

Résultat du calcul : x=-1332 et y=731

J'obtiens un décalage de -2px en y.

Djibouti : long=43.08° ; lat=11.35° ; h=3240 ; w=6480
Résultat du calcul : x=775 et y=205
J'obtiens un décalage de -3px en y et +2px en x.

George Island (Malouine) : long=-59.44° ; lat=52.22° ; h=3240 ; w=6480
Résultat du calcul : x=-1070 et y=-944
J'obtiens un décalage de -2px en x.

Wellington : long=174.46° ; lat=41.17° ; h=3240 ; w=6480
Résultat du calcul : x=3140 et y=-744
J'obtiens un décalage de -6px en x.

Pour Wellington, voici ce que j'obtiens depuis Google Earth.

Alors que le calcul nous amène en pleine mer (entre les deux îles de nouvelle Zélande ).

Toi qui a l'habitude de manipuler ce genre de carte, saurai-tu me dire d'où proviennent ces imprécisions ?

@Athanor68 Bonsoir, je viens de retirer des spams votre dernier message, si cela ce produit encore vous pouvez poster dans ce sujet Si votre message est considéré comme spam

-- Message supprimé --

-
Edité par Athanor68 10 janvier 2021 à 21:24:07

Est-ce que tu as bien converti les angles sexagésimaux en angles décimaux ?

Par exemple la longitude de Wellington est de 174°46' c'est-à-dire 174,77°. Je crois que ça pourrait expliquer les petites différences.

New York a une latitude de 40°43' soit 40,72°. Oui, on dirait bien que l'erreur est ici.

-
Edité par robun 11 janvier 2021 à 9:57:25

Bonsoir Robun.

Bien joué ! Effectivement tu as mis le doigt dessus, l'imprécision provenait du fait que le travaillait avec sur une base 60 au lieu d'une base 100. Du coup, j'obtiens précisément les mêmes positionnements que sur Google Earth.

Pour finaliser mon travail j'ai cherché à calculer des superficies à partir des pixels de la carte. Pour cela j'ai commencé par définir la superficie pour 1 pixel au niveau de l'équateur pour ne pas tenir compte de la déformation et j'obtiens :

w = 6480 pixel ; h = 3240 pixel ; Périmètre méridional = 40007.86 km ; Périmètre équatorial = 40075.02 km

Largeur d'un pixel = 40075.02 / w = 6.184417 km

Hauteur d'un pixel = 40007.86 / 2 / h  = 6.17405 km

Donc, à l'équateur un pixel vaut 6.184417 * 6.17405 = 38.18291 km²

Étant donné la déformation de la carte au niveau des pôles, il faut appliquer un coefficient pour compenser la déformation en X. Dans mes recherches j'ai découvert que la déformation se définissait par cos(latitude en rad). Comme l'image ci-dessous nous en informe à 50° de latitude, il y a une déformation de 0.642.

Donc, puisque la projection cylindrique équidistante déforme les surfaces en X en les étirant, j'ai fais le calcul suivant :

lp (longueur pixel) = 6.184417 ; hp (hauteur pixel) = 6.17405 ; cl (coefficient latitudinal)

Surface en km² = hp * (lp*cl)

Seulement, lorsque j'applique ces calculs mes résultats sembles très approximatifs. Par exemple, l'île "Gran Canaria" fait environ 1560 km² et sur la carte, elle est recouverte par 58 pixels. Pour simplifié l'exemple, j'ai calculer un coefficient latitudinal moyen de 0.883397 (27.94°N). D'autre part, j'ai repris la longueur d'un pixel valant 6.184417 km et sa hauteur 6.17405 km.

Donc La superficie calculée donne :

6.17405 * (6.184417 * 0.883397) * 58 = 1956.9 km²

Étant donné l'écart de plus de 25% que j'observe, je ne sais pas si cela est dû à une erreur de calcul ou du fait de la résolution trop grossière. En même temps, comme tu peux le constater sur l'image ci-contre, les pixels recouvrent l'île sans excès.

Que penses-tu de ma démarche ? Y aurait-il un problème au niveau de mes calculs ou du coefficient latitudinal peut-être ?

Merci !

-
Edité par Athanor68 13 janvier 2021 à 1:35:42

Quand tu disposes 58 petits carrés cote à cote, en forme plus ou moins circulaire, tu en as environ 31 qui sont à l'intérieur du disque, et le reste ( 27 environ) qui sont à la frontière.

On peut considérer que les 27 carrés qui forment le contour, ils sont pour moitié dans l'ile, et pour moitié sur la mer.

Le surface est donc de 31+27/2 pixels.  Et ça nous donne une surface pour l'île de 1501km².

Comment on peut trouver ce prorata (58 carrés en tout donnent 31 à l'intérieur, et 27 sur le contour).

Je pars du principe que la forme est plus ou moins circulaire. Plus on s'éloigne d'un cercle, plus c'est faux. Pour une longueur de périmètre imposée, le cercle donne la surface la plus grande qu'on puisse obtenir.

Un cercle qui a une surface de 58 carrés, il a un rayon de racine(58/3.14), soit 4.3 

Il a donc une circonférence de 2*3.14*4.3=27. Donc 27 carrés sur le contour.  Puis par différence, 31 sur l'intérieur.

Et donc (31+27/2) comme surface, si on estime que tous les pixels sur le contour doivent être comptés pour moitié.

C'est une approximation grossière. (mathémétiquement, j'ai des petits doutes, mais je pense que l'erreur faite nous arrange bien ; l'erreur, si erreur il y a , ferait un résultat un peu sous-estimé pour une forme parfaitement circulaire, mais peu surestimée pour une forme autre)

Bien évidemment, il faudrait une résolution plus fine pour améliorer le résultat. 

La superficie que tu as calculée, elle correspond à tout ce qui est en rouge sur l'image ? Si oui, la zone en rouge semble correspondre de très près à l'île. (Cela dit on voit mal les pixels sur cette image. On les voyait bien mieux sur l'image de la Nouvelle Zélande et, là, l'explication de tbc92 semble en effet expliquer la surestimation.)

Le principe du calcul m'a l'air correct à part un petit détail : plus haut tu disais que la latitude maximale était 89,59°. Tu voulais dire 89°59' ? Si c'est 89,59°, dans le calcul de la hauteur d'un pixel en km il ne faut pas utiliser le périmètre méridional mais une valeur un poil plus petite . (Si c'est 89°59' la correction existe aussi mais sera vraiment minuscule.)

Est-ce que tu as la possibilité de distinguer les pixels situés à l'intérieur et ceux situés à cheval sur la frontière ? Si oui, on pourrait encadrer la superficie par la somme des pixels intérieurs et la somme {pixels intérieurs + pixels à la frontière}.

-
Edité par robun 13 janvier 2021 à 13:28:40

Merci pour vos réponses, à toi tbc92 et surtout à toi robun ! Vôtre aide m'a permit de calculer tout ce dont j'ai besoin pour finaliser mon projet.


Vous avez tout deux mit dans le mille, en effet le problème de précision provenait des pixels en lisières d'autres surface qui, de fait gonflait les surfaces calculées.


> Robun : La mesure de 89°59' dans mes calculs était bien une erreur d'interprétation. Je n'avais par encore compris la conversion des angles sexagésimaux en décimaux à ce moment là. Google Earth génère bien une représentation de la terre jusqu'au pôle.

Pour finir en beauté, je propose un résumé des solutions que vous m'avez permit de découvrir qui servira peut-être à d'autre curieux de la cartographie

Encore merci pour tout !

------------------

Projection cylindrique équidistante dites quadratique

Carte de travail

Coordonnées latitude et longitude vers la coordonné d'un pixel

1. Passer les coordonnés latitude et longitude en négatif pour les latitudes Sud et les longitudes Est

2. Conversion du sexagésimal au décimal
Formule : lgs = longitude en ° sexagésimal ; lts = latitude en ° sexagésimal

longitude en ° décimal (lgd) = lgs+((lgs-arrondi(lgs))/60*100)

latitude en ° décimal (ltd) = lts+((lts-arrondi(lts))/60*100)

Exemple :    

Ville                   Lat °60        Long °60    Lat °100    Long °100
New York          40°43'        -74°             40,71         -74
Djibouti             11°35'        43°08'         11,58         43,13
George Island   -52°22'       -59°44'       -52,37        -59,73
Wellington        -41°17'       174°46'       -41,28        174,77

3. Calculer les coordonnées X et Y d'un pixel
Dans notre exemple la carte fait 6480 * 3240 pixels

Formules :    lgd = longitude en ° (décimal) ; ltd = latitude en ° (décimal) ; w = largeur de l'image = 6480 ; h = hauteur de l'image = 3240
    x = arrondi((lgd/180)*(w/2))
    y = arrondi(ltd/90*(h/2))

Exemple :

Ville                  x                 y    
New York          -1332         733
Djibouti             776            208
George Island  -1075         -943
Wellington        3146         -743

Note : Les positions x et y prennent pour point de départ le centre de la carte au coordonné (0°, 0°).
Attention : Toutes les cartes ne sont pas centrées sur les coordonnées 0°,0°, dans ce cas il faut faire la manipulation suivante :

Coordonnées d'un pixel vers une latitude et une longitude

1. Passer les coordonnés latitude et longitude en négatif pour les latitudes Sud et les longitudes Est
2. Calculer les coordonnées longitudinales et latitudinales décimal
Formule : w = largeur de l'image = 6480 ; h = hauteur de l'image = 3240
lgd = (x/(w/2))*180
ltd = (y/(h/2))*90

Exemple :

Ville                   x            y         long en °100    lat en °100
New York          -1332     733     -74                   40.7222
Djibouti             776        208     43.1111           11.5555
George Island   -1075    -943    -59.713            -52.3727
Wellington        3146      -743    174.7238         -41.2523
            
3. Conversion du décimal au sexagésimal
lgs = lgd+((lgd-arrondi(lgd))/100*60)
lts = ltd+((ltd-arrondi(ltd))/100*60)

Note  : Pour les valeurs négatives, pensez à arrondir vers le haut plutôt que vers le bas.

Exemple :

Ville                    x            y        long en 60    lat en 60
New York           -1332     733    -74               40.43
Djibouti              776        208    43.07           11.33
George Island    -1075    -943   -59.43          -52.22
Wellington         3146      -743   174.43         -41.2523

Calculer une surface

1. Calculer la surface d'un pixel par rapport à la superficie de la terre
Formules : w = 6480 pixel ; h = 3240 pixel ; Périmètre méridional = 40007.86 km ; Périmètre équatorial = 40075.02 km
Largeur d'un pixel = 40075.02 / w = 6.184417 km
Hauteur d'un pixel = 40007.86 / 2 / h  = 6.17405 km
Superficie d'un pixel = 6.184417 * 6.17405 = 38.18291 km²

2. Calculer la latitude en radian
Formule : ltd : latitude en ° (décimal)
ltrad = ltd*(Pi/180)

Exemple :     Ville            lat en °100        lat en rad
                    New York        40.71            0.71052354
            
3. Calculer la déformation en x selon la latitude
Formule : coefficientLatitude = COS(ltrad)

4. Calcul de la superficie
Note : les pixels ne permettent pas une précision suffisante pour délimiter parfaitement les changements de zone d'une carte, il faut prendre en compte ce paramètre pour obtenir des superficies proches de la réalité.

Exemple avec l'île "Gran Canaria" faisant au total 1560 km² :
Celle-ci est recouverte par 61 pixels, 30 à l'intérieur de la zone, et 31 en bordure. Le coefficient latitudinal moyen est de 0.88363 et un pixel au niveau de l'équateur représente une superficie de 38.18291 km².
La formule pour les 30 pixels de la surface interne (si) est donc :
    si = 30 * (38.18291*0.88363) = 1012.19 km²

Note : En condition réelle il faut calculer un coefficient latitudinal pour chaque coordonnées Y traitées. Le coefficient moyen sert uniquement à simplifier notre exemple.

Les pixels à la lisière de l'océan doivent être traité spécifiquement pour ne considérer qu'une fraction de la surface qu'ils représentent.
Sur l'image ci-contre, voici comment j'ai procédé :
- les pixels en gris claire possède 2 ou 3 côté en contact avec l'océan (surface bleu) : Je ne prendrais que 25% de la surface totale du pixel à cette latitude.
- les pixels en gris foncé possède 1 côté en contact avec l'océan : Je ne considérerai que 50% de la surface totale du pixel
- les pixels en noir possède un coin en contact avec l'océan : je considère pour ces derniers la totalité de la surface du pixel comme s'il s'agissait de pixels internes.

Nous obtenons donc : 8 pixels noir, 12 pixels gris foncé et 11 pixel gris claire.
La formule pour calculer la surface total est donc : si = surface interne ; st = surface totale
st = si + 8*(38.18291*0.88363) + 12*(38.18291*0.88363)*0.5 + 11*(38.18291*0.88363)*0.25 = 1552.02 km²

Attention : Il faut penser à attribuer aux surfaces connexes ce que l'on aura déduit en ne considérant qu'une fraction de la surface des pixels en bordure.

Note : Pour améliorer le calcul, il faut améliorer la résolution de la carte afin d'obtenir un grain plus fin.

Bon courage à tous dans vos projets !

-
Edité par Athanor68 15 janvier 2021 à 22:52:43

Super compte-rendu final ! 

L'idée de pondérer les pixels frontières est très bonne. Si un pixel frontière a 1 seul côté en contact avec la mer, ou 2 côtés en contact avec la mer, ce n'est pas la même chose, et c'est une bonne idée de ne pas appliquer le même coefficient

Je suis un peu jaloux de ne pas y avoir pensé avant.

Nos idées se complètes

J'ai toujours été sidéré de voir ce que sont capables les individus lorsqu'ils mettent leurs connaissances en commun.. vive la complémentarité !