Saltar al contenido

El lento vagar del iceberg A-68 visualizado con cartopy en Python

Etiquetas:

La barrera de hielo Larsen se encuentra en el mar de Weddell a lo largo de la península antártica. Esta barrera se divide a su vez en varias barreras menores de hielo que han ido sufriendo numerosas variaciones. En julio de 2017 se produjo una de estas variaciones, un gran cambio muy destacable, que tuvo lugar en la barrera menor llamada Larsen C.

¿Qué pasó? se desprendió gran parte del hielo de Larsen C dando lugar al iceberg A-68. Este iceberg es uno de los más grandes de los que conocemos.

El tamaño inicial de este iceberg fue de 5800 km², más grande que todas las islas Baleares juntas. Mallorca, la isla mayor, tiene un área de unos 3640 km².

Los nombres de los icebergs antárticos los pone el US ice center. Este centro monitoriza los icebergs a nivel global y de la Antártida y nombra a los icebergs antárticos que tienen un tamaño superior a 10 millas náuticas de la siguiente forma:

  • La letra inicial proviene del lugar donde se origina el iceberg. A es porque se ha originado entre longitud 0º y longitud 90º oeste.
  • El número es un número que se va incrementando con cada nuevo iceberg en cada uno de los orígenes.

Este iceberg lleva desde entonces vagando por el hemisferio sur y ahora se encuentra muy cerca de las islas Georgias del sur las cuales se encuentran a más de 1500 km de distancia de donde se originó el iceberg A68.

Si algún oso polar se quedó atrapado en el iceberg lleva más de tres años vagando por el mundo.

Pero este iceberg se va despedazando poco a poco. Cada trozo se renombra con una letra que se añade al final del nombre del iceberg. Hasta ahora se ha despedazado en cuatro trozos: A-68A, A-68B, A-68C, A-68D.

Vamos a ver un poco de su historia usando cartopy y Python.

Importamos unas pocas bibliotecas:

import datetime as dt
from pathlib import Path

import requests
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import get_named_colors_mapping
import cartopy.crs as ccrs
plt.style.use('bmh')

Lo siguiente que necesito es una función que me descargue datos del US ice center. Te puedes descargar datos semanales con información de los icebergs en formato csv:

def download_locations(start: dt.date, end: dt.date) -> None:
    # create path to save all weekly data
    p = Path('icebergs_path')
    p.mkdir(0o755, exist_ok=True)
    date = start
    while date <= end:
        year = date.year
        datestr = date.strftime('%Y%m%d')
        url = (
            'https://usicecenter.gov/File/DownloadProduct'
            '?products=%2Ficeberg%2F'
            f'{year}&fName=AntarcticIcebergs_{datestr}.csv'
        )
        r = requests.get(url)
        with open(p / f'icebergs_{datestr}.csv', 'w') as f:
            f.write(r.content.decode('utf-8'))
        print(datestr)
        date += dt.timedelta(days=7)

No hay mucho misterio en la función anterior, genera una serie de urls a partir de fechas y descarga los datos que guarda en una carpeta que se llama ‘icebergs_path’ y que se creará en el directorio desde donde ejecutes el código.

Ahora llamo a la función anterior con fechas desde mediados de julio del 2017 hasta la actualidad (la actualidad es la fecha de este artículo):

# Descargamos los datos de icebergs del US ice center        
download_locations(dt.date(2017, 7, 14), dt.date.today())

Con todos los ficheros que nos hemos descargado vamos a crear un dataframe con la información relacionada solo con el iceberg A-68, el resto lo desechamos.

# Creamos un DF con la info interesante
files = sorted(list(Path('.', 'icebergs_path').glob('*.csv')))

columns = (
    'Iceberg,Length (NM),Width (NM),Latitude,Longitude,Remarks,Last Update'
).split(',')

def to_date(datestr: str):
    m, d, y = map(int, datestr.split('/'))
    return dt.date(y, m, d)

df = pd.DataFrame(
    columns=columns
)

