davidfrantz / force

Framework for Operational Radiometric Correction for Environmental monitoring
GNU General Public License v3.0
172 stars 50 forks source link

No data pixels appear in Landsat-8/9 processing #191

Closed mhelleis closed 2 years ago

mhelleis commented 2 years ago

When processing Landsat-8 or 9 data using FORCE pixels with nodata values appear in the processed BOA and QAI image (they are not in the L1C data). nodata values are shown yellow in the image:

image

Things we know:

Things we tried turning on/off without changing

While we can simply fill these nodata values using interpolation it would be good to know what the cause is (they don't appear if we use tiling).

mhelleis commented 2 years ago

No insights on this issue?

davidfrantz commented 2 years ago

Hi @mhelleis ,

with "identical", do you mean that they appear on the same spot or is there some randomness involved?

Commonly, there are two main things that may cause this.

1) a gap in the DEM (we can apparently rule this out) 2) a 0-value in the input TOA image. Can you please check if this might be the case?

Cheers, David

mhelleis commented 2 years ago

Hi David,

thanks for the reply!

I just checked again:

  1. There is randomness involved, the nodata pixels appear at different locations for the Landsat-8 and Landsat-9 scene (from the same path/row). Sorry for the confusion.
  2. There is no gap in the DEM data at any of these locations.
  3. The values of the TOA images are not zero at these locations but have reasonable values.
  4. The nodata values appear also in the corresponding QA masks.

Also, we are at the LPS until Friday. If you find the time, we could meet for a coffee and I could show you how the result looks like. If not, no problem.

Best, Max

davidfrantz commented 2 years ago

Hi Max,

yes, we could try to meet up somewhen tomorrow in a coffee break. But we would need to do this spontaneously, I assume.

Cheers, David

mhelleis commented 2 years ago

Sure, you can just write me an email if you have time spontaneously. Or if you can name any break already: Any is fine by me.

Best, Max

mhelleis commented 2 years ago

I've just run another Landsat-8 scene from another area (near Los Angeles, whereas the other scenes were from eastern Australia).

The same artefacts appear also in this scene, but the amount is drastically lower. Most of the the nodata pixels appear here over the ocean, only some appear in areas close to the ocean and pretty much none in the rest of the scene (by looking only at the land are of this scene, I probably would not have notices them). The image below shows a small region of the scene containing the most artefacts (red, pixel value is -9999).

I've checked:

The scene I have processed is LC08_L1TP_041036_20220404_20220412_02_T1.

The output of the log is

/app/force/testdata/LosAngeles/LC08_L1TP_041036_20220404_20220412_02_T1. 
: dc:  66.35%. wc:   6.77%. sc:   0.61%. cc:  15.85%. AOD: 0.0902. # of targets: 242/78.  1 product(s) written. Success! Processing time: 01 mins 42 secs

Not sure if I'm using a wrong parameter somewhere, but I couldn't find anything.

Below is some additional information and the screenshot.

Expand to view used CopDEM30 file names (they cover a larger area than the actual scene)
Copernicus_DSM_COG_10_N32_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N32_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N32_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N33_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N33_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N33_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N33_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N34_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N34_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N34_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N34_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N34_00_W121_00_DEM.tif
Copernicus_DSM_COG_10_N35_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N35_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N35_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N35_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N35_00_W121_00_DEM.tif
Copernicus_DSM_COG_10_N35_00_W122_00_DEM.tif
Copernicus_DSM_COG_10_N36_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N36_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N36_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N36_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N36_00_W121_00_DEM.tif
Copernicus_DSM_COG_10_N36_00_W122_00_DEM.tif
Copernicus_DSM_COG_10_N37_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N37_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N37_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N37_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N37_00_W121_00_DEM.tif
Copernicus_DSM_COG_10_N37_00_W122_00_DEM.tif
Copernicus_DSM_COG_10_N37_00_W123_00_DEM.tif
Copernicus_DSM_COG_10_N38_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N38_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N38_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N38_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N38_00_W121_00_DEM.tif
Copernicus_DSM_COG_10_N38_00_W122_00_DEM.tif
Copernicus_DSM_COG_10_N38_00_W123_00_DEM.tif
Copernicus_DSM_COG_10_N38_00_W124_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W121_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W122_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W123_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W124_00_DEM.tif
Copernicus_DSM_COG_10_N39_00_W125_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W121_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W122_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W123_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W124_00_DEM.tif
Copernicus_DSM_COG_10_N40_00_W125_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W117_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W118_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W119_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W120_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W121_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W122_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W123_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W124_00_DEM.tif
Copernicus_DSM_COG_10_N41_00_W125_00_DEM.tif
Expand to view used parameter file
++PARAM_LEVEL2_START++

