Explorando datos de población a gran escala en formato raster

Explorando datos de población en formato raster a gran escala

Visualizando datos de población geoespacial en varias escalas usando Python: datos globales, por país y a nivel urbano

He visto con frecuencia mapas de población hermosos circulando en línea; sin embargo, a menudo me quedaba atascado en algunas partes técnicas, como visualizar otros segmentos de mapa que no se muestran en el tutorial, o convertir los datos de raster a formatos vectoriales más amigables para el cálculo. Supero algunas de estas limitaciones en este artículo con una guía práctica de dos fuentes principales de datos de población global.

También es importante destacar que, además de su valor estético, los datos de población y los mapas que los muestran son información fundamental y valiosa que se puede recopilar e incorporar para cualquier tarea de desarrollo urbano o inteligencia de ubicación. Son especialmente útiles en casos de uso como la planificación de nuevas comodidades, la selección de sitios y el análisis de alcance, la estimación de la escala de productos urbanos o el perfilado de diferentes vecindarios, por mencionar solo algunos.

1. Fuentes de datos

Me baso en las siguientes dos fuentes de datos de estimación de población detallada, de las cuales puedes descargar los archivos a través de los enlaces adjuntos (a la fecha de publicación):

• El GHSL (Global Human Settlement Layer) de la Comisión Europea mide el nivel de población de cada celda de la cuadrícula. Encuentra la descripción general aquí y el conjunto de datos específico que utilicé de su informe 2023 con una resolución espacial de 100m aquí.
• WorldPop hub, desde donde tomaré Alemania como ejemplo utilizando el conjunto de datos de países individuales restringidos con una resolución de 100m. Encuentra el listado aquí y los datos de Alemania aquí.

2. Visualizando la capa global de asentamientos humanos

2.1. Importar los datos

Me encontré por primera vez con este conjunto de datos en el tutorial de Datashader del Architecture Performance. Después de reproducir su visualización, me encontré con algunos obstáculos al ampliarla a un mapa global, lo que inició este trabajo, así que ahora te muestro las soluciones alternativas que encontré sobre cómo hacerlo.

• 9 LLMs y Agentes de Código Abierto para Observar
• Inteligencia Artificial Tradicional vs Inteligencia Artificial Gene...
• Construyendo una Búsqueda de Vectores de Grado de Producción

Primero, analiza el archivo raster utilizando el paquete xarray.

import rioxarrayfile_path  = "GHS_POP_E2030_GLOBE_R2023A_54009_100_V1_0/GHS_POP_E2030_GLOBE_R2023A_54009_100_V1_0.tif" data_array = rioxarray.open_rasterio(file_path, chunks=True, lock=False)data_array

La salida de esta celda es una descripción detallada del conjunto de datos:

2.2. Visualizando segmentos de los datos

Ya podemos ver que esta es una cantidad de datos bastante desafiante para la mayoría de las computadoras portátiles estándar. De todos modos, intentemos visualizarlo utilizando Datashader, una herramienta realmente útil para conjuntos de datos geoespaciales de esta magnitud:

# ADVERTENCIA: ESTE BLOQUE DE CÓDIGO PROBABLEMENTE PROVOQUE UN ERROR DE DESBORDAMIENTO DE MEMORIAimport datashader as dsimport xarray as xrfrom colorcet import palettefrom datashader import transfer_functions as tf# preparar para trazar los datosdata_array_p = xr.DataArray(data_array)[0]data_array_p = data_array_p.where(data_array_p > 0)data_array_p = data_array_p.compute()# obtener el tamaño de la imágenesize = 1200asp  = data_array_p.shape[0] / data_array_p.shape[1]# crear el lienzo de Datashadercvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))raster = cvs.raster(data_array_p)# dibujar la imagen en el lienzo de Datashadercmap = palette["fire"]img = tf.shade(raster, how="eq_hist", cmap=cmap)img

