Einen SVG-Globus rendern

In diesem Tutorial zeige ich Ihnen, wie Sie eine SVG-Karte als Vektor auf einen Globus projizieren. Um die mathematischen Transformationen durchzuführen, die zum Projizieren der Karte auf eine Kugel erforderlich sind, müssen wir Python-Skripte verwenden, um die Kartendaten zu lesen und in ein Bild eines Globus zu übersetzen. In diesem Tutorial wird davon ausgegangen, dass Sie Python 3.4 ausführen, das neueste verfügbare Python.

Inkscape hat eine Art Python-API, die für eine Vielzahl von Dingen verwendet werden kann. Da wir jedoch nur an der Transformation von Formen interessiert sind, ist es einfacher, einfach ein eigenständiges Programm zu schreiben, das SVG-Dateien selbst liest und druckt.

1. Formatieren Sie die Karte

Der gewünschte Kartentyp wird als äquirektanguläre Karte bezeichnet. In einer äquirektangulären Karte entsprechen der Längen- und Breitengrad eines Ortes seiner x- und y-Position auf der Karte. Eine äquirektanguläre Weltkarte ist auf Wikimedia Commons zu finden (hier eine Version mit US-Bundesstaaten).

SVG-Koordinaten können auf verschiedene Weise definiert werden. Sie können beispielsweise relativ zum zuvor definierten Punkt oder absolut vom Ursprung aus definiert sein. Um uns das Leben zu erleichtern, wollen wir die Koordinaten in der Karte in die absolute Form umwandeln. Inkscape kann dies tun. Gehen Sie zu den Inkscape-Einstellungen (im Menü Bearbeiten) und setzen Sie unter Eingabe/Ausgabe > SVG-Ausgabe das Pfadzeichenfolgenformat auf Absolut.

Inkscape konvertiert die Koordinaten nicht automatisch; Sie müssen eine Art Transformation auf den Pfaden durchführen, um dies zu erreichen. Der einfachste Weg, dies zu tun, besteht darin, einfach alles auszuwählen und es mit jeweils einem Druck auf die Auf- und Abwärtspfeile nach oben und unten zu verschieben. Speichern Sie dann die Datei erneut.

2. Starten Sie Ihr Python-Skript

Erstellen Sie eine neue Python-Datei. Importieren Sie die folgenden Module:

1

import sys

2

import re

3

import math

4

import time

5

import datetime

6

import numpy as np

7

8

import xml.etree.ElementTree as ET

Sie müssen NumPy installieren, eine Bibliothek, mit der Sie bestimmte Vektoroperationen wie Punktprodukt und Kreuzprodukt ausführen können.

3. Die Mathematik von den perspektivischen Projektion

Das Projizieren eines Punkts im dreidimensionalen Raum in ein 2D-Bild umfasst das Finden eines Vektors von der Kamera zum Punkt und das anschließende Aufteilen dieses Vektors in drei senkrechte Vektoren.

Die beiden Teilvektoren senkrecht zum Kameravektor (der Richtung, in die die Kamera gerichtet ist) werden die x- und y-Koordinaten eines orthogonal projizierten Bildes. Der zum Kameravektor parallele Teilvektor wird als z-Abstand des Punktes bezeichnet. Um ein orthogonales Bild in ein perspektivisches Bild umzuwandeln, teilen Sie jede x- und y-Koordinate durch den z-Abstand.

An dieser Stelle ist es sinnvoll, bestimmte Kameraparameter zu definieren. Zuerst müssen wir wissen, wo sich die Kamera im 3D-Raum befindet. Speichern Sie seine x-, y- und z-Koordinaten in einem Wörterbuch.

1

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

Der Globus befindet sich am Ursprung, daher ist es sinnvoll, die Kamera darauf auszurichten. Das bedeutet, dass der Kamerarichtungsvektor das Gegenteil der Kameraposition ist.

1

cameraForward = {'x': -1*camera['x'], 'y': -1*camera['y'], 'z': -1*camera['z']}