#INPUT/OUTPUT DIRECTORIES
#------------------------------------------------------------------------
#The file queue specifies, which images are to be processed. The full path
#to the file needs to be given. Do  not  paste  the content of the file queue
#into the parameter file. The file queue is mandatory for force-level2, but
#may be NULL for force-l2ps.
#Type: full file path
FILE_QUEUE = /app/force/queue_LA.txt
#This is the output directory where the Level 2 data will be stored. Note
#that data will be overwritten/mosaicked if you reprocess images. It is
#safe and recommended to use a single Level 2 data pool for different
#sensors (provided the same grid and projection is used). The higher-level
#programs of FORCE can handle different spatial resolutions (e.g. 30m
#Landsat and 10m Sentinel-2).
#Type: full directory path
DIR_LEVEL2 = /app/force/testdata/LosAngeles/output
#This is the directory where logfiles should be saved.
#Type: full directory path
DIR_LOG = /app/force/logs
#This is a temporary directory that is used to extract compressed images
#for force-level2. Note that images already need to be extracted when using
#force-l2ps directly. The extracted data will be deleted once they were
#processed. If you cancel processing, you may want to delete any left-overs
#in this directory. A file 'cpu-$TIME' is temporarily created in DIR_TEMP.
#This file can be modified to re-adjust the number of CPUs while(!) force-
#level2 is running. Note that the effect is not immediate, as the load is
#only adjusted after one of the running jobs (images) is finished. 
#Type: full directory path
DIR_TEMP = /app/force/tmp

#DIGITAL ELEVATION MODEL
#------------------------------------------------------------------------
#This file specifies the DEM. It is highly recommended to use a DEM. It is
#used for cloud / cloud shadow detection, atmospheric correction and topo-
#graphic correction. The DEM should be a mosaic that should completely cover
#the area you are preprocessing. If there are nodata values in the DEM, the
#Level 2 outputs will have holes, too. It is possible to process without a
#DEM (DEM = NULL). In this case, the surface is assumed flat @ z = 0m.
#Topographic correction cannot be used without a DEM. The quality of atmo-
#spheric correction and cloud /cloud shadow detection will suffer without
#a DEM.
#Type: full file path
FILE_DEM = /app/force/testdata/LosAngeles/dem/copernicus_mosaic.vrt
#Nodata value of the DEM.
#Type: Integer. Valid range: [-32768,32767]
DEM_NODATA = -32767

#DATA CUBES
#------------------------------------------------------------------------
#This indicates whether the images should be reprojected to the target
#coordinate system or if they should stay in their original UTM projection.
#If you want to work with force-higher-level routines, give TRUE.
#Type: Logical. Valid values: {TRUE,FALSE}
DO_REPROJ = FALSE
#This indicates whether the images should be gridded after processing.
#If TRUE, sub-directories for the tiles are generated in DIR_LEVEL2. 
#If FALSE, sub-directories for the original spatial reference systems
#are generated in DIR_LEVEL2. If you want to work with force-higher-level
#routines, give TRUE.
#Type: Logical. Valid values: {TRUE,FALSE}
DO_TILE = FALSE
#This is the tile allow-list. It is an appional file that holds all tiles
#that should be output. Tiles, which are not specified in this file are
#not written to disc. This paremeter is ignored if DO_TILE = FALSE.
#If no tile allow-list should be used, give FILE_TILE = NULL, in which
#case all tiles are output.
#Type: full file path
FILE_TILE = NULL
#This is the tile size (in target units, commonly in meters) of the 
#gridded output. tiles are square; not used if DO_TILE = FALSE. 
#Type: Double. Valid range: ]0,...
TILE_SIZE = NULL
#This is the block size (in target units, commonly in meters) of the 
#image chips. Blocks are stripes, i.e. they are as wide as the tile,
#and as high as specified here; not used if DO_TILE = FALSE or 
#OUTPUT_FORMAT = ENVI. The blocks are the primary processing unit of
#the force-hher-level routines.
#Type: Double. Valid range: ]0,TILE_SIZE]
BLOCK_SIZE = NULL
#This is the spatial resolution of Landsat output; not used if DO_REPROJ 
#= FALSE. Note that the tile and block sizes must be a multiple of the
#pixel resolution.
#Type: Double. Valid range: ]0,...
RESOLUTION_LANDSAT = NULL
#This is the spatial resolution of Sentinel-2 output; not used if DO_REPROJ 
#= FALSE. Note that the tile and block sizes must be a multiple of the
#pixel resolution.
#Type: Double. Valid range: ]0,...
RESOLUTION_SENTINEL2 = NULL
#These are the origin coordinates of the grid system in decimal degree
#(negative values for West/South). The upper left corner of tile 
#X0000_Y0000 represents this point. It is a good choice to use a coord-
#inate that is North-West of your study area – to avoid negative tile
#numbers. Not used if DO_TILE = FALSE.
#Type: Double. Valid range: [-90,90]
#Type: Double. Valid range: [-180,180]
ORIGIN_LON = NULL
ORIGIN_LAT = NULL
#This defines the target coordinate system. If DO_REPROJ = FALSE, the
#projection string can be NULL. The coordinate system must either be
#given as WKT string - or can be a predefined coordinate/grid system.
#If one of the predefined systems are used, TILE_SIZE, BLOCK_SIZE,
#ORIGIN_LAT, and ORIGIN_LON are ignored and internally replaced with
#predefined values. Currently, EQUI7 and GLANCE7 are availble. Both
#are globally defined sets of projections with a corresponding grid 
#system. EQUI7 consists of 7 Equi-Distant, continental projections,
#with a tile size of 100km. GLANCE7 consists of 7 Equal-Area, conti-
#nental projections, with a tile size of 150km. One datacube will be
#generated for each continent.
#Type: Character. Valid values: {,EQUI7,GLANCE7}
PROJECTION = PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+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"],AUTHORITY["EPSG","3857"]]
#This is the resampling appion for the reprojection; you can choose
#between Nearest Neighbor (NN), Bilinear (BL) and Cubic Convolution
#(CC); not used if DO_REPROJ = FALSE.
#Type: Character. Valid values: {NN,BL,CC}
RESAMPLING = BL

