4. Scripts de ejemplo

En esta sección se muestran algunos ejemplos de scripts de MATLAB para mostrar datos geoespaciales. Se utilizan las funciones descritas en la sección anterior para cargar datos, crear mapas y mostrar puntos en un mapa. Se recomienda revisar la documentación de cada función para obtener más información sobre su uso y opciones disponibles.

4.1. Datos para realizar las ejemplos

Enlace a los datos

Región

Formato

Fuente

Datos de sismos (2010-2023)

Global

csv

USGS

Temperatura superficial del mar (2015)

Bahía de La Paz

mat

ERDDAP

Temperatura superficial del mar (2018)

Sur del GC

mat

ERDDAP

Clorofila (2020)

Sur del GC

mat

ERDDAP

Clorofila (2015)

Sur del GC

nc

ERDDAP

Temperatura superficial del mar (2015)

Sur del GC

nc

ERDDAP

Temperatura superficial del mar (2022)

Sur del GC

nc

ERDDAP

Corrientes marinas CMEMS (2023)

Golfo de Mexico

nc

Copernicus

Corrientes geostróficas (2023)

Atlántico ecuatorial

nc

Copernicus

4.2. Con datos vectoriales

4.2.1. Tsunamis en el mundo y en México

Práctica #1

Crear varios tipos de mapas para representar los tsunamis almacenados en el archivo tsunamis.xlsx. Para este ejemplo, no se requiere la descarga de datos adicionales. Los que se utilizan estan disponibles con la instalación de Mapping Toolbox.

  1. Utiliza la función readtable para cargar los datos desde el archivo de Excel.

  2. Crea los ejes del mapa utilizando la función worldmap.

  3. Utiliza la función geoshow para mostrar los continentes (landareas).

  4. Muestra la ubicación de los tsunamis utilizando la función geoshow.

 1% Cargar los datos (tsunamis y continentes)
 2tsunamis = readtable('tsunamis.xlsx');
 3land = readgeotable('landareas.shp');
 4
 5% Crear los ejes del mapa
 6figure
 7worldmap('World')
 8
 9% Mostrar los continentes (polígonos) y tsunamis (puntos)
10geoshow(land, 'DisplayType', 'polygon')
11geoshow(tsunamis.Latitude, tsunamis.Longitude, 'DisplayType', 'point')

  1. Modifica el script para mejorar la apariencia del mapa (visualización de los tsunamis y los continentes).

  2. Agrega instrucciones al script para crear otro mapa, esta vez de los sismos ocurridos en México.

Nota

Sugerencia: Primero filtra los datos de los tsunamis para seleccionar solo aquellos que ocurrieron en México (por ejemplo, usando las funciones find y contains en la columna que almacena el nombre del país). Posteriormente la función worldmap para crear un mapa de México y la función geoshow para mostrar los puntos de los tsunamis en el mapa.

  1. Agrega un título al mapa utilizando la función title.

  2. Guarda el mapa como una imagen utilizando la función print, exportgraphics o saveas.

  3. Genera un mapa adicional (para todo el mundo), que represente las ubicaciones de los tsunamis (puntos) de tamaño diferente (dependiendo de la “magnitud”) y con la escala de color según la “causa”. Para esto, utiliza la función geobubble. También agrega un título y guarda el mapa como una imagen (.png).

  4. Guarda el script como xyz_00_tsunamis.mlx (xyz = tus iniciales) en la carpeta matlab.

  5. Sube el script y los mapas (.png) a la plataforma Teams.

4.2.2. Sismos

Sismos en México (2010-2024)

Sismos en México (2010-2024)

Práctica #2

