how to batch process avhrr data in seadas

By irenealabia Date 2011-12-08 04:02

Hi,

I have downloaded monthly global avhrr 4-km res. data from podacc from 1999-2008 (sample file: 200201.s04m3pfv50-bsst.hdf) and i need to subset the data to my study area (20N; 50N; 140W;-160E), resample them to 9km resolution and output them to readable GMT format for further analysis. I just don't have the idea how to get things started with seadas. Any suggestion will be highly appreciated. Many thanks in advance.

Regards,
irene

By msb65 Date 2011-12-08 14:31

Hi,

I am curious about how to do this as well.

Regarding loading a Pathfinder v5 file into SeaDAS, is it as simple as writing a SeaDAS batch file and utilizing the "load" function? Is there a way to also supply the corresponding flag file and flag threshold so that the SST can be quality controlled during the loading process (similar to loading a MODIS Level-2/SMI SST file)? Or will that have to be done either with band functions in SeaDAS, or afterwards in Matlab, IDL, etc. ?

Thanks!

Mike

By @norman Date 2011-12-08 18:25 Edited 2011-12-12 15:44

Hi Irene,

I can't answer Mike's SeaDAS questions, but the following
work for me in GMT.

# Get the input file from the PODAAC.
ncftpget ftp://podaac-ftp.jpl.nasa.gov/allData/avhrr/L3/pathfinder_v5/monthly/day/04km/2002/01/200201.s04m3pfv50-bsst.hdf

# Extract the BSST data.
hdp dumpsds -n bsst -d -b 200201.s04m3pfv50-bsst.hdf > 200201.s04m3pfv50-bsst.dat

# Convert to GMT .grd format.
xyz2grd 200201.s04m3pfv50-bsst.dat -G200201.s04m3pfv50-bsst.scaled.grd -I0.0439453125 -Rd -F -N0 -ZTLh

# Convert scaled integers to temperatures.
grdmath 200201.s04m3pfv50-bsst.scaled.grd 0.075 MUL 3 SUB = 200201.s04m3pfv50-bsst.grd

# Extract your region of interest and regrid to a 5-minute interval
# (our standard 9-kilometer interval).
grdsample 200201.s04m3pfv50-bsst.grd -G200201.s04m3pfv50-bsst.50N.20N.160E.140W.grd -I5m -R160/220/20/50

# Turn the grid file into an image using a color palette table (sst.cpt).
grdimage 200201.s04m3pfv50-bsst.50N.20N.160E.140W.grd -Csst.cpt -JB190/35/30/40/6 -P -X1 -Y1 -B10g5 \
| gs -sDEVICE=pngalpha -sOutputFile=- -q -r100 -g800x544 - -dGraphicsAlphaBits=4 -dTextAlphaBits=4 \
> 200201.s04m3pfv50-bsst.50N.20N.160E.140W.png

I attach the output image and color palette table I used below.

Norman

Attachment: sst.cpt (13.1k)

By @sean Date 2011-12-08 18:47

Irene, Mike,

To read the Pathfinder AVHRR files, yes, it is as simple as using the 'load' function.

However, as the flags are contained in a separate dataset, there is no way to mask
the data on load. You will have to read both datasets in and create a mask from the
flag data within SeaDAS (or Matlab).