#RADIOMETRIC CORRECTION appIONS
#------------------------------------------------------------------------
#This indicates if topographic correction should be performed. If TRUE,
#a DEM need to be given.
#Type: Logical. Valid values: {TRUE,FALSE}
DO_ATMO = TRUE
#This indicates if atmospheric correction should be performed. If TRUE,
#Bottom-of-Atmosphere reflectance is computed. If FALSE, only Top-of-Atmo-
#sphere reflectance is computed.
#Type: Logical. Valid values: {TRUE,FALSE}
DO_TOPO = TRUE
#This indicates if BRDF correction should be performed. If TRUE, output is
#nadir BRDF adjusted reflectance instead of BOA reflectance (the output is
#named BOA nonetheless).
#Type: Logical. Valid values: {TRUE,FALSE}
DO_BRDF = TRUE
#This indicates if adjacency effect correction should be performed.
#Type: Logical. Valid values: {TRUE,FALSE}
ADJACENCY_EFFECT = TRUE
#This indicates if multiple scattering (TRUE) or the single scattering
#approximation (FALSE) should be used in the radiative transfer calculations.
#Type: Logical. Valid values: {TRUE,FALSE}
MULTI_SCATTERING = TRUE

#WATER VAPOR CORRECTION appIONS
#------------------------------------------------------------------------
#This is the directory where the water vapor tables are located. Water
#vapor tables are not required for Sentinel-2, in this case DIR_WVPLUT
#may be NULL. For Landsat, it is recommended to use this functionality.
#As a minimum requirement, DIR_WVPLUT may be NULL and a global value
#for WATER_VAPOR needs to be specified. If a directory is given, 
#WATER_VAPOR is ignored. DIR_WVPLUT must contain water vapor tables.
#The 12 climatology tables must exist at least. They are used if the
#daily tables do not exist or if there is no valid daily value.
#Type: full directory path
DIR_WVPLUT = /app/force/water_vapour_database_2000_2020_monthly
#This specifies a global value for atmospheric water vapor content in
#g cm-2. This parameter can be a dummy value to quickly process an image
#without needing to generate a water vapor database. Note that especially
#Landsat-8 is relatively insensitive to atmospheric water vapor (depending
#on wavelength), and external water vapor is not needed to process
#Sentinel-2. The error in using a dummy value is significant for the TM
#sensors.
#Type: Float. Valid range: [0,15]
WATER_VAPOR = 1

#AEROSOL appICAL DEPTH appIONS
#------------------------------------------------------------------------
#This indicates whether the internal AOD estimation (TRUE) or externally
#generated AOD values should be used (FALSE).
#Type: Logical. Valid values: {TRUE,FALSE}
DO_AOD  = TRUE
#This is the directory where the aerosol appical depth look-up-tables are
#located. They can be used to input external AOD values. It is recom-
#mended to use the internal algorithm only. If a path is given, and 
#DO_ATMO = TRUE, internal AOD estimation is used and external AOD values
#are used as fallback appion.
#Type: full directory path
DIR_AOD  = NULL