Tomando como referencia el script anterior, crear 14 mapas con la ubicación de los sismos en México por año (2010-2024). Además, crear un mapa con todos los sismos en México en ese periodo.

  1. Descarga los “Datos de sismos” (2010-2023) de la sección de datos. Además, descarga de la página earthquake.usgs.gov los datos de los sismos (≥4 de magnitud) ocurridos durante el año 2024 en todo el mundo (en formato csv). Guargar este archivo como sismos-2024.csv. Transfiere estos 14 archivos csv a la carpeta datos/csv.

  2. Modificar el script para acceder a la carpeta datos/csv y leer todos los archivos csv (sismos de 2010-2024).

  3. Para cada archivo csv seleccionar los sismos ocurridos en México, cuya magnitud fue ≥ 4.

Nota

Sugerencia: Puedes utilizar la función dir para obtener una lista de los archivos en la carpeta y luego leer cada archivo en un bucle. Por ejemplo, puedes usar dir('../datos/csv/sismos-*.csv') para obtener la lista de archivos csv. Luego, usar un ciclo for para iterar sobre cada archivo y cargar los datos utilizando la función readtable. El procedimiento para filtrar los registros de interés es el mismo que en el ejercicio anterior (funciones find y contains).

  1. Generar los mapas de sismos (uno por año) utilizando la función geoscatter. Mostrar una barra de color que indique la magnitud de los sismos. En el título de cada mapa especificar el año y el número de sismos en el año. Guardar estos mapas en formato png a 300 dpi.

  2. Generar un mapa que incluya todos los sismos ocurridos en México (2010-2024), similar a los obtenidos según el punto anterior. En este caso utilizar la función geodensityplot para presentar el mapa. Utilizar grayterrain como mapa base. ¿Cuántos sismos ocurrieron en ese periodo? (mostrarlo en el título del mapa). Guardar el mapa en formato png a 300 dpi.

  3. Crear una gráfica que muestre la magnitud de los sismos con magnitud ≥ 7 grados en el periodo 2010-2024. ¿Cuándo (fecha y hora) ocurrieron esos sismos? (según la zona horaria 'America/Mexico_City').

Nota

Utilizar la función datetime para convertir las fechas y horas de los sismos a un formato de fechas, así como para especificar la zona horaria (con el argumento 'TimeZone', 'UTC'). Después de crear la variable de tipo datetime puedes configurar propiedad TimeZone para cambiar la zona horaria. Así, cambia la zona horaria a 'America/Mexico_City'.

  1. Documentar el script con los comentarios pertinentes de acuerdo con los datos y procedimientos utilizados.

  2. Como evidencia, adjuntar en Teams el script (xyz_01_sismos.pdf). Debe incluir todos los mapas y la gráfica. Además, subir el archivo png generado con la función geodensityplot y uno de los archivos png obtenido con la función geoscatter.

Ejemplo de script para leer todos los archivos en una carpeta, especificando la ruta a la carpeta, parte del nombre de los archivos y su extensión (en este caso csv):

 1% Recuperar los nombres de archivos CSV de sismos en la carpeta datos/csv
 2archivos = dir('../datos/csv/sismos-*.csv');
 3
 4% Inicializar tabla para almacenar todos los sismos
 5todos_sismos = table();
 6
 7% Ciclo para leer cada archivo CSV
 8for k = 1:length(archivos)
 9    % Leer el archivo k
10    sismos = readtable(archivos(k).name);
11
12    % Filtrar sismos en México con magnitud >= 4
13    rens = find(contains(sismos.place, 'Mexico') & sismos.mag >= 4);
14
15    % Crear tabla con sismos en México y magnitud >= 4
16    sismos_mex = sismos(rens, :);
17
18    % Concatenar a la tabla general
19    todos_sismos = [todos_sismos; sismos_mex];
20
21    % Crear mapa para cada año...
22end
23
24% Crear el mapa con los sismos (2010-2024) en México...

4.2.3. Leer archivos shp

Los archivos “shapefile” aun son ampliamente utilizados para almacenar datos geoespaciales de tipo vectorial. En este ejemplo, se muestra cómo leer un archivo shp utilizando la función shaperead. La función shaperead permite incluir condiciones para filtrar los registros. Devuelve una estructura tipo struct con los atributos almacenados en el “shapefile”.

