import numpy as np
from numpy import cos, sin
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import gdal
%matplotlib inline
%cd D:\Data\Sentinel\2A\S2A_tile_20150829_31TGK_0
D:\Data\Sentinel\2A\S2A_tile_20150829_31TGK_0
toa_refl = dn / quanitification factor
SMAC with SRTM DEM and GEOS-5 atmospheric conditions
Angles (band averaged) from S2A metadata.xml
Calculated using SPCTRL2, with SRTM DEM and GEOS-5 atmospheric conditions
def load_GeoRaster(filename, xoff=None, yoff=None, xsize=None, ysize=None):
"""Helper function to load a subset
Returns: data, extent, gt, proj"""
ds = gdal.Open(filename)
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
if xoff is None:
data = ds.ReadAsArray()
extent = [gt[0], gt[0] + (gt[1] * ds.RasterXSize), gt[3] + (gt[5] * ds.RasterYSize), gt[3]]
else:
data = ds.ReadAsArray(xoff, yoff, xsize, ysize)
extent = [gt[0] + (gt[1] * xoff), gt[0] + ((xoff+xsize) * gt[1]),
gt[3] + (gt[5] * (yoff+ysize)), gt[3] + (yoff * gt[5])]
ds = None
return data, extent, gt, proj
norm = lambda x,vmin=0,vmax=0.13: np.clip((x - vmin) / vmax, 0, 1)
roll = lambda x: np.rollaxis(x,0,3)
# reflectances
subset_props = dict(xoff=3000, yoff=2000, xsize=500, ysize=1000)
#subset_props = dict(xoff=4500, yoff=3500, xsize=500, ysize=1000)
bgr_boa, extent, gt, proj = load_GeoRaster(r'rgb_stack_boa_refl.tif', **subset_props)
rgb_boa = bgr_boa[::-1,...] # bgr to rgb
bgr_toa, extent, gt, proj = load_GeoRaster(r'rgb_stack.vrt', **subset_props)
rgb_toa = bgr_toa[::-1,...] / 10000 # # bgr to rgb & dn to toa refl
# direct and diffuse radiation from spectral
rdif, extent, gt, proj = load_GeoRaster(r'spctrl_stack_10m_dif_tilted.tif', **subset_props)
rdif = rdif[::-1,...]
rdir, extent, gt, proj = load_GeoRaster(r'spctrl_stack_10m_dir_tilted.tif', **subset_props)
rdir = rdir[::-1,...]
d2r = np.deg2rad(1)
slope, extent, gt, proj = load_GeoRaster(r'srtm_slope.tif', **subset_props)
aspect, extent, gt, proj = load_GeoRaster(r'srtm_aspect.tif', **subset_props)
slope *= d2r
aspect *= d2r
Subset somewhere in the French Alps, south-east of Grenoble:
https://www.google.com/maps/@44.886034,5.9197751,13z/data=!3m1!1e3
# prep data for plotting
rgb_toa_plot = norm(roll(rgb_toa))
rgb_boa_plot = norm(roll(rgb_boa))
# plot
fig, axs = plt.subplots(1,2, figsize=(12,12), subplot_kw=dict(xticks=[], yticks=[]))
fig.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0)
axs[0].imshow(rgb_toa_plot)
axs[0].text(0.01, .99, 'TOA refl', transform=axs[0].transAxes, va='top', color='k', weight='bold', size=20)
axs[0].axis('off')
axs[1].imshow(rgb_boa_plot)
axs[1].axis('off')
axs[1].text(0.01, .99, 'BOA refl', transform=axs[1].transAxes, va='top', color='k', weight='bold', size=20)
<matplotlib.text.Text at 0x75227f0>
sza = 37.36
saa = 156.16
sza *= d2r
saa *= d2r
# incedence angle
i = np.arccos(cos(sza)*cos(slope)+sin(sza)*sin(slope)*cos(saa-aspect))
rtot = rdir + rdif
ground_albedo = 0.2
fsky = (1 + cos(slope)) / 2
fground = (1 - cos(slope)) / 2
corr_factor = rtot / (rdir * (cos(i) / cos(sza)) + rdif * fsky + rtot * fground * ground_albedo)
rgb_boa_corr = rgb_boa * corr_factor
rgb_boa_corr_plot = norm(roll(rgb_boa_corr))
# plot
fig, axs = plt.subplots(1,2, figsize=(12,12), subplot_kw=dict(xticks=[], yticks=[]))
fig.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0)
axs[0].imshow(rgb_boa_plot)
axs[0].text(0.01, .99, 'BOA refl', transform=axs[0].transAxes, va='top', color='k', weight='bold', size=20)
axs[0].axis('off')
axs[1].imshow(rgb_boa_corr_plot, cmap=plt.cm.Greys_r)
axs[1].axis('off')
axs[1].text(0.01, .99, 'BOA refl corr', transform=axs[1].transAxes, va='top', color='k', weight='bold', size=20)
axs[1].add_artist(Ellipse(xy=[.46,.56], width=.25, height=.1, angle=10, transform=axs[1].transAxes,
facecolor='none', edgecolor='red', linewidth=2.))
axs[1].add_artist(Ellipse(xy=[.65,.77], width=.2, height=.1, angle=0, transform=axs[1].transAxes,
facecolor='none', edgecolor='red', linewidth=2.))
<matplotlib.patches.Ellipse at 0xd8b4a90>