Salut,

je viens de lire votre article “Envoyez nous les scripts que vous ne pigez pas”.
L’article date un peu, ptêt j’arrive trop tard.
Dans le cas contraire.. read on..

Voilà: je me ronge le cerveau sur un script python depuis un moment déjà, mais il reste aujourd’hui une lamentable énigme.

La situation:

On m’offre l’intégrale (26Go!) des cartes IGN, tombées du cyber-camion, accompagné de quelques lignes de python permettant de soit-disant trouver la carte qui nous intéresse par ses coordonnées GPS et de, une fois qu’on a une carte, trouver celles qui viennent se placer au Nord/Sud/Est/Ouest.

D’la balle, que j’me dis. J’y connais rien au python, mais le code est lisible, commenté, easy quoi.
Après quelques heures de lecture, le monde sera à moi.. C’était il y a deux mois.
Échec total de communication entre mes synapses et le serpent.
Jamais plus je ne toucherai à du Python me suis-dis alors.
Jusqu’à ce que l’espoir renaisse à lecture de votre billet qui-fut-pour-moi-comme-un-arc-en-ciel (et hop, léchage de cul qui fera, à n’en point douter, toute la différence).

Je m’en remets donc à vous, nobles aventuriers, et vous remets votre carte au trésor ci-dessous:

* find_map_by_coordinate.py
http://0bin.net/paste/8a658d3a8f3149ee3af36b7f9abb0e291d00d3f1#OYNAOG9Bhuk/cdC6SgsgVxoW+s5eADhCwjUGACNI4wk=

* find_close_maps.py
http://0bin.net/paste/5a1a9236cddfbea5abbc4c21459856acfd21fd17#MncBXtbHcQFaKW8l3aEVmyNc+0gqe4/iTU6aIatAK4w=

* ign25.py
http://0bin.net/paste/f50e61ad54b9e53d6dd99a655d9cf564eb7a5ad7#QWa6Fz9GejLNK5+niZNqc6/fTm0Pwvl2fBjtPGL8I5Y=

Et pour répondre à votre question, oui, je peux vous envoyer les précieuses cartes.
Reste simplement à trouver un moyen de vous faire passer 26Go de données..

Bon décodage!

Bien à vous,

(pour partouze c’est toujours OK ?)

Bon, d’accord, j’ai rajouté la dernière ligne.

Donc, j’ai lu les scripts, et je les ai commenté au fil de l’eau. Vous allez voir, ça se sent, c’est assez organique et on sent bien se développer mes réactions au fur et à mesure face à ce script.

Spoiler : cher ami lecteur, tes synapses ne sont pas en cause, l’auteur du script n’est pas un Dieu de la prog, et en fait je ne suis même pas sûr que son code marche. Même un petit peu. Auteur anonyme si tu nous lis, toutes mes condoléances.

Comme ça va être long, et conforme à la tradition :

find_map_by_coordinate.py

Là je découvre, c’est l’aventure, on sent l’enthousiasme.

#!/usr/bin/python


# Ceci est un script de recherche de coordonnées, il trouve une carte qui
# contient ces coordonnées dans la base de données des cartes. Le but est donc
# de prendre en entrée des coordonnées sous forme de latitudes et longitudes
# saisies par l'utilisateur et d'afficher le nom des cartes sur son terminal.
#
# Ex:
#
# $ python ./find_map_by_coordinate.py (5, 32) (3, 44)
# Trying to open database...
# ...OK
# Searching for (104, 32),  (201, 44)...
# ...found.
# ['/chemin/vers/carte1', '/chemin/vers/carte2']

# Comme d'habitude on a les imports en premier, ici l'auteur a choisi
# de mettre les imports de ses propres modules avant ceux de la lib
# standard, je recommande l'inverse.

from ign25 import load_db, find

import os
import sys
import ast

# Plutôt que de uniquement des prints, l'auteur choisit d'utiliser
# le module de logging, mais l'utilise au final pour faire des prints.
# Du coup je ne suis pas certains de l'interêt de la chose, peut être
# que l'intention est de le faire plus plus tard.
# Dans tous les cas il met le niveau de logging sur DEBUG, ce qui fait
# que tous les messages seront affichés.
import logging
logging.basicConfig(level=logging.DEBUG)