1% Importar datos de un shapefile (solo los estados de interés)
2estados_interes = shaperead("gadm41_MEX_1.shp", "UseGeoCoords", true, "Selector", ...
3    {@(name) any(contains(name, {'Baja California', 'Sonora', 'Sinaloa'})), 'NAME_1'});
4
5% Crear mapa
6figure
7worldmap("Mexico")
8geoshow(estados_interes, "DisplayType", "polygon", "FaceColor", [0.5 0.5 0.5])
9title("Estados de interés")

Con shaperead también podemos especificar la región (BoundingBox) que nos interesa extraer del “shapefile”, como se muestra a continuación:

Nota

En este caso, la región se especifica con una matriz de 2×2 de la forma [lonmin,latmin; lonmax,latmax] (para coordenadas geográficas). La función shaperead no recorta las entidades que intersecan parcialmente el rectángulo. El parámetro Àttributes permite seleccionar los atributos (columnas) que se desean incluir en la estructura de salida.

1% Importar datos de un shapefile (región específica)
2bbox = [-117, 20; -105, 32]; % [longitud, latitud]
3estados_interes = shaperead("gadm41_MEX_1.shp", "UseGeoCoords", true, "BoundingBox", bbox, "Attributes", {'NAME_1'});
4
5% Crear mapa
6figure
7worldmap([bbox(1,2),bbox(2,2)],[bbox(1,1), bbox(2,1)]) % ([limlat], [limlon])
8geoshow(estados_interes, "DisplayType", "polygon", "FaceColor", [0.5 0.5 0.5])
9title("Estados en la región especificada")

También se puede usar la función readgeotable para leer “shapefiles”. Esta función es una alternativa más reciente para acceder a archivos vectoriales. Devuelve un tipo de estructura table:

 1% Importar datos de un shapefile
 2estados = readgeotable("gadm41_MEX_1.shp");
 3
 4% Filtrar los estados de interés
 5estados_interes = estados(contains(estados.NAME_1, {'Baja California', 'Sonora', 'Sinaloa'}), :);
 6
 7% Crear mapa
 8figure
 9worldmap("Mexico")
10geoshow(estados_interes, "DisplayType", "polygon", "FaceColor", [0.5 0.5 0.5])
11title("Estados de interés")

4.3. Con datos ráster

4.3.1. Archivos mat

TSM en la bahía de La Paz (2015) y en el sur del GC (2018)

En este ejemplo, creamos un mapa que muestra la temperatura superficial del mar (TSM) en la Bahía de La Paz (2015). Estos datos están disponibles en la sección 4.1 en formato .mat y se pueden leer utilizando la función load. La variable de interés es jplG1SST.SST.

Nota

Nota que las variables están contenidas en una estructura (tipo struct) llamada jplG1SST. Los datos de TSM están en un arreglo (3D) con dimensiones [tiempo, latitud, longitud]. Para crear el mapa, primero es necesario obtener una referencia geográfica utilizando la función georefcells, posteriormente, seleccionar una capa específica (por ejemplo, la primera capa) y utilizar la función geoshow para mostrar los datos en un mapa.

 1% Importar datos de TSM
 2load jplG1SST_blap_2015.mat;
 3
 4% Obtener los datos de la variable (struct) jplG1SST
 5% (Verificar el tipo de datos)
 6SST = jplG1SST.SST;
 7lat = jplG1SST.latitude;
 8lon = jplG1SST.longitude;
 9
