Add stretch options for visualizing images

Question

Add stretch options for visualizing images

aazuspan opened this issue a year ago · 8 comments

Description

Hi folks! I'd like to suggest adding "Stretch" options to adjust the vis param range for images automatically, similar to what can be done in the Code Editor. I'm happy to implement this and make a PR.

Implementation

I would suggest we follow the lead of the Code Editor add a dropdown menu to select different percent and sigma stretches. When selected, reducers would be run on the selected image bands over the current map bounds, and values would be pulled client-side and processed to set the current range. I believe all the different percentile and sigma stretches can be calculated using min, max, and std dev reducers. As in the Code Editor, we could disable the range inputs when a stretch option is selected, adding a Custom option for user-defined values.

API

It looks like there's a lot of great refactoring work going on in geemap, so I don't want to interfere with ongoing plans! My tentative suggestion would be to add a dropdown widget to the range_hbox in the main Map class, with a listener that would call a new method Map.get_stretch_values() that returns the min and max range values for a given ee.Image, but I'm open to any feedback!

UI

It will be a little tricky to fit the range slider, range values, and stretch dropdown on a single row (see screenshot of a quick prototype I put together below).

We could try swapping the slider for two text boxes, more similar to the Code Editor, although that would require a more substantial rewrite. I'm open to any other ideas as well!

Thanks!

Answer 1 · 2023-07-07T12:53:17.000Z

@aazuspan Good suggestion. This would be a nice feature to add. I would suggest adding the stretch dropdown as a new row above the range slider. Similar to the Code Editor, the default value would be Custom, and the options are stretch std and percentages. When the stretch option is changed, the slider beneath is updated automatically.

import ipywidgets as widgets
options = ['Custom', '1σ', '2σ', '3σ', '90%', '98%', '100%']
stretch = widgets.Dropdown(description='Stretch:', value='Custom', options=options)
stretch

This will only affect the Map.create_vis_widget() method. It won't affect the refactoring. You can go ahead submit a PR. Thanks.
https://github.com/gee-community/geemap/blob/master/geemap/geemap.py#L1193

Answer 2 · 2023-07-07T16:35:26.000Z

@aazuspan - this is great!

I was thinking about this a while back. I didn't get anywhere with the UI components but worked through some ideas about the reductions in this draft function. No expectation to use any of it, just wanting to share the product of some thoughts I had. There is probably a wrapper for ee.Image.reduceRegions needed or **kwargs to avoid setting the shared arguments in each reduction (I was just copy/pasting during quick dev). It's intended to get arguments from the Map and the vis selector UI. I was also wanting a way to run it manually without the UI. It would be nice if a person can fetch the results and store them for later use - like once calculated, could be applied to other layers.

I was imagining that this function would be a method of e.g. EELeafletTileLayer. When called, it would fetch the tile layer's associated self.ee_object (if ee.Image) (not a self attribute currently), fetch args from Map and style UI, perform the calculation and then update self.vis_params and self.url_format (w/ _get_tile_url_format(ee_object, self.vis_params)), and then trigger tile redraw. A person could save the vis params by accessing the updated self.vis_params attribute for later use on additional layers.