for f in files:
    with open(f) as fi:
        for line in fi:
            if line.startswith('A68'):
                values = line[:-1].split(',')
                dd = {k:v for k, v in zip(columns, values)}
                df = df.append(dd, ignore_index=True)
df['Latitude'] = df.loc[:, 'Latitude'].astype(np.float)
df['Longitude'] = df.loc[:, 'Longitude'].astype(np.float)
df['Length (NM)'] = df.loc[:, 'Length (NM)'].astype(np.float)
df['Width (NM)'] = df.loc[:, 'Width (NM)'].astype(np.float)
df['Last Update'] = df.loc[:, 'Last Update'].apply(to_date)

Nuevamente, el código no tiene mucha historia, leer ficheros en formato csv, modificar algunas cosas para tener los datos más usables, etc.

En el nuevo dataframe que hemos creado podemos agrupar los datos en función del iceberg o “sub-iceberg” (A-68, A-68A, A-68B,….):

# Creamos grupos con la info de cada una de las fases del iceberg A68
groups = df.groupby('Iceberg')
icebergs = list(groups.groups.keys())
colors = list(get_named_colors_mapping().values())[:len(icebergs)]

Vamos a ver la evolución de la longitud a lo largo del tiempo:

fig, ax = plt.subplots(figsize=(20, 10))
for i, g in enumerate(groups):
    x = g[1]['Last Update'].values
    y = g[1]['Length (NM)'].values
    print(x)
    ax.plot_date(x, y, label=icebergs[i], color=colors[i], lw=3)
ax.legend()
# Lo siguiente está comentado porque depende 
# de si solo quieres ver la figura o guardarla
#fig.show()
#fig.savefig('un_nombre.png')

El anterior código debería mostraros un resultado como el siguiente:

El iceberg A-68 se desgajo en el A-68A (el más largo) y el A-68B y en los últimos tiempos el A-68A ha perdido parte porque se ha separado el trozo A-68D cerca de las islas Georgias del sur.

Si yo fuese un habitante de esas islas (viven solo unas decenas de personas) y viese algo tan grande (alrededor de 1,000,000,000,000 de toneladas de hielo) como la isla que habito dirigiéndose hacia la misma estaría bastante intranquilo… 😮

Veamos en un mapa con cartopy cómo se ha ido desplazando esta brutalidad natural:

# Dibujamos la información de cada uno de los grupos
fig, ax = plt.subplots(
    figsize=(10, 10), 
    subplot_kw={'projection': ccrs.PlateCarree()}
)
ax.stock_img()
ax.coastlines(resolution='50m', color='black', linewidth=1)

xmin = df.Longitude.min() - 5
xmax = df.Longitude.max() + 5
ymin = df.Latitude.min() - 5
ymax = df.Latitude.max() + 5
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

for i, g in enumerate(groups):
    y = g[1].Latitude
    x = g[1].Longitude
    sizes = g[1]['Width (NM)'] * g[1]['Length (NM)'] / 10
    ax.plot(x, y, label=icebergs[i], color=colors[i], lw=3)
    ax.scatter(x, y, color=colors[i], s=sizes)
ax.legend()
# Lo siguiente está comentado porque depende 
# de si solo quieres ver la figura o guardarla
#fig.show()
#fig.savefig('otro_nombre.png')

El resultado de lo anterior sería:

En lo anterior, el código más novedoso relacionado con cartopy sería el de definir la proyección. En este caso uso la proyección ‘Plate carrée’ (un caso especial de la proyección cilíndrica equidistante). Luego uso dos métodos especiales del GeoAxesSubplot de cartopy:

  • stock_img que me permite añadir imágenes geográficas de fondo.
  • coastlines que me permite añadir las líneas de costa a distinta resolución.

Los puntos de cada una de las trayectorias van en relación con el tamaño de cada uno de los ‘trozos’ del iceberg. Si os fijáis, en el código he calculado sizes como si fuese un rectángulo de Length y Width. Esas sizes se usan para dibujar el tamaño del símbolo en el scatter. No es muy científico, simplemente para resaltar un poco el tamaño del iceberg en cada fase y mostrar que el A-68A es el más grande.

