Advertisement
  1. Code
  2. Python

Renderizar un SVG Globe

Scroll to top
Read Time: 15 min

Spanish (Español) translation by Manuel (you can also view the original English article)

Final product imageFinal product imageFinal product image
What You'll Be Creating

En este tutorial, voy a mostrar cómo tomar un mapa SVG y proyectarlo en un globo terráqueo, como un vector. Para llevar a cabo las transformaciones matemáticas necesarias para proyectar el mapa en una esfera, debemos utilizar scripts de Python para leer los datos del mapa y traducirlos en una imagen de un globo terráqueo. Este tutorial asume que estás ejecutando Python 3.4, la última versión disponible de Python.

Inkscape tiene una especie de API de Python que se puede utilizar para hacer una variedad de cosas. Sin embargo, ya que sólo estamos interesados en la transformación de las formas, es más fácil escribir un programa independiente que lee e imprime los archivos SVG por su cuenta.

1. Formatear el mapa

El tipo de mapa que queremos se llama mapa equirectangular. En un mapa equirectangular, la longitud y latitud de un lugar corresponden a su posición x e y en el mapa. En Wikimedia Commons se puede encontrar un mapamundi equirectangular (aquí hay una versión con los estados de Estados Unidos).

Las coordenadas SVG pueden definirse de varias maneras. Por ejemplo, pueden ser relativas al punto previamente definido, o definidas absolutamente desde el origen. Para facilitarnos la vida, queremos convertir las coordenadas del mapa a la forma absoluta. Inkscape puede hacer esto. Ve a las preferencias de Inkscape (en el menú Editar) y en Entrada/Salida > Salida SVG, establece el formato de cadena de ruta como Absoluto.

Inkscape preferencesInkscape preferencesInkscape preferences

Inkscape no convertirá automáticamente las coordenadas; tiene que realizar algún tipo de transformación en los trazados para que eso ocurra. La forma más fácil de hacerlo es simplemente seleccionar todo y moverlo hacia arriba y hacia abajo con una pulsación de cada una de las flechas hacia arriba y hacia abajo. Luego vuelve a guardar el archivo.

2. Inicia tu script de Python

Crea un nuevo archivo de Python. Importa los siguientes módulos:

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

Necesitarás instalar NumPy, una biblioteca que te permite realizar ciertas operaciones vectoriales como el producto punto y el producto cruz.

3. Las matemáticas de la proyección en perspectiva

Proyectar un punto en un espacio tridimensional en una imagen 2D implica encontrar un vector desde la cámara hasta el punto, y luego dividir ese vector en tres vectores perpendiculares.

Los dos vectores parciales perpendiculares al vector de la cámara (la dirección en la que está orientada la cámara) se convierten en las coordenadas x e y de una imagen proyectada ortogonalmente. El vector parcial paralelo al vector de la cámara se convierte en algo llamado distancia z del punto. Para convertir una imagen ortogonal en una imagen en perspectiva, divide cada coordenada x e y por la distancia z.

En este punto, tiene sentido definir ciertos parámetros de la cámara. En primer lugar, necesitamos saber dónde se encuentra la cámara en el espacio 3D. Guarda sus coordenadas x, y y z en un diccionario.

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

El globo terráqueo estará situado en el origen, por lo que tiene sentido orientar la cámara hacia él. Eso significa que el vector de dirección de la cámara será el opuesto a la posición de la cámara.

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

No basta con determinar la dirección en la que está orientada la cámara, sino que también hay que fijar una rotación para la cámara. Para ello, define un vector perpendicular al vector cameraForward.

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

1. Definir funciones vectoriales útiles

Será muy útil tener definidas ciertas funciones vectoriales en nuestro programa. Definir una función de magnitud vectorial:

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

Necesitamos poder proyectar un vector sobre otro. Como esta operación implica un producto punto, es mucho más fácil utilizar la biblioteca NumPy. NumPy, sin embargo, toma los vectores en forma de lista, sin los identificadores explícitos 'x', 'y', 'z', por lo que necesitamos una función para convertir nuestros vectores en vectores NumPy.

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 bueno tener una función que nos dé un vector unitario en la dirección de un vector dado:

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 }

Por último, tenemos que ser capaces de tomar dos puntos y encontrar un vector entre ellos:

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. Definir los ejes de la cámara

Ahora sólo tenemos que terminar de definir los ejes de la cámara. Ya tenemos dos de estos ejes, cameraForward y cameraPerpendicular, correspondientes a la distancia z y a la coordenada x de la imagen de la cámara.

Ahora sólo necesitamos el tercer eje, definido por un vector que represente la coordenada y de la imagen de la cámara. Podemos encontrar este tercer eje tomando el producto cruzado de esos dos vectores, usando NumPy, np.cross(vectorToList(cameraForward), vectorToList(cameraPerpendicular)).

