terracotta icon indicating copy to clipboard operation
terracotta copied to clipboard

Add support for seamless layers

Open dionhaefner opened this issue 6 years ago • 34 comments

dionhaefner avatar May 29 '18 17:05 dionhaefner

The notable difference to non-seamless is that the dataset to look up needs to be determined from xyz bounds, so keys alone are not enough to identify a dataset. We would also need to implement merging from several data files to one tile for low zoom levels. However, this comes with its own challenges: For the lowest zoom level, we would have to read pixels from every single file in the dataset, so we might have to enforce a minimum zoom level.

A workaround would be to require the user to create a single (possibly gigantic) GeoTiff, or to add all data files as separate layers and display them on top of each other.

dionhaefner avatar Jul 02 '18 08:07 dionhaefner

I just learned about terracotta (really awesome tool!) and this feature would be great. I frequently have a collection of rasters that are tiled, and would like to serve them as a single layer without merging them into a gigantic single GeoTiff.

mccarthyryanc avatar Mar 14 '19 22:03 mccarthyryanc

Maybe the best solution to larger-than-GeoTIFF layers is to store the data not in GeoTIFF tiles (e.g. Sentinel 2 tiles) but in even smaller chunks (say 256x256) in a fast database, like GEE seems to be doing it. But that might be beyond the scope of Terracotta - would be a hell of a driver...

j08lue avatar Mar 15 '19 08:03 j08lue

I've fallen out of love with seamless layer support because there are so many micro-decisions to make for relatively little gain (mainly how to handle overlap). But we could look into supporting GDAL-style VRTs instead, and offload all the heavy lifting to GDAL / rasterio. That way we'd hopefully only need very little additional driver code.

dionhaefner avatar Mar 15 '19 10:03 dionhaefner

mainly how to handle overlap

If that is the main concern, we could require that there is no overlap (or undefined behavior in such case). But there are probably more?

supporting GDAL-style VRTs

VRTs over S3 files? - But why not? As far as I know, rasterio reads VRTs fine and might even be able to create them (the XML representation) by now.

j08lue avatar Mar 15 '19 12:03 j08lue

If that is the main concern, we could require that there is no overlap (or undefined behavior in such case). But there are probably more?

I'm mostly concerned about numerical roundoff and off-by-one errors. I'd expect it to take a lot of fiddling to get right.

dionhaefner avatar Mar 15 '19 15:03 dionhaefner

There is a vrt module in rasterio, for in-memory vrt.

mccarthyryanc avatar Mar 15 '19 15:03 mccarthyryanc

Yes, we use those extensively. I was thinking of the GDAL VRTs in XML form that allow you to combine datasets transparently ahead of time.

dionhaefner avatar Mar 15 '19 15:03 dionhaefner

@dionhaefner , ah I see. I've created VRT's with gdalbuildvrt and specified rasters on S3 making sure relativeToVRT="0" is set and that the source filename includes /vsis3/.... Do you think that would work with serving:

terracotta serve -r /{name}/{date}_{band}_{}.vrt

mccarthyryanc avatar Mar 15 '19 16:03 mccarthyryanc

It might actually work by accident if rasterio exposes the VRT exactly like a regular raster file. Could you try and let us know?

dionhaefner avatar Mar 15 '19 17:03 dionhaefner

Okay, gave it a shot, but it didn't work. Here was my setup. I already had rasters on S3 that were optimized via terracotta optimize-rasters.

Build list to build vrt:

aws s3 ls s3://bucketname/prefix/ | \
    rev | cut -d\  -f1 | rev | \
    sed 's|^|/vsis3/bucketname/prefix/|g' \
    > 20190315_testing_filelist

Build the VRT:

gdalbuildvrt -input_file_list 20190315_testing_filelist 20190315_testing.vrt

Serve the vrt (from EC2 instance):

terracotta serve --allow-all-ips --port 8787 -r ./{date}_{name}.vrt

Connect from local machine:

terracotta connect INSTANCE-IP:8787

