Created attachment 215 [details]
Radiance spectra from DIRSIG and PLEXUS
What are the DIRSIG image units? I assumed that they were Watts/cm^2/steridian/micron, but as the attached graph shows, the values from a radiance spectrum produced by DIRSIG are far below those produced by PLEXUS. Both of these graphs were produced under the same conditions: 1:00 p.m. on 1 June 2001 at 43N 77W.
What is the format of that attachment? It helps if screen grabs are things that are portable like JPEG, PNG, etc.
Your units are correct, maybe. It really depends on the type of response you are using, the spectral units, etc. If you are using a functional channel shape (Gaussian, Rectangular, etc.) defined with "microns", then your units are correct
as W/(cm^2 sr micron). But, if you used a tabulated channel (wavelength and response pairs), then it is not "spectral radiance" (not normalized by the area under the curve) and hence it would be W/(cm^2 sr).
Created attachment 216 [details]
.jpg version of the Plexus-Dirsig comparison
Sorry, the attachment was an Enhanced Meta File (.emf), which I guess is used by Microsoft rather than Apple. Here is the .jpg version.
Yes, I am using the AVIRIS detector array response: 224 Gaussian channels between .36 microns and 2.55 microns. (The channel spectral response graph have units of microns on the x-axis.)
The problem is, I don't know which is correct, but somebody is either cooked to a cinder or buried in ice . . .
Some possibilities: my atmospheric profile is only defined up to 80000 feet (24.5km). Could it be that, say, Plexus assumes the atmosphere above that (such as it is) is empty space (giving higher radiances) while DIRSIG gives it some set of default values (lower radiances)?
The profile altitudes are defined as AGL (I'm pretty sure), but DIRSIG might interpret them as MSL values (giving a thinner atmosphere at 300m, but that should make DIRSIG radiances higher, not lower.)
I'm amazed by the congruity between Plexus and Planck's equation. Plexus is clearly inserting the absorption bands, but shouldn't skylight be torquing the envelope?
How many top-of-the-atmosphere radiance spectra can MODTRAN choose from? Newkur.dat is one of these, supposedly, but are there others, and is DIRSIG using a different one?
It is certainly possible that my conversion from W/cm^2/str/cm^-1 (the Plexus output) to W/cm^2/str/micron is wrong, but if it is, nobody has said why.
radiance = [10000./(radiance(:,1)) radiance(:,2).*radiance(:,1)];
radiance = [radiance(:,1) radiance(:,2)./radiance(:,1)];
Thanks for your help.
Created attachment 218 [details]
ENVI spectrum of material type 16 using the MLS standard tp5 file.
Okay, here is some more information. I defined a material type with uniform emissivity of zero and placed it in the megascene image. This attachment shows the resulting radiance spectrum (taken from ENVI) after the image is created. As you can see, the radiance values are in line with results I obtain from PLEXUS. (Although I still have PLEXUS problems I am trying to reconcile.)
The next image was generated using a user-defined atmosphere in the .tp5 file.
Created attachment 219 [details]
ENVI spectrum of material type 16 using the user-defined atmosphere in the tp5 file.
This spectrum was taken from Megascene (material type 16) using a user-defined atmosphere (subsequent attachment). As you can see, while it's shape is similar, the values are an order of magnitude lower. This matches the DIRSIG results we showed earlier. Other aspects of the image generation (geometry, etc) are identical.
Created attachment 220 [details]
User-defined tp5 file used to generate the lower-magnitude spectrum.
Created attachment 221 [details]
Standard MLS tp5 file as created by the DIRSIG simulation editor.
Comment on attachment 221 [details]
Standard MLS tp5 file as created by the DIRSIG simulation editor.
I should note that when this tp5 file is fed to Modtran (v5) directly, it generates an error. I assume that this file is parsed by DIRSIG before its data is used by MODTRAN during the make_adb run.
Thanks for the updates, I will look into it this week.
I'm calling in the troops on this one.
Can we get your platform file?
Created attachment 222 [details]
Modified AVIRIS platform file
Attached is the platform file (or one very much like it; I'm presently on vacation and only have the copy of the file that was on my laptop). The only change I made it it was to reduce the image size from 1024x1024 to 512x512 in order to reduce execution time.
Created attachment 223 [details]
I recently put together a 3-slide Powerpoint presentation for Dr. Chris Borel that outlines the problems I'm having. Since we've recently brought in some other folks to look at it, maybe it will be useful
Thanks for the platform file.
Can you explain what the "Actual" is on your plots? Especially since it appears to have an atmosphere with no water vapor in it?
I am busy running a DIRSIG training this week, but we have a few of us looking into it right now. Thanks for your info.
Sorry, I should have said specifically that the graphics are described in detail in the "note" section of the powerpoint file. To answer your specific question, "actual" is a plot of newkur.dat, an irradiance file contained in the modtran release. It's use is specified in LSUNFL parameter of Card 1A. I should note here that I ran Modtran directly on several of the LSUNFL options with only trivial differences in the results, certainly nothing of the magnitude we're seeing here.
my 2 cents
in raw capture mode: the radiometric units are W/cm2/sr/um
in basic capture mode: the units are W/cm2/sr w/o channel normalization
the units are W/cm2/sr/um w/ channel normalization
my understanding (as a user) has been that the channel normalization
divides out by the area under the spectral channel response curve to get
units back into W/cm2/sr/um, which essentially treats the channel response then as an averaging function to some extent.
The first problem that we have is some confusion with your plot. You said that "Actual" is "newcur.dat" from the MODTRAN DAT directory. And that file is actually an exo-atmospheric (top of atmosphere) solar irradiance (not radiance). So plotting it with these other radiances isn't valid since they are not the same units (or radiometric concept).
Knowing that "Actual" is the top-of-atmosphere solar irradiance, the fact that your PLEXUS curve is greater makes us question the units of your PLEXUS curve. A quick calculation would be:
Radiance = Irradiance * 1.0 / PI * tau
where Irradiance is your exo-atm solar curve, the 1/PI term is the reflectance into a given direction from your unity reflectance Lambertian target, and the tau term is the atmospheric transmission (which is always less that 1, too). From that you can see that a reflected radiance could never be greater than the source irradiance. That PLEXUS curve doesn't add up (so to speak).
More importantly, we have done hand calculations which have given us confidence that the DIRSIG curve is consistent with the "Actual" irradiance curve and your unity Lambertian reflectance scenario.
Thank you for the time you are devoting to this problem. You are of course correct about the newkur.dat file, and I apologize for not having been more specific. Here is the full matlab code I used to condition the newkur.dat to generate the graph as it appears in my slide:
After loading the newkur data into a 49951x2 variable called "radiance":
radiance = [10000./(radiance(:,1)) radiance(:,2).*radiance(:,1)];
radiance = [radiance(:,1) radiance(:,2)./radiance(:,1)];
% Consider only AVIRIS wavelengths
radiance = radiance((radiance(:,1) > .36 & radiance(:,1) < 2.55),:);
% The reflected light is scattered over 2pi steridians
radiance = [radiance(:,1),radiance(:,2)/(2*pi)];
So basically I convert the units from W/cm^2/cm^-1 to W/cm^2/micron, and then assume lambertian reflectance over a hemisphere; hence dividing by 2pi steridians. I notice that you divided by only pi steridians; I'm not sure I follow why this is. Also, I set tau = 1 for the purpose of removing atmospherics.
Can you tell me what hand calculations did you use to verify the DIRSIG graph? And just to clarify, are you speaking of the graph when using one of my LEEDR atmospheres, or the graph when using the default mls atmosphere?
Created attachment 224 [details]
newkur irradiance and "scaled" version
This is Adam Goodenough, and I'm still trying to figure what exactly you are
trying to compare. Let's establish some upper limits first, since I'm stuck on
how your PLEXUS/Theoretical values can be what they are. Once we're on the same page, we can certainly talk about what DIRSIG is doing.
I just grabbed the newkur.dat file from Modtran and verified that you are
indeed plotting this exo-atmospheric spectral irradiance divided by two pi as
the "Actual" curve in your plot (I attached a plot with that same curve and the
exoatm irradiance as well).
So, my first question: Given that we agree about the exoatmospheric irradiance
magnitude and the units, how can either the Plexus or theoretical values be
correct? Ignoring reality for the moment, let's say that the sun is at a
perfect zenith angle such that the exoatmospheric irradiance is defined
perpendicular to the surface of the earth where the reflector is, assume there
is absolutely no atmospheric interaction (i.e. no losses due to extinction, no
scattering, etc...). Even then, the irradiance incident on the surface is
reflected diffusely into the hemisphere above the surface (we can discuss the
pi vs. 2pi issue later) -- in either case the magnitude of the retrieved
reflected radiance will be less than the magnitude of the incident irradiance
(we can hash out the details of the irradiance integral and Lambertian
reflectance definition later if you want as well).
So given the upper limit case, where we really don't have an atmosphere (which
introduces losses, and re-distributes energy off angle) and are purely
reflecting the raw exoatmospheric irradiance from an (unlikely) nadir incident
direction -- we can still only get a retrieved radiance that is smaller in
magnitude than the initial irradiance. Yet the PLEXUS and Theoretical curves
are higher than the values you computed yourself for this very simplified
newkur.dat irradiance reflectance and look to be on par with the raw
exoatmospheric irradiance values.
Any ideas for what is going on with those two curves? The fact that these are
the same magnitude as the exoatmospheric irradiance worries me a great deal...
Allan, sorry for misspelling your name in the last comment.
"I'm still trying to figure what exactly you are trying to compare."
Broadly speaking, I am attempting to reconcile the output between DIRSIG on one side and either PLEXUS or command-line MODTRAN (hereafter CCM) on the other side.
"[H]ow can either the Plexus or theoretical values be correct?"
The "theoretical" values were generated with the following code:
% First, find the "brightness" of the sun
T = 5800; % Sun temp (degK)
c = 2.998e8; % Speed of light (m/s)
kb = 1.381e-23; % Boltzmann's constant (J/K)
h = 6.626e-34; % Planck's constant (J s)
lambda = linspace(.36e-6, 2.55e-6, 1000); % Frequency range
B = 2*h*c^2./lambda.^5./(exp(h*c./(kb.*lambda.*T)) - 1); % W/m^2/str/m
B = B/10^10; % Convert to W/cm^2/str/micron
lambda = lambda*1e6; % Convert wavelengths from meters to microns
% Next, find the brightness of light reflected from speculon
r_sun = 696e3; % Radius of sun (km)
R_se = 152.1e6; % Mean distance from sun to earth (km)
B = B*(r_sun/R_se)^2; % * fraction of sky dome filled by sun, squared
I'm aware that this calculation waves away a lot of complexity: the scene geometry and skylight are the two obvious factors. And you are correct that if newkur.dat is used, then any values above it would be impossible. I include the graph only because my "theoretical" calculation were so damn close to the Plexus output, not counting the absorption bands, that I was convinced I'm on to something.
I'm going to try to go at this from a different direction. Can you post the build date of the dirsig version you used to generate the plot? (it's in the first few lines of output when you run dirsig). Also, can you verify that the actual platform file that you used has the per-channel normalize flag set to "true"? Thanks...
Build Date: Jan 8 2010 11:07:28
. . .
From the .platform file:
<channel normalize="true" shape="gaussian" name="Channel #1" >
<polarizer type="none" />
<dataoutput type="total" />
. . .
Created attachment 225 [details]
comparison between hand calc and dirsig radiance
Thank you for the build date and verification of channel normalization. There was a problem with the normalization flag not being read in correctly in an earlier, unreleased version and I wanted to make sure that we weren't dealing with that issue somehow. We are not, so I went back and did some quick calc comparisons with DIRSIG so we can discuss exactly what is going on.
I'm just going to describe what's happening in DIRSIG in terms of the radiometry of the solution. I'm still not sure what is going on with the PLEXUS/theoretical curves, but hopefully we can reach agreement on what the answer should be and what you are comparing to.
I've plotted the exo-atmospheric irradiance curves in the attached plot. The red curve shows the actual values used within DIRSIG for this particular run and the purple plot is the newkur.dat file (the conversion to wavelength units is shown). This just demonstrates that dirsig is using what we'd expect for exoatm data.
I've setup a simple scene consisting of a flat "Lambertian" panel with an emissivity of zero. The Lambertian property refers to how energy is emitted/scattered from a surface -- specifically, the intensity distribution follows Lambert's cosine law, I(theta)=I(0)*cos(theta). The consequence of this property is that the observable, reflected radiance, L, measured from any direction is constant, i.e. the observed radiance is invariant to observer location (let me know if you want the equations for this).
Lambert's cosine law does not say anything about the weighting of incident irradiance, just the distribution of exitant radiance (i.e. it is not a bi-directional reflectance distribution function). So we characterize Lambertian reflectance by the scalar ratio between the total incident irradiance (E) from all directions and the exitant irradiance (exitance, or M) -- a constant value (in this case 1) is sufficient to parameterize the reflectance at each wavelength. In other words, R = M/E.
Since we want radiance out, we can derive a relationship between exitance and exitant radiance. This works out the L=M/pi (this is where that single pi comes from and I'd be happy to show the math if you'd like). So if we know the total incident irradiance (E) then the output radiance (L) is: L = R*E / pi.
We also have to consider losses of energy to the atmosphere due to scattering/absorption which we can collectively characterize as a source->surface path transmission and a surface->observer path transmission (t1 and t2, respectively). This brings our quick calc equation to: L = t2*R*E/pi, where E incorporates t1 and the transmittance is taken from modtran.
The only thing that's left is to figure out what E, the total incident irradiance is. For simplicity, we can leave out the downwelled irradiance (i.e. solar irradiance scattered by the atmosphere towards the surface), since the contribution is much lower than direct, solar irradiance for most of the spectrum. We know the exoatm irradiance, Eo, onto a surface perpendicular to the sun-earth vector so the total (direct) incident irradiance is E = t1*Eo*cos( dec ), where dec is the declination angle of the sun (for my plots, this is about 66 degrees) and the cosine accounts for the relative orientation of our surface and the virtual plane perpendicular to the exoatm vector that defines it.
Plugging in, we get L = t1*t2*R*Eo*cos(dec)/pi, where t1, t2, and Eo are all values computed by modtran, R = 1, and dec is given by the position of the sun at the time of simulation.
To verify what dirsig was doing, I plotted a curve of the computed radiance using that exact equation with t1, t2, and Eo stolen from modtran output (the green curve in the attached plot). The result from dirsig is also plotted (the blue curve) and it matches the hand calc except for the following differences:
1) the hand calc is too low in the vis (particularly in the blue) since we are not incorporating downwelled irradiance from the sky in the hand calc.
2) the dirsig solution is much smoother which is because we are not incorporating the (normalized) spectral channel responses into our solution, which effectively convolves the data with a gaussian function.
Based on this, and other similar comparisons, I'm going to remove this bug from the open problems list. However, I'm completely open to keeping this discussion going if you want to keep sorting things out (and I'm definitely open to disagreement with what I just laid out as well).
I'm removing this bug from the open issue list for reasons given in the previous comment. However, please feel free to continue the discussion...
A couple of points:
- I would indeed like to see the formula by which we convert irradiance to radiance by dividing by pi. I will concede, however, that this is what modtran is doing.
- By "declination", do you mean the angle between "straight up" and the sun? Because 66 degrees seems awfully low in the sky for the sun to be at 1 p.m. in June! Plexus is giving me an angle of 23 degrees for this geometry.
That said, and considering these differences, I am satisfied that your DIRSIG curve comes out to be pretty close to my Planck/PLEXUS (MODEL=18) curve from the earlier slides. These are also identical to the results I am getting from both DIRSIG when using the default mls.tp5 file and from command-line modtran (CLM) when using a standard atmosphere with this geometry.
The problem lies with my user-defined atmosphere. The results I get from both DIRSIG and CLM are radically different, by orders of magnitude, not only from the results for a standard atmosphere, but from each other; in other words, given the same .tp5 file, CLM and DIRSIG produce different results. It is noteworthy that the CLM results loosely match the PLEXUS results for MODEL=9, even when PLEXUS is using a standard atmosphere.
So . . . my slender claim on your continued attention is that DIRSIG is parsing the .tp5 file differently than MODTRAN and/or making different assumptions that we don't know about. Scott and I went very carefully over the .tp5 file formatting a bug or two ago. We were satisfied at the time that it was done correctly; however, if you could have a look at it, I would appreciate any insight you could offer.
Created attachment 226 [details]
lambertian radiance/irradiance relationship
>- I would indeed like to see the formula by which we convert irradiance to
>radiance by dividing by pi. I will concede, however, that this is what modtran
Please see the attached pdf (a lot easier for me to write equations in latex than in here). I've derived this relationship based on the radiometric definitions and Lambert's cosine law.
>- By "declination", do you mean the angle between "straight up" and the sun?
>Because 66 degrees seems awfully low in the sky for the sun to be at 1 p.m. in
>June! Plexus is giving me an angle of 23 degrees for this geometry.
Yes, that would be a bit low for June :) Sorry, that run I posted was actually done at 13:00 in December. I was doing a few runs (to just check against the hand calc) and that happened to be the one I was checking at the time. Sorry for confusing the situation.
>however, if you could have a look at it, I would appreciate any
>insight you could offer.
I will take a look at the consistency of the internal tp5 model and let you know what I find....
There are definitely some formatting issues in both your tape 5 and in our code that reads it in. The trouble on dirsig's side appears to be the use of sscanf, which really doesn't like handling structured whitespace and some non-robust input formatting. This isn't normally my area, but it's not a terribly difficult (though tedious) fix on our end. At the same time, I'd check over your tp5 as I definitely saw a few spacing problems.
I've yet to confirm my results with my committee, but I have isolated the cause of the anamolous output of MODTRAN (and, my extension, PLEXUS).
It turns out that MODTRAN doesn't behave well if the atmospheric profile is defined for altitude values below that of the target. In my case, the first altitude of my profile has always been zero, while the target altitude (ground) has always been 300m for the MegaScene image.
This suggests a couple of ways to fix my immediate problem, obviously, but I am still concerned that DIRSIG is giving results different than CLM. DIRSIG parses its .tp5 file before feeding it to MODTRAN, and I want to be as sure as possible that the altitude values mean the same thing to MODTRAN whether they are coming from DIRSIG, PLEXUS, or directly from a .tp5 file.
I can tell you two things:
1) going through dirsig, you shouldn't need to worry about the first altitude being at zero as dirsig will automatically resample the data you give it (I was actually wondering why the tape5 output by dirsig for your tp5 reset the minimum altitude to 0.3 instead of zero).
2) i can currently guarantee that dirsig's parsing is currently messing up at least card 1A and possibly a few others. this is due to some non-robust handling of some of the new parameters that have been added (which may not be given in the input tape5 file). i've initiated a separate bug for this and i've discussed some long term solutions with scott that we are moving forward on (though probably not as quickly as you might need). my suggestion for the short term would be to explicitly set each parameter in the last version of modtran4 which I believe is the format that is being used. if you want to send me any tape5 you have, i can verify for you whether or not it will be read in the way you expect. in the meantime we'll be moving towards a better way of handling the changing card formats...
Could you clarify "resample"?
For instance, I'm comparing my .tp5 input file with DIRSIG's tape5 output file, and I notice the following:
- The values I give at 0km are removed.
- The values for pressure and mixing ratio I give at .5km are assigned to the altitude at .3km in the output.
- The value for temperature I give at .5km almost matches the temperature at .3km in the output, but not exactly.
I will test DIRSIG against a .tp5 file containing profile data with MSL altitudes, starting at .3km; however, you can kind of see how DIRSIG's handling of atmospheric levels below the target altitude isn't really consistent with either an MSL or AGL interpretation.
Perhaps I used the term resampled too loosely and didn't address your actual problem as well as I should. I only meant to indicate that the DIRSIG tape5 model was apparently making sure that the lowest altitude wasn't below the altitude of the target if you did in fact set it as such. Whether that correction is valid or not, I have no idea (it sounds like it might be a superficial "fix" at best). The best solution seems to be providing data starting at the target altitude for any run.
... a year later
Finally finished a flexible tape 5 model that is driven by an xml-based version format document (to allow easily changing between versions of modtran). I'm attaching the debug log file created when reading the griffiss50.tp5 file in. There are minor differences in the tape 5 file the is written out, but are just minor formatting changes (I'm using the format codes in the MOD4 manual, aside from switching the card 4 wavelength/wavenumber format from F10.0 to F10.4 to support wavelengths).
Created attachment 241 [details]
debug log from the new tape 5 model tool