Es reicht nicht nur zu bestimmen, in welche Richtung die Kamera zeigt – Sie müssen auch eine Drehung für die Kamera festlegen. Definieren Sie dazu einen Vektor senkrecht zum cameraForward-Vektor.

1

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

1. Definieren Sie nützliche Vektorfunktionen

Es wird sehr hilfreich sein, bestimmte Vektorfunktionen in unserem Programm definiert zu haben. Definieren Sie eine Vektorbetragsfunktion:

1

#magnitude of a 3D vector

2

def sumOfSquares(vector):

3

    return vector['x']**2 + vector['y']**2 + vector['z']**2

4

def magnitude(vector):

5

	return math.sqrt(sumOfSquares(vector))

Wir müssen in der Lage sein, einen Vektor auf einen anderen zu projizieren. Da es sich bei dieser Operation um ein Punktprodukt handelt, ist es viel einfacher, die NumPy-Bibliothek zu verwenden. NumPy nimmt Vektoren jedoch in Listenform an, ohne die expliziten Bezeichner ‚x‘, ‚y‘, ‚z‘, daher benötigen wir eine Funktion, um unsere Vektoren in NumPy-Vektoren umzuwandeln.

1

#converts dictionary vector to list vector

2

def vectorToList (vector):

3

    return [vector['x'], vector['y'], vector['z']]

1

#projects u onto v

2

def vectorProject(u, v):

3

    return np.dot(vectorToList (v), vectorToList (u))/magnitude(v)

Es ist schön, eine Funktion zu haben, die uns einen Einheitsvektor in Richtung eines gegebenen Vektors liefert:

1

#get unit vector

2

def unitVector(vector):

3

    magVector = magnitude(vector)

4

	return {'x': vector['x']/magVector, 'y': vector['y']/magVector, 'z': vector['z']/magVector }

Schließlich müssen wir in der Lage sein, zwei Punkte zu nehmen und einen Vektor zwischen ihnen zu finden:

1

#Calculates vector from two points, dictionary form

2

def findVector (origin, point):    

3

	return { 'x': point['x'] - origin['x'], 'y': point['y'] - origin['y'], 'z': point['z'] - origin['z'] }

2. Kameraachsen definieren

Jetzt müssen wir nur noch die Definition der Kameraachsen beenden. Wir haben bereits zwei dieser Achsen – cameraForward und cameraPerpendicular, die dem z-Abstand und der x-Koordinate des Kamerabilds entsprechen.

Jetzt brauchen wir nur noch die dritte Achse, definiert durch einen Vektor, der die y-Koordinate des Kamerabildes darstellt. Wir können diese dritte Achse finden, indem wir das Kreuzprodukt dieser beiden Vektoren bilden, indem wir NumPy verwenden – np.cross(vectorToList(cameraForward), vectorToList(cameraPerpendicular)).

Das erste Element im Ergebnis entspricht der x-Komponente; die zweite zur y-Komponente und die dritte zur z-Komponente, so dass der erzeugte Vektor gegeben ist durch:

1

#Calculates horizon plane vector (points upward)

2

cameraHorizon = {'x': np.cross(vectorToList(cameraForward) , vectorToList(cameraPerpendicular))[0], 'y': np.cross(vectorToList(cameraForward) , vectorToList(cameraPerpendicular))[1], 'z': np.cross(vectorToList(cameraForward) , vectorToList(cameraPerpendicular))[2] }

3. Auf Orthogonal projizieren

Um den orthogonalen x-, y- und z-Abstand zu ermitteln, suchen wir zunächst den Vektor, der die Kamera und den fraglichen Punkt verbindet, und projizieren ihn dann auf jede der drei zuvor definierten Kameraachsen:

1

def physicalProjection (point):

2

    pointVector = findVector(camera, point)

3

		#pointVector is a vector starting from the camera and ending at a point in question

4

	return {'x': vectorProject(pointVector, cameraPerpendicular), 'y': vectorProject(pointVector, cameraHorizon), 'z': vectorProject(pointVector, cameraForward)}

Ein Punkt (dunkelgrau) wird auf die drei Kameraachsen (grau) projiziert. x ist rot, y ist grün und z ist blau.

