Preparing an execution

In order to execute the model, some data manipulation is needed to prepare it for the model. It is strongly recommended to separate the data from the code by creating a new directory somewhere on your computer and placing all the relevant data within.

When it is done, executing the init_run.py script will copy the necessary files to the working directory. The parameter files can then be modified to contain the correct values for your experiment.

Downloading and preparing the required satellite images

ILLUMINA requires some satellite data to run properly, namely a digital elevation model and the nocturnal light emittance. Theses data also need to be projected in a suitable spatial reference system and clipped to the simulation domain.

Domain definition

Defining the simulation domain is crucial in the input preparation, as it directly affects everything afterwards. The first step is to define the location(s) where the simulation of the artificial sky radiance is desired. Then, the projection needs to be defined, as the model need to work with coordinates in meters instead of degrees.

The parameters file named domain_params.in defines the domain. It should contain the following parameters:

• latitude and longitude are the coordinates of the observer. These two parameters can also be lists (of the format [lat1, lat2, ...]) for multiple observing locations.
• srs is the spatial reference system used by the model. Setting it to auto' lets the model select an appropriate one automatically, but it can also be selected manually by writing the corresponding EPSG code (for instance epsg:3561'').
• The domain is defined by a series of overlapping layers with different scale. This allows to describe vast domains while keeping a high resolution close to the observer. The last 5 parameters describe theses layers. ⁽¹⁾
  • nb_layers is the number of these layers,
  • nb_pixel is width of the domain in number of pixel. It will be the same in each layer ⁽²⁾,
  • scale_min is the resolution in meters of the smallest layer, We strongly suggest not to use a value lower than 20 m. This is a hard coded limitation in the illumina kernel. If you select a smaller value, the model will make the computation but we cannot confirm that the results will be good.
  • scale_factor is the ratio between the dimension of pixel of successive layers and
  • buffer is the size of the buffer (in kilometers) around each layer to properly consider each possible optical path the light can use to reach the observer. This is mainly useful to allow the calculation of the interaction of photons reflecting on the ground or 2^nd order of scattering photons blocking by the elevation model. Since in the model, the 2^nd order of scattering is calculated up to 100 km away from the source and from the line of sight, a buffer of at least 100 km wide would be ideal. However for most resolutions such a huge buffer isn't realistic. The native size of the domain for any modelling scale is 512 x 512 pixels. If the buffer exceed this size (i.e. nb_pixel+2xbuffer>512), the maximal value of the buffer will correspond to the domain size, not the requested buffer.

The defineDomain.py script is then used to generate the domain.ini file containing the details of the defined layers. It will print the geometric properties of each layer so you can validate that the dimensions are reasonable. We suggest a largest domain size of around 300-600 km. The script also prints the coordinates of the south-west and the north-east corners. Theses are useful for bounding the domain to download only the relevant satellite imagery in the following steps.

⁽¹⁾ Note that defineDomain.py can be called as often as required until you are satisfied with the layers/domain properties.

⁽²⁾ As a rule of thumb, we suggest not to exceed 255 for that number.

VIIRS-DNB imagery

The night emittance is obtained from VIIRS-DNB imagery that can be found here. One should download the appropriate tile(s) for the period desired (year and month). We suggest the VCMSLCFG configuration in the monthly composite product because of the stray light correction, but the choice is left to the user. You will want to extract the 'avg_rade9.tif' file as it's the actual values, whereas the 'cf_cvg.tif' file contain information related to the quality of the image. The tif file(s) should be placed inside a subfolder named VIIRS-DNB inside your experiment directory.

It is also possible to use the VIIRS background method proposed by Coesfeld et al. (2020) for more accurate results. In that case, the VCMCFG product needs to be used instead, and the correction data needs to be downloaded from here and decompressed in the VIIRS-DNB folder.

Watermask

When used with VIIRS-DNB input, the model need a water mask to calculate properly the light fluxes. You can download it here and save it to the experiment folder.

SRTM data

The digital elevation model is made with the SRTM elevation data that can be found here. One should use the spatial filter to select only the required tiles, and then follow the download procedure. The extracted hgt files should be placed inside a subfolder named SRTM inside your experiment directory.

Processing the input images with Illuminutils.py

