Remote Sensing Studies of Urban Canopies: 3D Radiative Transfer Modeling

Abstract

Need for better understanding and more accurate estimation of radiative fluxes in urban environments, specifically urban surface albedo and exitance, motivates development of new remote sensing and three‐dimensional (3D) radiative transfer (RT) modeling methods. The discrete anisotropic radiative transfer (DART) model, one of the most comprehensive physically based 3D models simulating Earth/atmosphere radiation interactions, was used in combination with satellite data (e.g., Landsat‐8 observations) to better parameterize the radiative budget components of cities, such as Basel in Switzerland. After presenting DART and its recent RT modeling functions, we present a methodological concept for estimating urban fluxes using any satellite image data.

1. Introduction

In today's world, concerns about sustainable energy production, related environmental issues, and their impact on urban agglomerations and their citizens are of high importance. In this context, the functioning of urban environments, including climate‐alternating interactions between atmosphere and human civilization, is studied globally. Remote sensing is a powerful tool that is increasingly used in such studies due to improvements achieved in sensor technology, as well as modeling of remote sensing measurements with respect to three‐dimensional (3D) radiative and energy budget. One of these studies is the H2020 project URBANFLUXES (http://urbanfluxes.eu), which aims to develop methods estimating the anthropogenic heat flux (Q_F) of urban environments by employing remote sensing data [1]. In other words, its goal is to estimate the impact of human urban activities on the energy budget of city using satellite images. Important parts of the surface energy balance computation are the 3D radiative budget components specifically as urban surface albedo including thermal exitance. However, no remote sensing model able to simulate accurately spatial distribution of urban spectral albedo and exitance has been previously available. Three conditions must be fulfilled in order to achieve an acceptable solution of an urban radiative budget simulation:

1. – The model must consider explicitly the 3D architecture of urban environments and simulate radiance images and radiative budget of urban environment. Hence, apart from physical modeling considerations, the model must be able to work with urban databases, including spatial information on vegetation and digital elevation model.

2. – The model must work within any atmospheric conditions and possibly with air pollution of an urban environment. This requires to model radiative transfer of both the atmosphere above and the air among urban objects.

3. – An operational methodology must allow calibrating outputs of the remote sensing model in terms of 2D distribution of albedo and exitance (i.e., to produce image outputs). This calibration is important, because one cannot expect to have access to the optical properties of all urban surface elements, which vary in space (e.g., tiles of roofs have different reflectance values depending on their age) and time (e.g., wet and dry roofs will exhibit different anisotropy of their reflectance).

Figure 1.

DART calibration with a remote sensing image (Landsat‐8) for computation of the urban surface albedo over the city of Basel, Switzerland.

Here, we present a 3D radiative transfer model, DART, that fulfils these requirements and its recent improvements for studying urban and natural Earth landscapes with remote sensing acquisitions. We present also the approach that was recently designed and implemented to assess the spatial distribution of DART input parameters: optical properties of surface elements (e.g., roofs, streets, vegetation). Figure 1 summarizes this approach. It leads to DART simulated albedo and exitance maps that are calibrated with real‐time satellite acquisition. Ideally, these maps have a spatial resolution that is equal to that of satellite images that are used for the calibration.

2. The DART model

2.1. Presentation

DART computes radiation propagation in the three‐dimensional (3D) Earth/atmosphere system in the entire optical domain from the visible to the thermal infrared parts of the electromagnetic spectrum (EMS) [2–6]. The Earth surfaces and the atmosphere are simulated as a three‐dimensional (3D) medium (Figure 2). For any urban and natural landscapes, DART simulates the 3D radiative budget and acquisitions by satellite and airborne imaging radiometers and LIDAR scanners aboard of space and airborne platforms. The DART model, developed in the CESBIO Laboratory (www.cesbio.ups‐tlse.fr/fr/dart.htm) since 1992, can work with any 3D experimental landscape configuration (atmosphere, terrain geomorphology, forest stands, agricultural crops, angular solar illumination of any day, Earth/atmosphere curvature, etc.) and instrument specifications (spatial and spectral resolutions, sensor viewing directions, platform altitude, etc.).