4. Perspektivisch projizieren

Die perspektivische Projektion nimmt einfach das x und y der orthogonalen Projektion und teilt jede Koordinate durch den z-Abstand. Dies macht es so, dass weiter entfernte Dinge kleiner aussehen als Dinge, die näher an der Kamera sind.

Da eine Division durch z sehr kleine Koordinaten ergibt, multiplizieren wir jede Koordinate mit einem Wert, der der Brennweite der Kamera entspricht.

1

focalLength = 1000

1

# draws points onto camera sensor using xDistance, yDistance, and zDistance

2

def perspectiveProjection (pCoords):

3

    scaleFactor = focalLength/pCoords['z']

4

	return {'x': pCoords['x']*scaleFactor, 'y': pCoords['y']*scaleFactor}

5. Konvertieren Sie sphärische Koordinaten in rechteckige Koordinaten

Die Erde ist eine Kugel. Somit sind unsere Koordinaten – Breiten- und Längengrad – Kugelkoordinaten. Wir müssen also eine Funktion schreiben, die Kugelkoordinaten in rechtwinklige Koordinaten umwandelt (sowie einen Erdradius definieren und die π-Konstante angeben):

1

radius = 10

2

pi = 3.14159

1

#converts spherical coordinates to rectangular coordinates

2

def sphereToRect (r, a, b):

3

	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) }

Wir können eine bessere Leistung erzielen, indem wir einige Berechnungen speichern, die mehr als einmal verwendet werden:

1

#converts spherical coordinates to rectangular coordinates

2

def sphereToRect (r, a, b):

3

    aRad = math.radians(a)

4

	bRad = math.radians(b)

5

	r_sin_b = r*math.sin(bRad)

6

	return {'x': r_sin_b*math.cos(aRad), 'y': r_sin_b*math.sin(aRad), 'z': r*math.cos(bRad) }

Wir können einige zusammengesetzte Funktionen schreiben, die alle vorherigen Schritte in einer Funktion kombinieren – direkt von sphärischen oder rechteckigen Koordinaten zu perspektivischen Bildern:

1

#functions for plotting points

2

def rectPlot (coordinate):

3

    return perspectiveProjection(physicalProjection(coordinate))

4

def spherePlot (coordinate, sRadius):

5

	return rectPlot(sphereToRect(sRadius, coordinate['long'], coordinate['lat']))

4. Rendern in SVG

Unser Skript muss in der Lage sein, in eine SVG-Datei zu schreiben. Es sollte also beginnen mit:

1

f = open('globe.svg', 'w')

2

f.write('<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n<svg viewBox="0 0 800 800" version="1.1"\nxmlns="http://www.w3.org/2000/svg" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n')

Und enden Sie mit:

1

f.write('</svg>')

Erstellen einer leeren, aber gültigen SVG-Datei. In dieser Datei muss das Skript SVG-Objekte erstellen können, daher werden wir zwei Funktionen definieren, die es ermöglichen, SVG-Punkte und -Polygone zu zeichnen:

1

#Draws SVG circle object

2

def svgCircle (coordinate, circleRadius, color):

3

    f.write('<circle cx=\"' + str(coordinate['x'] + 400) + '\" cy=\"' + str(coordinate['y'] + 400) + '\" r=\"' + str(circleRadius) + '\" style=\"fill:' + color + ';\"/>\n')

4

#Draws SVG polygon node

5

def polyNode (coordinate):

6

	f.write(str(coordinate['x'] + 400) + ',' + str(coordinate['y'] + 400) + ' ')

Wir können dies testen, indem wir ein kugelförmiges Punktgitter rendern:

1

#DRAW GRID

2

for x in range(72):

3

    for y in range(36):

4

		svgCircle (spherePlot( { 'long': 5*x, 'lat': 5*y }, radius ), 1, '#ccc')

Dieses Skript sollte, wenn es gespeichert und ausgeführt wird, in etwa so aussehen:

5. Transformieren Sie die SVG-Kartendaten

