A veces, es muy útil poder hacer selecciones, geoprocesamientos y visualizaciones sencillos de datos vectoriales. En vez de usar algún producto al uso, se puede usar Python con OGR y Matplotlib.

OGR

OGR es una librería que permite abrir y trabajar con ficheros vectoriales en bastantes formatos GIS. Entre ellos está el Shapefile de ESRI, pero también hay otros. Entre otras cosas, integra geoprocesamientos, reproyección de datos y trabajo con los atributos. Desafortunadamente, está bastante poco documentado, y la interfaz de Python, aun menos. Sin entrar en mucho más detalle, listo aquí algunos de los procesados que se pueden hacer:

  • Intersecciones
  • Equivalencias
  • Test para ver si dos features se tocan
  • Solapes
  • Distancia
  • Convex Hull
  • Buffer
  • Unión
  • Intersección
  • Diferencia...

La manera habitual de tratar datos es la siguiente:


import ogr
g = ogr.Open( fname) #Se abre el fichero
L = g.GetLayer(0) #Cogemos el layer 1
feat = L.GetNextFeature() # Primer feature
field_count = L.GetLayerDefn().GetFieldCount() # ¿Cuántos campos?
for i in range(field_count): 
 print feat.GetFieldDefnRef().GetName(), ";", # Títulos de los campos

print
while feat is not None: #Por cada feature...
 for i in range (field_count):
 print feat.GetFieldAsString(i),";", #Valores de los atributos
 print
 geom = feat.GetGeometryRef() #Es la geometría asociada al feature...
 geom.XXXXXX #Geoprocesamientos: Union, Intersect, Touches....

geo = feat.GetGeometryRef()
 Este código nos permite abrir un fichero, seleccionar una capa, y recorrerla feature a feature, imprimiendo los atributos a la pantalla. También guardamos la geometría de cada uno de estos features, para usarlo (añadiendo un buffer por ejemplo).
Está claro que es bastante fácil :)
Matplotlib
Matplotlib es una potente librería gráfica, muy parecida en su concepción a los comandos gráficos de Matlab, con los que cualquiera puede hacer unas gráficas científicas más o menos apañadas.
Cutre-Cartografía
Dícese de la cartografía que hago yo :) En el siguiente código:
  • Se lee un shapefile (sí, world.shp).
  • Se itera por sus features.
  • Se extrae cada vértice. 
  • Se dibuja cada una de ellas con Matplotlib.
Conviene saber que para OGR cada feature poligonal tiene asociadas varias geometrías (polígonos), por lo que hay que hacer iteraciones dentro de la geometría para acceder a los puntos que caracterizan cada vértice.

import ogr
import pylab

g = ogr.Open ("/tmp/mundo.shp")
L = g.GetLayer(0)
N = 0
feat = L.GetNextFeature()
while feat is not None:
	field_count = L.GetLayerDefn().GetFieldCount()
	geo = feat.GetGeometryRef()
	if geo.GetGeometryCount()<2:
		g1 = geo.GetGeometryRef( 0 )
		x =[g1.GetX(i) for i in range(g1.GetPointCount()) ]
		y =[g1.GetY(i) for i in range(g1.GetPointCount()) ]
		pylab.plot(x,y,'-',hold=True)
	for cuenta in range( geo.GetGeometryCount()):
		anillo = geo.GetGeometryRef ( cuenta )
		for cnt in range( anillo.GetGeometryCount()):
			g1 = anillo.GetGeometryRef( cnt )
			x =[g1.GetX(i) for i in range(g1.GetPointCount()) ]
			y =[g1.GetY(i) for i in range(g1.GetPointCount()) ]
			pylab.plot(x,y,'-',hold=True)

	feat = L.GetNextFeature()
