Batch processing script, MODIS chl for inland lake

By allyn Date 2012-02-22 10:58 Edited 2012-02-23 08:55

Hello,

This forum was extremely helpful the last time I had a problem, so here we go again.

I am trying to get chlorophyll data for Lake Kivu in the African Great Lakes Region, following a procedure established in Horion et al (2010) for Lake Tanganyika. http://perso.fundp.ac.be/~jpdescy/Horion%20et%20al%202010_RSE%20114%20781-7911.pdf.

I found a script on this forum that processes from L1A to L2, so I used a modified version of that (below).
I ran my script last night for files in /mnt/hgfs/seadas_shared. For some of my original L1A files, I obtained the original .hdf, a GEO file, L1B, L1B.anc, L2, and text files. For others, it gave me the hdf, geo, geo.met, geo.pcf, l1b, l1b.anc, l2, text, and met files. Question one: Why would this be different? Is this an error? What should I obtain? Note: I ran out of storage space, so some files weren't processed at all - not sure if this played a role.

I displayed a few images, and I get a colorful display when looking at products from the original L1A file and the L1B, but for L2, I'm getting -32767 as the raw data value (same for geophysical data value for chl) for everything in the spreadsheet, whether it is a chlorophyll product or nLw product. Question two: What would cause this? Is there a general problem with using SeaDAS ocean color for Lake Kivu? Or was there a processing error or user error? Or could it just be clouds...? Note: I haven't yet figured out how to use and interpret the L2 Flags properly.