Um eine SVG-Datei zu lesen, muss ein Skript eine XML-Datei lesen können, da SVG ein XML-Typ ist. Deshalb haben wir xml.etree.ElementTree importiert. Mit diesem Modul können Sie die XML/SVG als verschachtelte Liste in ein Skript laden:

1

tree = ET.parse('BlankMap Equirectangular states.svg')

2

root = tree.getroot()

Über die Listenindizes können Sie zu einem Objekt im SVG navigieren (normalerweise müssen Sie sich den Quellcode der Kartendatei ansehen, um deren Struktur zu verstehen). In unserem Fall befindet sich jedes Land unter root[4][0][x][n], wobei x die Nummer des Landes ist, beginnend mit 1 und n die verschiedenen Unterpfade darstellt, die das Land umreißen. Die tatsächlichen Konturen des Landes werden im d-Attribut gespeichert, auf das über root[4][0][x][n].attrib['d'] zugegriffen werden kann.

1. Konstruiere Schleifen

Wir können diese Karte nicht einfach durchlaufen, da sie am Anfang ein „Dummy“-Element enthält, das übersprungen werden muss. Wir müssen also die Anzahl der „Land“-Objekte zählen und eins subtrahieren, um den Dummy loszuwerden. Dann durchlaufen wir die restlichen Objekte.

1

countries = len(root[4][0]) - 1

2

3

for x in range(countries):

4

    root[4][0][x + 1]

Einige Länderobjekte enthalten mehrere Pfade, weshalb wir dann jeden Pfad in jedem Land durchlaufen:

1

countries = len(root[4][0]) - 1

2

3

for x in range(countries):

4

    for path in root[4][0][x + 1]:

Innerhalb jedes Pfads gibt es disjunkte Konturen, die durch die Zeichen „Z M“ in der d-Zeichenfolge getrennt sind, also teilen wir die d-Zeichenfolge entlang dieses Trennzeichens auf und durchlaufen diese.

1

countries = len(root[4][0]) - 1

2

3

for x in range(countries):

4

    for path in root[4][0][x + 1]:

5

		for k in re.split('Z M', path.attrib['d']):

Wir teilen dann jede Kontur durch die Trennzeichen „Z“, „L“ oder „M“, um die Koordinaten jedes Punktes im Pfad zu erhalten:

1

for x in range(countries):

2

    for path in root[4][0][x + 1]:

3

		for k in re.split('Z M', path.attrib['d']):

4

			for i in re.split('Z|M|L', k):

Dann entfernen wir alle nicht numerischen Zeichen aus den Koordinaten und teilen sie entlang der Kommas in zwei Hälften, um die Breiten- und Längengrade zu erhalten. Wenn beide vorhanden sind, speichern wir sie in einem sphereCoordinates-Wörterbuch (in der Karte gehen die Breitengradkoordinaten von 0 bis 180°, aber wir möchten, dass sie von –90° bis 90° – Nord und Süd – gehen, also subtrahieren wir 90°).

1

for x in range(countries):

2

    for path in root[4][0][x + 1]:

3

		for k in re.split('Z M', path.attrib['d']):

4

			for i in re.split('Z|M|L', k):

5

				breakup = re.split(',', re.sub("[^-0123456789.,]", "", i))

6

				if breakup[0] and breakup[1]:

7

					sphereCoordinates = {}

8

					sphereCoordinates['long'] = float(breakup[0])

9

					sphereCoordinates['lat'] = float(breakup[1]) - 90

Wenn wir es dann testen, indem wir einige Punkte zeichnen (svgCircle(spherePlot(sphereCoordinates, radius), 1, '#333')), erhalten wir etwa Folgendes:

2. Auflösen für Okklusion

Dabei wird nicht zwischen Punkten auf der nahen Seite des Globus und Punkten auf der anderen Seite des Globus unterschieden. Wenn wir nur Punkte auf die sichtbare Seite des Planeten drucken möchten, müssen wir in der Lage sein, herauszufinden, auf welcher Seite des Planeten sich ein bestimmter Punkt befindet.

Wir können dies tun, indem wir die beiden Punkte auf der Kugel berechnen, an denen sich ein Strahl von der Kamera zu dem Punkt mit der Kugel schneiden würde. Diese Funktion implementiert die Formel zum Lösen der Entfernungen zu diesen beiden Punkten – dNear und dFar:

1

cameraDistanceSquare = sumOfSquares(camera)

2

    	#distance from globe center to camera

3

4

def distanceToPoint(spherePoint):

5

	point = sphereToRect(radius, spherePoint['long'], spherePoint['lat'])

6

	ray = findVector(camera,point)

7

	return vectorProject(ray, cameraForward)

1

def occlude(spherePoint):

2

    point = sphereToRect(radius, spherePoint['long'], spherePoint['lat'])

3

	ray = findVector(camera,point)

4

	d1 = magnitude(ray)

5

		#distance from camera to point

6

7

	dot_l = np.dot( [ray['x']/d1, ray['y']/d1, ray['z']/d1], vectorToList(camera) )

8

		#dot product of unit vector from camera to point and camera vector

9

10

	determinant = math.sqrt(abs( (dot_l)**2 - cameraDistanceSquare + radius**2 ))

11

	dNear = -(dot_l) + determinant

12

	dFar = -(dot_l) - determinant

Wenn die tatsächliche Entfernung zum Punkt d1 kleiner oder gleich diesen beiden Entfernungen ist, befindet sich der Punkt auf der nahen Seite der Kugel. Aufgrund von Rundungsfehlern ist in diese Operation ein wenig Spielraum eingebaut:

1

    if d1 - 0.0000000001 <= dNear and d1 - 0.0000000001 <= dFar :

2

		return True

3

	else:

4

		return False

Die Verwendung dieser Funktion als Bedingung sollte das Rendering auf Punkte auf der nahen Seite beschränken:

1

    				if occlude(sphereCoordinates):

2

						svgCircle(spherePlot(sphereCoordinates, radius), 1, '#333')

6. Solide Länder rendern

Natürlich sind die Punkte keine echten geschlossenen, gefüllten Formen – sie vermitteln nur die Illusion geschlossener Formen. Das Zeichnen tatsächlich gefüllter Länder erfordert etwas mehr Raffinesse. Zunächst müssen wir die Gesamtheit aller sichtbaren Länder drucken.

Wir können dies tun, indem wir einen Schalter erstellen, der jedes Mal aktiviert wird, wenn ein Land einen sichtbaren Punkt enthält, währenddessen die Koordinaten dieses Landes vorübergehend gespeichert werden. Wenn der Schalter aktiviert ist, wird das Land anhand der gespeicherten Koordinaten gezeichnet. Wir werden auch Polygone anstelle von Punkten zeichnen.

1

for x in range(countries):

2

    for path in root[4][0][x + 1]:

3

    	for k in re.split('Z M', path.attrib['d']):

4

5

			countryIsVisible = False

6

			country = []

7

			for i in re.split('Z|M|L', k):

8

9

				breakup = re.split(',', re.sub("[^-0123456789.,]", "", i))

10

				if breakup[0] and breakup[1]:

11

					sphereCoordinates = {}

12

					sphereCoordinates['long'] = float(breakup[0])

13

					sphereCoordinates['lat'] = float(breakup[1]) - 90

14

15

					#DRAW COUNTRY

16

17

					if occlude(sphereCoordinates):

18

						country.append([sphereCoordinates, radius])

19

20

						countryIsVisible = True

21

22

					else:

23

						country.append([sphereCoordinates, radius])

24

25

			if countryIsVisible:

26

				f.write('<polygon points=\"')

27

				for i in country:

28

					polyNode(spherePlot(i[0], i[1]))

29

				f.write('\" style="fill:#ff3092;stroke: #fff;stroke-width:0.3\" />\n\n')

Es ist schwer zu sagen, aber die Länder am Rande der Welt falten sich in sich zusammen, was wir nicht wollen (siehe Brasilien).

1. Verfolgen Sie die Erdscheibe

