I found the routines at http://celestrak.com/software/vallado-sw.asp, associated with Vallado's book helpful for obtaining the needed positions. For satellite ephemerides in ECI J2000, I could convert them to ECEF (WGS84 xyz for example), and there are routines that can convert those to geodetic latitude, Longitude, altitude; I used Vermielle's algorithm for that step (Journal of Geodesy (2002) 76:451-454). 

So, given the pointing vector from the ground to the sun, and the pointing vector from the ground to the spacecraft, one can determine the angle between the reflected (specular) vector and the view vector, assuming a flat surface. In reality, because there are waves on surfaces, a lot of ocean remote sensing stays pretty far away from the specular vector - something like more than 20 or 30 degrees. For those purposes especially it is more tolerant of small errors. For glint from small bodies of water with very little fetch (ponds, swimming pools) I suspect you may need to be more careful, but one could calculate using a few different algorithms (if in doubt) and determine the differences in the answers to se if it really matters. 

Thank you Marcos, it was very instructive. I will look into the references you sent to develop a specific solution to my problem.

I have done a lot of reading and coding in the last 24 hours, and now I understand a little bit more about the issue. I manage to combine the information on the U.S. Naval Observatory's Almanac for Computers 1991 with that provided by Hapgood (1992) to find what I was hoping to be a suitable formulation; the results are still a little different from SPICE, but perhaps that is because J2000 is actually ICRF (which has a close alignment, but not equal). I will see the routines you recommended.

I didn't had time yet to test the sunglint calculations, having spent a lot of hours to learn about time measurements and reference frames and trying some test codes. But shortly I will give some feedback on that; as a beginner to the subject, I expect to find some more cumbersome issues.

Well, I could not do it... I'm sending a numerical example. The sensor in question is TOMS/EP (Total Ozone Mapping Spectrometer - Earth Probe). I tried to use simple trigonometric relations with the solar zenith, solar azimuth, viewing zenith and viewing azimuth, but azimuth data do not follow the definitions on the documentation (at surface, degrees from north, clockwise positive, -180 to 180)... viewing azimuth is always 0 and solar azimuth has a clear discontinuity at the nadir (column 18 of the 35 columns in a scan row), not showing the expected smooth change of solar azimuth along a scan row. I could correct this adding 180 to columns 19 to 35, but could not correct the nadir value (so in the figure "Result_1" attached, the center of each swath have a strip with very different results).

Then, I tried the formulation on the Mailhe paper I sent before, but my results were even worse ("Result_2", attached). I'm sending a numerical example below. The points A and B are over the Atlantic ocean (WGS84), and I know, from TOMS L3, that point A is within the sunglint region, and point B is not, but my glint angle for both scenes are above 68º... Any ideas?

*Time of observation:

2001-01-01 12:12:53 UT

*Earth-Sun unit vector, geocentric (ECEF x, y, z):

-0.2577401, 0.8836004, -0.3909224

*Spacecraft geodetic position (alt(m), lon, lat):

766000, -10.743, -27.129

*Earth-Spacecraft unit vector, geocentric (ECEF x, y, z):

0.875469, -0.1661022, -0.4538326

*Surface geodetic position (lon, lat):

Point A: -10.01, -26.768

Point B: -2.79, -25.307

*Earth-Surface unit vectors, geocentric (ECEF x, y, z):

(no correction made for geoid height, +- 100m)

Point A: 0.8804386, -0.15540352, -0.4479706

Point B: 0.9040616, -0.04405782, -0.4251253

*Glint angle from Mailhe formulations:

Point A: 70.23º

Point B: 68.99º

I'm sending a L3 TOMS data that excludes sunglint regions for comparison with my calculations of sun glint. The portions without data (in a north south configuration) that center in the equator are inter-swath spaces; the sunglint areas also have a north-south orientation but are centered in the southern hemisphere (in the middle of each swath). "Result_1" is the sunglint angle calculations using at surface geometry provided by TOMS, and "Result_2" is using Mailhe formulations. In Result_1, a circular feature is present in the middle of the sun glint area, but is circular and does not have north-south orientation as in the TOMS L3 data; Result_2 seems far wrong...

By "earth-sun" unit vector and "earth-spacecraft" unit vector, you mean from the center of the earth to the sun or spacecraft?

Yes.

I noticed that in your earlier reply, you transformed from ECEF to geodetic; but it is not better to use ECEF as the center is the same for all latitudes? What I understood is that geodetic latitudes are the polar angle of the normal to the surface with the equator plane, and do not pass through the Earth's center, as the Earth is not a perfect sphere.

ECEF is better for all vector related items, and certainly for the glint direction calculation.  Geodetic is useful for some things.

I'll take a look at this a bit later today. 

I haven't put numbers in your equations, but could you verify that the steps listed in the accompanying document is what you are doing?

Dear Marcos,

Thank you very much for your help. The trigonometry of the problem was clearer from your formulation. Even though, a few attempts latter I was still arriving at wrong values - it turned out that there was some bugs in my code...; the Vallado's book did help a lot to solve those out. I'm sending a final (and at last right) result (the dots on the image are the two test points from the previous reply). I'm in your debit.

Best regards,