Lastly, because of storage capacity, my L1A files are stored in two directories on my host Windows OS, one on the external hard drive and one on an internal data drive. I re-ran the same script (and then just the first part of it) for /mnt/hgfs/seadas_shared-2/*. This time, I got the following error in the terminal (below). Question three: what happened here?

I haven't used linux or seadas before, and none of my colleagues have ever used seadas, so your help is greatly appreciated and truly indispensable to me completing this project.

Allyn

SCRIPT

source /home/ack/Desktop/seadas6.3/config/seadas.env

WEST=28.8
SOUTH=-2.8
EAST=29.8
NORTH=-1.6

#for FILE in /mnt/hgfs/seadas_shared/*.hdf
for FILE in /mnt/hgfs/seadas_shared-2/*.hdf

do
echo $FILE
#BASE= echo $FILE |awk -F. '{ print $1 }'
BASE=$FILE
GEOFILE=$BASE.GEO
L1ASUB=$BASE_sub.L1A_LAC
GEOSUB=$BASE_sub.GEO
L1BFILE=$BASE.L1B
L2FILE=$BASE.L2

echo "GENERATING NEW GEOFILE..."
/home/ack/Desktop/seadas6.3/run/scripts/modis_GEO.py $FILE -o $GEOFILE

#Rest of script was commented out...

done

TERMINAL

ack@ubuntu:/mnt/hgfs/seadas_shared-2$ /home/ack/Desktop/script.sh
/mnt/hgfs/seadas_shared-2/MYD01.A2008123.1135.005.2012039111623.hdf
/home/ack/Desktop/script.sh: line 22: //: is a directory
GENERATING NEW GEOFILE...
Traceback (most recent call last):
File "/home/ack/Desktop/seadas6.3/run/scripts/modis_GEO.py", line 115, in <module>
modis_env(m)
File "/home/ack/Desktop/seadas6.3/run/scripts/modules/modis_utils.py", line 248, in modis_env
self.start, self.stop, self.sat_name = modis_timestamp(self.file)
File "/home/ack/Desktop/seadas6.3/run/scripts/modules/modis_utils.py", line 126, in modis_timestamp
meta = readMetadata(arg)
File "/home/ack/Desktop/seadas6.3/run/scripts/modules/MetaUtils.py", line 73, in readMetadata
if re.search('PRODUCT',text[0]):
IndexError: list index out of range
/home/ack/Desktop/script.sh: line 26: //: is a directory
/home/ack/Desktop/script.sh: line 30: //: is a directory
/home/ack/Desktop/script.sh: line 35: //: is a directory

By @sean Date 2012-02-24 21:02

> I ran my script last night for files in /mnt/hgfs/seadas_shared. For some of my original L1A files, I obtained the original .hdf, a GEO file, L1B, L1B.anc, L2, and text files. For others, it gave me the hdf, geo, geo.met, geo.pcf, l1b, l1b.anc, l2, text, and met files.
> Question one: Why would this be different? Is this an error? What should I obtain? Note: I ran out of storage space, so some files weren't processed at all - not sure if this played a role.

Yes, running out of space can cause all sorts of odd behaviour.

> I displayed a few images, and I get a colorful display when looking at products from the original L1A file and the L1B, but for L2, I'm getting -32767 as the raw data value (same for geophysical data value for chl) for everything in the spreadsheet, whether it is a
> chlorophyll product or nLw product. Question two: What would cause this? Is there a general problem with using SeaDAS ocean color for Lake Kivu? Or was there a processing error or user error? Or could it just be clouds...? Note: I haven't yet figured out how
> to use and interpret the L2 Flags properly.

Before processing in batch, it is best to work the problem out on a single file. Lake Kivu is significantly above sea level. To process data over such a body of water, you will need to enable terrain height correction when computing the geolocation files. This will both provide the correct location (and thus viewing geometries necessary for proccessing) and account for the altiude (important when computing Rayleigh radiance). When calling modis_GEO.py, use the --enable-dem flag - but first make sure you have installed the Digital Elevation Maps.

You should also get familiar with the L2flags - they can be quite helpful in understanding the data.

> Lastly, because of storage capacity, my L1A files are stored in two directories on my host Windows OS, one on the external hard drive and one on an internal data drive. I re-ran the same script (and then just the first part of it) for /mnt/hgfs/seadas_shared-2/*.
> This time, I got the following error in the terminal (below). Question three: what happened here?

Not sure. Try modifying the script to include the --verbose option in the call to modis_GEO.py - there may be a clue there.

Regards,
Sean

By allyn Date 2012-02-29 11:29

Sean,

Your response was very helpful. Thank you. I had already installed the DEM files, and I added -d to the modis_GEO.py processing step and -v to all steps, so in theory, the terrain correction should have worked because I saw no errors in the terminal output. And yet, when I looked at the L2 Flags, I found a CHLFAIL of 100%.

I played around with it a bit and tested a batch process to see if this was just a problem with my two test images. Batch processing yielded in the terminal Percent valid data (0.00) is less than threshold (95.00) for the geogen step. This problem reappeared when I tried to reprocess the single file, so I followed Gwyn's advice in this post: http://oceancolor.gsfc.nasa.gov/forum/oceancolor/topic_show.pl?pid=18222#pid18222

1. I ran update_luts.py aqua --verbose. All good.

2. I checked uname -a ; ls -l $SEADAS/run/data/modis/static/de200.eos `which geogen_modis`
RESULT (which I don't understand):
Linux ubuntu 2.6.32-28-generic #55-Ubuntu SMP Mon Jan 10 21:21:01 UTC 2011 i686 GNU/Linux
-rwxr-xr-x 1 ack ack 7255484 2011-10-31 10:43 /home/ack/Desktop/seadas6.3/run/bin/linux/geogen_modis
-rw-r--r-- 1 ack ack 4764304 2012-02-21 04:07 /home/ack/Desktop/seadas6.3/run/data/modis/static/de200.eos

3. I deleted de200.eos.

When I reran the single file, thus rebuilding the ephemeris files for this (and further) reruns, I had solved the "percent valid" problem, but I was back to the same original issue of having a CHLFAIL of 100%.

Here are the Flags with non-zero results for one file:
Flag # 2:    LAND    8951 90.2319
Flag # 5:    HILT    9079 91.5222
Flag # 7:    COASTZ 969   9.7681
Flag # 9:    STRAYLIGHT    4013 40.4536
Flag #10:    CLDICE    2929 29.5262
Flag #16: CHLFAIL    9920 100.0000
Flag #28: SSTWARN 969   9.7681
Flag #29: SSTFAIL    8951 90.2319
Flag #31:    PRODFAIL    9920 100.0000

Obviously, my region is mostly over land, and I know that cloud cover could be an issue for most of my images, but is there anything else that could explain the consistent CHLFAIL flag of 100%?

My defined extract region (WEST=28.8 SOUTH=-2.8 EAST=29.8 NORTH=-1.6) actually does include all of Lake Kivu, right?

Thank you again for all of your help. If I end up with any results at the end, it will be thanks to you, your colleagues and this forum.

Best,
Allyn

By @sean Date 2012-02-29 14:31

Allyn,

I processed a recent granule over Lake Kivu (A2012036115000.L1A_LAC ).
To get retrievals over the lake, I had to modify the cloud_albedo threshold - seems the default of 0.027 is too strict for
this high altitude lake. I also unset the masking on straylight, given the small size of the lake, the straylight mask
set due to the proximity to land masked out much of the lake. So, try these parameters in your l2gen run:

proc_land=1 (this will adjust the Rayleigh correction for the terrain height)
maskstlight=0 (turns off the default masking on the flag stray light)
cloud_albedo=0.054 (double the default, MOST of the lake had a cloud_albedo of 0.032, or about 20% higher than the default)

Regards,
Sean

By allyn Date 2012-02-29 16:23

And we have some improvement!

My first test-file yielded a CHLFAIL flag of 93%, but since 90% of pixels are over-land, it looks like most of my over-water pixels should work. (I think?!) Next step: batch processing. I'll be sure to let you know if I run into any brick walls with that.

Thank you so much for your help Sean.

Allyn

By allyn Date 2012-03-02 17:22

Dear Sean,

For my next step, I need to batch process the exportation of my L2 Products as ascii files. This will allow me to continue my work in MatLab, for which I have much more support because all my colleagues use matlab, but I'm the first to try seadas.

Even after looking through this forum, I have only an inkling of how to do that. IE in terminal: seadas -em -b 'filename'. And in that file, do I need IDL or seadas-style commands? Can I still use For-do-done and echo, the way I did in bash? Or how do I modify it?

In my first attempt, I think I mix languages... and basically, it cannot possibly work, but here it is anyway:

for FILE in /mnt/hgfs/seadas_shared/*.L2 & $

do & $

echo $FILE & $
load $FILE ftype=$MOS & $

display & $
landmask, color=6 & $
out/display & $

out_ascii, fbuf=1 ll_flag=1 region=[28.8, -1.6, 29.8, -2.8] & $

vars=[1,2,3,4,5,6] & $
#Which variables do I really need?

prod_name='chlor_a,chl_oc2, chl_oc3, chl_giop, chl_gsm, cdom_index, K_490, nLw_443, nLw_488, nLw_547' & $
#I want to output each product so that I can play with them. I mostly want to determine which chl algorithm will best match my in-situ chl data.

write file & $
#Will this write a new file for each product?

clear_up & $

done & $

Thanks again and have a great weekend.

Allyn

By @sean Date 2012-03-02 18:25

Allyn,

You're close. There is a FAQ on how to do this.
The batch script needs to be IDL syntax, as it's using IDL. So, "for FILE in ..." would be "for i=0, n_elements(files-1 do begin"
The IDL variables do not use prepended '$'. Here is the example script described in the FAQ listed above:

for i=0,n_elements(filenames)-1 do begin & $
   filename=filenames(i) & $
   base=strmid(filename,0,14) & $
   ascii2=base+'_latlon.asc' & $
   print,ascii2 & $
   load,filename,prod_name='chlor_a' & $
   out_ascii, ascii2, band_no=1, region=[pt4,pt3], /line, vars=[3,4,2], $
   format='(F10.2, F10.2, F10.4)', titles=['Latitude','Longitude','CHL'], $
   /ll_flag & $
   clear_up & $
endfor

The "filenames" list can be obtained using the "file_search" function to get a
list of files from the operating system.

Sean

By allyn Date 2012-03-06 08:28 Edited 2012-03-06 12:21

First of all, thank you for the help on my IDL script, I'll work on that, but for now, I have another bug to solve in the first step for processing L1A --> L2.

The "test batch" of five files worked at the beginning but got held up later. I thought it was because processing took a while, and I needed to rerun update_luts. So I tried to integrate update_luts into my script by calling it at the beginning of my for loop. This didn't seem to solve the problem either. So I figured out how to set up update_luts as a cronjob (and remove it from the script), which I did successfully, but it did not solve the problem, so I'm not sure what is going on that is preventing the last 3 of 5 files from processing, except that the 3rd file, which processes first, takes way longer than it should to create the GEO file. Not sure why.

For the full batch, I redid (again) the whole step of deleting de200.eos so that SeaDAS regenerates the ephemeris and attitude files. Same results.

Ideas?

Thanks,
Allyn

Terminal error for test batch (The first error sometimes occurs, the 2nd occurs with or without the first)
Running validation test on geolocation file..
Percent valid data (0.00) is less than threshold (95.00)
/home/ack/Desktop/script.sh: line 28: //: is a directory
GENERATING SUBFILES SUB.L1A_LAC and SUB.GEO...
Locating pixel/line range ...
Error locating pixel/line range to extract.

First error in logstatus.hdf.geo file from test and full batches
PGS_EPH_EphemAttit():PGSEPH_E_NO_SC_EPHEM_FILE:44547
TOOLKIT error

Weird "bounding coordinates" in logreport.geo for test and full batches - not the coordinates I used:
Bounding coordinates
east: 3.141593
west: -3.141593
north: -1.570796
south: 1.570796

By @sean Date 2012-03-07 14:50

Allyn,

The update_luts.py script need only be run at most once a day - once a week or two is probably sufficient.

The de200.eos file needs to be created only once. It is a static file, but as it is an architecture dependent binary, we require the file to be created on the local system.
Unless you change OS flavors, you should not need to regenerate this file.

The error:

> PGS_EPH_EphemAttit():PGSEPH_E_NO_SC_EPHEM_FILE:44547

suggests a problem with the ephemeris file. It is possible the file was not sucessfully retreived.
Please make sure the download of attitude and ephemeris files was successful. If you run modis_GEO.py with the --verbose option
you will see the attitude and ephemeris files the program is using. You could compare the file size with the size presented on our
website for the files in question: http://oceandata.sci.gsfc.nasa.gov/Ancillary/Attitude-Ephemeris/<YYYY>/<DDD>/

If the file sizes differ, you can manually replace them or remove them and rerun modis_GEO.py with the --refreshDB flag.

Regards,
Sean

By allyn Date 2012-03-08 09:16

Sean,

The file is there, and when I redownloaded it, the size didn't change:

PM1ATTNR.P2002229.1000.003 of 452 kB

The only variability among all such files is that occasionally one ends in ".1200.001" instead of ".1000.003" (eg PM1ATTNR.P2002231.1200.003). The two that worked end as does the above. All are the same size.

I did as you suggested, but my results didn't change, even when adding -r and -v to the script. I must be doing something wrong. Is sufficient to add the option like that or does it have to literally be "-v, --verbose" in order to function? Is there some other script-writing or SeaDAS quirk that I should know?

Again, I very much appreciate your help and apologize for taking up so much of your time.

Allyn

Script command:

echo "GENERATING NEW GEOFILE..."
/home/ack/Desktop/seadas6.3/run/scripts/modis_GEO.py $FILE -o $GEOFILE -r -d -v

Terminal output:

GENERATING NEW GEOFILE...
Determining required attitude and ephemeris files...
Searching database: /home/ack/Desktop/seadas6.3/log/ancillary_data.db

Input file: /mnt/hgfs/seadas_shared/2002-test/MYD01.A2002229.1115.005.2012039114417.hdf
Start time: 2002229111500
End time: 2002229112000

Found: PM1EPHND.P2002228.1200.001
Found: PM1ATTNR.P2002229.1000.003

Creating MODIS geolocation file...
GEO version: 5.0.14 built on Oct 31 2011 (13:43:00)
scan: 0 out of 203 Thu Mar 8 00:37:07 2012
% ... The scan for this particular file takes forever, and it is the first file with which I started having problems.

By @sean Date 2012-03-08 17:20

Allyn,

The script accepts -r or --refreshDB (and -v or --verbose). You shouldn't put -r in the script, as it is only necessary if there is a need to refresh you local database of ancillary files.
I was able to sucessfully process our version of the L1A file causing you trouble (A2002229111500.L1A_LAC)

Try running the file through again, this time use the --log flag to retain the processing log files. Attach these files to you post along with the complete terminal output of the modis_GEO.py process.

Sean

By allyn Date 2012-03-13 09:38 Edited 2012-03-13 09:40

Sean,

I'm attaching the complete terminal, my script as run, and the Log files for the first file of three that I tried to process. Thank you for your help.

Best,
Allyn

Attachment: terminal (4.8k)

Attachment: LogReport.MYD01.A2002229.1115.005.2012039114417.hdf.GEO (1.0k)

Attachment: LogStatus.MYD01.A2002229.1115.005.2012039114417.hdf.GEO (212.0k)

Attachment: LogUser.MYD01.A2002229.1115.005.2012039114417.hdf.GEO (1.0k)

Attachment: script.sh (1.4k)

By gwyn Date 2012-03-14 12:35

Allyn -

The LogStatus files show lots of errors like:
Error returned by function PGS_CSC_quatRotate(sun vector)

We often see these when the Planetary Ephemeris file was built for the wrong platform. I see you've already tried rebuilding de200.eos; all I can think of is that something went wrong with the rebuild.

Do your platform and the Planetary Ephemeris binary match? Try:
uname -a; file $LIB3_BIN/ephtobin

By gwyn Date 2012-03-14 12:40

Also, you should remove the "//" lines from your script. That's a comment delimiter for C, but it's "#" for shell scripts.

By allyn Date 2012-03-19 10:45

Hi Gwyn,

I'm not sure what the results (below) mean, but in case this is relevant, I'm running Xubuntu 10.04 that I installed in a VA. It is the same OS as supplied for the seadas virtual machine, but I installed it separately after running into problems with my original seadas installation. (I did much better with the step-by-step instructions for installing seadas on a linux machine than I did with the instructions for installing seadas in a VA on a windows machine.)

ack@ubuntu:~/Desktop$ uname -a; file $LIB3_BIN/ephtobin
Linux ubuntu 2.6.32-28-generic #55-Ubuntu SMP Mon Jan 10 21:21:01 UTC 2011 i686 GNU/Linux
/home/ack/Desktop/seadas6.3/run/bin3/linux/ephtobin: ELF 32-bit LSB executable, Intel 80386, version 1 (SYSV), dynamically linked (uses shared libs), for GNU/Linux 2.6.9, not stripped

Next suggestion? Sorry, I'm a total novice, but I really appreciate your help.

Best,
Allyn

By gwyn Date 2012-03-19 13:09

I'm stumped. The executable and system are compatible, as they're both 32-bit.
Sean suggested that if you're short on disk space, partially-written files could cause a lot of problems.

By allyn Date 2012-03-19 14:00 Edited 2012-03-19 17:00

Alright. Well, I fixed the disk storage problem a while back by moving files into two different locations, each with 40-60 gigs of extra space. So I can't imagine that being a problem, unless copying&pasting the original .hdf files and/or deleting old files from aborted processing runs can somehow screw something up...? Could it be the files themselves?

I tried reinstalling seadas, and the process seems to have worked for the 1st file, which didn't work before. And by that, I mean that I can get it to display colors in the shape of Lake Kivu. And then it gave me 100% LAND error for subsequent images.

For all of them, even the first file, the LogStatus comes up with the same errors about altitude.

PGS_CSC_ECRtoGEO():PGSCSC_W_INVALID_ALTITUDE:35877
invalid altitude - probably indicates bad input data

I'm pretty determined to figure this out, despite my lack of computer expertise, so if you guys are willing to put up with more questions along the way, we might get this working at some point.

Thanks again for all your help.

Allyn

Attachment: LogStatus.MYD01.A2002229.1115.005.2012039114417.hdf.GEO (2.6k)

By @sean Date 2012-03-19 19:12

> Could it be the files themselves?

Yes. The geolocation process will modify the input L1A file. IF that process was somehow interrupted, the source L1A file can be left in a not so happy state.

When uncompressing the bz2 files, you could add the -k option which will keep the compressed file around, then once GEO is successful, you can delete the
compressed file (or process to L1B, then delete the uncompressed L1A - storing the compressed L1A should you wish to process again).

Sean

By allyn Date 2012-03-22 10:33 Edited 2012-03-22 10:40

I've never ran benchmarks, but the bathymetry image - at least with a landmask I can see the outline of the lakes.

After redownloading my original files and testing those, I got the same error in LogStatus as with the geo-process-modified files... and 100% land fail error on output.

PGS_CSC_ECRtoGEO():PGSCSC_W_INVALID_ALTITUDE:35877
invalid altitude - probably indicates bad input data

PGS_CSC_GrazingRay():PGSCSC_W_HIT_EARTH:35897
The line of sight strikes the Earth

Note:
1. I downloaded my original uncompressed L1A .hdf files from LAADSweb via FTP instead of the bz2 files from OCBG via HTTP. Does that make a difference?
2. Half of my files are stored on an external harddrive. I've seen that this can cause problems depending on permissions, but I didn't have problems at first, so I'd be surprised if changed somehow, unless the fact of reinstalling seadas... oh I really don't know.

What would you do to get this working - start again from zero?

Thanks again.

By @sean Date 2012-03-22 15:13

Allyn,

LAADS L1A files are the same as ours - we just rename them and compress them...to verify that I grabbed the
MYD01.A2002229.1115.005.2012039114417.hdf file from LAADS. It processed just fine. The logs did
show the same messages you report:

PGS_CSC_ECRtoGEO():PGSCSC_W_INVALID_ALTITUDE:35877
invalid altitude - probably indicates bad input data

PGS_CSC_GrazingRay():PGSCSC_W_HIT_EARTH:35897
The line of sight strikes the Earth

But since the file did produce a valid GEO (and subsequent L1B & L2 files), I'll assume that at least these are red-herring notices.

Since it's easy to check, you may want to verify permissions are not a problem.

Sean

By allyn Date 2012-03-30 15:52

Sean & Gwyn,

After taking a break from this due to deadlines for other work, I came back to it today, played around a little, and low and behold, we are in business! My 1st batch went well and the 2nd one seems to be going well too.

By playing around, I mean that I mounted the external harddrive containing the L1A files directly on the linux VA rather than going through a shared folder with the window's host. Not sure if this changed anything. I also retested the processing steps from the gui, rechecked the seadas installation, luts files, and altered my script a little to have it match exactly what the README examples had. (I think this was key.) It sure takes a while for a newbie like me to learn my way around Seadas, but I'm really looking forward to the results!

Thank you again for all of your help. You have been wonderful.

All the best,
Allyn

By gwyn Date 2012-03-30 18:20

Allyn -

Glad it's working! Those particular geolocation error messages are indeed red herrings.

Gwyn

By allyn Date 2012-04-04 11:47

...and my out_ascii script is up and running. I basically stuck to the example.

Quick question:
I saw the discussion about PRODFAIL for l2bin, and well... Can you just confirm for me that, if I have a high PRODFAIL or even 100% PRODFAIL for one or more of my products, it will not make me lose valid pixel data for a different product in the l2gen and out_ascii processes.

Thank you again.
Allyn

By @sean Date 2012-04-04 12:51

Only if you mask those data.
It is a flag that is set by l2gen if a product has a failure condition, but it does not prevent the processing of other products.
However, if set and if you use PRODFAIL in the flaguse parameter of the binner, - no product for the flagged pixels will be binned.

Sean

By allyn Date 2012-05-14 12:05 Edited 2012-05-14 15:35

Hello again,

New question: Some files run perfectly, but often, I get errors in the terminal:

>>

geogen_modis created 'geofilename'.GEO successfully!
Locating pixel/line range ...
Error locating pixel/line range to extract.

<<OR>>

PROCESSING /media/Iomega_/seadas_shared/2006-test/MYD01.A2006286.1130.005.2012039111823.hdf
GENERATING NEW GEOFILE...
Determining required attitude and ephemeris files...
Searching database: /home/ack/Desktop/seadas6.3/log/ancillary_data.db

Input file: /media/Iomega_/seadas_shared/2006-test/MYD01.A2006286.1130.005.2012039111823.hdf
Start time: 2006286113000
End time: 2006286113500

Traceback (most recent call last):
File "/home/ack/Desktop/seadas6.3/run/scripts/modis_GEO.py", line 118, in <module>
m.atteph()
File "/home/ack/Desktop/seadas6.3/run/scripts/modules/modis_GEO_utils.py", line 226, in atteph
self.attfile1 = os.path.basename(get.files['att1'])
KeyError: 'att1'
GENERATING SUBFILES SUB.L1A_LAC and SUB.GEO...

Locating pixel/line range ...
Error locating pixel/line range to extract.

>>

Again, this kind of thing happens maybe one-quarter of the time. Any idea what causes is? And is there any way around it? Thank you again for all of your help.

Best,
Allyn