Python y Sistemas de Información Geográfica

Python // Geospatial Analysis

Python y Sistemas de Información Geográfica

Herramientas fundamentales para el análisis geoespacial con Python

Este capítulo se enfoca en aplicar Python a funciones típicamente realizadas por un Sistema de Información Geográfica (GIS) como QGIS o ArcGIS. Estas funciones son el corazón del análisis geoespacial. Utilizaremos la menor cantidad posible de dependencias externas para que las herramientas sean lo más reutilizables posible en diferentes entornos.

Los temas presentados aquí son funciones centrales que sirven como bloques de construcción que puedes recombinar para resolver desafíos más allá de este libro. Cubriremos desde medir la Tierra hasta editar datos, crear mapas y escalar con multiprocesamiento.

Python UTM

OGR Library

PyShp

Fiona

PNGCanvas

Pillow/PIL

Folium

GeoPy

MÓDULO_01

Medición de Distancia

La esencia del análisis geoespacial es descubrir las relaciones de objetos en la Tierra. Los elementos que están más cerca tienden a tener una relación más fuerte que aquellos que están más lejos. Este concepto se conoce como la Primera Ley de Geografía de Tobler. Por lo tanto, medir la distancia es una función crítica del análisis geoespacial.

Cada mapa es un modelo de la Tierra y todos están equivocados en algún grado. Por esta razón, medir la distancia precisa entre dos puntos mientras estamos frente a una computadora es imposible. Para medir distancia, debemos considerar: ¿Qué estamos midiendo? ¿Cuánto estamos midiendo? ¿Cuánta precisión necesitamos?

Modelos de la Tierra

Plano Euclidiano

Geometría estándar para áreas pequeñas donde la curvatura es despreciable.

Esférico

Modelo que asume la Tierra como una esfera perfecta. Usa la fórmula Haversine.

Elipsoide

Modelo más preciso usando la fórmula de Vincenty para distancias exactas.

Distancia Euclidiana (Teorema de Pitágoras)

Fórmula Euclidiana

d = √[(x₂ - x₁)² + (y₂ - y₁)²]

pythagorean_distance.py

import math

def euclidean_distance(x1, y1, x2, y2):
    """Calcula la distancia euclidiana entre dos puntos"""
    dx = x2 - x1
    dy = y2 - y1
    return math.sqrt(dx**2 + dy**2)

# Ejemplo: Dos puntos en un mapa
point_a = (0, 0)
point_b = (3, 4)

distance = euclidean_distance(*point_a, *point_b)
print(f"Distancia: {distance} unidades")

python pythagorean_distance.py

Distancia: 5.0 unidades

Fórmula Haversine (Modelo Esférico)

Para distancias más largas sobre la superficie curva de la Tierra, usamos la fórmula Haversine que calcula la distancia del círculo máximo entre dos puntos.

haversine_distance.py

import math

def haversine(lon1, lat1, lon2, lat2):
    """Calcula la distancia entre dos puntos usando Haversine"""
    # Radio de la Tierra en kilómetros
    R = 6371
    
    # Convertir grados a radianes
    lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
    
    # Diferencias
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    
    # Fórmula Haversine
    a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
    c = 2 * math.asin(math.sqrt(a))
    
    return R * c

# Ejemplo: Nueva Orleans a Memphis
nola = (-90.07, 29.95)  # lon, lat
memphis = (-90.05, 35.15)

dist = haversine(*nola, *memphis)
print(f"Distancia: {dist:.2f} km")

python haversine_distance.py

Distancia: 577.98 km

La fórmula Haversine es suficientemente precisa para la mayoría de aplicaciones. Para máxima precisión geodésica, utiliza la fórmula de Vincenty que considera el elipsoide terrestre.

MÓDULO_02

Conversión de Coordenadas

Las coordenadas de mapas tradicionalmente se representaban como grados, minutos y segundos (DMS) para navegación marítima. Sin embargo, en GIS (basado en computadoras), la latitud y longitud se representan como números decimales conocidos como grados decimales.

Formatos de Coordenadas

DMS: 35° 8' 58.31" N

⟷

DD: 35.14953°

Conversión DMS a Grados Decimales

coordinate_conversion.py

import math
import re