Aunque este código parece técnicamente correcto, mi MacBook Pro M1 de 2021 con 16GB de RAM arroja un miserable error de desbordamiento de memoria. Así que, ¡recortemos la imagen para observar los datos! Para esto, seguiré el enfoque del Architecture Performance y me centraré en Europa, que seguiré por el momento, ya que definitivamente funciona.

Sin embargo, la pregunta principal que responderé más adelante es cómo podemos visualizar los datos de todo el globo a pesar de estas limitaciones de memoria, pero aún utilizando nuestra máquina local. ¡Espera por ello!

import datashader as dsimport xarray as xrfrom colorcet import palettefrom datashader import transfer_functions as tfimport numpy as np# recortar el arreglo de datosdata_array_c = data_array.rio.clip_box(minx=-1_000_000.0, miny=4_250_000.0, maxx=2_500_000.0, maxy=7_750_000.0)data_array_c = xr.DataArray(data_array_c)# preparar para graficardata_array_c = xr.DataArray(data_array_c)[0]data_array_c = data_array_c.where(data_array_c > 0)data_array_c = data_array_c.compute()data_array_c = np.flip(data_array_c, 0)# obtener el tamaño de la imagentamaño = 1200asp  = data_array_c.shape[0] / data_array_c.shape[1]# crear el lienzo de Data Shadercvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))raster = cvs.raster(data_array_c)# dibujar la imagencmap = palette["fire"]img = tf.shade(raster, how="eq_hist", cmap=cmap)img = tf.set_background(img, "black")img

Este bloque de código produce la siguiente visualización:

El uso del mapa de colores ‘fire’ parece ser un estándar de la industria por una buena razón; sin embargo, si deseas variar un poco, puedes encontrar otros esquemas de color aquí y aplicarlos a continuación:

# crear el lienzo de Data Shadercvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))raster = cvs.raster(data_array_c)# dibujar la imagencmap = palette["bmw"]img = tf.shade(raster, how="eq_hist", cmap=cmap)img = tf.set_background(img, "black")img

Este bloque de código produce la siguiente visualización:

2.3. Visualizar el globo entero

Entonces, los datos están ahí, pero ¿qué pasa si tienes una computadora normal a mano y aún quieres visualizar todo el mundo con una resolución de 100 m? La solución alternativa que te mostraré aquí es relativamente simple: divido la imagen de raster completa en alrededor de cien fragmentos más pequeños para que mi computadora los maneje uno por uno, luego uso algunos trucos de manipulación de imágenes para fusionarlos en un solo archivo de imagen.

Sin embargo, antes de continuar, una nota rápida. También hay una opción para reducir el tamaño de los arreglos de XArray de la siguiente manera, pero no pude encontrar una reducción adecuada que pudiera manejar todo el conjunto de datos. Además, no quería perder precisión y quería ver el conjunto de datos completo en toda su gloria.

# una forma rápida de reducir el tamañofactor_de_reducción = 20data_array_reducido = data_array.coarsen(x=factor_de_reducción, y=factor_de_reducción).mean()data_array_reducido

Y la salida, que vale la pena comparar con el data_array previamente graficado:

Para dividir la imagen de raster completa en segmentos de cuadrícula, primero obtén sus límites y define N como el tamaño del paso. Luego, crea la lista de límites de cada segmento de imagen.

minx = float(data_array.x.min().values)maxx = float(data_array.x.max().values)miny = float(data_array.y.min().values)maxy = float(data_array.y.max().values)N = 10xstep = (maxx-minx) / Nystep = (maxy-miny) / Nxsteps = list(np.arange(minx, maxx, xstep)) ysteps = list(np.arange(miny, maxy, ystep))

Ahora, iterar sobre cada paso x e y y crear cada segmento de imagen, donde cada archivo de imagen se llama según su posición en la cuadrícula original. Este bucle puede llevar un tiempo.