Glad I could help! You're very welcome. 

Alternate approach:  the surface normal from the glint point must intersect the line between the sun and the satellite; assume there is only one such point, call it X

Implement the following algorithm:

1) set ToSun = Sun location, set ToSat = satellite location

Start binary search loop:

2) set MP = midpoint between ToSun and ToSat (SPICE routines VSUM and VSCL)

3) SP = point on the surface closest to MP (SPICE routine NEARPT)

4) SN = MP - SP (= surface normal at SP; SPICE routine VSUB)

5) INC = angle between (ToSun-SP) and SN; EMI = angle between (ToSat-SP) and SN (SPICE routine VSEP)

6) if INC==EMI (or within some tolerance), SP is the glint point, exit loop

7) If INC < EMI, then X is between MP and ToSat, and set ToSun = MP; else set ToSat = MP (SPICE routine VEQU)

8) Return to step 3

All vectors except SN have earth center as their origin and are in earth body-fixed reference frame (SPICE routine SPKEZR with frame IAU_EARTH).

Steps 2 through 8 implement a binary search that will converge in logartihmic time.

Code will have as many or more lines loading SPICE kernels and getting the time input than solving the problem.

Thank you Brian for an approach using SPICE. I'm now testing the validity of a hypothesis on TOMS's data (for which first I must exclude sunglint areas...), so I'm using approximate calculations and an interpreted language (R); but if it shows promising results, I will have to migrate to C to gain efficiency, in which case SPICE and/or SOFA routines could be used.

I didn't had the time yet to explore the available routines in SPICE and SOFA; so I would appreciate your recommendation on which library to use.

I am a bit of a spice evangelist, and I never heard of SOFA, but my first question would be, what format are the TOM ephemerides in, and is it easier to get them into SOFA-compatible or SPICE compatible format?  I see SPICE can take TLEs and generate positions, but I am not sure of the center or reference frame of those results.

In looking back I realize you are not trying to get the specular point, but for each data pixel you want to know what is the angle between the emission vector (the vector from the pixel intersection with the surface to the instrument) with the vector that is specular reflection of the sun at that same surface point.  Then you will make that calculation for every pixel in each observation.  Is that correct?

On second thought maybe the specular point is useful:  for culling which pixels need to be tested, because any pixel pointing away from the specular point by more than 30deg cannot have a glint angle less than 30deg.

This repo:

  https://github.com/drbitboy/glint

does what you want in file glint.py.  Most of the code parses the input and sets up the variables.

Five SPICE calls at the end do all the work.

You need to put the two SPICE kernels in the same directory as glint.py

Yes Brian, the routine calculates for each pixel the angle between the specular reflection and view vectors. TOMS has a large IFOV (50x50 to 70x140 km) so the amount of data to process for this test is manageable; nevertheless a more economical approach can be useful if a large amount of data come to be processed. "SOFA" stands for Standards of Fundamental Astronomy, and is a library maintained by the International Astronomical Union (IAU), but I don't know much about it: http://www.iausofa.org/ .

I don't have the ephemeris but the position (altitude, longitude and latitude) of the spacecraft is given for each scan line in the TOMS's level 2 data; it seems that SOFA has routines for geodetic transformations, so it should not be a problem. I didn't saw those in SPICE, but again, I did not focused on the libraries.

I had not seen that you posted routines on github. Thank you very much for your time and effort! I will look into that. My best,

Alexandre,

Check the github repo again; I wrote a Python script (tag.py TOMS Add Glint) that

- reads the L2 .he5 HDF5 files

- calculates the glint angles for each pixel

- copies the data from the L2 files into new HDF5 files (_glint.he5)

- adds the glint angles to the new file (.../GeoLocation Fields/GlintAngle)

On my laptop it processes about four files per second, so it should take a few hours to process the entire data set.

At some point you should check my glint angle calculation.  I use a SPICE kernel to get the location of the Sun in the IAU_EARTH frame, so there may be an inconsistency.

Best regards,

Brian

  https://github.com/drbitboy/glint

I ran my python code on the 37,098 TOMS/EP .he5 files from the FTP site; it completed in under 4h.

(Whether my implementation of the glint angle calcuation is correct is an entirely different matter;-)

Diagram of my how I calculate the specular reflection angle, which is a prerequisite to calculating the glint angle:

https://docs.google.com/presentation/d/18KIez0HScqvvw2HCDe-clH9S9On9p45wRNUU9TO2THM

Dear Brian, thank you again for all your attention and involvement. I must say that I did not expected that much enthusiasm! But very much like you, whenever I have the time, I try to easy the way for someone not familiar with the tools I use. There are no "royal roads" to any knowledge, but as a beginner programmer, your codes will certainly help me to advance faster (for example, how to use SPICE library).

You are very welcome, and thank you for the challenge.  I'm a bit of a SPICE evangelist and while I am not part of NAIF at JPL, the NAIF group at JPL is always looking to reach out and I like to support them (as they have supported me).  I think this particular problem would make an excellent case study for using SPICE, and I will probably bring it to their attention at some point.

The funny part for me is that the code that solves the problem is usually just a few SPICE calls, but the steps to get to those calls is what takes up most of the code.

Brian