Figure 2.

DART cell matrix of the Earth/atmosphere system. The atmosphere has three vertical levels: upper (i.e., just layers), mid (i.e., cells of any size), and lower atmosphere (i.e., same cell size as the land surface). Land surface elements are simulated as the juxtaposition of facets and turbid cells.

DART has been successfully employed in various scientific applications, including development of inversion techniques for airborne and satellite reflectance images [7–9], simulation of airborne sensor images of vegetation and urban landscapes [10], design of satellite sensors (e.g., NASA DESDynl, CNES Pleiades, CNES LIDAR mission project [11]), among others. DART forward simulations of vegetation reflectance were successfully verified by real measurements [12] and also cross‐compared against a number of independently designed 3D reflectance models (e.g., FLIGHT [13], Sprint [14], Raytran [15]) in the context of the RAdiation transfer Model Intercomparison (RAMI) experiment [16, 17].

DART creates and manages 3D landscapes independently from the RT modeling (e.g., visible and thermal infrared spectroradiometers, LIDAR, radiative budget). This multi‐sensor functionality allows users to simulate efficiently radiative transfer products of the same landscape as being captured by various sensors. Major scene elements are as follows: urban features, trees, grass and crop canopies, and water bodies. A DART simulated tree is made of a trunk, optionally with branches created with solid facets, and crown foliage that is simulated either as a set of facets or as a set of turbid cells, with specific vertical and horizontal distributions of leaf volume density. Trees of several species with different geometric and optical properties can be located within the simulated scene of any user‐defined size randomly or based on exact coordinates. Urban objects (houses, roads, etc.) contain solid walls and roofs built from triangular facets. Finally, water bodies (rivers, lakes, etc.) are simulated as facets of appropriate optical properties. DART can use external libraries to import and to some extent also edit landscape elements, digital elevation models (DEMs), and digital surface models (DSM) produced by other software or measured in field (e.g., translation, homothetic and rotational transformations). Most importantly, the imported and DART‐created landscape objects can be combined into virtual Earth scenes of user‐defined complexity. This allows importation of whole cities from urban databases provided by city councils and urban planners. The optical properties of landscape elements and their geometry, as well as and properties of atmosphere, are specified and stored in adjacent SQL databases.

