I have a question,that is, a HDF file(e.g. SeaWiFS L3 Mapped data), I want to extract some part(e.g, lon=105E,lon=122E,lat=4N,lat=25N) of it and then compute the average value of it, I need the value. Has SeaDAS the function? if it has, how to do it? can you tell me the detail of methods?

Bluesar,

There are several ways to output regional data from
SeaDAS displayed image.

Firstly to get simple ASCII text output with lat/lon/pixel/line/geodata
format:

(1) Use SeaDAS Annotation function to draw a box (or multiple boxes)
cover interested areas with a color.

(2) Use SeaDAS Output :

Output -> Data -> ASCII -> Select Blotchs with certain color -> ASCII text file.

This output is a text file and can be used in spreadsheet program outside SeaDAS
for further analysis.

Long

Secondly for multiple images, user can use
SeaDAS L3 display GUI to select an interested
region.

Use Function ->Arithmetic Band Function
-> Simple Mean

Or use User Defined Band Operation Command
function to calculate average (see Examples in
the function).

Long

Long,

I would like to do an band average using scripts.

But I couldnīt uderstand how to load the bands and make the average at the same script.

I could make a script to load the bands, but how to put it together with the average script I couldnīt make it work.

I've tryed to make a bin image with 4 km of resolution from GAC images, but it seems that I have to use 9 km of resolution becouse of the images problems. And I wouldn't like to loose the 4 km resolution image.

Do you think that using average I can make it?

Thanks

Once the bands are loaded you should be able to simply call the band_avg SeaDAS command, e.g.:

   band_avg,  band_no=[1,2],  /geo

> I've tryed to make a bin image with 4 km of resolution from GAC images, but it seems that I have to use 9 km of resolution becouse of the images problems. And I wouldn't like to loose the 4 km resolution image.

By "images problems" do you mean using 4km resolution results in missing pixels in your images? Can you please describe what you mean?

Also can you describe why you are averaging the bands in case we can offer suggestions as to the best approach?

Mike,

I should I call the band average at the load script?

The problem in my images is that some stripes with no data appear in my 7 days composition. I asked at the forum and someone answered that I need to use the 9 km resolution insted of the original 4 km of the GAC images, becouse of the problem of the missing pixels that were expected. But I didn't want to loose the 4 km resolution.

I am going to use SeaWiFS images, a composition of 7 days, for the period of 2000-2003. My research area is the equatorial coast of Brazil. 

I am averaging the bands becouse there is a lot of clouds in my area, and also becouse I am going to use temporal series (7 days) of  images and other data like wind direction and river flow.

Gustavo.

> I should I call the band average at the load script?

Example of using band_avg:

   filename='A2005205180000.L2'
   load, filename, prod_name=['nLw_412']
   load, filename, prod_name=['nLw_443']
   band_avg, band_no=[1,2], /geo   ; calculates average and puts it into band 3
   display, fbuf=1, band_no=3
   out,'test.png',/display,fbuf=1

> The problem in my images is that some stripes with no data appear in my 7 days composition.

This is a fact of life.. sometimes there just isn't coverage. Instead of increasing your bin size which "creates" data, you could keep your resolution at 4km and use SeaDAS and/or IDL to detect and exclude the missing pixels during your analysis. Your final time series will then correctly display regions that were cloud-covered (or missing for any other reason) for the entire week. If you must fill in these missing pixels for your research, I would suggest coming up with an intelligent interpolation scheme and applying it to the weekly 4km composites to "create" missing data.

Gustavo,

If you can collect some SeaWiFS L2 LAC images (instead GAC to avoid subsampling
in the data), you may have improved resolution.

I used only 3 SeaWiFS L2 MLAC images for a quick run with SeaDAS's l2bin + 2km
. The result looks good. If you use more images and use flags to mask clouds, the
resuled image should be better looking.

Long

Long and Mike,