The Illuminutils.py script should be executed from the experiment directory containing the two data subfolders explained above.

Two files should be produced by this script:

• stable_lights.hdf5
• srtm.hdf5

VIIRS	SRTM
stable_lights.bin	srtm.bin

Sample files for Hawaii

The standard format used by ILLUMINA is HDF version 5. Theses can be visualized using tools like hdfcompass. We also provide convenience python functions in the hdftools package included with ILLUMINA.

Converting other datasets

Illuminutils.py can also be used to convert other datasets to the Illumina format for a specified simulation domain. As long as the domain.ini file is in the current working directory, the script can be called as

where NAME will be the name of the output file (without the extension) and file list is a list of one or multiple files (the use of bash wildcards is recommended here) to be warped. Note that all the files will be warped togheter so they should be tiles of the same dataset. The supported formats are the ones that can be processed by GDAL.

Making light inventory

In order to model the propagation of the light, the properties of the light sources must be defined. There are two way to do this for ILLUMNA: 1- using VIIRS-DNB spaceborne radiance monthly product, 2- using a point source inventory. Both of them can be used together, as long as they are not overlapping. There can not be poing sources where in a pixel already containing sources derived with VIIRS-DNB.

Using uniform overlapping circular zones

The first way is to define overlapping circular zones of uniform properties. Theses zones are defined by a their center point and a radius and specify the mix of lamps assumed to be present in that area (different by their photometry function or light output pattern (LOP), their spectrum, lamp height) as well as the average distance between obstacles, obstacles height and obstacle filling factor. Two or more zones may be in the same geographical region or partly overlapping. Each new zone overwrite the previous in case of intersection between the zones. All the points that are not included in a zone will be ignored. To create a zone you have to edit an ASCII file with a simple text editor like kwrite or gedit following the format shown below:

Sample inventory file for the Hawaii territory

This file can have any number of header lines as long as they are preceded by a '#' symbol. Anything on the same line following this symbol will not be considered by the model. Each data line contains several parameters:

1. lat: central latitude of the circular zone.
2. lon: central longitude of the circular zone.
3. R: radius of the circular zone (in kilometers).
4. hobs: subgrid averaged obstacles height (in meters)
5. dobs: averaged distance between subgrid obstacles (in meters)
6. fobs: obstacle filling factor i.e. probability for a photon to hit an obstacle (0. to 1.)
7. hlamp: averaged lamps height relative to the ground (in meters)
8. list of lamps characteristics

Each lamp characteristics is composed of three fields separated by '_'.

• The first field is the weight of the zone that is defined by the following two characteristics. The weight is later converted to a ratio by normalizing to the sum of all the weight for that zone.
• The second is a reference word ⁽³⁾ corresponding to the spectral power distribution of the lamp.
• The third is a reference word ⁽³⁾ corresponding to the angular power distribution of the lamp or light ouput pattern (LOP)⁽³⁾.

The weighting is applied on the luminous flux of the spectral power distribution of the lamp in lumen. This means that the spectrums are weighted by the photopic sensitivity curve.

The data referenced by the last two fields must be located in a subfolder named 'Lights'. This folder must contain the following files in addition to the ones used to define the lamp inventory :

• photopic.dat
• scotopic.dat
• viirs.dat

Theses files can be found in the ILLUMINA installation directory (Examples/Lights). Any additional file used to characterize the lamp must follow the following format :

• Angular light output pattern (LOP) files must have the extension '.lop'. They are made of two columns ASCII data where the first column is the relative intensity and the second is the zenital angle in degree. The lop file must contain 181 data starting at z=0 to z=180 at 1 deg. step.
• spectral files must have the extension '.spct'. They are two columns ASCII data files with a single line header. The first column contains the wavelength in nm and the second contains the relative intensity.

The normalization of all theses files is not important, as it will be done by the programs.

⁽³⁾In all cases, any characters preceding the first underscore (_) in the lop or spct file name is the reference word that must be written in the inventory file.

All .spct files must have the same wavelength scale. All LOP Files must share the same angle scale.

Using a discrete light source inventory

The second way to describe the lights is to directly specify their properties on a lamp-by-lamb basis. In this case, the file needs to have the following format:

where

1. lat: Latitude of the light source
2. lon: Longitude of the light source
3. pow: Intensity of the source in lumen
4. hobs: Averaged obstacles height (in meters)
5. dobs: Averaged distance between obstacles (in meters)
6. fobs: Obstacle filling factor i.e. probability for a photon to hit an obstacle (0. to 1.)
7. hlamp: Light source height relative to the ground (in meters)
8. spct: Spectral power distribution keyword
9. lop: Angular power distribution keyword

Defining the experiment

The execution mode

You may be interested to run ILLUMINA for many reasons. By default, ILLUMINA will calculate the artificial diffuse radiance, the part that is produced by the clouds, the direct radiance reaching the observer from a sight to the sources and the direct radiance coming from a sight to reflecting surfaces. If you are more interested by the direct radiances, it may be a good idea to increase to the maximum the resolution near the observer. The calculation of the direct radiance inside the mean free path toward obstacles will not experience any obstacle blocking. The blocking by obstacles only occur when the observer is farther than the mean free path to the ground. This parameter is defined when you specified the "subgrid obstacle properties" with the variable dobs. Actually dobs is twice the value of the mean free path. If you are not interested to obtain the sky or cloud radiances, but only the direct radiance, then you can speedup the calculation by setting off the scattering.

Create the input parameters file

The parameters used by the model for executing the experiment are contained in the input_params.in file, as described below:

The reflective surface types are ASTER files located in the `Lights` folder.

Run make_inputs script

Once all the data is obtained as the input parameter file is created, the make_inputs.py script is used to prepare the data for the model.

Quick file check

The script should produce a directory named 'Inputs'. It should contain:

• N * W fctem_wl_WWW_zon_NNN.dat files, where N is the number of zones and W the number of wavelength used. In the case of the example, both these numbers are 5, and so you should have 25 fctem files. This is the angular emission informations.
• You should also have lumlp files, one for each wavelength and zone combination plus one for each wavelength, giving the global view. lumlp files are giving the total lamp spectral flux for any pixel of a zone at the given wavelength. Note that this flux is not corrected for the atmospheric attenuation and obstacles blocking between sources and the satellite. Such a correction will be done later on while running the model.
• There also must be a file named exp_name_altlp.hdf5 . This is the lamp height relative to the ground.
• A file named exp_name_obstd.hdf5 . This is the distance between obstacles.
• A file named exp_name_obsth.hdf5 . This is the height of obstacles.
• A file named exp_name_obstf.hdf5 . This is the obstacle filling factor.
• A file named origin.hdf5 . This is the a flag array telling if the inventory was determine with VIIRS or from point sources.
• You should also find two .lst files.
• In addition to all that, you should see some symbolic links, one for integration_limits.dat, one for srtm.hdf5 and W .mie.out files, one for each wavelength W.

Example run at 605nm with contrast boosted

Total lumlp	Zone 1 = Oahu
Zone 2 = Molokai and Lanai	Zone 3 = Maui
Zone 4 = Big Island	Map of the zones made with free map tools

Zone 5, wich is the lava lakes, is all black because its light isn't considered in the model. To achieve this, its lamp inventory was empty. You can verify that in the 'inventory.txt' file.

Alternative scenarios

You may be interested in simulating alternative scenarios based on the current situation. For example, artificially replacing all light sources to a new photometry. This is done with the alt_scenario_maker.py script. Help on that script is available by calling

If used with an alternative zones inventory, a replacement inventory needs to be in the same directory and contain a set of lamp characteristics for each zone. For example,

If used with an alternative lamps inventory, a replacement inventory needs to be in the same directory and contain a set of lamp characteristics in the same format as the initial inventory.

The script will then generate a folder named Inputs_NAME containing the relevant data.

Submitting the calculations to a Linux cluster

To perform the calculations, we now connect to a 'Cluster'. In our case, we connected to 'Mammouth serial II' located at Université de Sherbrooke. The task scheduler used is Slurm. You may need to manually adjust some files to match the execution environment you are using.

Then it is necessary to recompile the ILLUMINA model using 'makeILLUMINA' or through the following command:

The `Inputs` directory created for each experiment in step 6 should now be transferred to the cluster interactive node via the scp protocol.

