**3. Computation of urban albedo: remote sensing calibration of DART**

Bottom of atmosphere irradiance and exitance are essential terms in the 3D radiative budget of urban and natural landscapes. In the short‐wavelength spectral domain, the ratio of exitance and irradiance is simply the albedo. DART can compute these terms for any time, date, and spectral band Δ*λ*, by using urban database inputs and actual atmosphere conditions that can be derived from satellite images or data from meteorological centers (e.g., ECMWF) or in situ measurements (e.g., AERONET network). The urban database, presented to DART in a \*.obj file format, must contain all buildings of the considered area, as well as the relevant information about urban vegetation in that area (i.e., geographic locations and physical dimensions of trees). It also requires distinction between roofs, walls, and streets, which means each facet in the database has to be registered in one of these groups. DART uses this registration to assign optical properties of the different scene elements. An example of a 3D urban model in nadir and an oblique view, which was built based on the urban database of the city of Basel, is shown in **Figure 7**.

**Figure 7.** DART 3D view of the city of Basel.

The accuracy of DART simulated albedo *ADART* ,Δ*<sup>λ</sup>* and exitance *MDART* ,Δ*<sup>λ</sup>* corresponds with accuracy of the presented optical properties. **Figure 8**, which shows various views of DART simulated images of Basel in natural colors, illustrates importance of specific optical properties. Roofs and walls appear with different colors. These colors depend on the roof and wall optical properties. A major problem is that it is practically impossible to know the spatial distribution of these properties accurately beforehand, due to their spatial variability. In addition, they vary also in time. In consequence, a major objective of the URBANFLUXES project is to derive maps of optical properties per type of Earth surface element. The approach relies on the use of atmospherically corrected remote sensing data (in our example, Landsat‐8 images), in order to compute the 3D radiative budget accurately. Two calibration methods were developed to compute products albedo more accurately: (1) a direct and (2) an iterative calibration.

**Figure 8.** DART simulated images of Basel, Switzerland: (a) pushbroom image with zoom, (b) airborne camera image, (c) satellite RGB image, (d) satellite TIR image.

### **3.1. DART calibration with satellite images**

*LDART,Δλ(xDART, yDART, Ωv)* the per‐pixel simulations of radiance *LDART,Δλ, n(xDART, yDART, Ωv)* and cross section *σn(xDART, yDART, Ωv)* images of each type of scene element n in discreet directions

**3. Computation of urban albedo: remote sensing calibration of DART**

Bottom of atmosphere irradiance and exitance are essential terms in the 3D radiative budget of urban and natural landscapes. In the short‐wavelength spectral domain, the ratio of exitance and irradiance is simply the albedo. DART can compute these terms for any time, date, and spectral band Δ*λ*, by using urban database inputs and actual atmosphere conditions that can be derived from satellite images or data from meteorological centers (e.g., ECMWF) or in situ measurements (e.g., AERONET network). The urban database, presented to DART in a \*.obj file format, must contain all buildings of the considered area, as well as the relevant information about urban vegetation in that area (i.e., geographic locations and physical dimensions of trees). It also requires distinction between roofs, walls, and streets, which means each facet in the database has to be registered in one of these groups. DART uses this registration to assign optical properties of the different scene elements. An example of a 3D urban model in nadir and an oblique view, which was built based on the urban database of the city of Basel, is shown

The accuracy of DART simulated albedo *ADART* ,Δ*<sup>λ</sup>* and exitance *MDART* ,Δ*<sup>λ</sup>* corresponds with accuracy of the presented optical properties. **Figure 8**, which shows various views of DART simulated images of Basel in natural colors, illustrates importance of specific optical properties. Roofs and walls appear with different colors. These colors depend on the roof and wall optical properties. A major problem is that it is practically impossible to know the spatial distribution of these properties accurately beforehand, due to their spatial variability. In addition, they vary also in time. In consequence, a major objective of the URBANFLUXES project is to derive maps of optical properties per type of Earth surface element. The approach relies on the use of atmospherically corrected remote sensing data (in our example, Landsat‐8 images), in order

along the sensor viewing direction *Ωv)*.

234 Sustainable Urbanization

in **Figure 7**.

**Figure 7.** DART 3D view of the city of Basel.

D*irect calibration* is a straightforward method that calibrates DART without assessing the optical properties of the surface elements per pixel [23], while *iterative calibration* uses an iterative procedure in order to spatially derive the optical properties of the scene elements. In theory, this approach gives more accurate results, but is computationally more complex. Both methods are demonstrated on the example of DART simulation of the city of Basel (Switzer‐ land) calibrated using a Landsat‐8 multispectral image, with band specifications presented in **Table 1**, corrected of atmospheric effects by DLR using the ATCOR atmosphere correction model [24].


**Table 1.** Landsat‐8 band specifications.