import osfoldout = 'world_map_image_segments'if not os.path.exists(foldout):    os.makedirs(foldout)for idx_x, x_coord in enumerate(xsteps):    for idx_y, y_coord in enumerate(ysteps):            if not os.path.exists(foldout+'/'+str(idx_x)+'_'+str(idx_y)+'.png'):                        data_array_c = data_array.rio.clip_box( minx=x_coord,  miny=y_coord,  maxx=x_coord+xstep, maxy=y_coord+ystep)            data_array_c = xr.DataArray(data_array_c)[0]            data_array_c = data_array_c.fillna(0)            data_array_c = data_array_c.where(data_array_c > 0)            data_array_c = data_array_c.compute()            data_array_c = np.flip(data_array_c, 0)            size = 2000            asp  = data_array_c.shape[0] / data_array_c.shape[1]            cvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))            raster = cvs.raster(data_array_c)            cmap = palette["fire"]            img = tf.shade(raster, how="eq_hist", cmap=cmap)            img = tf.set_background(img, "black")            pil_image = img.to_pil()            pil_image.save(foldout+'/'+str(idx_x)+'_'+str(idx_y)+ '.png')            print('GUARDADO: ', x_coord, y_coord, y_coord+xstep,y_coord+ystep)

Finalmente, si tenemos todos los segmentos de imagen, podemos ensamblarlos rápidamente con la siguiente función. Para este fragmento de código, también pedí algunos consejos a ChatGPT para acelerar las cosas, pero como de costumbre, también necesitó algunos ajustes manuales.

from PIL import Imagedef find_dimensions(image_dir):    max_x = 0    max_y = 0    for filename in os.listdir(image_dir):        if filename.endswith(".png"):            x, y = map(int, os.path.splitext(filename)[0].split("_"))            max_x = max(max_x, x)            max_y = max(max_y, y)    return max_x + 1, max_y + 1 image_dir = foldoutsegment_width = sizesegment_height = int(asp*size)# Determinar las dimensiones de la imagen grandelarge_image_width, large_image_height = find_dimensions(image_dir)# Crear una imagen grande vacía (fondo blanco)large_image = Image.new("RGB", (large_image_width * segment_width, large_image_height * segment_height), "black")# Recorrer los segmentos de imagen individuales y pegarlos en la imagen grandepara filename in sorted(os.listdir(image_dir)):    if filename.endswith(".png"):        x, y = map(int, os.path.splitext(filename)[0].split("_"))        segment_image = Image.open(os.path.join(image_dir, filename))        x_offset = x * segment_width        y_offset = large_image_height * segment_height-1*y * segment_height        large_image.paste(segment_image, (x_offset, y_offset))# Guardar la imagen grande fusionadalarge_image.save("global_population_map.png")     

Y el resultado final aquí, el mapa completo del mundo:

3. Visualizando y transformando los datos de WorldPop

La segunda fuente de datos que me gustaría mostrarles es la base de datos de población de WorldPop, que tiene continentes y países por separado en varias resoluciones. En este ejemplo, complementando los niveles de continente y global de la sección anterior, aquí me enfoco en el nivel de países y ciudades. Por ejemplo, elijo Alemania y una resolución de 100m curada en 2020, y también les muestro cómo recortar una ciudad del país completo y convertirla en un formato vectorial fácil de usar utilizando GeoPandas.

3.1. Visualizar datos de WorldPop

Usando los métodos anteriores, nuevamente podemos visualizar rápidamente este archivo raster:

# analizar el archivodata_file = 'deu_ppp_2020_constrained.tif'data_array = rioxarray.open_rasterio(data_file, chunks=True, lock=False)# preparar los datosdata_array = xr.DataArray(data_array)[0]data_array = data_array.where(data_array > 0)data_array = data_array.compute()data_array = np.flip(data_array, 0)# obtener el tamaño de la imágenesize = 1200asp  = data_array.shape[0] / data_array.shape[1]# crear el lienzo de sombreador de datoscvs = ds.Canvas(plot_width=size, plot_height=int(asp*size))raster = cvs.raster(data_array)# dibujar la imagencmap = palette["fire"]img = tf.shade(raster, how="eq_hist", cmap=cmap)img = tf.set_background(img, "black")img