def dd2dms(lat, lon):
    """Convierte grados decimales a DMS"""
    latf, latn = math.modf(lat)
    lonf, lonn = math.modf(lon)
    latd = int(latn)
    latm = int(latf * 60)
    lats = (lat - latd - latm / 60) * 3600.00
    lond = int(lonn)
    lonm = int(lonf * 60)
    lons = (lon - lond - lonm / 60) * 3600.00
    
    compass = {'lat': ('N','S'), 'lon': ('E','W')}
    lat_compass = compass['lat'][0 if latd >= 0 else 1]
    lon_compass = compass['lon'][0 if lond >= 0 else 1]
    
    return f"{abs(latd)}° {abs(latm)}' {abs(lats):.2f}\" {lat_compass}"

def dms2dd(lat, lon):
    """Convierte DMS a grados decimales"""
    lat_deg, lat_min, lat_sec, lat_dir = re.split(r'[^\d\.A-Z]+', lat)
    lon_deg, lon_min, lon_sec, lon_dir = re.split(r'[^\d\.A-Z]+', lon)
    
    lat_dd = float(lat_deg) + float(lat_min)/60 + float(lat_sec)/3600
    lon_dd = float(lon_deg) + float(lon_min)/60 + float(lon_sec)/3600
    
    if lat_dir == 'S': lat_dd *= -1
    if lon_dir == 'W': lon_dd *= -1
    
    return (lat_dd, lon_dd)

# Uso
print(dd2dms(35.14953, -90.04898))
print(dms2dd('''29° 56' 0.00" N''', '''90° 4' 12.36" W'''))

python coordinate_conversion.py

35° 8' 58.31" N
(29.933333, -90.0701)

MÓDULO_03

Reproyección de Datos Vectoriales

En GIS, la reproyección consiste en cambiar las coordenadas de un dataset de un sistema de coordenadas a otro. El módulo utm de Python funciona para conversión de sistemas de referencia, pero para una reproyección completa, necesitamos ayuda de la API OGR de Python.

Flujo de Reproyección

Lambert Conformal

―― OGR ――›

WGS84 Geographic

Reproyección con OGR

reproject_shapefile.py

from osgeo import ogr
from osgeo import osr
import os
import shutil

# Definir archivos fuente y destino
srcName = 'NYC_MUSEUMS_LAMBERT.shp'
tgtName = 'NYC_MUSEUMS_GEO.shp'

# Crear referencia espacial destino (WGS84)
tgt_spatRef = osr.SpatialReference()
tgt_spatRef.ImportFromEPSG(4326)

# Configurar el driver y abrir shapefile
driver = ogr.GetDriverByName('ESRI Shapefile')
src = driver.Open(srcName, 0)
srcLyr = src.GetLayer()
src_spatRef = srcLyr.GetSpatialRef()

# Crear transformación de coordenadas
trans = osr.CoordinateTransformation(src_spatRef, tgt_spatRef)

# Crear capa destino
tgt = driver.CreateDataSource(tgtName)
lyrName = os.path.splitext(tgtName)[0]
tgtLyr = tgt.CreateLayer(lyrName, geom_type=ogr.wkbPoint)
featDef = srcLyr.GetLayerDefn()

# Transformar cada feature
srcFeat = srcLyr.GetNextFeature()
while srcFeat:
    geom = srcFeat.GetGeometryRef()
    geom.Transform(trans)
    feature = ogr.Feature(featDef)
    feature.SetGeometry(geom)
    tgtLyr.CreateFeature(feature)
    srcFeat = srcLyr.GetNextFeature()

# Crear archivo .prj
tgt_spatRef.MorphToESRI()
with open(lyrName + '.prj', 'w') as prj:
    prj.write(tgt_spatRef.ExportToWkt())

El código EPSG 4326 representa el sistema WGS84 Geographic, el estándar mundial para GPS y mapas web.

MÓDULO_04

Calculando Área de Polígonos

La unidad más básica de GIS es un punto. Dos puntos forman una línea. Múltiples líneas que comparten puntos finales forman una polilínea, y las polilíneas pueden formar polígonos. Los polígonos representan todo, desde una casa hasta un país entero.

En GIS, los cálculos de área van más allá de la geometría básica. El polígono yace sobre la Tierra, que es una superficie curva. Afortunadamente, existe un módulo Python puro llamado area que maneja estas complicaciones.

calculate_area.py

from area import area