Atmospheric cells are used to simulate attenuation effects of satellite at‐sensor radiance and also to model influence of the atmosphere composition on radiative budget of Earth surfaces. The atmosphere can be treated just as an interface above the simulated Earth scene or as a light‐propagating medium above and also within the simulated Earth scene, with cell sizes inversely proportional to the particle density. These cells are characterized by their gas and aerosol contents and spectral properties (i.e., phase functions, vertical profiles, extinction coefficients, spherical albedo). These quantities can be defined manually or obtained automatically from an atmospheric database. DART contains a database that stores the properties of major atmospheric gases and aerosol parameters for wavelengths between 0.3 and 50 µm. In addition, external databases can be imported, for instance, from the AErosol RObotic NETwork (AERONET; http://aeronet.gsfc.nasa.gov/) or from the European Centre for Medium‐Range Weather Forecasts (ECMWF; http://ecmwf.int/). Atmospheric RT modeling includes the Earth/atmosphere radiative coupling (i.e., radiation that is emitted and/or scattered by the Earth and backscattered by the atmosphere towards the Earth). It can be simulated for any spectral band within the optical domain from the ultraviolet up to the thermal infrared part of the electromagnetic spectrum. The Earth/atmosphere coupling was cross‐compared and successfully validated [18, 19] with simulations of the MODTRAN atmosphere RT model [20].

2.2. Recent improvements of DART

Set of improvements had to be recently implemented in the DART model in order to provide the optimal products required by the URBANFLUXES project.

2.2.1. Work preparation

The urban information used for the URBANFLUXES project was provided in the form of urban databases (for an example, see Figure 7). The two following adaptations had to be introduced to interface the databases with DART: (1) the format of 3D objects (i.e., triangular facets) stored in a common file format “*. obj”, and (2) the information about the vegetation, which was simulated as a turbid medium. Turbid vegetation was created from a set of characteristics, which for each tree in the urban scene provided geographical coordinates, physical dimensions, and optical properties. Those characteristics were obtained partially from external sources or partially from databases already exiting in DART (e.g., from database of optical properties for urban elements and vegetation).

2.2.2. Modeling of in situ camera acquisitions

In the frame of the URBANFLUXES project, in situ cameras are used for acquiring images of urban canopy. The major objective is to better assess the properties (i.e., spectral reflectance/emissivity and thermodynamic temperature) of the urban elements that are observed. The images acquired by these sensors can be very useful for validating the approach that we devised in order to retrieve maps of optical properties from satellite images. In this context, we introduced the modeling of these in situ sensors into DART. The major difficulty was linked to the fact that these sensors (e.g., fish‐eye cameras) have a wide field of view (FOV). Hence, different elements of a scene are not viewed along the same direction, which strongly impacts the radiometry and geometry of the acquired images. It is interesting to note that most remote sensing models neglect this fact: They consider that sensors have an infinitely small FOV. In short, DART can know simulate in situ cameras at any location in the Earth landscape (natural or urban), with any view direction, either upward or downward. Some examples of simulated images are shown here: downward looking sensor (Figure 3) that is on top of an urban district (Toulouse, France) and sensors with upward and oblique view directions (Figure 4).

Figure 3.

DART simulation of a fish‐eye camera acquisition above an urban district of Toulouse (France). Left: Red‐green‐blue spectral composite in natural colors; right: a thermal infrared image.

Figure 4.

DART simulation of in situ cameras: (a) Trees with leaves simulated as facets viewed by a camera with an upward looking direction, (b) and (c) trees with leaves simulated as turbid material in an urban environment viewed by a camera with a horizontal view direction.

Work continues in order to generalize the simulation of in situ sensors in case of any Earth surface element: fluids, atmosphere, water, etc., for any Earth scene configuration (i.e., isolated scenes and infinitely repetitive scenes with/without continuity of local digital elevation model). In addition to help in understanding and calibrating remote sensing measurements, the simulation of in situ sensors can be very helpful for a number of applications: to determine the optimal location and view direction of in situ sensors and to assess local atmosphere and pollution impact on sensor acquisitions.

2.2.3. Atmosphere database

In relation to the in situ cameras implementation, the DART atmosphere database and its management in the model were improved to offer higher flexibility of DART when dealing with different atmospheric conditions, especially in case of available in situ and/or satellite measurements of atmosphere. This atmosphere database was originally derived from simulations of the MODTRAN atmosphere radiative transfer model. Its use with DART was already validated with MODTRAN simulations. However, the simulation of aerosols was not as accurate as for gases. The DART atmosphere database was therefore completed using transmittance spectra for scattering and absorption mechanisms, derived from the atmospheric model MODTRAN. An interesting feature of this improvement is a new possibility to specify gas and aerosol amounts within the urban scene, independently of the atmosphere characteristics above the considered environment. This will be of a great help for assessment of local atmospheric properties and pollution impact using in situ sensor acquisitions. This improvement was accomplished by importing HITRAN [21] line‐by‐line cross section database (with specified temperature and pressure) for thermal infrared spectral domain, as well as the MPI‐Mainz [22] cross section database for visible and near‐infrared spectral domains.

Figures 5 and 6 show comparisons of DART and MODTRAN simulations in the visible and near‐infrared and the thermal infrared wavelengths, respectively. Both results demonstrate that the DART update and the introduction of new atmospheric database, combined with an improved radiative transfer modeling approach, brings DART simulations of atmosphere very close to MODTRAN 5.1 simulations, which is very encouraging especially for simulating accurately the in situ sensors.

Figure 5.

DART (red) vs. MODTRAN 5.1 (blue) simulations in short wavelengths (UV, VIS, near infrared). Gas model: US standard, aerosol model: Rural, visibility: 23 km. (a) Sun irradiance, (b) bottom of atmosphere (BOA) radiance, (c) top of atmosphere (TOA) reflectance (${\rho }_{ground}=0.5$).

Figure 6.

DART (red) vs. MODTRAN (blue) in the long wavelengths (thermal infrared). Gas model: Tropical, aerosol model: rural, visibility: 23 km. (a) Path radiance calculated at top of atmosphere (TOA) of scattered + emitted fluxes from atmosphere. (b) Direct transmitted radiance from Earth to TOA. (c) Total TOA radiance ($T_{ground}=299.15 \text{K}$). (d) TOA brightness temperature ($T_{ground}=299.15 \text{K}$).

2.2.4. Decomposition of a sensor image into images per type of scene element

A common difficulty for analyzing in situ sensor images (i.e., radiance images) is assessment of radiance and area proportion per type of surface material (e.g., wall, roof, atmosphere) inside an image pixel. Knowledge of the different radiance components is valuable for the “iterative calibration” that calibrates DART with remote sensing images as presented in Section 3.1.2. Hence, DART has been improved to facilitate in addition to the original sensor radiance image L_DART,Δλ(x_DART, y_DART, Ω_v) the per‐pixel simulations of radiance L_{DART,Δλ, n}(x_DART, y_DART, Ω_v) and cross section σ_n(x_DART, y_DART, Ω_v) images of each type of scene element n in discreet directions along the sensor viewing direction Ω_v).

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 $\text{Δ}\lambda$, 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 $A_{DART,\text{Δ}\lambda }$ and exitance $M_{DART,\text{Δ}\lambda }$ 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