def image_vis_stats(image: ee.Image, params: Dict) -> Dict:
    """
    Calculate visualization statistics for an Earth Engine Image.

    Args:
        image: An Earth Engine Image.
        params: A dictionary containing the following fields:
            - bands: List of band names to use.
            - bounds: A geometry representing the area to compute statistics.
            - latitude: Latitude of the center of the map.
            - zoom: Zoom level of the map.
            - statistic: The statistic to compute. Must be one of:
                'stddev2', 'percentile', or 'minmax'

    Returns:
        A dictionary containing the following fields:
            - bands: List of band names used.
            - min: Minimum value of the computed statistic.
            - max: Maximum value of the computed statistic.
    """
    bands = params['bands']
    reducer = None
    select_bands = ee.Image(image.select(bands))
    geometry = params['bounds']
    scale = get_resolution(params['latitude'], params['zoom']) # should be able to use Map.get_scale()
    crs = 'EPSG:3857'
    best_effort = True

    if params['statistic'] == "stddev2":
        reducer = ee.Reducer.mean().combine(ee.Reducer.stdDev(), sharedInputs=True)
        stats = select_bands.reduceRegion(
            reducer=reducer,
            geometry=geometry,
            scale=scale,
            crs=crs,
            bestEffort=best_effort
        ).getInfo()
        stddev_from_mean = [stats[f'{band}_mean'] + stats[f'{band}_stdDev'] * 2 for band in bands]
        stddev_from_mean_sum = sum(stddev_from_mean)
        max = stddev_from_mean_sum / len(bands)
        min = -max

    elif params['statistic'] == "percentile":
        reducer = ee.Reducer.percentile([1, 99])
        stats = select_bands.reduceRegion(
            reducer=reducer,
            geometry=geometry,
            scale=scale,
            crs=crs,
            bestEffort=best_effort
        ).getInfo()
        max = sum([stats[f'{band}_p99'] for band in bands]) / len(bands)
        min = sum([stats[f'{band}_p1'] for band in bands]) / len(bands)

    elif params['statistic'] == "minmax":
        reducer = ee.Reducer.minMax()
        stats = select_bands.reduceRegion(
            reducer=reducer,
            geometry=geometry,
            scale=scale,
            crs=crs,
            bestEffort=best_effort
        ).getInfo()
        max = sum([stats[f'{band}_max'] for band in bands]) / len(bands)
        min = sum([stats[f'{band}_min'] for band in bands]) / len(bands)

    else:
        raise ValueError(f"Invalid statistic: {params['statistic']}. Must be one of 'stddev2', 'percentile', or 'minmax'")

    return {'bands': bands, 'min': min, 'max': max}

Answer 3 · 2023-07-07T17:38:15.000Z

I forgot to include the get_resolution() helper, but actually, you should be able to use Map.get_scale(). Except, I'm seeing that it does not account for latitude, it assumes 0.

import math

def get_resolution(latitude: float, zoom_level: int) -> float:
    """
    Calculate the scale (in meters) of a tile layer at a given latitude
    and zoom level. Equation is from: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Resolution_and_Scale

    Args:
        latitude: The latitude of the center of the map viewport.
        zoom_level: The zoom level of the map.

    Returns:
        The resolution of the tile layer in meters.

    Raises:
        ValueError: If `latitude` is not between -90 and 90, or if `zoom_level`
        is less than 0.
    """

    if not (-90 <= latitude <= 90):
        raise ValueError("latitude must be between -90 and 90")

    if zoom_level < 0:
        raise ValueError("zoom_level must be greater than or equal to 0")

    cos_latitude = math.cos(math.radians(latitude))
    return 156543.03 * cos_latitude / (2 ** zoom_level)

Answer 4 · 2023-07-07T19:29:13.000Z

The stats can probably be saved to the Map.ee_layer_dict attribute.

import ee
import geemap
Map = geemap.Map(center=(40, -100), zoom=4)
dem = ee.Image('USGS/SRTMGL1_003')
landsat7 = ee.Image('LANDSAT/LE7_TOA_5YEAR/1999_2003').select(
    ['B1', 'B2', 'B3', 'B4', 'B5', 'B7']
)
states = ee.FeatureCollection("TIGER/2018/States")
vis_params = {
    'min': 0,
    'max': 4000,
    'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5'],
}
Map.addLayer(dem, vis_params, 'SRTM DEM')
Map.addLayer(
    landsat7,
    {'bands': ['B4', 'B3', 'B2'], 'min': 20, 'max': 200, 'gamma': 2.0},
    'Landsat 7',
)
Map.addLayer(states, {}, "US States")
Map

Map.ee_layer_dict

Answer 5 · 2023-07-08T01:40:54.000Z

Thanks for the great feedback, guys!

@giswqs I like the UI suggestions. Adding a separate row for stretching options sounds like a good plan. Thanks also for the tip on Map.ee_layer_dict.

@jdbcode Very helpful notes - you've clearly put a lot of thought into this already! I'll have to familiarize myself a little more with EELeafletTileLayer, but decoupling the vis param calulation from the UI by putting it in there makes perfect sense, and should improve maintainability and testability over my suggestion.

I forgot to include the get_resolution() helper, but actually, you should be able to use Map.get_scale(). Except, I'm seeing that it does not account for latitude, it assumes 0.

