Using GPT's WriteImage Operator to output GeoTIFF

[simkins]

Hi Norman,

Thanks for your interest and information. My goal is to go from raw L0 VIIRS data to L2 SST and output this as a grayscale TIFF or GeoTIFF. I have to follow strict pixel sizing instructions (1100x1100m Mercator_1SP projection) so this TIFF can be used by an external script. I've tried this multiple ways and the most recent way has me using l2gen, l2bin, and then l3mapgen and finally creating a grayscale TIFF (I've attached it as a screenshot of it because it's too large to attach by itself)


     1. sh viiirs_L0.sh
     2. sh viirs_RDR.sh

     3. l2gen ifile=SVM01.h5 geofile=GMTCO.h5 ofile=SV_L2_file \
         l2prod1="sst qual_sst" resolution=1100 ctl_pt_incr=1 ctl_pt_incr=1 proc_ocean=1 \
         proc_sst=1 proc_land=0 atmocor=1 maskcloud=1 maskland=1 maskhilt=0 maskstlight=0

     4. l2bin infile=SV_L2_file ofile=SV_L2_binned resolve=2 suite=SST l3bprod=sst \
         prodtype=regional qual_max=4 verbose=1

     5. l3mapgen ifile=SV_L2_binned ofile=ao1.tiff product=sst projection=mercator \
         resolution=1100.1100 interp=area north=36.274 south=19.0 east=-72.735 \
         west=-87.75 oformat=TIFF palfile=/home/james/seadas-7.4/ocssw/run/data/common/palette/default-black0.pal \
         palette_dir=/Applications/seadas-7.4/ocssw/run/data/common/palette/

The way I initially tried to do this was using gpt operators reproject and mosaic to output a grayscale TIFF to my specifications, but that feature is unavailable until v7.5. Now I *think* I've gotten pretty close here. The last question I really had was how to make sure that the final product is 1100x1100m Pixel Size for the specified domain and making sure that the projection is Mercator_1SP. When using l2bin in step 4, what should I set resolve to if I want 1100x1100m resolution given that my only options for that are 1km or 2km?

Thanks!

[OB SeaDAS - knowles]

(NOTE: I edited this entry and changed l2bin suggestion to use qual_max=2)

It is definitely helpful to see the commands you are using.

For l2gen: you do NOT need the resolution field (it is only used currently for MODIS in order to support bands of different resolutions).  The resolution which is actually being used by l2gen for VIIRS is the native resolution (750m).

For l2bin: Since the native resolution is 750m, I would suggest trying resolve=1.  If there is too much peppering of missing data at higher scan angles for what you want, then go for resolve=2.  Using resolve=H, while preserving more of the native resolution feature, will definitely result in peppering of missing data.  The key is that when you bin (at least currently),  each source pixel of data will go into precisely one and only one bin (nearest), and hence the peppering effect for bins of higher resolution than the source.
      
          You probably want to use qual_max=2 (4 is the worst case flagged as invalid sst data)
          
            Note: for l2bin there are actually many resolution options (from usage):  
                resolution (string) (alias=resolve) (default=H) = bin resolution (H = 0.5km,  Q=250m, HQ=100m, HH=50m,
                1 = 1.1km, 2 = 2.3km, 4 = 4.6km,
                9 = 9.2km, 18 = 18.5km, QD = 0.25 degree,
                36 = 36km, HD = 0.5 degree, 1D = 1 degree)

