Out of memory when optimizing large rasters even with --no-in-memory option

[htejwani]

Command used:

terracotta optimize-rasters --resampling-method nearest Raster300WMTiled.tif --no-in-memory --output-folder 
optimized_cogs/

Here is the error output after 1.5 hours of the command execution. Size of the output file before the error occurs is ~14GB:

[-] Uncaught exception!                                                                                                                                                  
Traceback (most recent call last):
  File "/home/SVaid/workspace/terracotta/terracotta/terracotta/scripts/cli.py", line 52, in entrypoint
    cli(obj={})
  File "/home/SVaid/miniconda3/envs/terracotta/lib/python3.9/site-packages/click/core.py", line 829, in __call__
    return self.main(*args, **kwargs)
  File "/home/SVaid/miniconda3/envs/terracotta/lib/python3.9/site-packages/click/core.py", line 782, in main
    rv = self.invoke(ctx)
  File "/home/SVaid/miniconda3/envs/terracotta/lib/python3.9/site-packages/click/core.py", line 1259, in invoke
    return _process_result(sub_ctx.command.invoke(sub_ctx))
  File "/home/SVaid/miniconda3/envs/terracotta/lib/python3.9/site-packages/click/core.py", line 1066, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/home/SVaid/miniconda3/envs/terracotta/lib/python3.9/site-packages/click/core.py", line 610, in invoke
    return callback(*args, **kwargs)
  File "/home/SVaid/workspace/terracotta/terracotta/terracotta/scripts/optimize_rasters.py", line 269, in optimize_rasters
    dst.build_overviews(overviews, rs_method)
  File "rasterio/_io.pyx", line 1619, in rasterio._io.DatasetWriterBase.build_overviews
  File "rasterio/_err.pyx", line 190, in rasterio._err.exc_wrap_int
rasterio._err.CPLE_OutOfMemoryError: overview.cpp, 87: cannot allocate 6183031428 bytes

Here is the gdalinfo output of the tif file:

Driver: GTiff/GeoTIFF
Files: Raster300WM.tif
Size is 1516563, 32061
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-10580430.809917071834207,4852233.339643679559231)
Pixel Size = (1.302869663735731,-1.302869663735710)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=ZSTD
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-10580430.810, 4852233.340) ( 95d 2'44.26"W, 39d54'20.15"N)
Lower Left  (-10580430.810, 4810462.035) ( 95d 2'44.26"W, 39d37' 1.73"N)
Upper Right (-8604546.884, 4852233.340) ( 77d17'45.46"W, 39d54'20.15"N)
Lower Right (-8604546.884, 4810462.035) ( 77d17'45.46"W, 39d37' 1.73"N)
Center      (-9592488.847, 4831347.687) ( 86d10'14.86"W, 39d45'41.48"N)
Band 1 Block=1516563x1 Type=Float32, ColorInterp=Gray
  NoData Value=0

[dionhaefner]

It looks like you're running out of memory when building overviews. This is handled by rasterio and nothing we have control over.

Did you try to use GDAL directly?

[htejwani]

Thanks for your reply. Tried it with GDAL, the COG tif size is 37GB. When served with terracotta, the rendering on map fails intermittently or freezes.

[dionhaefner]

Can you be more specific? Which command did you use? What does "fails" mean? Do you see any errors or warnings in the logs?

[htejwani]

1. GDAL command used: gdal_translate input.tif output_cog.tif -of COG -co BIGTIFF=YES --config GDAL_CACHEMAX 6144 -co COMPRESS=ZSTD -co NUM_THREADS=8
2. The tif file size is 37GB
3. On leaflet, the layer initially takes more than 15 seconds to appear. When I zoom-in to certain parts, at times the tiles are rendered quickly but in most cases, the map freezes and does not render the tiles or takes minutes to render.
4. I could not find any error info in console logs except the traffic logs.

[dionhaefner]

The file size shouldn't be a problem per se, we've served files in that ballpark without issues. Right, @j08lue ?

I'm not sure how smart the new GDAL is. Can you run it through a COG validator to make sure that it's a valid COG? And can you post the output of gdalinfo again?

[htejwani]

Yes, I checked wit rio-cogeo library and gdal's vaidation script. Both confirm that it is a valid COG.
Here's the output from gdalinfo command:

Driver: GTiff/GeoTIFF
Files: Raster300WMGDALCOG.tif
Size is 1516563, 32061
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-10580430.809917071834207,4852233.339643679559231)
Pixel Size = (1.302869663735731,-1.302869663735710)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=ZSTD
  INTERLEAVE=BAND
  LAYOUT=COG
Corner Coordinates:
Upper Left  (-10580430.810, 4852233.340) ( 95d 2'44.26"W, 39d54'20.15"N)
Lower Left  (-10580430.810, 4810462.035) ( 95d 2'44.26"W, 39d37' 1.73"N)
Upper Right (-8604546.884, 4852233.340) ( 77d17'45.46"W, 39d54'20.15"N)
Lower Right (-8604546.884, 4810462.035) ( 77d17'45.46"W, 39d37' 1.73"N)
Center      (-9592488.847, 4831347.687) ( 86d10'14.86"W, 39d45'41.48"N)
Band 1 Block=512x512 Type=Float32, ColorInterp=Gray
  NoData Value=0
  Overviews: 758281x16030, 379140x8015, 189570x4007, 94785x2003, 47392x1001, 23696x500, 11848x250, 5924x125, 2962x62, 1481x31, 740x15, 370x7

[dionhaefner]

That looks mostly OK to me, the only smaller issue is Block=512x512. Terracotta uses 256x256 by default, so you lose some efficiency there. But the overviews seem alright.

You do ingest the raster before serving it, right? And how do you serve?

Edit: Another thing is INTERLEAVE=BAND, where we use INTERLEAVE=PIXEL (I don't know if this makes a difference in single-banded images)

[j08lue]

2\. The tif file size is 37GB

we've served files in that ballpark without issues. Right, @j08lue?

In daily use, our images are mostly in the order of a 10th of that size on disk (1-5 GB). But size on disk is one thing with compression and all. If your raster is very sparsely filled or has a lot of areas of same value, it will compress very well. From your gdalinfo, 1516563 x 32061 pixels in 32 bit float is actually a 200 GB array (without compression). Do we have any sense of how large those GeoTIFF headers with the overviews are getting?