Una visualización de cuando el iceberg se separá del banco Larsen C se muestra a continuación:

Animación obtenida de Wikipedia.

Aquí tenéis otra animación mucho más completa.

El código completo lo puedes ver a continuación:

import datetime as dt
from pathlib import Path

import requests
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import get_named_colors_mapping
import cartopy.crs as ccrs
plt.style.use('bmh')

def download_locations(start: dt.date, end: dt.date) -> None:
    # create path to save all weekly data
    p = Path('icebergs_path')
    p.mkdir(0o755, exist_ok=True)
    date = start
    while date <= end:
        year = date.year
        datestr = date.strftime('%Y%m%d')
        url = (
            'https://usicecenter.gov/File/DownloadProduct'
            '?products=%2Ficeberg%2F'
            f'{year}&fName=AntarcticIcebergs_{datestr}.csv'
        )
        r = requests.get(url)
        with open(p / f'icebergs_{datestr}.csv', 'w') as f:
            f.write(r.content.decode('utf-8'))
        print(datestr)
        date += dt.timedelta(days=7)

# Descargamos los datos de icebergs del US ice center        
download_locations(dt.date(2017, 7, 14), dt.date.today())

# Creamos un DF con la info interesante
files = sorted(list(Path('.', 'icebergs_path').glob('*.csv')))

columns = (
    'Iceberg,Length (NM),Width (NM),Latitude,Longitude,Remarks,Last Update'
).split(',')

def to_date(datestr: str):
    m, d, y = map(int, datestr.split('/'))
    return dt.date(y, m, d)

df = pd.DataFrame(
    columns=columns
)

for f in files:
    with open(f) as fi:
        for line in fi:
            if line.startswith('A68'):
                values = line[:-1].split(',')
                dd = {k:v for k, v in zip(columns, values)}
                df = df.append(dd, ignore_index=True)
df['Latitude'] = df.loc[:, 'Latitude'].astype(np.float)
df['Longitude'] = df.loc[:, 'Longitude'].astype(np.float)
df['Length (NM)'] = df.loc[:, 'Length (NM)'].astype(np.float)
df['Width (NM)'] = df.loc[:, 'Width (NM)'].astype(np.float)
df['Last Update'] = df.loc[:, 'Last Update'].apply(to_date)

# Creamos grupos con la info de cada una de las fases del iceberg A68
groups = df.groupby('Iceberg')
icebergs = list(groups.groups.keys())
colors = list(get_named_colors_mapping().values())[:len(icebergs)]

fig, ax = plt.subplots(figsize=(20, 10))
for i, g in enumerate(groups):
    x = g[1]['Last Update'].values
    y = g[1]['Length (NM)'].values
    print(x)
    ax.plot_date(x, y, label=icebergs[i], color=colors[i], lw=3)
ax.legend()
# Lo siguiente está comentado porque depende 
# de si solo quieres ver la figura o guardarla
#fig.show()
#fig.savefig('un_nombre.png')

# Dibujamos la información de cada uno de los grupos
fig, ax = plt.subplots(
    figsize=(10, 10), 
    subplot_kw={'projection': ccrs.PlateCarree()}
)
ax.stock_img()
ax.coastlines(resolution='50m', color='black', linewidth=1)

xmin = df.Longitude.min() - 5
xmax = df.Longitude.max() + 5
ymin = df.Latitude.min() - 5
ymax = df.Latitude.max() + 5
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

for i, g in enumerate(groups):
    y = g[1].Latitude
    x = g[1].Longitude
    sizes = g[1]['Width (NM)'] * g[1]['Length (NM)'] / 10
    ax.plot(x, y, label=icebergs[i], color=colors[i], lw=3)
    ax.scatter(x, y, color=colors[i], s=sizes)
ax.legend()
# Lo siguiente está comentado porque depende 
# de si solo quieres ver la figura o guardarla
#fig.show()
#fig.savefig('otro_nombre.png')

(continuará…)

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

− five = four

Pybonacci