Rendre un globe SVG

Ce que vous allez créer

Dans ce tutoriel, je vais vous montrer comment utiliser une carte SVG et la projeter sur un globe, en tant que vecteur. Pour effectuer les transformations mathématiques nécessaires à la projection de la carte sur une sphère, nous devons utiliser un script Python pour lire les données de la carte et les traduire en image d'un globe. Ce didacticiel suppose que vous utilisez Python 3.4, le dernier logiciel Python disponible..

Inkscape dispose d'une sorte d'API Python qui peut être utilisée pour effectuer diverses tâches. Cependant, comme nous ne sommes intéressés que par la transformation de formes, il est plus facile d'écrire un programme autonome qui lit et imprime des fichiers SVG de manière autonome..

1. Formater la carte

Le type de carte que nous voulons s'appelle une carte équirectangulaire. Sur une carte équirectangulaire, la longitude et la latitude d’un lieu correspondent à sa X et y position sur la carte. Une carte mondiale équirectangulaire est disponible sur Wikimedia Commons (voici une version avec les États-Unis)..

Les coordonnées SVG peuvent être définies de différentes manières. Par exemple, ils peuvent être relatifs au point précédemment défini ou définis absolument à partir de l'origine. Pour rendre nos vies plus faciles, nous voulons convertir les coordonnées de la carte en forme absolue. Inkscape peut le faire. Accédez aux préférences d'Inkscape (sous l'onglet modifier menu) et sous Entrée sortie > Sortie SVG, ensemble Format de chaîne de chemin à Absolu.

Inkscape ne convertira pas automatiquement les coordonnées. pour que cela se produise, il faut transformer les chemins. La meilleure façon de le faire est de tout sélectionner et de le déplacer de haut en bas en appuyant une fois sur chacune des flèches haut et bas. Puis ré-enregistrez le fichier.

2. Lancez votre script Python

Créez un nouveau fichier Python. Importez les modules suivants:

import sys import re import math temps import math import import datetime import numpy en tant que np import xml.etree.ElementTree en tant que ET

Vous devrez installer NumPy, une bibliothèque qui vous permet d’effectuer certaines opérations vectorielles telles que le produit point et le produit croisé..

3. La projection mathématique de la perspective

Pour projeter un point dans un espace tridimensionnel dans une image 2D, vous devez rechercher un vecteur de la caméra au point, puis diviser ce vecteur en trois vecteurs perpendiculaires.. 

Les deux vecteurs partiels perpendiculaires au vecteur de la caméra (la direction dans laquelle la caméra fait face) deviennent le X et y coordonnées d'une image projetée orthogonalement. Le vecteur partiel parallèle au vecteur caméra devient ce que l’on appelle le z distance du point. Pour convertir une image orthogonale en une image en perspective, divisez chaque X et y coordonner par le z distance.

À ce stade, il est judicieux de définir certains paramètres de la caméra. Premièrement, nous devons savoir où se trouve la caméra dans un espace 3D. Stocker ses X, y, et z coordonnées dans un dictionnaire.

camera = 'x': -15, 'y': 15, 'z': 30

Le globe sera situé à l'origine, il est donc judicieux d'orienter la caméra vers lui. Cela signifie que le vecteur direction de la caméra sera l'opposé de la position de la caméra.

cameraForward = 'x': -1 * caméra ['x'], 'y': -1 * caméra ['y'], 'z': -1 * caméra ['z']

Il ne suffit pas de déterminer la direction dans laquelle la caméra est dirigée; vous devez également effectuer une rotation de la caméra. Faites-le en définissant un vecteur perpendiculaire à la cameraForward vecteur.

cameraPerpendicular = 'x': cameraForward ['y'], 'y': -1 * cameraForward ['x'], 'z': 0

1. Définir les fonctions vectorielles utiles

Il sera très utile de définir certaines fonctions vectorielles dans notre programme. Définir une fonction de magnitude vectorielle:

#magnitude d'un vecteur 3D def sumOfSquares (vecteur): vecteur de retour ['x'] ** 2 + vecteur ['y'] ** 2 + vecteur ['z'] ** 2 valeur de magnitude (vecteur): mathématiques de retour .sqrt (sumOfSquares (vector))

Nous devons pouvoir projeter un vecteur sur un autre. Comme cette opération implique un produit scalaire, il est beaucoup plus facile d'utiliser la bibliothèque NumPy. NumPy, cependant, prend les vecteurs sous forme de liste, sans les identificateurs explicites 'x', 'y', 'z', nous avons donc besoin d'une fonction pour convertir nos vecteurs en vecteurs NumPy.

#convertit le vecteur dictionnaire à la liste vecteur def vectorToList (vecteur): retour [vecteur ['x'], vecteur ['y'], vecteur ['z']]
#projects u sur v def vectorProject (u, v): retourne np.dot (vecteurToList (v), vecteurToList (u)) / magnitude (v)

C'est bien d'avoir une fonction qui nous donnera un vecteur unitaire dans la direction d'un vecteur donné:

#get unité vecteur def unitVector (vector): magVector = magnitude (vector) return 'x': vector ['x'] / magVector, 'y': vector ['y'] / magVector, 'z': vector [ 'z'] / magVector

Enfin, nous devons pouvoir prendre deux points et trouver un vecteur entre eux:

#Calcule le vecteur à partir de deux points, forme du dictionnaire def findVector (origine, point): return 'x': point ['x'] - origine ['x'], 'y': point ['y'] - origine [ 'y'], 'z': point ['z'] - origine ['z']

2. Définir les axes de la caméra

Il ne reste plus qu’à finir de définir les axes de la caméra. Nous avons déjà deux de ces axes-cameraForward et cameraPerpendiculaire, correspondant à la z la distance et X coordonnée de l'image de la caméra. 

Maintenant, nous avons juste besoin du troisième axe, défini par un vecteur représentant la y coordonnée de l'image de la caméra. Nous pouvons trouver ce troisième axe en prenant le produit croisé de ces deux vecteurs, en utilisant NumPy-np.cross (vectorToList (cameraForward), vectorToList (cameraPerpendicular)).

Le premier élément du résultat correspond à la X composant; la seconde à la y composant, et le troisième à la z composant, donc le vecteur produit est donné par:

#Calcule le vecteur dans le plan de l'horizon (points en haut) cameraHorizon = 'x': np.cross (vectorToList (cameraForward), vectorToList (cameraPerpendicular)) [0], 'y': np.cross (vectorToList (cameraForward), vectorToList (cameraForward), vectorToList (cameraPerpendicular) ))) [1], 'z': np.cross (vecteurToList (cameraForward), vecteurToList (cameraPerpendiculaire)) [2]

3. Projet d'orthogonal

Pour trouver l'orthogonal X, y, et z distance, nous trouvons d’abord le vecteur reliant la caméra au point en question, puis nous le projetons sur chacun des trois axes de la caméra définis précédemment:

def physicalProjection (point): pointVector = findVector (caméra, point) #pointVector est un vecteur partant de la caméra et se terminant à un point de la question renvoyée 'x': vectorProject (pointVector, cameraPerpendicular), 'y': vectorProject (pointVector , cameraHorizon), 'z': vecteurProject (pointVector, cameraForward)

Un point (gris foncé) projeté sur les trois axes de la caméra (gris). X est rouge, y est vert et z est bleu.

4. Projet à perspective

La projection en perspective prend simplement la X et y de la projection orthogonale et divise chaque coordonnée par le z distance. Cela fait en sorte que les choses plus éloignées aient l'air plus petites que celles qui sont plus proches de la caméra. 

Parce que diviser par z donne de très petites coordonnées, on multiplie chaque coordonnée par une valeur correspondant à la focale de la caméra.

longueur focale = 1000
# dessine des points sur le capteur de la caméra à l'aide de xDistance, yDistance et zDistance def perspectiveProjection (pCoords): scaleFactor = focalLength / pCoords ['z'] return 'x': pCoords ['x'] * scaleFactor, 'y': pCoords [ 'y'] * scaleFactor

5. Convertir les coordonnées sphériques en coordonnées rectangulaires

La Terre est une sphère. Ainsi, nos coordonnées-latitude et longitude-sont des coordonnées sphériques. Nous avons donc besoin d’écrire une fonction qui convertit les coordonnées sphériques en coordonnées rectangulaires (définisse également un rayon de la Terre et fournisse la π constant):

rayon = 10 pi = 3.14159
#convertit les coordonnées sphériques en coordonnées rectangulaires def sphereToRect (r, a, b): return 'x': r * math.sin (b * pi / 180) * math.cos (a * pi / 180), 'y' : r * math.sin (b * pi / 180) * math.sin (a * pi / 180), 'z': r * math.cos (b * pi / 180)

Nous pouvons obtenir de meilleures performances en stockant certains calculs utilisés plusieurs fois:

#convertit les coordonnées sphériques en coordonnées rectangulaires def sphereToRect (r, a, b): aRad = math.radians (a) bRad = math.radians (b) r_sin_b = r * math.sin (bRad) return 'x': r_sin_b * math.cos (aRad), 'y': r_sin_b * math.sin (aRad), 'z': r * math.cos (bRad)

Nous pouvons écrire des fonctions composites qui combineront toutes les étapes précédentes en une seule fonction allant directement des coordonnées sphériques ou rectangulaires aux images en perspective:

#fonctions pour le traçage des points def rectPlot (coordonnée): retour perspectiveProjection (physicalProjection (coordonnée)) def spherePlot (coordonnée, sRadius): retour rectPlot (sphereToRect (sRadius, coordonnée ['long'], coordonnée ['lat']])

4. Rendu en SVG

Notre script doit pouvoir écrire dans un fichier SVG. Donc, il devrait commencer par:

f = open ('globe.svg', 'w') f.write ('\ n\ n ')

Et finissez avec:

f.write ('')

Produire un fichier SVG vide mais valide. Dans ce fichier, le script doit pouvoir créer des objets SVG. Nous allons donc définir deux fonctions qui lui permettront de dessiner des points et des polygones SVG:

#Draws SVG objet cercle def svgCircle (coordonnée, circleRadius, couleur): f.write ('\ n ') # Dessine un noeud de polygone SVG def polyNode (coordonnée): f.write (str (coordonnée [' x '] + 400) +', '+ str (coordonnée [' y '] + 400) + ")

Nous pouvons le tester en générant une grille sphérique de points:

GRILLE DE TRAVAIL pour x dans la plage (72): pour y dans la plage (36): svgCircle (spherePlot ('long': 5 * x, 'lat': 5 ​​* y, rayon), 1, '#ccc' )

Ce script, une fois enregistré et exécuté, devrait produire quelque chose comme ceci:


5. Transformer les données cartographiques SVG

Pour lire un fichier SVG, un script doit pouvoir lire un fichier XML, car SVG est un type de XML. C'est pourquoi nous avons importé xml.etree.ElementTree. Ce module vous permet de charger XML / SVG dans un script sous forme de liste imbriquée:

tree = ET.parse ('BlankMap equirectangular states.svg') racine = tree.getroot ()

Vous pouvez naviguer jusqu'à un objet dans le SVG à travers les index de la liste (vous devez généralement examiner le code source du fichier de carte pour en comprendre la structure). Dans notre cas, chaque pays est situé à racine [4] [0] [X] [n], où X est le numéro du pays, en commençant par 1, et n représente les divers sous-chemins qui décrivent le pays. Les contours réels du pays sont stockés dans le attribut, accessible par racine [4] [0] [X] [n] .attrib ['d'].

1. Construire des boucles

Nous ne pouvons pas simplement parcourir cette carte car elle contient au début un élément «factice» qui doit être ignoré. Il faut donc compter le nombre d’objets «pays» et en soustraire un pour éliminer le mannequin. Puis nous parcourons les objets restants.

countries = len (racine [4] [0]) - 1 pour x dans la plage (pays): racine [4] [0] [x + 1]

Certains objets pays incluent plusieurs chemins. C'est pourquoi nous parcourons chaque chemin dans chaque pays:

countries = len (racine [4] [0]) - 1 pour x dans la plage (pays): pour le chemin dans la racine [4] [0] [x + 1]:

Dans chaque chemin, il y a des contours disjoints séparés par les caractères 'Z M' dans le chaîne, donc nous avons divisé le chaîne le long de ce délimiteur et parcourir à travers ceux.

countries = len (racine [4] [0]) - 1 pour x dans la plage (pays): pour chemin dans la racine [4] [0] [x + 1]: pour k dans re.split ('Z M', path.attrib ['d']):

Nous divisons ensuite chaque contour par les délimiteurs 'Z', 'L' ou 'M' pour obtenir la coordonnée de chaque point du tracé:

pour x dans la plage (pays): pour chemin dans la racine [4] [0] [x + 1]: pour k dans re.split ('Z M', chemin.attrib ['d']): pour i dans re .split ('Z | M | L', k):

Ensuite, nous supprimons tous les caractères non numériques des coordonnées et les scindons en deux le long des virgules, en donnant les latitudes et les longitudes. Si les deux existent, nous les stockons dans un sphèreCoordonnées dictionnaire (sur la carte, les coordonnées de latitude vont de 0 à 180 °, mais nous voulons qu’elles se situent entre -90 ° et 90 ° au nord et au sud, nous soustrayons donc à 90 °).

pour x dans la plage (pays): pour chemin dans la racine [4] [0] [x + 1]: pour k dans re.split ('Z M', chemin.attrib ['d']): pour i dans re .split ('Z | M | L', k): breakup = re.split (',', re.sub ("[^ - 0123456789.,]", "", i)) si breakup [0] et breakup [1]: sphereCoordinates =  sphereCoordinates ['long'] = float (breakup [0]) sphereCoordinates ['lat'] = float (breakup [1]) - 90

Ensuite, si nous le testons en traçant quelques points (svgCircle (spherePlot (sphereCoordinates, radius), 1, '# 333')), nous obtenons quelque chose comme ceci:

2. Résoudre pour occlusion

Cela ne fait pas la distinction entre les points situés à proximité du globe et ceux situés à l'autre bout du globe. Si nous voulons simplement imprimer des points sur le côté visible de la planète, nous devons pouvoir déterminer de quel côté de la planète se trouve un point donné.. 

Nous pouvons le faire en calculant les deux points de la sphère où un rayon de la caméra au point intersecterait la sphère. Cette fonction implémente la formule pour résoudre les distances à ces deux points-dNear et dFar:

cameraDistanceSquare = sumOfSquares (camera) #distance du centre du globe à la résolution de la caméra ray, cameraForward)
def occlude (spherePoint): point = sphereToRect (rayon, spherePoint ['long'], sphèrePoint ['lat']) ray = findVector (caméra, point) d1 = magnitude (rayon) #distance de la caméra au point dot_l = np. point ([rayon ['x'] / d1, rayon ['y'] / d1, rayon ['z'] / d1], vectorToList (appareil photo)) # produit du vecteur unitaire de la caméra au point et déterminant du vecteur de la caméra = math.sqrt (abs ((dot_l) ** 2 - cameraDistanceSquare + radius ** 2)) dNear = - (dot_l) + déterminant dFar = - (dot_l) - déterminant

Si la distance réelle au point, d1, est inférieur ou égal à tous les deux de ces distances, alors le point est du côté proche de la sphère. En raison des erreurs d'arrondis, cette opération comporte un peu de marge de manœuvre:

 si d1 - 0.0000000001 <= dNear and d1 - 0.0000000001 <= dFar : return True else: return False

Utiliser cette fonction comme condition devrait limiter le rendu aux points proches du côté:

 si occlude (sphereCoordinates): svgCircle (spherePlot (sphereCoordinates, radius), 1, '# 333')

6. Rendre les pays solides

Bien sûr, les points ne sont pas de vraies formes fermées, mais remplies - ils ne donnent que l'illusion de formes fermées. Dessiner des pays réellement remplis demande un peu plus de sophistication. Tout d'abord, nous devons imprimer l'intégralité de tous les pays visibles. 

Nous pouvons le faire en créant un commutateur qui s'active chaque fois qu'un pays contient un point visible, tout en stockant temporairement les coordonnées de ce pays. Si le commutateur est activé, le pays est tracé en utilisant les coordonnées enregistrées. Nous allons aussi dessiner des polygones au lieu de points.

pour x dans la plage (pays): pour le chemin dans la racine [4] [0] [x + 1]: pour le k dans re.split ('Z M', chemin.attrib ['d']): countryIsVisible = False pays = [] pour i dans re.split ('Z | M | L', k): décomposition = re.split (',', re.sub ("[^ - 0123456789.,]", "", i) ) si rupture [0] et rupture [1]: sphereCoordinates =  sphèreCoordinates ['long'] = float (rupture [0]) sphereCoordinates ['lat'] = float (rupture [1]) - 90 #DRAW COUNTRY if occlude (sphèreCoordinates): country.append ([sphèreCoordinates, rayon]) countryIsVisible = True else: country.append ([sphèreCoordinates, rayon]) si countryIsVisible: f.write ('\ n \ n ')

C'est difficile à dire, mais les pays du bout du monde se replient sur eux-mêmes, ce que nous ne voulons pas (regardez le Brésil).

1. Tracez le disque de la terre

Pour que les pays s'affichent correctement sur les bords du globe, nous devons d'abord tracer le disque du globe avec un polygone (le disque que vous voyez à partir des points est une illusion d'optique). Le disque est délimité par le bord visible du globe, un cercle. Les opérations suivantes calculent le rayon et le centre de ce cercle, ainsi que la distance du plan contenant le cercle de la caméra et le centre du globe..

#TRACE LIMB limbRadius = math.sqrt (rayon ** 2 - rayon ** 4 / cameraDistanceSquare) cx = caméra ['x'] * rayon ** 2 / cameraDistanceSquare cy = caméra ['y'] * rayon ** 2 / cameraDistanceSquare cz = camera ['z'] * rayon ** 2 / cameraDistanceSquare planeDistance = magnitude (caméra) * (1 - rayon ** 2 / cameraDistanceSquare) planeDisplacement = math.sqrt (cx ** 2 + cy ** 2 + cz ** 2)

La terre et la caméra (point gris foncé) vus d'en haut. La ligne rose représente le bord visible de la terre. Seul le secteur ombré est visible à la caméra.

Ensuite, pour tracer un cercle dans ce plan, nous construisons deux axes parallèles à ce plan:

#trade & negate x et y pour obtenir un vecteur perpendiculaire unitVectorCamera = unitVector (caméra) aV = unitVector ('x': -unitVectorCamera ['y'], 'y': unitVectorCamera ['x'], 'z': 0) bV = np.cross (vecteurToList (aV), vecteurToList (unitVectorCamera))

Ensuite, nous allons simplement tracer graphiquement sur ces axes par incréments de 2 degrés pour tracer un cercle dans ce plan avec ce rayon et ce centre (voir cette explication pour le calcul):

pour t dans la plage (180): thêta = math.radians (2 * t) cosT = math.cos (thêta) sinT = math.sin (thêta) limbPoint = 'x': cx + limbRadius * (cosT * aV [ 'x'] + sinT * bV [0]), 'y': cy + limbRadius * (cosT * aV ['y'] + sinT * bV [1]), 'z': cz + limbRadius * (cosT * aV ['z'] + sinT * bV [2])

Ensuite, nous encapsulons simplement tout cela avec du code de dessin de polygone:

f.write ('')

Nous créons également une copie de cet objet à utiliser ultérieurement en tant que masque d'écrêtage pour tous nos pays:

f.write ('')

Cela devrait vous donner ceci:

2. Coupure sur le disque

En utilisant le disque nouvellement calculé, nous pouvons modifier notre autre déclaration dans le code de traçage du pays (pour quand les coordonnées sont sur le côté caché du globe) pour tracer ces points quelque part en dehors du disque:

 sinon: tangentscale = (radius + planeDisplacement) / (pi * 0.5) rr = 1 + abs (math.tan ((distanceToPoint (sphereCoordinates) - planeDistance) / tangentscale)) country.append ([sphereCoordinates, rayon * rr])

Cela utilise une courbe tangente pour soulever les points cachés au-dessus de la surface de la Terre, donnant l’apparence qu’ils sont répartis autour de celle-ci:

Ce n’est pas totalement mathématique (ça tombe en panne si la caméra n’est pas pointée au centre de la planète), mais c’est simple et ça marche la plupart du temps. Puis en ajoutant simplement clip-path = "url (#clipglobe)" au code de dessin du polygone, nous pouvons couper les pays au bord du globe:

 si countryIsVisible: f.write ('

J'espère que vous avez apprécié ce tutoriel! Amusez-vous avec vos globes de vecteur!