It is important to note that the image of the 6th Landsat‐8 band was omitted, because of its spectral coincidence with the atmospheric absorption bands. However, this omission has no impact on presented methods as they are designed to work with any type of satellite data, with a number of spectral bands ≥2 and with any spatial resolution.

## *3.1.1. Direct method*

The direct method consists of the following five steps:

**Step 1.** Reflectance images corresponding to the satellite bands used for calibration are simulated. For these simulations, we set the optical properties to a realistic, but not exact value, as this piece of information is unavailable. The spatial resolution of images is equivalent to the size (*xDART* , *yDART* ) of the DART voxels (in this example to 2.5 m).

	- **◦** The corrected image corresponds to the so‐called direct sun illumination. Hence, bottom of atmosphere (BOA) irradiance is direct, without diffuse component. The reflectance is noted ρdd,BOA.
	- **◦** The corrected image corresponds to the actual BOA sun (direct) and atmosphere (diffuse) irradiance. The reflectance is then noted ρsd,BOA. It would be ρhd,BOA (hemispheric‐direct) if the atmosphere irradiance was isotropic.

Here, the atmospherically corrected satellite image corresponds to the case ρsd,BOA. The corresponding DART simulated term of interest is noted:

$$\rho\_{DART,\lambda} \left( \mathbf{x}\_{DART}, \mathbf{y}\_{DART}, \boldsymbol{\Omega}\_{\boldsymbol{S}}, \boldsymbol{E}\_{\boldsymbol{S},\boldsymbol{B}\boldsymbol{M}} \left( \boldsymbol{\Omega}\_{\boldsymbol{S}} \right), \boldsymbol{L}\_{\boldsymbol{a}\boldsymbol{m}} \left( \boldsymbol{\Omega} \right), \boldsymbol{\Omega}\_{\boldsymbol{s}\boldsymbol{a}}, t\_{\boldsymbol{s}\boldsymbol{a}} \right),$$

where—*xDART* , *yDART* : coordinates of a given point in the DART simulated scene,


**Step 2.** Obtained DART image *ρDART, Δλ(xDART, yDART, ΩS, ES, BOA(ΩS), Latm(Ω), Ωsat, tsat)* was georeferenced and resampled to the spatial resolution of the Landsat‐8 image (*xsat*, *ysat*) to ensure exact correspondence in size and geolocation.

**Step 3.** A calibration factor *K*Δ*λ*(*xsat*, *ysat*, *tsat*) is computed per DART pixel (*xsat*, *ysat*). The objective of this factor is to mitigate the use of approximate optical properties. We use the following equation:

$$K\_{\Delta\lambda} \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sat}}, t\_{\text{sat}} \right) = \frac{\rho\_{\text{sat}\Delta\text{ur}, \ \lambda} \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sat}}, \mathbf{\Omega}\_{\text{S}}, \mathbf{\Omega}\_{\text{sat}} \right)}{\rho\_{\text{std}\Delta\text{RRT}, \ \lambda} \left( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sat}}, \mathbf{\Omega}\_{\text{S}}, \mathbf{\Omega}\_{\text{sat}}, t\_{\text{sat}} \right)}. \tag{1}$$

**Step 4.** The map of urban spectral albedo *A*Δ*<sup>λ</sup>*, is calculated per spectral band Δ*λ* with:

$$\begin{split} \boldsymbol{A}\_{\boldsymbol{\Delta}\boldsymbol{\lambda}} \Big( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, \boldsymbol{\Omega}\_{\text{S}}, \boldsymbol{E}\_{\text{S,BOL,\lambda\lambda}} \Big( \boldsymbol{\Omega}\_{\text{S}} \Big), \boldsymbol{L}\_{\text{ann},\text{\lambda}\boldsymbol{\lambda}} \Big( \boldsymbol{\Omega} \Big), \boldsymbol{t}\_{\text{sur}} \Big) &= \boldsymbol{K}\_{\boldsymbol{\Delta}\boldsymbol{\lambda}} \Big( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, \boldsymbol{t}\_{\text{sur}} \Big), \\ \boldsymbol{A}\_{\boldsymbol{\Delta}\boldsymbol{\mathcal{R}}\boldsymbol{\Gamma},\boldsymbol{\lambda}\boldsymbol{\lambda}} \Big( \mathbf{x}\_{\text{sur}}, \mathbf{y}\_{\text{sur}}, \boldsymbol{\Omega}\_{\text{S}}, \boldsymbol{E}\_{\text{S,BOL,\lambda}} \Big( \boldsymbol{\Omega}\_{\text{S}} \Big), \boldsymbol{L}\_{\text{ann},\text{\lambda}} \Big( \boldsymbol{\Omega} \Big), \boldsymbol{t}\_{\text{sur}} \Big), \end{split} \tag{2}$$

