Ocean Color Forum - Not logged in
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.
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. ?
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 \
-sDEVICE=pngalpha -sOutputFile=- -q -r100 -g800x544 - -dGraphicsAlphaBits=4 -dTextAlphaBits=4 \
I attach the output image and color palette table I used below.
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=
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
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:
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=, rep_factor=1, titles=['SST'] & $
out, fbuf=2, out3, /display, ftype='PNG', /cbar & $
clear_up & $
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?
go through the script again, figured it out. mapimg option should be bands not band_no.
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.
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.
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
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
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.
Thanks a lot, it got me started.
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?
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
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.)
Hi Norman, bunch of thanks. its working now.
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
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
reference to how to specify foreground, background, and not-a-number
I personally have not used GMT tools to create color palette tables --
preferring to generate them with perl scripts that I have written instead.
Thanks for the information and links you've provided. You have been very helpful,