I am trying to get higher than 1000m resolution output file and have set resolution=250 in l1mapgen, l1brsgen, the  but this does not seem to change the output resolution. Running version v7.5.
What am I doing wrong?

find . -name \*.par -exec grep -l resolution {} \;
./share/modis/l1brsgen_defaults.par
./share/modis/l1mapgen_defaults.par
./share/modis/l2bin_defaults_GIBS.par
./share/modis/l2bin_defaults_SFREFL.par
./share/modis/l3mapgen_defaults_GIBS.par
./share/modis/l3mapgen_defaults_SFREFL.par
./share/modis/msl12_defaults.par
./share/viirs/l2bin_defaults_GIBS.par
./share/viirs/l2bin_defaults_SFREFL.par
./share/viirs/l3mapgen_defaults_GIBS.par
./share/viirs/l3mapgen_defaults_SFREFL.par

And I tried setting resolution=250 in the *.par files, but its not generating higher resolution files.

Also I noticed the following output messages and I'm wondering if this is related to the resolution issue.

Thanks,
-Stephen Woodbridge

Input file MOD00.P2020054.1900_1.PDS.L1B does not have standard name; must include string "1KM" or "LAC" in order to find "QKM" file.
Input file MOD00.P2020054.1900_1.PDS.L1B is MODIS Aqua Level-1B HDF-EOS product.

found in the output below:

l1mapgen ifile=MOD00.P2020054.1900_1.PDS.L1B ofile=MOD00.P2020054.1900_1.PDS.tif oformat=tiff geofile=MOD00.P2020054.1900_1.PDS.GEO east=-73 west=-105 north=46 south=24 width=2467 resolution=250 stype=1

l1brsgen ifile=MOD00.P2020054.1900_1.PDS.L1B ofile=MOD00.P2020054.1900_1.PDS-alt.tif oformat=tiff geofile=MOD00.P2020054.1900_1.PDS.GEO east=-73 west=-105 north=46 south=24 width=2467 resolution=250 stype=1
Input file MOD00.P2020054.1900_1.PDS.L1B does not have standard name; must include string "1KM" or "LAC" in order to find "QKM" file.
Input file MOD00.P2020054.1900_1.PDS.L1B is MODIS Aqua Level-1B HDF-EOS product.
Loading default parameters from /u/oceandata/share/common/msl12_defaults.par
Loading default parameters from /u/oceandata/share/common/l1brsgen_defaults.par
Loading default parameters for MODISA from /u/oceandata/share/modis/msl12_defaults.par
Loading default parameters for MODISA from /u/oceandata/share/modis/l1brsgen_defaults.par
Loading default sub-sensor parameters for MODISA from /u/oceandata/share/modis/aqua/msl12_defaults.par
Loading parameters for suite OC from /u/oceandata/share/modis/msl12_defaults_OC.par
Loading command line parameters

Loading user parameters for MODISA

...

OK, so maybe I making some progress on this, but It would be nice to get some confirmation or redirection if I'm approaching this wrong. After reading the source code is looks like I have to use _LAC in my filenames so I changed the script to do this and this is my current processing:

modis_L1A.py MOD00.P2020055.2120_1.PDS.L0 -o MOD00.P2020055.2120_1.PDS.L1A_LAC

modis_atteph.py -v MOD00.P2020055.2120_1.PDS.L1A_LAC --refreshDB

modis_GEO.py MOD00.P2020055.2120_1.PDS.L1A_LAC -o MOD00.P2020055.2120_1.PDS.GEO

getanc.py --refreshDB MOD00.P2020055.2120_1.PDS.L1A_LAC

modis_L1B.py MOD00.P2020055.2120_1.PDS.L1A_LAC MOD00.P2020055.2120_1.PDS.GEO -o MOD00.P2020055.2120_1.PDS.L1B_LAC -q MOD00.P2020055.2120_1.PDS.L1B_QKM -k MOD00.P2020055.2120_1.PDS.L1B_HKM

l2gen ifile=MOD00.P2020055.2120_1.PDS.L1B_LAC geofile=MOD00.P2020055.2120_1.PDS.GEO ofile=MOD00.P2020055.2120_1.PDS.L2_LAC_OC l2prod=chlor_a par=MOD00.P2020055.2120_1.PDS.L1A_LAC.anc resolution=250

