I have finally figured out that l3mapgen has a small problem (or I'm missing something). When I generate a TIFF image with l3mapgen with:
Code: Select all
l3mapgen interp=area projection='EPSG:3857' threshold=0.5 ofile=SST-l3mapgen-out.tif ifile=SNPP_VIIRS.20240619T204200.sst.L3b.nc oformat=tiff fudge=5.0 suite=SST product=sst apply_pal=0 quiet=0 resolution=qkm
Loading default parameters from /u/ocssw//share/viirs/npp/l3mapgen_defaults.par
Loading default parameters from /u/ocssw//share/common/l3mapgen_defaults_SST.par
l3mapgen 2.4.0-T2024.20 (Jun 6 2024 21:11:08)
ifile : SNPP_VIIRS.20240619T204200.sst.L3b.nc
ofile : SST-l3mapgen-out.tif
oformat : tiff
projection : EPSG:3857
resolution : 289.895m
product : sst
qual_prod :
north : 25.974
south : 2.589
east : -100.127
west : -125.586
image size : 9342 x 9778
99% complete
actual data min : -5.000000
actual data max : 32.980000
num filled pixels : 16901924
percent filled pixels : 18.50%
I'm not clear on what values should exist to indicate not actual data values, but it would be convenient to be able to specify something like nodata=-32768 to override the default output. This would make it much easier to work with the gdal tools.
Below is a gdalinfo report of the tif.
Code: Select all
$ gdalinfo -stats -hist SST-l3mapgen-out.tif
Driver: GTiff/GeoTIFF
Files: SST-l3mapgen-out.tif
Size is 9778, 9342
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-13980489.748659253120422,2996145.823402185924351)
Pixel Size = (289.894507274305511,-289.894507274305511)
Metadata:
AREA_OR_POINT=Area
OBPG_version=Unspecified
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-13980489.749, 2996145.823) (125d35'19.95"W, 25d58'34.68"N)
Lower Left (-13980489.749, 287951.336) (125d35'19.95"W, 2d35' 9.00"N)
Upper Right (-11145901.257, 2996145.823) (100d 7'31.20"W, 25d58'34.68"N)
Lower Right (-11145901.257, 287951.336) (100d 7'31.20"W, 2d35' 9.00"N)
Center (-12563195.503, 1642048.580) (112d51'25.58"W, 14d35'25.71"N)
Band 1 Block=9778x9342 Type=Float32, ColorInterp=Gray
Minimum=-32767.000, Maximum=32.980, Mean=-26699.336, StdDev=12734.114
0...10...20...30...40...50...60...70...80...90...100 - done.
256 buckets from -32831.3 to 97.2937:
74444152 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 16901924
Metadata:
STATISTICS_MAXIMUM=32.979999542236
STATISTICS_MEAN=-26699.336323204
STATISTICS_MINIMUM=-32767
STATISTICS_STDDEV=12734.11371322
STATISTICS_VALID_PERCENT=100
Ideally, I want to use gdal_translate to set NODATA=-32768 but since there are other values near that like -32767 and maybe others this doesn't work. I was able to work around this issue with the following to set all values < -5 to -32768 and then asserting the NoDataValue=-32768.
Code: Select all
gdal_calc.py -A SST-l3mapgen-out.tif --outfile=SST-out.tif --calc="where(A<-5,-32768,A)" --NoDataValue=-32768 --overwrite
Now gdalinfo on SST-out.tif report (truncated) is much more accurate WRT the actual data values:
Code: Select all
...
Corner Coordinates:
Upper Left (-13980489.749, 2996145.823) (125d35'19.95"W, 25d58'34.68"N)
Lower Left (-13980489.749, 287951.336) (125d35'19.95"W, 2d35' 9.00"N)
Upper Right (-11145901.257, 2996145.823) (100d 7'31.20"W, 25d58'34.68"N)
Lower Right (-11145901.257, 287951.336) (100d 7'31.20"W, 2d35' 9.00"N)
Center (-12563195.503, 1642048.580) (112d51'25.58"W, 14d35'25.71"N)
Band 1 Block=9778x1 Type=Float32, ColorInterp=Gray
Min=-5.000 Max=32.980
Minimum=-5.000, Maximum=32.980, Mean=4.728, StdDev=9.999
NoData Value=-32768
Metadata:
STATISTICS_MAXIMUM=32.979999542236
STATISTICS_MEAN=4.728431679632
STATISTICS_MINIMUM=-5
STATISTICS_STDDEV=9.9994284060557
STATISTICS_VALID_PERCENT=100