For l3mapgen: Use resolution=1100  (note: you've chosen mercator which is a cylindrical conformal projection and not an equal-area projection so pixels cannot be 1100x1100 except perhaps at the longitude of natural origin).
                The resolution field of l3mapgen sets the longitudinal width.  The specifics of the projection itself will dictate the latitudinal width for each pixel.
                When you use resolution=1100.1100 this is NOT doing what you think, as the resolution field expects a single integer or decimal value.

                Regarding the possible peppering effect caused by setting l2bin:resolve=1,  if it is undesired, you can fill in some or all of this in the mapped image by        using the l3mapgen parameter:               
                fudge=1.25 (or some other factor greater than unity)

                Note: when you specify  projection=mercator   this is actually a shortcut (which is in fact Mercator_1SP):
                Alternatively you could type in proj4 string with parameters, for instance:
                   projection="+proj=merc +lon_0=<central_meridian>"
             (see http://geotiff.maptools.org/proj_list/mercator_1sp.html)

This explanation is a little messy as it covers many programs, topics and parameters, but I think it addresses everything specific to what you are trying to do.

Danny

[OB WebDev - norman]

Here is a Mercator map of your region with length scales overlaid
to help reinforce what Danny said about the projection not being
equal-area.  Only at 30 degrees North in the attached image are the
pixels 1100x1100 meters.  North of that latitude pixels cover smaller
areas while pixels cover larger areas to the south.  Another way of
saying this is that the 1100-kilometer length scale is 1000 pixels
long at 30 North but shorter south of there and longer north of there.

Norman

[OB WebDev - norman]

P.S.
I note that our forum software shrinks the image a bit to fit
in the forum page frame.  Select "View Image" or "Save Image As ..."
in your browser, or go here to see the image at full resolution.
https://oceancolor.gsfc.nasa.gov/oceanforum/attach/69/33569/mercator_map_with_scales.png

Norman

[OB SeaDAS - knowles]

I would like to add that you should really consider running l2bin with qual_max=2 instead of a higher number for qual_max.  Be aware that if you use qual_max = 3 or qual_max=4 this will NOT mask out cloud, straylight, and other potentially bad pixels.  It really depends on whether you are looking for reliable quality science data or some sort of high retrieval imagery.

Danny

[simkins]

Thanks for all of the help Danny and Norman, you've really helped me out a lot.

Regarding the pixel sizing with Mercator Projection...this is something I also realized after creating the imagery and seeing that the pixelsizeX and pixelsizeY that were in my reproject.xml. I know that I'm roughly supposed to follow 1100m x 1100m, but have the exact pixelsizing and coordinates. Is there a way to calculate the approximate resolution (maybe from the mid point since it's not an equal area projection?) from pixel sizes and lat/lon coordinates?

Regarding the qual_max...hopefully the case study I've been working on is a bad pass, because l2bin wouldn't run when qual_max=2.

New Question > how do I properly specify the palette I want to use? I'm trying to use the sst-black0.pal palette but even when l3mapgen can find the palette without returning an error, it doesn't use the palette. Instead, it just uses the default colored palette. Below is my code and attached is a screenshot of the output tiff.

sh viirs_L0.sh $1
sh viirs_RDR.sh $1

l2gen ifile=SVM01.h5 geofile=GMTCO.h5 ofile=SV_L2_file \
l2prod1="sst qual_sst l2_flags chl_oc3" \
ctl_pt_incr=1 ctl_pt_incr=1 proc_ocean=1 \
proc_sst=1 proc_land=0 atmocor=1 maskcloud=1 \
maskland=1 maskhilt=0 maskstlight=0

l2bin infile=SV_L2_file ofile=SV_L2_binned \
resolve=1 \
suite=SST \
l3bprod=sst \
prodtype=regional \
qual_max=3 \
verbose=1

l3mapgen ifile=SV_L2_binned ofile="ao1.sst.$1.tiff" \
product=sst \
projection=mercator \
resolution=1100 \
interp=area \
north=36.274 \
south=19.0 \
east=-72.735 \
west=-87.75 \
oformat=TIFF \
palfile=/home/james/seadas-7.4/ocssw/run/data/common/palette/default-black0.pal

[OB SeaDAS - knowles]

To make things easier to diagnose, can you use this file for testing your l2bin and l3mapgen steps?
https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/V2015045190600.L2_SNPP_SST.nc

I know nothing about your source file, but it is clear from your resultant image that there are clouds in it.  By using qual_max=3 you are producing bogus sst results on pixels which are actually cloud and not water surface.

The palfile you mention "default-black0.pal" is not grey scale.   "black0" means that the ends are black.

Also the correct syntax is (without path and extension)
parfile=default-black0

Try instead:
apply_pal=no

Or if you want a gray-scale palette you can create a gray.pal file

There is no gray.pal file there currently but we'll have one there in the next release.

Danny