# On a ici l'idiome connu qui permet de n'exécuter le bloc que si
# le script et lancé directement (et pas importé). Dans ce cas c'est inutile,
# en effet il n'y a rien dans ce script qu'on puisse vouloir importer, c'est
# donc un script qui est TOUJOURS destiné à être lancé. Encore une fois
# c'est peut être un indice d'une intention future de l'auteur.
if __name__ == "__main__":

    # Le parsing des arguments est fait à la main. Pour un petit script
    # comme ça ça a du sens (on a pas toujours envie de sortir l'artillerie
    # lourde avec argparse)
    # En gros, si il y a les deux arguments nécessaires (sys.argv contient
    # le nom du script et chaque argument passé au script), on les traite,
    # sinon on affiche un message indiquant comment utiliser le script.
    # Le traitement des arguments est une évaluation du code comme si c'était
    # du Python. Encore une fois je ne recommande pas cela, demander les
    # degrés et les minutes sous la forme naturelle (34°27') et le parser
    # avec des regex est plus propre ou de passer des longitudes et latitudes
    # via 4 paramètres au lieu de deux ici.
    if len(sys.argv) >= 3:
        x = ast.literal_eval(sys.argv[1])
        y = ast.literal_eval(sys.argv[2])
    else:
        print("find map for given coordiantes")
        print("")
        print("usage: %s [lat] [lon]" % sys.argv[0])
        print("       with [lat,lon] in degrees or in (degrees, minutes)")
        sys.exit(0)

    # Ici on fait appel à la fonction load_db() qui se connecte à la
    # base de données
    db = load_db()

    # On lance la recherche dans la base de données avec la fonction find()
    # qui est aussi issue du module ign25, et on affiche les résultats.
    logging.info("Searching for %s, %s..." % (x, y))
    r = find(db, x, y)
    logging.info("...found.")
    print(r)

find_close_maps.py

Je tique un peu. Il y a des trucs bizarres. Mais bon, c’est peut être moi.

#!/usr/bin/python

# Même principe, mais pour les cartes proche d'une autre carte. Je ne vais
# donc que commenter les différences.
# Il attend un fichier .map en paramètre et va lister toutes les cartes de
# la base de données qui sont à 10km de la carte dont le fichier a été passé
# en paramètres.


import os
import os.path
import sys
import math

import logging
logging.basicConfig(level=logging.DEBUG)

from ign25 import *


if __name__ == "__main__":

    if len(sys.argv) < 2:
        print "find maps around a given map"
        print ""
        print "usage: %s [reference map]" % sys.argv[0]
        sys.exit(0)

    db = load_db()

    # On ouvre le fichier de carte, et on récupère les coordonnées des 4
    # extrémités du carré que forme la carte sous le forme d'un objet de type
    # MapInfo qui contiendra les attributs nw, ne, se, sw (Nord-Ouest, Nord-Est,
    # Sud-Est, Sud-Est)
    try:
        reference = get_info(open(sys.argv[1] + ".map"))
    except:
        print "map %s not found" % sys.argv[1]
    maxDistance = 10 #in km
    # la distance maximale est écrite en dure, on pourrait la passer en paramètre

    # Ensuite la recherche se fait avec encore une fois une fonction de
    # ign25, donc la logique de ce script est assez simple puisque la
    # complexité est dans ign25.py. On récupère une liste de carte.
    closemaps = gen_closemaps(db, reference, maxDistance)

    # La liste est affichée mais on remplace la carte par une image de la carte
    # au format PNG (qui je suppose est fournie avec la carte ?)
    for i in closemaps:
        print os.path.abspath("../png/" + i.filename + ".png")

ign25.py

... yet so far

#!/usr/bin/python

"""Module to work with my set of IGN TOP25 maps

"""

# Toute la logique des scripts précédents est définie ici.
# Encore une fois la docstring est un peu légère et des commentaires
# en plus ne serait pas de refus car du coup pour comprendre ce que fais
# le script juste en les recevant par email, c'est moins facile, forcément ^^

# Encore une fois, je vous invite à ne pas faire plein d'imports sur une ligne
# si ils ne sont pas du même module. D'autant qu'ici les modules shutils et
# math ne sont pas utilisé dans le script, et os.path est importé deux fois.
# Cher auteur de script, loin de moi l'idée d'insulter ton travail, j'ai
# conscience que tu es sans aucun doute un géographe et non un développeur,
# ne m'en veut donc pas, je ne fais que commenter le script à des fins
# pédagogiques et pas du tout pour te descendre. 
import os, re, math, shutil, os.path, sys
import ast
import logging
import os.path
import math

