15

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.

16
+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)
  • 1
    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 '16 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 '17 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 '17 at 14:38
7

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.

1

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
  • 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 '16 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 '16 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 '16 at 16:01
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).

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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