where *ADART* ,Δ*<sup>λ</sup>* is the albedo computed by DART, for spectral interval Δ*λ*:

It is important to note that the image of the 6th Landsat‐8 band was omitted, because of its spectral coincidence with the atmospheric absorption bands. However, this omission has no impact on presented methods as they are designed to work with any type of satellite data, with

**Step 1.** Reflectance images corresponding to the satellite bands used for calibration are simulated. For these simulations, we set the optical properties to a realistic, but not exact value, as this piece of information is unavailable. The spatial resolution of images is equivalent to the

**•** The DART reflectance image of interest *ρDART* ,Δ*λ*(*xDART* , *yDART* , Ω*sat*) is simulated in the viewing direction of the available satellite image. It is important to note that this approach can also work with an image which has not been atmospherically corrected. In this case, the simulated image corresponds to a "top of atmosphere" (TOA) data in a "direct‐direct" configuration, for which the incident and exiting radiative fluxes are following a single direction (ρdd,TOA). Classical atmosphere correction methods lead to two types of images:

**◦** The corrected image corresponds to the so‐called direct sun illumination. Hence, bottom of atmosphere (BOA) irradiance is direct, without diffuse component. The reflectance is

**◦** The corrected image corresponds to the actual BOA sun (direct) and atmosphere (diffuse) irradiance. The reflectance is then noted ρsd,BOA. It would be ρhd,BOA (hemispheric‐direct) if

Here, the atmospherically corrected satellite image corresponds to the case ρsd,BOA. The

**Step 2.** Obtained DART image *ρDART, Δλ(xDART, yDART, ΩS, ES, BOA(ΩS), Latm(Ω), Ωsat, tsat)* was georeferenced and resampled to the spatial resolution of the Landsat‐8 image (*xsat*, *ysat*) to

( *xy E L t* , ,W W WW , , ( ), ,, ( ) *sat sat*)

*DART DART DART S S BOA S atm* ,Δ

where—*xDART* , *yDART* : coordinates of a given point in the DART simulated scene,

a number of spectral bands ≥2 and with any spatial resolution.

size (*xDART* , *yDART* ) of the DART voxels (in this example to 2.5 m).

The direct method consists of the following five steps:

the atmosphere irradiance was isotropic.

r

— *L atm*(Ω): atmosphere radiance,

— *tsat*: time of the satellite acquisition.

corresponding DART simulated term of interest is noted:

l

— Ω*<sup>S</sup>* , Ω*sat*: sun and satellite view direction, respectively,

ensure exact correspondence in size and geolocation.

— *ES* ,*BOA*: sun irradiance at the bottom of the atmosphere (BOA),

*3.1.1. Direct method*

236 Sustainable Urbanization

noted ρdd,BOA.