# Déclaration de constantes. Mettez les constantes en majuscules car la notion
# de constante n'existe pas en Python et seule cette convention permet
# de savoir qu'il ne faut pas y toucher.

easting2km = 0.73847
northing2km = 0.53848
northAdjust = easting2km / northing2km

# Une classe vide qui va servir de conteneur. Je comprends que si l'on vient
# d'un langage avec des structs, c'est une tendance naturelle. Néanmoins,
# si il n'y a aucune logique associée à vos données, préférez un dictionnaire
# ou un namedtuple.
class MapInfo:
    """Store map information"""
    pass

# Ah, un peu de doc, ça fait plaisir :-)

# Routines to read .map files

def get_coordinates(line):
    """Read the coordiantes of a point

    line:       line describing the point
    returns:    a tuple of tuples ((xdeg, xmin), (ydeg, ymin))"""
    try:
        # La fonction analyse une ligne de texte formattée ainsi:
        # champ1,    champ2    , champ3    ,    champ4...
        # et en extrait les degrés et les minutes avec un split sur
        # le caractère virgule entouré par un nombre indéterminé de caractères
        # non imprimables (espaces, tabulations...)
        fields = re.split( "\s*,\s*", line )
        # adjust for East, West. We are always in the northern hemisphere.

        # Detecte si on est dans l'hémisphèse Nord (je suppose que le onzième
        # champ est égale à "W" dans ce cas, et "E" dans l'autre cas mais
        # je n'ai pas le format sous les yeux pour le vérifier).
        # Bref, si on est dans l'hémisphère Nord, on donne aux champs une valeur
        # négative. Je suppose que ça parle aux amateurs de GIS, comme je
        # suis incompétent dans le domaine, je ne vais pas faire mon malin
        # et tenter d'expliquer pourquoi.
        if fields[11] == "W":
            fields[9] = "-" + fields[9]
            fields[10] = "-" + fields[10]

        # La fonction retourne ensuite un tuple de deux tuples avec les valeurs
        # ou None si il a eu une erreur.
        # Je ne sais pas où il peu y avoir une ValueError là dedans, donc
        # difficile de dire quelle erreur l'auteur cherche à gérer.

        return ((int(fields[6]), float(fields[7])), (int(fields[9]), float(fields[10])))
    except ValueError:
        return None

# La fonction utilisée dans find_close_maps.py, qui lit un fichier carte
# (attend un file like object en paramètre) et retourne un objet MapInfo().
# Comme je vous l'ai dis dans ce cas il vaudrait mieux retourner un dictionnaire
# ou un namedtuple.
def get_info(file):
    """Read map information in `file`

    file:       file object to read
    returns:    MapInfo"""
    i = MapInfo()
    # Je pense que le but recherche est d'extraire le nom du fichier sans
    # extention. os.path.splitext(file.name)[0] aurait été préférable
    i.filename = file.name[:-4]
    # Pour chaque ligne du fichier carte, on va vérifier si la ligne
    # contient la chaîne Point0x, extraire les coordonnées de la ligne
    # et les stocker dans un attribut représentant un des coins de la carte
    # de l'objet MapInfo().
    #
    # Pour les améliorations possibles :
    # - mettre le return dans la boucle for
    # - re.match(r'regex', line) suffit
    # - retourner None plutôt que False
    for line in file:
        if re.compile(r'Point03').match(line):
            i.nw = get_coordinates(line)
        elif re.compile(r'Point02').match(line):
            i.ne = get_coordinates(line)
        elif re.compile(r'Point04').match(line):
            i.sw = get_coordinates(line)
        elif re.compile(r'Point01').match(line):
            i.se = get_coordinates(line)
    if hasattr(i, 'nw'):
        return i
    else:
        return False
# Une manière concise d'écrire la même fonction aurait pu être:
# data = {'type': os.path.splitext(file.name)[0]}
# for line in file:
#     res = re.search('Point0(?P1)?(?P2)?(?P3)?(?P3)?', line)
#     if res:
#         res = (name for name, val in res.groupdict().items() if val).next()
#         data[name] = get_coordinates(line)
#         if len(data) == 4
#             return data
# return None
# Mais bon, la c'est de l'enculage de mouche, l'important c'est que ça marche.