Este bloque de código muestra la siguiente visualización:

3.2. Transformando los datos de WorldPop

Después de visualizar el planeta completo, el continente de Europa y el país de Alemania, me gustaría trabajar más de cerca con la ciudad de Berlín y mostrarte cómo transformar datos ráster en formato vectorial y manipularlos fácilmente con GeoPandas. Para esto, accedo a los límites administrativos de Berlín en formato geojson aquí.

Este archivo admin contiene los distritos de Berlín, por lo que primero los fusiono en la ciudad como un todo.

from shapely.ops import cascaded_unionimport geopandas as gpdadmin = gpd.read_file('tufts-berlin-bezirke-boroughs01-geojson.json')admin = gpd.GeoDataFrame(cascaded_union(admin.geometry.to_list()), columns = ['geometry']).head(1)admin.plot()

Este bloque de código muestra la siguiente visualización:

Ahora convierte el xarray en un DataFrame de Pandas, extrae la información de la geometría y crea un GeoDataFrame de GeoPandas. Una forma de hacer esto es:

import pandas as pddf_berlin = pd.DataFrame(data_array.to_series(), columns = ['population']).dropna()

Ahora, crea un GeoDataFrame a partir de esto, enfocándote en Berlín:

from shapely.geometry import Point# encuentra la caja delimitadora para una selección más fácil de las coordenadasminx, miny, maxx, maxy = admin.bounds.T[0].to_list()points = []population = df_berlin.population.to_list()indicies   = list(df_berlin.index)# crea geometrías de punto a partir de los puntos que caen dentro de esta caja delimitadorageodata = []for ijk, (lon, lat) in enumerate(indicies):    if minx <= lat <= maxx and miny <= lon <= maxy:           geodata.append({'geometry' : Point(lat, lon), 'population' : population[ijk]})        # construye un GeoDataFramegdf_berlin = gpd.GeoDataFrame(geodata)gdf_berlin = gpd.overlay(gdf_berlin, admin)

Luego, visualiza la población como datos vectoriales:

import matplotlib.pyplot as pltf, ax = plt.subplots(1,1,figsize=(15,15))admin.plot(ax=ax, color = 'k', edgecolor = 'orange', linewidth = 3)gdf_berlin.plot(column = 'population',                 cmap = 'inferno',                 ax=ax,                 alpha = 0.9,                 markersize = 0.25)ax.axis('off')f.patch.set_facecolor('black')

Este bloque de código muestra la siguiente visualización:

Finalmente, aquí tenemos un GeoDataFrame estándar con niveles de población de resolución de 100m asignados a cada geometría de punto correspondiente a cada píxel en el archivo ráster.

Resumen

En este artículo, exploré dos conjuntos de datos de población global basados en varias aproximaciones, que estiman niveles de población combinando diferentes métodos de aproximación, medición y modelado a una notable resolución espacial de 100m utilizando cuadrículas ráster. Este tipo de información es muy valiosa para una amplia gama de aplicaciones en desarrollo urbano e inteligencia de ubicación, como planificación de infraestructura, selección de sitios, perfilado de vecindarios y más. Desde el lado técnico, mostré ejemplos en tres niveles espaciales, cubriendo el planeta entero, luego acercándome a los países y finalmente a las ciudades. Si bien la metodología puede manejar incluso resoluciones más pequeñas, todo esto sucedió en una sola computadora portátil utilizando poderosas bibliotecas de Python como Xarray, DataShader y GeoPandas.

• Data
• Geospatial
• GIS
• Maps
• Python