How to efficiently access files with GDAL from an S3 bucket using VSIS3?

Question

So, GDAL has recently added a new feature that allows random reading of S3 bucket files. I am looking to crop GDAL images from multiple tiles of an image without having to download the whole file. I've only seen very sparse documentation on how to configure and access an S3 bucket though GDAL and am a little confused on how to begin? Would someone be kind enough to provide an extremely short example/tutorial on how one would go about setting the virtual filesystem for GDAL in order to accomplish this goal? Bonus pts if your solution allows it to be scripted via Python!

To clarify: We already have it done in Python. The issue with Python is that you have to download the whole image to operate it with it. The newest version of GDAL has support for mounting the S3 bucket so that if we need to say a crop a small portion of the image, we can operate directly on that smaller portion. Alas, as the feature only was released on the stable branch in January, I haven't found any documentation on it. So the solution should use the VSI3 system in the newest release of GDAL or otherwise smartly uses the system to prevent the user from needing to download the entire image to an EBS drive to operate on it.

That is to say the bounty will be awarded to answer that uses the VSI APIs found in the newest versions of GDAL so that the whole file does not need to be read into memory or disk. Also, we the buckets we use are not always public so many of the HTTP tricks being posted won't work in many of our situations.

Answer 1

I've found when something isn't particularly well documented in GDAL, that looking through their tests can be useful.

The /vsis3 test module has some simple examples, though it doesn't have any examples of actually reading chunks.

I've cobbled together the code below based on the test module, but I'm unable to test as GDAL /vsis3 requires credentials and I don't have an AWS account.

"""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

Since /vsis3/ is implemented in GDAL you can also use rasterio to read Windows of S3 datasets. This requires either your credentials to be set up for boto or using rasterios AWS session handler.

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.

See also rasterios windowed-rw and VSI docs.

Answer 3

Try using an XML file to store the WMS info in, more details are at the GDAL WMS documentation.

Here's an example WMS XML file to retrieve data from Mapzen's Elevation API:

<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>

You can then clip to a bounding box like so:

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

Answer 4

To read from S3, I used '/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

I don't know much about S3 buckets but it seems that it's a cloud storage drive with authentication using http REST services. i.e. could be used as an ordinary mounting point, with an associated uri.

If you are looking for cropping parts of images/raster then the file needs to be in an appropriated format.

Take a look at the TMS specification http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification

(Perhaps netCDF could also do the trick.)

GDAL also reads and writes TMS formats. Basically it's only a standard directory structure with some metadata files.

Now, the trick is to create on the fly the url with the geographic extent parameters through the TMS driver.

Take a look at the OpenLayers TMS driver documentation: http://dev.openlayers.org/docs/files/OpenLayers/Layer/TMS-js.html To see how it handles the requests based on location, zoom and extents.

Of course it can be done in Python. You need first to create the appropriate "mounting-point" (or path) URI with viscurl (according to the documentation) and then, once is mounted go to the specific tile according to the TMS specification (which is an extension of the path).

Answer 6

In my case I tried without any success all suggestions here :'(. What ended up working for me:

from pathlib import Path

from uuid import uuid4
GDAL_DATASET = gdal.Dataset
import boto3
from pydantic import BaseSettings
from osgeo import gdal

class S3Configuration(BaseSettings):
    """
    S3 configuration class
    """
    s3_access_key_id: str = ''
    s3_secret_access_key: str = ''
    s3_region_name: str = ''
    s3_endpoint_url: str = ''
    s3_bucket_name: str = ''
    s3_use: bool = False

S3_CONF = S3Configuration()
S3_STR = 's3'
S3_SESSION = boto3.session.Session()
S3 = S3_SESSION.resource(
    service_name=S3_STR,
    aws_access_key_id=S3_CONF.s3_access_key_id,
    aws_secret_access_key=S3_CONF.s3_secret_access_key,
    endpoint_url=S3_CONF.s3_endpoint_url,
    region_name=S3_CONF.s3_region_name,
    use_ssl=True,
    verify=True,
)
BUCKET = S3_CONF.s3_bucket_name
ZIP_EXT = '.zip'


def get_gdal_data(file_path: Path, s3_use: S3_CONF.s3_use) -> GDAL_DATASET:
    """
    Retrieves the tif content with gdal from the passed file_path. Do so either locally or from S3
    """
    if s3_use:
        return get_in_memory_tile(get_s3_object(file_path))
    return get_tif_tile(file_path)


def get_s3_object(file_path: Path) -> bytes:
    """
    Retrieve as bytes the content associated to the passed file_path
    """
    return S3.Object(bucket_name=BUCKET, key=forge_key(file_path)).get()['Body'].read()


def forge_key(file_path: Path) -> str:
    """
    Edit this code at your convenience to forge the bucket key out of the passed file_path
    """
    return str(file_path.relative_to(*file_path.parts[:2]))


def get_tif_tile(file_path):
    """
    Retrieves tif content with gdal from local file
    """
    return gdal.Open(str(file_path))

def get_in_memory_tile(in_memory_data: bytes):
    """
     Retrieves tif content with gdal from vsimem
    """
    filename = _in_memory_filename()
    gdal.FileFromMemBuffer(filename, in_memory_data)
    tile = gdal.Open(filename)
    unlinked = gdal.Unlink(filename)
    return tile

def _in_memory_filename():
    """
    Generate a random unique name to avoid concurrent accesses if multiple reading from vsimem
    """
    return f'/vsimem/{get_new_uuid()}'


def get_new_uuid() -> str:
    """
    Generate uuid4 strings
    """
    return str(uuid4())

There is quite a lot of code, but the first part is very helpful to easily use a minio proxy for local devs (just have to change the .env).

The key to solve the issue for me was the use of gdal not so well documented (in my opinion) but life saver (in my case :)) FileFromMemBuffer in conjunction with vsimem