Spectra: comparison between plan-parallel and spherical
In this example, we present how you can calculate a reflectance spectrum with htrdr-planets from a set of 1D data. Here, we compare between plan-parallel and spherical geometries.
Imports
We start by importing the necessary libraries.:
import numpy as np
import matplotlib.pyplot as plt
import htrdrPy as htrdr
import scipy.constants as cst
import json
Data
The first step is to generate an Data object that will contain the
optical properties. This object aims at handling the data in order to generate
the input files. In this example, the optical properties are loaded from an
.npz file that you can find here. We first create an instance of the
Data object with the correct planet radius:
radius = 2.575e6 # m
dataSp= htrdr.Data(radius = radius, nTheta=80, nPhi=80, name="Sphere")
dataPP = htrdr.Data(radius = 1e9, name="PP")
Note
In plan parallel, the radius is used to define the x and y extension of the atmosphere. Make sure this value is large enough, but also note that a too large value can result in issues within htrdr. For this example, we found that any value larger the 1e9 would result in inconsistent results.
Note
At the creation of each Data instance, an inputs_{name}/ repository is created and will contain all the input files necessary for htrdr. A outputs_{name}/ repository is also created and will contain the files created by htrdr.
Then, we read the data and convert it into a python disctionnary. The
Data object will work with the dictionnary, which therefore needs
to have the correct keys. In this example, the keys used to generate the
.npz file are already the right ones. Each corresponds to a table with a
given shape. All these informations (keys name and array shapes) are provided
further down in this documentation (c.f. htrdrPy.data module).
data = np.load("Example_spectra_data/1D_spetral_data.npz")
data = dict(data)
for key, array in data.items():
print(key, ", shape : ", array.shape)
Additional information, such as the dimensions (number of angles in the phase function, number of wavelength, etc.) are also required to be present in the dictionnary. We therefore add the necessary information:
nWavelength, nLevel, nCoeff = data["absorption (m-1)"].shape
nAngle = len(data["angles (°)"])
wavelengths = data["wavelength"]
data.update({
'nLevel' : nLevel,
'nWavelength' : nWavelength,
'nCoeff' : nCoeff,
'nAngle' : nAngle,
})
print(data.keys())
The next step is to provide those data to the Data object.
Depending on the context, different methods exist to pass the information. In
this example, we use the Data.makeMixture() method, which handles the
case where the data provided are those of the atmosphere mixture comprising all
the gases and aerosols. The dim parameter indicates the dimension of the
provided data: 1 means that we provide a column, 2 for a slice along altitude
latitude and 3 for a fully 3D heteorgeneous atmosphere. Obviously, the shapes of
the tables provided through the dictionnary change accordingly. dim=0 sets a
plan-parallel atmosphere.
dataSp.makeMixture(data, dim=1)
dataPP.makeMixture(data, dim=0)
Then, we generate the surface by providing the temperature map (a single float
in 1D) and the brdf dictionnary to the Data instance through
the Data.makeGround() methods. The method to use depends on the
dimension of the input data.
brdf = {
"kind": 'lambertian',
"albedo": data['surface albedo'],
"wavelengths": wavelengths
}
dataPP.makeGroundFrom1D_PP(91, brdf)
dataSp.makeGroundFrom1D(91, brdf)
Generating input files
We can now generate the input files with the command
Data.writeInputs():
dataPP.writeInputs()
dataSp.writeInputs()
We can also generate the VTK files with the command
Data.writeVTKfiles(). VTK files are not necessary and only serves to
visualize the data that will be send to htrdr-planets.
dataPP.writeVTKfiles()
dataSp.writeVTKfiles()
Observation geometry
Once the input files are generated (this can actually be done before), we have
to define the observation geometry. These informations concerns the camera (such
as its position or the target point observed), the image (the (x,y) definition
and the number of sample per pixel) and the source (size, radius, distance,
longitude, latitude, etc.). We start by creating an instance of
Geometry (two actually, one for the plan-parallel
calculation and one for the spherical calulation):
geomPP = htrdr.Geometry(case='PP')
geomSp = htrdr.Geometry(case='Sphere')
Note
A geometry can be used multiple times, with different
Data and Script object. Here, we have
to create one insatnce for each because the target won’t be at the same
position in both case: in the spherical scenario, the planet radius adds a
shift to the scene.
In the current scenario, we aim to reproduce some observation with constraints
being the incidence, emergence, phase and azimut angles. To calculate the camera
and source position, the Geometry objects has a method
Geometry.makeGeomFromAPIE() (a routine that makes you happy!) that
will automatically generate the geometry from the constraints on the observation
angles. Those information are read from a file (contained in the data repository
downloaded earlier). Additional information are requirted by the routine, such
as the distance between the camera and target point and the source properties
(distance, size and temperature).
with open("Example_spectra_data/observation.json", 'r') as f:
obs = json.loads(f.read())
cameraDist = 1000 * cst.kilo # m
solDist = 10 * cst.au # m
solRad = 7e8 # m
solTemp = 5800 # K
geomPP.setImage([64,64], 10)
geomPP.makeGeomFromAPIE(obs, 0, cameraDist, solDist, solRad, srcTemp=solTemp)
geomPP.exportGeometry()
geomSp.setImage([64,64], 10)
geomSp.makeGeomFromAPIE(obs, radius, cameraDist, solDist, solRad, srcTemp=solTemp)
geomSp.exportGeometry()
Note
The source temperature will be used by htrdr to determine the spectral
distribution of the incoming stellar radiation, but note that is it also
pausible to directly use a spectrum (c.f.
Geometry.setSource() documentation).
Note
The image information must be provided separately through
Geometry.setImage() or when creating the
Geometry instance.
Note
Geometry.exportGeometry() will produce a file within the
geometries repository, containing the source, camera and image
parameters.
Scripts
Now, we create an instance of Script that handles the last details
of the calulation and the call to htrdr. The different methods of the
Script object simplifies the call to htrdr via a bunch of
predefined scirpts (c.f. htrdrPy.script module for a full review of the
possibilities). Here, we start a reflectance spectrum (I/F) and we therefore use
the Script.reflectanceSpectrum() method, which requires the instance
of Geometry previously created, the type of calculation (“sw” or
“lw”) and the list of wavlengths:
scriptPP = htrdr.Script(case='Spectrum',
MPIcmd="mpirun -np 4 --map-by socket:PE=9")
scriptPP.reflectanceSpectrum(geomPP, "sw", wavelengths)
scriptSp = htrdr.Script(case='Spectrum',
MPIcmd="mpirun -np 4 --map-by socket:PE=9")
scriptSp.reflectanceSpectrum(geomSp, "sw", wavelengths)
Then, we can start the calculation by calling the instance of
Script with the previously created instance of Data:
scriptPP(dataPP)
scriptSp(dataSp)
Post-process
Finally, a call to Postprocess will treat the raw output of
htrdr to generate the required results, here it will be spectra, stored in
the results_Sphere/ and the results_PP/ repositories:
htrdr.Postprocess(scriptPP)
htrdr.Postprocess(scriptSp)
We can now recover and plot the calculated spectra along with the observed specrtum as well as additonal data, all provided in the spectrum.txt file previously downloaded:
with open("results_PP/reflectance_spectrum_Spectrum.json", 'r') as f:
resPP = json.loads(f.read())
with open("results_Sphere/reflectance_spectrum_Spectrum.json", 'r') as f:
resSp= json.loads(f.read())
with open("Example_spectra_data/spectrum.txt", 'r') as f:
f.readline()
wvl = []
IF_obs = []
IF_SHDOM = []
IF_MCC = []
for line in f:
l = line.split()
wvl .append(l[0])
IF_obs .append(l[2])
IF_SHDOM.append(l[3])
IF_MCC .append(l[4])
wvl = np.array(wvl , dtype=float)
IF_obs = np.array(IF_obs , dtype=float)
IF_SHDOM = np.array(IF_SHDOM, dtype=float)
IF_MCC = np.array(IF_MCC , dtype=float)
print(wvl)
print(IF_obs)
fig, ax = plt.subplots()
ax.plot(wvl, IF_obs , label="Observed", ls="", marker="o")
ax.plot(wvl, IF_SHDOM, label="SHDOMPP", ls="", marker="o")
ax.plot(wvl, IF_MCC , label="MCCSPHE", ls="", marker="o")
ax.errorbar(np.array(resPP['wavelength'])/cst.micro,
np.array(resPP['reflectance spectrum']),
yerr = 3 * np.array(resPP['reflectance std deviation']),
label = "plan-parallel")
ax.errorbar(np.array(resSp['wavelength'])/cst.micro,
np.array(resSp['reflectance spectrum']),
yerr = 3 * np.array(resSp['reflectance std deviation']),
label = "spherical")
ax.set_ylabel(r'$\frac{I}{F}$', rotation='horizontal', fontsize=15)
ax.set_xlabel('Wavlength [µm]', fontsize=15)
ax.legend()
fig.savefig('Spectrum.png')