El primer elemento del resultado corresponde a la componente x; el segundo a la componente y, y el tercero a la componente z, por lo que el vector producido viene dado por:

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. Proyectar a Ortogonal

Para encontrar la distancia ortogonal x, y y z, primero encontramos el vector que une la cámara y el punto en cuestión, y luego lo proyectamos sobre cada uno de los tres ejes de la cámara definidos anteriormente:

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)}
A point being projected onto the three camera axesA point being projected onto the three camera axesA point being projected onto the three camera axes

Un punto (gris oscuro) que se proyecta sobre los tres ejes de la cámara (gris). x es rojo, y es verde y z es azul.

4. Proyecto de perspectiva

La proyección en perspectiva simplemente toma la x y la y de la proyección ortogonal, y divide cada coordenada por la distancia z. Esto hace que las cosas que están más lejos se vean más pequeñas que las que están más cerca de la cámara.

Como al dividir por z se obtienen coordenadas muy pequeñas, multiplicamos cada coordenada por un valor correspondiente a la distancia focal de la cámara.

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. Convertir coordenadas esféricas en coordenadas rectangulares

La Tierra es una esfera. Por tanto, nuestras coordenadas, latitud y longitud, son coordenadas esféricas. Así que tenemos que escribir una función que convierta las coordenadas esféricas en coordenadas rectangulares (además de definir un radio de la Tierra y proporcionar la constante π):

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

Podemos conseguir un mejor rendimiento almacenando algunos cálculos utilizados más de una vez:

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

Podemos escribir algunas funciones compuestas que combinen todos los pasos anteriores en una sola función, pasando directamente de las coordenadas esféricas o rectangulares a las imágenes en perspectiva:

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. Renderización a SVG

Nuestro script tiene que ser capaz de escribir en un archivo SVG. Así que debe comenzar con:

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

Y terminar con:

1
f.write('</svg>')

Produciendo un archivo SVG vacío pero válido. Dentro de ese archivo el script tiene que ser capaz de crear objetos SVG, por lo que definiremos dos funciones que le permitirán dibujar puntos y polígonos SVG:

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) + ' ')

Podemos probarlo representando una cuadrícula esférica de puntos:

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

Este script, cuando se guarda y se ejecuta, debería producir algo como esto:

A dot sphere rendered with perspectiveA dot sphere rendered with perspectiveA dot sphere rendered with perspective

5. Transformar los datos del mapa SVG

Para leer un archivo SVG, un script debe ser capaz de leer un archivo XML, ya que SVG es un tipo de XML. Por eso importamos xml.etree.ElementTree. Este módulo permite cargar el XML/SVG en un script como una lista anidada:

1
tree = ET.parse('BlankMap Equirectangular states.svg')
2
root = tree.getroot()

Se puede navegar hasta un objeto en el SVG a través de los índices de la lista (normalmente hay que echar un vistazo al código fuente del archivo del mapa para entender su estructura). En nuestro caso, cada país se encuentra en root[4][0][x][n], donde x es el número del país, empezando por el 1, y n representa los distintos subcaminos que perfilan el país. Los contornos reales del país se almacenan en el atributo d, accesible a través de root[4][0][x][n].attrib['d'].

1. Construir bucles

No podemos simplemente iterar a través de este mapa porque contiene un elemento "ficticio" al principio que debe ser saltado. Así que tenemos que contar el número de objetos "país" y restar uno para deshacernos del ficticio. A continuación, hacemos un bucle con los objetos restantes.

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

2


3
for x in range(countries):

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

Algunos objetos de países incluyen múltiples caminos, por lo que a continuación iteramos por cada camino en cada país:

1
countries = len(root[4][0]) - 1
2
3
for x in range(countries):
4
    for path in root[4][0][x + 1]:

Dentro de cada camino, hay contornos disjuntos separados por los caracteres 'Z M' en la cadena d, así que dividimos la cadena d a lo largo de ese delimitador e iteramos a través de ellos.

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']):

A continuación, dividimos cada contorno por los delimitadores "Z", "L" o "M" para obtener la coordenada de cada punto de la trayectoria:

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

