While playing with the Sentinel Playground, trying to trace the Strait of Magellan, I came across a very nice cloud-free Sentinel-2 image of Monte Sarmiento in Tierra del Fuego. Its double peak rises 2,246 m above the dark cold waters of the surrounding fjords and canals. Couple that with the picturesque Schiaparelli glacier that extends down its slopes and almost reaches the sea, but stops just short to form a beautiful milk-watered proglacial lake on the shore, and it’s easy to see why Darwin called this “the most sublime spectacle in Tierra del Fuego”. Darwin watched it from an entirely different perspective, but the satellite view fits his description just as well.
The default Sentinel Hub rendering shown in the image below has some issues, though: the waters and coastal vegetation are too dark, the snow-covered peaks look over-exposed and are all flat white, and there is a bluish tint across the whole image.
I set out to improve the rendering and create something similar to images produced by Pierre Markuse, Antti Lipponen and other masters of the craft. I wanted to achieve a nice color-balanced image using only the Sentinel Hub’s scripting capability, instead of downloading the image and using desktop image processing tools, which are usually employed to achieve such improvements (see guides by Robert Simmon on Planet and Landsat images).
Start with sRGB output
I first choose the date of May 5th, 2016, and then set cloud cover filter to 100% to make sure that I get the image of the full swath, even if some of the neighboring tiles are cloudy. After switching to ‘custom script’ rendering mode, I get the following script:
The script maps reflectances from bands 4, 3 and 2 into red, green and blue components of the output color. The default gain factor is set to 2.5, which corresponds to white-point reflectance of 40% (a pixel with 40% or more reflectance in all three bands will appear white in the resulting image). I will deal with white-point selection later, so I set this factor to 1 for now.
Sentinel Hub renderer converts the returned array into colors of the output image by multiplying the values with
255. Since color components will be interpreted by the web browser as non-linear sRGB numbers and I want to keep working in linear sRGB space, I add the sRGB encoding transfer function as the last step of the calculation:
If you’re wondering why we need this complication: color models are hard. The sRGB color space (standard Red Green Blue) is the most widely used color space, and the default for all web browsers and other internet image exchange. Non-linear transfer function is a peculiarity of the sRGB standard, introduced to make better use of the 256 values that are available for encoding the intensity of light in a given component into a byte. With nonlinear encoding, the browser can render more different shades of dark gray, which the human eye is more sensitive to, at the cost of having fewer shades of white, where we don’t notice the difference that much.
When working with satellite data, the numbers usually represent reflectance, which is directly proportional to the intensity of light, hence the default mapping should be to the linear version of the sRGB color space in which values also represent intensities of light, lest we introduce unintended color transformations.
Subtract the atmosphere for a nicely saturated image
Encoding into proper sRGB values reveals the fact that the image was taken by a sensor orbiting 800 km above ground, looking through about 10 km of atmosphere. The bluish veil that covers the scene is caused by scattering of sunlight by atmospheric particles — it is the same blue thing that we normally call the sky but this time we’re looking at it in the opposite direction.
The effect can be largely removed by simply subtracting some small constant from the reflectance values in each band. To determine the optimal amount and prevent me from darkening the image too much, I added some code to highlight pixels whose components become negative after adjustment.
Using this rendering, I choose the offsets for atmospheric adjustment in such a way that very few pixels get highlighted, and I get as many different colors as possible. As I increase the offsets, the image becomes noticeably less hazy, showing vivid colors and enhanced contrasts.
I want to keep a non-zero blue component over water, so I don’t mind if highlights are mostly red and green there. In the image that I’m working with, appropriate offsets seem to be around 4 %, 7 % and 12 % for red, green and blue bands, respectively. Other images will have different atmospheric conditions, so other offsets should be chosen.
I can simplify the selection of offsets a bit by considering the fact that most of the atmospheric contribution comes from Rayleigh scattering, the amount of which is inversely proportional to the fourth power of wavelength. Applying this on the spectral sensitivities of the three Sentinel-2 bands that I’m using gives ratios of (1.0, 2.0, 3.25) for the (r, g, b) channels, which is pretty close to the ratios I determined empirically in the script above (1.0, 1.75, 3.0).
Atmospheric correction can now be performed by finding the value of a single parameter, with hardly noticeable difference in results — highlights are different, of course, but the difference is barely noticeable in the rest of the image. 🌐
Balance color components to get the perfect shade of white
Once I’ve compensated for the offset in the measured radiance, I add some multiplicative factors to each component, to ensure that white things in the image appear truly white. Again I start by highlighting the pixels that have some component larger than 1 with the following:
As for the selection of factors, I can get pretty nice results by assuming that they are related to the offset
c0 that I’ve determined already: the surface is illuminated by all the light coming from the sun except the part that was scattered and reflected back into space. The illumination of the surface is therefore proportional to
(1 - c0). Once the light is reflected and travels back towards the satellite, some of it gets lost again, which means that the total flux has to be multiplied by
(1 - c0) again. The resulting full atmospheric adjustment code is then simply:
The image shows the highlighted colors to be quite nicely distributed, especially considering the very crude atmospheric model that I’ve derived with some hand-waving and wishful thinking*.
*DISCLAIMER: Please note that any atmospheric model described in this article is very inaccurate (if not entirely wrong) and as such shouldn’t be used for anything other than pretty pictures (and even that is done at your own risk). If you want to learn more about atmospheric effects in satellite images, you can start with the excellent series of blog posts on atmospheric correction by Olivier Hagolle.
As a side note — both the black-point offset
c0 and the white-point factor
c1 can be determined for each channel individually, without assuming anything about the atmosphere. Code that would do that is very similar, but it has 6 parameters to play with instead of just one: 🌐
Stretch contrast to bring out the sunlit slopes
The image above shows many highlighted pixels, where reflectance values are higher than 100%. This often happens on snow-covered sun-facing slopes, especially when the sun’s angle is low, which is usually true for any image taken around the poles, such as this one.
To prevent flattening of these very bright peaks, I add another constant,
max, which I use to scale reflectance values when converting to output color components. Thus I can fine-tune the contrast so that only the very brightest pixels appear white. The resulting image, after setting
max = 3, starts too look really nice — the scaling brings out lovely details in the snowy slopes of the mountain and the cracks and ridges of the glaciers.
Adjust the contrast curve to improve contrast in the dark areas
The stretched image is compelling, but a lot of the contrast was lost in the darker coastal areas, especially in the shadows. To bring that back, I first modify the scaling function and split the linear adjustment function into two parts, forming a piecewise linear function which increases contrast in the dark areas without decreasing it too much in the bright regions.
The result looks nice, but it appears that the sharp break in the adjustment function creates some strange unnatural shading in the snow-covered slopes.
To try and improve that, I replace the straight-line piecewise function with a smoother version, based on a function that is the quotient of two polynomials and matches the slopes of the original at both extreme ends of the reflectance range. The parameters of this function are the same, but the output is much smoother.
The resulting image retains most of the contrast across the brightness range and removes the artifacts caused by the sharp break, which means that I’m almost finished.
Note that other adjustment functions could be used for the same purpose, e.g. a simpler quotient, cubic splines, logarithm, power (gamma) functions etc. I chose this one because it is easy to understand graphically and offers a simple way to fine-tune the slope around 0, where most of the color enhancement is needed.
Increase saturation to make the image prettier
Any nonlinear adjustment curve applied on the individual r, g, b components will inevitably change the saturation of the rendered colors. Concave curves, such as the one I used, decrease the saturation (because larger components are effectively multiplied by a smaller contrast factor), and that makes the image a bit dull.
This can be corrected with a very crude method of saturation enhancement which increases the distance of the three color components from their average. I select a moderate amount of correction to accentuate the light blue color of the glacial ice, while making sure that the rocks and vegetation don’t look too artificial.
The final image
It took me a while to get here, but I think the final image looks pretty good. And the script is not too complicated, so it should be quite useful for lots of other scenes and potentially adjusted to render something other than glaciers.
I recommend using Sentinel Playground for viewing the image in full detail. Make sure to look around and check out the view at different zoom levels — there are plenty of gems to find in that part of the world.
The following script performs all of the steps mentioned in this article, and supports different options of atmospheric compensation:
Bonus for readers that made it this far
Since I started with Darwin, it seems appropriate to finish with him —my colleague Matej Batič suggested to try and render this scene in 3D, to see if we can get closer to Darwin’s viewpoint and see what he saw all those years ago.
Thank you for reading. The final script is also available in our repository of custom scripts. Let me know if you find the scripts and approaches useful, or if you think of ways to improve them. I’m looking forward to seeing pretty images that you create.
I leave you with the full quote from The Voyage of H.M.S. Beagle:
This mountain, which is one of the highest in Tierra del Fuego, has an altitude of 6800 feet. Its base, for about an eighth of its total height, is clothed by dusky woods, and above this a field of snow extends to the summit. These vast piles of snow, which never melt, and seem destined to last as long as the world holds together, present a noble and even sublime spectacle. The outline of the mountain was admirably clear and defined. Owing to the abundance of light reflected from the white and glittering surface, no shadows were cast on any part; and those lines which intersected the sky could alone be distinguished: hence the mass stood out in the boldest relief. Several glaciers descended in a winding course from the upper great expanse of snow to the sea-coast: they may be likened to great frozen Niagaras; and perhaps these cataracts of blue ice are full as beautiful as the moving ones of water. [John van Wyhe, ed. 2002-. The Complete Work of Charles Darwin Online http://darwin-online.org.uk/]