l2mapgen ifile=MOD00.P2020055.2120_1.PDS.L2_LAC_OC ofile=MOD00.P2020055.2120_1.PDS.L2_OC.tif prod=chlor_a apply_pal stype=0 threshold=0.5 east=-107 west=-136 north=38 south=16 width=10148 outmode=tiff

Where width should be set something like this:

$width = int(($bounds{EAST} - $bounds{WEST}) * floor(111120.0 / $resolution) * cos($bounds{NORTH}/57.2957795) +0.5);

Steve,
Your approach to l2mapgen is reasonable, but there is a better way.   The l1/l2mapgen codes, while still functional, are effectively deprecated. 
A better method is to bin the L2 data at a high resolution, then use l3mapgen to generate a mapped output at whatever resolution you desire.
An added bonus is you can also map to projections other than plate carrée (which is all you get with the older mappers).

We have an initial stab at a replacement processing script using l2bin/l3mapgen.  See $OCSSWROOT/scripts/1mpagen.py  and $OCSSWROOT/scripts/l2mapgen.py ...while not perfect (yet) and in need of some more polish before we advertise it more widely, they should give you an idea of how to string the processors together to have a better mousetrap.

Sean

BTW, the l1/2brsgen codes do not map the data, and have no concept of resolution.  In the default configuration, they just produce a subsampled depiction of the input data.

Hi Sean,

Thank you! I am trying to get l2mapgen.py to work but running into a problem(s). First off, a few minor issues to report.

