EPSG and its effect on COG and Tiles Opacity

geoffreynyaga commented 2 years ago

Great work with the library.

I have an issue with opacity/clarity of tiles, thumbnails and overall Geotiffs/COGs as well as COGs not getting tiled.

Opacity Issue: When using the swagger API to get a thumbnail of an existing COG, if I include "EPSG:3857" on the request, I get a washed out image

However, if I omit the EPSG Code, I get a clear image thumbnail. I have tested with multiple COGs and raw GeoTiffs

In Django Admin, the same Geotiffs and COGs are also washed out

I have also used a separate Leaflet html page to access the /api/maps/{id}/{z}/{x}/{y}.png and I get the same washed out images

<html>
  <head>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no" />
    <meta charset="utf-8" />
    <style>
      html,
      body {
        height: 100%;
        margin: 0;
        padding: 0;
      }
      #map {
        height: 100%;
      }
    </style>
    <link
      rel="stylesheet"
      href="https://unpkg.com/leaflet@1.8.0/dist/leaflet.css"
      integrity="sha512-hoalWLoI8r4UszCkZ5kL8vayOGVae1oxXe/2A4AO6J9+580uKHDO3JdHb7NzwwzK5xr/Fs0W40kiNHxM9vyTtQ=="
      crossorigin=""
    />
    <script
      src="https://unpkg.com/leaflet@1.8.0/dist/leaflet.js"
      integrity="sha512-BB3hKbKWOc9Ez/TAwyWxNXeoV9c1v6FIeYiBieIWkpLjauysF18NzgR1MBNBXf8/KABdlkX68nAhlwcDFLGPCQ=="
      crossorigin=""
    ></script>
  </head>
  <body>
    <div id="map"></div>
    <script>
      var map = L.map("map").setView(
        [-0.24511062205544606, 34.87011481076246],
        16
      );
      L.tileLayer("https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png", {
        attribution:
          '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors',
        maxZoom: 18,
      }).addTo(map);
      L.tileLayer(
        "http://localhost:8000/api/maps/1/tiles/{z}/{x}/{y}.png?projection=EPSG%3A3857",
        {
          attribution:
            '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors',
          maxZoom: 23,
          opacity: 1,
        }
      ).addTo(map);

    </script>
  </body>
</html>

However, in the above example, if I do not append the EPSG code on the url in html, then I get an error (no map rendered) and a 400 bad request

As for individual tiles, if I include the EPSG Code in the request body in swagger UI, i get a washed out tile

However, if I omit the EPSG code, I get a bad request error

An example COG metadata is

{
  "geospatial": true,
  "levels": 8,
  "sizeX": 32465,
  "sizeY": 28586,
  "sourceLevels": 8,
  "sourceSizeX": 32465,
  "sourceSizeY": 28586,
  "tileWidth": 256,
  "tileHeight": 256,
  "bounds": {
    "ll": {
      "x": 34.85392278135806,
      "y": -0.23501171685313998
    },
    "ul": {
      "x": 34.85392278135806,
      "y": -0.2273338589976692
    },
    "lr": {
      "x": 34.86258419062559,
      "y": -0.23501171685313998
    },
    "ur": {
      "x": 34.86258419062559,
      "y": -0.2273338589976692
    },
    "srs": "EPSG:4326",
    "xmin": 34.85392278135806,
    "xmax": 34.86258419062559,
    "ymin": -0.23501171685313998,
    "ymax": -0.2273338589976692
  },
  "sourceBounds": {
    "ll": {
      "x": 34.85392278135806,
      "y": -0.23501171685313998
    },
    "ul": {
      "x": 34.85392278135806,
      "y": -0.2273338589976692
    },
    "lr": {
      "x": 34.86258419062559,
      "y": -0.23501171685313998
    },
    "ur": {
      "x": 34.86258419062559,
      "y": -0.2273338589976692
    },
    "srs": "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs",
    "xmin": 34.85392278135806,
    "xmax": 34.86258419062559,
    "ymin": -0.23501171685313998,
    "ymax": -0.2273338589976692
  },
  "bands": {
    "1": {
      "min": 55,
      "max": 255,
      "mean": 193.05740298087738,
      "stdev": 72.13546456636203,
      "units": "metre",
      "interpretation": "red",
      "maskband": 4
    },
    "2": {
      "min": 63,
      "max": 255,
      "mean": 194.14802446569178,
      "stdev": 70.6158617103163,
      "units": "metre",
      "interpretation": "green",
      "maskband": 4
    },
    "3": {
      "min": 41,
      "max": 255,
      "mean": 182.4484673790776,
      "stdev": 84.39825325395805,
      "units": "metre",
      "interpretation": "blue",
      "maskband": 4
    },
    "4": {
      "min": 0,
      "max": 255,
      "mean": 109.17366774465692,
      "stdev": 126.17604981162934,
      "units": "metre",
      "interpretation": "alpha"
    }
  },
  "magnification": null,
  "mm_x": 2.969893339423529,
  "mm_y": 2.969893339423529,
  "proj4": [
    "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"
  ],
  "frames": false
}

banesullivan commented 2 years ago

Thank you for filing such a detailed issue! This is indeed very odd. I'm curious if you would be able to to share the raster itself for me to reproduce or at least try sharing a "cleaned" version of the file by using: rio-faux https://github.com/cogeotiff/rio-faux

banesullivan commented 2 years ago

Please note the interpretation in band 4.