Direct 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 (Switzerland) 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].

Band	Central wavelength [nm]	Lower wavelength [nm]	Upper wavelength [nm]	Spatial resolution [m]
1	443.0	435	451	30
2	482.6	453	512	30
3	561.3	533	590	30
4	654.6	636	673	30
5	864.6	636	673	30
6	1373.5	1363	1384	30
7	1609.1	1566	1651	30
8	2201.2	2107	2294	30

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 $\ge 2$ and with any spatial resolution.

3.1.1. Direct method

The direct method consists of the following five steps:

1. 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 ($x_{DART},y_{DART}$) of the DART voxels (in this example to 2.5 m).

   • The DART reflectance image of interest ${\rho }_{DART,\text{Δ}\lambda }\left(x_{DART},y_{DART}, {\text{Ω}}_{sat}\right)$ 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 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,\text{Δ}\lambda }\left(x_{DART},y_{DART}, {\text{Ω}}_S, E_{S,BOA}\left({\text{Ω}}_S\right),L_{atm}\left(\text{Ω}\right),{\text{Ω}}_{sat},t_{sat}\right)$E101

   where—$x_{DART},y_{DART}:$ coordinates of a given point in the DART simulated scene,

   1. — ${\text{Ω}}_S, {\text{Ω}}_{sat}$: sun and satellite view direction, respectively,

   2. — $E_{S,BOA}$: sun irradiance at the bottom of the atmosphere (BOA),

   3. — $L_{atm}(\text{Ω})$: atmosphere radiance,

   4. — $t_{sat}$: time of the satellite acquisition.

2. Step 2. Obtained DART image ρ_DART, _Δλ(x_DART, y_DART, Ω_S, E_{S, BOA}(Ω_S), L_atm(Ω), Ω_sat, t_sat) was georeferenced and resampled to the spatial resolution of the Landsat‐8 image ($x_{sat},y_{sat}$) to ensure exact correspondence in size and geolocation.