# Definir polígono en formato GeoJSON
polygon = {
    "type": "Polygon",
    "coordinates": [[
        [-89.324, 30.312],
        [-89.326, 30.31],
        [-89.322, 30.31],
        [-89.321, 30.311],
        [-89.321, 30.312],
        [-89.324, 30.312]
    ]]
}

# Calcular área en metros cuadrados
a = area(polygon)
print(f"Área: {round(a, 2)} m²")

python calculate_area.py

Área: 80235.14 m²

MÓDULO_05

Edición de Shapefiles

Los shapefiles son uno de los formatos de datos más comunes en GIS, tanto para intercambio de datos como para análisis. Para editar shapefiles, nos enfocamos principalmente en dos tipos de archivos:

.SHP

Contiene la geometría (puntos, líneas, polígonos)

.DBF

Contiene los atributos de la geometría correspondiente

.PRJ

Define el sistema de coordenadas (opcional)

Lectura de Atributos

read_shapefile.py

import shapefile

# Abrir shapefile
sf = shapefile.Reader("cities.shp")

# Obtener nombres de campos
fields = [field[0] for field in sf.fields[1:]]
print("Campos:", fields)

# Leer registros
for rec in sf.records():
    print(rec)

# Leer geometría
for shape in sf.shapes():
    print("Tipo:", shape.shapeType)
    print("Puntos:", shape.points)
    print("BBox:", shape.bbox)

Creación de Shapefiles

create_shapefile.py

import shapefile

# Crear writer para puntos
w = shapefile.Writer(shapefile.POINT)

# Definir campos
w.field('NAME', 'C', 40)
w.field('POPULATION', 'N', 10)

# Agregar registros con geometría
w.point(-90.07, 29.95)
w.record('New Orleans', 390144)

w.point(-90.05, 35.15)
w.record('Memphis', 651073)

# Guardar
w.save('ciudades')

En PyShp, cada registro de geometría debe coincidir con un registro de atributos. El orden de creación importa: primero geometría, luego registro.

MÓDULO_06

Selecciones y Consultas

Las selecciones son operaciones fundamentales para extraer subconjuntos de datos basados en criterios espaciales o de atributos.

Point-in-Polygon (Algoritmo Ray Casting)

point_in_polygon.py

def point_in_poly(x, y, poly):
    """Determina si un punto está dentro de un polígono"""
    n = len(poly)
    inside = False
    
    p1x, p1y = poly[0]
    for i in range(1, n + 1):
        p2x, p2y = poly[i % n]
        if y > min(p1y, p2y):
            if y <= max(p1y, p2y):
                if x <= max(p1x, p2x):
                    if p1y != p2y:
                        xinters = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
                    if p1x == p2x or x <= xinters:
                        inside = not inside
        p1x, p1y = p2x, p2y
    
    return inside

# Ejemplo
polygon = [(0,0), (10,0), (10,10), (0,10)]
print(point_in_poly(5, 5, polygon))  # True
print(point_in_poly(15, 5, polygon)) # False

Selección por Atributos con Fiona

attribute_selection.py

import fiona

# Seleccionar áreas urbanas con población < 5000
with fiona.open('MS_UrbanAreas.shp') as sf:
    # Filtrar por atributo
    filtered = filter(lambda f: f['properties']['POP'] < 5000, sf)
    
    # Obtener metadatos
    drv = sf.driver
    crs = sf.crs
    schm = sf.schema
    
    # Escribir subset
    with fiona.open('MS_Urban_Subset.shp', 'w',
                     driver=drv, crs=crs, schema=schm) as w:
        for rec in filtered:
            w.write(rec)

Las selecciones espaciales son computacionalmente costosas por los cálculos de punto flotante. Siempre intenta usar selección por atributos primero para reducir el conjunto de datos.

MÓDULO_07

Creación de Mapas Temáticos

¡Ahora pasamos de cálculos y edición de datos a algo visual! Crearemos dos tipos específicos de mapas temáticos usando Python puro.

Dot Density

Muestra concentraciones de datos usando puntos distribuidos aleatoriamente dentro de un área, con una proporción fija a través del dataset. Comúnmente usado para mapas de densidad de población.

Choropleth

Usa sombreado, colores o símbolos para mostrar un valor promedio o cantidad dentro de un área. Facilita visualizar grandes cantidades de datos como resumen.

Mapa Dot Density con PNGCanvas

dot_density_map.py

import shapefile
import pngcanvas
import random

