26.2 Un mapa como un ejemplo de campo

El siguiente ejemplo presenta un muestreo irregular de puntos, tal como lo podría almacenar un LIDAR, y su interpolación para ubicar las medidas en celdas regulares con el objeto de representar en puntos \(x, y\), la altura z del terreno.

Fuente de lo siguiente: https://towardsdatascience.com/plotting-regional-topographic-maps-from-scratch-in-python-8452fd770d9d

Se cargó desde https://topex.ucsd.edu/cgi-bin/get_data.cgi la información geográfica entre \(-4\) y \(8\) grados de latitud y \(-74.8\) y \(-73.8\) grados de longitud. Un buen pedazo del territorio nacional.

data = np.loadtxt('../Datos/Colombia_topografia.txt')
data

## array([[  280.975 ,    12.2637, -3481.    ],
##        [  280.9916,    12.2637, -3488.    ],
##        [  281.0083,    12.2637, -3490.    ],
##        ...,
##        [  293.475 ,    -4.138 ,    63.    ],
##        [  293.4916,    -4.138 ,    61.    ],
##        [  293.5083,    -4.138 ,    63.    ]])

Se observa que se trata de una matriz con tres columnas, que debemos imaginar son latitud, longitud y elevación.

from scipy.interpolate import griddata
Long = data[:,0]; Lat = data[:,1]; Elev = data[:,2]; #Variables
puntos=1000000;

Se definen tres vectores, cada uno con una de las columnas de la matriz presentada.

En el siguiente código, la primera línea hace uso de la función numpy.meshgrid() que se utiliza para crear una cuadrícula rectangular a partir de dos vectores (longitud y latitud). Esta cuadrícula sera de \(1000 \times 1000\) (definida a partir del objeto puntos).

Una base de datos topográfica como la que se cargó se corresponde probablemente con una grilla rectilínea, en donde se ha muestreado la elevación de las coordenadas a espaciados no uniformes. En las llanuras se puede muestrear la altura con bastante distancia debido a su uniformidad, en cambio, en las montañas, se debe hacer frecuentemente para captar el cambio de elevación.

La segunda línea hace uso de la función spyci.griddata() que interpola los puntos de la grilla a partir del vector de elevación. ¿Interpolar para qué? No coincide el espaciado regular de la cuadrícula de \(1000 \times 1000\) con los datos topográficos. La intención es ajustar el espaciado regular de la cuadrícula a partir del espaciado irregular de los datos.

[x,y]=np.meshgrid(np.linspace(np.min(Long), np.max(Long), round(np.sqrt(puntos))),
                  np.linspace(np.min(Lat), np.max(Lat), round(np.sqrt(puntos))));
z = griddata((Long, Lat), Elev, (x, y), method='linear');

En este momento los valores \(x, y\) y \(z\) representan Longitud, Latitud y Elevación, pero no son los valores cargados desde la base de datos, sino los valores ubicados en la cuadrícula de \(1000\times 1000\). Las siguientes tres líneas de código aplanan (flatten) cada matriz. Es decir, cada matriz de \(n\times n\) la convierte en un vector de longitud \(n\times n\) con el objeto de utilizar scatter de matplotlib:

x = np.matrix.flatten(x) # Longitud
y = np.matrix.flatten(y) # Latitud
z = np.matrix.flatten(z) # Elevación

import matplotlib.colors as mcolors
plt.gca().set_aspect('equal')
plt.scatter(x, y, 1, z, cmap="terrain")

## <matplotlib.collections.PathCollection object at 0x0000012D869390A0>

plt.colorbar(label='Elevación respecto al nivel del mar')

## <matplotlib.colorbar.Colorbar object at 0x0000012D86971940>

plt.xlabel('Longitud [°]')

## Text(0.5, 0, 'Longitud [°]')

plt.ylabel('Latitud [°]')

## Text(0, 0.5, 'Latitud [°]')

La longitud está medida respecto al meridiano de Greenwich, pero sobre \(360\)°.