28

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.

4

6 Answers 6

26
+50

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)
3
  • 3
    Woot works like a charm! Here is a cropping example from command line btw: gdal_translate --config AWS_REGION "some_region" --config AWS_ACCESS_KEY_ID "KEY_ID" --config AWS_SECRET_ACCESS_KEY "ACCESS_KEY" \ -srcwin 000 000 1000 1000 \ "/vsis3/bucket/file.ext" from_s3.tif
    – Skylion
    Jul 20, 2016 at 15:29
  • What do those values you hid look like? I think KEY_ID is a short text string, like a username. What is ACCESS_KEY? It seems like it is what is in a pem file but that is around 1000 characters, so it must be something else.
    – Solx
    Feb 7, 2017 at 21:31
  • Those will be just strings with numbers and letters kinda like a username and password. You can obtain those strings by setting IAM roles in AWS
    – RutgerH
    Aug 23, 2017 at 14:38
11

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.

4

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
3
  • While this is a useful answer, we already cache the metadata in a similar manner, but we want to know how to use VSI API so we can quickly crop small portions of large images.
    – Skylion
    Jul 19, 2016 at 14:43
  • I'm not certain if it's because the Mapzen API endpoint is a tiled WMS but the above code ran for me in under a minute, are you sure that the VSI API will be any faster?
    – clhenrick
    Jul 19, 2016 at 15:45
  • We are working with VERY large rasters and large raster datasets the bottleneck is definitely IO. Also the buckets we use are private and require credentials meaning using the S3 http API won't work in our case. It's not that we have to read every image, it's that we know we have to crap a small portion of a very large image.
    – Skylion
    Jul 19, 2016 at 16:01
1

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

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

1
  • I just added some clarifications to differentiate it from just using the S3 interface in Python.
    – Skylion
    Jul 18, 2016 at 18:10
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.