Raster image stretch e.g. min max stretching

[jordancaraballo]

Hi, has there been any work or examples on stretching the rasters when used with localtileserver.get_leaflet_tile_layer? Working with uint16 raster at the time and imagery is extremely dark. I have seen this example using the ee library, https://github.com/giswqs/earthengine-py-notebooks/blob/master/Visualization/image_stretch.py, but I have not seen any available args in order to work on stretching the image (https://girder.github.io/large_image/tilesource_options.html#style). Thanks!

[banesullivan]

This should be doable, and perhaps there is a terminology difference between EE and large-image. With large-image, I think the min/max options are what you need.

Would you please provide some sample data and perhaps an example (screenshot) of the desired output so that I can try to put an example together for you?

[giswqs]

TiTiler uses a parameter called rescale, so you can set something like rescale=0,1000 (same for all bands) or rescale=["164,223","130,211","99,212"] (different values for different bands). See https://github.com/giswqs/leafmap/pull/299

[giswqs]

Here are two examples for showing the differences. The leafmap.add_cog_layer uses TiTiler while leafmap.add_raster() uses localtileserver. It would be great if localtileserver can support this as well.

m = leafmap.Map()
url = 'https://opendata.digitalglobe.com/events/california-fire-2020/post-event/2020-08-14/pine-gulch-fire20/10300100AAC8DD00.tif'
m.add_cog_layer(url)
m

m = leafmap.Map()
url = 'https://opendata.digitalglobe.com/events/california-fire-2020/post-event/2020-08-14/pine-gulch-fire20/10300100AAC8DD00.tif'
m.add_cog_layer(url, rescale='0, 255')
m

[banesullivan]

Hm. I'm not actually needing to rescale that image when using localtileserver/large-image as the GeoTIFF properly labels the RGB channels. However, I can show how to rescale with this image (or any image) using the style specification from large-image (reference https://localtileserver.banesullivan.com/user-guide/rgb.html and https://girder.github.io/large_image/tilesource_options.html#style):

from localtileserver import TileClient, get_leaflet_tile_layer, examples
from ipyleaflet import Map

url = 'https://opendata.digitalglobe.com/events/california-fire-2020/post-event/2020-08-14/pine-gulch-fire20/10300100AAC8DD00.tif'

client = TileClient(url)

# Using https://girder.github.io/large_image/tilesource_options.html#style
style = {
  'bands': [
    {'band': 1, 'palette': '#f00', 'min': 75, 'max': 150},
    {'band': 2, 'palette': '#0f0', 'min': 25, 'max': 200},
    {'band': 3, 'palette': '#00f', 'min': 150, 'max': 175},
  ]
}

t = get_leaflet_tile_layer(client, style=style)
m = Map(center=client.center(), zoom=client.default_zoom)
m.add_layer(t)
m

@jordancaraballo, does this accomplish what you are looking for?

Please also note there is a vmin/vmax API within localtileserver as well but I've found working directly with the large-image styling specification to be the best solution.

[jordancaraballo]

@banesullivan thank you very much for the example, that is definitely going in the right direction and it changes the colors extremely well. This can definitely go in one of the example notebooks, I will try to create one with what I have if you think it is appropriate to include. The downside for this I guess is the manual lookup of the min and max values which I just gathered with rasterio on the fly.

Following the same line of thought, the example you provided is this min/max stretch, do you know of similar ways to mimic then mean/standard deviation stretching? large_image does not seem to provide such option, besides the precomputed gamma example from their site.

[jordancaraballo]

I think I answered my own question. Simple example below. I think this can be made as an option of the localtileserver get_leaflet_tile_layer function. Something like image_stretch with min/max, mean/std, etc. options. I can look at the source and add something cleaner if you think that would be a good addition. Ideally we can calculate the metadata on the fly and pass the attribute to the large_image style parameter.

PS: Cannot share the data example since this is Maxar WorldView data under a signed agreement, but I will find another example that can actually be shared.

import rioxarray as rxr
import numpy as np

filename = '/home/jovyan/efs/projects/3sl/data/Tappan/Tappan02_WV02_20120218_M1BS_103001001077BE00_data.tif'

sigma = 2
x = rxr.open_rasterio(filename)
x = x.where(x>-10001,np.nan)

print(x.min().values, x.max().values)

b7mean = x[6, :, :].mean().values
b3mean = x[2, :, :].mean().values
b2mean = x[1, :, :].mean().values

b7std = x[6, :, :].std().values
b3std = x[2, :, :].std().values
b2std = x[1, :, :].std().values

b7newmin = b7mean - (b7std * sigma)
b7newmax = b7mean + (b7std * sigma)

b3newmin = b3mean - (b3std * sigma)
b3newmax = b3mean + (b3std * sigma)

b2newmin = b2mean - (b2std * sigma)
b2newmax = b2mean + (b2std * sigma)

print(b7newmin, b7newmax) # used in style
print(b3newmin, b3newmax) # used in style
print(b2newmin, b2newmax) # used in style

[banesullivan]

Fantastic example!

I'll see if I can add some helper methods to build the styling parameters for different types of stretches like this. If I implement, @jordancaraballo, would you be willing to review the work?

[jordancaraballo]

Yes, I would be willing to take a look at it. Also, just noticed the .rasterio option from the TileClient has all the needed statistics, so it will be even easier to implement.

raster_client = TileClient(in_raster)
print("RASTER_CLIENT", raster_client.rasterio.statistics(7))
RASTER_CLIENT Statistics(min=718.0, max=8111.0, mean=2116.8099705569166, std=365.11115105581723)

[banesullivan]

This is reassuring me that #111 and https://github.com/girder/large_image/pull/927 would be hugely beneficial... Would be great if I could refactor localtileserver to be rasterio centric and only do the tile serving through large-image... maybe in time!

---

I'll get working on these stretching options when I have some time!

[jordancaraballo]

Indeed, #111 would make a lot of more interactions with rasterio easier. This is my quick bandaid for the dashboard I am building, hopefully I will remove it after your implementation. Thanks again for looking into this!

        # create TileClient object
        raster_client = TileClient(in_raster)

        style_list = []
        for bid, pid in zip(data_bands, ['#f00', '#0f0', '#00f']):
            band_stats = raster_client.rasterio.statistics(bid)
            newmin = band_stats.mean - (band_stats.std * sigma)
            newmax = band_stats.mean + (band_stats.std * sigma)
            style_list.append(
                {'band': bid, 'palette': pid, 'min': newmin, 'max': newmax})

        self.add_layer(
            get_leaflet_tile_layer(
                raster_client, show=False, band=data_bands,
                cmap=cmap, max_zoom=self.default_zoom,
                max_native_zoom=self.default_zoom,
                name=layer_name, scheme='linear',
                dtype='uint16', style={'bands': style_list}
            )
        )
        self.center = raster_client.center()  # center Map around raster
        self.zoom = raster_client.default_zoom  # zoom to raster center