# Cette fonction n'est pas appelée dans les autres scripts.
# Elle retourne un générateur qui retourne, pour chaque fichier passé
# dans l'itérable en paramètre (par exemple une liste de chemins de fichier)
# un générateur donc chaque élément est un objet MapInfo contenant les
# infos de la carte correspondant à chaque nom de fichier
def gen_info(filenames):
    for filename in filenames:
        try:
            logging.info("Reading %s" % filename)
            f = open(filename)
            info = get_info(f)
            if info:
                yield info
        except IOError:
            pass

# Pareil, pas utilisé. Ca fabrique une liste en listant les fichiers d'un
# dossier nomme "map" (on va supposé qu'il est déjà là ?), on passant
# cette liste à gen_info() qu'on vient de voir plus haut, et en transformant
# le générateur obtenu en liste.
def generate_database():
    if os.path.isdir('map'):
        mapdir = 'map'
    else:
        mapdir = '.'
    fnames = os.listdir (mapdir)
    return list(gen_info(fnames))

# Routine to print db

# Idem, pas utilisé. On passe une base de données en paramètres
# Et on print un listing de tous les fichiers de cartes de la BDD et
# leurs coordonnées.
def print_database(db):
    print 'filename,', 'nw,', 'ne,', 'sw,', 'se,'
    for i in db:
        # La seule partie compliquée : on définie une fonction
        # dans le corps de la boucle for et on l'utilise cash pistache.
        # Une manière plus propre aurait été :
        # vals = (i.filename, i.nw, i.ne, i.sw, i.se)
        # print ' '.join(str(x).replace(' ','') for x in vals)
        # Si on avait un namedtuple plutôt qu'un objet, une ligne aurait suffit
        def rms(s):
            return str(s).replace(' ','')
        print i.filename, rms(i.nw), rms(i.ne), rms(i.sw), rms(i.se)

# Itou, pas utilisé. On prend chaque carte dd'un dossier, et on écrit
# toutes les infos obtenues avec les fonctions précédentes dans un fichier.
# Encore une fois notez que si on avait utilisé un namedtupple, l'opération
# aurait pris 4 lignes.
# Il y a néanmoins 3 choses qui m'interpellent dans ce code :
# - pourquoi appeler ça une base de données ? C'est juste un listing texte
# - pourquoi ne pas avoir choisit le format CSV (il y a un module pour ça)
# - pourquoi writelineS() et pas writeline() ? Est-ce un bug ?
def save_database(db, path):
    f = open(path, 'w')
    f.writelines(('filename ', 'nw ', 'ne ', 'sw ', 'se ', '\n'))
    for i in db:
        def rms(s):
            return str(s).replace(' ','')
        f.writelines((i.filename, ' ',
                      rms(i.nw), ' ',
                      rms(i.ne), ' ',
                      rms(i.sw), ' ',
                      rms(i.se), '\n'))

# Routine to import db L'inverse de la fonction précédente. Chaque tout le
# fichier et retourne une liste d'objet MapInfos(). Je pense que l'auteur ne
# connais pas la fonction float()
def import_database(path):
    db = []
    f = open(path)
    f.readline()
    for line in f:
        fields = line.split()
        i = MapInfo()
        i.filename = str(fields[0])
        i.nw = ast.literal_eval(fields[1])
        i.ne = ast.literal_eval(fields[2])
        i.sw = ast.literal_eval(fields[3])
        i.se = ast.literal_eval(fields[4])
        db.append(i)
    return db

# Ca appelle import_database et si ça ne marche pas, ça
# génère une "base de données".
def load_db():
    try:
        logging.info("Trying to open database...")
        db = import_database("map.db")
        logging.info("...OK.")
    except IOError:
        logging.info("...failed. Reading .map files...")
        db = generate_database()
        logging.info("Saving database...")
        save_database(db, "map.db")
    return db