3. Step 3. A calibration factor $K_{\text{Δ}\lambda }\left(x_{sat},y_{sat},t_{sat}\right)$ is computed per DART pixel ($x_{sat},y_{sat}$). The objective of this factor is to mitigate the use of approximate optical properties. We use the following equation:

   $K_{\text{Δ}\lambda }\left(x_{sat},y_{sat},t_{sat}\right)=\frac{{\rho }_{sd,sat,\text{Δ}\lambda }(x_{sat},y_{sat}, {\text{Ω}}_S,{\text{Ω}}_{sat})}{{\rho }_{sd,DART,\text{Δ}\lambda }(x_{sat},y_{sat}, {\text{Ω}}_S,{\text{Ω}}_{sat},t_{sat})}$E1

4. Step 4. The map of urban spectral albedo $A_{\text{Δ}\lambda }$, is calculated per spectral band $\text{Δ}\lambda$ with:

   $A_{\text{Δ}\lambda }\left(x_{sat},y_{sat}, {\text{Ω}}_S, E_{S,BOA,\text{Δ}\lambda }\left({\text{Ω}}_S\right),L_{atm,\text{Δ}\lambda }\left(\text{Ω}\right),t_{sat}\right)=K_{\text{Δ}\lambda }\left(x_{sat},y_{sat},t_{sat}\right).A_{DART,\text{Δ}\lambda }\left(x_{sat},y_{sat}, {\text{Ω}}_S, E_{S,BOA,\text{Δ}\lambda }\left({\text{Ω}}_S\right),L_{atm,\text{Δ}\lambda }\left(\text{Ω}\right),t_{sat}\right)$E2

   where $A_{DART,\text{Δ}\lambda }$ is the albedo computed by DART, for spectral interval $\text{Δ}\lambda$:

   $A_{DART,\text{Δ}\lambda }\left(x_{DART},y_{DART}, {\text{Ω}}_S, E_{S,BOA}\left({\text{Ω}}_S\right),L_{atm}\left(\text{Ω}\right),t_{sat}\right)= \frac{{\rho }_{dh}{E}_{S,BOA,\text{Δ}\lambda }+{{\displaystyle \int }}^{\text{}}{\rho }_{dh,\text{Δ}\lambda }\left(\text{Ω}\right)L_{atm,\text{Δ}\lambda }\left(\text{Ω}\right)\cos \left(\theta \right)d\text{Ω}}{E_{S,BOA,\text{Δ}\lambda }+{{\displaystyle \int }}^{\text{}}{L}_{atm,\text{Δ}\lambda }\left(\text{Ω}\right)\cos \left(\theta \right)d\text{Ω}}$E3

   In (3), $\theta$ is the zenith angle corresponding to the viewing direction $\text{Ω}$.

5. Step 5. Desired urban surface albedo $A$ is calculated as the integral of all the spectral albedos $A_{\text{Δ}\lambda }\left(x_{sat},y_{sat}, {\text{Ω}}_S, E_{S,BOA,\text{Δ}\lambda }\left({\text{Ω}}_S\right),L_{atm,\text{Δ}\lambda }\left(\text{Ω}\right),t_{sat}\right)$ across the entire spectrum, weighted by the BOA irradiance $E_{S,BOA,\text{Δ}\lambda }\left({\text{Ω}}_S\right)$.

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 spatial resolution 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

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:

1. Step 1. Regular grid is laid over the satellite image, with a mesh size equal to M M satellite pixels and with $M>\sqrt{N}$. This ensures that each cell of the mesh contains a number of satellite pixels equal or larger than N.

