¿Cómo acceder eficientemente a archivos con GDAL desde un bucket S3 usando VSIS3?

Question

Así, GDAL ha añadido recientemente un nueva función que permite la lectura aleatoria de archivos de cubo S3. Estoy buscando para recortar imágenes GDAL de múltiples azulejos de una imagen sin tener que descargar todo el archivo. Sólo he visto documentación muy escasa sobre cómo configurar y acceder a un cubo de S3 a través de GDAL y estoy un poco confundido sobre cómo empezar? ¿Sería alguien tan amable de proporcionar un ejemplo/tutorial extremadamente corto sobre cómo uno iría sobre la configuración del sistema de archivos virtual para GDAL con el fin de lograr este objetivo? ¡Pts de bonificación si su solución permite que sea de secuencias de comandos a través de Python!

Para aclarar: ya lo tenemos hecho en Python. El problema con Python es que hay que descargar toda la imagen para poder operar con ella. La versión más reciente de GDAL tiene soporte para montar el cubo S3 de modo que si necesitamos, por ejemplo, recortar una pequeña porción de la imagen, podemos operar directamente sobre esa pequeña porción. Por desgracia, como la característica sólo se publicó en la rama estable en enero, no he encontrado ninguna documentación al respecto. Así que la solución debe utilizar el sistema VSI3 en la versión más reciente de GDAL o de otra manera inteligente utiliza el sistema para evitar que el usuario tenga que descargar toda la imagen a una unidad de EBS para operar en él.

Es decir, la recompensa se otorgará a la respuesta que utilice las API de VSI que se encuentran en las versiones más recientes de GDAL para que no sea necesario leer todo el archivo en memoria o en disco. Además, los buckets que utilizamos no siempre son públicos, por lo que muchos de los trucos HTTP que se publican no funcionarán en muchas de nuestras situaciones.

Answer 1

He descubierto que cuando algo no está particularmente bien documentado en GDAL, puede ser útil echar un vistazo a sus pruebas.

En /vsis3 módulo de prueba tiene algunos ejemplos sencillos, aunque no tiene ningún ejemplo de lectura real de trozos.

He improvisado el siguiente código basado en el módulo de prueba, pero no puedo probar como GDAL /vsis3 requiere credenciales y no tengo una cuenta de AWS.

"""This should read from the Sentinal-2 public dataset
   More info - http://sentinel-pds.s3-website.eu-central-1.amazonaws.com"""

from osgeo import gdal
import numpy as np

# These only need to be set if they're not already in the environment,
# ~/.aws/config, or you're running on an EC2 instance with an IAM role.
gdal.SetConfigOption('AWS_REGION', 'eu-central-1')
gdal.SetConfigOption('AWS_SECRET_ACCESS_KEY', 'MY_AWS_SECRET_ACCESS_KEY')
gdal.SetConfigOption('AWS_ACCESS_KEY_ID', 'MY_AWS_ACCESS_KEY_ID')
gdal.SetConfigOption('AWS_SESSION_TOKEN', 'MY_AWS_SESSION_TOKEN')

# 'sentinel-pds' is the S3 bucket name
path = '/vsis3/sentinel-pds/tiles/10/S/DG/2015/12/7/0/B01.jp2'
ds = gdal.Open(path)

band = ds.GetRasterBand(1)

xoff, yoff, xcount, ycount = (0, 0, 10, 10)
np_array = band.ReadAsArray(xoff, yoff, xcount, ycount)

Answer 2

Desde /vsis3/ está implementado en GDAL también puede utilizar rasterio para leer ventanas de conjuntos de datos S3. Esto requiere que sus credenciales estén configuradas para boto o utilizando rasterios Gestor de sesiones de AWS .

import rasterio

with rasterio.open('s3://landsat-pds/L8/139/045/LC81390452014295LGN00/LC81390452014295LGN00_B1.TIF') as ds:
    window = ds.read(window=((0, 100), (0, 100)))  # read a 100 by 100 window in the upper left corner.

Véase también rasterios windowed-rw et VSI docs.

Answer 3

Pruebe a utilizar un archivo XML para almacenar la información de WMS, más detalles en el Documentación del SGA GDAL .

A continuación se muestra un archivo XML WMS de ejemplo para recuperar datos de la Elevation API de Mapzen:

<GDAL_WMS>
  <Service name="TMS">
    <ServerUrl>https://s3.amazonaws.com/elevation-tiles-prod/geotiff/${z}/${x}/${y}.tif</ServerUrl>
  </Service>
  <DataWindow>
    <UpperLeftX>-20037508.34</UpperLeftX>
    <UpperLeftY>20037508.34</UpperLeftY>
    <LowerRightX>20037508.34</LowerRightX>
    <LowerRightY>-20037508.34</LowerRightY>
    <TileLevel>14</TileLevel>
    <TileCountX>1</TileCountX>
    <TileCountY>1</TileCountY>
    <YOrigin>top</YOrigin>
  </DataWindow>
  <Projection>EPSG:3857</Projection>
  <BlockSizeX>512</BlockSizeX>
  <BlockSizeY>512</BlockSizeY>
  <BandsCount>1</BandsCount>
  <DataType>Int16</DataType>
  <ZeroBlockHttpCodes>403,404</ZeroBlockHttpCodes>
  <DataValues>
    <NoData>-32768</NoData>
  </DataValues>
  <Cache/>
</GDAL_WMS>

A continuación, puede recortar a un cuadro delimitador de la siguiente manera:

gdalwarp -of "GTiff" -te -13648825.0817 4552130.7825 -13627575.5878 4565507.2624 mapzen_wms.xml test.tif

Answer 4

Para leer desde S3, utilicé '/vsis3_streaming/ (AWS S3 files: streaming)'.

( https://gdal.org/user/virtual_file_systems.html#vsis3-streaming-aws-s3-files-streaming ).

# first set s3 client
session = boto3.Session(profile_name='default')
s3_client = session.client('s3')
s3_bucket_name = 'my-data-lake'
s3 = boto3.resource('s3')
my_bucket  = s3.Bucket(s3_bucket_name)

# reading from vsis3_streaming
bucket_tif = 'sandbox/user01/field00/myfile.tif'
ds = gdal.Open(str('vsis3_streaming/' + s3_bucket_name + '/' + bucket_tif))

Answer 5

No sé mucho sobre cubos S3 pero parece que es una unidad de almacenamiento en la nube con autenticación mediante servicios http REST. Es decir, se podría usar como un punto de montaje normal y corriente, con una uri asociada.

Si desea recortar partes de imágenes/raster, el archivo debe estar en un formato adecuado.

Consulte el pliego de condiciones de TMS http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification

(Quizá netCDF también podría servir).

GDAL también lee y escribe formatos TMS. Básicamente es sólo una estructura de directorios estándar con algunos archivos de metadatos.

Ahora, el truco consiste en crear sobre la marcha la url con los parámetros de extensión geográfica a través del controlador TMS.

Eche un vistazo a la documentación del controlador OpenLayers TMS: http://dev.openlayers.org/docs/files/OpenLayers/Layer/TMS-js.html Para ver cómo gestiona las solicitudes en función de la ubicación, el zoom y la extensión.

Por supuesto que se puede hacer en Python. En primer lugar, es necesario crear el "punto de montaje" apropiado (o ruta) URI con viscurl (de acuerdo con la documentación) y luego, una vez que se monta ir a la baldosa específica de acuerdo con la especificación TMS (que es una extensión de la ruta).