Long, I've tried to use LAC images but the coverage isn't good even using 7 days compositions. It is so bad to use GAC images? The values aren't so good?

Mike, I've searched on the scripts available for SEADAS but at the msl12 there ins't the option to don't use bad values, and at the l3-space bin there is an option to use user defined product mins, but it isn't working.

Iīve tried to use the bands average but the image that is resulting only have 64 columns and 100 lines, that's very strange. I couldn't find an option to choose the lines and columns of the final image.
Here is the scritp that I wrote with your help.

filename1='S1999360141038.L2'
filename2='S1999360154934.L2'
filename3='S1999361145444.L2'
filename4='S1999361163714.L2'
filename5='teste_media.hdf'

load, filename1, prod_name=['chl_oc4']
load, filename2, prod_name=['chl_oc4']
load, filename3, prod_name=['chl_oc4']
load, filename4, prod_name=['chl_oc4']

band_avg, band_no=[1,2,3,4], /geo, min=0, max=64

out, filename5, /data, band=5, ftype='HDF', win="MAIN", /geo

Thanks to you all

Gustavo Molleri

> It is so bad to use GAC images?

GAC data is ok as long as the 4km resolution is high enough for you.

> Mike, I've searched on the scripts available for SEADAS but at the msl12 there ins't the option to don't use bad values

I don't know what you mean.. msl12 automatically detects bad pixels and sets them to a predefined missing value.

> and at the l3-space bin there is an option to use user defined product mins, but it isn't working.

l2bin and l3bin don't have options for product mins. Where are you seeing this option?

> band_avg, band_no=[1,2,3,4], /geo, min=0, max=64

When you use band_avg all the bands must have the same dimensions. This means you have to project all of the Level 2 files to maps of the same x/y dimension before you use band_ave. To project SeaWiFS Level 2 data with the GUI use either Process->SeaWiFS->Other->bl2map or Utilities->Data Manipulation->Map Projection. To project the data with commands use either bl2map or map_img. You could also process the Level 2 data to Level 3 and work directly with the data bins (no mapping or display involved), or project the Level 3 data with bl3map or smigen.

A question you may want to ask yourself is whether to use Level 2 mapped data or Level 3 data for your research. Level 3 data are equal-area, stored as floats, and include statistical parameters unlike mapped data. Ideally the binning process will preserve both the central tendency (e.g., average) and the variation in the data points that contribute to a bin. So I'm guessing you might want to take a close look at the L3 format. Some documentation to get you started is here.

mike,

I did as you told me so. I mapped the l2 file.

Now the load program can't recognize that the mapped image is a seawifs image.

here is my script:

filename1='/dados2/gustavo/seawifs/GAC/media/S1999360141038.L2_mapped'
filename2='/dados2/gustavo/seawifs/GAC/media/S1999360154934.L2_mapped'
filename3='/dados2/gustavo/seawifs/GAC/media/S1999361145444.L2_mapped'
filename4='/dados2/gustavo/seawifs/GAC/media/S1999361163714.L2_mapped'
filename5='/dados2/gustavo/seawifs/GAC/media/teste_media.hdf'

load, filename1, ftype='SeaWiFS', prod_name=['chl_oc4']
load, filename2, ftype='SeaWiFS', prod_name=['chl_oc4']
load, filename3, ftype='SeaWiFS', prod_name=['chl_oc4']
load, filename4, ftype='SeaWiFS', prod_name=['chl_oc4']

band_avg, band_no=[1,2,3,4], /geo, min=0, max=64

out, filename5, /data, band=5, ftype='HDF', win="MAIN", /geo

here is error message:

[LOAD] - Invalid SeaWiFS file : /dados2/gustavo/seawifs/GAC/media/S1999360141038.L2_mapped
[LOAD] - Invalid SeaWiFS file : /dados2/gustavo/seawifs/GAC/media/S1999360154934.L2_mapped
[LOAD] - Invalid SeaWiFS file : /dados2/gustavo/seawifs/GAC/media/S1999361145444.L2_mapped
[LOAD] - Invalid SeaWiFS file : /dados2/gustavo/seawifs/GAC/media/S1999361163714.L2_mapped
SDP_VAL_BAND: No band available for band_no : 1
SDP_VAL_BAND: No band available for band_no : 2
SDP_VAL_BAND: No band available for band_no : 3
SDP_VAL_BAND: No band available for band_no : 4
% Subscript range values of the form low:high must be >= 0, < size, with low <= high: BAND_INFO (SDS_BAND).
[MAIN]: Error executing the command: band_avg, band_no=[1,2,3,4], /geo, min=0, max=64

I couldn't find were is the error. The S in front of the name of the file doesn't mean the sensor!
Coud be the  ".L2_mapped" that is wrong?

Thanks

Gustavo

Gustavo,

Once you "map" the image, it is no longer ftype='SeaWiFS', but ftype='MAPPED'.

Mark

Mark,

Thanks!
A very small problem, but it makes a huge difference!

Now I am dealing with another problem!

I am trying to load a lot of images at the same time to make an average, but the script that I am using isn't working. It seems that the load program doesn't recognize the others files, just the first. I think that I have to make the load program loop before going to the band average program, but couldn't find a way to do it.

Here is my script:

file=findfile('/dados2/gustavo/seawifs/GAC/2000/2000_01/S*.L2_mapped')
n = n_elements(file)
for i=0,n-1 do begin & $
base=strmid(file(i),33,7) & $

load, file(i), FTYPE="mapped", prod_name=["chl_oc4"] & $

band_avg, band_no=[1,2,3,4,5,6,7,8,9,10,11,12,13], /geo, min=0.01, max=64.0, replace=-1, /exclusive

out_ascii, ascii, band_no=14, vars=[2], /ll_flag, rep_factor=[695], format="(F11.5)"

end

Thanks

Gustavo Molleri

Just put the band_avg and out_ascii commands outside of the loop at the end of the script like so:

file=findfile('/dados2/gustavo/seawifs/GAC/2000/2000_01/S*.L2_mapped')
n = n_elements(file)
for i=0,n-1 do begin & $
  base=strmid(file(i),33,7) & $
  load, file(i), FTYPE="mapped", prod_name=["chl_oc4"] & $
end
band_avg, band_no=[1,2,3,4,5,6,7,8,9,10,11,12,13], /geo, min=0.01, max=64.0, replace=-1, /exclusive
out_ascii, ascii, band_no=14, vars=[2], /ll_flag, rep_factor=[695], format="(F11.5)"

Mike,

I've tried what you said, but the same message appears.

Error message:

IDL Version 6.3 (linux x86 m32). (c) 2006, Research Systems, Inc.

% Embedded IDL: NASA GSFC SeaDAS Development, SeaDAS.
% Embedded IDL: NASA GSFC SeaDAS Development, SeaDAS.
SeaDAS Version 5.0.5 (pid = 22225)
% Attempt to subscript FILE with I is out of range.
[MAIN]: Error executing the command:  load, file(i), FTYPE="mapped", prod_name=["chl_oc4"] & end
MAIN: Error executing command in batch file

I think that the LOAD program doesn't work with this kind  "file(i)" of command!

Thanks Gustavo Molleri

Gustavo,

I need to see your full script again.  The LOAD command should work fine with "file(i)".  Something else is probably wrong.  The "& end" at the end of the LOAD command looks suspect in the error message.  That's why I need to see your full script.

Mark

Mark,

here is my script:

file=findfile('/dados2/gustavo/seawifs/GAC/2000/2000_01/S*.L2_mapped')
n = n_elements(file)
for i=0,n-1 do begin & $
base=strmid(file(i),33,7) & $