#CLOUD DETECTION appIONS
#------------------------------------------------------------------------
#If this parameter is enabled, confident cloud detections will be erased in the 
#reflectance product, i.e. pixels are set to nodata. The cloud flag in the QAI 
#product will still mark these pixels as clouds. Use this appion is disk space 
#is of concern.
#Type: Logical. Valid values: {TRUE,FALSE}
ERASE_CLOUDS = FALSE
#This parameter cancels the processing of images that exceed the given
#threshold. The processing will be canceled after cloud detection.
#Type: Integer. Valid range: ]0,100]
MAX_CLOUD_COVER_FRAME = 100
#This parameter works on a tile basis. It suppresses the output for chips
#(tiled image) that exceed the given threshold.
#Type: Integer. Valid range: ]0,100]
MAX_CLOUD_COVER_TILE  = 100
#Buffer sizes (radius in meters) for cloud, cloud shadow and snow masks.
#Type: Float. Valid range: [0,10000]
CLOUD_BUFFER  = 300
SHADOW_BUFFER = 90
SNOW_BUFFER   = 30
#These are the main thresholds of the Fmask algorithm.
#Type: Float. Valid range: [0,1]
CLOUD_THRESHOLD  = 0.225
SHADOW_THRESHOLD = 0.02

#RESOLUTION MERGING
#------------------------------------------------------------------------
#This parameter defines the method used for improving the spatial reso-
#lution of Sentinel-2’s 20 m bands to 10 m. Pixels flagged as cloud or
#shadow will be skipped. Following methods are available: IMPROPHE uses
#the ImproPhe code in a spectral-only setup; REGRESSION uses a multi-
#parameter regression (results are expected to be best, but processing
#time is significant); STARFM uses a spectral-only setup of the Spatial
#and Temporal Adaptive Reflectance Fusion Model (prediction artifacts
#may occur between land cover boundaries); NONE disables resolution merge;
#in this case, 20m bands are quadrupled.
#Type: Character. Valid values: {IMPROPHE,REGRESSION,STARFM,NONE}
RES_MERGE = IMPROPHE

#CO-REGISTRATION appIONS
#------------------------------------------------------------------------
#This parameter only applies for Sentinel-2 data. This parameter defines
#the path to a directory that contains monthly Landsat NIR base images.
#If given, a co-registration is attempted. If it fails (no tie points),
#the image won't be processed.
#Type: full directory path
DIR_COREG_BASE = NULL
#This parameter defines the nodata values of the coregistration base images.
#Type: Integer. Valid values: [-32768,32767]
COREG_BASE_NODATA = -9999

#MISCELLANEOUS appIONS
#------------------------------------------------------------------------
#This parameter defines if impulse noise should be removed. Ony applies
#to 8bit input data.
#Type: Logical. Valid values: {TRUE,FALSE}
IMPULSE_NOISE = TRUE
#This parameter defines if nodata pixels should be buffered by 1 pixel.
#Type: Logical. Valid values: {TRUE,FALSE}
BUFFER_NODATA = FALSE

#TIER LEVEL
#------------------------------------------------------------------------
#This parameter specifies the acceptable tier level of Landsat Level 1 data.
#For pre-collection data, TIER = 1 will only accept L1T images, TIER = 2
#will also accept L1Gt and L1G images. For collection data, TIER = 1 will
#only accept L1TP images, TIER = 2 will also accept T2 images, TIER = 3
#will additionaly accept RT images.
#Type: Integer. Valid range: [1,3]
TIER = 3

#PARALLEL PROCESSING
#------------------------------------------------------------------------
#Multiprocessing appions (NPROC, DELAY) only apply when using the batch
#utility force-level2. They are not used by the core function force-l2ps.
#------------------------------------------------------------------------
#This module is using hybrid parallelization, i.e. a mix of multiprocessing
#and multithreading. Each input image is one process, each process may use
#multiple threads. In general, it is most efficient to use as much multi-
#processing as possible (i.e. NTHREAD = 1 or 2). However, if you only have
#a small number of images - or if your system does not have enough RAM,
#it is adviced to use less processes and more threads per process. The
#number of processes and threads is given by following parameters.
#Type: Integer. Valid range: [1,...
NPROC = 32
NTHREAD = 2
#This parameter controls whether the individual bands of the Level 1 input
#images are read sequentially or in parallel. Note that we observed two kinds
#of GDAL installation: (1) the JPEG driver reads each band parallely, but 
#separated images in sequence - we recommend to disable PARALLEL_READS in this
#case (for Sentinel-2). (2) The GDAL JPEG drived does not do anything in 
#parallel - use PARALLEL_READ to speed up the work (also use it for Landsat).
#Type: Logical. Valid values: {TRUE,FALSE}
PARALLEL_READS = FALSE
#This parameter sets a delay before starting a new process. This can be help-
#ful to prevent I/O jams when using a lot of processes. The delay is given
#in seconds.
#Type: Integer. Valid range: [0,...
DELAY = 3
#This parameter sets a timeout for unzipping the Level 1 data (only applies when
#images are still in zip/tar.gz format. Only applies for force-level2).
#The timeout is given in seconds.
#Type: Integer. Valid range: [0,...
TIMEOUT_ZIP = 30