10% Obtener los límites de latitud y longitud
11latlim = [lat(1), lat(end)];
12lonlim = [lon(1), lon(end)];
13
14% Consular el tamaño del arreglo (3D)
15[t, rens, cols] = size(SST);
16
17% Crear la referencia geográfica
18R = georefcells(latlim, lonlim, [rens, cols]);
19
20% Seleccionar la primera capa de datos
21sst = squeeze(SST(1,:,:));
22
23% Crear el mapa
24figure
25worldmap(latlim, lonlim) % Verificar el tipo de las variables
26geoshow(sst, R, 'DisplayType', 'surface')
27colorbar
28colormap('jet')
29title('TSM en la Bahía de La Paz')

Práctica #3

Tomando como referencia el script revisado en clase, utiliza los archivos jplG1SST_blap_2015.mat, sst_2018.mat y gadm41_MEX_2.shp para generar los siguientes mapas y gráficas:

  1. Mapa de TSM en la bahía de La Paz para un día específico (2015).

  2. Mapa de la TSM promedio en todo el año (2015).

  3. Gráfico Hovmöller, (pcolor) mostrando en el eje x las longitudes y en el eje y las fechas (meses).

  4. Gráfica (plot) mostrando la variación de la temperatura diaria promedio a lo largo del año (2015).

  5. Mapa de TSM en el sur del golfo de California para un día específico (2018).

  6. Mapa de TSM en la bahía de La Paz, recortada de los datos del 2018 (geocrop), para el mismo día elegido en el punto anterior.

  7. Gráfica para comparar la variación de la TSM diaria promedio a lo largo del 2015 y del 2018 en la bahía de La Paz. Para ello, primero se deben cortar (geocrop) las capas del 2018 según la cobertura de los datos del 2015.

  8. ¿Cuál es la cobertura, resolución temporal y espacial de los datos de TSM (2015 y 2018)? Incluye la respuesta en el script mencionando cómo obtuviste esa información.

  9. Documentar el script con los comentarios pertinentes de acuerdo con los datos y procedimientos utilizados.

  10. Como evidencia, adjuntar en Teams el script (xyz_03_tsm_blap.pdf). Debe incluir todos los mapas y las gráficas.

Nota

Los mapas deben incluir una escala de color y de distancia, un título apropiado con la fecha de los datos y la región del shapefile correspondiente al municipio o a los estados en la cobertura. Las gráficas deben incluir título y etiquetas en los ejes, y cuando se requiera, escala de color.

Representación de la estructura de los datos de TSM

SST Bahía de La Paz (2015)

4.3.2. Archivos .nc (NetCDF)

NetCDF es un formato de archivo ampliamente utilizado para almacenar datos científicos, especialmente en el ámbito de la oceanografía y la meteorología.

TSM en el sur del GC (2015)

Los archivos .nc se pueden leer utilizando la función ncread para acceder a las variables y ncinfo para obtener información sobre el contenido del archivo. La creación de la referencia geográfica se realiza con la función georefcells, como en el ejemplo anterior.

1% Importar datos de TSM
2SST = ncread('sst_2022.nc', 'analysed_sst');
3lat = ncread('sst_2022.nc', 'latitude');
4lon = ncread('sst_2022.nc', 'longitude');
5
6% Crear referencia geográfica
7
8% Crear mapa

Puesto que los comandos para visualizar los datos son similares en los ejemplos revisados, es conveniente implementar una función que permita crear este tipo de mapas. Por ejemplo:

Nota

Puedes agregar argumentos adicionales a la función para personalizar el mapa, como un shapefile, color de los polígonos, leyenda de escala de color, condicionantes para el rango, etc.

 1% Crea una mapa a partir de una capa raster
 2% data   - capa ráster
 3% r      - referencia geográfica
 4% rango  - límites de los valores [a, b]
 5% cmap   - escala de colores ('jet', 'turbo',etc)
 6% titulo - título del mapa
 7function [] = surfaceMap(data, r, rango, cmap, titulo)
 8    figure
 9    worldmap(r.LatitudeLimits, r.LongitudeLimits)
10    geoshow(zeros(size(data)), r, 'CData', data, 'DisplayType', 'surface')
11    title(titulo)
12    colormap(cmap)
13    colorbar;
14    clim(rango);