At this point I get the preview app open within a browser, and I can see the dataset listed. When I put the cursor over that row, I get the correct outlined area on the map (plus the low-res sample image). But when I select it nothing is displayed, however I do get the correct limits set on the "Adjust contrast" slider bar. When I pan/zoom it looks like the server is getting requests, just not displaying anything.

rio info gives valid info on both the VRT and any of the individual raster tiles. Hopefully this help, happy to run other tests.

Server output:


Ingesting raster files: 100%|███████████████████████████| 1/1 [00:00<00:00, 3407.23it/s]

 * Serving Flask app "terracotta.server" (lazy loading)
 * Environment: production
   WARNING: Do not use the development server in a production environment.
   Use a production WSGI server instead.
 * Debug mode: off
 * Running on http://0.0.0.0:8787/ (Press CTRL+C to quit)
LOCALIP - - [15/Mar/2019 22:35:25] "GET /keys HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:25] "GET /swagger.json HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /colormap?colormap=viridis&stretch_range=[0,1]&num_values=100 HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /colormap?colormap=greys_r&stretch_range=[0,1]&num_values=100 HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /colormap?colormap=rdbu_r&stretch_range=[0,1]&num_values=100 HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /colormap?colormap=ylgn&stretch_range=[0,1]&num_values=100 HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /colormap?colormap=bugn&stretch_range=[0,1]&num_values=100 HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /colormap?colormap=magma&stretch_range=[0,1]&num_values=100 HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /colormap?colormap=gist_earth&stretch_range=[0,1]&num_values=100 HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /colormap?colormap=ocean&stretch_range=[0,1]&num_values=100 HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /keys HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:35:28] "GET /datasets?limit=16&page=0 HTTP/1.1" 200 -
 [!] /home/ubuntu/src/terracotta/terracotta/drivers/sqlite.py:318: PerformanceWarning: Raster file /tmp/20190315_testing.vrt is not a valid cloud-optimized GeoTIFF. Any interaction with it will be significantly slower. Consider optimizing it through `terracotta optimize-rasters` before ingestion.
  metadata = self.compute_metadata(filepath[keys], max_shape=self.LAZY_LOADING_MAX_SHAPE)

LOCALIP - - [15/Mar/2019 22:37:10] "GET /metadata/20190315/testing HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/preview.png?tile_size=[128,128] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/0/2.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/3/0.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/3/1.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/0/3.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/0/0.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/0/1.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/1/1.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/2/1.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/2/1/2.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/4/5/4.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/4/6/6.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/4/4/6.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/4/3/6.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/4/7/6.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/5/10/11.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/5/7/12.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/5/11/13.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/5/8/11.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:07] "GET /singleband/20190315/testing/5/9/12.png?colormap=greys_r&stretch_range=[0,1] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:45:23] "GET /singleband/20190315/testing/10/269/407.png?colormap=greys_r&stretch_range=[172.37033081054688,320.1650085449219] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:47:37] "GET /singleband/20190315/testing/10/268/407.png?colormap=greys_r&stretch_range=[172.37033081054688,320.1650085449219] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:47:37] "GET /singleband/20190315/testing/10/269/406.png?colormap=greys_r&stretch_range=[172.37033081054688,320.1650085449219] HTTP/1.1" 200 -
LOCALIP - - [15/Mar/2019 22:47:37] "GET /singleband/20190315/testing/10/270/407.png?colormap=greys_r&stretch_range=[172.37033081054688,320.1650085449219] HTTP/1.1" 200 -

Connection output:

 * Serving Flask app "terracotta.client" (lazy loading)
 * Environment: production
   WARNING: Do not use the development server in a production environment.
   Use a production WSGI server instead.
 * Debug mode: off
 * Running on http://127.0.0.1:5100/ (Press CTRL+C to quit)
[12825:12825:0315/153527.026940:ERROR:sandbox_linux.cc(364)] InitializeSandbox() called with multiple threads in process gpu-process.
127.0.0.1 - - [15/Mar/2019 15:35:27] "GET / HTTP/1.1" 200 -
[12751:12806:0315/153527.050398:ERROR:browser_process_sub_thread.cc(209)] Waited 14 ms for network service
Opening in existing browser session.