$$\begin{split} \boldsymbol{A}\_{\text{DART},\boldsymbol{\lambda}}\left(\mathbf{x}\_{\text{DART}},\mathbf{y}\_{\text{DART}},\boldsymbol{\Omega}\_{\text{S}},\boldsymbol{E}\_{\text{S}\_{\text{BOM}}}\left(\boldsymbol{\Omega}\_{\text{S}}\right),\boldsymbol{L}\_{\text{ann}}\left(\boldsymbol{\Omega}\right),\boldsymbol{t}\_{\text{sat}}\right) &= \\ \frac{\rho\_{\text{d}\boldsymbol{\theta}}\boldsymbol{E}\_{\text{S},\boldsymbol{\text{BOM}},\boldsymbol{\lambda}}\,+\,\big[\rho\_{\text{d}\boldsymbol{\theta}\text{-}\boldsymbol{\Lambda}\boldsymbol{\lambda}}\big(\boldsymbol{\Omega}\big)\boldsymbol{L}\_{\text{ann\\_\text{AL}}}\left(\boldsymbol{\Omega}\right)\boldsymbol{\cos}\left(\boldsymbol{\theta}\right)d\boldsymbol{\Omega}\boldsymbol{\Omega}}{\boldsymbol{E}\_{\text{S},\boldsymbol{\text{BOM}},\boldsymbol{\lambda}}\,+\,\big[\boldsymbol{L}\_{\text{ann\\_\text{AL}}}\big(\boldsymbol{\Omega}\big)\boldsymbol{\cos}\left(\boldsymbol{\theta}\right)d\boldsymbol{\Omega}}.\end{split}\tag{3}$$

In (3), *θ* is the zenith angle corresponding to the viewing direction Ω.

**Step 5.** Desired urban surface albedo *A* is calculated as the integral of all the spectral albe‐ dos *A*Δ*λ*(*xsat*, *ysat*, Ω*<sup>S</sup>* , *ES* ,*BOA*,Δ*λ*(Ω*<sup>S</sup>* ), *L atm*,Δ*λ*(Ω), *tsat*) across the entire spectrum, weighted by the BOA irradiance *ES* ,*BOA*,Δ*λ*(Ω*<sup>S</sup>* ).

**Figure 9.** Landsat‐8 atmospherically corrected image for an NIR band (864.6 nm) over the city of Basel, acquired on the 24th April, 2015.

Here, we consider an application of the method to the city of Basel. **Figure 9** shows the atmospherically corrected Landsat‐8 band 5 (**Table 1**) image.

The original DART image of the same spectral band without calibration, at the simulation resolution of 2.5 m, is shown in **Figure 10**. This image has been georeferenced using the spatial information retrieved from the Landsat‐8 satellite image.

**Figure 10.** DART reflectance image for an NIR band (864.6 nm), over the city of Basel, corresponding to a Landsat‐8 image, with erroneous (but realistic) optical properties used as input.

Image in **Figure 10** was resampled to the Landsat‐8 native spatialresolution using ILWIS image processing software. After computation of the calibration factor (step 3) and the albedo per spectral band, we integrated the final broadband albedo product as shown in **Figure 11**. In **Figure 11**, the color bar indicates the albedo values, and the axes simply indicate the pixels positions.

**Figure 11.** Urban surface albedo image computed using the direct calibration method for the city of Basel, with a Landsat‐8 multispectral image.

## *3.1.2. Iterative method*

The original DART image of the same spectral band without calibration, at the simulation resolution of 2.5 m, is shown in **Figure 10**. This image has been georeferenced using the spatial

**Figure 10.** DART reflectance image for an NIR band (864.6 nm), over the city of Basel, corresponding to a Landsat‐8

Image in **Figure 10** was resampled to the Landsat‐8 native spatialresolution using ILWIS image processing software. After computation of the calibration factor (step 3) and the albedo per spectral band, we integrated the final broadband albedo product as shown in **Figure 11**. In **Figure 11**, the color bar indicates the albedo values, and the axes simply indicate the pixels

**Figure 11.** Urban surface albedo image computed using the direct calibration method for the city of Basel, with a

information retrieved from the Landsat‐8 satellite image.

image, with erroneous (but realistic) optical properties used as input.

positions.

238 Sustainable Urbanization

Landsat‐8 multispectral image.

Iterative calibration method computes the actual optical properties for each type of element in the scene per pixel of the satellite image or per group of pixels, which provides the desirable spatial variability in optical properties of present elements. If N types of elements are present in a single pixel, the reflectance of this pixel is formed by a number of different optical properties, that is, by N unknowns. Therefore, a set of equations is required to resolve the contribution of every surface element in a final reflectance value, which involves the image per type of scene elements described in Section 2.2. Taking into account a number of pixels, one obtains a system of equations, where a number of equations are equal or greater than the number of unknowns.

Consecutive steps of this approach are presented below, with N being the number of types of scene elements present in the studied urban area and parameter k being the order of the iteration initialized at k = 1:

**Step 1.** Regular grid is laid over the satellite image, with a mesh size equal to M M satellite pixels and with *M* > *N* . This ensures that each cell of the mesh contains a number of satel‐ lite pixels equal or larger than N.

**Step 2.** DART simulates the total radiance image *L DART* ,Δ*λ*(*xDART* , *yDART* , Ω*sat*) of the scene and also the radiance images per scene element n *L DART* ,Δ*λ*,*n*(*xDART* , *yDART* , Ω*sat*) in the satellite viewing direction Ω*sat*. Consequently, each surface element n viewed by a DART pixel d along the satellite viewing direction Ω*sat* has a reflectance value equal to *ρn*,*<sup>d</sup>* ,Δ*<sup>λ</sup> <sup>k</sup>* (*xDART* , *yDART* ).

It is important to note that reflectance and radiance are proportional terms, linked by the BOA sun irradiance:

$$\left(\rho\_{DART\lambda}\left(\mathbf{x}\_{DART},\mathbf{y}\_{DART},\boldsymbol{\Omega}\_{\text{sat}}\right) = \frac{\pi \boldsymbol{L}\_{DART,\text{\AA}\boldsymbol{\ell}}\left(\mathbf{x}\_{DART},\mathbf{y}\_{DART},\boldsymbol{\Omega}\_{\text{sat}}\right)}{\boldsymbol{E}\_{DART,\text{ROA},\text{\AA}\boldsymbol{\ell}}}.\tag{4}$$

The first iteration k = 1 assumes that the reflectance of an urban elements *ρn*,*<sup>d</sup>* ,Δ*<sup>λ</sup> <sup>k</sup>* (*xDART* , *yDART* ) is constant across the entire DART scene. Its plausible value is provided by a reasonable first guess, taken, for instance, from a recent airborne spectral data. Next step is to compute the reflectance value *ρn*,*<sup>d</sup>* ,Δ*<sup>λ</sup> <sup>k</sup>* <sup>+</sup><sup>1</sup> (*xDART* , *yDART* ) of the urban surface elements, which we will be used by DART in the follow‐up iteration k + 1. Objective of each iteration is to bring the DART reflectance closer to the reflectance of a real satellite image. In a DART pixel d, mean irradi‐ ance of an element type n is given by:

$$E\_{\text{DART},\lambda,\eta,d}\left(\mathbf{x}\_{\text{DART}},\mathbf{y}\_{\text{DART}}\right) = \frac{\pi \, L\_{\text{DART},\lambda,\eta}\left(\mathbf{x}\_{\text{DART}},\mathbf{y}\_{\text{DATT}}\right)}{\rho\_{\text{s}}^{\text{k}}\left(\mathbf{x}\_{\text{DATT}},\mathbf{y}\_{\text{DATT}}\right)} \cdot \frac{\Delta\mathbf{x}\_{\text{DATT}} \,\Delta\mathbf{y}\_{\text{DATT}} \,\cos\theta\_{\text{sat}}}{\sigma\_{\text{s},d}\left(\mathbf{x}\_{\text{DATT}},\mathbf{y}\_{\text{DATT}},\,\,\Omega\_{\text{sat}}\right)},\tag{5}$$

where *θsat* is the zenith angle of the satellite viewing direction and *σn*,*<sup>d</sup>* (*xDART* , *yDART* , Ω*sat*) is the cross section of the element of type n, which is viewed by the DART pixel d along the satellite viewing direction Ω*sat*. It is also an output of the DART model. The ratio *∆xDART* .*∆ yDART* .*cosθsat <sup>σ</sup>n*,*<sup>d</sup>* (*xDART* , *yDART* , <sup>Ω</sup>*sat*) had to be introduced because the radiance of any DART pixel is simulat‐ ed per effective square meter of the pixel area *∆xDART* .*∆yDART* .*cosθsat*, while the radiance of the element n is expressed per effective square meter of the area of the surface element itself *σn*,*<sup>d</sup>* (*xDART* , *yDART* , Ω*sat*). Both radiances are equal, if only a single element is present in the considered pixel.

**Step 3.** DART radiance *L DART* ,Δ*λ*(Ω*sat*) and element *L DART* ,Δ*λ*,*n*(Ω*sat*) images are resampled in the satellite image spatial resolution (Δ*xsat*, Δ*ysat*). If selecting an appropriate spatial resolu‐ tion of aDARTsimulation, one can make the assumption that each satellite pixel m is composed of an integer number *D* <sup>2</sup> of DART pixels. Therefore, for any given satellite pixel m, where d is the index ofthe DART pixels in m (*d* ∈ 1 *D* <sup>2</sup> ) and n is the index ofthe elementtype (*n* ∈ 1 *N* ), one can define:

**•** Radiance as:

$$L\_{\text{DARTA},\lambda,m} \left( \mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}}, \ \ \ \ \ \ \ \ \mathbf{D}\_{\text{sat}} \right) = \frac{\sum\_{d=1}^{D^2} L\_{\text{DARTA},\lambda,d}}{D^2},\tag{6}$$

**•** Radiance of scene element n as:

$$L\_{\text{DART}\downarrow\lambda,\eta,m} \left( \chi\_{\text{sart}}, \mathcal{Y}\_{\text{sag}\,\Omega\_{\text{an}}} \right) = \frac{\sum\_{d=1}^{D^2} L\_{\text{DART}\downarrow\lambda,\eta,d}}{D^2},\tag{7}$$

**•** Irradiance of elements of type n as:

$$E\_{\text{DART},\lambda,\eta,d} \left(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}}\right) = \frac{\sum\_{d=1}^{D^2} E\_{\text{DART},\lambda,\eta,d} \left(\mathbf{x}\_{\text{DART}}, \mathbf{y}\_{\text{DART}}\right) \sigma\_{\eta,d} \left(\mathbf{x}\_{\text{DATT}}, \mathbf{y}\_{\text{DATT}}\right)}{\sum\_{d=1}^{D^2} \sigma\_{\eta,d} \left(\mathbf{x}\_{\text{DATT}}, \mathbf{y}\_{\text{DATT}}\right)},\tag{8}$$