pylab.xlabel("Longitud")
pylab.ylabel("Latitud")
pylab.grid(True)
pylab.title("Paises del mundo")
pylab.show()
La imagen de salida es esta:
Usando PostGIS
PostGIS son las extensiones espaciales para la base de datos PostgreSQL (en  otra página  se ha discutido cómo instalar PostGIS en Kubuntu). Permite almacenar datos que generalmente se tienen en formato Shapefile. Además de poder hacer consultas usando los campos alfanuméricos de la base de datos (usando el lenguaje SQL), permite hacer consultas espaciales (así que se podría hacer una consulta similar a esta: "dáme todos los registros que tengan el campo 'población' mayor a 40 millones y que estén a menos de 5000 km de Kazakhastan"). Las funciones de consulta de PostGIS son bastante completas, y podemos encontrar uniones geométricas, intersecciones, buffer, etc.
En el presente ejemplo, vamos a ver cómo usar OGR, Python y Matplotlib para extraer una serie de datos de PostGIS, reproyectarlos y generar una cartografía. La idea en seleccionar los países que tienen su centroide dentro de una caja que representa la extensión de Europa, proyectar sus geometrías (originalmente en longitude latitude, WGS84) a la proyección sancionada por la directiva INSPIRE para Europa, ETRS-LAEA (Lambert Equal Area, datum geodético ETRS89, elipsoide GRS80). Para darle más emoción, cada polígono se va a dibujar como tal, y los vamos a colorear de acuerdo a un campo de la tabla. El código fuente del ejemplo está aquí (con sintaxis coloreada. Todo un lujo). En principio, es muy parecido al ejemplo anterior, con un par de salvedades.
La conexión OGR-PostGIS
La conexión con PostGIS se hace dándole al método Open de OGR el nombre de usuario, la base de datos que nos interesa y la dirección IP de la base de datos: 

g=ogr.Open("PG:host=127.0.0.1 dbname=db_gis user=gis")


donde

PG:
Indica base de datos PostGIS.
host=127.0.0.1
Indica la dirección IP de la base de datos (en este caso, localhost).
dbname=db_gis
Indica qué base de datos vamos a usar.
user
Indica el usuario de la base de datos.
La consulta espacial
La consulta espacial se delega en PostGIS. Nuestro programa únicamente creará la consulta en lenguaje SQL, se la pasará a PostGIS, y PostGIS nos devolverá las filas que la cumplan. Nuestra consulta va a escoger una serie de campos de la tabla world que 

  1. Tengan el centroide dentro de una caja geográfica determinada.
  2. Devuelva los datos ya proyectados.
El código es


europe_extent="GeomFromText('POLYGON((-25 35, 40 35, 40 85, -25 85, -25 35))',4326)" 

#Codigo PostGIS para LLAC-EEA: 3035 

L = g.ExecuteSQL("SELECT gid,cntry_name,color_map,transform(the_geom,3035) as the_geom 
 from world 
 where centroid(the_geom)@%s "%europe_extent )
 Vemos cómo la función transform(the_geom,SRID> transforma el campo geometría a nuestra proyección (en este caso, la 3035, LLAC-EEA). La función centroid(the_geom) calcula el centroide de una geometría, y el operador @ indica que se escojan aquellos features que caen dentro de otra geometría (en este caso, la geometría la hemos definido como una caja con límites inferior derecho cerca de Israel, y superior izquierdo cerca de Groenlandia). Esta consulta devuelve alrededor de 50 países.

Dibujando los polígonos
Para dibujar los polígonos, se usa los patches de Matplotlib. Este módulo define una serie de objetos útiles par representar formas cerradas (círculos, polígonos...). El código que usamos es:


m=transpose(r_[[x],[y]]) #Vértices del polígono
poligono = matplotlib.patches.Polygon( m ) #Definimos el polígono
#Le añadimos color al polígono, basado en un campo de la base de datos
poligono.set_facecolor ( color_map[feat.GetFieldAsInteger('color_map')-1])
#Se añade el polígono a la figura
ax.add_patch ( poligono )
El resultado final de nuestro ejemplo es el siguiente:
(también lo podemos sacar en EPS y pasarlo a PDF. Así quedaría en el PDF.
Conclusión
La combinación de python como lenguaje de programación, OGR como capa de enlace a datos geoespaciales, PostGIS como base de datos espacial y Matplotlib+scipy como herramientas de visualización y análisis es una herramienta para poder generar análisis espaciales complejos muy versátil. Y todo esto, sin usar software privativo, sin pagar un duro en licencias, y de manera bastante sencilla.