ascii='/dados2/gustavo/seawifs/GAC/2000/ascii/chl_oc4/'+base(0)+'.chl_oc4.ascii' 
png='/dados2/gustavo/seawifs/GAC/2000/png/chl_oc4/'+base(0)+'.chl_oc4.png'
mapped='/dados2/gustavo/seawifs/GAC/2000/mapped/chl_oc4/'+base(0)+'chl_oc4.mapped'

load, file(i), FTYPE="mapped", prod_name=["chl_oc4"] & $
end

band_avg, band_no=[1,2,3,4,5,6,7,8,9,10,11,12,13], /geo, min=0.01, max=64.0, replace=-1, /exclusive

out_ascii, ascii, band_no=14, vars=[2], /ll_flag, rep_factor=[695], format="(F11.5)"

out,mapped,/data, band=14, ftype='mapped', /geo

display, band_no=14, smin=0.01, smax=64.00, stype="LOG", /GEO
loadpal, "/dados2/gustavo/scripts/chl_oc4.txt"
load_graph, "/dados2/gustavo/scripts/graph_695x695_1b.flat", fbuf=fbuf, merge=merge, /OVER
out, png

endfor

Thanks

Gustavo Molleri

Gustavo, you need to put "& $" at the end of every line of code in the loop, e.g.:

for i=0,n-1 do begin & $
base=strmid(file(i),33,7) & $

ascii='/dados2/gustavo/seawifs/GAC/2000/ascii/chl_oc4/'+base(0)+'.chl_oc4.ascii' & $
png='/dados2/gustavo/seawifs/GAC/2000/png/chl_oc4/'+base(0)+'.chl_oc4.png' & $
mapped='/dados2/gustavo/seawifs/GAC/2000/mapped/chl_oc4/'+base(0)+'chl_oc4.mapped' & $

load, file(i), FTYPE="mapped", prod_name=["chl_oc4"] & $
end

Mike,

I could make the mapped file! Thanks!

But one thing that I notice is that the mapped image doesn't have navigation information, I couldn't insert a coastline or a grid becouse of this. Is that right?

thanks!

Gustavo Molleri

Gustavo,

Geo info should be retained as SeaDAS mapped file output
even after 'band_avg' calculation.

Long

Long,

I don't know if I understood what you said.
But the images that I am using to make the averages are Seadas Mapped.
The problem occours when I open the mapped images generated from the average.
Those images doesn't have spatial information. Is that normal? I am going to have to map the images again?

Thanks

Gustavo Molleri

Gustavo, I just did a test in a new SeaDAS session:

  load, 'A2002185205500_chlor_a_map.hdf', FTYPE="mapped", prod_name=["Mapped - chlor_a"]
  load, 'A2002185205500_chlor_a_map2.hdf', FTYPE="mapped", prod_name=["Mapped - chlor_a"]
  band_avg, band_no=[1,2], /geo, min=0.01, max=64.0, replace=-1, /exclusive
  out,'mean.hdf',/data, band=3, ftype='mapped', /geo
  load, 'mean.hdf', FTYPE="mapped"
  display,band=4
  coast

...this results in a new SeaDAS mapped file called mean.hdf that does have geolocation information. So can you please try to do the same as above and make sure it works for you.. (your mapped products might be called something other than "Mapped - chlor_a")

Gustavo,

Yes, Mike and I both did tests for you with either SeaWiFS or MODIS mapped files to
create a new band with 'band_avg'. Coastline and other geo-info are still working.

You may try two test files in :

ftp://samoa.gsfc.nasa.gov/seadas/pub/gustavo/S2000123173200_chlor_a_map.hdf
ftp://samoa.gsfc.nasa.gov/seadas/pub/gustavo/S2000132173043_chlor_a_map.hdf

With simple SeaDAS command like:
    band_avg,band_no=[1,2],/geo,min=0,max=5

Also suggest you try to first average only a few files (2-3) to see the results.

Long