2. Step 2. DART simulates the total radiance image $L_{DART,\text{Δ}\lambda }\left(x_{DART},y_{DART}, {\text{Ω}}_{sat}\right)$ of the scene and also the radiance images per scene element n $L_{DART,\text{Δ}\lambda ,n}\left(x_{DART},y_{DART}, {\text{Ω}}_{sat}\right)$ in the satellite viewing direction ${\text{Ω}}_{sat}$. Consequently, each surface element n viewed by a DART pixel d along the satellite viewing direction ${\text{Ω}}_{sat}$ has a reflectance value equal to ${\rho }_{n,d,\text{Δ}\lambda }^k(x_{DART},y_{DART})$.

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

   ${\rho }_{DART,\text{Δ}\lambda }\left(x_{DART},y_{DART}, {\text{Ω}}_{sat}\right)=\frac{\pi .L_{DART,\Delta \lambda }\left(x_{DART},y_{DART}, {\Omega }_{sat}\right)}{E_{DART,BOA, \Delta \lambda }}$E4

   The first iteration k = 1 assumes that the reflectance of an urban elements ${\rho }_{n,d,\text{Δ}\lambda }^k(x_{DART},y_{DART})$ 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 ${\rho }_{n,d,\text{Δ}\lambda }^{k+1}(x_{DART},y_{DART})$ 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 irradiance of an element type n is given by:

   $E_{DART,\text{Δ}\lambda ,n,d}\left(x_{DART},y_{DART}\right)=\frac{\pi .L_{DART,\text{Δ}\lambda ,n}(x_{DART},y_{DART})}{{\rho }_n^k(x_{DART},y_{DART})}.\frac{∆x_{DART}.∆y_{DART}.cos{\theta }_{sat}}{{\sigma }_{n,d}\left(x_{DART}, y_{DART},{\text{ Ω}}_{sat}\right)}$E5

   where ${\theta }_{sat}$ is the zenith angle of the satellite viewing direction and ${\sigma }_{n,d}\left(x_{DART}, y_{DART},{\text{ Ω}}_{sat}\right)$ is the cross section of the element of type n, which is viewed by the DART pixel d along the satellite viewing direction ${\text{Ω}}_{sat}$. It is also an output of the DART model. The ratio $\frac{∆x_{DART}.∆y_{DART}.cos{\theta }_{sat}}{{\sigma }_{n,d}\left(x_{DART}, y_{DART},{\text{ Ω}}_{sat}\right)}$ had to be introduced because the radiance of any DART pixel is simulated per effective square meter of the pixel area $∆x_{DART}.∆y_{DART}.cos{\theta }_{sat}$, while the radiance of the element n is expressed per effective square meter of the area of the surface element itself ${\sigma }_{n,d}\left(x_{DART}, y_{DART},{\text{ Ω}}_{sat}\right)$. Both radiances are equal, if only a single element is present in the considered pixel.

3. Step 3. DART radiance $L_{DART,\text{Δ}\lambda }({\text{Ω}}_{sat})$ and element $L_{DART,\text{Δ}\lambda ,n}({\text{Ω}}_{sat})$ images are resampled in the satellite image spatial resolution $(\text{Δ}{x}_{sat},\text{Δ}{y}_{sat})$. If selecting an appropriate spatial resolution of a DART simulation, one can make the assumption that each satellite pixel m is composed of an integer number $D^2$ of DART pixels. Therefore, for any given satellite pixel m, where d is the index of the DART pixels in m ($d\in [1 D^2]$) and n is the index of the element type ($n\in [1 N]$), one can define:

   • Radiance as:

     $L_{DART,\text{Δ}\lambda ,m}\left(x_{sat}, y_{sat},{\text{ Ω}}_{sat}\right)=\frac{{{\displaystyle \sum }}_{d=1}^{D^2}{L}_{DART,\text{Δ}\lambda ,d}}{D^2}$E6

   • Radiance of scene element n as:

     $L_{DART,\text{Δ}\lambda ,n,m}\left(x_{sat}, y_{sat,{\text{ Ω}}_{sat}}\right)=\frac{{{\displaystyle \sum }}_{d=1}^{D^2}{L}_{DART,\text{Δ}\lambda ,n,d}}{D^2}$E7

   • Irradiance of elements of type n as:

     $E_{DART,\text{Δ}\lambda ,n,d}\left(x_{sat},y_{sat}\right)=\frac{{{\displaystyle \sum }}_{d=1}^{D^2}{E}_{DART,\text{Δ}\lambda ,n,d}\left(x_{DART},y_{DART}\right).{\sigma }_{n,d}\left(x_{DART}, y_{DART}\right)}{{{\displaystyle \sum }}_{d=1}^{D^2}{\sigma }_{n,d}\left(x_{DART}, y_{DART}\right)}$E8
     and

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

     ${\sigma }_{n,m}\left(x_{sat}, y_{sat}\right)= {\displaystyle \sum }_{d=1}^{D^2}{\sigma }_{n,d}\left(x_{DART}, y_{DART}\right)$E9

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

     $E_{DART,\text{Δ}\lambda ,n,m}\left(x_{sat}, y_{sat}\right)=\frac{\pi .L_{DART,\text{Δ}\lambda ,n}(x_{sat}, y_{sat},{\text{ Ω}}_{sat})}{{\rho }_{n,m}^k(x_{sat}, y_{sat})}.\frac{∆x_{sat}.∆y_{sat}.cos{\theta }_{sat}}{{\sigma }_{n,m}\left(x_{sat}, y_{sat},{\text{ Ω}}_{sat}\right)}$E10

     In this expression, ${\rho }_{n,m}^k\left(x_{sat}, y_{sat}\right)$ 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.