A continuación, eliminamos todos los caracteres no numéricos de las coordenadas y las dividimos por la mitad a lo largo de las comas, obteniendo las latitudes y las longitudes. Si existen ambas, las almacenamos en un diccionario sphereCoordinates (en el mapa, las coordenadas de latitud van de 0 a 180°, pero nosotros queremos que vayan de -90° a 90°, al norte y al sur, así que restamos 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

Entonces, si lo probamos trazando algunos puntos (svgCircle(spherePlot(sphereCoordinates, radius), 1, '#333')), obtenemos algo así:

A dot rendering of national and state bordersA dot rendering of national and state bordersA dot rendering of national and state borders

2. Resolver la oclusión

Esto no distingue entre puntos en el lado cercano del planeta y puntos en el lado lejano del planeta. Si solo queremos imprimir puntos en el lado visible del planeta, tenemos que ser capaces de averiguar en qué lado del planeta se encuentra un punto determinado.

Podemos hacerlo calculando los dos puntos de la esfera en los que un rayo desde la cámara hasta el punto se cruzaría con la esfera. Esta función implementa la fórmula para resolver las distancias a esos dos puntos, dNear y 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

Si la distancia real al punto, d1, es menor o igual que estas dos distancias, entonces el punto está en el lado cercano de la esfera. Debido a los errores de redondeo, en esta operación hay un pequeño margen de maniobra:

1
    if d1 - 0.0000000001 <= dNear and d1 - 0.0000000001 <= dFar :
2
		return True
3
	else:
4
		return False

El uso de esta función como condición debería restringir el renderizado a los puntos cercanos:

1
    				if occlude(sphereCoordinates):
2
						svgCircle(spherePlot(sphereCoordinates, radius), 1, '#333')
Dot globe only displaying points on the near side of the planetDot globe only displaying points on the near side of the planetDot globe only displaying points on the near side of the planet

6. Renderizar países sólidos

Por supuesto, los puntos no son verdaderas formas cerradas y rellenas, solo dan la ilusión de formas cerradas. Dibujar los países rellenos reales requiere un poco más de sofisticación. En primer lugar, tenemos que imprimir la totalidad de los países visibles.

Podemos hacerlo creando un interruptor que se active cada vez que un país contenga un punto visible, almacenando temporalmente las coordenadas de ese país. Si se activa el interruptor, se dibuja el país, utilizando las coordenadas almacenadas. También dibujaremos polígonos en lugar de puntos.

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')
Solid rendering of the entirety of all visible countriesSolid rendering of the entirety of all visible countriesSolid rendering of the entirety of all visible countries

Es difícil decirlo, pero los países del borde del globo se repliegan sobre sí mismos, cosa que no queremos (miren a Brasil).

1. Trazar el disco de la Tierra

Para que los países se representen correctamente en los bordes del globo, primero tenemos que trazar el disco del globo con un polígono (el disco que se ve por los puntos es una ilusión óptica). El disco está delimitado por el borde visible del globo terráqueo: un círculo. Las siguientes operaciones calculan el radio y el centro de este círculo, así como la distancia del plano que contiene el círculo desde la cámara, y el centro del globo.

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)
2D analogy of finding the edge of the visible disk2D analogy of finding the edge of the visible disk2D analogy of finding the edge of the visible disk

La Tierra y la cámara (punto gris oscuro) vistas desde arriba. La línea rosa representa el borde visible de la tierra. Solo el sector sombreado es visible para la cámara.

Entonces, para graficar una circunferencia en ese plano, construimos dos ejes paralelos a ese plano:

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

Entonces solo tenemos que graficar en esos ejes por incrementos de 2 grados para trazar un círculo en ese plano con ese radio y centro (ver esta explicación para las matemáticas):

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

Entonces encapsulamos todo eso con código de dibujo de polígonos:

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\" />')

También creamos una copia de ese objeto para usarla después como máscara de recorte para todos nuestros países:

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

Eso debería darte esto:

Rendering of the visible diskRendering of the visible diskRendering of the visible disk

2. Recorte en el disco

Usando el disco recién calculado, podemos modificar nuestra sentencia else en el código de trazado de países (para cuando las coordenadas están en el lado oculto del globo) para trazar esos puntos en algún lugar fuera del disco:

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

Utiliza una curva tangente para elevar los puntos ocultos por encima de la superficie de la Tierra, dando la apariencia de que están repartidos alrededor de ella:

Lifted portions of countries on the far side of the globeLifted portions of countries on the far side of the globeLifted portions of countries on the far side of the globe

Esto no es del todo correcto desde el punto de vista matemático (se rompe si la cámara no apunta aproximadamente al centro del planeta), pero es sencillo y funciona la mayoría de las veces. Entonces, simplemente añadiendo clip-path="url(#clipglobe)" al código de dibujo del polígono, podemos recortar limpiamente los países al borde del globo:

1
    		if countryIsVisible:
2
				f.write('<polygon clip-path="url(#clipglobe)" points=\"')
Final product with all countries clipped to the visible diskFinal product with all countries clipped to the visible diskFinal product with all countries clipped to the visible disk

Espero que hayas disfrutado de este tutorial. ¡Diviértete con tus globos vectoriales!

Another globe renderingAnother globe renderingAnother globe rendering
Advertisement
Did you find this post useful?
Want a weekly email summary?
Subscribe below and we’ll send you a weekly email summary of all new Code tutorials. Never miss out on learning about the next big thing.
Advertisement
Looking for something to help kick start your next project?
Envato Market has a range of items for sale to help get you started.