# Lat/lon to filename
# Une fonction pour te rendre bien deg.
# Elle converti un tuple ou une liste Latitude / Longitude en une valeur
# en degré. Malheureusement elle tue complètement le duck typing (elle pourrait
# accepter n'importe quel iterable). Une version pythonique serait :
# try:
#   deg, min = x # marche avec TOUT itérable
# except TypeError:
#   deg, min = x, 0
# return float(deg) + float(min)/60. # accepte tout ce qui est castable en float
#
def to_deg(x):
    # isinstance est rarement une bonne idée en Python. Très rarement.
    if isinstance(x, tuple) or isinstance(x, list):
        xdeg, xmin = x # un bon petit unpacking !
    else:
        xdeg = x
        xmin = 0
    return xdeg + xmin/60.

# Trouve une carte qui contienne ces coordonnées, retourne le filename
# plutôt que l'objet MapInfo() (pourquoi ?).
def find(db, x, y):
    target_xdeg = to_deg(x)
    target_ydeg = to_deg(y)
    #print "target", target_xdeg, target_ydeg
    for i in db:
        nw_x, nw_y = i.nw
        se_x, se_y = i.se
        nw_xdeg = to_deg(nw_x)
        nw_ydeg = to_deg(nw_y)
        se_xdeg = to_deg(se_x)
        se_ydeg = to_deg(se_y)
        # est-ce que les bordures de la carte contiennent les coordonnées ?
        if nw_xdeg <= target_xdeg and target_xdeg <= se_xdeg:
            #print "candidate x", nw_xdeg, se_xdeg
            if nw_ydeg >= target_ydeg and target_ydeg >= se_ydeg:
                #print "candidate y", nw_ydeg, se_ydeg
                return i.filename

# Close maps
# Vérifie qu'une carte est proche d'une autre carte. Le PEP8 recommanderait
# de l'appeler is_close().
def isclose(candidate, reference, max_distance):
    if candidate.filename == reference.filename:
        return False
    # put everything in minutes
    # On repasse toutes les coordonnées en minutes
    infoEastings = to_deg(candidate.nw[0])*60
    infoNorthings = to_deg(candidate.nw[1])*60*northAdjust
    referenceEastings = to_deg(reference.nw[0])*60
    referenceNorthings = to_deg(reference.nw[1])*60*northAdjust

    # Pythagore, les maths de 6eme, ça vous revient ?
    distance = ( math.sqrt(
        ( infoEastings - referenceEastings )**2 + ( infoNorthings - referenceNorthings )**2 )
          * easting2km )
    if distance < max_distance:
        return True
    else:
        return False

# un générateur qui yield toutes les cartes proches d'une distance
def gen_closemaps(db, reference, max_distance):
    for i in db:
        if isclose(i, reference, max_distance):
            yield i

# Tile
# Une fonction pas utilisée qui convertie les degré en..., heu... lat / long ?
# Merci le module math de Python qui permet de s'amuser avec la trigonométrie
def deg2num(zoom, lat_deg, lon_deg):
    """Lon./lat. to tile numbers"""
    lat_rad = math.radians(lat_deg)
    n = 2.0 ** zoom
    xtile = int((lon_deg + 180.0) / 360.0 * n)
    ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
    return (xtile, ytile)

# Pareil, pas utilisé. C'est la fête du slip. On fait des fonctions pour
# le fun, si ça se trouve on les importera un jour dans le shell, un dimanche
# de pluie
def num2deg(zoom, xtile, ytile):
    """Tile numbers to lon./lat."""
    n = 2.0 ** zoom
    lon_deg = xtile / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    lat_deg = math.degrees(lat_rad)
    return (lat_deg, lon_deg)

# obtenir une les map qui contienne les coordonnées suivantes
# Une tile est un morceau de carte qu'on compte utiliser pour faire une
# carte plus grosse en assemblant les parties
def get_tile(z, x, y):
    lon_deg, lat_deg = num2deg(z, x, y)
    logging.debug("conv deg: %f, %f" % (lon_deg, lat_deg))
    db = load_db()
    filename = find(db, lon_deg, lat_deg)
    return filename

# Une version plus simple de find_map_by_coordinate.py qui utilise
# if __name__ == "__main__" à bon escient
if __name__ == "__main__":

    if len(sys.argv) >= 3:
        x = ast.literal_eval(sys.argv[1])
        y = ast.literal_eval(sys.argv[2])
    else:
        x = (48, 41.29)
        y = (2, 22)

    db = load_db()

    logging.info("Searching %f, %f..." % (x, y))
    r = find(db, x, y)
    logging.info("...found: %s" % r)

En conclusion, "tout est bien qui part à point sous les ponts".