and

**•** Cross section of element n viewed by pixel m as:

Remote Sensing Studies of Urban Canopies: 3D Radiative Transfer Modeling http://dx.doi.org/10.5772/63886 241

$$
\sigma\_{n,m} \left( \mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}} \right) = \sum\_{d=1}^{D^2} \sigma\_{n,d} \left( \mathbf{x}\_{DART}, \mathbf{y}\_{DART} \right). \tag{9}
$$

For a first approximation, Eq. (8) can be written similarly as Eq. (5):

( ) ( )

r

p

,Δ , ,

*∆xDART* .*∆ yDART* .*cosθsat*

240 Sustainable Urbanization

considered pixel.

one can define:

**•** Radiance as:

of an integer number *D* <sup>2</sup>

**•** Radiance of scene element n as:

**•** Irradiance of elements of type n as:

,Δ , ,

and

l

l

*DART n d DART DART k*

,Δ ,

*L xy x y cos E xy*

l

( ) ( )

D D <sup>=</sup> <sup>W</sup> (5)

q

*DART n DART DART DART DART sat*

 s

*xy xy*

. , . . , ., , , ,

where *θsat* is the zenith angle of the satellite viewing direction and *σn*,*<sup>d</sup>* (*xDART* , *yDART* , Ω*sat*) is the cross section of the element of type n, which is viewed by the DART pixel d along the satellite viewing direction Ω*sat*. It is also an output of the DART model. The ratio