Preparing the batch execution

Now we need to run the program makeBATCH.py that will prepare the execution folder for each calculation directory on the cluster. This must be done from the Inputs folder(s). The documentation of the function is available by calling makeBATCH.py -h. Note that files with names conflicting with the batch name provided either at the command line or in the 'input_params.in' file will be removed prior to executing. If you want to prepare multiple execution, make sure that they have different batch names.

To execute the calculations, simply execute the bash file(s) produced by the makeBATCH.py script.

Find failed calculations

In many cases you will probably have a lot of calculations to be done to complete your modeling experiment. Each calculation going to a given core and/or node (if you run on a cluster). Then for some reasons there is some chance that some of your calculations can fail. Finding the failed calculations can be a difficult task. For that reason we provided a script called find-failed-runs.py. All you need is to wait for all calculations to finish and then go to the experiment folder and run the script. The script will show the path of folders containing failed runs. If you run it with the -e option, the scritp will generate the code to launch the failed runs. You should probably want to store it in a file and then run it as a bash script.

Then simply start the aborted runs by running this script

Extracting results

ILLUMINA generates two different output per calculation:

1. An image file showing the relative contribution of each pixel to the calculated diffuse radiance. We call this image the PCL file or the contribution map.
2. The numerical values for the Total diffuse radiance, the cloud radiance, the direct radiance from sources and the direct radiance from reflection surfaces.

To extract the data, the extract-output-data.py script is used. It can extract either the value of the diffuse radiance or all available components. It can also extract the contribution maps. Moreover, filters are available to only process certain parameter values. Documentation is available by calling

The script will output the data directly, so it should be redirected to a file with

If you are not only interested in the total diffuse radiance (clouds + preceeding atmosphere) and want also to extract the cloud contribution to the radiance and the direct and direct reflected radiance, you will neet to run the script in the full mode.

The script will output the data directly, so it should be redirected to a file with

There will be a column for each radiometric value.

As stated in the documentation, the contribution map can be extracted using the `-c` flag.

Analysing the results

The analysis can be done with your favorite tools. We strongly recommend the use of python and provide convenience functions in the pytools and hdftools packages provided with ILLUMINA.

Transforming to magnitude per arc sq seconds (for astronomers...)

Transforming diffuse radiances to sky brightnesses (SB) in units of mag/sqarcsec is not an simple task. First of all you have to consider that illumina is only dealing with the artificial component of the SB. If you are using illumina in a relatively dark site, the artificial SB can represents only a small part of the total SB. To transform radiance to total SB, you will need a relevant estimate of the natural contribution to the total SB. The natural SB is highly variable with time, altitude, season, observing direction etc. It is composed of many sources like the zodiacal light, the starlight, the sky glow, the Milky Way etc. Given that complexity, we suggest to determine it experimentally for the modeled site and period you are interested in. To do it, you need an in situ measurement of the total SB from which you will be able to extract the natural component and then eventually consider this component as a constant natural contribution to the SB for your specific site and period, no matter the viewing angle or light inventory, obstacles properties etc. Lets call the radiance responsible for that natural contribution the background radiance (R_{bg}). Lets assume that you have an in situ measurement of the total Johnson-Cousins SB. You need to accomplish the following steps to convert your artificial modeled radiance to total SB.

Integrate your radiances {R_a} and according to the Johnson-Cousins filter. {R_a} being the modeled radiance you want to convert to SB . The sensitivity curve of Johnson-Cousins filters are provided in the Example/Lights folder (e.g. JC_V.dat).

R_{bg} the radiance corresponding to the natural level of the sky brightness (SB_{bg}) can be estimated using the SB data provided by Benn & Elison 1998, except for the R band that was taken from La Palma P99 ASTMON (2018-2019) measurements. We simply added 0.03 from the B and V values provided for La Palma (as recommended by Benn & Elison 1998) and then used the formula below to calculate the radiance.

For a value of modeled artificial radiance ({R_a}), use the following formulae to convert to total SB:

R_0 are derived from zero points obtained by Bessell 1979 calibration (DOI 10.1086/130542) and given in the table below.

Quick git using guide for developers

To publish your locally modified codes on the server

To update your local version of the codes