"4": {
      "min": 0,
      "max": 255,
      "mean": 109.17366774465692,
      "stdev": 126.17604981162934,
      "units": "metre",
      "interpretation": "alpha"
    }

from the metadata about the bands. It appears as though band 4 is always interpreted as alpha/opacity when it may not actually be opacity.

The best way to handle this would be to update the source raster file to change that band's interpretation, or perhaps there is a bug in large-image itself where band 4 is forcibly interpreted as alpha (hard to tell without the source file). cc @manthey

A workaround for this is to specify the bands when viewing these particular images via the style parameter on the endpoints:

// See https://girder.github.io/large_image/tilesource_options.html#style
var style = {
  bands: [
    {band: 1, palette: ['#000', '#f00']},  // red
    {band: 2, palette: ['#000', '#0f0']},  // green
    {band: 3, palette: ['#000', '#00f']}   // blue
  ]
};
var styleEncoded = encodeURIComponent(JSON.stringify(style))
var thumbnailUrl = `http://localhost:8000/api/image-file/${imageId}/data/thumbnail.png?style=${styleEncoded}`;

Here is the stringified value for that which you could copy/paste as the style parameter of these endpoints:

{"bands":[{"band":1,"palette":["#000","#f00"]},{"band":2,"palette":["#000","#0f0"]},{"band":3,"palette":["#000","#00f"]}]}

banesullivan commented 2 years ago

@manthey, perhaps we could add an option in large-image to ignore the alpha channel when getting thumbnails/tiles?

geoffreynyaga commented 2 years ago

@banesullivan @manthey I have uploaded sample dataset on the link below

https://drive.google.com/drive/folders/1JS6hFnjaGTYJe8_R7-tEzojf0_qHMMgX?usp=sharing

I believe its an alpha channel issue

I have tried converting the COGs with multiple tools and command line GDAL but I get same result.

However, if I "don't convert background colour to transparent" and leave it to white or black, the COGs are visible clearly at all zoom levels; we just now have an unwanted white bounding box around the raster. I also get the same exact resolution in swagger playground, regardless whether I include the EPSG Code or not.

But if I make the background void tiles transparent, the washed out issue arises.

In the folder, I have uploaded both instances of a raster COG with transparent background, one with a white background and a third COG for comparison.

banesullivan commented 2 years ago

Alpha channel

This is quite odd. The alpha channel seems to be a mask of data vs no data (which is good):

import large_image
import rasterio as rio
import matplotlib.pyplot as plt
import numpy as np
import PIL.Image
import io

ds = rio.open(path)
band = ds.read(4)
plt.hist(band.ravel())

download

>>> np.unique(band.ravel())
array([  0, 255], dtype=uint8)

And we can see this clearly when getting a JPEG thumbnail:

src = large_image.open(path, projection='EPSG:3857', style={'bands': [{'band': 4, 'palette': ['#000', '#fff']}]})
src.getThumbnail(encoding='JPEG')[0]

download

However, the PNG still looks off (look at with dark/black background to see features):

src.getThumbnail(encoding='PNG')[0]

download

The PNG encoded thumbnail is very wrong. You can see structure from other bands and thus it's not actually giving a thumbnail of just that band but some sort of composite (ignoring the style) -- this might be a bug in large_image

I have no idea where the mid-range opacity/alpha values are coming from... thoughts, @manthey ?

thumb = src.getThumbnail(encoding='PNG')[0]
im = PIL.Image.open(io.BytesIO(thumb))
pix = np.array(im)
plt.hist(pix[:,:,-1].ravel())

download

@geoffreynyaga, we will have to continue looking into this.

FYI, these might not be COGs

@geoffreynyaga, what did you use to convert your images to COGs? According to the GDAL check, these are not totally valid COGs:

>>> import large_image
>>> src = large_image.open('4-ortho-rgb-cog.tif')
>>> src.validateCOG()
---------------------------------------------------------------------------
TileSourceInefficientError                Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 src.validateCOG()

File ~/Software/kw/large_image/sources/gdal/large_image_source_gdal/__init__.py:1227, in GDALFileTileSource.validateCOG(self, check_tiled, full_check, strict, warn)
   1221 warnings, errors, details = validate(
   1222     self._largeImagePath,
   1223     check_tiled=check_tiled,
   1224     full_check=full_check
   1225 )
   1226 if errors:
-> 1227     raise TileSourceInefficientError(errors)
   1228 if strict and warnings:
   1229     raise TileSourceInefficientError(warnings)

TileSourceInefficientError: ['The offset of the main IFD should be 8. It is 455085854 instead', 'The offset of the first block of overview of index 3 should be after the one of the overview of index 4', 'The offset of the first block of overview of index 2 should be after the one of the overview of index 3', 'The offset of the first block of overview of index 1 should be after the one of the overview of index 2', 'The offset of the first block of overview of index 0 should be after the one of the overview of index 1', 'The offset of the first block of the main resolution image should be after the one of the overview of index 4']

FWIW, I have written up a few things on converting COGs:

manthey commented 2 years ago

This turns out to be a typo in the large-image-gdal-source module where there is code to apply fallback alpha information if there is no alpha channel. Due to a typo, this was being invoked even though there was an alpha channel. See https://github.com/girder/large_image/pull/909.

banesullivan commented 2 years ago

Thanks for digging into this, @manthey!

@geoffreynyaga, this will be fixed on the next release of large-image

girder / django-large-image

EPSG and its effect on COG and Tiles Opacity #39

Alpha channel

FYI, these might not be COGs