* parser.add_argument('-scaletype', nargs=1, choices=['linear','log','arctan','DEFAULT'],type=str,default=(["DEFAULT"]), ...
  The type is float and needs to be changed to str otherwise it fails if you use -scaletype log
* the help strings need a little work
  "The argument-list is a set of keyword=value pairs."  this is not ture, it requires -option value
* I'm not totally clear on what I should use for ifile and parfile. I assumed ifile would be the same as ifile got l2gen and parfile would be the *.anc file.

Below is my run and how it blows up.

Thanks,
  -Steve

modis_L1A.py MOD00.P2020058.1700_1.PDS -o MOD00.P2020058.1700_1.PDS.L1A_LAC

getanc.py --refreshDB MOD00.P2020058.1700_1.PDS.L1A_LAC

modis_GEO.py MOD00.P2020058.1700_1.PDS.L1A_LAC -o MOD00.P2020058.1700_1.PDS.GEO

modis_L1B.py MOD00.P2020058.1700_1.PDS.L1A_LAC MOD00.P2020058.1700_1.PDS.GEO -o MOD00.P2020058.1700_1.PDS.L1B_LAC -q MOD00.P2020058.1700_1.PDS.L1B_QKM -k MOD00.P2020058.1700_1.PDS.L1B_HKM

# ls
MOD00.P2020058.1700_1.PDS
MOD00.P2020058.1700_1.PDS.GEO
MOD00.P2020058.1700_1.PDS.L1A_LAC
MOD00.P2020058.1700_1.PDS.L1A_LAC.anc
MOD00.P2020058.1700_1.PDS.L1B_HKM
MOD00.P2020058.1700_1.PDS.L1B_LAC
MOD00.P2020058.1700_1.PDS.L1B_QKM

l2mapgen.py -ifile MOD00.P2020058.1700_1.PDS.L1B_LAC -ofile MOD00.P2020058.1700_1.PDS.L2_K490.tif -product Kd_490 -resolution 250m -oformat tiff -fullrange -projection platecarree -scaletype log -threshold 0.5 -apply_pal -palfile /u/oceandata/share/common/palette/universal_bluered-black0.pal MOD00.P2020058.1700_1.PDS.L1A_LAC.anc
{'projection': ['platecarree'], 'met1': ['/u/oceandata/var/anc/2020/058/N202005812_MET_NCEP_6h.hdf'], 'threshold': [0.5], 'ofile': ['MOD00.P2020058.1700_1.PDS.L2_K490.tif'], 'fullrange': True, 'ifile': ['MOD00.P2020058.1700_1.PDS.L1B_LAC'], 'west': ['DEFAULT'], 'oformat': ['tiff'], 'east': ['DEFAULT'], 'fudge': [1.0], 'met2': ['/u/oceandata/var/anc/2020/058/N202005812_MET_NCEP_6h.hdf'], 'met3': ['/u/oceandata/var/anc/2020/058/N202005812_MET_NCEP_6h.hdf'], 'product': ['Kd_490'], 'central_meridian': ['DEFAULT'], 'north': ['DEFAULT'], 'apply_pal': True, 'gibs': False, 'datamax': ['DEFAULT'], 'quiet': False, 'palfile': ['/u/oceandata/share/common/palette/universal_bluered-black0.pal'], 'sstfile': ['/u/oceandata/var/anc/2020/057/N2020057_SST_OIV2AV_24h.nc'], 'icefile': ['/u/oceandata/var/anc/2020/057/N202005700_SEAICE_NSIDC_24h.hdf'], 'parfile': 'MOD00.P2020058.1700_1.PDS.L1A_LAC.anc', 'datamin': ['DEFAULT'], 'resolution': ['250m'], 'south': ['DEFAULT'], 'scaletype': ['log']}
['l2bin', 'flaguse=ATMFAIL,CLDICE,BOWTIEDEL', 'prodtype=regional', 'ofile=tmp.l3bin', 'infile=MOD00.P2020058.1700_1.PDS.L1B_LAC', 'resolve=HQ', 'l3bprod=Kd_490']
l2bin 5.0.0 (Aug 16 2019 12:47:25)
-E- Can not look up sensor ID for MOD00.P2020058.1700_1.PDS.L1B_LAC.
Input listing file: "" not found.
ncdump: MOD00.P2020058.1700_1.PDS.L1B_LAC: NetCDF: Unknown file format
Traceback (most recent call last):
  File "/u/oceandata/scripts/l2mapgen.py", line 285, in <module>
    if __name__ == "__main__": main()
  File "/u/oceandata/scripts/l2mapgen.py", line 234, in main
    geoExtent = getGeoExtent(dict_args['ifile'][0])
  File "/u/oceandata/scripts/l2mapgen.py", line 58, in getGeoExtent
    geoExtent['center_latitude'] = (geoExtent['northernmost_latitude'] + geoExtent['southernmost_latitude']) / 2.0
KeyError: 'northernmost_latitude'

Related to this, my goal is to generate sst, chloro_a, Kd_490 and True Color. Knowing that I want all of these, can these all be computed in a single pass? and how would I do that. Somewhere I thought I saw that you could put a comma separated list of products you want and I think it puts results as different layers in the output file. What was not sure of was how extract them into separate tif files in the later processing steps. I assume this would be much faster than what I'm currently doing by running l2gen multiple times. I also presume that I have to do something different for truecolor as it does not seem to be a product that I can reference.

-Steve

...I did say it needs more polish   Admittedly, I haven't tried it in  a long time, so I'm not surpised that it's got issues.  I'll see if I can find time to stare at it soon.
As for your l2gen question, yes, you can output more than one product per run, e.g.:  l2prod=sst,chlor_a,Kd_490,rhos_vvv

The rhos_vvv will get you surface reflectance for all visible bands of the sensor.  You can be more specific, as each mission has a set of preferred bands for making true color images.
MODIS would be rhos_469,rhos_555,rhos_645

In the workflow defined in the l1/l2mapgen.py scripts, the L2 files are binned and mapped using l2bin/l3mapgen.
You could bin the products separately, and probably should - especially the rhos for TC output and SST as the masking/flagging would be different than the ocean color products.

The l3mapgen code will use these to produce a TC image:
   use_rgb (boolean) (default=no) = should we use product_rgb to make a
        pseudo-true color image
   product_rgb (string) (default=rhos_670,rhos_555,rhos_412) =
        Three products to use for RGB.  Default is sensor-specific.


Sean

Sean,

I'm slowly starting to figure this stuff out, with your advice.
So ignore my error above, I didn't generate L2 file and was passing an L1B file to l2mapgen.py, so that issue is solved.

If I understand correctly I should be able to generate an L2 file using l2gen once that contains multiple products like you suggested with l2prod=sst,chlor_a,Kd_490,rhos_vvv and then extract each product by calling l2mapgen.py (or the equivalent processing) to generate the various individual files I need. I'll write some equivalent processing to l2mapgen.py using it as an example as it still has a few issues.

Some issues are:
* those mentioned above
* the scale_type parameter is messed up and not getting translated to an appropriate int value, so the command blows up if it is set
* resolution parameter needs some work, if you set it to 250m that gets translated to HQ for l2bin where Q=250m and HQ=100m and HH=50m so the output image it stippled rather than smooth.

Ok, I got it working using the following steps. One question I have is on resolution.

If l2bin resolution=250m, then I get holes in the final output image. The only way I could get a smooth output image was to set l2bin resolution=1km. Am I missing a parameter that would allow me to generate higher resolution images that are smooth? This seems like I could generate a 1km resolution output image, than use gdal_translate to to stretch the resolution to 250m/pixel with bilinear interpolation to get the same results and it would probably be much faster. Any thoughts on this?

modis_L1A.py
getanc.py
modis_GEO.py
modis_L1B.py
l2gen l2prod=sst,chlor_a,Kd_490,rhos_vvv
# then follow what l2magen is doing:
l2bin l3bprod=Kd_490 resolution=1
l3mapgen resolution=250m
l2bin l3bprod=sst resolution=1
l3mapgen resolution=250m
l2bin l3bprod=chlor_a resolution=1
l3mapgen resolution=250m
l2bin l3bprod=rhos_469,rhos_555,rhos_645 resolution=1
l3mapgen resolution=250m

I'm using the default palettes and mask at the moment, but want to not mask out the land for the true color, so I need to fiddle with that and I might want to play with the palettes a little, but so far so good.

Thanks for the help Sean.

Remember that most of the MODIS bands used to calculate the products are 1KM,  so you need to decide if 250m resolution is appropriate.  This it what I would change to get the highest possible resolution from MODIS.

modis_L1A.py
getanc.py
modis_GEO.py
modis_L1B.py
l2gen resolution=250 l2prod=sst,chlor_a,Kd_490,rhos_vvv

The problem with the above is that chl will cause a pixel to fail for SST and rhos.  If you have the time and want as many pixels output, split the l2gen into 2 runs:

l2gen resolution=250 l2prod=chlor_a,Kd_490 ofile=chl.L2.nc
l2gen resolution=250 l2prod=sst,rhos_vvv ofile=sst.L2.nc

now bin the two different L2 files.  The area weighting parameter gets rid of the holes and is a better way to bin.

l2bin l3bprod=chlor_a,Kd_490, resolution=Q area_weighting=1 ofile=chl.L3b.nc
l2bin l3bprod=sst,rhos_vvv resolution=Q area_weighting=1 ofile=sst.L3b.nc

now use l3mapgen to produce your pictures.  You might be able to get away with interp=nearest, but see if you like interp=area better.

l3mapgen resolution=qkm ifile=chl.L3b.nc product=chlor_a
l3mapgen resolution=qkm ifile=chl.L3b.nc product=Kd_490

l3mapgen resolution=qkm ifile=sst.L3b.nc product=sst
l3mapgen resolution=qkm ifile=sst.L3b.nc use_rgb=yes product_rgb=rhos_645,rhos_555,rhos_469

don

Don, Thanks! Sounds like a great plan. I appreciate the additional education on this, as I'm more of a programmer/user of the data than scientist when it comes to the ins/outs of this. Thanks again!

I'm not finding area_weighting as a parameter to any seadas commands. Did you mean average=area which seems to be available in l2bin?

Steve,
The *next* release of SeaDAS will have the area_weighting option in l2bin.  Don thought we'd already put it out.
There is, however, an area weighting available in l3mapgen, which is almost as good

Set interp=area for the l3mapgen calls.  You can squeeze a little more out by also adding fudge=3 .

Sean

Ok, thanks and fudge=3 really helped to eliminate the holes in the image.

I have a question on the scaling, what is the math that is used for the linear and the log scaling?
For example, is it based on the default data_min, data_max for a given product? such that if multiple shots are processed they all get the same scaling?

If so, then I'm assuming that the math for linear scaling is something like this to scale it into 0-255 image palette:
pixel = (value - data_min)/(data_max - data_min)*(255 - 0) + 0

So then is log scaling done like:
pixel = (log(value) - log(data_min))/(log(data_max) - log(data_min))*(255 - 0) + 0

For the Kd_490 band I think data_min=0.01 and data_max=6.4 and the units are m**-1, so using the blue-red palette larger numbers are more red, which is what I'm getting. I can't quite get my head around how to interpret/visualize inverse meters (m**-1). I get particle reflectance, so more particles more reflected, clear water nothing reflected. Is that the right way to think about this?

Yes, your scaling math is appropriate

Kd_490 is a diffuse attenuation, the larger the value the more light is attenuated (per meter) with depth (along the path of light).   So, it's not that there are more particles (and thus there is more reflectance) it's that there is more stuff attenuating light along the path. This means less light penetration.  Since attenuation includes contribution from both scattering and absorption, it could be that increased Kd is due to increased particles concentration, which would result in increased backscattering, and thus increased reflectance.  But it could also be increased absorption, so less backscattering and less reflectance.
That said, the algorithm currently used employs a band ratio of remote sensing reflectance

Sean

Am I correct in noting the scaling ONLY seems to be applied when applying a pallet? ie: apply_pal=1. If I generate a grayscale image with apply_pal=0 like:

l3mapgen interp=area 'projection=+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs' threshold=0.5 ofile=OCL.tif north=42.42578 south=23.89974 west=-90.39094 east=-69.06195 ifile=MOD00.P2020058.1835_1.PDS.chl.L3b.nc oformat=tiff fudge=3.0 product=chlor_a apply_pal=0 quiet=0 resolution=qkm scale_type=linear

I seem to get the same output regardless of setting scale_type= linear or log.
This seems makes sense as you are scaling a potentially larger range of values into the fixed 255 bucket palette. where as the gray scale is outputting float32 value so there is not need to scale that.

Sorry for all the questions, but this has been extremely useful. Thank you very much.

Steve,
Correct, scaling is only applied when a color palette is applied.  As you noticed, when grayscale is output, the TIFF is written as a float, so no scaling is applied.

Sean

More questions on true color generation.

The l3mapgen help page shows and example of using: rhos_670,rhos_555,rhos_412
above says to use: rhos_469,rhos_555,rhos_645  (this generates a buddy brown that overlays the chlor_a, and sst data, but no clouds,water or land)

I assume the numbers are frequencies so lower is more red and higher is more blue.

what are the best to get a reasonable true color image?
What should I use for flags if I want to see clouds, land, water, etc, basically just ignore bad/questionable pixels.
Do I need the apply_pal=1, since the output image is not paletted?
What value do nodata pixels get set to?

I'm currently using:

l2bin prodtype=regional ofile=$name.sst.L3b.nc infile=$name.sst.L2.nc resolution=Q average=area l3bprod=sst,rhos_469,rhos_555,rhos_645

l3mapgen interp=area 'projection=$proj' threshold=0.5 ofile=$name.L2_TC.tif north=$north south=$south west=$west east=$east ifile=$name.sst.L3b.nc oformat=tiff fudge=3.0 use_rgb=1 product_rgb=rhos_469,rhos_555,rhos_645 apply_pal=1 quiet=0 resolution=qkm

I presume, that I'll need to run l2bin once for sst and again for the rhos_vvv set using different flags which I'm ok with. from the l2bin Help page:
flaguse (string) (default=ATMFAIL,LAND,HILT,HISATZEN,STRAYLIGHT,CLDICE,COCCOLITH,LOWLW,CHLFAIL,CHLWARN,NAVWARN,ABSAER,MAXAERITER,ATMWARN,HISOLZEN,NAVFAIL,FILTER) = flags masked [see /SENSOR/l2bin_defaults.par]
and ./share/modis/l2bin_defaults.par has
flaguse=ATMFAIL,LAND,HILT,HISATZEN,STRAYLIGHT,CLDICE,COCCOLITH,LOWLW,CHLWARN,CHLFAIL,NAVWARN,MAXAERITER,ATMWARN,HISOLZEN,NAVFAIL,FILTER,HIGLINT
resolve=4
rowgroup=270
noext=1
qual_max=2

Thanks,
-Steve

Steve,

The numbers refer to the wavelength of the band, related to but not frequency.  So, smaller (shorter) wavelengths tend toward blue, longer red.
469nm is blue, 645nm is red.

SeaDAS comes with some default files for l2bin and l3mapgen that make it easy to generate a true color image.  Simply use the suite SFREFL, e.g.
l2bin suite=SFREFL ifile=<input L2> ofile=<output L3b>
l3mapgen suite=SFREFL ifile=<input L3b> ofile=<output mapped TC>

The suite just defines parameter so you don't have to remember them all:

$ cat l2bin_defaults_SFREFL.par
flaguse=ATMFAIL
resolution=1
$ cat l3mapgen_defaults_SFREFL.par
resolution=1km
interp=area
fudge=3
oformat=png
use_rgb=1
product_rgb=rhos_645,rhos_555,rhos_469


As you can see, your calls are pretty close.  So yes, what you were missing was the flaguse parameter.  The default flaguse parameter excludes a lot of things...like land and clouds...so you get not so true color images

Thanks for all the help, I have everything working great now. You guys are great!
Best regards,