*<sup>σ</sup>n*,*<sup>d</sup>* (*xDART* , *yDART* , <sup>Ω</sup>*sat*) had to be introduced because the radiance of any DART pixel is simulat‐ ed per effective square meter of the pixel area *∆xDART* .*∆yDART* .*cosθsat*, while the radiance of the element n is expressed per effective square meter of the area of the surface element itself *σn*,*<sup>d</sup>* (*xDART* , *yDART* , Ω*sat*). Both radiances are equal, if only a single element is present in the

**Step 3.** DART radiance *L DART* ,Δ*λ*(Ω*sat*) and element *L DART* ,Δ*λ*,*n*(Ω*sat*) images are resampled in the satellite image spatial resolution (Δ*xsat*, Δ*ysat*). If selecting an appropriate spatial resolu‐ tion of aDARTsimulation, one can make the assumption that each satellite pixel m is composed

the index ofthe DART pixels in m (*d* ∈ 1 *D* <sup>2</sup> ) and n is the index ofthe elementtype (*n* ∈ 1 *N* ),

,Δ , <sup>2</sup> , , ,

,Δ , , , Ω <sup>2</sup> , , *sat*

( ) ( ) ( )

,. , , ,

*DART n d DART DART n d DART DART d*

*n d DART DART d*

, 1

s

2 ,Δ,, , 1

=

2

*L*

*D*

2

*L*

*D*

,Δ , 1

,Δ,, 1

( )

s

,

*x y*

l

*DART n d d*

l

<sup>=</sup> W = å (6)

<sup>=</sup> <sup>=</sup> å (7)

<sup>å</sup> (8)

*DART d d*

( )

*L xy <sup>D</sup>*

( )

*L xy <sup>D</sup>*

*E xy xy E xy*

l

*DART m sat sat sat*

*DART n m sat sat*

2

=

*D*

l

*DART n d sat sat D*

**•** Cross section of element n viewed by pixel m as:

= å

l

,

of DART pixels. Therefore, for any given satellite pixel m, where d is

*n DART DART n d DART DART sat*

$$E\_{\text{DART},\lambda,\text{\tiny{u}},\text{\tiny{u}}}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}}\right) = \frac{\pi \, L\_{\text{DART},\lambda,\text{\tiny{u}}}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}},\,\mathbf{\Omega}\_{\text{sat}}\right)}{\rho\_{\text{u},\text{u}}^{\text{k}}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}}\right)} \cdot \frac{\Delta\mathbf{x}\_{\text{sat}}\,\Delta\mathbf{y}\_{\text{sat}}\,\cos\theta\_{\text{sat}}}{\sigma\_{\text{u},\text{u}}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}},\,\mathbf{\Omega}\_{\text{sat}}\right)}.\tag{10}$$

In this expression, *ρn*,*<sup>m</sup> <sup>k</sup>* (*xsat*, *ysat*) is the reflectance of the scene element of type n in the satellite pixel m at iteration k, as computed in the last iteration or equal to the initial value of k = 1. It assumes that the reflectance value of each DART pixel d corresponding to the satellite pixel m is constant.

**Step 4.** Resampled DART radiance is compared with the satellite radiance of all considered image pixels. If the difference between the two is acceptable, the procedure is terminated, and the desired urban albedo is the one computed by DART in the iteration of its simulation. Similarly to the direct calibration method, the albedo product is computed from all the spectral albedos. However, if the difference between the two radiances is according to user require‐ ments too large, the computation enters in a next step.

**Step 5.** The actual calibration of DART for each cell of the mesh defined in step 1 is conduct‐ ed for all *M* <sup>2</sup> satellite pixels and in each cell u of the mesh. Since *M* <sup>2</sup> > *N* , the number of pixels m analyzed in each cell u is larger than the number of different surface element types. Therefore, a deconvolution could be applied to retrieve the optical properties of each surface element present in the studied cell. Each cell u contains *M* <sup>2</sup> satellite pixels m, leading to a system of *M* <sup>2</sup> equations verifying if the DART image and the satellite image are equal:

$$L\sum\_{\mathfrak{u}=1}^{N}L\_{\text{DART},\mathfrak{u},\mathfrak{u},\mathfrak{u}}\left(\mathbf{x}\_{\text{sar}},\mathbf{y}\_{\text{sar}},\boldsymbol{\Omega}\_{\text{sar}}\right) = L\_{\text{sar},\mathfrak{u},\mathfrak{u}}\left(\mathbf{x}\_{\text{sar}},\mathbf{y}\_{\text{sar}},\boldsymbol{\Omega}\_{\text{sar}}\right), \forall m \in \left[1\,M^2\right].\tag{11}$$

Obviously, the verification is negative if the two images differ. At iteration k, the DART radiance values *L DART* ,Δ*λ*,*n*,*m*(*xsat*, *ysat*, Ω*sat*) are computed with inaccurate approximated optical properties *ρn*,*<sup>u</sup> <sup>k</sup>* (*xsat*, *ysat*). If we define *<sup>L</sup> DART* ,Δ*λ*,*n*,*<sup>m</sup>* ' (*xsat*, *ysat*, Ω*sat*) as the DART radiance value computed from *ρn*,*<sup>u</sup> k* +1 (*xsat*, *ysat*), and if we accept the fact that radiance values are proportional to reflectance values [Eq. (10)], then we can write:

$$\frac{L\_{\text{DARTA},\text{A},\text{u},\text{u}}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}},\boldsymbol{\Omega}\_{\text{sat}}\right)}{\boldsymbol{\rho}^{k}\_{\text{u},\text{u}}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}}\right)} = \frac{L^{\top}\_{\text{DARTA},\text{A},\text{u},\text{u}}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}},\boldsymbol{\Omega}\_{\text{sat}}\right)}{\boldsymbol{\rho}^{k\cdot\text{1}}\_{\text{u},\text{u}}\left(\mathbf{x}\_{\text{sat}},\mathbf{y}\_{\text{sat}}\right)}.\tag{12}$$

Consequently, we can rewrite the system of Eq. (11) to a system of *M* <sup>2</sup> equations, in which the unknowns are the expected reflectance values *ρn*,*<sup>u</sup> k* +1 (*xsat*, *ysat*):

$$\sum\_{n=1}^{N} \frac{\rho\_{n,\boldsymbol{u}}^{k+1}(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}})}{\rho\_{n,\boldsymbol{u}}^{k}(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}})} \cdot L\_{\text{DART},\lambda,\mathbf{u},\mathbf{u}}(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}}) = L\_{\text{auța }\lambda,\mathbf{u}}(\mathbf{x}\_{\text{sat}}, \mathbf{y}\_{\text{sat}}), \forall m \in \left[1M^{2}\right]. \tag{13}$$

Resolving this system for each cell u of the mesh leads to the new desired reflectance values *ρn*,*<sup>u</sup> k* +1 (*xsat*, *ysat*), which is used in the follow‐up iteration starting in step 2. This system will find no solutions, if one of the satellite pixels of the cell u contains N types of elements, but all other pixels contain just a single unique element of the same type. In this case, the solution needs a new adaptation of the mesh grid size.

The iterative method has been executed on the same case as the direct method, but so far we retrieved the new optical properties after only a single iteration. We show in **Figure 12** the preliminary result assessing the operational functioning of the method. On the left, we see the spectral Landsat‐8 image used for the calibration, for the 5th spectral band (864.6 nm). In the center is the corresponding initial DART reflectance image, resampled to the satellite resolu‐ tion. On the bottom is the new DART reflectance image, with updated optical properties for each element, on the new grid.

**Figure 12.** (Top left) Landsat‐8 image used for the calibration (NIR band). (Top right) DART image with initial optical properties (NIR band), resampled in the satellite resolution. (Bottom) DART image with updated optical properties (NIR band).

After one iteration, the computed mean relative error between the DART image and the Landsat‐8 image went from *rel* =0.3519 to *rel*,*new* =0.0998. It improved by a factor of 3.5. A clear advantage of this method compared to the direct one is the actual computation of new optical properties, which we are able to use in the simulation of other satellite images, whereas the direct calibration only works for a single image. A better accuracy is also expected after several iterations.

#### **3.2. Albedo derived from climatological data**

Consequently, we can rewrite the system of Eq. (11) to a system of *M* <sup>2</sup>

,Δ , , ,Δ , 1

l

*k DART n m sat sat sat m sat sat n*

( ) ( ) ( ) <sup>1</sup>

*k* +1

, . , , , <sup>1</sup> . ,

<sup>=</sup> <sup>=</sup> " Î é ù å ë û (13)

 l