def world2screen(bbox, w, h, x, y):
    """Convierte coordenadas geoespaciales a píxeles"""
    minx, miny, maxx, maxy = bbox
    xratio = w / (maxx - minx)
    yratio = h / (maxy - miny)
    px = int(w - ((maxx - x) * xratio))
    py = int((maxy - y) * yratio)
    return (px, py)

# Abrir shapefile de censo
inShp = shapefile.Reader('GIS_CensusTract_poly')
iwidth, iheight = 600, 400

# Encontrar índice del campo de población
pop_index = None
for i, f in enumerate(inShp.fields):
    if f[0] == 'POPULATION':
        pop_index = i - 1

# Calcular densidad (1 punto por cada 100 personas)
dots = []
for sr in inShp.shapeRecords():
    population = sr.record[pop_index]
    density = population // 100
    found = 0
    while found < density:
        minx, miny, maxx, maxy = sr.shape.bbox
        x = random.uniform(minx, maxx)
        y = random.uniform(miny, maxy)
        if point_in_poly(x, y, sr.shape.points):
            dots.append((x, y))
            found += 1

# Crear imagen PNG
c = pngcanvas.PNGCanvas(iwidth, iheight)
c.color = (255, 0, 0, 0xff)  # Rojo

for d in dots:
    x, y = world2screen(inShp.bbox, iwidth, iheight, *d)
    c.filled_rectangle(x-1, y-1, x+1, y+1)

with open('DotDensity.png', 'wb') as img:
    img.write(c.dump())

Mapa Choropleth con PIL

choropleth_map.py

import math
import shapefile
from PIL import Image, ImageDraw

# Abrir shapefile
inShp = shapefile.Reader('GIS_CensusTract_poly')
iwidth, iheight = 600, 400

# Crear imagen PIL
img = Image.new('RGB', (iwidth, iheight), (255,255,255))
draw = ImageDraw.Draw(img)

# Obtener índices de campos
pop_index, area_index = None, None
for i, f in enumerate(inShp.fields):
    if f[0] == 'POPULATION': pop_index = i - 1
    elif f[0] == 'AREASQKM': area_index = i - 1

# Dibujar polígonos con color basado en densidad
for sr in inShp.shapeRecords():
    density = sr.record[pop_index] / sr.record[area_index]
    weight = min(math.sqrt(density / 80.0), 1.0) * 50
    
    # Calcular color (más oscuro = más denso)
    R = int(205 - weight)
    G = int(215 - weight)
    B = int(245 - weight)
    
    pixels = [world2screen(inShp.bbox, iwidth, iheight, x, y) 
              for x, y in sr.shape.points]
    draw.polygon(pixels, outline=(255,255,255), fill=(R,G,B))

img.save('choropleth.png')

MÓDULO_08

Creación de Heat Maps

Un heat map muestra el agrupamiento geográfico de datos usando una imagen raster que muestra densidad. Usaremos la biblioteca Folium para crear mapas web interactivos basados en Leaflet.

heat_map.py

import folium
from folium.plugins import HeatMap

# Leer datos de avistamientos
with open('sightings.csv', 'r') as f:
    lines = f.readlines()
    lines.pop(0)  # Quitar header

# Parsear coordenadas
data = [list(map(float, l.strip().split(','))) for l in lines]

# Crear mapa base
m = folium.Map(
    [32.75, -89.52],
    tiles='stamentonerbackground',
    zoom_start=7
)

# Agregar capa de calor
HeatMap(data, max_zoom=16, radius=22, 
        min_opacity=1, blur=30).add_to(m)

# Guardar como HTML
m.save('heatmap.html')

El archivo HTML generado puede abrirse en cualquier navegador web y proporciona interactividad completa con zoom y pan.

MÓDULO_09

Uso de Datos GPS

El formato GPS más común actualmente es Garmin GPX, un formato XML. Sin embargo, existe otro tipo de datos GPS más antiguo llamado NMEA (National Marine Electronics Association), que son sentencias ASCII diseñadas para streaming.

parse_nmea.py

from pynmea.streamer import NMEAStream

# Ejemplo de sentencias NMEA
# $GPRMC,012417.859,V,1856.599,N,15145.602,W,12.0,7.27,020713,,E*4F
# $GPGGA,012418.859,1856.599,N,15145.602,W,0,00,,,M,,M,,*54