Damit die Länder an den Rändern des Globus richtig gerendert werden, müssen wir zuerst die Scheibe des Globus mit einem Polygon nachzeichnen (die Scheibe, die Sie an den Punkten sehen, ist eine optische Täuschung). Die Scheibe wird vom sichtbaren Rand der Weltkugel – einem Kreis – umrissen. Die folgenden Operationen berechnen den Radius und den Mittelpunkt dieses Kreises sowie den Abstand der Ebene, die den Kreis enthält, von der Kamera und dem Mittelpunkt des Globus.

1

#TRACE LIMB

2

limbRadius = math.sqrt( radius**2 - radius**4/cameraDistanceSquare )

3

4

cx = camera['x']*radius**2/cameraDistanceSquare

5

cy = camera['y']*radius**2/cameraDistanceSquare

6

cz = camera['z']*radius**2/cameraDistanceSquare

7

8

planeDistance = magnitude(camera)*(1 - radius**2/cameraDistanceSquare)

9

planeDisplacement = math.sqrt(cx**2 + cy**2 + cz**2)

Die Erde und die Kamera (dunkelgrauer Punkt) von oben gesehen. Die rosa Linie stellt den sichtbaren Rand der Erde dar. Nur der schattierte Sektor ist für die Kamera sichtbar.

Um dann einen Kreis in dieser Ebene zu zeichnen, konstruieren wir zwei Achsen parallel zu dieser Ebene:

1

#trade & negate x and y to get a perpendicular vector

2

unitVectorCamera = unitVector(camera)

3

aV = unitVector( {'x': -unitVectorCamera['y'], 'y': unitVectorCamera['x'], 'z': 0} )

4

bV = np.cross(vectorToList(aV), vectorToList( unitVectorCamera ))

Dann zeichnen wir diese Achsen einfach in Schritten von 2 Grad auf, um einen Kreis in dieser Ebene mit diesem Radius und Mittelpunkt zu zeichnen (siehe diese Erklärung für die Mathematik):

1

for t in range(180):

2

    theta = math.radians(2*t)

3

    cosT = math.cos(theta)

4

    sinT = math.sin(theta)

5

6

    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]) }

Dann kapseln wir das alles einfach mit Polygon-Zeichnungscode:

1

f.write('<polygon id=\"globe\" points=\"')

2

for t in range(180):

3

    theta = math.radians(2*t)

4

	cosT = math.cos(theta)

5

	sinT = math.sin(theta)

6

7

	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]) }

8

9

	polyNode(rectPlot(limbPoint))

10

11

f.write('\" style="fill:#eee;stroke: none;stroke-width:0.5\" />')

Wir erstellen auch eine Kopie dieses Objekts, um sie später als Schnittmaske für alle unsere Länder zu verwenden:

1

f.write('<clipPath id=\"clipglobe\"><use xlink:href=\"#globe\"/></clipPath>')

Das sollte dir das geben:

2. Clipping auf die Festplatte

Mit der neu berechneten Scheibe können wir unsere else-Anweisung im Länder-Plot-Code ändern (für Koordinaten auf der verborgenen Seite des Globus), um diese Punkte irgendwo außerhalb der Scheibe zu zeichnen:

1

    				else:

2

						tangentscale = (radius + planeDisplacement)/(pi*0.5)

3

						rr = 1 + abs(math.tan( (distanceToPoint(sphereCoordinates) - planeDistance)/tangentscale ))

4

						country.append([sphereCoordinates, radius*rr])

Dies verwendet eine Tangentenkurve, um die verborgenen Punkte über die Erdoberfläche zu heben, was den Anschein erweckt, dass sie um sie herum verteilt sind:

Das ist mathematisch nicht ganz schlüssig (es bricht zusammen, wenn die Kamera nicht grob auf das Zentrum des Planeten ausgerichtet ist), aber es ist einfach und funktioniert die meiste Zeit. Durch einfaches Hinzufügen von clip-path="url(#clipglobe)" zum Polygon-Zeichnungscode können wir die Länder sauber an den Rand des Globus beschneiden:

1

    		if countryIsVisible:

2

				f.write('<polygon clip-path="url(#clipglobe)" points=\"')

Ich hoffe, Ihnen hat dieses Tutorial gefallen! Viel Spaß mit Ihren Vektorgloben!