How to get higher than 1000m/pixel resolution?

[woodbri]

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

...

[woodbri]

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);

[OB.DAAC - SeanBailey]

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.

[woodbri]

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'

[woodbri]

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

[OB.DAAC - SeanBailey]

...I did say it needs more polish :razz:  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

[woodbri]

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.

[woodbri]

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.

[OB SeaDAS - dshea]

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

[woodbri]

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!