Use the load_graph function to easily apply a mask of flags in the range of 4 to 7
from the sst flag field (in this case loaded band #2) to the displayed sst band (fbuf=1):

load_graph,2,ftype='band', fbuf=1, range=[4,7], color=[5]

Use the mapimg function to reproject the data. You'd need to know the dimensions of your region to down-res the data properly.
For the full globe, Cylindrical projection of dimensions 4096x2048 will yeild approximatly 9km resolution (or half the input approx 4km)

Use the output function to output the remapped, masked data to one of several formats.

See our FAQ for examples of how to automate this process.

Sean

By irenealabia Date 2011-12-09 03:44

Hi Sean,

I started with the first guess sst (bsst) at 4km. do i still need to apply a mask of flags for these? i tried to work around with the sample script you have provided and i'm getting error at the mapimg and saving the output files. this is the sample script i used:

files=findfile('*.hdf')
n=n_elements(files)

FOR i=0,n-1 do begin & $
fname=files(i) & $
base=strmid(fname,0,6) & $
out1=base + '.bin' & $
out2=base + '.asc' & $
out3=base + '.png' & $
print, out1,out2,out3 & $
load, fname, ftype='PFV5', pixel=8192, line=4096, skipbytes=0, slope=0.075, intercept=-3, latlon=[90,-180,-90,180], stype=0 & $
display, band_no=1, fbuf=1, smin=-3, smax=35, STYPE='LIN', /GEO & $
mapimg, band_no=1,/CYLIN, limit=[20,130,50,180], xsize=200, ysize=120 & $
display, band_no=2, fbuf=2, smin=-3, smax=35, STYPE='LIN', /GEO & $
landmask, color=7 & $
loadpal, 'RAINBOW' & $
out, fbuf=2, out1,/data, ftype='FLAT', /GEO & $
out_ascii, fbuf=2, out2, vars=[2], rep_factor=1, titles=['SST'] & $
out, fbuf=2, out3, /display, ftype='PNG', /cbar & $
clear_up & $
ENDFOR

thanks, irene

By irenealabia Date 2011-12-09 04:07

hi norman,

is the ncftpget command included in the gmt package or it has to be separately installed? i am using wget to batch avhrr data via ftp and the primitive way i do a batch download is to create a list of all the wget ftp_links and run it as a sh file.txt. is there a better way to do it?

thanks, irene

By irenealabia Date 2011-12-09 04:30

hi sean,

go through the script again, figured it out. mapimg option should be bands not band_no.

regards, irene

By @norman Date 2011-12-09 13:53

Hi Irene,

There is nothing wrong with the way you download
data using wget. I just tend to use ncftpget for FTP
sites because its default behaviour is better suited
to FTP (at least with our firewall set-up). I usually
use wget for HTTP sites for the same reason.

Norman

By irenealabia Date 2011-12-09 17:43

hi Norman,

oh i see, thanks. back to gmt again, i want to make a script to automate the series of gmt commands on your example, any suggestion on the best way to do it. thanks.

regards, irene

By @norman Date 2011-12-09 18:08 Edited 2011-12-12 15:43

Hi Irene,

There are lots of ways to do this. You may wish
to spend some time reading about scripting on
whichever operating system you use.

If you are on a unix/linux type of system and
using a Cshell, you could do something like
the following.

foreach f (*.s04m3pfv50-bsst.hdf)

hdp dumpsds -n bsst -d -b $f > $f:r.dat

xyz2grd $f:r.dat -G$f:r.scaled.grd -I0.0439453125 -Rd -F -N0 -ZTLh

grdmath $f:r.scaled.grd 0.075 MUL 3 SUB = $f:r.grd

grdsample $f:r.grd -G$f:r.50N.20N.160E.140W.grd -I5m -R160/220/20/50

grdimage $f:r.50N.20N.160E.140W.grd -Csst.cpt -JB190/35/30/40/6 -P -X1 -Y1 -B10g5 | gs -sDEVICE=pngalpha -sOutputFile=- -q -r100 -g800x544 - -dGraphicsAlphaBits=4 -dTextAlphaBits=4 > $f:r.50N.20N.160E.140W.png

end

Regards,
Norman

By WhiteG Date 2011-12-11 14:36

If you are going to use GMT then I'm not sure SeaDAS 6 has much to offer. Others have suggested approaches that do not rely on SeaDAS. Are you aware of the new GHRSST PFV52 daily data? We are planning to use that for a current pilot project, but although the reprocessing was done using (a presumably modified) SeaDAS, the .nc files aren't recognized by SeaDAS 6.3 (but can be used with BEAM, so presumably with SeaDAS 7). For PFV50 daily data, we have used SeaDAS batch mode (with embedded IDL) to load each pair (sst, qual) of files, apply quality criteria using the mband command with a simple script, apply a map projection, and then save a SeaDAS mapped hdf file for further processing with SeaDAS.

By irenealabia Date 2011-12-12 01:37

Hi Norman,

Thanks a lot, it got me started.

Best regards,
Irene

By irenealabia Date 2011-12-12 01:40

Hi,

Thanks for the info. haven't heard of the GHRSST PFV52 as of yet, i'd like to check on it. Are this data available at nodc website?

regards, irene

By irenealabia Date 2011-12-12 07:19

Hi Norman,

I did try to run your code using bash and change some command options in gmt. it did work except when i arrived at the grdimage syntax i was getting an error message: GPL Ghostscript 8.71: Unrecoverable error, exit code 1. what could have been wrong.

grdsample $i.grd -G"$i".npgrd -I15m -R130/180/20/50
grdimage $i.npgrd -Csst.cpt -JM6i -P -X1 -Y1 -B10g5 | gs -sDEVICE=pngalpha -sOutputFile=- -q -r100 -g800x544 - -dGraphics AlphaBits=4 -dTextGraphicsAlphaBits=4 > "$i".png

many thanks,irene

By @norman Date 2011-12-12 15:36 Edited 2011-12-12 15:45

Irene,

The "-dGraphicsAlphaBits=4", which turns on antialiasing during
image rasterization, should not contain any white space. The corresponding
option for text rasterization should be "-dTextAlphaBits=4". (Sorry, that
second option was my typo.)

Norman

By irenealabia Date 2011-12-13 01:23

Hi Norman, bunch of thanks. its working now.

cheers,irene

By irenealabia Date 2011-12-13 05:59

Hi Norman,

Sorry i'm back again with a question. I am just wondering how you were able to generate a nice color palette. i tried to make one using gmt command grd2cpt option with the options below. However, it seemed not doing a mask on pixels with no data. How will i make it do the masking?

grd2cpt $i.fin.grd -Crainbow.cpt -Z > test.cpt

many thanks,irene

By @norman Date 2011-12-13 14:04

Hi Irene,

Since your questions are now primarily concerning GMT, you might
find it useful to subscribe to the GMT mailing list where the GMT gurus
answer questions similar to the way we do for ocean color.

The GMT web site also hosts descriptions of the CPT format including
reference to how to specify foreground, background, and not-a-number
colors.

I personally have not used GMT tools to create color palette tables --
preferring to generate them with perl scripts that I have written instead.

Regards,
Norman