mccarthyryanc avatar Mar 15 '19 22:03 mccarthyryanc

@mccarthyryanc can you post here the VRT file you created?

j08lue avatar Mar 16 '19 07:03 j08lue

At this point I get the preview app open within a browser, and I can see the dataset listed. When I put the cursor over that row, I get the correct outlined area on the map (plus the low-res sample image). But when I select it nothing is displayed, however I do get the correct limits set on the "Adjust contrast" slider bar. When I pan/zoom it looks like the server is getting requests, just not displaying anything.

That is very reassuring. If the preview image works, tile retrieval should work, too. And correct limits and footprint means the metadata computation works, too. I'm just a bit confused about all the &stretch_range=[0,1] in your queries. Could be the reason why you don't see anything.

dionhaefner avatar Mar 16 '19 11:03 dionhaefner

@j08lue I can't provide the exact VRT, contains private info. But here is an example VRT that is pretty much the same. It was built with plain gdalbuildvrt:

<VRTDataset rasterXSize="61000" rasterYSize="56000">
  <SRS>PROJCS["NAD83(2011) / UTM zone 16N",GEOGCS["NAD83(2011)",DATUM["NAD83_National_Spatial_Reference_System_2011",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree"$
  <GeoTransform>  6.1400000000000000e+05,  1.0000000000000000e+00,  0.0000000000000000e+00,  3.8120000000000000e+06,  0.0000000000000000e+00, -1.0000000000000000e+00</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>-9999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">/vsis3/bucketname/TESTING/20190314_rastername_1.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1000" RasterYSize="1000" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
      <DstRect xOff="0" yOff="25000" xSize="1000" ySize="1000" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">/vsis3/bucketname/TESTING/20190314_rastername_2.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1000" RasterYSize="1000" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
      <DstRect xOff="1000" yOff="26000" xSize="1000" ySize="1000" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">/vsis3/bucketname/TESTING/20190314_rastername_3.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1000" RasterYSize="1000" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
      <DstRect xOff="60000" yOff="41000" xSize="1000" ySize="1000" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">/vsis3/bucketname/TESTING/20190314_rastername_4.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1000" RasterYSize="1000" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
      <DstRect xOff="60000" yOff="40000" xSize="1000" ySize="1000" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
  </VRTRasterBand>
  <MaskBand>
    <VRTRasterBand dataType="Byte">
      <SimpleSource>
        <SourceFilename relativeToVRT="0">/vsis3/bucketname/TESTING/20190314_rastername_1.tif</SourceFilename>
        <SourceBand>mask,1</SourceBand>
        <SourceProperties RasterXSize="1000" RasterYSize="1000" DataType="Byte" BlockXSize="256" BlockYSize="256" />
        <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
        <DstRect xOff="0" yOff="25000" xSize="1000" ySize="1000" />
      </SimpleSource>
      <SimpleSource>
        <SourceFilename relativeToVRT="0">/vsis3/bucketname/TESTING/20190314_rastername_2.tif</SourceFilename>
        <SourceBand>mask,1</SourceBand>
        <SourceProperties RasterXSize="1000" RasterYSize="1000" DataType="Byte" BlockXSize="256" BlockYSize="256" />
        <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
        <DstRect xOff="1000" yOff="26000" xSize="1000" ySize="1000" />
      </SimpleSource>
      <SimpleSource>
        <SourceFilename relativeToVRT="0">/vsis3/bucketname/TESTING/20190314_rastername_3.tif</SourceFilename>
        <SourceBand>mask,1</SourceBand>
        <SourceProperties RasterXSize="1000" RasterYSize="1000" DataType="Byte" BlockXSize="256" BlockYSize="256" />
        <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
        <DstRect xOff="60000" yOff="41000" xSize="1000" ySize="1000" />
      </SimpleSource>
      <SimpleSource>
        <SourceFilename relativeToVRT="0">/vsis3/bucketname/TESTING/20190314_rastername_4.tif</SourceFilename>
        <SourceBand>mask,1</SourceBand>
        <SourceProperties RasterXSize="1000" RasterYSize="1000" DataType="Byte" BlockXSize="256" BlockYSize="256" />
        <SrcRect xOff="0" yOff="0" xSize="1000" ySize="1000" />
        <DstRect xOff="60000" yOff="40000" xSize="1000" ySize="1000" />
      </SimpleSource>
    </VRTRasterBand>
  </MaskBand>

mccarthyryanc avatar Mar 18 '19 22:03 mccarthyryanc

@dionhaefner and @j08lue, something strange happened and it actually did work... However, it took probably 30min before the tiles were rendered. I tried to retest today, saw that nothing was rendering and then accidentally left the preview app open on a browser and the server running. After something like 30min tiles started popping up!

So, maybe this is something to do with running terracotta server on an EC2 instance and then terracotta connect on a local machine?

~~My other guess would be slow connection speed, but I tested again while monitoring network bandwidth and saw hardly any download/upload for the local machine/EC2 instance.~~ Actually, after monitoring network activity again it looks like the EC2 instance running terracotta server is maxed out on download trying to pull rasters from S3.

Perhaps terracotta is reverting to non cloud optimized geotiffs reads when given a VRT file, even though the underlying rasters are cloud optimized?

mccarthyryanc avatar Mar 18 '19 22:03 mccarthyryanc

Your log even says that:

[!] /home/ubuntu/src/terracotta/terracotta/drivers/sqlite.py:318: PerformanceWarning: Raster file /tmp/20190315_testing.vrt is not a valid cloud-optimized GeoTIFF. Any interaction with it will be significantly slower. Consider optimizing it through terracotta optimize-rasters before ingestion.

That's certainly bad news; maybe there is something we can do to optimize the VRT file, otherwise we might have to take it up with Rasterio / GDAL people.

But you said the preview loaded fast?

dionhaefner avatar Mar 19 '19 09:03 dionhaefner

@dionhaefner , it loaded faster than the tiles on the map, but still slower than serving a large geotiff.

mccarthyryanc avatar Mar 19 '19 14:03 mccarthyryanc

From the GDAL docs:

in the general case, the VRT bands themselves will not expose overviews.

But there seem to be ways to build overviews for VRTs (.vrt.ovr). Don't know whether you can get to the full COG specs, though.

j08lue avatar Mar 19 '19 18:03 j08lue

Maybe there is some neat way to create a scale-dependent combining internal (GeoTIFF) overviews and VRT-level ones, like mentioned here:

build a few internal overview levels (2 4 8) and then create a new physical file from .vrt with pixel size of the next overview (16). Create internal overviews for the subsampled image and finally combine all together into a scale dependent group.

GDAL black magic...

j08lue avatar Mar 19 '19 18:03 j08lue

I took the discussion here: https://rasterio.groups.io/g/main/topic/30751150#167

If everything works as it should, the performance issues should only be present for low zoom levels, for which we don't have overviews.

dionhaefner avatar Mar 25 '19 15:03 dionhaefner

@dionhaefner, I'm not sure I follow. Based on the groups.io post the workflow is:

  1. Create COGs for individual rasters via terracotta optimize-rasters
  2. Create VRT for rasters via gdalbuildvrt
  3. Create external overviews for VRT via gdaladdo -ro flie.vrt OVERVIEWLEVELS

Does step 3 require the same overiew levels as provided during the COG creation step?

mccarthyryanc avatar Apr 10 '19 22:04 mccarthyryanc

Sorry, haven't pursued this any further. I think you got it right, just two points:

  • Terracotta tells GDAL to disregard external auxiliary files, so you will have to disable that to see any effect of the VRT level overviews. Replacing _RIO_ENV_KEYS here with an empty dict should do the trick:

    https://github.com/DHI-GRAS/terracotta/blob/b7d2fce976895586ccfe6fa545cf2402903247a2/terracotta/drivers/raster_base.py#L47

  • There is no need to put overview levels into the external overview file that are already covered by the internal dataset overviews. So if your VRT is made up of 4 x 4 rasters, it should be sufficient to add 2 additional levels in your last step.

All of this is pretty experimental and I don't have any idea whether it will work. But it would be cool if you could play with it and report back.

dionhaefner avatar Apr 13 '19 23:04 dionhaefner

@dionhaefner sorry for the long wait, but I finally found some time to test this out! TL;DR: It works until I zoom in too far on the map.

I created a new branch and conda env to keep things isolated:

# New conda env
conda create -n tc-vrt python=3.6 -c conda-forge
conda activate tc-vrt
pip install git+https://github.com/mccarthyryanc/terracotta.git@RIO_ENV_KEYS

I had a collection of 873 tiles that I added to a text file vrt_filelist, eg:

/vsis3/bucket/prefix/tile1.tif
/vsis3/bucket/prefix/tile2.tif
...
/vsis3/bucket/prefix/tile873.tif

Built the VRT and overviews:

# Build VRT: This step took 1 min and 41 sec to finish
gdalbuildvrt -input_file_list vrt_filelist ovr_testing.vrt
# Build overviews: This step took 4 mins ans 52 sec to finish. OVR file was 385M
gdaladdo -ro --config COMPRESS_OVERVIEW DEFLATE ovr_testing.vrt 2 4 8 16

Serve and connect like usual:

# Serve from EC2 instance
terracotta serve --allow-all-ips --port 8788 -r ./{type}_{name}.vrt
# Connect from local instance
terracotta connect INSTANCE-IP:8788

At this point I can pan and zoom without issues, but subjectively it does seem a little slower than when using a single large GeoTiff. If I zoom in too far, I see no tiles and the server spits out the following error:

Server Error

[2019-05-06 22:58:20,339] ERROR in app: Exception on /singleband/ovr/testing/17/34397/52244.png [GET]
Traceback (most recent call last):
  File "rasterio/_io.pyx", line 698, in rasterio._io.DatasetReaderBase._read
  File "rasterio/shim_rasterioex.pxi", line 133, in rasterio._shim.io_multi_band
  File "rasterio/_err.pyx", line 182, in rasterio._err.exc_wrap_int
rasterio._err.CPLE_AppDefinedError: IReadBlock failed at X offset 0, Y offset 0: '/vsis3/bucket/prefix/tile330.tif' does not exist in the file system, and is not recognized as a supported dataset name.

During handling of the above exception, another exception occurred:

Traceback (most recent call last): File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/flask/app.py", line 2292, in wsgi_app response = self.full_dispatch_request() File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/flask/app.py", line 1815, in full_dispatch_request rv = self.handle_user_exception(e) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/flask/app.py", line 1718, in handle_user_exception reraise(exc_type, exc_value, tb) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/flask/_compat.py", line 35, in reraise raise value File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/flask/app.py", line 1813, in full_dispatch_request rv = self.dispatch_request() File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/flask/app.py", line 1799, in dispatch_request return self.view_functionsrule.endpoint File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/terracotta/server/flask_api.py", line 50, in inner return fun(*args, **kwargs) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/terracotta/server/singleband.py", line 121, in get_singleband return _get_singleband_image(keys, tile_xyz) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/terracotta/server/singleband.py", line 166, in _get_singleband_image image = singleband(parsed_keys, tile_xyz=tile_xyz, **options) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/contextlib.py", line 52, in inner return func(*args, **kwds) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/terracotta/handlers/singleband.py", line 45, in singleband tile_size=tile_size, preserve_values=preserve_values File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/terracotta/xyz.py", line 46, in get_tile_data preserve_values=preserve_values, asynchronous=asynchronous File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/terracotta/drivers/base.py", line 20, in inner return fun(self, *args, **kwargs) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/terracotta/drivers/raster_base.py", line 566, in get_raster_tile return task() File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/cachetools/init.py", line 87, in wrapper v = method(self, *args, **kwargs) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/contextlib.py", line 52, in inner return func(*args, **kwds) File "/home/ubuntu/miniconda/envs/tc-vrt/lib/python3.6/site-packages/terracotta/drivers/raster_base.py", line 525, in _get_raster_tile 1, resampling=resampling_enum, window=out_window, out_shape=tile_size File "rasterio/_warp.pyx", line 937, in rasterio._warp.WarpedVRTReaderBase.read File "rasterio/_io.pyx", line 349, in rasterio._io.DatasetReaderBase.read File "rasterio/_io.pyx", line 701, in rasterio._io.DatasetReaderBase._read rasterio.errors.RasterioIOError: Read or write failed. IReadBlock failed at X offset 0, Y offset 0: '/vsis3/bucket/prefix/tile330.tif' does not exist in the file system, and is not recognized as a supported dataset name.

I imagine the error is the result of my choice in overview levels. The tiles did not form a nice grid. I did try this with a subset of tiles that made a 4x4 and built overview levels at 2 and 4 (like you suggested). This worked without any issues.

What I didn't try was uploading the VRT/OVR files to S3 and inserting an entry to MySQL via terracotta driver to make it work from Lambda. But at this point I feel like the data has jumped through a few too many hoops for me to use this in production... I think would rather wait (willing to help where I can) until terracotta supported seamless layers.

mccarthyryanc avatar May 06 '19 23:05 mccarthyryanc

Thanks for reporting back!

rasterio.errors.RasterioIOError: Read or write failed. IReadBlock failed at X offset 0, Y offset 0: '/vsis3/bucket/prefix/tile330.tif' 

Does this file exist / contain a valid GTiff?

dionhaefner avatar May 07 '19 06:05 dionhaefner

Yep, file exists and is a valid GeoTiff.

mccarthyryanc avatar May 07 '19 15:05 mccarthyryanc

Very odd, the error message reads almost like https://github.com/DHI-GRAS/terracotta/issues/139. If we could only get hold of a small example that demonstrates this bug...

dionhaefner avatar May 08 '19 11:05 dionhaefner

I've tried recreating this issue with a small set of rasters that I can share, but haven't had any luck.

@dionhaefner in your original post you mention another workaround: adding all rasters as separate layers and displaying them on top of eachother. Do you have any advice on trying that? I'm assuming this would be part of any viewing app. For example, adding a ton of XYZ TileLayers to a leaftlet map?

mccarthyryanc avatar May 10 '19 16:05 mccarthyryanc

I've tried recreating this issue with a small set of rasters that I can share, but haven't had any luck.

Thanks for trying! I think I might ask someone from mapbox to have a look at this, maybe they have an idea what's going wrong.

For example, adding a ton of XYZ TileLayers to a leaftlet map?

Yes, exactly, you would do that in the client. But it does become infeasible at some point (I wouldn't do it for hundreds of layers), because then Terracotta gets bombarded for requests for tiles that are out of range. I think the only fully supported solution is to create one or a handful of big GeoTIFFs 😕

dionhaefner avatar May 11 '19 10:05 dionhaefner

👋 Because I've never been able to work with VRT (and that you mostly need to create external overviews) I've been working on a simple (but fast) solution to get mosaic tiles from multiple COGs, the result is explained in https://medium.com/devseed/cog-talk-part-2-mosaics-bbbf474e66df

The general idea is to create a json files linking the COGs together (the gzip .json file should be really small, but it could also be a database item) and then use https://github.com/cogeotiff/rio-tiler-mosaic to merge the tile arrays together (with a pixel selection method)

We have an implementation example over https://github.com/developmentseed/cogeo-mosaic

mosaicJSON is still a WIP (work in progress) and we are looking for feedback 😄 .

Please feel free to ping me if you have any questions.

vincentsarago avatar May 27 '19 14:05 vincentsarago

Do I understand that correctly that you don't have overviews for the entire mosaic? So if I zoom out rerally far it touches all of the COGs?

dionhaefner avatar May 28 '19 10:05 dionhaefner