Cuando se utilizan datos en formato NetCDF, es común que se requiera especificar regiones para observar cambios en las variables en un periodo determinado. Por ejemplo, se podría utilizar un polígono irregular para delimitar una región de interés, luego aplicar una máscara a los datos ráster para mostrar solo la región seleccionada y proyectar esa selección en el cubo de datos para generar una serie de tiempo.

El siguiente script muestra un ejemplo de cómo crear una o varias máscaras binarias a partir de un archivo kmz o kml con polígono(s) dentro de la cobertura del ráster.

 1% Crea una o varias máscaras binaras (ráster) para la selección de ROIs
 2% Variables de entrada:
 3%   kmzName: nombre del archivo kml o kmz [polígono(s) dentro de la ROI (R)]
 4%   r: matriz de referencia geográfica del ráster al que se aplicará la máscara (tipo GeographicCellsReference)
 5% Regresa la matriz de la máscara:
 6%   mask: mascara binaria (1 | NaN)
 7% Si el archivo vectorial contiene n polígonos, regresa n mascaras [mask(:,:,n)]
 8
 9function mask = creaMask(kmzName, r)
10    gt = readgeotable(kmzName);
11    tb = geotable2table(gt,["lat","lon"]);
12    n = height(tb);
13    mask = zeros(r.RasterSize(1), r.RasterSize(2), n);
14    for j=1:n
15        [x, y] = geographicToIntrinsic(r, cell2mat(tb.lat(j)), cell2mat(tb.lon(j)));
16        mask(:,:,j) = double(poly2mask(x, y, r.RasterSize(1), r.RasterSize(2)));
17    end
18    mask(mask == 0) = NaN;
19end

Práctica #4

Descarga datos diarios de TSM para un año de tu área de estudio (NetCDF) con una cobertura aproximada de 4 x 5º. Si es necesario, también descarga un shapefile de esa región (similar al gadm41_MEX_1.shp). También digitaliza en Google Earth Pro al menos tres puntos y tres polígonos en la cobertura del archivo de TSM (guardarlos archivos kml o kmz, un archivo para los puntos y otro para los polígonos).

Tomando como referencia el script revisado en clase y los archivos mencionados en el párrafo anterior, generar los siguientes mapas y gráficas:

  1. Mapa de TSM en tu área de estudio para un día específico. Sobreponer al mapa los puntos y los polígonos digitalizados, mostrando las etiquetas para los puntos.

  2. Gráfica (plot) mostrando la variación de la temperatura diaria promedio a lo largo del año, en la ubicación de los polígonos digitalizados.

  3. Mapas que muestren la TSM (de un día específico) recortada con los polígonos digitalizados (generar un mapa por cada polígono).

  4. Gráfica (plot) que muestre la variación de la temperatura diaria promedio en la cobertura de los puntos digitalizados, a lo largo del año. Suavizar las series de tiempo (movmean).

  5. Mapas de los promedios mensuales de TSM en tu área de estudio (para todo el año). Los 12 mapas deben tener la misma escala de temperatura, en el título de cada mapa, incluir el mes y el año. Guardar estos mapas en formato png a 300 dpi.

  6. Documentar el script con los comentarios pertinentes de acuerdo con los datos y procedimientos utilizados.

  7. Como evidencia, adjuntar en Teams el script (xyz_04_tsm_aaaa_nomArea.pdf, aaaa = año al que corresponden los datos y nomArea = Nombre del área en una palabra). El archivo pdf debe incluir todos los mapas y las gráficas. Además, incluir los 12 mapas (png) de los promedios mensuales.

Nota

Los mapas deben incluir escala de color y de distancia, título apropiado con la fecha de los datos y el shapefile correspondiente al estado o región de la cobertura. Las gráficas deben incluir título y etiquetas en los ejes, y cuando se requiera, escala de color o leyenda.

⁂