AWS tiene un repositorio de datos abiertos. Tienes datos sobre observación de la tierra, meteorología, secuenciación del genoma, climatología, procesamiento del lenguaje natural,…
Estos datos están almacenados en Simple Storage Service (S3), que es un servicio de almacenamiento en la nube de Amazon Web Services (AWS).
En el artículo de hoy vamos a ver como podemos acceder a estos datos libres almacenados y distribuidos mediante S3 en AWS usando Python y la biblioteca s3fs.
Los recursos a los que queremos acceder se almacenan en algo que AWS llama buckets. Un bucket lo podríamos ver como una unidad de almacenamiento en la nube (como un disco duro) y la biblioteca s3fs nos ayudará precisamente a eso, a ver y tratar ese bucket como si fuera una unidad de disco (s3fs viene de S3 y de File System (o sistema de ficheros)).
Para el ejercicio de hoy he elegido los datos sobre teselas del terreno (Terrain Tiles). El bucket se encuentra en elevation-tiles-prod. Podéis ver un mapa de estos datos aquí.
¿Empezamos a trastear? Como siempre, vamos a empezar importando algunas cosas que vamos a usar:
|
1 2 3 4 5 6 7 8 |
from pathlib import Path import numpy as np import s3fs import xarray as xr import matplotlib.pyplot as plt import rasterio from pyproj import Transformer |
Vamos a ‘crear’ el sistema de ficheros. Como para acceder al bucket que vamos a usar no hace falta estar identificados lo haremos de forma anónima.
|
1 |
fs = s3fs.S3FileSystem(anon=True) |
Vamos a var qué cosas útiles podemos tener en el objeto fs (instancia de S3FileSystem):
|
1 2 3 |
for i in dir(fs): if not i.startswith('_') and callable(getattr(fs, i)): print(i) |
Lo anterior te debería de mostrar en pantalla algo parecido a lo siguiente:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 |
bulk_delete cat cat_file checksum chmod clear_instance_cache connect copy copy_basic copy_managed cp cp_file created current delete disk_usage download du end_transaction exists expand_path find from_json get get_delegated_s3pars get_file get_mapper get_tags getxattr glob head info invalidate_cache isdir isfile listdir ls makedir makedirs merge metadata mkdir mkdirs modified move mv object_version_info open pipe pipe_file put put_file put_tags read_block rename rm rm_file rmdir setxattr sign size split_path start_transaction stat tail to_json touch ukey upload url walk |
Algunas de las cosas que se ven más arriba son nombres que nos pueden sonar similar a herramientas que usamos desde nuestra terminal: ls, glob, makedir, touch,…
Vamos a hacer, por ejemplo, un ls del bucket, recuerda que el bucket se llamaba elevation-tiles-prod. Uso la palabra clave detail para obtener un poco más de información:
|
1 2 3 |
stuff = fs.ls('elevation-tiles-prod', detail=True) for s in stuff: print(f'{s["type"]: <10s}: {s["Key"]}') |
Lo anterior te debería mostrar en pantalla algo parecido a:
|
1 2 3 4 5 6 7 8 9 10 |
file : elevation-tiles-prod/index.html file : elevation-tiles-prod/main.js directory : elevation-tiles-prod/docs directory : elevation-tiles-prod/geotiff directory : elevation-tiles-prod/lib directory : elevation-tiles-prod/logs directory : elevation-tiles-prod/normal directory : elevation-tiles-prod/skadi directory : elevation-tiles-prod/terrarium directory : elevation-tiles-prod/v2 |
Vemos que hay varios ficheros y directorios. Para ir un directorio de los indicados más arriba simplemente ponemos la ruta completa. Voy a hacerlo usando el directorio “geotiff”:
|
1 2 3 |
more_stuff = fs.ls('elevation-tiles-prod/geotiff/') for s in more_stuff: print(s) |
Lo anterior debería mostrar algo parecido a:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
elevation-tiles-prod/geotiff/0 elevation-tiles-prod/geotiff/1 elevation-tiles-prod/geotiff/10 elevation-tiles-prod/geotiff/11 elevation-tiles-prod/geotiff/12 elevation-tiles-prod/geotiff/13 elevation-tiles-prod/geotiff/14 elevation-tiles-prod/geotiff/2 elevation-tiles-prod/geotiff/3 elevation-tiles-prod/geotiff/4 elevation-tiles-prod/geotiff/5 elevation-tiles-prod/geotiff/6 elevation-tiles-prod/geotiff/7 elevation-tiles-prod/geotiff/8 elevation-tiles-prod/geotiff/9 |
Creo que vais entendiendo como funciona esto. Voy a descargar unos cuantos GeoTiffs para ver cómo funciona lo de descargar la información que nos interese.
Los GeoTiffs están organizados en lo que se conoce como una pirámide de teselas. Podéis leer más sobre ello aquí. Básicamente, tenemos para cada nivel de zoom una carpeta para las ‘x’ y un fichero para cada ‘y’ dentro de esa ‘x’. Si quiero ir a la imagen que está en zoom=7, x=4 e y=10 la ruta será 7/4/10.tif.
En el listado resultado del código anterior vemos que tenemos rutas que acaban en números. Esos números corresponden al zoom o ‘z’.
Una vez leída la mini explicación anterior voy a descargar algunas imágenes a la máxima resolución (mayor zoom) en un directorio (local) que llamaré geotiffs.
Primero creo la carpeta geotiffs dentro del directorio desde donde estoy ejecutando el código:
|
1 2 |
p = Path('.', 'geotiffs') p.mkdir(mode=0o755, exist_ok=True) |
Y ahora las descargo. Para descargarlas usamos el método get con la ruta en remoto (lo que me quiero descargar) y la ruta en local (donde lo quiero descargar):
|
1 2 3 4 5 6 7 8 9 |
z = 14 x = range(11670, 11679) y = range(6440, 6447) for xi in x: for yi in y: remote = f'elevation-tiles-prod/geotiff/{z}/{xi}/{yi}.tif' local = str(p / f'zoom{z}_x{xi}_y{yi}.tif') fs.get(remote, local) |
Y, si no ha fallado nada, tendremos los ficheros descargados en nuestro PC.
Y, ya que hemos llegado hasta aquí, vamos a hacer algo con la información descargada. La información descargada muestra la elevación alrededor del K2. En este caso, este artículo es un triste y humilde homenaje a un alpinista que acaba de perder la vida en la montaña salvaje, Sergi Mingote.
Pintamos la información usando xarray y combinando todas las teselas que hemos descargado:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
files = list(p.glob('*')) ds_list = [] crs_dest = None for f in files: if crs_dest is None: _ = xr.open_rasterio(f) crs_dest = _.crs _ = _.to_dataset(name='elev') else: _ = xr.open_rasterio(f).to_dataset(name='elev') ds_list.append(_) combined = xr.combine_by_coords(ds_list) |
El gráfico con la información (además de la localización del K2 en la imagen) lo podemos obtener usando el código siguiente:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
fig, ax = plt.subplots(figsize=(18, 14)) ax.contourf( combined.x.values, combined.y.values, combined.elev.values[0,...], levels=np.arange(2000, 9500, 100), cmap=plt.cm.Greens_r ) ax.contour( combined.x.values, combined.y.values, combined.elev.values[0,...], levels=np.arange(2000, 9500, 500), colors='gray', linewidths=1 ) lat = 35.88111 lon = 76.51333 t = Transformer.from_crs('epsg:4326', crs_dest[6:]) x, y = t.transform(lat, lon) ax.scatter([x], [y], marker='^', c='k', s=300) ax.text(x + 2000, y, 'Mt. K2', fontsize=20, horizontalalignment='right') fig.savefig('terrain_tiles_k2.png') |
El anterior código os debería haber creado una figura llamada “terrain_tiles_k2.png” en la misma carpeta donde estés ejecutando el código.
El resultado final es:
Sed felices y mis condolencias a la familia y amigos de Sergi Mingote.
El instructivo es maravilloso y súper práctico. Pero lo mejor de todo es el homenaje a Sergi Mingote.
Te mando un abrazo y muchas gracias
Gracias a ti por tu comentario.