4. 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 requirements too large, the computation enters in a next step.

5. Step 5. The actual calibration of DART for each cell of the mesh defined in step 1 is conducted for all $M^2$ satellite pixels and in each cell u of the mesh. Since $M^2>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^2$ satellite pixels m, leading to a system of $M^2$ equations verifying if the DART image and the satellite image are equal:

   ${\displaystyle \sum }_{n=1}^N L_{DART,\text{Δ}\lambda ,n,m}(x_{sat}, y_{sat},{\text{Ω}}_{sat})=L_{sat,\text{Δ}\lambda ,m}\left(x_{sat}, y_{sat},{\text{Ω}}_{sat}\right), \forall m\in [1 M^2]$E11

   Obviously, the verification is negative if the two images differ. At iteration k, the DART radiance values $L_{DART,\text{Δ}\lambda ,n,m}(x_{sat}, y_{sat},{\text{Ω}}_{sat})$ are computed with inaccurate approximated optical properties ${\rho }_{n,u}^k\left(x_{sat}, y_{sat}\right)$. If we define $L_{DART,\text{Δ}\lambda ,n,m}^{\text{'}}(x_{sat}, y_{sat},{\text{Ω}}_{sat})$ as the DART radiance value computed from ${\rho }_{n,u}^{k+1}\left(x_{sat}, y_{sat}\right)$, and if we accept the fact that radiance values are proportional to reflectance values [Eq. (10)], then we can write:

   $\frac{L_{DART,\text{Δ}\lambda ,n,m}\left(x_{sat}, y_{sat},{\text{Ω}}_{sat}\right)}{{\rho }_{n,u}^k\left(x_{sat}, y_{sat}\right)}=\frac{L_{DART,\text{Δ}\lambda ,n,m}^{\text{'}}\left(x_{sat}, y_{sat},{\text{Ω}}_{sat}\right)}{{\rho }_{n,u}^{k+1}\left(x_{sat}, y_{sat}\right)}$E12

   Consequently, we can rewrite the system of Eq. (11) to a system of $M^2$ equations, in which the unknowns are the expected reflectance values ${\rho }_{n,u}^{k+1}\left(x_{sat}, y_{sat}\right)$:

   ${\displaystyle \sum }_{n=1}^N\frac{{\rho }_{n,u}^{k+1}\left(x_{sat}, y_{sat}\right)}{{\rho }_{n,u}^k\left(x_{sat}, y_{sat}\right)} .L_{DART,\text{Δ}\lambda ,n,m}(x_{sat}, y_{sat})=L_{sat,\text{Δ}\lambda ,m}(x_{sat}, y_{sat}), \forall m\in [1 M^2]$E13

   Resolving this system for each cell u of the mesh leads to the new desired reflectance values ${\rho }_{n,u}^{k+1}\left(x_{sat}, y_{sat}\right)$, 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 resolution. 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

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 $E_{S,BOA,\text{Δ}\lambda }$ and diffuse sky irradiance $E_{atm,\text{Δ}\lambda }$. This information can be obtained from several sources including ECMWF (www.ecmwf.int), Meteo France (www.meteofrance.com), and other in situ meteorological sensor networks.

In this alternative approach, DART calculates the so‐called white sky and black sky albedos. The white sky albedo $A_{DART,white sky, \text{Δ}\lambda }(x_{sat},y_{sat},E_{S,BOA}\left({\text{Ω}}_S\right)=0,L_{atm}=\frac{E_{atm}}{\pi },t_{sat})$ is computed for an irradiance coming solely from the atmosphere, that is, without any direct sunlight contribution. The black sky albedo A_{DART, white sky,  Δλ}(x_sat, y_sat, E_{S, BOA}(Ω_S)=0, L_atm=0, t_sat) is computed for the direct sunlight only. The desired final albedo product is then computed as a combination of both components:

$A_{\text{Δ}\lambda }\left(x_{sat},y_{sat}, {\text{Ω}}_S, E_{s,BOA}\left({\text{Ω}}_S\right),E_{atm},t_{sat}\right)=E_{s,BOA}\left({\text{Ω}}_S\right). A_{\text{DART},\text{black sky},\text{Δ}\lambda }+E_{atm}. A_{\text{DART},\text{white sky},\text{Δ}\lambda }$E14

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 spatially 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 $ϵ\times L\left(T\right)$, where $L\left(T\right)$ is Planck's law for temperature $\left(T\right)$. 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 for these 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.

Thermal; band [nm]	Central; wavelength [nm]	Lower; wavelength [nm]	Upper; wavelength [nm]	Spatial; resolution [m]
1	10,895	10,600	11,190	100
2	12,005	11,500	12,510	100

Table 2.

Landsat‐8 thermal infrared band specifications.

4. Conclusions and future perspectives

An innovative methodology using a 3D RT modeling and satellite reflectance data to derive urban surface albedo images has been designed and implemented. The method was tested using an atmospherically corrected Landsat‐8 multispectral image of the city of Basel. Although the preliminary results are encouraging, the newly presented methods have to be still properly validated. Different validation approaches are foreseen for testing their robustness:

• Application of the methodology on other satellite images acquired over the city of Basel. Using images of other satellite platforms with a different spatial resolution within relatively short time intervals (e.g., Landsat‐8 images acquired earlier or later and/or images from Sentinel‐2) will allow us to check how much the derived optical properties vary in time (urban surfaces elements’ optical properties are expected not to vary significantly in time, whereas optical properties of vegetation should be more variable).

• Application of the methodology over additional cities, in particular the cities of Heraklion and London. These two other cities are indeed of scientific interest in the URBANFLUXES project. This will allow us to test the robustness of the method for different climatic conditions and urban structures, and to study the impact of those differences in terms of method accuracy as well as technical constraints.

A comprehensive comparison of the direct and iterative method technical performances is also required. It is especially important for future transformation of the methodology into a standardized operational approach.

Another important development lies in finalization of the thermal infrared spectral domain calibration. Combination of both calibrations in the visible and thermal part of the electromagnetic spectrum would lead to a complete 3D radiative budget produced by DART at the same time as the other RT products.

To summarize, an original methodology for deriving urban surface albedo images from remote sensing data without a need of in situ measurements was implemented, although ground measurements could improve the final results. Both calibrations are computationally operational, but the results for the iterative method are considered as preliminary, due to the limited implementation of only a single iteration. Iterative method is expected to be more accurate than for the simpler direct one. Nevertheless, the direct calibration might prove to be easier to use and more robust in case of inaccurate urban datasets, which depends on spatial resolution and geometric co‐registration of DART and satellite imagery. Finally, both methods rely on the fact that DART is a reliable RT model simulating both satellite images and the radiative budget of cities with an acceptable accuracy. In a perspective of increasing remote sensing data availability, arising from new satellite systems such as Sentinel‐2, these approaches are expected to play an important role in future surveys of cities and their dynamic climate.