#OUTPUT appIONS
#------------------------------------------------------------------------
#Output format, which is either uncompressed flat binary image format aka
#ENVI Standard or GeoTiff. GeoTiff images are compressed with LZW and hori-
#zontal differencing; BigTiff support is enabled; the Tiff is structured 
#with striped blocks according to the TILE_SIZE (X) and BLOCK_SIZE (Y) speci-
#fications. Metadata are written to the ENVI header or directly into the Tiff
#to the FORCE domain. If the size of the metadata exceeds the Tiff's limit,
#an external .aux.xml file is additionally generated.
#Type: Character. Valid values: {ENVI,GTiff}
OUTPUT_FORMAT = GTiff
#Output the cloud/cloud shadow/snow distance output? Note that this is NOT
#the cloud mask (which is sitting in the mandatory QAI product). This pro-
#duct can be used in force-level3; no other higher-level FORCE module is
#using this.
#Type: Logical. Valid values: {TRUE,FALSE}
OUTPUT_DST = FALSE
#Output Aerosol appical Depth map for the green band? No higher-level FORCE
#module is using this.
#Type: Logical. Valid values: {TRUE,FALSE}
OUTPUT_AOD = FALSE
#Output the Water Wapor map? No higher-level FORCE module is using this.
#Type: Logical. Valid values: {TRUE,FALSE}
OUTPUT_WVP = FALSE
#Output the view zenith map? This product can be used in force-level3; no
#other higher-level FORCE module is using this.
#Type: Logical. Valid values: {TRUE,FALSE}
OUTPUT_VZN = FALSE
#Output the  Haze appimzed Transformation output? This product can be
#used in force-level3; no other higher-level FORCE module is using this.
#Type: Logical. Valid values: {TRUE,FALSE}
OUTPUT_HOT = FALSE
#Output overview thumbnails? These are jpegs at reduced spatial resolution,
#which feature an RGB overview + quality information overlayed (pink: cloud,
#red: cirrus, cyan: cloud shadow, yellow: snow, orange: saturated, green:
#subzero reflectance). No higher-level FORCE module is using this.
#Type: Logical. Valid values: {TRUE,FALSE}
OUTPUT_OVV = FALSE

++PARAM_LEVEL2_END++
image

Cheers, Max

davidfrantz commented 2 years ago

Hi @mhelleis,

I just tested with the LA image. This is what happens:

The cirrus band (B9) seems to be the culprit.

In order to convert the DNs to reflectance, you have to apply additive and multiplicative factors.

ref = (meta->cal[b].radd + meta->cal[b].rmul*dn_[b][p]) / sun_[g];

The factors are REFLECTANCE_ADD_BAND_9 and REFLECTANCE_MULT_BAND_9 from the metadata.

5000 is the DN where the reflectance becomes negative.

> -0.100000 + 4999*2.0000E-05
[1] -2e-05
> -0.100000 + 5000*2.0000E-05
[1] 0
> -0.100000 + 5001*2.0000E-05
[1] 2e-05

As we cannot recover from negative reflectance, I masked that out.

I have to think a bit about whether we can relax this for the cirrus band, or if this will have some other unwanted follow-up effects.

I'll keep you posted.

Cheers, David

mhelleis commented 2 years ago

Hi David,

thanks a lot for the info! We were not aware of this issue and - looking at their website - the USGS only mentions this problem over water bodies and doesn't mention it at all for land areas. However, now I found this document, where the same behaviour was observed for the cirrus band over land.

Regarding FORCE: Should we not be able to set these pixels to value zero, as it should be a problem rooted in the calibration values? But if that does cause follow-up problems, one could simply add a pixel to the QA mask so that this behaviour is documented there.

davidfrantz commented 2 years ago

Hi Max,

I have implemented this now (included in develop). 0 reflectance in cirrus band is now allowed. I believe this should not have any negative side effects.

I will close this issue now :)

Thanks, David