# Abrir archivo NMEA
nmeaFile = open('nmea.txt')
nmea_stream = NMEAStream(stream_obj=nmeaFile)

# Parsear objetos
nmea_objects = []
next_data = nmea_stream.get_objects()
while next_data:
    nmea_objects += next_data
    next_data = nmea_stream.get_objects()

# Extraer coordenadas
for nmea_ob in nmea_objects:
    if hasattr(nmea_ob, 'lat'):
        print(f'Lat/Lon: ({nmea_ob.lat}, {nmea_ob.lon})')

Las coordenadas NMEA se almacenan en formato "grados minutos decimales". Por ejemplo, 4533.35 = 45° 33.35' = 45° 33' 21".

MÓDULO_10

Geocodificación

La geocodificación es el proceso de convertir una dirección postal en latitud y longitud. Esta operación es crítica para sistemas de navegación y sitios web de direcciones de manejo.

◉

geocoder

Múltiples servicios (Google, Bing, Yahoo)

◉

geopy

OpenStreetMap, reverse geocoding

geocoding.py

import geocoder
from geopy.geocoders import Nominatim

# Usando geocoder (Google Maps)
g = geocoder.google('1403 Washington Ave, New Orleans, LA 70130')
print(g.wkt)
# POINT(-90.08421849999999 29.9287839)

# Usando geopy (OpenStreetMap)
geolocator = Nominatim(user_agent="myapp")
location = geolocator.geocode('88360 Diamondhead Dr E, Diamondhead, MS')

print(f"Lat: {location.latitude}, Lon: {location.longitude}")

# Reverse geocoding
rev = geolocator.reverse(f"{location.latitude},{location.longitude}")
print(rev.address)

MÓDULO_11

Multiprocesamiento

Los datasets geoespaciales son muy grandes. Procesarlos puede tomar horas o incluso días. El módulo multiprocessing de Python puede generar múltiples procesos para aprovechar todos los procesadores disponibles.

Pool de Procesadores

Lista de Ciudades

―→

Pool (CPU cores)

―→

Resultados Paralelos

multiprocessing_geocode.py

from geopy.geocoders import Nominatim
import multiprocessing as mp

# Crear geocoder
g = Nominatim(user_agent="myapp")

def gcode(address):
    """Geocodifica una dirección individual"""
    location = g.geocode(address)
    print(f"Geocodificando: {address}")
    return location

# Lista de ciudades a procesar
cities = [
    "New Orleans, LA",
    "Biloxi, MS",
    "Memphis, TN",
    "Atlanta, GA",
    "Little Rock, AR",
    "Destin, FL"
]

# Crear pool de procesadores
pool = mp.Pool(processes=mp.cpu_count())

# Mapear ciudades a la función de geocodificación
results = pool.map(gcode, cities)

# Imprimir resultados
for r in results:
    print(r)

El multiprocesamiento funciona mejor con operaciones que pueden dividirse en cálculos discretos. No todo tipo de procesamiento puede paralelizarse de esta manera.

Resumen de Módulos

01 // Medición de Distancia

Fórmulas Euclidiana, Haversine y Vincenty para diferentes precisiones y escalas.

02 // Conversión de Coordenadas

Transformación entre DMS y grados decimales para compatibilidad de datos.

03 // Reproyección

Cambio de sistemas de coordenadas usando OGR/OSR para integración de datos.

04 // Cálculo de Área

Área de polígonos considerando la curvatura terrestre con el módulo area.

05 // Edición de Shapefiles

Lectura, escritura y modificación de archivos shapefile con PyShp y Fiona.

06 // Selecciones

Consultas espaciales y por atributos para extracción de subconjuntos de datos.

07 // Mapas Temáticos

Visualización con Dot Density y Choropleth usando PNGCanvas y PIL.

08 // Heat Maps

Mapas web interactivos con Folium y Leaflet para clustering geográfico.

09 // Datos GPS

Parsing de formatos GPX y NMEA para datos de campo y tracking.

10 // Geocodificación

Conversión de direcciones a coordenadas con geocoder y geopy.

11 // Multiprocesamiento

Aceleración de procesos usando múltiples CPUs para datasets grandes.

∞ // Próximo Paso

Sensores Remotos: procesamiento de imágenes satelitales y aéreas.

Archivo

Categorias

Visitas al Blog

Acerca de Mi

jueves, 22 de enero de 2026