, 2

Resolving this system for each cell u of the mesh leads to the new desired reflectance values

The iterative method has been executed on the same case as the direct method, but so far we retrieved the new optical properties after only a single iteration. We show in **Figure 12** the preliminary result assessing the operational functioning of the method. On the left, we see the spectral Landsat‐8 image used for the calibration, for the 5th spectral band (864.6 nm). In the center is the corresponding initial DART reflectance image, resampled to the satellite resolu‐ tion. On the bottom is the new DART reflectance image, with updated optical properties for

**Figure 12.** (Top left) Landsat‐8 image used for the calibration (NIR band). (Top right) DART image with initial optical properties (NIR band), resampled in the satellite resolution. (Bottom) DART image with updated optical properties

(*xsat*, *ysat*), which is used in the follow‐up iteration starting in step 2. This system will find no solutions, if one of the satellite pixels of the cell u contains N types of elements, but all other pixels contain just a single unique element of the same type. In this case, the solution needs a

*x y L xy L xy m M*

(*xsat*, *ysat*):

unknowns are the expected reflectance values *ρn*,*<sup>u</sup>*

( )

*n u sat sat*

*x y*

,

+

r

242 Sustainable Urbanization

r

*ρn*,*<sup>u</sup> k* +1

(NIR band).

*<sup>k</sup> <sup>N</sup> n u sat sat*

new adaptation of the mesh grid size.

each element, on the new grid.

equations, in which the

Alternative methodology to derive urban albedo images has been adapted for a case when no satellite data of the area of interest are available. DART images are, in this case, simulated with optical properties of surface elements calibrated from the last available satellite image and with available atmospheric illumination conditions of direct sun irradiance *ES* ,*BOA*,Δ*<sup>λ</sup>* and diffuse sky irradiance *Eatm*,Δ*<sup>λ</sup>*. This information can be obtained from several sources including ECMWF (www.ecmwf.int), Meteo France (www.meteofrance.com), and other in situ mete‐ orological sensor networks.

In this alternative approach, DART calculates the so‐called white sky and black sky albedos.

The white sky albedo *ADART* ,*white sky*, <sup>Δ</sup>*λ*(*xsat*, *ysat*, *ES* ,*BOA*(Ω*<sup>S</sup>* ) =0, *<sup>L</sup> atm* <sup>=</sup> *Eatm <sup>π</sup>* , *tsat*) is computed for an irradiance coming solely from the atmosphere, that is, without any direct sunlight contri‐ bution. The black sky albedo *ADART, white sky, Δλ(xsat, ysat, ES, BOA(ΩS)=0, Latm=0, tsat)* is computed for the direct sunlight only. The desired final albedo product is then computed as a combination of both components:

$$A\_{\Delta \mathcal{k}} \left( \chi\_{\rm aer}, \mathcal{y}\_{\rm aer}, \mathfrak{Q}\_{\rm S}, E\_{s, \rm BOA} \left( \mathfrak{Q}\_{\rm S} \right), E\_{\rm am}, t\_{\rm air} \right) = E\_{s, \rm BOA} \left( \mathfrak{Q}\_{\rm S} \right). \tag{14}$$

$$A\_{\rm DART, \rm bluck \ sky, \rm A\mathcal{J}} + E\_{\rm am} \cdot A\_{\rm DART, \rm white \ sky, \rm A\mathcal{J}}$$

The black sky albedo must be pre‐computed for a set of sun directions such that the black sky albedo of an actual configuration (i.e., sun direction/date) can be derived by interpolation on pre‐computed black sky albedos.

#### **3.3. Computation of thermal exitance images**

The calibration of the DART model in the thermal infrared domain, using satellite thermal infrared images, follows the same kind of approach. Indeed, the method will rely on spatial‐ ly setting up systems of equations over groups of pixels in the satellite image, to solve for the desired parameters. However, an additional difficulty comes from the fact that there are now 2 unknowns for each single measurement: the temperature *T* and the emissivity of the urban surface elements. Furthermore, the treatment of satellite pixels from thermal imagery only give information on the value of the product × *L* (*T* ), where *L* (*T* ) is Planck's law for temperature (*T* ). Hence, separating the variables by simply taking more pixels to create equations for the system is not adequate. The adopted strategy is to consider 2 thermal infrared images that correspond to 2 close spectral bands, with the assumption that emissivity is the same forthese 2 spectral bands. This approach leads to the determination of the thermodynamic temperature and emissivity per surface element of type n. Then by considering the other spectral bands, if the satellite image has more than 2 thermal infrared bands, it leads to the determination of the emissivity of each type of surface element in those spectral bands. In the URBANFLUXES project, we will work with satellite images that have several thermal infrared bands, as shown in **Table 2**.


**Table 2.** Landsat‐8 thermal infrared band specifications.