The dynamic selection of scale based on zoom is interesting! Poking into the Code Editor implementation, I notice that they just run the reducer at a fixed scale=1 and rely on bestEffort=True and maxPixels=10000 to keep the computation fast. Any thoughts on how this approach would compare? It adds a bit of complexity, but maybe would tend to achieve more accurate stats or better performance?

statistic: The statistic to compute. Must be one of:
'stddev2', 'percentile', or 'minmax'

From your implementation, it looks like you planned to use these instead of the stretch and sigma options from the Code Editor? Want me to follow that lead?

When called, it would fetch the tile layer's associated self.ee_object (if ee.Image) (not a self attribute currently), fetch args from Map and style UI, perform the calculation and then update self.vis_params and self.url_format (w/ _get_tile_url_format(ee_object, self.vis_params)), and then trigger tile redraw.

A question for both of you, probably: I'm a bit confused by the potential duplication between the properties stored in Map.ee_layer_dict and in the EELeafletTileLayer itself. Both seem to store vis_params and ee_object. Would every update to those need to be applied in both places? Just want to make sure I understand the connection between those!

perform the calculation and then update self.vis_params and self.url_format (w/ _get_tile_url_format(ee_object, self.vis_params)

Would it make sense for EELeafletTileLayer.url_format to be converted to a cached_property? By caching it based on the vis_params, we could avoid having to manually set it when vis_params is changed as it would re-calculate automatically for new parameters.

A couple other general questions while I'm thinking about this:

Any thoughts on testing? Unit testing server-side operations is always tricky...
You mentioned adding this functionality to EELeafletTileLayer - do we need to duplicate it for EEFoliumTileLayer to support both mapping backends? If so, any way we can abstract that out to a shared base class?

Answer 6 · 2023-07-08T02:35:10.000Z

Map.ee_layer_dict was introduced a long time ago, and it is used significantly in the source code for accessing the list of ee layers for other purposes. EELeafletTileLayer and EEFoliumTileLayer were introduced by the code refactoring recently. The vis_params attribute of two TilerLayer classes seems unused. We can probably remove it. Also, the vis_params is only tied to a specific layer, unlike Map.ee_layer_dict that gets the entire list of all layers. Users can access the vis_params attribute of any layer through Map.ee_layer_dict[layer_name]['vis_params'].

A question for both of you, probably: I'm a bit confused by the potential duplication between the properties stored in Map.ee_layer_dict and in the EELeafletTileLayer itself. Both seem to store vis_params and ee_object. Would every update to those need to be applied in both places? Just want to make sure I understand the connection between those!

Answer 7 · 2023-07-08T02:43:29.000Z

The vis_params attribute is being removed from EEFoliumTileLayer and EELeafletTileLayer. #1619

Answer 8 · 2023-07-08T02:57:06.000Z

The url_format attribute is the tile URL, which does not need to be changed.

Would it make sense for EELeafletTileLayer.url_format to be converted to a cached_property? By caching it based on the vis_params, we could avoid having to manually set it when vis_params is changed as it would re-calculate automatically for new parameters.

What we need to do is when users invoke the layer visualize GUI, we compute the statistics of the selected layer and add the stats key to the selected layer under Map.ee_layer_dict. Something like this:

{'SRTM DEM': {'ee_object': <ee.image.Image at 0x7f228dd4e170>,
  'ee_layer': EELeafletTileLayer(attribution='Google Earth Engine', max_zoom=24, name='SRTM DEM', options=['attribution', 'bounds', 'detect_retina', 'max_native_zoom', 'max_zoom', 'min_native_zoom', 'min_zoom', 'no_wrap', 'tile_size', 'tms', 'zoom_offset'], url='https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/maps/82302bf76471cbc4d71eafdc4e3c615a-29d3130e1529b0a52f048bf2007db785/tiles/{z}/{x}/{y}'),
  'vis_params': {'min': 0,
   'max': 4000,
   'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5']},
   'stats': {'min': 0.0, 'max': 4000.0, 'sigma1': [a1, b1], 'sigma2': [a2, b2], 'sigma3': [a3, b3], 'pct90': [a4, b4], 'pct98': [a5, b5], 'pct100': [a6, b6]}},
   }

I am not sure if we want to re-compute the stats every time the map zoom level and view region are changed. It is overkilled, making the map unresponsive.