From cc67769decc1702051e0c5a9cce08ba1c69bc5a7 Mon Sep 17 00:00:00 2001 From: Paolo Veglio <paolo.veglio@ssec.wisc.edu> Date: Thu, 7 Sep 2023 19:00:09 +0000 Subject: [PATCH] running pre-commit for the first time --- .flake8 | 3 +- .pre-commit-config.yaml | 38 + mvcm/__init__.py | 2 +- mvcm/c_tools/bilinearInterpSST.c | 13 +- mvcm/c_tools/cithr.c | 14 +- mvcm/c_tools/get_GEOS.c | 1 - mvcm/c_tools/get_geos_times.c | 2 +- mvcm/c_tools/get_granule_times.c | 2 +- mvcm/c_tools/getcoord.c | 70 +- .../include/mask_processing_constants.h | 2 +- mvcm/c_tools/swap_bytes.c | 1 - mvcm/cli_mvcm.py | 170 +-- mvcm/conf.py | 39 +- mvcm/main.py | 509 ++++---- mvcm/preprocess_thresholds.py | 604 ++++++---- mvcm/read_data.py | 546 ++++++--- mvcm/restoral.py | 342 ++++-- mvcm/scene.py | 725 +++++++----- mvcm/spectral_tests.py | 1046 ++++++++++------- mvcm/utility_functions.py | 95 +- mvcm/utils.py | 156 +-- mvcm/write_output.py | 37 +- pyproject.toml | 3 +- tests/test_conf.py | 14 +- tests/test_read_data.py | 131 ++- tests/test_restoral.py | 82 ++ tests/test_spectral_tests.py | 254 ++-- 27 files changed, 3016 insertions(+), 1885 deletions(-) create mode 100644 .pre-commit-config.yaml create mode 100644 tests/test_restoral.py diff --git a/.flake8 b/.flake8 index 6deafc2..8dd399a 100644 --- a/.flake8 +++ b/.flake8 @@ -1,2 +1,3 @@ [flake8] -max-line-length = 120 +max-line-length = 88 +extend-ignore = E203 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..c0c9c98 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,38 @@ +# See https://pre-commit.com for more information +# See https://pre-commit.com/hooks.html for more hooks +exclude: '^$' +fail_fast: false +repos: +# - repo: https://github.com/PyCQA/flake8 +# rev: 6.1.0 +# hooks: +# - id: flake8 +# additional_dependencies: [flake8-docstrings, flake8-debugger, flake8-bugbear, mccabe] +# args: [--max-complexity, "10"] + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v3.2.0 + hooks: + - id: trailing-whitespace + - id: end-of-file-fixer + - id: check-yaml + - id: check-added-large-files + - repo: https://github.com/pycqa/isort + rev: 5.12.0 + hooks: + - id: isort + language_version: python3 + # - repo: https://github.com/pre-commit/mirrors-mypy + # rev: 'v1.5.1' # Use the sha / tag you want to point at + # hooks: + # - id: mypy + # additional_dependencies: + # - types-docutils + # - types-pkg-resources + # - types-PyYAML + # - types-requests + # args: ["--python-version", "3.8", "--ignore-missing-imports"] + - repo: https://github.com/psf/black + rev: 23.7.0 # Replace by any tag/version: https://github.com/psf/black/tags + hooks: + - id: black + language_version: python3 # Should be a command that runs python3.6+ diff --git a/mvcm/__init__.py b/mvcm/__init__.py index df9144c..485f44a 100644 --- a/mvcm/__init__.py +++ b/mvcm/__init__.py @@ -1 +1 @@ -__version__ = '0.1.1' +__version__ = "0.1.1" diff --git a/mvcm/c_tools/bilinearInterpSST.c b/mvcm/c_tools/bilinearInterpSST.c index cea548d..2b26c77 100644 --- a/mvcm/c_tools/bilinearInterpSST.c +++ b/mvcm/c_tools/bilinearInterpSST.c @@ -1,6 +1,6 @@ #include <stdlib.h> #include <stdio.h> -#include <math.h> +#include <math.h> #include "sst.h" @@ -55,8 +55,8 @@ int return(-1); if ( i < 0 || i > (npoints_x -1)) return(-1); - - ldi = lon - i; + + ldi = lon - i; if (ldi >= 0.5) { i1 = i + 1; if ( i == (npoints_x - 1)) @@ -67,7 +67,7 @@ int if ( i == 0) i1 = npoints_x - 1; } - + ldj = lat - (j - 90.0); if (ldj >= 0.5) { @@ -82,9 +82,9 @@ int } p = lon - (i1 + 0.5); - if ( p < -1.0 ) + if ( p < -1.0 ) p = p + 360.0; - if ( p > 1.0) + if ( p > 1.0) p = p - 360.0; ip0 = fabsf( p / dx); ip1 = dx - ip0; @@ -109,4 +109,3 @@ int return (0); } - diff --git a/mvcm/c_tools/cithr.c b/mvcm/c_tools/cithr.c index 93fbd93..5ab53e8 100644 --- a/mvcm/c_tools/cithr.c +++ b/mvcm/c_tools/cithr.c @@ -1,7 +1,7 @@ /********************************************************************** Description: - Float function cithr.c + Float function cithr.c Computes the bi-dimensional linear or quadratic interpolation (lagrange form) of a given value to a table of scan angle and 11 um brightness temperature dependent 11-12 micron brightness temperature @@ -70,8 +70,8 @@ float cithr(int key, float sec_vza, float bt11) .520,.405,.305,.194,.199,.294,.409,.715,1.141,2.591,4.433,7.389,11.381, .491,.391,.319,.238,.257,.367,.467,.729,1.083,2.140,3.866,6.720,10.567, .450,.370,.340,.300,.340,.470,.550,.750,1.000,1.500,3.060,5.770, 9.410 }; - - + + /* Initializations */ lt0 = 0.0; @@ -98,7 +98,7 @@ float cithr(int key, float sec_vza, float bt11) j1 = 0; j2 = 0; jj = 0; - + /* Bounds check. */ u = sec_vza; @@ -202,7 +202,7 @@ float cithr(int key, float sec_vza, float bt11) t0 = ttab[j0]; t1 = ttab[j1]; t2 = ttab[j2]; - + /* lagrange polynomials */ lu0 = (u-u1) * (u-u2) / (u0-u1) / (u0-u2); lu1 = (u-u0) * (u-u2) / (u1-u0) / (u1-u2); @@ -210,12 +210,12 @@ float cithr(int key, float sec_vza, float bt11) lt0 = (t-t1) * (t-t2) / (t0-t1) / (t0-t2); lt1 = (t-t0) * (t-t2) / (t1-t0) / (t1-t2); lt2 = (t-t0) * (t-t1) / (t2-t0) / (t2-t1); - + /* interpolating polynomials for the first dimension */ p0 = tab[i0][j0] * lu0 + tab[i1][j0] * lu1 + tab[i2][j0] * lu2; p1 = tab[i0][j1] * lu0 + tab[i1][j1] * lu1 + tab[i2][j1] * lu2; p2 = tab[i0][j2] * lu0 + tab[i1][j2] * lu1 + tab[i2][j2] * lu2; - + /* interpolating polynomial for second dimension */ p = p0*lt0 + p1*lt1 + p2*lt2; diff --git a/mvcm/c_tools/get_GEOS.c b/mvcm/c_tools/get_GEOS.c index 2e0692a..9f84b86 100644 --- a/mvcm/c_tools/get_GEOS.c +++ b/mvcm/c_tools/get_GEOS.c @@ -203,4 +203,3 @@ void get_GEOS(float *lat, float *lon, int res, char *granule_start_time, char *a // return ( return_code); } - diff --git a/mvcm/c_tools/get_geos_times.c b/mvcm/c_tools/get_geos_times.c index 447fc5b..17a107e 100644 --- a/mvcm/c_tools/get_geos_times.c +++ b/mvcm/c_tools/get_geos_times.c @@ -48,7 +48,7 @@ int get_geos_times(char s1[64], char s2[64], int *geoshr1, int *geoshr2) /******************************************************************************/ // Get GEOS start hours. - + // printf("get time 1: %s\n", s1); // printf("get time 2: %s\n", s2); diff --git a/mvcm/c_tools/get_granule_times.c b/mvcm/c_tools/get_granule_times.c index ef97765..a2a74d4 100644 --- a/mvcm/c_tools/get_granule_times.c +++ b/mvcm/c_tools/get_granule_times.c @@ -10,7 +10,7 @@ Input arguments: Output arguments: targhr decimal granule start hour - (target hour for time interplation of GDAS + (target hour for time interplation of GDAS variables) Function output: diff --git a/mvcm/c_tools/getcoord.c b/mvcm/c_tools/getcoord.c index 8b5872b..273ea85 100644 --- a/mvcm/c_tools/getcoord.c +++ b/mvcm/c_tools/getcoord.c @@ -1,16 +1,16 @@ -/* +/* !C ********************************************************************** !Description: - GIHLS211 transforms latitude, longitude coordinates to line, sample + GIHLS211 transforms latitude, longitude coordinates to line, sample Renamed GIHLS211 to GETCOORD !! - coordinates for an image in the Goode's Interrupted Homolosine - projection. This routine was complied and run using the C compiler on + coordinates for an image in the Goode's Interrupted Homolosine + projection. This routine was complied and run using the C compiler on SunOS 4.2 (UNIX). Results were accurate at the time of testing, but they are by no means guaranteed! -!Input Parameters: +!Input Parameters: lt latitude (double precision) -ln longitude (double precision) +ln longitude (double precision) !Output Parameters: lne line number of 1km mask for given latitude/longitude @@ -50,7 +50,7 @@ smp element number of 1km maks for given latitude/longitude extern void getcoord_(double *lt, double *ln, double *lne, double *smp); /*void getcoord_(double *lt, double *ln, double *lne, double *smp);*/ void bad_input_parms(); -void goode_init(double r); +void goode_init(double r); int goode_forward(double lon, double lat, double* x, double* y); /* void ptitle(char* A); @@ -85,7 +85,7 @@ static double feast[12]; /* False easting, one for each region */ /* Transformation routine -------------------------------*/ -void getcoord_(double *lt, double *ln, double *lne, double *smp) +void getcoord_(double *lt, double *ln, double *lne, double *smp) { float pixsiz; double x,y; @@ -109,7 +109,7 @@ pixsiz = 1000.0F; /* Report parameters and image geometric characteristics to the user -----------------------------------------------------------------*/ -/*printf("Converting latitude, longitude to line, sample coordinates\n\n"); +/*printf("Converting latitude, longitude to line, sample coordinates\n\n"); printf("Pixel size is %f km\n\n", pixsiz/1000.0); */ /* Modified by JC Guu 08/26/96 */ @@ -126,7 +126,7 @@ if (fabs(pixsiz-1000.0F) < 0.000001F) { nl = 17347; ns = 40031; */ - } + } else { ul_x = -20011500.0; ul_y = 8669500.0; @@ -134,7 +134,7 @@ else { nl = 2168; ns = 5004; */ - } + } /*printf("Image size is %d lines by %d samples with an upper left\n",nl,ns); printf("corner of UL_X = %lf and UL_Y = %lf meters.\n", ul_x, ul_y); */ @@ -165,7 +165,7 @@ goode_init(6370997.0); line = (ul_y - y) / pixsiz + 1.0; samp = (x - ul_x) / pixsiz + 1.0; -/* printf("%12.6f %12.6f\n",lat, lon); +/* printf("%12.6f %12.6f\n",lat, lon); printf("%12.6f %12.6f\n",line, samp); */ *lne = line; @@ -174,7 +174,7 @@ goode_init(6370997.0); *ln = lon; return; -} +} @@ -212,7 +212,7 @@ printf(" pixsize in km = 1 or 8\n\n");*/ /* is changed to the modern form and the */ /* return type is given as void. */ -void goode_init(double r) +void goode_init(double r) { /* Modified by JC Guu 08/30/96 */ @@ -224,9 +224,9 @@ void goode_init(double r) !Description: Initialize the Goode`s Homolosine projection. Place parameters in static storage for common use. !Input Parameters: r Radius of the earth (sphere) -!Output Parameters: lon_center[12] Central meridians, one for +!Output Parameters: lon_center[12] Central meridians, one for each region - feast[12]; False easting, one for each + feast[12]; False easting, one for each region !Revision History: !Team-unique Header: @@ -270,7 +270,7 @@ feast[11] = R * 2.44346095279; /* Report parameters to the user -----------------------------*/ -/*ptitle("Goode`s Homolosine Equal Area"); +/*ptitle("Goode`s Homolosine Equal Area"); radius(r): */ /* Modified by JC Guu 08/26/96 */ @@ -294,12 +294,12 @@ int goode_forward(double lon, double lat, double* x, double* y) /* !C -!Description:Goode`s Homolosine forward equations--mapping +!Description:Goode`s Homolosine forward equations--mapping lat,long to x,y !Input Parameters: lon Longitude lat Latitude !Output Parameters: x X projection coordinate - y Y projection coordinate + y Y projection coordinate !Revision History: !Team-unique Header: !END @@ -366,7 +366,7 @@ else delta_theta = -(theta + sin(theta) - constant) / (1.0 + cos(theta)); theta += delta_theta; if (fabs(delta_theta) < EPSLN) break; - if (i >= 30) + if (i >= 30) {giherror("Iteration failed to converge","Goode-forward");return(ERROR);} } theta /= 2.0; @@ -389,7 +389,7 @@ return(OK); /* the K&R style to the modern form. */ /* -void ptitle(char* A) +void ptitle(char* A) */ /* Modified by JC Guu 08/30/96 */ @@ -414,7 +414,7 @@ void ptitle(char* A) /* -void radius(double A) +void radius(double A) */ /* Modified by JC Guu 08/30/96 */ @@ -448,7 +448,7 @@ void radius(double A) /* is changed to the modern form and the */ /* return type is given as void. */ -void giherror(char* what, char* where) +void giherror(char* what, char* where) /* Modified by JC Guu 08/30/96 */ /* The required prolog components are added */ @@ -490,11 +490,11 @@ void getcoord_sincos(double val, double* sin_val, double* cos_val) /* !C -!Description: Function to calculate the sine and cosine in - one call. Some computer systems have implemented +!Description: Function to calculate the sine and cosine in + one call. Some computer systems have implemented this function, resulting in a faster implementation - than calling each function separately. It is - provided here for those computer systems which + than calling each function separately. It is + provided here for those computer systems which don`t implement this function. !Input Parameters: val the value to be operated on @@ -518,7 +518,7 @@ void getcoord_sincos(double val, double* sin_val, double* cos_val) /* is changed to the modern form and the */ /* return type is given as int. */ -int getcoord_sign(double x) +int getcoord_sign(double x) { if (x < 0.0) return(-1); else return(1);} /* Modified by JC Guu 08/30/96 */ @@ -550,7 +550,7 @@ double getcoord_adjust_lon(double x) /* !C !Description: Function to adjust longitude to -180 to 180 -!Input Parameters: x +!Input Parameters: x !Output Parameters: !Revision History: !Team-unique Header: @@ -569,7 +569,7 @@ double getcoord_adjust_lon(double x) /* from the K&R style to the modern form for */ /* the following functions. */ -double getcoord_e0fn(double x) +double getcoord_e0fn(double x) /* Modified by JC Guu 08/30/96 */ /* The required prolog components are added */ @@ -579,7 +579,7 @@ double getcoord_e0fn(double x) !C !Description: Function to compute e0 !Input Parameters: x -!Output Parameters: +!Output Parameters: !Revision History: !Team-unique Header: !END @@ -599,7 +599,7 @@ double getcoord_e1fn(double x) !C !Description: Function to compute e1 !Input Parameters: x -!Output Parameters: +!Output Parameters: !Revision History: !Team-unique Header: !END @@ -607,9 +607,9 @@ double getcoord_e1fn(double x) {return(0.375*x*(1.0+0.25*x*(1.0+0.46875*x)));} - -double getcoord_e2fn(double x) + +double getcoord_e2fn(double x) /* Modified by JC Guu 08/30/96 */ /* The required prolog components are added */ @@ -629,7 +629,7 @@ double getcoord_e2fn(double x) -double getcoord_mlfn(double e0,double e1,double e2,double phi) +double getcoord_mlfn(double e0,double e1,double e2,double phi) /* Modified by JC Guu 08/30/96 */ /* The required prolog components are added */ diff --git a/mvcm/c_tools/include/mask_processing_constants.h b/mvcm/c_tools/include/mask_processing_constants.h index bc1b6da..c4b349c 100644 --- a/mvcm/c_tools/include/mask_processing_constants.h +++ b/mvcm/c_tools/include/mask_processing_constants.h @@ -7,7 +7,7 @@ Revision History: 09/2012 R. Frey Original version 11/2012 R. Frey Changed 'badsza_data' to 'badgeo_data'. - 08/2013 R. Frey Eliminated 'max_vza' and + 08/2013 R. Frey Eliminated 'max_vza' and added 'MODIS_max_vza' and 'VIIRS_max_vza'. 08/2015 R. Frey Added "min" and "max" macros. diff --git a/mvcm/c_tools/swap_bytes.c b/mvcm/c_tools/swap_bytes.c index 47fd705..e000555 100644 --- a/mvcm/c_tools/swap_bytes.c +++ b/mvcm/c_tools/swap_bytes.c @@ -123,4 +123,3 @@ fbyte4_(void *buffer, int *num) cbuf += 4; } } - diff --git a/mvcm/cli_mvcm.py b/mvcm/cli_mvcm.py index 123e74a..9965f13 100644 --- a/mvcm/cli_mvcm.py +++ b/mvcm/cli_mvcm.py @@ -1,84 +1,88 @@ import argparse -from glob import glob +# from glob import glob from pkg_resources import get_distribution from mvcm.main import main -_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires' -_fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0] -_fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1312.001.*.uwssec.nc')[0] -_fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0] -_fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1312.001.*.uwssec.nc')[0] +_datapath = "/ships19/hercules/pveglio/mvcm_viirs_hires" +# _fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0] +# _fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1312.001.*.uwssec.nc')[0] +# _fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0] +# _fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1312.001.*.uwssec.nc')[0] # thresholds: -_threshold_file = '/home/pveglio/mvcm/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml' +_threshold_file = "/home/pveglio/mvcm/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml" # ancillary files: -_geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4' -_geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4' -_geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4' -_geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4' -_geos_constants = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4' -_ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf' -_sst_file = 'oisst.20220622' -_eco_file = 'goge1_2_img.v1' +_geos_atm_1 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4" +_geos_atm_2 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4" +_geos_land = "GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4" +_geos_ocean = "GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4" +_geos_constants = "GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4" +_ndvi_file = "NDVI.FM.c004.v2.0.WS.00-04.177.hdf" +_sst_file = "oisst.20220622" +_eco_file = "goge1_2_img.v1" -_VERSION = get_distribution('mvcm').version +_VERSION = get_distribution("mvcm").version def mvcm(): + parser = argparse.ArgumentParser(prog="MVCM", description="") - parser = argparse.ArgumentParser(prog='MVCM', description='') - - parser.add_argument('--satellite', - help='satellite name, not case-sensitive. Available options: [snpp, ]', - choices=['snpp', ]) - parser.add_argument('--sensor', - help='sensor name, not case-sensitive. Available options: [viirs, ]', - choices=['viirs', ]) - parser.add_argument('--path', - help='path where the data is stored') - parser.add_argument('--l1b', - help='level 1b file name') - parser.add_argument('--geolocation', - help='geolocation file name') - parser.add_argument('--hires-l1b', - help='VIIRS IMG02 file name') - parser.add_argument('--hires-geo', - help='VIIRS IMG03 file name') - parser.add_argument('-t', '--threshold', - help='thresholds file name') - parser.add_argument('--atmos-1', - help='file name of the first GEOS-IT file for atmospheric parameters') - parser.add_argument('--atmos-2', - help='file name of the second GEOS-IT file for atmospheric parameters') - parser.add_argument('--land', - help='file name of GEOS-IT land parameters') - parser.add_argument('--ocean', - help='file name of GEOS-IT ocean parameters') - parser.add_argument('--constants', - help='file name of GEOS-IT constants') - parser.add_argument('--ndvi', - help='NDVI file name') - parser.add_argument('--sst', - help='Sea surface temperature file name') - parser.add_argument('--eco', - help='Ecosystem file') - parser.add_argument('-o', '--out', - help='output file name') - parser.add_argument('-V', '--version', action='version', version=_VERSION, - help='print version and exit') - parser.add_argument('-v', '--verbose', action='store_true', - help='print verbose information') + parser.add_argument( + "--satellite", + help="satellite name, not case-sensitive. Available options: [snpp, ]", + choices=[ + "snpp", + ], + ) + parser.add_argument( + "--sensor", + help="sensor name, not case-sensitive. Available options: [viirs, ]", + choices=[ + "viirs", + ], + ) + parser.add_argument("--path", help="path where the data is stored") + parser.add_argument("--l1b", help="level 1b file name") + parser.add_argument("--geolocation", help="geolocation file name") + parser.add_argument("--hires-l1b", help="VIIRS IMG02 file name") + parser.add_argument("--hires-geo", help="VIIRS IMG03 file name") + parser.add_argument("-t", "--threshold", help="thresholds file name") + parser.add_argument( + "--atmos-1", + help="file name of the first GEOS-IT file for atmospheric parameters", + ) + parser.add_argument( + "--atmos-2", + help="file name of the second GEOS-IT file for atmospheric parameters", + ) + parser.add_argument("--land", help="file name of GEOS-IT land parameters") + parser.add_argument("--ocean", help="file name of GEOS-IT ocean parameters") + parser.add_argument("--constants", help="file name of GEOS-IT constants") + parser.add_argument("--ndvi", help="NDVI file name") + parser.add_argument("--sst", help="Sea surface temperature file name") + parser.add_argument("--eco", help="Ecosystem file") + parser.add_argument("-o", "--out", help="output file name") + parser.add_argument( + "-V", + "--version", + action="version", + version=_VERSION, + help="print version and exit", + ) + parser.add_argument( + "-v", "--verbose", action="store_true", help="print verbose information" + ) args = parser.parse_args() - satellite = args.satellite or 'snpp' - sensor = args.sensor or 'viirs' + satellite = args.satellite or "snpp" + sensor = args.sensor or "viirs" data_path = args.path or _datapath - mod02 = args.l1b or _fname_mod02 - mod03 = args.geolocation or _fname_mod03 + mod02 = args.l1b or "" # _fname_mod02 + mod03 = args.geolocation or "" # _fname_mod03 img02 = args.hires_l1b or None img03 = args.hires_geo or None threshold_file = args.threshold or _threshold_file @@ -90,26 +94,26 @@ def mvcm(): ndvi_file = args.ndvi or _ndvi_file sst_file = args.sst or _sst_file eco_file = args.eco or _eco_file - out_file = args.out or 'test_out.nc' + out_file = args.out or "test_out.nc" verbose = args.verbose or False - main(satellite=satellite, - sensor=sensor, - data_path=data_path, - mod02=mod02, - mod03=mod03, - img02=img02, - img03=img03, - threshold_file=threshold_file, - geos_atm_1=geos_atm_1, - geos_atm_2=geos_atm_2, - geos_land=geos_land, - geos_ocean=geos_ocean, - geos_constants=constants, - ndvi_file=ndvi_file, - sst_file=sst_file, - eco_file=eco_file, - out_file=out_file, - verbose=verbose) - - + main( + satellite=satellite, + sensor=sensor, + data_path=data_path, + mod02=mod02, + mod03=mod03, + img02=img02, + img03=img03, + threshold_file=threshold_file, + geos_atm_1=geos_atm_1, + geos_atm_2=geos_atm_2, + geos_land=geos_land, + geos_ocean=geos_ocean, + geos_constants=constants, + ndvi_file=ndvi_file, + sst_file=sst_file, + eco_file=eco_file, + out_file=out_file, + verbose=verbose, + ) diff --git a/mvcm/conf.py b/mvcm/conf.py index 98b826c..6b5969a 100644 --- a/mvcm/conf.py +++ b/mvcm/conf.py @@ -1,8 +1,7 @@ import numpy as np -def conf_test_new(rad: np.ndarray, - thr: np.ndarray) -> np.ndarray: +def conf_test_new(rad: np.ndarray, thr: np.ndarray) -> np.ndarray: """Assuming a linear function between min and max confidence level, the plot below shows how the confidence (y axis) is computed as function of radiance (x axis). This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be @@ -46,7 +45,7 @@ def conf_test_new(rad: np.ndarray, flipped[hicut < locut] = 1 # Rad between alpha and beta - range_ = 2. * (beta - alpha) + range_ = 2.0 * (beta - alpha) s1 = (rad - alpha) / range_ idx = np.nonzero((rad <= beta) & (flipped == 0)) confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx]) @@ -54,7 +53,7 @@ def conf_test_new(rad: np.ndarray, confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx]) # Rad between beta and gamma - range_ = 2. * (beta - gamma) + range_ = 2.0 * (beta - gamma) s1 = (rad - gamma) / range_ idx = np.nonzero((rad > beta) & (flipped == 0)) confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx]) @@ -73,7 +72,7 @@ def conf_test_new(rad: np.ndarray, def conf_test(rad, thr): - ''' + """ Assuming a linear function between min and max confidence level, the plot below shows how the confidence (y axis) is computed as function of radiance (x axis). This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be @@ -91,7 +90,7 @@ def conf_test(rad, thr): e 0________/ | alpha --------- radiance ----------> - ''' + """ radshape = rad.shape rad = rad.reshape(np.prod(radshape)) @@ -117,7 +116,7 @@ def conf_test(rad, thr): flipped[hicut < locut] = 1 # Rad between alpha and beta - range_ = 2. * (beta - alpha) + range_ = 2.0 * (beta - alpha) s1 = (rad - alpha) / range_ idx = np.nonzero((rad <= beta) & (flipped == 0)) confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx]) @@ -125,7 +124,7 @@ def conf_test(rad, thr): confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx]) # Rad between beta and gamma - range_ = 2. * (beta - gamma) + range_ = 2.0 * (beta - gamma) s1 = (rad - gamma) / range_ idx = np.nonzero((rad > beta) & (flipped == 0)) confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx]) @@ -189,23 +188,23 @@ def conf_test_dble(rad, coeffs): # Inner region fails test # Value is within range of lower set of limits - range_ = 2. * (beta1 - alpha1) + range_ = 2.0 * (beta1 - alpha1) s1 = (rad - alpha1) / range_ idx = np.nonzero((rad <= alpha1) & (rad >= beta1) & (alpha1 - gamma1 > 0)) confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx]) - range_ = 2. * (beta1 - gamma1) + range_ = 2.0 * (beta1 - gamma1) s1 = (rad - gamma1) / range_ idx = np.nonzero((rad >= gamma1) & (rad < beta1) & (alpha1 - gamma1 > 0)) confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx]) # Value is within range of upper set of limits - range_ = 2. * (beta2 - alpha2) + range_ = 2.0 * (beta2 - alpha2) s1 = (rad - alpha2) / range_ idx = np.nonzero((rad > alpha1) & (rad <= beta2) & (alpha1 - gamma1 > 0)) confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx]) - range_ = 2. * (beta2 - gamma2) + range_ = 2.0 * (beta2 - gamma2) s1 = (rad - gamma2) / range_ idx = np.nonzero((rad > alpha1) & (rad > beta2) & (alpha1 - gamma1 > 0)) confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx]) @@ -223,23 +222,31 @@ def conf_test_dble(rad, coeffs): # Value is within range of lower set of limits range_ = 2 * (beta1 - alpha1) s1 = (rad - alpha1) / range_ - idx = np.nonzero((rad > alpha1) & (rad <= gamma1) & (rad <= beta1) & (alpha1 - gamma1 <= 0)) + idx = np.nonzero( + (rad > alpha1) & (rad <= gamma1) & (rad <= beta1) & (alpha1 - gamma1 <= 0) + ) confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx]) range_ = 2 * (beta1 - gamma1) s1 = (rad - gamma1) / range_ - idx = np.nonzero((rad > alpha1) & (rad <= gamma1) & (rad > beta1) & (alpha1 - gamma1 <= 0)) + idx = np.nonzero( + (rad > alpha1) & (rad <= gamma1) & (rad > beta1) & (alpha1 - gamma1 <= 0) + ) confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx]) # Values is within range of upper set of limits range_ = 2 * (beta2 - alpha2) s1 = (rad - alpha2) / range_ - idx = np.nonzero((rad > gamma2) & (rad < alpha2) & (rad >= beta2) & (alpha1 - gamma1 <= 0)) + idx = np.nonzero( + (rad > gamma2) & (rad < alpha2) & (rad >= beta2) & (alpha1 - gamma1 <= 0) + ) confidence[idx] = 1.0 - coeff[idx] * np.power(s1[idx], power[idx]) range_ = 2 * (beta2 - gamma2) s1 = (rad - gamma2) / range_ - idx = np.nonzero((rad > gamma2) & (rad < alpha2) & (rad < beta2) & (alpha1 - gamma1 <= 0)) + idx = np.nonzero( + (rad > gamma2) & (rad < alpha2) & (rad < beta2) & (alpha1 - gamma1 <= 0) + ) confidence[idx] = coeff[idx] * np.power(s1[idx], power[idx]) confidence[(alpha1 - gamma1 <= 0) & ((rad > gamma1) | (rad < gamma2))] = 0 diff --git a/mvcm/main.py b/mvcm/main.py index 7ffc5ac..cb56848 100644 --- a/mvcm/main.py +++ b/mvcm/main.py @@ -4,43 +4,41 @@ except ImportError: import ruamel_yaml as yml import functools -import time import logging +import time +from copy import copy +from glob import glob import numpy as np import xarray as xr -from glob import glob -from copy import copy - import mvcm.read_data as rd -import mvcm.spectral_tests as tst import mvcm.scene as scn +import mvcm.spectral_tests as tst import mvcm.utility_functions as utils - from mvcm.restoral import Restoral # #################################################################### # # TEST CASE # data: -_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires' +_datapath = "/ships19/hercules/pveglio/mvcm_viirs_hires" # _fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0] # _fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1312.001.*.uwssec.nc')[0] # _fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0] # _fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1312.001.*.uwssec.nc')[0] # thresholds: -_threshold_file = '/home/pveglio/mvcm/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml' +_threshold_file = "/home/pveglio/mvcm/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml" # ancillary files: -_geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4' -_geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4' -_geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4' -_geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4' -_geos_constants = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4' -_ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf' -_sst_file = 'oisst.20220622' -_eco_file = 'goge1_2_img.v1' +_geos_atm_1 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4" +_geos_atm_2 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4" +_geos_land = "GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4" +_geos_ocean = "GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4" +_geos_constants = "GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4" +_ndvi_file = "NDVI.FM.c004.v2.0.WS.00-04.177.hdf" +_sst_file = "oisst.20220622" +_eco_file = "goge1_2_img.v1" # #################################################################### # @@ -52,30 +50,33 @@ def timer(func): value = func(*args, **kwargs) end_time = time.time() run_time = end_time - start_time - print(f'Finished {func.__name__!r} in {run_time:.2f} secs\n') + print(f"Finished {func.__name__!r} in {run_time:.2f} secs\n") return value + return wrapper_timer @timer -def main(satellite: str = 'snpp', - sensor: str = 'viirs', - data_path: str = _datapath, - mod02: str = '', # _fname_mod02, - mod03: str = '', # _fname_mod03, - img02: str = '', # _fname_img02, - img03: str = '', # _fname_img03, - threshold_file: str = _threshold_file, - geos_atm_1: str = _geos_atm_1, - geos_atm_2: str = _geos_atm_2, - geos_land: str = _geos_land, - geos_ocean: str = _geos_ocean, - geos_constants: str = _geos_constants, - ndvi_file: str = _ndvi_file, - sst_file: str = _sst_file, - eco_file: str = _eco_file, - out_file: str = 'temp_out.nc', - verbose: bool = False) -> None: +def main( + satellite: str = "snpp", + sensor: str = "viirs", + data_path: str = _datapath, + mod02: str = "", # _fname_mod02, + mod03: str = "", # _fname_mod03, + img02: str = "", # _fname_img02, + img03: str = "", # _fname_img03, + threshold_file: str = _threshold_file, + geos_atm_1: str = _geos_atm_1, + geos_atm_2: str = _geos_atm_2, + geos_land: str = _geos_land, + geos_ocean: str = _geos_ocean, + geos_constants: str = _geos_constants, + ndvi_file: str = _ndvi_file, + sst_file: str = _sst_file, + eco_file: str = _eco_file, + out_file: str = "temp_out.nc", + verbose: bool = False, +) -> None: """Main function for the MVCM. Parameters @@ -125,20 +126,21 @@ def main(satellite: str = 'snpp', else: logging.basicConfig(level=logging.WARNING) - file_names = {'MOD02': f'{mod02}', - 'MOD03': f'{mod03}', - 'IMG02': f'{img02}', - 'IMG03': f'{img03}', - 'GEOS_atm_1': f'{geos_atm_1}', - 'GEOS_atm_2': f'{geos_atm_2}', - 'GEOS_land': f'{geos_land}', - 'GEOS_ocean': f'{geos_ocean}', - 'GEOS_constants': f'{geos_constants}', - 'NDVI': f'{ndvi_file}', - 'SST': f'{sst_file}', - 'ECO': f'{eco_file}', - 'ANC_DIR': f'{data_path}/ancillary' - } + file_names = { + "MOD02": f"{mod02}", + "MOD03": f"{mod03}", + "IMG02": f"{img02}", + "IMG03": f"{img03}", + "GEOS_atm_1": f"{geos_atm_1}", + "GEOS_atm_2": f"{geos_atm_2}", + "GEOS_land": f"{geos_land}", + "GEOS_ocean": f"{geos_ocean}", + "GEOS_constants": f"{geos_constants}", + "NDVI": f"{ndvi_file}", + "SST": f"{sst_file}", + "ECO": f"{eco_file}", + "ANC_DIR": f"{data_path}/ancillary", + } with open(threshold_file) as f: text = f.read() @@ -146,7 +148,7 @@ def main(satellite: str = 'snpp', ################# ################# - if (img02 is None or img03 is None): + if img02 is None or img03 is None: use_hires = False else: use_hires = True @@ -154,7 +156,7 @@ def main(satellite: str = 'snpp', ################# viirs_data = rd.get_data(satellite, sensor, file_names, hires=use_hires) - sunglint_angle = thresholds['Sun_Glint']['bounds'][3] + sunglint_angle = thresholds["Sun_Glint"]["bounds"][3] scene_flags = scn.find_scene(viirs_data, sunglint_angle) Restore = Restoral(data=viirs_data, thresholds=thresholds, scene_flag=scene_flags) @@ -162,7 +164,30 @@ def main(satellite: str = 'snpp', cmin_3d = np.ones((18, viirs_data.M11.shape[0], viirs_data.M11.shape[1])) ########################################################## - _bitlist = ['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12', '13', '14', '15', '16'] + _bitlist = [ + "01", + "02", + "03", + "04", + "05", + "06", + "07", + "08", + "09", + "10", + "11", + "12", + "13", + "14", + "15", + "16", + "r01", + "r02", + "r03", + "r04", + "r05", + "r06", + ] # bits, bit, qabit = {}, {}, {} bits = {} # for b in _bitlist: @@ -170,26 +195,43 @@ def main(satellite: str = 'snpp', # qabit[b] = np.zeros(viirs_data.M11.shape) for b in _bitlist: - bits[b] = {'test': np.zeros(viirs_data.M11.shape), - 'qa': np.zeros(viirs_data.M11.shape)} + bits[b] = { + "test": np.zeros(viirs_data.M11.shape), + "qa": np.zeros(viirs_data.M11.shape), + } i = 0 # temp_bit = np.zeros((18, viirs_data.M11.shape[0], viirs_data.M11.shape[1])) - for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast', - 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean', - 'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert', - 'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Land_Night', - 'Polar_Night_Land', 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']: - print(f'Processing {scene_name}') + for scene_name in [ + "Land_Day", + "Land_Day_Coast", + "Land_Day_Desert", + "Land_Day_Desert_Coast", + "Ocean_Day", + "Ocean_Night", + "Polar_Day_Ocean", + "Polar_Night_Ocean", + "Polar_Day_Land", + "Polar_Day_Coast", + "Polar_Day_Desert", + "Polar_Day_Desert_Coast", + "Polar_Day_Snow", + "Land_Night", + "Polar_Night_Land", + "Polar_Night_Snow", + "Day_Snow", + "Night_Snow", + ]: + print(f"Processing {scene_name}") if np.all(viirs_data[scene_name].values == 0): - logging.info('Skipping, no pixels in scene.') + logging.info("Skipping, no pixels in scene.") continue - MyScene = tst.CloudTests(data=viirs_data, - scene_name=scene_name, - thresholds=thresholds) + MyScene = tst.CloudTests( + data=viirs_data, scene_name=scene_name, thresholds=thresholds + ) # Initialize the confidence arrays for the various test groups cmin_G1 = np.ones(viirs_data.M11.shape) @@ -200,84 +242,162 @@ def main(satellite: str = 'snpp', # cmin_temp = np.ones(viirs_data.M11.shape) if use_hires is True: - m15_name = 'M15hi' + m15_name = "M15hi" else: - m15_name = 'M15' - - cmin_G1, bits['01'] = MyScene.test_11um(m15_name, cmin_G1, bits['01']) - cmin_G1, bits['02'] = MyScene.surface_temperature_test('', viirs_data, cmin_G1, bits['02']) - cmin_G1, bits['03'] = MyScene.sst_test('M15', 'M16', cmin_G1, bits['03']) - - cmin_G2, bits['04'] = MyScene.bt_diff_86_11um('M14-M15', cmin_G2, bits['04']) - cmin_G2, bits['05'] = MyScene.test_11_12um_diff('M15-M16', cmin_G2, bits['05']) - cmin_G2, bits['06'] = MyScene.variability_11um_test(m15_name, cmin_G2, bits['06']) - cmin_G2, bits['07'] = MyScene.bt_difference_11_4um_test_ocean('M15-M12', cmin_G2, bits['07']) - # # cmin_G2, bits['08'] = MyScene.bt_difference_11_4um_test_land('M15-M12', cmin_G2), bits['08'] - cmin_G2, bits['09'] = MyScene.oceanic_stratus_11_4um_test('M15-M12', cmin_G2, bits['09']) - - cmin_G3, bits['10'] = MyScene.nir_reflectance_test('M07', cmin_G3, bits['10']) + m15_name = "M15" + + cmin_G1, bits["01"] = MyScene.test_11um(m15_name, cmin_G1, bits["01"]) + cmin_G1, bits["02"] = MyScene.surface_temperature_test( + "", viirs_data, cmin_G1, bits["02"] + ) + cmin_G1, bits["03"] = MyScene.sst_test("M15", "M16", cmin_G1, bits["03"]) + + cmin_G2, bits["04"] = MyScene.bt_diff_86_11um("M14-M15", cmin_G2, bits["04"]) + cmin_G2, bits["05"] = MyScene.test_11_12um_diff("M15-M16", cmin_G2, bits["05"]) + cmin_G2, bits["06"] = MyScene.variability_11um_test( + m15_name, cmin_G2, bits["06"] + ) + cmin_G2, bits["07"] = MyScene.bt_difference_11_4um_test_ocean( + "M15-M12", cmin_G2, bits["07"] + ) + cmin_G2, bits["08"] = MyScene.bt_difference_11_4um_test_land( + "M15-M12", cmin_G2, bits["08"] + ) + cmin_G2, bits["09"] = MyScene.oceanic_stratus_11_4um_test( + "M15-M12", cmin_G2, bits["09"] + ) + + cmin_G3, bits["10"] = MyScene.nir_reflectance_test("M07", cmin_G3, bits["10"]) # I need to make sure that this test is performing as intended. Right now it doesn't # seem to have an impact in any way whatsoever - cmin_G3, bits['11'] = MyScene.vis_nir_ratio_test('M07-M05ratio', cmin_G3, bits['11']) - cmin_G3, bits['12'] = MyScene.test_16_21um_reflectance('M10', cmin_G3, bits['12']) + cmin_G3, bits["11"] = MyScene.vis_nir_ratio_test( + "M07-M05ratio", cmin_G3, bits["11"] + ) + cmin_G3, bits["12"] = MyScene.test_16_21um_reflectance( + "M10", cmin_G3, bits["12"] + ) # this test seems to be almost working, but not 100% there yet - cmin_G3, bits['13'] = MyScene.visible_reflectance_test('M128', cmin_G3, bits['13']) - cmin_G3, bits['11'] = MyScene.gemi_test('GEMI', cmin_G3, bits['11']) + cmin_G3, bits["13"] = MyScene.visible_reflectance_test( + "M128", cmin_G3, bits["13"] + ) + cmin_G3, bits["11"] = MyScene.gemi_test("GEMI", cmin_G3, bits["11"]) - cmin_G4, bits['15'] = MyScene.test_1_38um_high_clouds('M09', cmin_G4, bits['15']) + cmin_G4, bits["15"] = MyScene.test_1_38um_high_clouds( + "M09", cmin_G4, bits["15"] + ) - cmin_G5, bits['16'] = MyScene.thin_cirrus_4_12um_BTD_test('M13-M16', cmin_G5, bits['16']) + cmin_G5, bits["16"] = MyScene.thin_cirrus_4_12um_BTD_test( + "M13-M16", cmin_G5, bits["16"] + ) bit = {} + restoral_bit = {} for b in _bitlist: - bit[b] = bits[f'{b}']['test'] + bit[b] = bits[f"{b}"]["test"] + restoral_bit[b] = utils.restoral_flag(bits[f"{b}"]) # if utils.group_count(bit) != 0: - ''' + """ cmin_3d[i, :, :] = np.fmin(cmin_temp, np.power(cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5, 1/utils.group_count(bit))) - ''' - cmin = np.power(cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5, 1/utils.group_count(bit)) - - sunglint_bits = bit['01'] * bit['03'] * bit['05'] * bit['15'] - sh_water_bits = {'ir': bit['01'] * bit['05'], - 'nir_1': bit['10'], - 'nir_2': bit['15'], - 'sst': bit['03'] - } - land_bits = bit['15'] * bit['05'] + """ + cmin = np.power( + cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5, 1 / utils.group_count(bit) + ) + + sunglint_bits = ( + restoral_bit["01"] + * restoral_bit["03"] + * restoral_bit["05"] + * restoral_bit["15"] + ) + sh_water_bits = { + "ir": restoral_bit["01"] * restoral_bit["05"], + "nir_1": restoral_bit["10"], + "nir_2": restoral_bit["15"], + "sst": restoral_bit["03"], + } + desert_bits = restoral_bit["15"] * restoral_bit["05"] # make sure that the land_bits are calculated properly - idx = np.nonzero(scene_flags['desert'] == 0) - land_bits[idx] = land_bits[idx] * bit['10'][idx] * bit['07'][idx] * bit['11'][idx] - coast_bits = bit['05'] - land_night_bits = bit['16'] - - idx = np.nonzero((scene_flags['uniform'] == 1) & (scene_flags['water'] == 1) & - (scene_flags['ice'] == 0) & (cmin >= 0.05) & (cmin <= 0.99)) - cmin[idx] = Restore.spatial_variability(copy(cmin))[idx] - - idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['water'] == 1) & - (scene_flags['uniform'] == 1) & (scene_flags['sunglint'] == 1) & - (cmin <= 0.95)) - cmin[idx] = Restore.sunglint(sunglint_bits, copy(cmin))[idx] - - idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['water'] == 1) & - (scene_flags['ice'] == 0)) - cmin[idx] = Restore.shallow_water(sh_water_bits, copy(cmin))[idx] - - idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) & - (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & - (cmin <= 0.95)) - cmin[idx] = Restore.land(land_bits, copy(cmin))[idx] - - idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) & - (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & - (scene_flags['coast'] == 1)) - cmin = Restore.coast(coast_bits, copy(cmin)) - - idx = np.nonzero((scene_flags['night'] == 1) & (scene_flags['land'] == 1) & - (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & - (cmin <= 0.95)) - # cmin[idx] = Restore.land_night(land_night_bits, copy(cmin))[idx] + idx = np.nonzero(scene_flags["desert"] == 0) + land_bits = np.zeros(desert_bits.shape) + land_bits[idx] = ( + desert_bits[idx] + * restoral_bit["10"][idx] + * restoral_bit["08"][idx] + * restoral_bit["09"][idx] + * restoral_bit["11"][idx] + ) + coast_bits = restoral_bit["05"] + land_night_bits = restoral_bit["16"] + + idx = np.nonzero( + (scene_flags["uniform"] == 1) + & (scene_flags["water"] == 1) + & (scene_flags["ice"] == 0) + & (cmin >= 0.05) + & (cmin <= 0.99) + ) + cmin[idx], bits["r01"]["test"][idx] = Restore.spatial_variability( + cmin, idx, bits["r01"] + ) + + idx = np.nonzero( + (scene_flags["day"] == 1) + & (scene_flags["water"] == 1) + & (scene_flags["uniform"] == 1) + & (scene_flags["sunglint"] == 1) + & (cmin <= 0.95) + ) + cmin[idx], bits["r02"]["test"][idx] = Restore.sunglint( + sunglint_bits, copy(cmin), idx, bits["r02"] + ) + + idx = np.nonzero( + (scene_flags["day"] == 1) + & (scene_flags["water"] == 1) + & (scene_flags["ice"] == 0) + ) + cmin[idx], bits["r03"]["test"][idx] = Restore.shallow_water( + sh_water_bits, copy(cmin), idx, bits["r03"] + ) + + idx = np.nonzero( + (scene_flags["day"] == 1) + & (scene_flags["land"] == 1) + & (scene_flags["snow"] == 0) + & (scene_flags["ice"] == 0) + & (cmin <= 0.95) + & (scene_flags["desert"] == 0) + ) + cmin[idx], bits["r04"]["test"][idx] = Restore.land( + land_bits, copy(cmin), idx, bits["r04"] + ) + idx = np.nonzero( + (scene_flags["day"] == 1) + & (scene_flags["desert"] == 1) + & (scene_flags["snow"] == 0) + & (scene_flags["ice"] == 0) + & (cmin <= 0.95) + ) + cmin[idx], bits["r04"]["test"][idx] = Restore.land( + desert_bits, copy(cmin), idx, bits["r04"] + ) + + idx = np.nonzero( + (scene_flags["day"] == 1) + & (scene_flags["land"] == 1) + & (scene_flags["snow"] == 0) + & (scene_flags["ice"] == 0) + & (scene_flags["coast"] == 1) + ) + cmin[idx], bits["r05"]["test"][idx] = Restore.coast( + coast_bits, copy(cmin), idx, bits["r05"] + ) + + # idx = np.nonzero((scene_flags['night'] == 1) & (scene_flags['land'] == 1) & + # (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & + # (cmin <= 0.95)) + # cmin[idx], bits['r06']['test'][idx] = Restore.land_night(land_night_bits, copy(cmin), bits['r06'])[idx] cmin_3d[i, :, :] = cmin i += 1 @@ -285,7 +405,7 @@ def main(satellite: str = 'snpp', cmin_final = np.min(cmin_3d, axis=0) # bit['15'] = np.max(temp_bit, axis=0) - ''' + """ idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['water'] == 1) & (scene_flags['ice'] == 0)) cmin[idx] = Restore.shallow_water(sh_water_bits, cmin)[idx] @@ -304,7 +424,7 @@ def main(satellite: str = 'snpp', (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) & (cmin <= 0.95)) cmin[idx] = Restore.land_night(land_night_bits, cmin)[idx] - ''' + """ # bits translation MVCM-python -> MVCM-C # 01 13 test_11um # 02 27 surface_temperature_test @@ -336,63 +456,47 @@ def main(satellite: str = 'snpp', # # cmin = MVCM.clear_sky_restoral(cmin) # # # return cmin - d = {'latitude': {'dims': ('x', 'y'), - 'data': viirs_data.latitude.values}, - 'longitude': {'dims': ('x', 'y'), - 'data': viirs_data.longitude.values}, - 'confidence': {'dims': ('x', 'y'), - 'data': cmin_final}, - 'ndvi': {'dims': ('x', 'y'), - 'data': viirs_data.ndvi.values}, - 'sunglint_bits': {'dims': ('x', 'y'), - 'data': sunglint_bits}, - 'eco': {'dims': ('x', 'y'), - 'data': viirs_data.ecosystem.values}, - 'bit1': {'dims': ('x', 'y'), - 'data': bit['01']}, - 'bit2': {'dims': ('x', 'y'), - 'data': bit['02']}, - 'bit3': {'dims': ('x', 'y'), - 'data': bit['03']}, - 'bit4': {'dims': ('x', 'y'), - 'data': bit['04']}, - 'bit5': {'dims': ('x', 'y'), - 'data': bit['05']}, - 'bit6': {'dims': ('x', 'y'), - 'data': bit['06']}, - 'bit7': {'dims': ('x', 'y'), - 'data': bit['07']}, - 'bit8': {'dims': ('x', 'y'), - 'data': bit['08']}, - 'bit9': {'dims': ('x', 'y'), - 'data': bit['09']}, - 'bit10': {'dims': ('x', 'y'), - 'data': bit['10']}, - 'bit11': {'dims': ('x', 'y'), - 'data': bit['11']}, - 'bit12': {'dims': ('x', 'y'), - 'data': bit['12']}, - 'bit13': {'dims': ('x', 'y'), - 'data': bit['13']}, - 'bit14': {'dims': ('x', 'y'), - 'data': bit['14']}, - 'bit15': {'dims': ('x', 'y'), - 'data': bit['15']}, - 'bit16': {'dims': ('x', 'y'), - 'data': bit['16']}, - # 'thr': {'dims': ('x', 'y'), - # 'data': thr_temp}, - # 'shallow_water_bits': {'dims': ('x', 'y'), - # 'data': sh_water_bits}, - # 'b15': {'dims': ('x', 'y'), - # 'data': viirs_data.M15.values}, - # 'b16': {'dims': ('x', 'y'), - # 'data': viirs_data.M16.values}, - 'GEMI': {'dims': ('x', 'y'), - 'data': viirs_data.GEMI.values}, - } - - sf = {k: {'dims': ('x', 'y'), 'data': scene_flags[k]} for k in scene_flags} + d = { + "latitude": {"dims": ("x", "y"), "data": viirs_data.latitude.values}, + "longitude": {"dims": ("x", "y"), "data": viirs_data.longitude.values}, + "confidence": {"dims": ("x", "y"), "data": cmin_final}, + "ndvi": {"dims": ("x", "y"), "data": viirs_data.ndvi.values}, + "sunglint_bits": {"dims": ("x", "y"), "data": sunglint_bits}, + "eco": {"dims": ("x", "y"), "data": viirs_data.ecosystem.values}, + "bit1": {"dims": ("x", "y"), "data": bit["01"]}, + "bit2": {"dims": ("x", "y"), "data": bit["02"]}, + "bit3": {"dims": ("x", "y"), "data": bit["03"]}, + "bit4": {"dims": ("x", "y"), "data": bit["04"]}, + "bit5": {"dims": ("x", "y"), "data": bit["05"]}, + "bit6": {"dims": ("x", "y"), "data": bit["06"]}, + "bit7": {"dims": ("x", "y"), "data": bit["07"]}, + "bit8": {"dims": ("x", "y"), "data": bit["08"]}, + "bit9": {"dims": ("x", "y"), "data": bit["09"]}, + "bit10": {"dims": ("x", "y"), "data": bit["10"]}, + "bit11": {"dims": ("x", "y"), "data": bit["11"]}, + "bit12": {"dims": ("x", "y"), "data": bit["12"]}, + "bit13": {"dims": ("x", "y"), "data": bit["13"]}, + "bit14": {"dims": ("x", "y"), "data": bit["14"]}, + "bit15": {"dims": ("x", "y"), "data": bit["15"]}, + "bit16": {"dims": ("x", "y"), "data": bit["16"]}, + "bit_r1": {"dims": ("x", "y"), "data": bit["r01"]}, + "bit_r2": {"dims": ("x", "y"), "data": bit["r02"]}, + "bit_r3": {"dims": ("x", "y"), "data": bit["r03"]}, + "bit_r4": {"dims": ("x", "y"), "data": bit["r04"]}, + "bit_r5": {"dims": ("x", "y"), "data": bit["r05"]}, + "bit_r6": {"dims": ("x", "y"), "data": bit["r06"]}, + # 'thr': {'dims': ('x', 'y'), + # 'data': thr_temp}, + # 'shallow_water_bits': {'dims': ('x', 'y'), + # 'data': sh_water_bits}, + # 'b15': {'dims': ('x', 'y'), + # 'data': viirs_data.M15.values}, + # 'b16': {'dims': ('x', 'y'), + # 'data': viirs_data.M16.values}, + "GEMI": {"dims": ("x", "y"), "data": viirs_data.GEMI.values}, + } + + sf = {k: {"dims": ("x", "y"), "data": scene_flags[k]} for k in scene_flags} out_dict = d | sf ds_out = xr.Dataset.from_dict(out_dict) @@ -411,14 +515,19 @@ def main(satellite: str = 'snpp', def test_main(): - - rad1 = [[255, 260, 265, 248, 223], - [278, 285, 270, 268, 256], - [275, 273, 266, 254, 259]] - rad2 = [[270, 273, 271, 268, 265], - [277, 286, 275, 277, 269], - [280, 281, 272, 270, 267]] - thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml' + rad1 = [ + [255, 260, 265, 248, 223], + [278, 285, 270, 268, 256], + [275, 273, 266, 254, 259], + ] + rad2 = [ + [270, 273, 271, 268, 265], + [277, 286, 275, 277, 269], + [280, 281, 272, 270, 267], + ] + thresh_file = ( + "/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml" + ) with open(thresh_file) as f: text = f.read() @@ -430,11 +539,11 @@ def test_main(): thresholds = yml.safe_load(text) - confidence[0, :, :] = tests.test_11um(rad1, thresholds['Daytime_Ocean']) - confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thresholds['Daytime_Ocean']) + confidence[0, :, :] = tests.test_11um(rad1, thresholds["Daytime_Ocean"]) + confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thresholds["Daytime_Ocean"]) - print(f'Confidence[0,:,:]: \n {confidence[0, :, :]}') - print(f'Confidence[1,:,:]: \n {confidence[1, :, :]}') + print(f"Confidence[0,:,:]: \n {confidence[0, :, :]}") + print(f"Confidence[1,:,:]: \n {confidence[1, :, :]}") return confidence diff --git a/mvcm/preprocess_thresholds.py b/mvcm/preprocess_thresholds.py index 4bd4e92..bdd57b5 100644 --- a/mvcm/preprocess_thresholds.py +++ b/mvcm/preprocess_thresholds.py @@ -1,50 +1,50 @@ +from typing import Dict + import numpy as np import xarray as xr +from numpy.lib.stride_tricks import sliding_window_view import ancillary_data as anc import mvcm.utils as utils -from numpy.lib.stride_tricks import sliding_window_view -from typing import Dict - -_dtr = np.pi/180 -_DTR = np.pi/180 +_dtr = np.pi / 180 +_DTR = np.pi / 180 -def prepare_11_12um_thresholds(thresholds: np.ndarray, - dim1: int) -> Dict: - +def prepare_11_12um_thresholds(thresholds: np.ndarray, dim1: int) -> Dict: coeff_values = np.empty((dim1, 2)) - coeff_values[:, 0] = np.full(dim1, thresholds['coeffs'][0]) - coeff_values[:, 1] = np.full(dim1, thresholds['coeffs'][1]) - cmult_values = np.full(dim1, thresholds['cmult']) - adj_values = np.full(dim1, thresholds['adj']) - if 'bt1' in list(thresholds): - bt1 = np.full(dim1, thresholds['bt1']) + coeff_values[:, 0] = np.full(dim1, thresholds["coeffs"][0]) + coeff_values[:, 1] = np.full(dim1, thresholds["coeffs"][1]) + cmult_values = np.full(dim1, thresholds["cmult"]) + adj_values = np.full(dim1, thresholds["adj"]) + if "bt1" in list(thresholds): + bt1 = np.full(dim1, thresholds["bt1"]) else: bt1 = np.full(dim1, -999) - if 'lat' in list(thresholds): - lat = np.full(dim1, thresholds['lat']) + if "lat" in list(thresholds): + lat = np.full(dim1, thresholds["lat"]) else: lat = np.full(dim1, -999) - thr_dict = {'coeffs': coeff_values, - 'cmult': cmult_values, - 'adj': adj_values, - 'bt1': bt1, - 'lat': lat, - } + thr_dict = { + "coeffs": coeff_values, + "cmult": cmult_values, + "adj": adj_values, + "bt1": bt1, + "lat": lat, + } return thr_dict -def thresholds_11_12um(data: xr.Dataset, - thresholds: np.ndarray, - scene: str, - scene_idx: np.ndarray) -> np.ndarray: +def thresholds_11_12um( + data: xr.Dataset, thresholds: np.ndarray, scene: str, scene_idx: np.ndarray +) -> np.ndarray: cosvza = np.cos(data.sensor_zenith.values[scene_idx].ravel() * _DTR) schi = np.full(cosvza.shape, 99.0) - schi[np.abs(cosvza) > 0] = 1.0/cosvza - schi = np.array(schi, dtype=np.float32) # this is because the C function expects a float + schi[np.abs(cosvza) > 0] = 1.0 / cosvza + schi = np.array( + schi, dtype=np.float32 + ) # this is because the C function expects a float m15 = np.array(data.M15.values[scene_idx].ravel(), dtype=np.float32) latitude = data.latitude.values[scene_idx].ravel() @@ -52,84 +52,120 @@ def thresholds_11_12um(data: xr.Dataset, thr_dict = prepare_11_12um_thresholds(thresholds, m15.shape[0]) midpt = np.full(m15.shape[0], thr) - idx = np.nonzero((thr < 0.1) | (np.abs(schi-99.0) < 0.0001)) - midpt[idx] = thr_dict['coeffs'][idx, 0] - locut = midpt + (thr_dict['cmult'] * midpt) - - if scene in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast', - 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean']: - hicut = midpt - thr_dict['adj'] - - elif scene in ['Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert', - 'Polar_Day_Desert_Coast', 'Polar_Day_Snow']: - hicut = midpt - (thr_dict['adj'] * midpt) - - elif scene in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']: - _coeffs = {'Land_Night': 0.3, 'Polar_Night_Land': 0.3, 'Polar_Night_Snow': 0.3, - 'Day_Snow': 0.0, 'Night_Snow': 0.3} + idx = np.nonzero((thr < 0.1) | (np.abs(schi - 99.0) < 0.0001)) + midpt[idx] = thr_dict["coeffs"][idx, 0] + locut = midpt + (thr_dict["cmult"] * midpt) + + if scene in [ + "Land_Day", + "Land_Day_Coast", + "Land_Day_Desert", + "Land_Day_Desert_Coast", + "Ocean_Day", + "Ocean_Night", + "Polar_Day_Ocean", + "Polar_Night_Ocean", + ]: + hicut = midpt - thr_dict["adj"] + + elif scene in [ + "Polar_Day_Land", + "Polar_Day_Coast", + "Polar_Day_Desert", + "Polar_Day_Desert_Coast", + "Polar_Day_Snow", + ]: + hicut = midpt - (thr_dict["adj"] * midpt) + + elif scene in [ + "Land_Night", + "Polar_Night_Land", + "Polar_Night_Snow", + "Day_Snow", + "Night_Snow", + ]: + _coeffs = { + "Land_Night": 0.3, + "Polar_Night_Land": 0.3, + "Polar_Night_Snow": 0.3, + "Day_Snow": 0.0, + "Night_Snow": 0.3, + } midpt = midpt - (_coeffs[scene] * locut) - if scene in ['Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']: + if scene in ["Polar_Night_Land", "Polar_Night_Snow", "Night_Snow"]: hicut = np.full(m15.shape, midpt - 1.25) - idx = np.nonzero(m15 < thr_dict['bt1']) + idx = np.nonzero(m15 < thr_dict["bt1"]) hicut[idx] = midpt[idx] - (0.2 * locut[idx]) - elif scene in ['Land_Night']: + elif scene in ["Land_Night"]: hicut = np.full(m15.shape, 1.25) - idx = np.nonzero((m15 < thr_dict['bt1']) & (latitude > thr_dict['lat'])) - hicut[idx] = -0.1 - np.power(90.0 - np.abs(latitude[idx])/60, 4) * 1.15 - elif scene in ['Day_Snow']: - hicut = locut - (thr_dict['cmult'] * locut) + idx = np.nonzero((m15 < thr_dict["bt1"]) & (latitude > thr_dict["lat"])) + hicut[idx] = -0.1 - np.power(90.0 - np.abs(latitude[idx]) / 60, 4) * 1.15 + elif scene in ["Day_Snow"]: + hicut = locut - (thr_dict["cmult"] * locut) else: - raise ValueError('Scene name not valid') + raise ValueError("Scene name not valid") else: - raise ValueError('Scene name not valid') + raise ValueError("Scene name not valid") - thr_out = np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape))) + thr_out = np.dstack( + (locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape)) + ) return np.squeeze(thr_out.T) def thresholds_NIR(data, thresholds, scene, test_name, scene_idx): - sza = data.solar_zenith.values[scene_idx].ravel() # NOTE: the visud condition in the C code is equivalent to having sza <= 85 # For the time being the visud filtering is not implemented - c = np.array(thresholds[scene][test_name]['coeffs']) - if test_name == 'NIR_Reflectance_Test': + c = np.array(thresholds[scene][test_name]["coeffs"]) + if test_name == "NIR_Reflectance_Test": band_n = 2 - vzcpow = thresholds['VZA_correction']['vzcpow'][0] - elif test_name == '1.6_2.1um_NIR_Reflectance_Test': + vzcpow = thresholds["VZA_correction"]["vzcpow"][0] + elif test_name == "1.6_2.1um_NIR_Reflectance_Test": band_n = 7 - vzcpow = thresholds['VZA_correction']['vzcpow'][2] + vzcpow = thresholds["VZA_correction"]["vzcpow"][2] else: - raise ValueError('Test name not valid') + raise ValueError("Test name not valid") refang = data.sunglint_angle.values[scene_idx].ravel() - sunglint_thresholds = thresholds['Sun_Glint'] + sunglint_thresholds = thresholds["Sun_Glint"] sunglint_flag = utils.sunglint_scene(refang, sunglint_thresholds) nir_thresh = thresholds[scene][test_name] - if test_name == 'NIR_Reflectance_Test': - hicut = c[0] + c[1]*sza + c[2]*np.power(sza, 2) + c[3]*np.power(sza, 3) - elif test_name == '1.6_2.1um_NIR_Reflectance_Test': - hicut = c[0] + c[1]*sza + c[2]*np.power(sza, 2) + c[3]*np.power(sza, 3) + c[4]*np.power(sza, 4) + if test_name == "NIR_Reflectance_Test": + hicut = c[0] + c[1] * sza + c[2] * np.power(sza, 2) + c[3] * np.power(sza, 3) + elif test_name == "1.6_2.1um_NIR_Reflectance_Test": + hicut = ( + c[0] + + c[1] * sza + + c[2] * np.power(sza, 2) + + c[3] * np.power(sza, 3) + + c[4] * np.power(sza, 4) + ) else: pass - hicut = (hicut * 0.01) + nir_thresh['adj'] - hicut = (hicut * nir_thresh['bias']) - midpt = hicut + (nir_thresh['midpt_coeff'] * nir_thresh['bias']) - locut = midpt + (nir_thresh['locut_coeff'] * nir_thresh['bias']) - thr = np.array([locut, midpt, hicut, nir_thresh['thr'][3]*np.ones(refang.shape)]) + hicut = (hicut * 0.01) + nir_thresh["adj"] + hicut = hicut * nir_thresh["bias"] + midpt = hicut + (nir_thresh["midpt_coeff"] * nir_thresh["bias"]) + locut = midpt + (nir_thresh["locut_coeff"] * nir_thresh["bias"]) + thr = np.array([locut, midpt, hicut, nir_thresh["thr"][3] * np.ones(refang.shape)]) - cosvza = np.cos(data.sensor_zenith.values[scene_idx]*_DTR).ravel() + cosvza = np.cos(data.sensor_zenith.values[scene_idx] * _DTR).ravel() corr_thr = np.zeros((4, refang.shape[0])) - corr_thr[:3, sunglint_flag == 0] = thr[:3, sunglint_flag == 0] * (1./np.power(cosvza[sunglint_flag == 0], vzcpow)) + corr_thr[:3, sunglint_flag == 0] = thr[:3, sunglint_flag == 0] * ( + 1.0 / np.power(cosvza[sunglint_flag == 0], vzcpow) + ) corr_thr[3, sunglint_flag == 0] = thr[3, sunglint_flag == 0] for flag in range(1, 4): if len(refang[sunglint_flag == flag]) > 0: - dosgref = utils.get_sunglint_thresholds(refang, sunglint_thresholds, band_n, flag, thr) - corr_thr[:3, sunglint_flag == flag] = dosgref[:3, sunglint_flag == flag] * \ - (1./np.power(cosvza[sunglint_flag == flag], vzcpow)) + dosgref = utils.get_sunglint_thresholds( + refang, sunglint_thresholds, band_n, flag, thr + ) + corr_thr[:3, sunglint_flag == flag] = dosgref[:3, sunglint_flag == flag] * ( + 1.0 / np.power(cosvza[sunglint_flag == flag], vzcpow) + ) corr_thr[3, sunglint_flag == flag] = dosgref[3, sunglint_flag == flag] return corr_thr @@ -137,10 +173,10 @@ def thresholds_NIR(data, thresholds, scene, test_name, scene_idx): def thresholds_surface_temperature(data, thresholds, scene_idx): # def preproc_surf_temp(data, thresholds): - thr_sfc1 = thresholds['desert_thr'] - thr_sfc2 = thresholds['regular_thr'] - thr_df1 = thresholds['channel_diff_11-12um_thr'] - thr_df2 = thresholds['channel_diff_11-4um_thr'] + thr_sfc1 = thresholds["desert_thr"] + thr_sfc2 = thresholds["regular_thr"] + thr_df1 = thresholds["channel_diff_11-12um_thr"] + thr_df2 = thresholds["channel_diff_11-4um_thr"] max_vza = 70.13 # This values is set based on sensor. Check mask_processing_constants.h for MODIS value df1 = (data.M15 - data.M16).values[scene_idx].ravel() @@ -148,21 +184,26 @@ def thresholds_surface_temperature(data, thresholds, scene_idx): desert_flag = data.Desert.values[scene_idx].ravel() thresh = np.ones(df1.shape) * thr_sfc1 - idx = np.where((df1 >= thr_df1[0]) | ((df1 < thr_df1[0]) & ((df2 <= thr_df2[0]) | (df2 >= thr_df2[1])))) + idx = np.where( + (df1 >= thr_df1[0]) + | ((df1 < thr_df1[0]) & ((df2 <= thr_df2[0]) | (df2 >= thr_df2[1]))) + ) thresh[idx] = thr_sfc2 idx = np.where(desert_flag == 1) thresh[idx] == thr_sfc1 midpt = thresh idx = np.where(df1 >= thr_df1[1]) - midpt[idx] = thresh[idx] + 2.0*df1[idx] + midpt[idx] = thresh[idx] + 2.0 * df1[idx] - corr = np.power(data.sensor_zenith.values[scene_idx].ravel()/max_vza, 4) * 3.0 + corr = np.power(data.sensor_zenith.values[scene_idx].ravel() / max_vza, 4) * 3.0 midpt = midpt + corr locut = midpt + 2.0 hicut = midpt - 2.0 - thr_out = np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape))) + thr_out = np.dstack( + (locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape)) + ) return np.squeeze(thr_out.T) @@ -173,11 +214,17 @@ def preproc_sst(data, thresholds): m32c = data.M16 - 273.16 m31c_m32c = m31c - m32c sstc = data.geos_sfct - 273.16 - cosvza = np.cos(data.sensor_zenith*_DTR) + cosvza = np.cos(data.sensor_zenith * _DTR) - a = thresholds['coeffs'] + a = thresholds["coeffs"] - modsst = 273.16 + a[0] + a[1]*m31c + a[2]*m31c_m32c*sstc + a[3]*m31c_m32c*((1/cosvza) - 1) + modsst = ( + 273.16 + + a[0] + + a[1] * m31c + + a[2] * m31c_m32c * sstc + + a[3] * m31c_m32c * ((1 / cosvza) - 1) + ) sfcdif = data.geos_sfct - modsst return sfcdif @@ -186,9 +233,11 @@ def preproc_sst(data, thresholds): def var_11um(data, thresholds): rad = data.M15.values var = np.zeros((rad.shape[0], rad.shape[1], 9)) - var_thr = thresholds['Daytime_Ocean_Spatial_Variability']['dovar11'] + var_thr = thresholds["Daytime_Ocean_Spatial_Variability"]["dovar11"] - test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3)) + test = sliding_window_view( + np.pad(rad, [1, 1], mode="constant"), (3, 3) + ) - np.expand_dims(rad, (2, 3)) var[np.abs(test).reshape(rad.shape[0], rad.shape[1], 9) < var_thr] = 1 var = var.sum(axis=2) @@ -203,16 +252,16 @@ def get_b1_thresholds(data, thresholds, scene_idx): # this is hardcoded in the function delta_ndvi_bin = 0.1 - des_ndvi = thresholds['Misc']['des_ndvi'] - thr_adj_fac_desert = thresholds['Misc']['adj_fac_desert'] - thr_adj_fac_land = thresholds['Misc']['adj_fac_land'] - ndvi_bnd1 = thresholds['Misc']['ndvi_bnd1'] - ndvi_bnd2 = thresholds['Misc']['ndvi_bnd2'] - fill_ndvi = thresholds['Misc']['fill_ndvi'] + des_ndvi = thresholds["Misc"]["des_ndvi"] + thr_adj_fac_desert = thresholds["Misc"]["adj_fac_desert"] + thr_adj_fac_land = thresholds["Misc"]["adj_fac_land"] + ndvi_bnd1 = thresholds["Misc"]["ndvi_bnd1"] + ndvi_bnd2 = thresholds["Misc"]["ndvi_bnd2"] + fill_ndvi = thresholds["Misc"]["fill_ndvi"] - coeff1 = np.array(thresholds['Coeffs_Band1_land_thresh']).reshape(10, 3, 4) + coeff1 = np.array(thresholds["Coeffs_Band1_land_thresh"]).reshape(10, 3, 4) coeff2 = np.zeros((10, 3, 4)) - coeff2[:3, :, :] = np.array(thresholds['Coeffs_Band8_land_thresh']).reshape(3, 3, 4) + coeff2[:3, :, :] = np.array(thresholds["Coeffs_Band8_land_thresh"]).reshape(3, 3, 4) coeff = np.stack((coeff1, coeff2)) indvi = np.zeros(ndvi.shape) @@ -223,12 +272,12 @@ def get_b1_thresholds(data, thresholds, scene_idx): # this is equivalent to interp=1 in the C code idx = np.nonzero((ndvi >= ndvi_bnd1) & (ndvi < ndvi_bnd2)) - indvi[idx] = (ndvi[idx]/delta_ndvi_bin) - 0.5 + indvi[idx] = (ndvi[idx] / delta_ndvi_bin) - 0.5 indvi[idx] = np.trunc(np.clip(indvi[idx], 0, None)) - x1 = delta_ndvi_bin*indvi + delta_ndvi_bin/2.0 + x1 = delta_ndvi_bin * indvi + delta_ndvi_bin / 2.0 x2 = x1 + delta_ndvi_bin - x[idx] = (ndvi[idx] - x1[idx])/(x2[idx] - x1[idx]) + x[idx] = (ndvi[idx] - x1[idx]) / (x2[idx] - x1[idx]) x[idx] = np.clip(x[idx], 0, 1) indvi = np.array(indvi, dtype=np.int32) @@ -240,54 +289,73 @@ def get_b1_thresholds(data, thresholds, scene_idx): idxdes = np.nonzero((ndvi >= ndvi_bnd1) & (ndvi < ndvi_bnd2) & (ndvi < des_ndvi)) for i in range(3): - y1 = coeff[0, indvi, i, 0] + \ - coeff[0, indvi, i, 1]*sctang + \ - coeff[0, indvi, i, 2]*sctang**2 + \ - coeff[0, indvi, i, 3]*sctang**3 - y1[des] = coeff[1, indvi[des], i, 0] + \ - coeff[1, indvi[des], i, 1]*sctang[des] + \ - coeff[1, indvi[des], i, 2]*sctang[des]**2 + \ - coeff[1, indvi[des], i, 3]*sctang[des]**3 - - y2[idx] = coeff[0, indvi[idx], i, 0] + \ - coeff[0, indvi[idx], i, 1]*sctang[idx] + \ - coeff[0, indvi[idx], i, 2]*sctang[idx]**2 + \ - coeff[0, indvi[idx], i, 3]*sctang[idx]**3 - y2[idxdes] = coeff[1, indvi[idxdes], i, 0] + \ - coeff[1, indvi[idxdes], i, 1]*sctang[idxdes] + \ - coeff[1, indvi[idxdes], i, 2]*sctang[idxdes]**2 + \ - coeff[1, indvi[idxdes], i, 3]*sctang[idxdes]**3 + y1 = ( + coeff[0, indvi, i, 0] + + coeff[0, indvi, i, 1] * sctang + + coeff[0, indvi, i, 2] * sctang**2 + + coeff[0, indvi, i, 3] * sctang**3 + ) + y1[des] = ( + coeff[1, indvi[des], i, 0] + + coeff[1, indvi[des], i, 1] * sctang[des] + + coeff[1, indvi[des], i, 2] * sctang[des] ** 2 + + coeff[1, indvi[des], i, 3] * sctang[des] ** 3 + ) + + y2[idx] = ( + coeff[0, indvi[idx], i, 0] + + coeff[0, indvi[idx], i, 1] * sctang[idx] + + coeff[0, indvi[idx], i, 2] * sctang[idx] ** 2 + + coeff[0, indvi[idx], i, 3] * sctang[idx] ** 3 + ) + y2[idxdes] = ( + coeff[1, indvi[idxdes], i, 0] + + coeff[1, indvi[idxdes], i, 1] * sctang[idxdes] + + coeff[1, indvi[idxdes], i, 2] * sctang[idxdes] ** 2 + + coeff[1, indvi[idxdes], i, 3] * sctang[idxdes] ** 3 + ) thr[:, i] = ((1.0 - x) * y1) + (x * y2) thr_adj[ndvi < des_ndvi, i] = thr[ndvi < des_ndvi, i] * thr_adj_fac_desert thr_adj[ndvi >= des_ndvi, i] = thr[ndvi >= des_ndvi, i] * thr_adj_fac_land - hicut = ((thr[:, 0] + thr_adj[:, 0])/100) - midpt = ((thr[:, 1] + thr_adj[:, 1])/100) - locut = ((thr[:, 2] + thr_adj[:, 2])/100) + hicut = (thr[:, 0] + thr_adj[:, 0]) / 100 + midpt = (thr[:, 1] + thr_adj[:, 1]) / 100 + locut = (thr[:, 2] + thr_adj[:, 2]) / 100 idx = np.nonzero((ndvi >= fill_ndvi[0]) | (ndvi <= fill_ndvi[1])) hicut[idx] = -999 midpt[idx] = -999 locut[idx] = -999 -# out_thr = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(data.ndvi.shape), -# np.full(data.ndvi.shape, 2))), -# dims=('number_of_lines', 'number_of_pixels', 'z')) -# -# return out_thr + # out_thr = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(data.ndvi.shape), + # np.full(data.ndvi.shape, 2))), + # dims=('number_of_lines', 'number_of_pixels', 'z')) + # + # return out_thr return locut, midpt, hicut def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx): - thresholds = thresholds[scene] - if ((test_name == '4-12um_BTD_Thin_Cirrus_Test') and (scene in ['Land_Night', 'Night_Snow']) or - (test_name == '7.3-11um_BTD_Mid_Level_Cloud_Test') and (scene == 'Land_Night')): - locut = thresholds[test_name]['thr'][0] * np.ones(data.M15.values[scene_idx].ravel()) - midpt = thresholds[test_name]['thr'][1] * np.ones(data.M15.values[scene_idx].ravel()) - hicut = thresholds[test_name]['thr'][2] * np.ones(data.M15.values[scene_idx].ravel()) - power = thresholds[test_name]['thr'][3] * np.ones(data.M15.values[scene_idx].ravel()) + if ( + (test_name == "4-12um_BTD_Thin_Cirrus_Test") + and (scene in ["Land_Night", "Night_Snow"]) + or (test_name == "7.3-11um_BTD_Mid_Level_Cloud_Test") + and (scene == "Land_Night") + ): + locut = thresholds[test_name]["thr"][0] * np.ones( + data.M15.values[scene_idx].ravel() + ) + midpt = thresholds[test_name]["thr"][1] * np.ones( + data.M15.values[scene_idx].ravel() + ) + hicut = thresholds[test_name]["thr"][2] * np.ones( + data.M15.values[scene_idx].ravel() + ) + power = thresholds[test_name]["thr"][3] * np.ones( + data.M15.values[scene_idx].ravel() + ) # out_thr = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(data.ndvi.shape), power)), # dims=('number_of_lines', 'number_of_pixels', 'z')) @@ -296,7 +364,7 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx): return out_thr.T rad = data.M15.values[scene_idx].ravel() - bt_bounds = thresholds[test_name]['bt11_bounds'] + bt_bounds = thresholds[test_name]["bt11_bounds"] locut, midpt = np.empty(rad.shape), np.empty(rad.shape) hicut, power = np.empty(rad.shape), np.empty(rad.shape) @@ -305,56 +373,63 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx): conf_range = np.empty(rad.shape) idx = np.nonzero(rad < bt_bounds[0]) - locut[idx] = thresholds[test_name]['low'][0] - midpt[idx] = thresholds[test_name]['low'][1] - hicut[idx] = thresholds[test_name]['low'][2] - power[idx] = thresholds[test_name]['low'][3] + locut[idx] = thresholds[test_name]["low"][0] + midpt[idx] = thresholds[test_name]["low"][1] + hicut[idx] = thresholds[test_name]["low"][2] + power[idx] = thresholds[test_name]["low"][3] idx = np.nonzero(rad > bt_bounds[3]) - locut[idx] = thresholds[test_name]['high'][0] - midpt[idx] = thresholds[test_name]['high'][1] - hicut[idx] = thresholds[test_name]['high'][2] - power[idx] = thresholds[test_name]['high'][3] + locut[idx] = thresholds[test_name]["high"][0] + midpt[idx] = thresholds[test_name]["high"][1] + hicut[idx] = thresholds[test_name]["high"][2] + power[idx] = thresholds[test_name]["high"][3] # # # # # - idx = np.nonzero((rad >= bt_bounds[0]) & (rad <= bt_bounds[3]) & - (bt_bounds[1] == 0) & (bt_bounds[2] == 0)) - lo[idx] = thresholds[test_name]['bt11_bounds'][0] - hi[idx] = thresholds[test_name]['bt11_bounds'][3] - lo_thr[idx] = thresholds[test_name]['mid1'][0] - hi_thr[idx] = thresholds[test_name]['mid1'][1] - power[idx] = thresholds[test_name]['mid1'][3] - conf_range[idx] = thresholds[test_name]['mid1'][2] + idx = np.nonzero( + (rad >= bt_bounds[0]) + & (rad <= bt_bounds[3]) + & (bt_bounds[1] == 0) + & (bt_bounds[2] == 0) + ) + lo[idx] = thresholds[test_name]["bt11_bounds"][0] + hi[idx] = thresholds[test_name]["bt11_bounds"][3] + lo_thr[idx] = thresholds[test_name]["mid1"][0] + hi_thr[idx] = thresholds[test_name]["mid1"][1] + power[idx] = thresholds[test_name]["mid1"][3] + conf_range[idx] = thresholds[test_name]["mid1"][2] idx = np.nonzero((rad >= bt_bounds[0]) & (rad < bt_bounds[1])) - lo[idx] = thresholds[test_name]['bt11_bounds'][0] - hi[idx] = thresholds[test_name]['bt11_bounds'][1] - lo_thr[idx] = thresholds[test_name]['mid1'][0] - hi_thr[idx] = thresholds[test_name]['mid1'][1] - power[idx] = thresholds[test_name]['mid1'][3] - conf_range[idx] = thresholds[test_name]['mid1'][2] + lo[idx] = thresholds[test_name]["bt11_bounds"][0] + hi[idx] = thresholds[test_name]["bt11_bounds"][1] + lo_thr[idx] = thresholds[test_name]["mid1"][0] + hi_thr[idx] = thresholds[test_name]["mid1"][1] + power[idx] = thresholds[test_name]["mid1"][3] + conf_range[idx] = thresholds[test_name]["mid1"][2] idx = np.nonzero((rad >= bt_bounds[1]) & (rad < bt_bounds[2])) - lo[idx] = thresholds[test_name]['bt11_bounds'][1] - hi[idx] = thresholds[test_name]['bt11_bounds'][2] - lo_thr[idx] = thresholds[test_name]['mid2'][0] - hi_thr[idx] = thresholds[test_name]['mid2'][1] - power[idx] = thresholds[test_name]['mid2'][3] - conf_range[idx] = thresholds[test_name]['mid2'][2] + lo[idx] = thresholds[test_name]["bt11_bounds"][1] + hi[idx] = thresholds[test_name]["bt11_bounds"][2] + lo_thr[idx] = thresholds[test_name]["mid2"][0] + hi_thr[idx] = thresholds[test_name]["mid2"][1] + power[idx] = thresholds[test_name]["mid2"][3] + conf_range[idx] = thresholds[test_name]["mid2"][2] idx = np.nonzero((rad >= bt_bounds[2]) & (rad < bt_bounds[3])) - lo[idx] = thresholds[test_name]['bt11_bounds'][2] - hi[idx] = thresholds[test_name]['bt11_bounds'][3] - lo_thr[idx] = thresholds[test_name]['mid3'][0] - hi_thr[idx] = thresholds[test_name]['mid3'][1] - power[idx] = thresholds[test_name]['mid3'][3] - conf_range[idx] = thresholds[test_name]['mid3'][2] - - idx = np.nonzero(((rad >= bt_bounds[0]) & (rad < bt_bounds[3])) | - (bt_bounds[1] == 0.0) | (bt_bounds[2] == 0)) - - a = (rad[idx] - lo[idx])/(hi[idx] - lo[idx]) - midpt[idx] = lo_thr[idx] + (a*(hi_thr[idx] - lo_thr[idx])) + lo[idx] = thresholds[test_name]["bt11_bounds"][2] + hi[idx] = thresholds[test_name]["bt11_bounds"][3] + lo_thr[idx] = thresholds[test_name]["mid3"][0] + hi_thr[idx] = thresholds[test_name]["mid3"][1] + power[idx] = thresholds[test_name]["mid3"][3] + conf_range[idx] = thresholds[test_name]["mid3"][2] + + idx = np.nonzero( + ((rad >= bt_bounds[0]) & (rad < bt_bounds[3])) + | (bt_bounds[1] == 0.0) + | (bt_bounds[2] == 0) + ) + + a = (rad[idx] - lo[idx]) / (hi[idx] - lo[idx]) + midpt[idx] = lo_thr[idx] + (a * (hi_thr[idx] - lo_thr[idx])) hicut[idx] = midpt[idx] - conf_range[idx] locut[idx] = midpt[idx] + conf_range[idx] @@ -373,40 +448,52 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx): # get_nl_thresholds def land_night_thresholds(data, threshold, coast=True): if coast is False: - lo_val = threshold['bt_diff_bounds'][0] - hi_val = threshold['bt_diff_bounds'][1] - lo_val_thr = threshold['thr_mid'][0] - hi_val_thr = threshold['thr_mid'][1] - conf_range = threshold['thr_mid'][2] - - a = (data['M15-M16'].values - lo_val) / (hi_val - lo_val) - midpt = lo_val_thr + a*(hi_val_thr - lo_val_thr) + lo_val = threshold["bt_diff_bounds"][0] + hi_val = threshold["bt_diff_bounds"][1] + lo_val_thr = threshold["thr_mid"][0] + hi_val_thr = threshold["thr_mid"][1] + conf_range = threshold["thr_mid"][2] + + a = (data["M15-M16"].values - lo_val) / (hi_val - lo_val) + midpt = lo_val_thr + a * (hi_val_thr - lo_val_thr) hicut = midpt - conf_range locut = midpt + conf_range - power = np.full(midpt.shape, threshold['thr_mid'][3]) - - idx = np.nonzero(data['M15-M16'].values > threshold['bt_diff_bounds'][0]) - locut[idx] = threshold['thr_low'][0] - midpt[idx] = threshold['thr_low'][1] - hicut[idx] = threshold['thr_low'][2] - power[idx] = threshold['thr_low'][3] - - idx = np.nonzero(data['M15-M16'].values < threshold['bt_diff_bounds'][1]) - locut[idx] = threshold['thr_hi'][0] - midpt[idx] = threshold['thr_hi'][1] - hicut[idx] = threshold['thr_hi'][2] - power[idx] = threshold['thr_hi'][3] - - out_thr = np.dstack((locut.ravel(), midpt.ravel(), hicut.ravel(), - np.ones(locut.ravel().shape), power.ravel())) + power = np.full(midpt.shape, threshold["thr_mid"][3]) + + idx = np.nonzero(data["M15-M16"].values > threshold["bt_diff_bounds"][0]) + locut[idx] = threshold["thr_low"][0] + midpt[idx] = threshold["thr_low"][1] + hicut[idx] = threshold["thr_low"][2] + power[idx] = threshold["thr_low"][3] + + idx = np.nonzero(data["M15-M16"].values < threshold["bt_diff_bounds"][1]) + locut[idx] = threshold["thr_hi"][0] + midpt[idx] = threshold["thr_hi"][1] + hicut[idx] = threshold["thr_hi"][2] + power[idx] = threshold["thr_hi"][3] + + out_thr = np.dstack( + ( + locut.ravel(), + midpt.ravel(), + hicut.ravel(), + np.ones(locut.ravel().shape), + power.ravel(), + ) + ) return np.squeeze(out_thr.T) else: - b0 = threshold['coeffs'][0] - b1 = threshold['coeffs'][1] - b2 = threshold['coeffs'][2] - - thr = b0 + threshold['int_adj'] + b1*data.geos_tpw.values + b2*data.geos_tpw.values**2 + b0 = threshold["coeffs"][0] + b1 = threshold["coeffs"][1] + b2 = threshold["coeffs"][2] + + thr = ( + b0 + + threshold["int_adj"] + + b1 * data.geos_tpw.values + + b2 * data.geos_tpw.values**2 + ) hicut = np.empty((2, thr.shape[0], thr.shape[1])) midpt = np.empty((2, thr.shape[0], thr.shape[1])) locut = np.empty((2, thr.shape[0], thr.shape[1])) @@ -417,15 +504,24 @@ def land_night_thresholds(data, threshold, coast=True): locut[0, :] = hicut[0, :] - 1 locut[1, :] = hicut[1, :] + 1 - out_thr = np.dstack((hicut[0, :], midpt[0, :], locut[0, :], locut[1, :], midpt[1, :], hicut[1, :], - np.ones(locut[0, :].shape), np.ones(locut[0, :].shape))) + out_thr = np.dstack( + ( + hicut[0, :], + midpt[0, :], + locut[0, :], + locut[1, :], + midpt[1, :], + hicut[1, :], + np.ones(locut[0, :].shape), + np.ones(locut[0, :].shape), + ) + ) return np.moveaxis(out_thr, -1, 0) def vis_refl_thresholds(data, thresholds, scene, scene_idx): - locut, midpt, hicut = get_b1_thresholds(data, thresholds, scene_idx) - bias_adj = thresholds[scene]['Visible_Reflectance_Test']['adj'] + bias_adj = thresholds[scene]["Visible_Reflectance_Test"]["adj"] ndvi = data.ndvi.values[scene_idx].ravel() m01 = data.M05.values[scene_idx].ravel() m02 = data.M07.values[scene_idx].ravel() @@ -437,8 +533,8 @@ def vis_refl_thresholds(data, thresholds, scene, scene_idx): b1_midpt = midpt * bias_adj b1_hicut = hicut * bias_adj - if scene in ['Land_Day_Desert', 'Land_Day_Desert_Coast']: - ndvi_desert_thr = thresholds[scene]['Visible_Reflectance_Test']['ndvi_thr'] + if scene in ["Land_Day_Desert", "Land_Day_Desert_Coast"]: + ndvi_desert_thr = thresholds[scene]["Visible_Reflectance_Test"]["ndvi_thr"] idx = np.nonzero(ndvi < ndvi_desert_thr) b1_locut[idx] = locut[idx] @@ -449,19 +545,19 @@ def vis_refl_thresholds(data, thresholds, scene, scene_idx): b1_power = np.full(b1_locut.shape, 2) idx = np.nonzero(locut == -999) - b1_locut[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][0] - b1_midpt[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][1] - b1_hicut[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][2] - b1_power[idx] = thresholds[scene]['Visible_Reflectance_Test']['thr'][3] - if scene in ['Land_Day_Desert', 'Land_Day_Desert_Coast']: + b1_locut[idx] = thresholds[scene]["Visible_Reflectance_Test"]["thr"][0] + b1_midpt[idx] = thresholds[scene]["Visible_Reflectance_Test"]["thr"][1] + b1_hicut[idx] = thresholds[scene]["Visible_Reflectance_Test"]["thr"][2] + b1_power[idx] = thresholds[scene]["Visible_Reflectance_Test"]["thr"][3] + if scene in ["Land_Day_Desert", "Land_Day_Desert_Coast"]: m128[idx] = m02[idx] cosvza = np.cos(data.sensor_zenith.values[scene_idx] * _DTR).ravel() - vzcpow = thresholds['VZA_correction']['vzcpow'][0] - b1_locut = (b1_locut * (1.0 / np.power(cosvza, vzcpow))) - b1_midpt = (b1_midpt * (1.0 / np.power(cosvza, vzcpow))) - b1_hicut = (b1_hicut * (1.0 / np.power(cosvza, vzcpow))) + vzcpow = thresholds["VZA_correction"]["vzcpow"][0] + b1_locut = b1_locut * (1.0 / np.power(cosvza, vzcpow)) + b1_midpt = b1_midpt * (1.0 / np.power(cosvza, vzcpow)) + b1_hicut = b1_hicut * (1.0 / np.power(cosvza, vzcpow)) # out_thr = xr.DataArray(data=np.dstack((b1_locut, b1_midpt, b1_hicut, np.ones(data.ndvi.shape), # b1_power.reshape(data.ndvi.shape))), @@ -477,11 +573,11 @@ def gemi_thresholds(data, thresholds, scene_name, scene_idx): gemi_thr = np.ones((ndvi.shape[0], 5)) idx = np.nonzero(ndvi < 0.1) - gemi_thr[idx, :3] = thresholds['gemi0'][:3] + gemi_thr[idx, :3] = thresholds["gemi0"][:3] idx = np.nonzero((ndvi >= 0.1) & (ndvi < 0.2)) - gemi_thr[idx, :3] = thresholds['gemi1'][:3] + gemi_thr[idx, :3] = thresholds["gemi1"][:3] idx = np.nonzero((ndvi >= 0.2) & (ndvi < 0.3)) - gemi_thr[idx, :3] = thresholds['gemi2'][:3] + gemi_thr[idx, :3] = thresholds["gemi2"][:3] # thr_out = xr.DataArray(data=np.dstack((gemi_thr[:, :, 0], gemi_thr[:, :, 1], gemi_thr[:, :, 2], # np.ones(gemi_thr[:, :, 0].shape), @@ -492,49 +588,70 @@ def gemi_thresholds(data, thresholds, scene_name, scene_idx): def bt_diff_11_4um_thresholds(data, threshold, scene_idx): - c = threshold['coeffs'] + c = threshold["coeffs"] tpw = data.geos_tpw.values[scene_idx].ravel() - thr = c[0] + threshold['corr'] + c[1]*tpw + c[2]*np.power(tpw, 2) - - hicut0 = thr + threshold['hicut_coeff'][0] - hicut1 = thr + threshold['hicut_coeff'][1] - midpt0 = hicut0 + threshold['midpt_coeff'][0] - midpt1 = hicut1 + threshold['midpt_coeff'][1] - locut0 = hicut0 + threshold['locut_coeff'][0] - locut1 = hicut1 + threshold['locut_coeff'][1] - - thr_out = np.dstack([locut0, midpt0, hicut0, hicut1, midpt1, locut1, - np.ones(hicut0.shape), np.ones(hicut0.shape)]) + thr = c[0] + threshold["corr"] + c[1] * tpw + c[2] * np.power(tpw, 2) + + hicut0 = thr + threshold["hicut_coeff"][0] + hicut1 = thr + threshold["hicut_coeff"][1] + midpt0 = hicut0 + threshold["midpt_coeff"][0] + midpt1 = hicut1 + threshold["midpt_coeff"][1] + locut0 = hicut0 + threshold["locut_coeff"][0] + locut1 = hicut1 + threshold["locut_coeff"][1] + + thr_out = np.dstack( + [ + locut0, + midpt0, + hicut0, + hicut1, + midpt1, + locut1, + np.ones(hicut0.shape), + np.ones(hicut0.shape), + ] + ) return np.squeeze(thr_out.T) def thresholds_1_38um_test(data, thresholds, scene_name, scene_idx): sza = data.solar_zenith.values[scene_idx] vza = data.sensor_zenith.values[scene_idx] - thresh = thresholds[scene_name]['1.38um_High_Cloud_Test'] - c = thresh['coeffs'] - vzcpow = thresholds['VZA_correction']['vzcpow'][1] - - hicut = c[0] + c[1]*sza + c[2]*np.power(sza, 2) + c[3]*np.power(sza, 3) + \ - c[4]*np.power(sza, 4) + c[5]*np.power(sza, 5) - hicut = hicut*0.01 + (np.maximum((sza/90.)*thresh['szafac']*thresh['adj'], thresh['adj'])) + thresh = thresholds[scene_name]["1.38um_High_Cloud_Test"] + c = thresh["coeffs"] + vzcpow = thresholds["VZA_correction"]["vzcpow"][1] + + hicut = ( + c[0] + + c[1] * sza + + c[2] * np.power(sza, 2) + + c[3] * np.power(sza, 3) + + c[4] * np.power(sza, 4) + + c[5] * np.power(sza, 5) + ) + hicut = hicut * 0.01 + ( + np.maximum((sza / 90.0) * thresh["szafac"] * thresh["adj"], thresh["adj"]) + ) midpt = hicut + 0.001 locut = midpt + 0.001 - cosvza = np.cos(vza*_DTR) + cosvza = np.cos(vza * _DTR) - locut = locut * (1/np.power(cosvza, vzcpow)) - midpt = midpt * (1/np.power(cosvza, vzcpow)) - hicut = hicut * (1/np.power(cosvza, vzcpow)) + locut = locut * (1 / np.power(cosvza, vzcpow)) + midpt = midpt * (1 / np.power(cosvza, vzcpow)) + hicut = hicut * (1 / np.power(cosvza, vzcpow)) # out_thr = xr.DataArray(data=np.dstack((locut, midpt, hicut, np.ones(data.ndvi.shape), # np.ones(data.ndvi.shape))), # dims=('number_of_lines', 'number_of_pixels', 'z')) - out_thr = np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape))) + out_thr = np.dstack( + (locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape)) + ) return np.squeeze(out_thr.T) + # NOTE: 11-12um Cirrus Test # hicut is computed in different ways depending on the scene # 1. midpt - adj @@ -579,4 +696,3 @@ def thresholds_1_38um_test(data, thresholds, scene_name, scene_idx): # 2. # - Ocean_Day # - Polar_Ocean_Day - diff --git a/mvcm/read_data.py b/mvcm/read_data.py index a67c2ff..fdd3157 100644 --- a/mvcm/read_data.py +++ b/mvcm/read_data.py @@ -1,23 +1,42 @@ -import xarray as xr +import logging +import os +from datetime import datetime as dt +from typing import Dict + import numpy as np import numpy.typing as npt +import xarray as xr +from attrs import Factory, define, field, validators import ancillary_data as anc import mvcm.scene as scn -import os -import logging -from datetime import datetime as dt -from typing import Dict -from attrs import define, field, validators, Factory - -_DTR = np.pi/180. -_RTD = 180./np.pi +_DTR = np.pi / 180.0 +_RTD = 180.0 / np.pi _bad_data = -999.0 -_datapath = '/ships19/hercules/pveglio/mvcm_test_data' -_band_names = ['um0_41', 'um0_44', 'um0_49', 'um0_55', 'um0_64h', 'um0_67', 'um0_74', 'um0_86', - 'um0_86h', 'um1_24', 'um1_37', 'um1_61', 'um2_25', 'um3_70', 'um3_70h', 'um4_05', - 'um8_55', 'um10_76', 'um11_45h', 'um12_01'] +_datapath = "/ships19/hercules/pveglio/mvcm_test_data" +_band_names = [ + "um0_41", + "um0_44", + "um0_49", + "um0_55", + "um0_64h", + "um0_67", + "um0_74", + "um0_86", + "um0_86h", + "um1_24", + "um1_37", + "um1_61", + "um2_25", + "um3_70", + "um3_70h", + "um4_05", + "um8_55", + "um10_76", + "um11_45h", + "um12_01", +] logger = logging.getLogger(__name__) # logging.basicConfig(level=logging.INFO, format='%(name)s - %(levelname)s - %(message)s') # logging.basicConfig(level=logging.INFO, filename='logfile.log', 'filemode='w', @@ -54,31 +73,80 @@ class CollectInputs(object): dims [optional]: str name of the dimensions for the arrays in the xarray.Dataset output """ - file_name_geo: str = field(default=f'{_datapath}/VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc', - validator=[validators.instance_of(str), ]) - file_name_l1b: str = field(default=f'{_datapath}/VNP02MOD.A2022173.1312.001.2022174011547.uwssec.nc', - validator=[validators.instance_of(str), ]) - ancillary_dir: str = field(default=f'{_datapath}/ancillary', - validator=[validators.instance_of(str), ]) - sst_file: str = field(default='oisst.20220622', - validator=[validators.instance_of(str), ]) - eco_file: str = field(default='goge1_2_img.v1', - validator=[validators.instance_of(str), ]) - ndvi_file: str = field(default='NDVI.FM.c004.v2.0.WS.00-04.177.hdf', - validator=[validators.instance_of(str), ]) - geos_file_1: str = field(default='GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4', - validator=[validators.instance_of(str), ]) - geos_file_2: str = field(default='GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4', - validator=[validators.instance_of(str), ]) - geos_land: str = field(default='GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4', - validator=[validators.instance_of(str), ]) - geos_ocean: str = field(default='GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4', - validator=[validators.instance_of(str), ]) - geos_constants: str = field(default='GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4', - validator=[validators.instance_of(str), ]) - - dims: tuple = field(default=('number_of_lines', 'number_of_pixels'), - validator=[validators.instance_of(tuple), ]) + + file_name_geo: str = field( + default=f"{_datapath}/VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc", + validator=[ + validators.instance_of(str), + ], + ) + file_name_l1b: str = field( + default=f"{_datapath}/VNP02MOD.A2022173.1312.001.2022174011547.uwssec.nc", + validator=[ + validators.instance_of(str), + ], + ) + ancillary_dir: str = field( + default=f"{_datapath}/ancillary", + validator=[ + validators.instance_of(str), + ], + ) + sst_file: str = field( + default="oisst.20220622", + validator=[ + validators.instance_of(str), + ], + ) + eco_file: str = field( + default="goge1_2_img.v1", + validator=[ + validators.instance_of(str), + ], + ) + ndvi_file: str = field( + default="NDVI.FM.c004.v2.0.WS.00-04.177.hdf", + validator=[ + validators.instance_of(str), + ], + ) + geos_file_1: str = field( + default="GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4", + validator=[ + validators.instance_of(str), + ], + ) + geos_file_2: str = field( + default="GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4", + validator=[ + validators.instance_of(str), + ], + ) + geos_land: str = field( + default="GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4", + validator=[ + validators.instance_of(str), + ], + ) + geos_ocean: str = field( + default="GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4", + validator=[ + validators.instance_of(str), + ], + ) + geos_constants: str = field( + default="GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4", + validator=[ + validators.instance_of(str), + ], + ) + + dims: tuple = field( + default=("number_of_lines", "number_of_pixels"), + validator=[ + validators.instance_of(tuple), + ], + ) @define(slots=True, kw_only=True) @@ -94,14 +162,16 @@ class ReadData(CollectInputs): hires: bool high resolution data flag """ - satellite: str = field(validator=[validators.instance_of(str), - validators.in_(['snpp'])]) - sensor: str = field(validator=[validators.instance_of(str), - validators.in_(['viirs'])]) - hires: bool = field(default=False, - validator=[validators.instance_of(bool)]) - logger.debug('Instance of ReadData created') + satellite: str = field( + validator=[validators.instance_of(str), validators.in_(["snpp"])] + ) + sensor: str = field( + validator=[validators.instance_of(str), validators.in_(["viirs"])] + ) + hires: bool = field(default=False, validator=[validators.instance_of(bool)]) + + logger.debug("Instance of ReadData created") def read_viirs_geo(self) -> xr.Dataset: """Read VIIRS geolocation data and generate additional angles @@ -116,24 +186,28 @@ class ReadData(CollectInputs): geo_data xarray.Dataset dataset containing all geolocation data """ - logger.debug(f'Reading {self.file_name_geo}') - geo_data = xr.open_dataset(self.file_name_geo, group='geolocation_data') + logger.debug(f"Reading {self.file_name_geo}") + geo_data = xr.open_dataset(self.file_name_geo, group="geolocation_data") - relazi = self.relative_azimuth_angle(geo_data.sensor_azimuth.values, geo_data.solar_azimuth.values) - sunglint = self.sun_glint_angle(geo_data.sensor_zenith.values, geo_data.solar_zenith.values, relazi) - scatt_angle = self.scattering_angle(geo_data.solar_zenith.values, geo_data.sensor_zenith.values, - relazi) + relazi = self.relative_azimuth_angle( + geo_data.sensor_azimuth.values, geo_data.solar_azimuth.values + ) + sunglint = self.sun_glint_angle( + geo_data.sensor_zenith.values, geo_data.solar_zenith.values, relazi + ) + scatt_angle = self.scattering_angle( + geo_data.solar_zenith.values, geo_data.sensor_zenith.values, relazi + ) - geo_data['relative_azimuth'] = (self.dims, relazi) - geo_data['sunglint_angle'] = (self.dims, sunglint) - geo_data['scattering_angle'] = (self.dims, scatt_angle) + geo_data["relative_azimuth"] = (self.dims, relazi) + geo_data["sunglint_angle"] = (self.dims, sunglint) + geo_data["scattering_angle"] = (self.dims, scatt_angle) - logger.debug('Geolocation file read correctly') + logger.debug("Geolocation file read correctly") return geo_data - def read_viirs_l1b(self, - solar_zenith: npt.NDArray[float]) -> xr.Dataset: + def read_viirs_l1b(self, solar_zenith: npt.NDArray[float]) -> xr.Dataset: """Read VIIRS L1b data Parameters @@ -143,33 +217,42 @@ class ReadData(CollectInputs): solar_zenith: np.ndarray solar zenith angle derived from the geolocation file """ - logger.debug(f'Reading {self.file_name_l1b}') - l1b_data = xr.open_dataset(self.file_name_l1b, group='observation_data', decode_cf=False) + logger.debug(f"Reading {self.file_name_l1b}") + l1b_data = xr.open_dataset( + self.file_name_l1b, group="observation_data", decode_cf=False + ) rad_data = xr.Dataset() for band in list(l1b_data.variables): - if 'reflectance' in l1b_data[band].long_name: - if hasattr(l1b_data[band], 'VCST_scale_factor'): - scale_factor = l1b_data[band].VCST_scale_factor * l1b_data[band].bias_correction + if "reflectance" in l1b_data[band].long_name: + if hasattr(l1b_data[band], "VCST_scale_factor"): + scale_factor = ( + l1b_data[band].VCST_scale_factor + * l1b_data[band].bias_correction + ) else: scale_factor = l1b_data[band].scale_factor - rad_data[band] = (self.dims, l1b_data[band].values * scale_factor / np.cos(solar_zenith*_DTR)) - - elif 'radiance' in l1b_data[band].long_name: - bt_lut = f'{band}_brightness_temperature_lut' - rad_data[band] = (self.dims, - l1b_data[bt_lut].values[l1b_data[band].values]) + rad_data[band] = ( + self.dims, + l1b_data[band].values * scale_factor / np.cos(solar_zenith * _DTR), + ) + + elif "radiance" in l1b_data[band].long_name: + bt_lut = f"{band}_brightness_temperature_lut" + rad_data[band] = ( + self.dims, + l1b_data[bt_lut].values[l1b_data[band].values], + ) else: pass - logger.debug('L1b file read correctly') + logger.debug("L1b file read correctly") return rad_data - def preprocess_viirs(self, - geo_data: xr.Dataset, - viirs: xr.Dataset, - hires_data: xr.Dataset = None) -> xr.Dataset: + def preprocess_viirs( + self, geo_data: xr.Dataset, viirs: xr.Dataset, hires_data: xr.Dataset = None + ) -> xr.Dataset: """Create band combinations (e.g. ratios, differences) that are used by some spectral tests. Parameters @@ -192,21 +275,41 @@ class ReadData(CollectInputs): if hires_data is not None: viirs_out = xr.Dataset() - mod_bands = ['M01', 'M02', 'M03', 'M04', 'M06', 'M08', 'M09', 'M11', 'M13', 'M14', 'M15', 'M16'] + mod_bands = [ + "M01", + "M02", + "M03", + "M04", + "M06", + "M08", + "M09", + "M11", + "M13", + "M14", + "M15", + "M16", + ] for bnd in mod_bands: - viirs_out[bnd] = (self.dims, viirs[bnd].values.repeat(2, 0).repeat(2, 1)) + viirs_out[bnd] = ( + self.dims, + viirs[bnd].values.repeat(2, 0).repeat(2, 1), + ) - if (('I01' in hires_data) and ('I02' in hires_data)): - viirs_out['M05'] = (self.dims, hires_data.I01.values) - viirs_out['M07'] = (self.dims, hires_data.I02.values) + if ("I01" in hires_data) and ("I02" in hires_data): + viirs_out["M05"] = (self.dims, hires_data.I01.values) + viirs_out["M07"] = (self.dims, hires_data.I02.values) - if (('M05' in viirs) and ('M07' in viirs)): + if ("M05" in viirs) and ("M07" in viirs): m01 = viirs_out.M05.values m02 = viirs_out.M07.values - r1 = 2.0 * (np.power(m02, 2.0) - np.power(m01, 2.0)) + (1.5 * m02) + (0.5 * m01) + r1 = ( + 2.0 * (np.power(m02, 2.0) - np.power(m01, 2.0)) + + (1.5 * m02) + + (0.5 * m01) + ) r2 = m01 + m02 + 0.5 r3 = r1 / r2 - gemi = r3 * (1.0 - 0.25*r3) - ((m01 - 0.125) / (1.0 - m01)) + gemi = r3 * (1.0 - 0.25 * r3) - ((m01 - 0.125) / (1.0 - m01)) else: gemi = np.full((viirs_out.M16.shape), _bad_data) @@ -235,37 +338,49 @@ class ReadData(CollectInputs): """ # Compute channel differences and ratios that are used in the tests - if (('M05' in viirs) and ('M07' in viirs)): - viirs_out['M07-M05ratio'] = (self.dims, viirs_out.M07.values / viirs_out.M05.values) + if ("M05" in viirs) and ("M07" in viirs): + viirs_out["M07-M05ratio"] = ( + self.dims, + viirs_out.M07.values / viirs_out.M05.values, + ) else: - viirs_out['M07-M05ratio'] = (self.dims, np.full(viirs_out.M16.shape, _bad_data)) - viirs_out['M13-M16'] = (self.dims, viirs_out.M13.values - viirs_out.M16.values) - viirs_out['M14-M15'] = (self.dims, viirs_out.M14.values - viirs_out.M15.values) + viirs_out["M07-M05ratio"] = ( + self.dims, + np.full(viirs_out.M16.shape, _bad_data), + ) + viirs_out["M13-M16"] = (self.dims, viirs_out.M13.values - viirs_out.M16.values) + viirs_out["M14-M15"] = (self.dims, viirs_out.M14.values - viirs_out.M15.values) if hires_data is not None: - viirs_out['M15-M12'] = (self.dims, hires_data.I05.values - hires_data.I04.values) + viirs_out["M15-M12"] = ( + self.dims, + hires_data.I05.values - hires_data.I04.values, + ) else: - viirs_out['M15-M12'] = (self.dims, viirs_out.M15.values - viirs_out.M12.values) - viirs_out['M15-M13'] = (self.dims, viirs_out.M15.values - viirs_out.M13.values) - viirs_out['M15-M16'] = (self.dims, viirs_out.M15.values - viirs_out.M16.values) - viirs_out['GEMI'] = (self.dims, gemi) + viirs_out["M15-M12"] = ( + self.dims, + viirs_out.M15.values - viirs_out.M12.values, + ) + viirs_out["M15-M13"] = (self.dims, viirs_out.M15.values - viirs_out.M13.values) + viirs_out["M15-M16"] = (self.dims, viirs_out.M15.values - viirs_out.M16.values) + viirs_out["GEMI"] = (self.dims, gemi) if hires_data is not None: - viirs_out['M10'] = (self.dims, hires_data.I03.values) - viirs_out['M12'] = (self.dims, hires_data.I04.values) - viirs_out['M15hi'] = (self.dims, hires_data.I05.values) + viirs_out["M10"] = (self.dims, hires_data.I03.values) + viirs_out["M12"] = (self.dims, hires_data.I04.values) + viirs_out["M15hi"] = (self.dims, hires_data.I05.values) # temp value to force the code to work - viirs_out['M128'] = (self.dims, np.zeros(viirs_out.M16.shape)) + viirs_out["M128"] = (self.dims, np.zeros(viirs_out.M16.shape)) viirs_out.update(geo_data) - viirs_out = viirs_out.set_coords(['latitude', 'longitude']) + viirs_out = viirs_out.set_coords(["latitude", "longitude"]) - logger.debug('Viirs preprocessing completed successfully.') + logger.debug("Viirs preprocessing completed successfully.") return viirs_out - def relative_azimuth_angle(self, - sensor_azimuth: npt.NDArray[float], - solar_azimuth: npt.NDArray[float]) -> npt.NDArray[float]: + def relative_azimuth_angle( + self, sensor_azimuth: npt.NDArray[float], solar_azimuth: npt.NDArray[float] + ) -> npt.NDArray[float]: """Computation of the relative azimuth angle Parameters @@ -279,16 +394,18 @@ class ReadData(CollectInputs): ------- relative_azimuth: np.ndarray """ - rel_azimuth = np.abs(180. - np.abs(sensor_azimuth - solar_azimuth)) + rel_azimuth = np.abs(180.0 - np.abs(sensor_azimuth - solar_azimuth)) - logger.debug('Relative azimuth calculated successfully.') + logger.debug("Relative azimuth calculated successfully.") return rel_azimuth - def sun_glint_angle(self, - sensor_zenith: npt.NDArray[float], - solar_zenith: npt.NDArray[float], - rel_azimuth: npt.NDArray[float]) -> npt.NDArray[float]: + def sun_glint_angle( + self, + sensor_zenith: npt.NDArray[float], + solar_zenith: npt.NDArray[float], + rel_azimuth: npt.NDArray[float], + ) -> npt.NDArray[float]: """Computation of the sun glint angle Parameters @@ -304,19 +421,22 @@ class ReadData(CollectInputs): ------- sunglint_angle: np.ndarray """ - cossna = (np.sin(sensor_zenith*_DTR) * np.sin(solar_zenith*_DTR) * np.cos(rel_azimuth*_DTR) + - np.cos(sensor_zenith*_DTR) * np.cos(solar_zenith*_DTR)) + cossna = np.sin(sensor_zenith * _DTR) * np.sin(solar_zenith * _DTR) * np.cos( + rel_azimuth * _DTR + ) + np.cos(sensor_zenith * _DTR) * np.cos(solar_zenith * _DTR) cossna[cossna > 1] = 1 sunglint_angle = np.arccos(cossna) * _RTD - logger.debug('Sunglint generated') + logger.debug("Sunglint generated") return sunglint_angle - def scattering_angle(self, - solar_zenith: npt.NDArray[float], - sensor_zenith: npt.NDArray[float], - relative_azimuth: npt.NDArray[float]) -> npt.NDArray[float]: + def scattering_angle( + self, + solar_zenith: npt.NDArray[float], + sensor_zenith: npt.NDArray[float], + relative_azimuth: npt.NDArray[float], + ) -> npt.NDArray[float]: """Computation of the scattering angle Parameters @@ -332,13 +452,16 @@ class ReadData(CollectInputs): ------- scattering_angle: np.ndarray """ - cos_scatt_angle = -1. * (np.cos(solar_zenith*_DTR) * np.cos(sensor_zenith*_DTR) - - np.sin(solar_zenith*_DTR) * np.sin(sensor_zenith*_DTR) * - np.cos(relative_azimuth*_DTR)) + cos_scatt_angle = -1.0 * ( + np.cos(solar_zenith * _DTR) * np.cos(sensor_zenith * _DTR) + - np.sin(solar_zenith * _DTR) + * np.sin(sensor_zenith * _DTR) + * np.cos(relative_azimuth * _DTR) + ) scatt_angle = np.arccos(cos_scatt_angle) * _RTD - logger.debug('Scattering angle calculated successfully') + logger.debug("Scattering angle calculated successfully") return scatt_angle @@ -356,18 +479,32 @@ class ReadAncillary(CollectInputs): resolution: int flag to switch between MOD (1) and IMG (2) resolution """ - latitude: npt.NDArray[float] = field(validator=[validators.instance_of(np.ndarray), ]) - longitude: npt.NDArray[float] = field(validator=[validators.instance_of(np.ndarray), ]) - resolution: int = field(validator=[validators.instance_of(int), - validators.in_([1, 2]), ]) - out_shape: tuple = field(init=False, - default=Factory(lambda self: self.latitude.shape, takes_self=True)) + latitude: npt.NDArray[float] = field( + validator=[ + validators.instance_of(np.ndarray), + ] + ) + longitude: npt.NDArray[float] = field( + validator=[ + validators.instance_of(np.ndarray), + ] + ) + resolution: int = field( + validator=[ + validators.instance_of(int), + validators.in_([1, 2]), + ] + ) + + out_shape: tuple = field( + init=False, default=Factory(lambda self: self.latitude.shape, takes_self=True) + ) def get_granule_time(self): """Get granule timestamp and format it for MVCM""" - vnp_time = '.'.join(os.path.basename(self.file_name_l1b).split('.')[1:3]) - return dt.strftime(dt.strptime(vnp_time, 'A%Y%j.%H%M'), '%Y-%m-%d %H:%M:%S.000') + vnp_time = ".".join(os.path.basename(self.file_name_l1b).split(".")[1:3]) + return dt.strftime(dt.strptime(vnp_time, "A%Y%j.%H%M"), "%Y-%m-%d %H:%M:%S.000") def get_sst(self) -> npt.NDArray[float]: """Read Reynolds SST file @@ -381,11 +518,17 @@ class ReadAncillary(CollectInputs): array containing the Reynolds SST interpolated at the sensor's resolution """ if not os.path.isfile(os.path.join(self.ancillary_dir, self.sst_file)): - logger.error('SST file not found') + logger.error("SST file not found") sst = np.empty(self.out_shape, dtype=np.float32).ravel() - sst = anc.py_get_Reynolds_SST(self.latitude.ravel(), self.longitude.ravel(), self.resolution, - self.ancillary_dir, self.sst_file, sst) - logger.debug('SST file read successfully') + sst = anc.py_get_Reynolds_SST( + self.latitude.ravel(), + self.longitude.ravel(), + self.resolution, + self.ancillary_dir, + self.sst_file, + sst, + ) + logger.debug("SST file read successfully") return sst.reshape(self.out_shape) def get_ndvi(self) -> npt.NDArray[float]: @@ -400,11 +543,17 @@ class ReadAncillary(CollectInputs): NDVI interpolated at the sensor's resolution """ if not os.path.isfile(os.path.join(self.ancillary_dir, self.ndvi_file)): - logger.error('NDVI file not found') + logger.error("NDVI file not found") ndvi = np.empty(self.out_shape, dtype=np.float32).ravel() - ndvi = anc.py_get_NDVI_background(self.latitude.ravel(), self.longitude.ravel(), self.resolution, - self.ancillary_dir, self.ndvi_file, ndvi) - logger.debug('NDVI file read successfully') + ndvi = anc.py_get_NDVI_background( + self.latitude.ravel(), + self.longitude.ravel(), + self.resolution, + self.ancillary_dir, + self.ndvi_file, + ndvi, + ) + logger.debug("NDVI file read successfully") return ndvi.reshape(self.out_shape) def get_eco(self) -> npt.NDArray[float]: @@ -419,9 +568,14 @@ class ReadAncillary(CollectInputs): Olson ecosystem type interpolated at the sensor's resolution """ eco = np.empty(self.out_shape, dtype=np.ubyte).ravel() - eco = anc.py_get_Olson_eco(self.latitude.ravel(), self.longitude.ravel(), self.resolution, - self.ancillary_dir, eco) - logger.debug('Olson ecosystem file read successfully') + eco = anc.py_get_Olson_eco( + self.latitude.ravel(), + self.longitude.ravel(), + self.resolution, + self.ancillary_dir, + eco, + ) + logger.debug("Olson ecosystem file read successfully") return eco.reshape(self.out_shape) def get_geos(self) -> Dict: @@ -436,29 +590,47 @@ class ReadAncillary(CollectInputs): dictionary containing all quantities required by MVCM (see geos_variables here below) """ if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_file_1)): - logger.error('GEOS-5 file 1 not found') + logger.error("GEOS-5 file 1 not found") if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_file_2)): - logger.error('GEOS-5 file 2 not found') + logger.error("GEOS-5 file 2 not found") if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_land)): - logger.error('GEOS-5 land file not found') + logger.error("GEOS-5 land file not found") if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_ocean)): - logger.error('GEOS-5 ocean file not found') + logger.error("GEOS-5 ocean file not found") if not os.path.isfile(os.path.join(self.ancillary_dir, self.geos_constants)): - logger.error('GEOS-5 constants file not found') - - geos_variables = ['tpw', 'snow_fraction', 'ice_fraction', 'ocean_fraction', - 'land_ice_fraction', 'surface_temperature'] - geos_data = {var: np.empty(self.out_shape, dtype=np.float32).ravel() for var in geos_variables} - - geos_data = anc.py_get_GEOS(self.latitude.ravel(), self.longitude.ravel(), self.resolution, - self.get_granule_time(), self.ancillary_dir, self.geos_file_1, - self.geos_file_2, self.geos_land, self.geos_ocean, - self.geos_constants, geos_data) + logger.error("GEOS-5 constants file not found") + + geos_variables = [ + "tpw", + "snow_fraction", + "ice_fraction", + "ocean_fraction", + "land_ice_fraction", + "surface_temperature", + ] + geos_data = { + var: np.empty(self.out_shape, dtype=np.float32).ravel() + for var in geos_variables + } + + geos_data = anc.py_get_GEOS( + self.latitude.ravel(), + self.longitude.ravel(), + self.resolution, + self.get_granule_time(), + self.ancillary_dir, + self.geos_file_1, + self.geos_file_2, + self.geos_land, + self.geos_ocean, + self.geos_constants, + geos_data, + ) for var in list(geos_variables): geos_data[var] = geos_data[var].reshape(self.out_shape) - logger.debug('GEOS data read successfully') + logger.debug("GEOS data read successfully") return geos_data def pack_data(self) -> xr.Dataset: @@ -473,21 +645,20 @@ class ReadAncillary(CollectInputs): dataset containing all the ancillary data """ ancillary_data = xr.Dataset() - ancillary_data['sst'] = (self.dims, self.get_sst()) - ancillary_data['ecosystem'] = (self.dims, self.get_eco()) - ancillary_data['ndvi'] = (self.dims, self.get_ndvi()) + ancillary_data["sst"] = (self.dims, self.get_sst()) + ancillary_data["ecosystem"] = (self.dims, self.get_eco()) + ancillary_data["ndvi"] = (self.dims, self.get_ndvi()) geos_tmp = self.get_geos() for var in list(geos_tmp.keys()): - ancillary_data[f'geos_{var}'] = (self.dims, geos_tmp[var]) + ancillary_data[f"geos_{var}"] = (self.dims, geos_tmp[var]) return ancillary_data -def get_data(satellite: str, - sensor: str, - file_names: Dict, - hires: bool = False) -> xr.Dataset: +def get_data( + satellite: str, sensor: str, file_names: Dict, hires: bool = False +) -> xr.Dataset: """Store every variable in one xarray dataset Parameters @@ -506,42 +677,49 @@ def get_data(satellite: str, input_data: xr.Dataset dataset containing all the satellite and ancillary data """ - Viirs = ReadData(file_name_geo=file_names['MOD03'], - file_name_l1b=file_names['MOD02'], - satellite=satellite, - sensor=sensor) + Viirs = ReadData( + file_name_geo=file_names["MOD03"], + file_name_l1b=file_names["MOD02"], + satellite=satellite, + sensor=sensor, + ) if hires is True: - Viirs_hires = ReadData(file_name_geo=file_names['IMG03'], - file_name_l1b=file_names['IMG02'], - satellite=satellite, - sensor=sensor, - hires=True) + Viirs_hires = ReadData( + file_name_geo=file_names["IMG03"], + file_name_l1b=file_names["IMG02"], + satellite=satellite, + sensor=sensor, + hires=True, + ) geo_data_mod = Viirs.read_viirs_geo() viirs_data_mod = Viirs.read_viirs_l1b(geo_data_mod.solar_zenith.values) if hires is True: geo_data = Viirs_hires.read_viirs_geo() viirs_data_img = Viirs_hires.read_viirs_l1b(geo_data.solar_zenith.values) - viirs_data = Viirs_hires.preprocess_viirs(geo_data, viirs_data_mod, viirs_data_img) + viirs_data = Viirs_hires.preprocess_viirs( + geo_data, viirs_data_mod, viirs_data_img + ) res = 2 else: viirs_data = Viirs.preprocess_viirs(geo_data_mod, viirs_data_mod) geo_data = geo_data_mod res = 1 - Ancillary = ReadAncillary(latitude=geo_data.latitude.values, - longitude=geo_data.longitude.values, - resolution=res, - ancillary_dir=file_names['ANC_DIR'], - geos_file_1=file_names['GEOS_atm_1'], - geos_file_2=file_names['GEOS_atm_2'], - geos_land=file_names['GEOS_land'], - geos_ocean=file_names['GEOS_ocean'], - geos_constants=file_names['GEOS_constants'], - ndvi_file=file_names['NDVI'], - sst_file=file_names['SST'], - eco_file=file_names['ECO'] - ) + Ancillary = ReadAncillary( + latitude=geo_data.latitude.values, + longitude=geo_data.longitude.values, + resolution=res, + ancillary_dir=file_names["ANC_DIR"], + geos_file_1=file_names["GEOS_atm_1"], + geos_file_2=file_names["GEOS_atm_2"], + geos_land=file_names["GEOS_land"], + geos_ocean=file_names["GEOS_ocean"], + geos_constants=file_names["GEOS_constants"], + ndvi_file=file_names["NDVI"], + sst_file=file_names["SST"], + eco_file=file_names["ECO"], + ) viirs_data.update(Ancillary.pack_data()) @@ -550,10 +728,10 @@ def get_data(satellite: str, scene_xr = xr.Dataset() for s in scn._scene_list: - scene_xr[s] = (('number_of_lines', 'number_of_pixels'), scene[s]) + scene_xr[s] = (("number_of_lines", "number_of_pixels"), scene[s]) input_data = xr.Dataset(viirs_data, coords=scene_xr) - input_data.drop_vars(['latitude', 'longitude']) + input_data.drop_vars(["latitude", "longitude"]) return input_data diff --git a/mvcm/restoral.py b/mvcm/restoral.py index 8da1cf4..a218f46 100644 --- a/mvcm/restoral.py +++ b/mvcm/restoral.py @@ -1,14 +1,16 @@ +import logging +from typing import Dict + import numpy as np import xarray as xr - -from typing import Dict -from attrs import define, field, validators, Factory +from attrs import Factory, define, field, validators from numpy.lib.stride_tricks import sliding_window_view + +import mvcm.utility_functions as utils + # from copy import deepcopy -import logging -import mvcm.utility_functions as utils # import importlib _bad_data = -999.0 @@ -75,11 +77,23 @@ logger = logging.getLogger(__name__) @define(kw_only=True, slots=True) class Restoral(object): - """ - """ - data: xr.Dataset = field(validator=[validators.instance_of(xr.Dataset), ]) - thresholds: str = field(validator=[validators.instance_of(Dict), ]) - scene_flag: Dict = field(validator=[validators.instance_of(Dict), ]) + """ """ + + data: xr.Dataset = field( + validator=[ + validators.instance_of(xr.Dataset), + ] + ) + thresholds: str = field( + validator=[ + validators.instance_of(Dict), + ] + ) + scene_flag: Dict = field( + validator=[ + validators.instance_of(Dict), + ] + ) # confidence: np.ndarray = field(validator=[validators.instance_of(np.ndarray), ]) # VIIRS does not have water vapor channels (6.7um and 7.2um) for this test @@ -87,39 +101,56 @@ class Restoral(object): pass # chk_spatial_var - def spatial_variability(self, - confidence: np.ndarray) -> np.ndarray: - - threshold = self.thresholds['Daytime_Ocean_Spatial_Variability'] + def spatial_variability( + self, confidence: np.ndarray, scn_idx: tuple, bits: Dict + ) -> np.ndarray: + threshold = self.thresholds["Daytime_Ocean_Spatial_Variability"] bt = self.data.M15.values reflectance = self.data.M07.values - var_bt = utils.spatial_var(bt, threshold['var_11um']) - var_reflectance = utils.spatial_var(reflectance, threshold['var_0_86um']) + var_bt = utils.spatial_var(bt, threshold["var_11um"]) + var_reflectance = utils.spatial_var(reflectance, threshold["var_0_86um"]) - idx = np.nonzero((confidence > 0.95) & (self.scene_flag['day'] == 1) & - (var_bt == 1) & (var_reflectance == 1)) + idx = np.nonzero( + (confidence > 0.95) + & (self.scene_flag["day"] == 1) + & (var_bt == 1) + & (var_reflectance == 1) + ) confidence[idx] = 1 - - idx = np.nonzero((confidence > 0.66) & (confidence <= 0.95) & ( - ((self.scene_flag['night'] == 1) & (var_bt == 1)) | - ((self.scene_flag['day'] == 1) & (var_bt == 1) & (var_reflectance == 1)))) + bits["test"][idx] = 1 + + idx = np.nonzero( + (confidence > 0.66) + & (confidence <= 0.95) + & ( + ((self.scene_flag["night"] == 1) & (var_bt == 1)) + | ( + (self.scene_flag["day"] == 1) + & (var_bt == 1) + & (var_reflectance == 1) + ) + ) + ) confidence[idx] = 0.96 + bits["test"][idx] = 1 idx = np.nonzero((confidence >= 0.05) & (confidence <= 0.66) & (var_bt == 1)) confidence[idx] = 0.67 + bits["test"][idx] = 1 - logger.info('spatial_variability restoral test performed successfully') + logger.info("spatial_variability restoral test performed successfully") - return confidence + return confidence[scn_idx], bits["test"][scn_idx] # chk_sunglint - def sunglint(self, - bit: np.ndarray, - confidence: np.ndarray) -> np.ndarray: - - threshold = self.thresholds['Sun_Glint'] - variability_threshold = self.thresholds['Daytime_Ocean_Spatial_Variability']['var_11um'] + def sunglint( + self, bit: np.ndarray, confidence: np.ndarray, scn_idx: tuple, bits: Dict + ) -> np.ndarray: + threshold = self.thresholds["Sun_Glint"] + variability_threshold = self.thresholds["Daytime_Ocean_Spatial_Variability"][ + "var_11um" + ] m02 = self.data.M07.values m09 = self.data.M02.values m20 = self.data.M12.values @@ -130,132 +161,225 @@ class Restoral(object): bt_var = utils.spatial_var(m31, variability_threshold) diff37_11 = m20 - m31 - idx = np.nonzero((bt_var == 1) & (bit == 1) & (m09 < threshold['m09']) & - (diff37_11 >= threshold['btdiff_37_11um'])) + idx = np.nonzero( + (bt_var == 1) + & (bit == 1) + & (m09 < threshold["m09"]) + & (diff37_11 >= threshold["btdiff_37_11um"]) + ) confidence[idx] = 0.67 - - idx = np.nonzero((bt_var == 1) & (bit == 1) & (m09 < threshold['m09']) & - (diff37_11 >= threshold['btdiff_37_11um']) & - (reg_var_mean*reg_std < threshold['m02'])) + bits["test"][idx] = 1 + + idx = np.nonzero( + (bt_var == 1) + & (bit == 1) + & (m09 < threshold["m09"]) + & (diff37_11 >= threshold["btdiff_37_11um"]) + & (reg_var_mean * reg_std < threshold["m02"]) + ) confidence[idx] = 0.96 + bits["test"][idx] = 1 - logger.info('sunglint restoral test performed successfully') + logger.info("sunglint restoral test performed successfully") - return confidence + return confidence[scn_idx], bits["test"][scn_idx] # chk_shallow_water - def shallow_water(self, - bits: Dict, - confidence: np.ndarray) -> np.ndarray: - - threshold = self.thresholds['Ocean_NDVI_and_Shallow_Water'] + def shallow_water( + self, bits: Dict, confidence: np.ndarray, scn_idx: tuple, bits_test: Dict + ) -> np.ndarray: + threshold = self.thresholds["Ocean_NDVI_and_Shallow_Water"] ndvi = utils.compute_ndvi(self.data.M05.values, self.data.M07.values) std_m05 = utils.local_nxn_standard_deviation(self.data.M08.values) - idx = np.nonzero((bits['ir'] == 1) & (bits['nir_1'] == 1) & (bits['nir_2'] == 1) & - (bits['sst'] == 0) & (self.scene_flag['sh_ocean'] == 0) & - ((ndvi <= threshold['ndvi_boundaries'][0]) | (ndvi >= threshold['ndvi_boundaries'][1])) & - (confidence < 1.0)) + idx = np.nonzero( + (bits["ir"] == 1) + & (bits["nir_1"] == 1) + & (bits["nir_2"] == 1) + & (bits["sst"] == 0) + & (self.scene_flag["sh_ocean"] == 0) + & ( + (ndvi <= threshold["ndvi_boundaries"][0]) + | (ndvi >= threshold["ndvi_boundaries"][1]) + ) + & (confidence < 1.0) + ) confidence[idx] = 1 - - idx = np.nonzero((bits['ir'] == 1) & (bits['nir_2'] == 1) & (self.scene_flag['sh_ocean'] == 1) & - (std_m05 < threshold['std_124']) & (confidence < 1.0)) + bits_test["test"][idx] = 1 + + idx = np.nonzero( + (bits["ir"] == 1) + & (bits["nir_2"] == 1) + & (self.scene_flag["sh_ocean"] == 1) + & (std_m05 < threshold["std_124"]) + & (confidence < 1.0) + ) confidence[idx] = 1 + bits_test["test"][idx] = 1 - logger.info('shallow_water restoral test performed successfully') + logger.info("shallow_water restoral test performed successfully") - return confidence + return confidence[scn_idx], bits_test["test"][scn_idx] # chk_land - def land(self, - bit: np.ndarray, - confidence: np.ndarray) -> np.ndarray: - + def land( + self, bit: np.ndarray, confidence: np.ndarray, scn_idx: tuple, bits: Dict + ) -> np.ndarray: rad = self.data.M15.values - m05m04ratio = self.data.M08.values/self.data.M04.values + m05m04ratio = self.data.M08.values / self.data.M04.values m20m22diff = self.data.M12.values - self.data.M13.values m22m31diff = self.data.M13.values - self.data.M15.values - threshold = self.thresholds['Land_Restoral'] + threshold = self.thresholds["Land_Restoral"] - thr11um = np.full((rad.ravel().shape[0], 3), np.array(threshold['bt11'])) - thr11um[self.data.ecosystem.values.ravel() == 8, :] = np.array(threshold['bt11_desert']) + thr11um = np.full((rad.ravel().shape[0], 3), np.array(threshold["bt11"])) + thr11um[self.data.ecosystem.values.ravel() == 8, :] = np.array( + threshold["bt11_desert"] + ) for i in range(3): - thr11um[:, i] = thr11um[:, i] - utils.bt11_elevation_correction(self.data.height.values).ravel() + thr11um[:, i] = ( + thr11um[:, i] + - utils.bt11_elevation_correction(self.data.height.values).ravel() + ) thr11um = thr11um.reshape((rad.shape[0], rad.shape[1], 3)) - thr_ratio = np.full((rad.ravel().shape[0], ), np.array(threshold['m05m04_ratio'])) - thr_ratio[self.data.ecosystem.values.ravel() == 8] = np.array(threshold['m05m04_ratio_desert']) + thr_ratio = np.full( + (rad.ravel().shape[0],), np.array(threshold["m05m04_ratio"]) + ) + thr_ratio[self.data.ecosystem.values.ravel() == 8] = np.array( + threshold["m05m04_ratio_desert"] + ) thr_ratio = thr_ratio.reshape(rad.shape) - idx = np.nonzero((bit == 1) & (self.scene_flag['day'] == 1) & (self.scene_flag['land'] == 1) & - (self.scene_flag['snow'] == 0) & (self.scene_flag['ice'] == 0) & - (confidence <= 0.95) & (rad > thr11um[:, :, 2])) + idx = np.nonzero( + (bit == 1) + & (self.scene_flag["day"] == 1) + & (self.scene_flag["land"] == 1) + & (self.scene_flag["snow"] == 0) + & (self.scene_flag["ice"] == 0) + & (confidence <= 0.95) + & (rad > thr11um[:, :, 2]) + ) confidence[idx] = 1 - idx = np.nonzero((bit == 1) & (self.scene_flag['day'] == 1) & (self.scene_flag['land'] == 1) & - (self.scene_flag['snow'] == 0) & (self.scene_flag['ice'] == 0) & - (confidence <= 0.95) & (rad > thr11um[:, :, 1])) + bits["test"][idx] = 1 + idx = np.nonzero( + (bit == 1) + & (self.scene_flag["day"] == 1) + & (self.scene_flag["land"] == 1) + & (self.scene_flag["snow"] == 0) + & (self.scene_flag["ice"] == 0) + & (confidence <= 0.95) + & (rad > thr11um[:, :, 1]) + & (rad <= thr11um[:, :, 2]) + ) confidence[idx] = 0.96 - idx = np.nonzero((bit == 1) & (self.scene_flag['day'] == 1) & (self.scene_flag['land'] == 1) & - (self.scene_flag['snow'] == 0) & (self.scene_flag['ice'] == 0) & - (confidence <= 0.95) & (rad > thr11um[:, :, 0])) + bits["test"][idx] = 1 + idx = np.nonzero( + (bit == 1) + & (self.scene_flag["day"] == 1) + & (self.scene_flag["land"] == 1) + & (self.scene_flag["snow"] == 0) + & (self.scene_flag["ice"] == 0) + & (confidence <= 0.95) + & (rad > thr11um[:, :, 0]) + & (rad <= thr11um[:, :, 1]) + ) confidence[idx] = 0.95 + # bits['test'][idx] = 1 # Perform additional clear-sky spectral tests - idx = np.nonzero((bit == 1) & (self.scene_flag['day'] == 1) & (self.scene_flag['land'] == 1) & - (self.scene_flag['snow'] == 0) & (self.scene_flag['ice'] == 0) & - (confidence <= 0.95) & (m20m22diff < threshold['m20m22diff']) & - (m22m31diff < threshold['m22m31diff']) & (m05m04ratio > thr_ratio)) + idx = np.nonzero( + (bit == 1) + & (self.scene_flag["day"] == 1) + & (self.scene_flag["land"] == 1) + & (self.scene_flag["snow"] == 0) + & (self.scene_flag["ice"] == 0) + & (confidence <= 0.95) + & (m20m22diff < threshold["m20m22diff"]) + & (m22m31diff < threshold["m22m31diff"]) + & (m05m04ratio > thr_ratio) + ) confidence[idx] = 0.96 + bits["test"][idx] = 1 - logger.info('land restoral test performed successfully') + logger.info("land restoral test performed successfully") - return confidence + return confidence[scn_idx], bits["test"][scn_idx] # chk_coast - def coast(self, - bit: np.ndarray, - confidence: np.ndarray) -> np.ndarray: - + def coast( + self, bit: np.ndarray, confidence: np.ndarray, scn_idx: tuple, bits: Dict + ) -> np.ndarray: ndvi = utils.compute_ndvi(self.data.M05.values, self.data.M07.values) - threshold = self.thresholds['Coastal_NDVI_Thresholds'] - - idx = np.nonzero((self.scene_flag['day'] == 1) & (self.scene_flag['land'] == 1) & - (self.scene_flag['snow'] == 0) & (self.scene_flag['ice'] == 0) & - (self.scene_flag['coast'] == 1) & - (bit == 1) & ((ndvi < threshold['coast_ndvi'][0]) | (ndvi >= threshold['coast_ndvi'][1]))) + threshold = self.thresholds["Coastal_NDVI_Thresholds"] + + idx = np.nonzero( + (self.scene_flag["day"] == 1) + & (self.scene_flag["land"] == 1) + & (self.scene_flag["snow"] == 0) + & (self.scene_flag["ice"] == 0) + & (self.scene_flag["coast"] == 1) + & (bit == 1) + & ( + (ndvi < threshold["coast_ndvi"][0]) + | (ndvi >= threshold["coast_ndvi"][1]) + ) + ) confidence[idx] = 1 + bits["test"][idx] = 1 - logger.info('coast restoral test performed successfully') + logger.info("coast restoral test performed successfully") - return confidence + return confidence[scn_idx], bits["test"][scn_idx] # chk_land_night - def land_night(self, - bit: np.ndarray, - confidence: np.ndarray) -> np.ndarray: - + def land_night( + self, bit: np.ndarray, confidence: np.ndarray, bits: Dict + ) -> np.ndarray: rad = self.data.M15.values - threshold = self.thresholds['Land_Restoral'] - thr11um = np.full((self.data.height.values.ravel().shape[0], 3), np.array(threshold['bt11'])) + threshold = self.thresholds["Land_Restoral"] + thr11um = np.full( + (self.data.height.values.ravel().shape[0], 3), np.array(threshold["bt11"]) + ) for i in range(3): - thr11um[:, i] = threshold['bt11_land_night'][i] - utils.bt11_elevation_correction(self.data.height.values).ravel() + thr11um[:, i] = ( + threshold["bt11_land_night"][i] + - utils.bt11_elevation_correction(self.data.height.values).ravel() + ) thr11um = thr11um.reshape((rad.shape[0], rad.shape[1], 3)) - idx = np.nonzero((self.scene_flag['night'] == 1) & (self.scene_flag['snow'] == 0) & - (self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) & - (rad > thr11um[:, :, 0])) + idx = np.nonzero( + (self.scene_flag["night"] == 1) + & (self.scene_flag["snow"] == 0) + & (self.scene_flag["ice"] == 0) + & (bit == 1) + & (confidence <= 0.95) + & (rad > thr11um[:, :, 0]) + ) confidence[idx] = 0.95 - idx = np.nonzero((self.scene_flag['night'] == 1) & (self.scene_flag['snow'] == 0) & - (self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) & - (rad > thr11um[:, :, 1])) + bits["test"][idx] = 1 + idx = np.nonzero( + (self.scene_flag["night"] == 1) + & (self.scene_flag["snow"] == 0) + & (self.scene_flag["ice"] == 0) + & (bit == 1) + & (confidence <= 0.95) + & (rad > thr11um[:, :, 1]) + ) confidence[idx] = 0.96 - idx = np.nonzero((self.scene_flag['night'] == 1) & (self.scene_flag['snow'] == 0) & - (self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) & - (rad > thr11um[:, :, 2])) + bits["test"][idx] = 1 + idx = np.nonzero( + (self.scene_flag["night"] == 1) + & (self.scene_flag["snow"] == 0) + & (self.scene_flag["ice"] == 0) + & (bit == 1) + & (confidence <= 0.95) + & (rad > thr11um[:, :, 2]) + ) confidence[idx] = 1 + bits["test"][idx] = 1 - logger.info('land night restoral test performed successfully') + logger.info("land night restoral test performed successfully") return confidence @@ -303,7 +427,7 @@ def compute_ndvi(b01: np.ndarray, ''' -''' +""" def sunglint(viirs_data, threshold, bit, conf): m09 = viirs_data.M02.values @@ -413,4 +537,4 @@ def coast(viirs_data, threshold, scene, conf): conf[idx] = 1 return conf -''' +""" diff --git a/mvcm/scene.py b/mvcm/scene.py index 26dcda3..2e7beef 100644 --- a/mvcm/scene.py +++ b/mvcm/scene.py @@ -1,38 +1,82 @@ import numpy as np + try: from ruamel import yaml as yml except ImportError: import ruamel_yaml as yml import logging -# from glob import glob -import mvcm.read_data as rd import ancillary_data as anc +import mvcm.read_data as rd + +# from glob import glob + # lsf: land sea flag -_scene_list = ['Ocean_Day', 'Ocean_Night', 'Land_Day', 'Land_Night', 'Day_Snow', 'Night_Snow', 'Coast_Day', - 'Land_Day_Desert', 'Antarctic_Day', 'Polar_Day_Snow', 'Polar_Day_Desert', 'Polar_Day_Ocean', - 'Polar_Day_Desert_Coast', 'Polar_Day_Coast', 'Polar_Day_Land', 'Polar_Night_Snow', - 'Polar_Night_Land', 'Polar_Night_Ocean', 'Land_Day_Desert_Coast', 'Land_Day_Coast', 'Desert', - 'Australia'] -_flags = ['day', 'night', 'land', 'coast', 'sh_lake', 'sh_ocean', 'water', 'polar', 'sunglint', - 'greenland', 'high_elevation', 'antarctica', 'desert', 'visusd', 'vrused', 'map_snow', 'map_ice', - 'ndsi_snow', 'snow', 'ice', 'new_zealand', 'uniform', 'australia'] +_scene_list = [ + "Ocean_Day", + "Ocean_Night", + "Land_Day", + "Land_Night", + "Day_Snow", + "Night_Snow", + "Coast_Day", + "Land_Day_Desert", + "Antarctic_Day", + "Polar_Day_Snow", + "Polar_Day_Desert", + "Polar_Day_Ocean", + "Polar_Day_Desert_Coast", + "Polar_Day_Coast", + "Polar_Day_Land", + "Polar_Night_Snow", + "Polar_Night_Land", + "Polar_Night_Ocean", + "Land_Day_Desert_Coast", + "Land_Day_Coast", + "Desert", + "Australia", +] +_flags = [ + "day", + "night", + "land", + "coast", + "sh_lake", + "sh_ocean", + "water", + "polar", + "sunglint", + "greenland", + "high_elevation", + "antarctica", + "desert", + "visusd", + "vrused", + "map_snow", + "map_ice", + "ndsi_snow", + "snow", + "ice", + "new_zealand", + "uniform", + "australia", +] # temp value, need to verify what the actual bad_data value is in the C code _bad_data = -999.0 -_dtr = np.pi/180. -_rtd = 180./np.pi +_dtr = np.pi / 180.0 +_rtd = 180.0 / np.pi logger = logging.getLogger(__name__) # I'm defining here the flags for difference scenes. Eventually I want to find a better way of doing this land = 1 # coast = .2 -sh_lake = .3 -sh_ocean = .4 +sh_lake = 0.3 +sh_ocean = 0.4 water = 5 polar = 60 sunglint = 70 @@ -42,41 +86,41 @@ night = 200 # #################################################################### # # TEST CASE # data: -#datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires' -#fname_mod02 = glob(f'{datapath}/VNP02MOD.A2022173.1448.001.*.uwssec*.nc')[0] -#fname_mod03 = glob(f'{datapath}/VNP03MOD.A2022173.1448.001.*.uwssec.nc')[0] -#fname_img02 = glob(f'{datapath}/VNP02IMG.A2022173.1448.001.*.uwssec*.nc')[0] -#fname_img03 = glob(f'{datapath}/VNP03IMG.A2022173.1448.001.*.uwssec.nc')[0] +# datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires' +# fname_mod02 = glob(f'{datapath}/VNP02MOD.A2022173.1448.001.*.uwssec*.nc')[0] +# fname_mod03 = glob(f'{datapath}/VNP03MOD.A2022173.1448.001.*.uwssec.nc')[0] +# fname_img02 = glob(f'{datapath}/VNP02IMG.A2022173.1448.001.*.uwssec*.nc')[0] +# fname_img03 = glob(f'{datapath}/VNP03IMG.A2022173.1448.001.*.uwssec.nc')[0] # thresholds: -#threshold_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml' +# threshold_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml' # ancillary files: -geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4' -geos_atm_2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4' -geos_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1430.V01.nc4' -geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1430.V01.nc4' -geos_constants = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4' -ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf' -sst_file = 'oisst.20220622' -eco_file = 'goge1_2_img.v1' +geos_atm_1 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4" +geos_atm_2 = "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4" +geos_land = "GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1430.V01.nc4" +geos_ocean = "GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1430.V01.nc4" +geos_constants = "GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4" +ndvi_file = "NDVI.FM.c004.v2.0.WS.00-04.177.hdf" +sst_file = "oisst.20220622" +eco_file = "goge1_2_img.v1" # #################################################################### # def test_scene(): - - ancillary_file_names = {'GEOS_atm_1': geos_atm_1, - 'GEOS_atm_2': geos_atm_2, - 'GEOS_land': geos_land, - 'GEOS_ocean': geos_ocean, - 'GEOS_constants': geos_constants, - 'NDVI': ndvi_file, - 'SST': sst_file, - 'ANC_DIR': f'{datapath}/ancillary' - } - - viirs_data = rd.read_data('viirs', f'{fname_mod02}', f'{fname_mod03}') + ancillary_file_names = { + "GEOS_atm_1": geos_atm_1, + "GEOS_atm_2": geos_atm_2, + "GEOS_land": geos_land, + "GEOS_ocean": geos_ocean, + "GEOS_constants": geos_constants, + "NDVI": ndvi_file, + "SST": sst_file, + "ANC_DIR": f"{datapath}/ancillary", + } + + viirs_data = rd.read_data("viirs", f"{fname_mod02}", f"{fname_mod03}") viirs_data = rd.read_ancillary_data(ancillary_file_names, viirs_data) @@ -84,7 +128,7 @@ def test_scene(): text = f.read() thresholds = yml.safe_load(text) - sunglint_angle = thresholds['Sun_Glint']['bounds'][3] + sunglint_angle = thresholds["Sun_Glint"]["bounds"][3] scene = find_scene(viirs_data, sunglint_angle) @@ -92,29 +136,31 @@ def test_scene(): def find_scene(data, sunglint_angle): - eco = np.array(data['ecosystem'].values, dtype=np.uint8) + eco = np.array(data["ecosystem"].values, dtype=np.uint8) # logger.warning('"eco" has been renamed "ecosystem" in the most recent version of the read_data function') # eco = np.array(data['eco'].values, dtype=np.uint8) - snowfr = data['geos_snow_fraction'].values - icefr = data['geos_ice_fraction'].values - lsf = np.array(data['land_water_mask'].values, dtype=np.uint8) - lat = data['latitude'].values - lon = data['longitude'].values - sza = data['solar_zenith'].values - vza = data['sensor_zenith'].values - raz = data['relative_azimuth'].values - b065 = data['M05'].values - b086 = data['M07'].values - elev = data['height'].values - ndvibk = data['ndvi'].values - sfctmp = data['geos_surface_temperature'].values + snowfr = data["geos_snow_fraction"].values + icefr = data["geos_ice_fraction"].values + lsf = np.array(data["land_water_mask"].values, dtype=np.uint8) + lat = data["latitude"].values + lon = data["longitude"].values + sza = data["solar_zenith"].values + vza = data["sensor_zenith"].values + raz = data["relative_azimuth"].values + b065 = data["M05"].values + b086 = data["M07"].values + elev = data["height"].values + ndvibk = data["ndvi"].values + sfctmp = data["geos_surface_temperature"].values dim1 = data.latitude.shape[0] dim2 = data.latitude.shape[1] day = np.zeros((dim1, dim2)) day[sza <= 85] = 1 - cos_refang = np.sin(vza*_dtr) * np.sin(sza*_dtr) * np.cos(raz*_dtr) + np.cos(vza*_dtr) * np.cos(sza*_dtr) + cos_refang = np.sin(vza * _dtr) * np.sin(sza * _dtr) * np.cos(raz * _dtr) + np.cos( + vza * _dtr + ) * np.cos(sza * _dtr) refang = np.arccos(cos_refang) * _rtd # tmp = np.ones((dim1, dim2)) @@ -123,97 +169,129 @@ def find_scene(data, sunglint_angle): scene_flag = {flg: np.zeros((dim1, dim2)) for flg in _flags} - scene_flag['day'][sza <= 85] = 1 - scene_flag['visusd'][sza <= 85] = 1 - scene_flag['night'][sza > 85] = 1 + scene_flag["day"][sza <= 85] = 1 + scene_flag["visusd"][sza <= 85] = 1 + scene_flag["night"][sza > 85] = 1 - scene_flag['polar'][np.abs(lat) > 60] = 1 + scene_flag["polar"][np.abs(lat) > 60] = 1 # ################# need to pass refang (once I figure out what it is) and sunglint_angle. The latter # comes from the thresholds file. In the C code is snglnt_bounds[3] - idx = np.nonzero((scene_flag['day'] == 1) & (refang <= sunglint_angle)) - scene_flag['sunglint'][idx] = 1 + idx = np.nonzero((scene_flag["day"] == 1) & (refang <= sunglint_angle)) + scene_flag["sunglint"][idx] = 1 # Force consistency between lsf and ecosystem type for water idx = np.nonzero((lsf == 0) | ((lsf >= 5) & (lsf < 7))) eco[idx] = 14 # start by defining everything as land - scene_flag['land'] = np.ones((dim1, dim2)) - scene_flag['water'] = np.zeros((dim1, dim2)) + scene_flag["land"] = np.ones((dim1, dim2)) + scene_flag["water"] = np.zeros((dim1, dim2)) # Fix-up for missing ecosystem data in eastern Greenland and north-eastern Siberia. # Without this, these regions become completely "coast". idx = np.nonzero((lsf != 255) & ((lsf == 1) | (lsf == 4))) - scene_flag['land'][idx] = 1 + scene_flag["land"][idx] = 1 # idx = np.nonzero((lsf != 255) & (eco == 14)) idx = np.nonzero(((lsf == 1) | (lsf == 4)) & (eco == 14) & (lat < 64.0)) - scene_flag['coast'][idx] = 1 - - idx = np.nonzero((eco == 14) & ((lsf == 1) | (lsf == 4)) & - ((lat >= 67.5) & (lon < -40.0) & (lon > -168.6))) - scene_flag['coast'][idx] = 1 - idx = np.nonzero((eco == 14) & ((lsf == 1) | (lsf == 4)) & - ((lat >= 67.5) & (lon > -12.5))) - scene_flag['coast'][idx] = 1 - - idx = np.nonzero((eco == 14) & ((lsf == 1) | (lsf == 4)) & - ((lat >= 64.0) & (lat < 67.5) & (lon < -40.0) & (lon > -168.5))) - scene_flag['coast'][idx] = 1 - idx = np.nonzero((eco == 14) & ((lsf == 1) | (lsf == 4)) & - ((lat >= 64.0) & (lat < 67.5) & (lon > -30.0))) - scene_flag['coast'][idx] = 1 + scene_flag["coast"][idx] = 1 + + idx = np.nonzero( + (eco == 14) + & ((lsf == 1) | (lsf == 4)) + & ((lat >= 67.5) & (lon < -40.0) & (lon > -168.6)) + ) + scene_flag["coast"][idx] = 1 + idx = np.nonzero( + (eco == 14) & ((lsf == 1) | (lsf == 4)) & ((lat >= 67.5) & (lon > -12.5)) + ) + scene_flag["coast"][idx] = 1 + + idx = np.nonzero( + (eco == 14) + & ((lsf == 1) | (lsf == 4)) + & ((lat >= 64.0) & (lat < 67.5) & (lon < -40.0) & (lon > -168.5)) + ) + scene_flag["coast"][idx] = 1 + idx = np.nonzero( + (eco == 14) + & ((lsf == 1) | (lsf == 4)) + & ((lat >= 64.0) & (lat < 67.5) & (lon > -30.0)) + ) + scene_flag["coast"][idx] = 1 idx = np.nonzero(lsf == 2) - scene_flag['coast'][idx] = 1 - scene_flag['land'][idx] = 1 + scene_flag["coast"][idx] = 1 + scene_flag["land"][idx] = 1 idx = np.nonzero(lsf == 3) - scene_flag['land'][idx] = 1 - scene_flag['sh_lake'][idx] = 1 + scene_flag["land"][idx] = 1 + scene_flag["sh_lake"][idx] = 1 idx = np.nonzero((lsf == 0) | ((lsf >= 5) & (lsf <= 7))) - scene_flag['water'][idx] = 1 - scene_flag['land'][idx] = 0 + scene_flag["water"][idx] = 1 + scene_flag["land"][idx] = 0 - scene_flag['sh_ocean'][lsf == 0] = 1 + scene_flag["sh_ocean"][lsf == 0] = 1 # Need shallow lakes to be processed as "coast" for day, but not night idx = np.nonzero((lsf == 3) & (day == 1)) - scene_flag['coast'][idx] = 1 + scene_flag["coast"][idx] = 1 # if land/sea flag is missing, then calculate visible ratio to determine if land or water. - idx = np.nonzero((lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086/b065 > 0.9)) - scene_flag['land'][idx] = 1 + idx = np.nonzero( + (lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086 / b065 > 0.9) + ) + scene_flag["land"][idx] = 1 - idx = np.nonzero((lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086/b065 <= 0.9)) - scene_flag['land'][idx] = 0 - scene_flag['water'][idx] = 1 + idx = np.nonzero( + (lsf == 255) & (b065 != _bad_data) & (b086 != _bad_data) & (b086 / b065 <= 0.9) + ) + scene_flag["land"][idx] = 0 + scene_flag["water"][idx] = 1 # Check surface elevation # First, define "Greenland". - idx = np.nonzero((scene_flag['land'] == 1) & - (lat >= 60.0) & (lat < 67.0) & (lon >= -60.0) & (lon < -30.0)) - scene_flag['greenland'][idx] = 1 - - idx = np.nonzero((scene_flag['land'] == 1) & - (lat >= 67.0) & (lat < 75.0) & (lon >= -60.0) & (lon < -10.0)) - scene_flag['greenland'][idx] = 1 - - idx = np.nonzero((scene_flag['land'] == 1) & - (lat >= 75.0) & (lon >= -70.0) & (lon < -10.0)) - scene_flag['greenland'][idx] = 1 - - scene_flag['high_elevation'][elev > 2000] = 1 - idx = np.nonzero((elev > 200) & (scene_flag['greenland'] == 1) & (scene_flag['land'] == 1)) - scene_flag['high_elevation'][idx] = 1 - - idx = np.nonzero((lat >= 75.0) & (lat <= 79.0) & (lon >= -73.0) & (lon <= -50.0) & - (scene_flag['land'] == 1)) - scene_flag['high_elevation'][idx] = 1 - - scene_flag['antarctica'][lat < -60.0] = 1 + idx = np.nonzero( + (scene_flag["land"] == 1) + & (lat >= 60.0) + & (lat < 67.0) + & (lon >= -60.0) + & (lon < -30.0) + ) + scene_flag["greenland"][idx] = 1 + + idx = np.nonzero( + (scene_flag["land"] == 1) + & (lat >= 67.0) + & (lat < 75.0) + & (lon >= -60.0) + & (lon < -10.0) + ) + scene_flag["greenland"][idx] = 1 + + idx = np.nonzero( + (scene_flag["land"] == 1) & (lat >= 75.0) & (lon >= -70.0) & (lon < -10.0) + ) + scene_flag["greenland"][idx] = 1 + + scene_flag["high_elevation"][elev > 2000] = 1 + idx = np.nonzero( + (elev > 200) & (scene_flag["greenland"] == 1) & (scene_flag["land"] == 1) + ) + scene_flag["high_elevation"][idx] = 1 + + idx = np.nonzero( + (lat >= 75.0) + & (lat <= 79.0) + & (lon >= -73.0) + & (lon <= -50.0) + & (scene_flag["land"] == 1) + ) + scene_flag["high_elevation"][idx] = 1 + + scene_flag["antarctica"][lat < -60.0] = 1 ################################## # somewhere here I need to add # @@ -221,37 +299,48 @@ def find_scene(data, sunglint_angle): ################################## # this is a temporary variable for the 11um elevation correction - elev_correction = elev/1000.0 * 5.0 + elev_correction = elev / 1000.0 * 5.0 ## Get surface temperature from NWP and SST fields ## if it's land get it from GDAS/GEOS5 - #sfctmp[scene_flag['land'] == 1] = sfct + # sfctmp[scene_flag['land'] == 1] = sfct ## otherwise use the ReynoldsSST - #sfctmp[scene_flag['land'] == 0] = reynSST + # sfctmp[scene_flag['land'] == 0] = reynSST # Use background NDVI to define "desert" - idx = np.nonzero((scene_flag['land'] == 1) & (ndvibk < 0.3)) - scene_flag['desert'][idx] = 1 - idx = np.nonzero((scene_flag['land'] == 1) & (lat < -69.0)) - scene_flag['desert'][idx] = 1 - - idx = np.nonzero((eco == 2) | (eco == 8) | (eco == 11) | (eco == 40) | (eco == 41) | (eco == 46) | - (eco == 51) | (eco == 52) | (eco == 59) | (eco == 71) | (eco == 50)) - scene_flag['vrused'] = np.ones((dim1, dim2)) - scene_flag['vrused'][idx] = 0 - - snow_fraction = data['geos_snow_fraction'] - perm_ice_fraction = data['geos_land_ice_fraction'] - ice_fraction = data['geos_ice_fraction'] + idx = np.nonzero((scene_flag["land"] == 1) & (ndvibk < 0.3)) + scene_flag["desert"][idx] = 1 + idx = np.nonzero((scene_flag["land"] == 1) & (lat < -69.0)) + scene_flag["desert"][idx] = 1 + + idx = np.nonzero( + (eco == 2) + | (eco == 8) + | (eco == 11) + | (eco == 40) + | (eco == 41) + | (eco == 46) + | (eco == 51) + | (eco == 52) + | (eco == 59) + | (eco == 71) + | (eco == 50) + ) + scene_flag["vrused"] = np.ones((dim1, dim2)) + scene_flag["vrused"][idx] = 0 + + snow_fraction = data["geos_snow_fraction"] + perm_ice_fraction = data["geos_land_ice_fraction"] + ice_fraction = data["geos_ice_fraction"] idx = tuple(np.nonzero((snow_fraction > 0.10) & (snow_fraction <= 1.0))) - scene_flag['map_snow'][idx] = 1 + scene_flag["map_snow"][idx] = 1 idx = tuple(np.nonzero((perm_ice_fraction > 0.10) & (perm_ice_fraction <= 1.0))) - scene_flag['map_snow'][idx] = 1 + scene_flag["map_snow"][idx] = 1 idx = tuple(np.nonzero((ice_fraction > 0.10) & (ice_fraction <= 1.0))) - scene_flag['map_ice'][idx] = 1 + scene_flag["map_ice"][idx] = 1 # need to define this function and write this block better # if day == 1: @@ -260,198 +349,310 @@ def find_scene(data, sunglint_angle): # # scene_flag['ndsi_snow'] = run_snow_mask() # scene_flag['ndsi_snow'] = np.zeros((dim1, dim2)) - idx = np.nonzero((day == 1) & (water == 1) & (lat >= -60.0) & (lat <= 25.0) & - (scene_flag['map_snow'] == 1) & (scene_flag['ndsi_snow'] == 1)) - scene_flag['ice'][idx] = 1 - - idx = np.nonzero((day == 1) & (water == 1) & (lat < -60.0) & - (scene_flag['ndsi_snow'] == 1)) - scene_flag['ice'][idx] = 1 - - idx = np.nonzero((day == 1) & (water == 1) & (lsf == 3) | (lsf == 5) & - (scene_flag['ndsi_snow'] == 1)) - scene_flag['ice'][idx] = 1 - - idx = np.nonzero((day == 1) & (water == 1) & - (scene_flag['map_ice'] == 1) & (scene_flag['ndsi_snow'] == 1)) - scene_flag['ice'][idx] = 1 + idx = np.nonzero( + (day == 1) + & (water == 1) + & (lat >= -60.0) + & (lat <= 25.0) + & (scene_flag["map_snow"] == 1) + & (scene_flag["ndsi_snow"] == 1) + ) + scene_flag["ice"][idx] = 1 + + idx = np.nonzero( + (day == 1) & (water == 1) & (lat < -60.0) & (scene_flag["ndsi_snow"] == 1) + ) + scene_flag["ice"][idx] = 1 + + idx = np.nonzero( + (day == 1) & (water == 1) & (lsf == 3) + | (lsf == 5) & (scene_flag["ndsi_snow"] == 1) + ) + scene_flag["ice"][idx] = 1 + + idx = np.nonzero( + (day == 1) + & (water == 1) + & (scene_flag["map_ice"] == 1) + & (scene_flag["ndsi_snow"] == 1) + ) + scene_flag["ice"][idx] = 1 # Define New Zealand region which receives snow but snow map does not show it. - idx = np.nonzero((day == 1) & (land == 1) & - (lat >= 48.0) & (lat <= -34.0) & (lon >= 165.0) & (lon <= 180.0)) - scene_flag['new_zealand'][idx] = 1 - - idx = np.nonzero((day == 1) & (land == 1) & (lat >= -60.0) & (lat <= 25.0) & - (scene_flag['map_snow'] == 1) & (scene_flag['ndsi_snow'] == 1) | - (scene_flag['new_zealand'] == 1)) - scene_flag['snow'][idx] = 1 + idx = np.nonzero( + (day == 1) + & (land == 1) + & (lat >= 48.0) + & (lat <= -34.0) + & (lon >= 165.0) + & (lon <= 180.0) + ) + scene_flag["new_zealand"][idx] = 1 + + idx = np.nonzero( + (day == 1) + & (land == 1) + & (lat >= -60.0) + & (lat <= 25.0) + & (scene_flag["map_snow"] == 1) + & (scene_flag["ndsi_snow"] == 1) + | (scene_flag["new_zealand"] == 1) + ) + scene_flag["snow"][idx] = 1 idx = np.nonzero((day == 1) & (land == 1) & (lat < -60.0)) - scene_flag['snow'][idx] = 1 + scene_flag["snow"][idx] = 1 - idx = np.nonzero((day == 1) & (land == 1) & (scene_flag['ndsi_snow'] == 1)) - scene_flag['snow'][idx] = 1 + idx = np.nonzero((day == 1) & (land == 1) & (scene_flag["ndsi_snow"] == 1)) + scene_flag["snow"][idx] = 1 - idx = np.nonzero((day == 0) & (scene_flag['map_snow'] == 1) & - (sfctmp > 280.0) & (elev < 500.0)) - scene_flag['snow'][idx] = 0 + idx = np.nonzero( + (day == 0) & (scene_flag["map_snow"] == 1) & (sfctmp > 280.0) & (elev < 500.0) + ) + scene_flag["snow"][idx] = 0 - idx = np.nonzero((day == 0) & (scene_flag['map_snow'] == 1) & - (sfctmp > 280.0) & (elev < 500.0)) - scene_flag['ice'][idx] = 0 + idx = np.nonzero( + (day == 0) & (scene_flag["map_snow"] == 1) & (sfctmp > 280.0) & (elev < 500.0) + ) + scene_flag["ice"][idx] = 0 idx = np.nonzero((day == 0) & (lat > 86.0)) - scene_flag['ice'][idx] = 1 + scene_flag["ice"][idx] = 1 idx = np.nonzero((lat <= -11) & (lat > -40) & (lon <= 155) & (lon > 110)) - scene_flag['australia'][idx] = 1 + scene_flag["australia"][idx] = 1 # Check regional uniformity # Check for data border pixels # eco_out = np.array(np.empty(eco.shape), dtype=np.uint8) - scene_flag['uniform'] = anc.py_check_reg_uniformity(eco, eco, snowfr, icefr, lsf)['loc_uniform'] + scene_flag["uniform"] = anc.py_check_reg_uniformity(eco, eco, snowfr, icefr, lsf)[ + "loc_uniform" + ] # print(scene_flag['uniform'][:10, :10]) # NOTE: These three lines need to be verified. They seem correct but I need to make sure - scene_flag['water'][scene_flag['uniform'] < 0] = 0 - scene_flag['coast'][scene_flag['uniform'] < 0] = 1 - scene_flag['land'][scene_flag['uniform'] < 0] = 1 + scene_flag["water"][scene_flag["uniform"] < 0] = 0 + scene_flag["coast"][scene_flag["uniform"] < 0] = 1 + scene_flag["land"][scene_flag["uniform"] < 0] = 1 # TEMP VALUES FOR DEBUGGING - scene_flag['lat'] = lat - scene_flag['lon'] = lon + scene_flag["lat"] = lat + scene_flag["lon"] = lon return scene_flag def scene_id(scene_flag): - - dim1, dim2 = scene_flag['day'].shape[0], scene_flag['day'].shape[1] + dim1, dim2 = scene_flag["day"].shape[0], scene_flag["day"].shape[1] scene = {scn: np.zeros((dim1, dim2)) for scn in _scene_list} # Ocean Day - idx = np.nonzero((scene_flag['water'] == 1) & (scene_flag['day'] == 1) & - ((scene_flag['ice'] == 0) | (scene_flag['snow'] == 0))) # & - # (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0) & - # (scene_flag['coast'] == 0) & (scene_flag['desert'] == 0)) - scene['Ocean_Day'][idx] = 1 + idx = np.nonzero( + (scene_flag["water"] == 1) + & (scene_flag["day"] == 1) + & ((scene_flag["ice"] == 0) | (scene_flag["snow"] == 0)) + ) # & + # (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0) & + # (scene_flag['coast'] == 0) & (scene_flag['desert'] == 0)) + scene["Ocean_Day"][idx] = 1 # Ocean Night - idx = np.nonzero((scene_flag['water'] == 1) & (scene_flag['night'] == 1) & - ((scene_flag['ice'] == 0) | (scene_flag['snow'] == 0))) # & - # (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0) & - # (scene_flag['coast'] == 0) & (scene_flag['desert'] == 0)) - scene['Ocean_Night'][idx] = 1 + idx = np.nonzero( + (scene_flag["water"] == 1) + & (scene_flag["night"] == 1) + & ((scene_flag["ice"] == 0) | (scene_flag["snow"] == 0)) + ) # & + # (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0) & + # (scene_flag['coast'] == 0) & (scene_flag['desert'] == 0)) + scene["Ocean_Night"][idx] = 1 # Land Day - idx = np.nonzero((scene_flag['land'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['ice'] == 0) & (scene_flag['snow'] == 0) & - (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0) & - (scene_flag['coast'] == 0) & (scene_flag['desert'] == 0)) - scene['Land_Day'][idx] = 1 + idx = np.nonzero( + (scene_flag["land"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["ice"] == 0) + & (scene_flag["snow"] == 0) + & (scene_flag["polar"] == 0) + & (scene_flag["antarctica"] == 0) + & (scene_flag["coast"] == 0) + & (scene_flag["desert"] == 0) + ) + scene["Land_Day"][idx] = 1 # Land Night - idx = np.nonzero((scene_flag['land'] == 1) & (scene_flag['night'] == 1) & - (scene_flag['ice'] == 0) & (scene_flag['snow'] == 0) & - (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0) & - (scene_flag['coast'] == 0)) - scene['Land_Night'][idx] = 1 + idx = np.nonzero( + (scene_flag["land"] == 1) + & (scene_flag["night"] == 1) + & (scene_flag["ice"] == 0) + & (scene_flag["snow"] == 0) + & (scene_flag["polar"] == 0) + & (scene_flag["antarctica"] == 0) + & (scene_flag["coast"] == 0) + ) + scene["Land_Night"][idx] = 1 # Snow Day - idx = np.nonzero((scene_flag['day'] == 1) & - ((scene_flag['ice'] == 1) | (scene_flag['snow'] == 1)) & - (scene_flag['polar'] == 1) & (scene_flag['antarctica'] == 1)) - scene['Day_Snow'][idx] = 1 + idx = np.nonzero( + (scene_flag["day"] == 1) + & ((scene_flag["ice"] == 1) | (scene_flag["snow"] == 1)) + & (scene_flag["polar"] == 1) + & (scene_flag["antarctica"] == 1) + ) + scene["Day_Snow"][idx] = 1 # Snow Night - idx = np.nonzero((scene_flag['night'] == 1) & - ((scene_flag['ice'] == 1) | (scene_flag['snow'] == 1)) & - (scene_flag['polar'] == 1) & (scene_flag['antarctica'] == 1)) - scene['Night_Snow'][idx] = 1 + idx = np.nonzero( + (scene_flag["night"] == 1) + & ((scene_flag["ice"] == 1) | (scene_flag["snow"] == 1)) + & (scene_flag["polar"] == 1) + & (scene_flag["antarctica"] == 1) + ) + scene["Night_Snow"][idx] = 1 # Land Day Coast - idx = np.nonzero((scene_flag['land'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['coast'] == 1) & (scene_flag['desert'] == 0) & - (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0)) - scene['Land_Day_Coast'][idx] = 1 + idx = np.nonzero( + (scene_flag["land"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["coast"] == 1) + & (scene_flag["desert"] == 0) + & (scene_flag["polar"] == 0) + & (scene_flag["antarctica"] == 0) + ) + scene["Land_Day_Coast"][idx] = 1 # Land Day Desert - idx = np.nonzero((scene_flag['land'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['desert'] == 1) & (scene_flag['coast'] == 0) & - (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0)) - scene['Land_Day_Desert'][idx] = 1 + idx = np.nonzero( + (scene_flag["land"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["desert"] == 1) + & (scene_flag["coast"] == 0) + & (scene_flag["polar"] == 0) + & (scene_flag["antarctica"] == 0) + ) + scene["Land_Day_Desert"][idx] = 1 # Land Day Desert Coast - idx = np.nonzero((scene_flag['land'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['desert'] == 1) & (scene_flag['coast'] == 1) & - (scene_flag['polar'] == 0) & (scene_flag['antarctica'] == 0)) - scene['Land_Day_Desert_Coast'][idx] = 1 + idx = np.nonzero( + (scene_flag["land"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["desert"] == 1) + & (scene_flag["coast"] == 1) + & (scene_flag["polar"] == 0) + & (scene_flag["antarctica"] == 0) + ) + scene["Land_Day_Desert_Coast"][idx] = 1 # Antarctic Day - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['antarctica'] == 1) & (scene_flag['land'] == 1)) - scene['Antarctic_Day'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["antarctica"] == 1) + & (scene_flag["land"] == 1) + ) + scene["Antarctic_Day"][idx] = 1 # Polar Day Snow - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['day'] == 1) & - ((scene_flag['snow'] == 1) | (scene_flag['ice'] == 1)) & - (scene_flag['antarctica'] == 0)) - scene['Polar_Day_Snow'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["day"] == 1) + & ((scene_flag["snow"] == 1) | (scene_flag["ice"] == 1)) + & (scene_flag["antarctica"] == 0) + ) + scene["Polar_Day_Snow"][idx] = 1 # Polar Day Desert - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['land'] == 1) & (scene_flag['desert'] == 1) & - (scene_flag['coast'] == 0) & (scene_flag['antarctica'] == 0) & - (scene_flag['ice'] == 0) & (scene_flag['snow'] == 0)) - scene['Polar_Day_Desert'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["land"] == 1) + & (scene_flag["desert"] == 1) + & (scene_flag["coast"] == 0) + & (scene_flag["antarctica"] == 0) + & (scene_flag["ice"] == 0) + & (scene_flag["snow"] == 0) + ) + scene["Polar_Day_Desert"][idx] = 1 # Polar Day Desert Coast - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['land'] == 1) & (scene_flag['desert'] == 1) & - (scene_flag['coast'] == 1) & (scene_flag['antarctica'] == 0) & - (scene_flag['ice'] == 0) & (scene_flag['snow'] == 0)) - scene['Polar_Day_Desert_Coast'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["land"] == 1) + & (scene_flag["desert"] == 1) + & (scene_flag["coast"] == 1) + & (scene_flag["antarctica"] == 0) + & (scene_flag["ice"] == 0) + & (scene_flag["snow"] == 0) + ) + scene["Polar_Day_Desert_Coast"][idx] = 1 # Polar Day Coast - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['land'] == 1) & (scene_flag['coast'] == 1) & - (scene_flag['desert'] == 0) & (scene_flag['antarctica'] == 0) & - (scene_flag['ice'] == 0) & (scene_flag['snow'] == 0)) - scene['Polar_Day_Coast'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["land"] == 1) + & (scene_flag["coast"] == 1) + & (scene_flag["desert"] == 0) + & (scene_flag["antarctica"] == 0) + & (scene_flag["ice"] == 0) + & (scene_flag["snow"] == 0) + ) + scene["Polar_Day_Coast"][idx] = 1 # Polar Day Land - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['land'] == 1) & (scene_flag['coast'] == 0) & - (scene_flag['desert'] == 0) & (scene_flag['antarctica'] == 0) & - (scene_flag['ice'] == 0) & (scene_flag['snow'] == 0)) - scene['Polar_Day_Land'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["land"] == 1) + & (scene_flag["coast"] == 0) + & (scene_flag["desert"] == 0) + & (scene_flag["antarctica"] == 0) + & (scene_flag["ice"] == 0) + & (scene_flag["snow"] == 0) + ) + scene["Polar_Day_Land"][idx] = 1 # Polar Day Ocean - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['day'] == 1) & - (scene_flag['water'] == 1) & (scene_flag['coast'] == 0) & - (scene_flag['desert'] == 0) & (scene_flag['antarctica'] == 0) & - (scene_flag['ice'] == 0) & (scene_flag['snow'] == 0)) - scene['Polar_Day_Ocean'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["day"] == 1) + & (scene_flag["water"] == 1) + & (scene_flag["coast"] == 0) + & (scene_flag["desert"] == 0) + & (scene_flag["antarctica"] == 0) + & (scene_flag["ice"] == 0) + & (scene_flag["snow"] == 0) + ) + scene["Polar_Day_Ocean"][idx] = 1 # Polar Night Snow - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['night'] == 1) & - ((scene_flag['snow'] == 1) | (scene_flag['ice'] == 1))) - scene['Polar_Night_Snow'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["night"] == 1) + & ((scene_flag["snow"] == 1) | (scene_flag["ice"] == 1)) + ) + scene["Polar_Night_Snow"][idx] = 1 # Polar Night Land - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['night'] == 1) & - (scene_flag['land'] == 1)) - scene['Polar_Night_Land'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["night"] == 1) + & (scene_flag["land"] == 1) + ) + scene["Polar_Night_Land"][idx] = 1 # Polar Night Ocean - idx = np.nonzero((scene_flag['polar'] == 1) & (scene_flag['night'] == 1) & - (scene_flag['water'] == 1)) - scene['Polar_Night_Ocean'][idx] = 1 - - idx = np.nonzero(scene_flag['desert'] == 1) - scene['Desert'][idx] = 1 - - idx = np.nonzero(scene_flag['australia'] == 1) - scene['Australia'][idx] = 1 + idx = np.nonzero( + (scene_flag["polar"] == 1) + & (scene_flag["night"] == 1) + & (scene_flag["water"] == 1) + ) + scene["Polar_Night_Ocean"][idx] = 1 + + idx = np.nonzero(scene_flag["desert"] == 1) + scene["Desert"][idx] = 1 + + idx = np.nonzero(scene_flag["australia"] == 1) + scene["Australia"][idx] = 1 return scene diff --git a/mvcm/spectral_tests.py b/mvcm/spectral_tests.py index f4816a9..a4d56aa 100644 --- a/mvcm/spectral_tests.py +++ b/mvcm/spectral_tests.py @@ -1,29 +1,46 @@ +import functools +import importlib +import logging +from typing import Dict + import numpy as np import xarray as xr - +from attrs import Factory, define, field, validators from numpy.lib.stride_tricks import sliding_window_view -from attrs import define, field, validators, Factory -from typing import Dict - -import functools -import logging # import utils import mvcm.conf as conf -# import dev.conf_xr as conf_xr -import mvcm.scene as scn import mvcm.preprocess_thresholds as preproc import mvcm.restoral as restoral +import mvcm.scene as scn + # import mvcm.utility_functions as utils -import importlib -_RTD = 180./np.pi -_DTR = np.pi/180 -_scene_list = ['Ocean_Day', 'Ocean_Night', 'Land_Day', 'Land_Night', 'Day_Snow', 'Night_Snow', 'Coast_Day', - 'Land_Day_Desert', 'Antarctic_Day', 'Polar_Day_Snow', 'Polar_Day_Desert', 'Polar_Day_Ocean', - 'Polar_Day_Desert_Coast', 'Polar_Day_Coast', 'Polar_Day_Land', 'Polar_Night_Snow', - 'Polar_Night_Land', 'Polar_Night_Ocean', 'Land_Day_Desert_Coast', 'Land_Day_Coast', ] +_RTD = 180.0 / np.pi +_DTR = np.pi / 180 +_scene_list = [ + "Ocean_Day", + "Ocean_Night", + "Land_Day", + "Land_Night", + "Day_Snow", + "Night_Snow", + "Coast_Day", + "Land_Day_Desert", + "Antarctic_Day", + "Polar_Day_Snow", + "Polar_Day_Desert", + "Polar_Day_Ocean", + "Polar_Day_Desert_Coast", + "Polar_Day_Coast", + "Polar_Day_Land", + "Polar_Night_Snow", + "Polar_Night_Land", + "Polar_Night_Ocean", + "Land_Day_Desert_Coast", + "Land_Day_Coast", +] # # 'Desert', 'Australia'] # this is used for testing, eventually we want to remove it @@ -35,8 +52,8 @@ logger = logging.getLogger(__name__) def restoral_flag(bits): - qa_bit = bits['qa_bit'] - test_bit = bits['test_bit'] + qa_bit = bits["qa_bit"] + test_bit = bits["test_bit"] # the logic of this is convoluted as hell, but basically the point is: # if test_bit == 0 it means that a cloud it's likely been detected. # The return value is written this way to be consistent with the way @@ -46,188 +63,220 @@ def restoral_flag(bits): @define(kw_only=True, slots=True) class CloudTests(object): - """ - """ - data: xr.Dataset = field(validator=[validators.instance_of(xr.Dataset), ]) - scene_name: str = field(validator=[validators.instance_of(str), - validators.in_(_scene_list)]) - thresholds: str = field(validator=[validators.instance_of(Dict), ]) - scene_idx: tuple = field(init=False, - default=Factory(lambda self: tuple(np.nonzero(self.data[self.scene_name] == 1)), - takes_self=True), - validator=[validators.instance_of(tuple), ]) - pixels_in_scene: bool = field(init=False, - default=Factory(lambda self: self.scene_idx[0].shape[0] != 0, - takes_self=True), - validator=[validators.instance_of(bool), ]) + """ """ + + data: xr.Dataset = field( + validator=[ + validators.instance_of(xr.Dataset), + ] + ) + scene_name: str = field( + validator=[validators.instance_of(str), validators.in_(_scene_list)] + ) + thresholds: str = field( + validator=[ + validators.instance_of(Dict), + ] + ) + scene_idx: tuple = field( + init=False, + default=Factory( + lambda self: tuple(np.nonzero(self.data[self.scene_name] == 1)), + takes_self=True, + ), + validator=[ + validators.instance_of(tuple), + ], + ) + pixels_in_scene: bool = field( + init=False, + default=Factory(lambda self: self.scene_idx[0].shape[0] != 0, takes_self=True), + validator=[ + validators.instance_of(bool), + ], + ) def run_if_test_exists_for_scene(test_name): def decorate(func): - @functools.wraps(func) def wrapper(self, *args, **kwargs): if test_name not in self.thresholds[self.scene_name]: # test_bit = np.zeros(args[-1].shape) return args[-2], args[-1] else: - kwargs['confidence'] = np.ones(self.data.latitude.shape) + kwargs["confidence"] = np.ones(self.data.latitude.shape) # kwargs['qa_bit'] = np.zeros(self.data.latitude.shape) # kwargs['test_bit'] = np.zeros(self.data.latitude.shape) - kwargs['thresholds'] = self.thresholds[self.scene_name][test_name] + kwargs["thresholds"] = self.thresholds[self.scene_name][test_name] return func(self, *args, **kwargs) return wrapper - return decorate - def compute_test_bits(self, - rad: np.ndarray, - thr: np.ndarray, - ineq: str, - scene_idx=None) -> np.ndarray: + return decorate + def compute_test_bits( + self, rad: np.ndarray, thr: np.ndarray, ineq: str, scene_idx=None + ) -> np.ndarray: scene_idx = scene_idx or self.scene_idx rad = rad.ravel() - if ineq == '<=': - idx = np.nonzero((rad <= thr) & (self.data[self.scene_name].values[scene_idx] == 1)) - if ineq == '<': - idx = np.nonzero((rad < thr) & (self.data[self.scene_name].values[scene_idx] == 1)) - if ineq == '>': - idx = np.nonzero((rad > thr) & (self.data[self.scene_name].values[scene_idx] == 1)) - if ineq == '>=': - idx = np.nonzero((rad >= thr) & (self.data[self.scene_name].values[scene_idx] == 1)) + if ineq == "<=": + idx = np.nonzero( + (rad <= thr) & (self.data[self.scene_name].values[scene_idx] == 1) + ) + if ineq == "<": + idx = np.nonzero( + (rad < thr) & (self.data[self.scene_name].values[scene_idx] == 1) + ) + if ineq == ">": + idx = np.nonzero( + (rad > thr) & (self.data[self.scene_name].values[scene_idx] == 1) + ) + if ineq == ">=": + idx = np.nonzero( + (rad >= thr) & (self.data[self.scene_name].values[scene_idx] == 1) + ) test_bits = np.zeros(rad.shape) test_bits[idx] = 1 return test_bits - @run_if_test_exists_for_scene('11um_Test') - def test_11um(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - + @run_if_test_exists_for_scene("11um_Test") + def test_11um( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: # logger.info(f'Running test 11um_Test on {self.scene_name}') - threshold = kwargs['thresholds'] + threshold = kwargs["thresholds"] - if (threshold['perform'] is True and self.pixels_in_scene is True): + if threshold["perform"] is True and self.pixels_in_scene is True: logger.info(f'Running test 11um_Test on "{self.scene_name}"\n') - bits['qa'][self.scene_idx] = 1 + bits["qa"][self.scene_idx] = 1 rad = self.data[band].values[self.scene_idx] # idx = np.nonzero((self.data[band].values >= threshold['thr'][1]) & # (self.data[self.scene_name] == 1)) - bits['test'][self.scene_idx] = self.compute_test_bits(rad, threshold['thr'][1], '>=') - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad, threshold['thr']) + bits["test"][self.scene_idx] = self.compute_test_bits( + rad, threshold["thr"][1], ">=" + ) + kwargs["confidence"][self.scene_idx] = conf.conf_test_new( + rad, threshold["thr"] + ) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) - logger.info(f'Test 11um_Test performed on {self.scene_name}') - return cmin, bits # restoral_flag(kwargs) + logger.info(f"Test 11um_Test performed on {self.scene_name}") + return cmin, bits # restoral_flag(bits) # THIS FUNCTION DOES NOT HAVE THE TEST_BIT - @run_if_test_exists_for_scene('Surface_Temperature_Test') - def surface_temperature_test(self, - band: str, - viirs_data: xr.Dataset, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 - logger.info(f'Running test Surface_Temperature_Test on "{self.scene_name}"\n') + @run_if_test_exists_for_scene("Surface_Temperature_Test") + def surface_temperature_test( + self, band: str, viirs_data: xr.Dataset, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 + logger.info( + f'Running test Surface_Temperature_Test on "{self.scene_name}"\n' + ) # print(f'Testing "{self.scene_name}"\n') rad = self.data[band].values[self.scene_idx] sfcdif = viirs_data.geos_sfct.values[self.scene_idx] - rad - thr = preproc.thresholds_surface_temperature(viirs_data, threshold, self.scene_idx) + thr = preproc.thresholds_surface_temperature( + viirs_data, threshold, self.scene_idx + ) # idx = np.nonzero((sfcdif < thr[1, :]) & (self.data[self.scene_name] == 1)) # kwargs['test_bit'][idx] = 1 - bits['test'][self.scene_idx] = self.compute_test_bits(sfcdif, thr[:, 1], '<') - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(sfcdif, thr) + bits["test"][self.scene_idx] = self.compute_test_bits( + sfcdif, thr[:, 1], "<" + ) + kwargs["confidence"][self.scene_idx] = conf.conf_test_new(sfcdif, thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('SST_Test') - def sst_test(self, - band31: str, - band32: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: + @run_if_test_exists_for_scene("SST_Test") + def sst_test( + self, band31: str, band32: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): + if threshold["perform"] is True and self.pixels_in_scene is True: logger.info(f'Running test SST_Test on "{self.scene_name}"\n') - bits['qa'][self.scene_idx] = 1 + bits["qa"][self.scene_idx] = 1 m31 = self.data[band31].values - 273.16 m32 = self.data[band32].values - 273.16 bt_diff = m31 - m32 # bt_diff = self.data[band31].values - self.data[band32].values sst = self.data.sst.values - 273.16 - cosvza = np.cos(self.data.sensor_zenith.values*_DTR) + cosvza = np.cos(self.data.sensor_zenith.values * _DTR) - c = threshold['coeffs'] + c = threshold["coeffs"] - modsst = 273.16 + c[0] + c[1]*m31 + c[2]*bt_diff*sst + c[3]*bt_diff*((1.0/cosvza) - 1.0) + modsst = ( + 273.16 + + c[0] + + c[1] * m31 + + c[2] * bt_diff * sst + + c[3] * bt_diff * ((1.0 / cosvza) - 1.0) + ) sfcdif = self.data.sst.values - modsst - # idx = np.nonzero((sfcdif < threshold['thr'][1]) & (self.data[self.scene_name] == 1)) + # idx = np.nonzero((sfcdif < threshold['thr'][1]) & + # (self.data[self.scene_name] == 1)) # kwargs['test_bit'][idx] = 1 - bits['test'][self.scene_idx] = self.compute_test_bits(sfcdif[self.scene_idx], - threshold['thr'][1], '<') - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(sfcdif[self.scene_idx], threshold['thr']) + bits["test"][self.scene_idx] = self.compute_test_bits( + sfcdif[self.scene_idx], threshold["thr"][1], "<" + ) + kwargs["confidence"][self.scene_idx] = conf.conf_test_new( + sfcdif[self.scene_idx], threshold["thr"] + ) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) - logger.info(f'Test 11um_Test performed on {self.scene_name}') + logger.info(f"Test 11um_Test performed on {self.scene_name}") return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('8.6-11um_Test') - def bt_diff_86_11um(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: + @run_if_test_exists_for_scene("8.6-11um_Test") + def bt_diff_86_11um( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 # print(f'Testing "{self.scene_name}"\n') logger.info(f'Running test 8.6-11um_Test on "{self.scene_name}"\n') rad = self.data[band].values[self.scene_idx] # idx = np.nonzero((self.data[band].values < threshold['thr'][1]) & # (self.data[self.scene_name] == 1)) # kwargs['test_bit'][idx] = 1 - bits['test'][self.scene_idx] = self.compute_test_bits(rad, threshold['thr'][1], '<') - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad, threshold['thr']) + bits["test"][self.scene_idx] = self.compute_test_bits( + rad, threshold["thr"][1], "<" + ) + kwargs["confidence"][self.scene_idx] = conf.conf_test_new( + rad, threshold["thr"] + ) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('11-12um_Cirrus_Test') - def test_11_12um_diff(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 - thr = preproc.thresholds_11_12um(self.data, threshold, self.scene_name, self.scene_idx) + @run_if_test_exists_for_scene("11-12um_Cirrus_Test") + def test_11_12um_diff( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 + thr = preproc.thresholds_11_12um( + self.data, threshold, self.scene_name, self.scene_idx + ) # print(f'Testing "{self.scene_name}"\n') rad = self.data[band].values[self.scene_idx] logger.info(f'Running test 11-12um_Cirrus_Test on "{self.scene_name}"\n') @@ -235,73 +284,83 @@ class CloudTests(object): # tmp_bit = kwargs['test_bit'][self.scene_idx] # tmp_bit[idx] = 1 # kwargs['test_bit'][self.scene_idx] = tmp_bit - if self.scene_name == 'Ocean_Day': - comparison = '<' + if self.scene_name == "Ocean_Day": + comparison = "<" else: - comparison = '<=' - bits['test'][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], comparison) - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad, thr) + comparison = "<=" + bits["test"][self.scene_idx] = self.compute_test_bits( + rad, thr[1, :], comparison + ) + kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr) - cmin = np.fmin(cmin, kwargs['confidence']) - bits['test'] = bits['test'].reshape(self.data[band].shape) + cmin = np.fmin(cmin, kwargs["confidence"]) + bits["test"] = bits["test"].reshape(self.data[band].shape) return cmin, bits # restoral_flag(kwargs) # THIS NEEDS TO BE TESTED - @run_if_test_exists_for_scene('11-4um_BT_Difference_Test_Ocean') - def bt_difference_11_4um_test_ocean(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 + @run_if_test_exists_for_scene("11-4um_BT_Difference_Test_Ocean") + def bt_difference_11_4um_test_ocean( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 # print(f'Testing "{self.scene_name}"\n') - logger.info(f'Running test 11-4um_BT_Difference_Test_Ocean on "{self.scene_name}"\n') - thr = preproc.bt_diff_11_4um_thresholds(self.data, threshold, self.scene_idx) + logger.info( + f'Running test 11-4um_BT_Difference_Test_Ocean on "{self.scene_name}"\n' + ) + thr = preproc.bt_diff_11_4um_thresholds( + self.data, threshold, self.scene_idx + ) rad = self.data[band].values[self.scene_idx] # idx = np.nonzero((rad >= thr[1, :]) & (self.data[self.scene_name] == 1)) # tmp_bit = kwargs['test_bit'][self.scene_idx] # tmp_bit[idx] = 1 # kwargs['test_bit'][self.scene_idx] = tmp_bit - bits['test'][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], '<') - kwargs['confidence'][self.scene_idx] = conf.conf_test_dble(rad, thr) + bits["test"][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], "<") + kwargs["confidence"][self.scene_idx] = conf.conf_test_dble(rad, thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) # THIS NEEDS TO BE TESTED - @run_if_test_exists_for_scene('11-4um_BT_Difference_Test_Land') - def bt_difference_11_4um_test_land(self, - band1: str, - band2: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - + @run_if_test_exists_for_scene("11-4um_BT_Difference_Test_Land") + def bt_difference_11_4um_test_land( + self, band1: str, band2: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: # confidence = np.ones(self.data.M15.shape) # qa_bit = np.zeros(self.data[band1].shape) # test_bit = np.zeros(self.data[band1].shape) - threshold = kwargs['thresholds'] - test_name = '11-4um_BT_Difference_Test_Land' # this is temporary + threshold = kwargs["thresholds"] + test_name = "11-4um_BT_Difference_Test_Land" # this is temporary - scene_flag = scn.find_scene(self.data, self.thresholds['Sun_Glint']['bounds'][3]) + scene_flag = scn.find_scene( + self.data, self.thresholds["Sun_Glint"]["bounds"][3] + ) - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 # print(f'Testing "{self.scene_name}"\n') - logger.info(f'Running test 11-4um_BT_Difference_Test_Land on "{self.scene_name}"\n') - if self.scene_name == 'Land_Night': + logger.info( + f'Running test 11-4um_BT_Difference_Test_Land on "{self.scene_name}"\n' + ) + # some scenes are not implemented... + if self.scene_name == "Land_Night": australia = self.data.Australia.values ndvi = self.data.ndvi.values - coast = scene_flag['coast'] - sh_lake = scene_flag['sh_lake'] - idx = np.where(((australia == 0) & (ndvi < threshold['ndvi']) & (coast == 0)) | - ((australia == 1) & (ndvi < threshold['ndvi_australia']) & (coast == 0))) + coast = scene_flag["coast"] + sh_lake = scene_flag["sh_lake"] + idx = np.where( + ((australia == 0) & (ndvi < threshold["ndvi"]) & (coast == 0)) + | ( + (australia == 1) + & (ndvi < threshold["ndvi_australia"]) + & (coast == 0) + ) + ) not_idx = np.ones(australia.shape, dtype=bool) not_idx[idx] = False @@ -313,52 +372,76 @@ class CloudTests(object): thr[:3, idx_lake] = thr[:3, idx_lake] + 2 temp_thr = np.squeeze(thr[:, idx[0]]) temp_rad = rad11_4[idx] - bits['test'][idx] = self.compute_test_bits(temp_rad[idx], thr[1, idx], '<=') - kwargs['confidence'][idx] = conf.conf_test_new(temp_rad, temp_thr) + bits["test"][idx] = self.compute_test_bits( + temp_rad[idx], thr[1, idx], "<=" + ) + kwargs["confidence"][idx] = conf.conf_test_new(temp_rad, temp_thr) thr = preproc.land_night_thresholds(self.data, threshold, coast=True) temp_thr = np.squeeze(thr[:, not_idx]) temp_rad = rad11_4[not_idx] - thr[0:3, idx_lake[0], idx_lake[1]] = thr[0:3, idx_lake[0], idx_lake[1]] - 1 - thr[3:6, idx_lake[0], idx_lake[1]] = thr[3:6, idx_lake[0], idx_lake[1]] + 1 - - kwargs['confidence'][not_idx] = conf.conf_test_dble(rad11_4[not_idx], thr[:, not_idx]) - - if self.scene_name in ['Polar_Night_Land', 'Polar_Day_Snow', 'Day_Snow']: + thr[0:3, idx_lake[0], idx_lake[1]] = ( + thr[0:3, idx_lake[0], idx_lake[1]] - 1 + ) + thr[3:6, idx_lake[0], idx_lake[1]] = ( + thr[3:6, idx_lake[0], idx_lake[1]] + 1 + ) + + kwargs["confidence"][not_idx] = conf.conf_test_dble( + rad11_4[not_idx], thr[:, not_idx] + ) + + if self.scene_name in ["Polar_Night_Land", "Polar_Day_Snow", "Day_Snow"]: # rad11_12 = self.data[band1].values[self.scene_idx] rad11_4 = self.data[band2].values[self.scene_idx] - thr = preproc.polar_night_thresholds(self.data, threshold, self.scene_name, - test_name, self.scene_idx) - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad11_4, thr) + thr = preproc.polar_night_thresholds( + self.data, threshold, self.scene_name, test_name, self.scene_idx + ) + kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad11_4, thr) - if self.scene_name in ['Polar_Night_Snow']: + if self.scene_name in ["Polar_Night_Snow"]: # rad11_12 = self.data[band1].values rad11_4 = self.data[band2].values - idx_top = np.nonzero(self.data.M15.values > threshold['bt_thr']) - thr_top = preproc.polar_night_thresholds(self.data, threshold['top'], self.scene_name, - test_name, idx_top) - idx_bottom = np.nonzero((self.data.M15.values > threshold['bt_thr']) & - (self.data.geos_tpw > threshold['tpw_thr'])) - thr_bottom = preproc.polar_night_thresholds(self.data, threshold['bottom'], self.scene_name, - test_name, idx_top) - - kwargs['confidence'][idx_top] = conf.conf_test_new(rad11_4[idx_top], thr_top) - kwargs['confidence'][idx_bottom] = conf.conf_test_new(rad11_4[idx_bottom], thr_bottom) - - if self.scene_name in ['Night_Snow']: + idx_top = np.nonzero(self.data.M15.values > threshold["bt_thr"]) + thr_top = preproc.polar_night_thresholds( + self.data, threshold["top"], self.scene_name, test_name, idx_top + ) + idx_bottom = np.nonzero( + (self.data.M15.values > threshold["bt_thr"]) + & (self.data.geos_tpw > threshold["tpw_thr"]) + ) + thr_bottom = preproc.polar_night_thresholds( + self.data, threshold["bottom"], self.scene_name, test_name, idx_top + ) + + kwargs["confidence"][idx_top] = conf.conf_test_new( + rad11_4[idx_top], thr_top + ) + kwargs["confidence"][idx_bottom] = conf.conf_test_new( + rad11_4[idx_bottom], thr_bottom + ) + + if self.scene_name in ["Night_Snow"]: # rad11_12 = self.data[band1].values rad11_4 = self.data[band2].values - idx_top = np.nonzero(self.data.M15.values > threshold['bt_thr']) - thr_top = threshold['thr'] - idx_bottom = np.nonzero((self.data.M15.values > threshold['bt_thr']) & - (self.data.geos_tpw > threshold['tpw_thr'])) - thr_bottom = preproc.polar_night_thresholds(self.data, threshold['bottom'], self.scene_name, - test_name, idx_top) - - kwargs['confidence'][idx_top] = conf.conf_test_new(rad11_4[idx_top], thr_top) - kwargs['confidence'][idx_bottom] = conf.conf_test_new(rad11_4[idx_bottom], thr_bottom) - - cmin = np.fmin(cmin, kwargs['confidence']) + idx_top = np.nonzero(self.data.M15.values > threshold["bt_thr"]) + thr_top = threshold["thr"] + idx_bottom = np.nonzero( + (self.data.M15.values > threshold["bt_thr"]) + & (self.data.geos_tpw > threshold["tpw_thr"]) + ) + thr_bottom = preproc.polar_night_thresholds( + self.data, threshold["bottom"], self.scene_name, test_name, idx_top + ) + + kwargs["confidence"][idx_top] = conf.conf_test_new( + rad11_4[idx_top], thr_top + ) + kwargs["confidence"][idx_bottom] = conf.conf_test_new( + rad11_4[idx_bottom], thr_bottom + ) + + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) @@ -372,105 +455,135 @@ class CloudTests(object): def water_vapor_cloud_test(): pass - @run_if_test_exists_for_scene('11um_Variability_Test') - def variability_11um_test(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] + @run_if_test_exists_for_scene("11um_Variability_Test") + def variability_11um_test( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] var = np.zeros([self.data.M15.shape[0], self.data.M15.shape[1], 9]) radshape = self.data.M15.shape - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 # print(f'Testing "{self.scene_name}"\n') logger.info(f'Running test 11um_Variability_Test on "{self.scene_name}"\n') rad = self.data[band].values - test = (sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - - np.expand_dims(rad, (2, 3))) - var[np.abs(test).reshape(radshape[0], radshape[1], 9) < threshold['variability']] = 1 + test = sliding_window_view( + np.pad(rad, [1, 1], mode="constant"), (3, 3) + ) - np.expand_dims(rad, (2, 3)) + var[ + np.abs(test).reshape(radshape[0], radshape[1], 9) + < threshold["variability"] + ] = 1 var = var.sum(axis=2) idx = np.where(var == 9) - bits['test'][idx] = 1 - thr = np.array(threshold['thr']) - kwargs['confidence'][idx] = conf.conf_test_new(rad[idx], thr) + bits["test"][idx] = 1 + thr = np.array(threshold["thr"]) + kwargs["confidence"][idx] = conf.conf_test_new(rad[idx], thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) # THIS FUNCTION NEEDS THE TEST_BIT DOUBLE CHECKED - @run_if_test_exists_for_scene('11-4um_Oceanic_Stratus_Test') - def oceanic_stratus_11_4um_test(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 + @run_if_test_exists_for_scene("11-4um_Oceanic_Stratus_Test") + def oceanic_stratus_11_4um_test( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 # print(f'Testing "{self.scene_name}"\n') - logger.info(f'Running test 11-4um_Oceanic_Stratus_Test on "{self.scene_name}"\n') + logger.info( + f'Running test 11-4um_Oceanic_Stratus_Test on "{self.scene_name}"\n' + ) # rad = self.data[band].values[self.scene_idx] rad = self.data[band].values - thr = threshold['thr'] - - if self.scene_name in ['Ocean_Day']: - scene_flag = scn.find_scene(self.data, self.thresholds['Sun_Glint']['bounds'][3]) - idx = tuple(np.nonzero((self.data[band].values >= threshold['thr'][1]) & - (self.data[self.scene_name] == 1) & - (scene_flag['sunglint'] == 0))) - bits['test'][idx] = 1 - idx = tuple(np.nonzero((self.data[self.scene_name] == 1) & - (scene_flag['sunglint'] == 0))) - kwargs['confidence'][idx] = conf.conf_test_new(rad[idx], thr) - - elif self.scene_name in ['Land_Day', 'Land_Day_Coast']: - idx = tuple(np.nonzero((self.data[band].values >= threshold['thr'][1]) & - (self.data[self.scene_name] == 1))) - bits['test'][idx] = 1 + thr = threshold["thr"] + + if self.scene_name in ["Ocean_Day"]: + scene_flag = scn.find_scene( + self.data, self.thresholds["Sun_Glint"]["bounds"][3] + ) + idx = tuple( + np.nonzero( + (self.data[band].values >= threshold["thr"][1]) + & (self.data[self.scene_name] == 1) + & (scene_flag["sunglint"] == 0) + ) + ) + bits["test"][idx] = 1 + idx = tuple( + np.nonzero( + (self.data[self.scene_name] == 1) + & (scene_flag["sunglint"] == 0) + ) + ) + kwargs["confidence"][idx] = conf.conf_test_new(rad[idx], thr) + + elif self.scene_name in ["Land_Day", "Land_Day_Coast"]: + idx = tuple( + np.nonzero( + (self.data[band].values >= threshold["thr"][1]) + & (self.data[self.scene_name] == 1) + ) + ) + bits["test"][idx] = 1 idx = tuple(np.nonzero(self.data[self.scene_name] == 1)) - kwargs['confidence'][idx] = conf.conf_test_new(rad[idx], thr) - - elif self.scene_name in ['Land_Day_Desert', 'Land_Day_Desert_Coast', 'Polar_Day_Desert', - 'Polar_Day_Desert_Coast']: - idx = tuple(np.nonzero((self.data[band].values >= threshold['thr'][4]) & - (self.data[band].values <= threshold['thr'][1]) & - (self.data.M15.values < threshold['bt_cutoff']) & - (self.data[self.scene_name] == 1))) - bits['test'][idx] = 1 - kwargs['confidence'][self.scene_idx] = conf.conf_test_dble(rad[self.scene_idx], thr) + kwargs["confidence"][idx] = conf.conf_test_new(rad[idx], thr) + + elif self.scene_name in [ + "Land_Day_Desert", + "Land_Day_Desert_Coast", + "Polar_Day_Desert", + "Polar_Day_Desert_Coast", + ]: + idx = tuple( + np.nonzero( + (self.data[band].values >= threshold["thr"][4]) + & (self.data[band].values <= threshold["thr"][1]) + & (self.data.M15.values < threshold["bt_cutoff"]) + & (self.data[self.scene_name] == 1) + ) + ) + bits["test"][idx] = 1 + kwargs["confidence"][self.scene_idx] = conf.conf_test_dble( + rad[self.scene_idx], thr + ) else: - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad[self.scene_idx], thr) - - # these scenes have not been implemented yet... and it might not be necessary - # elif self.scene_name in ['Land_Night', 'Polar_Night_Land', 'Polar_Day_Snow', 'Polar_Night_Snow', + kwargs["confidence"][self.scene_idx] = conf.conf_test_new( + rad[self.scene_idx], thr + ) + + # these scenes have not been implemented yet... and it might not be + # necessary + # elif self.scene_name in ['Land_Night', 'Polar_Night_Land', + # 'Polar_Day_Snow', 'Polar_Night_Snow', # 'Day_Snow', 'Night_Snow', 'Antarctic_Day']: # pass - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('NIR_Reflectance_Test') - def nir_reflectance_test(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 - thr = preproc.thresholds_NIR(self.data, self.thresholds, self.scene_name, - 'NIR_Reflectance_Test', self.scene_idx) + @run_if_test_exists_for_scene("NIR_Reflectance_Test") + def nir_reflectance_test( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 + thr = preproc.thresholds_NIR( + self.data, + self.thresholds, + self.scene_name, + "NIR_Reflectance_Test", + self.scene_idx, + ) # print(f'Testing "{self.scene_name}"\n') logger.info(f'Running test NIR_Reflectance_Test on "{self.scene_name}"\n') rad = self.data[band].values[self.scene_idx] @@ -479,168 +592,200 @@ class CloudTests(object): # tmp_bit = kwargs['test_bit'][self.scene_idx] # tmp_bit[idx] = 1 # kwargs['test_bit'][self.scene_idx] = tmp_bit - bits['test'][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], '<=') - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad, thr) + bits["test"][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], "<=") + kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('Vis/NIR_Ratio_Test') - def vis_nir_ratio_test(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] + @run_if_test_exists_for_scene("Vis/NIR_Ratio_Test") + def vis_nir_ratio_test( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] temp_conf = np.ones(self.data[band].shape) - if (threshold['perform'] is True and self.pixels_in_scene is True): - # I'm not using the self.scene_idx here because it messes up with the indexing of the - # confidence array and I don't want to think of a solution at the moment. I will need - # to figure out the logic to make it work once I'm at the stage where I want to optimize + if threshold["perform"] is True and self.pixels_in_scene is True: + # I'm not using the self.scene_idx here because it messes up with the + # indexing of the confidence array and I don't want to think of a solution + # at the moment. I will need to figure out the logic to make it work once + # I'm at the stage where I want to optimize # the code logger.info(f'Running test Vis/NIR_Ratio_Test on "{self.scene_name}"\n') - bits['qa'][self.scene_idx] = 1 + bits["qa"][self.scene_idx] = 1 rad = self.data[band].values # [self.scene_idx] solar_zenith = self.data.solar_zenith.values # [self.scene_idx] - sunglint = scn.find_scene(self.data, self.thresholds['Sun_Glint']['bounds'][3])['sunglint'] + sunglint = scn.find_scene( + self.data, self.thresholds["Sun_Glint"]["bounds"][3] + )["sunglint"] idx = np.nonzero((solar_zenith <= 85) & (sunglint == 1)) - thr_no_sunglint = np.array([threshold['thr'][i] for i in range(8)]) - thr_sunglint = np.array([self.thresholds['Sun_Glint']['snglnt'][i] for i in range(8)]) + thr_no_sunglint = np.array([threshold["thr"][i] for i in range(8)]) + thr_sunglint = np.array( + [self.thresholds["Sun_Glint"]["snglnt"][i] for i in range(8)] + ) # print(thr_no_sunglint) # print(thr_sunglint) # tmp = self.thresholds['Sun_Glint']['snglnt'] - # thr_sunglint = np.array([tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], 1]) - # kwargs['confidence'][sunglint == 0] = conf.conf_test_dble(rad[sunglint == 0], thr_no_sunglint) + # thr_sunglint = np.array( + # [tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], 1] + # ) + # kwargs['confidence'][sunglint == 0] = conf.conf_test_dble( + # rad[sunglint == 0], thr_no_sunglint + # ) # kwargs['confidence'][idx] = conf.conf_test_dble(rad[idx], thr_sunglint) - temp_conf[sunglint == 0] = conf.conf_test_dble(rad[sunglint == 0], thr_no_sunglint) + temp_conf[sunglint == 0] = conf.conf_test_dble( + rad[sunglint == 0], thr_no_sunglint + ) temp_conf[idx] = conf.conf_test_dble(rad[idx], thr_sunglint) - kwargs['confidence'][self.scene_idx] = temp_conf[self.scene_idx] - - idx = np.nonzero(((rad < thr_no_sunglint[1]) | (rad > thr_no_sunglint[4])) & - (self.data[self.scene_name].values == 1) & (sunglint == 0)) - bits['test'][tuple(idx)] = 1 - - idx = np.nonzero(((rad < thr_sunglint[1]) | (rad > thr_sunglint[4])) & - ((self.data[self.scene_name].values == 1) & - (solar_zenith <= 85) & (sunglint == 1))) - bits['test'][tuple(idx)] = 1 - - cmin = np.fmin(cmin, kwargs['confidence']) + kwargs["confidence"][self.scene_idx] = temp_conf[self.scene_idx] + + idx = np.nonzero( + ((rad < thr_no_sunglint[1]) | (rad > thr_no_sunglint[4])) + & (self.data[self.scene_name].values == 1) + & (sunglint == 0) + ) + bits["test"][tuple(idx)] = 1 + + idx = np.nonzero( + ((rad < thr_sunglint[1]) | (rad > thr_sunglint[4])) + & ( + (self.data[self.scene_name].values == 1) + & (solar_zenith <= 85) + & (sunglint == 1) + ) + ) + bits["test"][tuple(idx)] = 1 + + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('1.6_2.1um_NIR_Reflectance_Test') - def test_16_21um_reflectance(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 - thr = preproc.thresholds_NIR(self.data, self.thresholds, self.scene_name, - '1.6_2.1um_NIR_Reflectance_Test', self.scene_idx) + @run_if_test_exists_for_scene("1.6_2.1um_NIR_Reflectance_Test") + def test_16_21um_reflectance( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 + thr = preproc.thresholds_NIR( + self.data, + self.thresholds, + self.scene_name, + "1.6_2.1um_NIR_Reflectance_Test", + self.scene_idx, + ) # print(f'Testing "{self.scene_name}"\n') - logger.info(f'Running test 1.6_2.1um_NIR_Reflectance_Test on "{self.scene_name}"\n') + logger.info( + f'Running test 1.6_2.1um_NIR_Reflectance_Test on "{self.scene_name}"\n' + ) rad = self.data[band].values[self.scene_idx] # idx = np.nonzero((rad <= thr[1, :]) & # (self.data[self.scene_name].values[self.scene_idx] == 1)) # tmp_bit = kwargs['test_bit'][self.scene_idx] # tmp_bit[idx] = 1 # kwargs['test_bit'][self.scene_idx] = tmp_bit - bits['test'][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], '<=') - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad, thr) + bits["test"][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], "<=") + kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) # this test needs some TLC - @run_if_test_exists_for_scene('Visible_Reflectance_Test') - def visible_reflectance_test(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - logger.warning('The pytest seems to be working now but it used to fail for this function. ' + - 'If it fails again, go check the test before turning this function inside out') - - if (threshold['perform'] is True and self.pixels_in_scene is True): + @run_if_test_exists_for_scene("Visible_Reflectance_Test") + def visible_reflectance_test( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + logger.warning( + "The pytest seems to be working now but it used to fail for this function. " + + "If it fails again, go check the test before turning this " + + "function inside out" + ) + + if threshold["perform"] is True and self.pixels_in_scene is True: # print(f'Testing "{self.scene_name}"\n')a - bits['qa'][self.scene_idx] = 1 - logger.info(f'Running test Visible_Reflectance_Test on "{self.scene_name}"\n') - thr, rad = preproc.vis_refl_thresholds(self.data, self.thresholds, self.scene_name, - self.scene_idx) + bits["qa"][self.scene_idx] = 1 + logger.info( + f'Running test Visible_Reflectance_Test on "{self.scene_name}"\n' + ) + thr, rad = preproc.vis_refl_thresholds( + self.data, self.thresholds, self.scene_name, self.scene_idx + ) - # idx = np.nonzero((rad <= thr[1, :]) & (self.data[self.scene_name].values[self.scene_idx] == 1)) + # idx = np.nonzero((rad <= thr[1, :]) & + # (self.data[self.scene_name].values[self.scene_idx] == 1)) # tmp_bit = kwargs['test_bit'][self.scene_name] # tmp_bit[idx] = 1 # kwargs['test_bit'][self.scene_idx] = tmp_bit - bits['test'][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], '<=') - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad, thr) + bits["test"][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], "<=") + kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('GEMI_Test') - def gemi_test(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - scene_idx = tuple(np.nonzero((self.data[self.scene_name] == 1) & (self.data['ndvi'] >= 0.1))) - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][scene_idx] = 1 - thr = preproc.gemi_thresholds(self.data, threshold, self.scene_name, scene_idx) + @run_if_test_exists_for_scene("GEMI_Test") + def gemi_test( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + scene_idx = tuple( + np.nonzero((self.data[self.scene_name] == 1) & (self.data["ndvi"] >= 0.1)) + ) + + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][scene_idx] = 1 + thr = preproc.gemi_thresholds( + self.data, threshold, self.scene_name, scene_idx + ) rad = self.data[band].values[scene_idx] logger.info(f'Running test GEMI_Test on "{self.scene_name}"\n') - # idx = np.nonzero((rad > thr[1, :]) & (self.data[self.scene_name].values[self.scene_idx] == 1)) + # idx = np.nonzero((rad > thr[1, :]) & + # (self.data[self.scene_name].values[self.scene_idx] == 1)) # tmp_bit = kwargs['test_bit'][self.scene_name] # tmp_bit[idx] = 1 # kwargs['test_bit'][self.scene_idx] = tmp_bit - bits['test'][scene_idx] = self.compute_test_bits(rad, thr[1, :], '>', scene_idx=scene_idx) - kwargs['confidence'][scene_idx] = conf.conf_test_new(rad, thr) + bits["test"][scene_idx] = self.compute_test_bits( + rad, thr[1, :], ">", scene_idx=scene_idx + ) + kwargs["confidence"][scene_idx] = conf.conf_test_new(rad, thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('1.38um_High_Cloud_Test') - def test_1_38um_high_clouds(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - scene_flag = scn.find_scene(self.data, self.thresholds['Sun_Glint']['bounds'][3]) - if self.scene_name in ['Ocean_Day', 'Polar_Day_Ocean']: - bits['qa'][self.scene_idx] = 1 - thr = preproc.thresholds_1_38um_test(self.data, self.thresholds, self.scene_name, - self.scene_idx) + @run_if_test_exists_for_scene("1.38um_High_Cloud_Test") + def test_1_38um_high_clouds( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + + if threshold["perform"] is True and self.pixels_in_scene is True: + scene_flag = scn.find_scene( + self.data, self.thresholds["Sun_Glint"]["bounds"][3] + ) + if self.scene_name in ["Ocean_Day", "Polar_Day_Ocean"]: + bits["qa"][self.scene_idx] = 1 + thr = preproc.thresholds_1_38um_test( + self.data, self.thresholds, self.scene_name, self.scene_idx + ) scene_idx = self.scene_idx else: # return cmin, kwargs['test_bit'] - scene_idx = tuple(np.nonzero((self.data[self.scene_name] == 1) & - (scene_flag['high_elevation'] == 0))) - bits['qa'][scene_idx] = 1 - thr = np.full((len(scene_idx[0]), 4), threshold['thr'][:4]).T + scene_idx = tuple( + np.nonzero( + (self.data[self.scene_name] == 1) + & (scene_flag["high_elevation"] == 0) + ) + ) + bits["qa"][scene_idx] = 1 + thr = np.full((len(scene_idx[0]), 4), threshold["thr"][:4]).T # print(f'Testing "{self.scene_name}"\n') logger.info(f'Running test 1.38um_High_Cloud_Test on "{self.scene_name}"\n') @@ -650,48 +795,50 @@ class CloudTests(object): # tmp_bit = kwargs['test_bit'][self.scene_idx] # tmp_bit[idx] = 1 # kwargs['test_bit'][self.scene_idx] = tmp_bit - bits['test'][scene_idx] = self.compute_test_bits(rad, thr[1, :], '<=', scene_idx=scene_idx) - kwargs['confidence'][scene_idx] = conf.conf_test_new(rad, thr) + bits["test"][scene_idx] = self.compute_test_bits( + rad, thr[1, :], "<=", scene_idx=scene_idx + ) + kwargs["confidence"][scene_idx] = conf.conf_test_new(rad, thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) - @run_if_test_exists_for_scene('4-12um_BTD_Thin_Cirrus_Test') - def thin_cirrus_4_12um_BTD_test(self, - band: str, - cmin: np.ndarray, - bits: Dict, - **kwargs) -> np.ndarray: - - threshold = kwargs['thresholds'] - - if (threshold['perform'] is True and self.pixels_in_scene is True): - bits['qa'][self.scene_idx] = 1 - thr = preproc.polar_night_thresholds(self.data, self.thresholds, self.scene_name, - '4-12um_BTD_Thin_Cirrus_Test', self.scene_idx) + @run_if_test_exists_for_scene("4-12um_BTD_Thin_Cirrus_Test") + def thin_cirrus_4_12um_BTD_test( + self, band: str, cmin: np.ndarray, bits: Dict, **kwargs + ) -> np.ndarray: + threshold = kwargs["thresholds"] + + if threshold["perform"] is True and self.pixels_in_scene is True: + bits["qa"][self.scene_idx] = 1 + thr = preproc.polar_night_thresholds( + self.data, + self.thresholds, + self.scene_name, + "4-12um_BTD_Thin_Cirrus_Test", + self.scene_idx, + ) rad = self.data[band].values[self.scene_idx] - logger.info(f'Running test 4-12um_BTD_Thin_Cirrus_Test on "{self.scene_name}"\n') + logger.info( + f'Running test 4-12um_BTD_Thin_Cirrus_Test on "{self.scene_name}"\n' + ) # idx = np.nonzero((rad <= thr[1, :]) & # (self.data[self.scene_name].values[self.scene_idx] == 1)) # tmp_bit = kwargs['test_bit'][self.scene_idx] # tmp_bit[idx] = 1 # kwargs['test_bit'][self.scene_idx] = tmp_bit - bits['test'][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], '<=') - kwargs['confidence'][self.scene_idx] = conf.conf_test_new(rad, thr) + bits["test"][self.scene_idx] = self.compute_test_bits(rad, thr[1, :], "<=") + kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad, thr) - cmin = np.fmin(cmin, kwargs['confidence']) + cmin = np.fmin(cmin, kwargs["confidence"]) return cmin, bits # restoral_flag(kwargs) class ComputeTests(CloudTests): - - def __init__(self, - data: xr.Dataset, - scene_name: str, - thresholds: Dict) -> None: + def __init__(self, data: xr.Dataset, scene_name: str, thresholds: Dict) -> None: super().__init__(data, scene_name, thresholds) def run_tests(self): @@ -702,49 +849,74 @@ class ComputeTests(CloudTests): cmin_G5 = np.ones(self.data.M01.shape) # Group 1 - cmin_G1 = self.test_11um('M15', cmin_G1, test_name='11um_Test') - cmin_G1 = self.sst_test('M15', 'M16', cmin_G1, test_name='SST_Test') + cmin_G1 = self.test_11um("M15", cmin_G1, test_name="11um_Test") + cmin_G1 = self.sst_test("M15", "M16", cmin_G1, test_name="SST_Test") # Group 2 - cmin_G2 = self.bt_diff_86_11um('M14-M15', cmin_G2, test_name='8.6-11um_Test') - cmin_G2 = self.test_11_12um_diff('M15-M16', cmin_G2, test_name='11-12um_Cirrus_Test') - cmin_G2 = self.oceanic_stratus_11_4um_test('M15-M13', cmin_G2, - test_name='11-4um_Oceanic_Stratus_Test') + cmin_G2 = self.bt_diff_86_11um("M14-M15", cmin_G2, test_name="8.6-11um_Test") + cmin_G2 = self.test_11_12um_diff( + "M15-M16", cmin_G2, test_name="11-12um_Cirrus_Test" + ) + cmin_G2 = self.oceanic_stratus_11_4um_test( + "M15-M13", cmin_G2, test_name="11-4um_Oceanic_Stratus_Test" + ) # Group 3 - cmin_G3 = self.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test') - cmin_G3 = self.vis_nir_ratio_test('M07-M05ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test') - cmin_G3 = self.nir_reflectance_test('M10', cmin_G3, test_name='1.6_2.1um_NIR_Reflectance_Test') - cmin_G3 = self.visible_reflectance_test('M128', cmin_G3, test_name='Visible_Reflectance_Test') - cmin_G3 = self.gemi_test('GEMI', cmin_G3, test_name='GEMI_Test') + cmin_G3 = self.nir_reflectance_test( + "M07", cmin_G3, test_name="NIR_Reflectance_Test" + ) + cmin_G3 = self.vis_nir_ratio_test( + "M07-M05ratio", cmin_G3, test_name="Vis/NIR_Ratio_Test" + ) + cmin_G3 = self.nir_reflectance_test( + "M10", cmin_G3, test_name="1.6_2.1um_NIR_Reflectance_Test" + ) + cmin_G3 = self.visible_reflectance_test( + "M128", cmin_G3, test_name="Visible_Reflectance_Test" + ) + cmin_G3 = self.gemi_test("GEMI", cmin_G3, test_name="GEMI_Test") # Group 4 - cmin_G4 = self.test_1_38um_high_clouds('M09', cmin_G4, test_name='1.38um_High_Cloud_Test') + cmin_G4 = self.test_1_38um_high_clouds( + "M09", cmin_G4, test_name="1.38um_High_Cloud_Test" + ) # Group 5 - cmin_G5 = self.thin_cirrus_4_12um_BTD_test('M13-M16', cmin_G5, - test_name='4-12um_BTD_Thin_Cirrus_Test') + cmin_G5 = self.thin_cirrus_4_12um_BTD_test( + "M13-M16", cmin_G5, test_name="4-12um_BTD_Thin_Cirrus_Test" + ) cmin = cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5 return cmin - def clear_sky_restoral(self, - cmin: np.ndarray) -> np.ndarray: - + def clear_sky_restoral(self, cmin: np.ndarray) -> np.ndarray: total_bit = np.full(self.data.M01.shape, 3) - sunglint_angle = self.thresholds['Sun_Glint']['bounds'][3] + sunglint_angle = self.thresholds["Sun_Glint"]["bounds"][3] scene_flags = scn.find_scene(self.data, sunglint_angle) cmin_tmp = cmin + 0 - idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['ice'] == 0) & - (scene_flags['uniform'] == 1) & (cmin <= 0.99) & (cmin >= 0.05)) - cmin[idx] = restoral.spatial(self.data, self.thresholds['Sun_Glint'], scene_flags, cmin_tmp)[idx] + idx = np.nonzero( + (scene_flags["water"] == 1) + & (scene_flags["ice"] == 0) + & (scene_flags["uniform"] == 1) + & (cmin <= 0.99) + & (cmin >= 0.05) + ) + cmin[idx] = restoral.spatial( + self.data, self.thresholds["Sun_Glint"], scene_flags, cmin_tmp + )[idx] cmin_tmp = cmin + 0 - idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['sunglint'] == 1) & - (scene_flags['uniform'] == 1) & (cmin <= 0.95)) - cmin[idx] = restoral.sunglint(self.data, self.thresholds['Sun_Glint'], total_bit, cmin_tmp)[idx] + idx = np.nonzero( + (scene_flags["water"] == 1) + & (scene_flags["sunglint"] == 1) + & (scene_flags["uniform"] == 1) + & (cmin <= 0.95) + ) + cmin[idx] = restoral.sunglint( + self.data, self.thresholds["Sun_Glint"], total_bit, cmin_tmp + )[idx] """ idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) & @@ -759,15 +931,3 @@ class ComputeTests(CloudTests): """ return cmin - - - - - - - - - - - - diff --git a/mvcm/utility_functions.py b/mvcm/utility_functions.py index c53f5a8..0306b54 100644 --- a/mvcm/utility_functions.py +++ b/mvcm/utility_functions.py @@ -1,66 +1,87 @@ -import numpy as np +from typing import Dict +import numpy as np from numpy.lib.stride_tricks import sliding_window_view -from typing import Dict - -def local_nxn_mean(arr: np.ndarray, - x: int = 3, - y: int = 3) -> np.ndarray: +def local_nxn_mean(arr: np.ndarray, x: int = 3, y: int = 3) -> np.ndarray: """Returns the mean value of each NxM (default 3) pixel region for a given array""" - local_mean = sliding_window_view(np.pad(arr, [x-2, y-2], mode='constant'), - (x, y)).reshape(arr.shape[0], arr.shape[1], x*y).mean(axis=2) + local_mean = ( + sliding_window_view(np.pad(arr, [x - 2, y - 2], mode="constant"), (x, y)) + .reshape(arr.shape[0], arr.shape[1], x * y) + .mean(axis=2) + ) return local_mean -def local_nxn_standard_deviation(arr: np.ndarray, - x: int = 3, - y: int = 3) -> np.ndarray: +def local_nxn_standard_deviation(arr: np.ndarray, x: int = 3, y: int = 3) -> np.ndarray: """Returns the standard deviation of each NxM pixel region for a given array""" - local_std = sliding_window_view(np.pad(arr, [x-2, y-2], mode='constant'), - (x, y)).reshape(arr.shape[0], arr.shape[1], x*y).std(axis=2) + local_std = ( + sliding_window_view(np.pad(arr, [x - 2, y - 2], mode="constant"), (x, y)) + .reshape(arr.shape[0], arr.shape[1], x * y) + .std(axis=2) + ) return local_std -def spatial_var(rad: np.ndarray, - threshold: Dict) -> np.ndarray: +def spatial_var(rad: np.ndarray, threshold: Dict) -> np.ndarray: """Compute 3x3 spatial variability and returns fractional number of pixels within the given threshold""" - test = sliding_window_view(np.pad(rad, [1, 1], mode='constant'), (3, 3)) - np.expand_dims(rad, (2, 3)) + test = sliding_window_view( + np.pad(rad, [1, 1], mode="constant"), (3, 3) + ) - np.expand_dims(rad, (2, 3)) test = np.abs(test.reshape(rad.shape[0], rad.shape[1], 9)) var = np.zeros(test.shape) var[test <= threshold] = 1 - return var.sum(axis=2)/9 + return var.sum(axis=2) / 9 -def compute_ndvi(b01: np.ndarray, - b02: np.ndarray) -> np.ndarray: +def compute_ndvi(b01: np.ndarray, b02: np.ndarray) -> np.ndarray: """Compute Normalized Difference Vegetation Index (NDVI)""" - ndvi = (b02 - b01)/(b01 + b02) + ndvi = (b02 - b01) / (b01 + b02) return ndvi def bt11_elevation_correction(height: np.ndarray) -> np.ndarray: """Compute the elevation correction for the 11um brightness temperature thresholds""" - return 5*height/1000 + return 5 * height / 1000 def group_count(bits: Dict) -> int: - g1 = (np.any(bits['01'] == 1) or np.any(bits['02'] == 1) or np.any(bits['03'] == 1)) + 0 - - g2 = (np.any(bits['04'] == 1) or np.any(bits['05'] == 1) or np.any(bits['06'] == 1) or - np.any(bits['07'] == 1) or np.any(bits['08'] == 1) or np.any(bits['09'] == 1)) + 0 - - g3 = (np.any(bits['10'] == 1) or np.any(bits['11'] == 1) or np.any(bits['12'] == 1) or - np.any(bits['13'] == 1) or np.any(bits['14'] == 1)) + 0 - - g4 = np.any(bits['15'] == 1) + 0 - - g5 = np.any(bits['16'] == 1) + 0 - - return np.fmax(1, g1+g2+g3+g4+g5) - - - + g1 = ( + np.any(bits["01"] == 1) or np.any(bits["02"] == 1) or np.any(bits["03"] == 1) + ) + 0 + + g2 = ( + np.any(bits["04"] == 1) + or np.any(bits["05"] == 1) + or np.any(bits["06"] == 1) + or np.any(bits["07"] == 1) + or np.any(bits["08"] == 1) + or np.any(bits["09"] == 1) + ) + 0 + + g3 = ( + np.any(bits["10"] == 1) + or np.any(bits["11"] == 1) + or np.any(bits["12"] == 1) + or np.any(bits["13"] == 1) + or np.any(bits["14"] == 1) + ) + 0 + + g4 = np.any(bits["15"] == 1) + 0 + + g5 = np.any(bits["16"] == 1) + 0 + + return np.fmax(1, g1 + g2 + g3 + g4 + g5) + + +def restoral_flag(bits: np.ndarray) -> np.ndarray: + qa_bit = bits["qa"] + test_bit = bits["test"] + # the logic of this is convoluted as hell, but basically the point is: + # if test_bit == 0 it means that a cloud it's likely been detected. + # The return value is written this way to be consistent with the way + # this is defined in the C code, at some point I might change it + return np.logical_not(np.logical_and(qa_bit, np.logical_not(test_bit))) diff --git a/mvcm/utils.py b/mvcm/utils.py index 96c15dd..26ea2ab 100644 --- a/mvcm/utils.py +++ b/mvcm/utils.py @@ -14,52 +14,54 @@ _test_thr = [5, 10, 15, 1, 1] def sunglint_scene(refang, sunglint_thr): sunglint_flag = np.zeros(refang.shape) refang = np.abs(refang) - sunglint_flag[(refang > sunglint_thr['bounds'][2]) & (refang <= sunglint_thr['bounds'][3])] = 1 - sunglint_flag[(refang > sunglint_thr['bounds'][1]) & (refang <= sunglint_thr['bounds'][2])] = 2 - sunglint_flag[refang <= sunglint_thr['bounds'][1]] = 3 + sunglint_flag[ + (refang > sunglint_thr["bounds"][2]) & (refang <= sunglint_thr["bounds"][3]) + ] = 1 + sunglint_flag[ + (refang > sunglint_thr["bounds"][1]) & (refang <= sunglint_thr["bounds"][2]) + ] = 2 + sunglint_flag[refang <= sunglint_thr["bounds"][1]] = 3 return sunglint_flag def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr): - - band = f'band{band_n}' + band = f"band{band_n}" sunglint_thr = np.zeros((4, len(refang))) -# if refang > thresholds['bounds'][3]: -# sunglint = sunglint -# # dosgref[2] = hicnf -# # dosgref[0] = locnf -# # dosgref[1] = mdcnf -# # sunglint[3] = doref2[3] + # if refang > thresholds['bounds'][3]: + # sunglint = sunglint + # # dosgref[2] = hicnf + # # dosgref[0] = locnf + # # dosgref[1] = mdcnf + # # sunglint[3] = doref2[3] # if refang <= thresholds['bounds'][1]: if sunglint_flag == 3: - sunglint_thr = np.full((len(refang), 4), thresholds[f'{band}_0deg']).T + sunglint_thr = np.full((len(refang), 4), thresholds[f"{band}_0deg"]).T # sunglint_thr[:, :] = thresholds[f'{band}_0deg'] else: - # if (refang > thresholds['bounds'][1] and refang <= thresholds['bounds'][2]): if sunglint_flag == 2: - lo_ang = thresholds['bounds'][1] - hi_ang = thresholds['bounds'][2] - lo_ang_val = thresholds[f'{band}_10deg'][0] - hi_ang_val = thresholds[f'{band}_10deg'][1] - power = thresholds[f'{band}_10deg'][3] - conf_range = thresholds[f'{band}_10deg'][2] + lo_ang = thresholds["bounds"][1] + hi_ang = thresholds["bounds"][2] + lo_ang_val = thresholds[f"{band}_10deg"][0] + hi_ang_val = thresholds[f"{band}_10deg"][1] + power = thresholds[f"{band}_10deg"][3] + conf_range = thresholds[f"{band}_10deg"][2] # elif (refang > thresholds['bounds'][2] and refang <= thresholds['bounds'][3]): elif sunglint_flag == 1: - lo_ang = thresholds['bounds'][2] - hi_ang = thresholds['bounds'][3] - lo_ang_val = thresholds[f'{band}_20deg'][0] + lo_ang = thresholds["bounds"][2] + hi_ang = thresholds["bounds"][3] + lo_ang_val = thresholds[f"{band}_20deg"][0] hi_ang_val = thr[1] - power = thresholds[f'{band}_20deg'][3] - conf_range = thresholds[f'{band}_20deg'][2] + power = thresholds[f"{band}_20deg"][3] + conf_range = thresholds[f"{band}_20deg"][2] else: raise ValueError("Wrong sunglint flag") a = (refang - lo_ang) / (hi_ang - lo_ang) - midpt = lo_ang_val + a*(hi_ang_val - lo_ang_val) + midpt = lo_ang_val + a * (hi_ang_val - lo_ang_val) sunglint_thr[1, :] = midpt sunglint_thr[2, :] = midpt - conf_range sunglint_thr[0, :] = midpt + conf_range @@ -69,7 +71,7 @@ def get_sunglint_thresholds(refang, thresholds, band_n, sunglint_flag, thr): def conf_test(rad=_test_rad, thr=_test_thr): - ''' + """ Assuming a linear function between min and max confidence level, the plot below shows how the confidence (y axis) is computed as function of radiance (x axis). This case illustrates alpha < gamma, obviously in case alpha > gamma, the plot would be @@ -87,7 +89,7 @@ def conf_test(rad=_test_rad, thr=_test_thr): e 0________/ | alpha --------- radiance ----------> - ''' + """ coeff = np.power(2, (thr[3] - 1)) hicut = thr[0] @@ -106,18 +108,18 @@ def conf_test(rad=_test_rad, thr=_test_thr): flipped = True # Rad between alpha and beta - range_ = 2. * (beta - alpha) + range_ = 2.0 * (beta - alpha) s1 = (rad[rad <= beta] - alpha) / range_ if flipped is False: confidence[rad <= beta] = coeff * np.power(s1, power) if flipped is True: - confidence[rad <= beta] = 1. - (coeff * np.power(s1, power)) + confidence[rad <= beta] = 1.0 - (coeff * np.power(s1, power)) # Rad between beta and gamma - range_ = 2. * (beta - gamma) + range_ = 2.0 * (beta - gamma) s1 = (rad[rad > beta] - gamma) / range_ if flipped is False: - confidence[rad > beta] = 1. - (coeff * np.power(s1, power)) + confidence[rad > beta] = 1.0 - (coeff * np.power(s1, power)) if flipped is True: confidence[rad > beta] = coeff * np.power(s1, power) @@ -166,12 +168,11 @@ def conf_test_dble(rad, coeffs): coeff = np.power(2, (power - 1)) radshape = rad.shape - rad = rad.reshape((rad.shape[0]*rad.shape[1])) + rad = rad.reshape((rad.shape[0] * rad.shape[1])) c = np.zeros(rad.shape) # Find if interval between inner cutoffs passes or fails test - if (alpha1 - gamma1 > 0): - + if alpha1 - gamma1 > 0: # Value is within range of lower set of limits range_ = 2 * (beta1 - alpha1) s1 = (rad[(rad <= alpha1) & (rad >= beta1)] - alpha1) / range_ @@ -196,7 +197,6 @@ def conf_test_dble(rad, coeffs): c[(rad < gamma1) | (rad > gamma2)] = 1 else: - # Value is withing range of lower set of limits range_ = 2 * (beta1 - alpha1) s1 = (rad[(rad <= gamma1) & (rad <= beta1)] - alpha1) / range_ @@ -239,23 +239,23 @@ def get_b1_threshold(thresholds, ndvi_thresholds): coeffs = np.full((10, 3, 4, len(pxin_ndvibk))) - indvi = np.zeros((len(pxin_ndvibk), )) - x = np.zeros((len(pxin_ndvibk), )) - y1 = np.zeros((len(pxin_ndvibk), )) - y2 = np.zeros((len(pxin_ndvibk), )) - interp = np.zeros((len(pxin_ndvibk), )) + indvi = np.zeros((len(pxin_ndvibk),)) + x = np.zeros((len(pxin_ndvibk),)) + y1 = np.zeros((len(pxin_ndvibk),)) + y2 = np.zeros((len(pxin_ndvibk),)) + interp = np.zeros((len(pxin_ndvibk),)) # Check for special cases and fill values ('ndvibk'=32.767) # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) & # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) & # (pxin_ndvibk < ndvi_thresholds['ndvi_bnd1'])) - idx = np.nonzero(pxin_ndvibk < ndvi_thresholds['ndvi_bnd1']) + idx = np.nonzero(pxin_ndvibk < ndvi_thresholds["ndvi_bnd1"]) # all gets set to zero, so no need to do anything for this idx # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) & # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) & # (pxin_ndvibk >= ndvi_thresholds['ndvi_bnd2'])) - idx = np.nonzero(pxin_ndvibk >= ndvi_thresholds['ndvi_bnd2']) + idx = np.nonzero(pxin_ndvibk >= ndvi_thresholds["ndvi_bnd2"]) indvi[idx] = 9 # else @@ -263,12 +263,14 @@ def get_b1_threshold(thresholds, ndvi_thresholds): # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) & # (pxin_ndvibk >= ndvi_thresholds['ndvi_bnd1']) & # (pxin_ndvibk < ndvi_thresholds['ndvi_bnd2'])) - idx = np.nonzero((pxin_ndvibk >= ndvi_thresholds['ndvi_bnd1']) & - (pxin_ndvibk < ndvi_thresholds['ndvi_bnd2'])) + idx = np.nonzero( + (pxin_ndvibk >= ndvi_thresholds["ndvi_bnd1"]) + & (pxin_ndvibk < ndvi_thresholds["ndvi_bnd2"]) + ) interp[idx] = 1 # this is not really needed anymore, maybe indvi = (pxin_ndvibk[idx] / delta_ndvi_bin) - 0.5 indvi[indvi < 0] = 0 - x = (pxin_ndvibk - delta_ndvi_bin*indvi + delta_ndvi_bin/2) / delta_ndvi_bin + x = (pxin_ndvibk - delta_ndvi_bin * indvi + delta_ndvi_bin / 2) / delta_ndvi_bin x[x < 0] = 0 x[x > 1] = 1 @@ -276,22 +278,30 @@ def get_b1_threshold(thresholds, ndvi_thresholds): # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) & # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1])) - y1[:, i] = (coeffs[indvi, i, 0] + - coeffs[indvi, i, 1]*pxin_sctang + - coeffs[indvi, i, 2]*pxin_sctang**2 + - coeffs[indvi, i, 3]*pxin_sctang**3) + y1[:, i] = ( + coeffs[indvi, i, 0] + + coeffs[indvi, i, 1] * pxin_sctang + + coeffs[indvi, i, 2] * pxin_sctang**2 + + coeffs[indvi, i, 3] * pxin_sctang**3 + ) # idx = np.nonzero((pxin_ndvibk < ndvi_thresholds['fill_ndvi'][0]) & # (pxin_ndvibk > ndvi_thresholds['fill_ndvi'][1]) & interp == 1) - y2[idx, i] = (coeffs[indvi[idx]+1, i, 0] + - coeffs[indvi[idx]+1, i, 1]*pxin_sctang[idx] + - coeffs[indvi[idx]+1, i, 2]*pxin_sctang[idx]**2 + - coeffs[indvi[idx]+1, i, 3]*pxin_sctang[idx]**3) + y2[idx, i] = ( + coeffs[indvi[idx] + 1, i, 0] + + coeffs[indvi[idx] + 1, i, 1] * pxin_sctang[idx] + + coeffs[indvi[idx] + 1, i, 2] * pxin_sctang[idx] ** 2 + + coeffs[indvi[idx] + 1, i, 3] * pxin_sctang[idx] ** 3 + ) thr[:, i] = ((1.0 - x) * y1) + (x * y2) - thr_adj[pxin_ndvibk < des_ndvi, i] = thr[pxin_ndvibk < des_ndvi, 0]*ndvi_thresholds['adj_fac_desert'] - thr_adj[pxin_ndvibk >= des_ndvi, i] = thr[pxin_ndvibk >= des_ndvi, 0]*ndvi_thresholds['adj_fac_land'] + thr_adj[pxin_ndvibk < des_ndvi, i] = ( + thr[pxin_ndvibk < des_ndvi, 0] * ndvi_thresholds["adj_fac_desert"] + ) + thr_adj[pxin_ndvibk >= des_ndvi, i] = ( + thr[pxin_ndvibk >= des_ndvi, 0] * ndvi_thresholds["adj_fac_land"] + ) thr[:, 3] = 2.0 @@ -299,23 +309,24 @@ def get_b1_threshold(thresholds, ndvi_thresholds): def find_scene(data): - scene = {'ocean_day': 0, - 'ocean_night': 0, - 'land_day': 0, - 'land_night': 0, - 'snow_day': 0, - 'coast_day': 0, - 'desert_day': 0, - 'antarctic_day': 0, - 'polar_day_snow': 0, - 'polar_day_desert': 0, - 'polar_day_desert_coast': 0, - 'polar_day_coast': 0, - 'polar_day_land': 0, - 'polar_night_snow': 0, - 'polar_night_land': 0, - 'polar_ocean_night': 0, - } + scene = { + "ocean_day": 0, + "ocean_night": 0, + "land_day": 0, + "land_night": 0, + "snow_day": 0, + "coast_day": 0, + "desert_day": 0, + "antarctic_day": 0, + "polar_day_snow": 0, + "polar_day_desert": 0, + "polar_day_desert_coast": 0, + "polar_day_coast": 0, + "polar_day_land": 0, + "polar_night_snow": 0, + "polar_night_land": 0, + "polar_ocean_night": 0, + } return scene @@ -327,4 +338,3 @@ def snow_mask(viirs_data): # pred_ndvi = prd_ndvi_const[0] + 3.0*ndsi pass - diff --git a/mvcm/write_output.py b/mvcm/write_output.py index 7cd3f55..b93762d 100644 --- a/mvcm/write_output.py +++ b/mvcm/write_output.py @@ -2,29 +2,26 @@ import xarray as xr def save_output(fname): - d = { - 'varname_1': { - 'dims': (), - 'data': 0, - 'attrs': { - 'units': '', - 'long_name': '', - 'description': '', - }, + "varname_1": { + "dims": (), + "data": 0, + "attrs": { + "units": "", + "long_name": "", + "description": "", }, - - 'varname_2': { - 'dims': (), - 'data': 0, - 'attrs': { - 'units': '', - 'long_name': '', - 'description': '', - }, + }, + "varname_2": { + "dims": (), + "data": 0, + "attrs": { + "units": "", + "long_name": "", + "description": "", }, - - } + }, + } out_data = xr.Dataset.from_dict(d) out_data.to_netcdf(fname) diff --git a/pyproject.toml b/pyproject.toml index a2037e8..c2b4d79 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,6 +24,7 @@ dependencies = [ 'numpy', 'xarray', 'attrs', + 'pre-commit', ] dynamic = ['version'] @@ -47,5 +48,3 @@ build-backend = 'hatchling.build' addopts = [ '--import-mode=importlib', ] - - diff --git a/tests/test_conf.py b/tests/test_conf.py index 33a193c..3f6fc54 100644 --- a/tests/test_conf.py +++ b/tests/test_conf.py @@ -1,6 +1,6 @@ import os -import numpy as np +import numpy as np import pytest import mvcm.conf as conf @@ -8,7 +8,7 @@ import mvcm.conf as conf @pytest.fixture def fixturepath(): - return os.path.join(os.path.dirname(__file__), 'fixtures') + return os.path.join(os.path.dirname(__file__), "fixtures") @pytest.fixture @@ -28,12 +28,12 @@ def double_threshold(): @pytest.fixture def reference_data(fixturepath): - return os.path.join(fixturepath, 'ref_conf.npz') + return os.path.join(fixturepath, "ref_conf.npz") def test_single_threshold(rad, single_threshold, reference_data): - ref_confidence = np.load(reference_data)['ref_sgl'] - ref_confidence_flipped = np.load(reference_data)['ref_sgl_flipped'] + ref_confidence = np.load(reference_data)["ref_sgl"] + ref_confidence_flipped = np.load(reference_data)["ref_sgl_flipped"] c = conf.conf_test_new(rad, single_threshold) assert np.all(c == ref_confidence) @@ -44,8 +44,8 @@ def test_single_threshold(rad, single_threshold, reference_data): def test_double_threshold(rad, double_threshold, reference_data): - ref_confidence = np.load(reference_data)['ref_dbl'] - ref_confidence_flipped = np.load(reference_data)['ref_dbl_flipped'] + ref_confidence = np.load(reference_data)["ref_dbl"] + ref_confidence_flipped = np.load(reference_data)["ref_dbl_flipped"] c = conf.conf_test_dble(rad, double_threshold) print(c) diff --git a/tests/test_read_data.py b/tests/test_read_data.py index 7f99a65..6e3db30 100644 --- a/tests/test_read_data.py +++ b/tests/test_read_data.py @@ -1,25 +1,29 @@ -import numpy as np -import xarray as xr import os +import numpy as np import pytest +import xarray as xr import mvcm.read_data as rd @pytest.fixture def fixturepath(): - return os.path.join(os.path.dirname(__file__), 'fixtures') + return os.path.join(os.path.dirname(__file__), "fixtures") @pytest.fixture def l1b_file(fixturepath): - return os.path.join(fixturepath, 'VNP02MOD.A2022173.1312.001.2022174011547.uwssec.nc') + return os.path.join( + fixturepath, "VNP02MOD.A2022173.1312.001.2022174011547.uwssec.nc" + ) @pytest.fixture def geo_file(fixturepath): - return os.path.join(fixturepath, 'VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc') + return os.path.join( + fixturepath, "VNP03MOD.A2022173.1312.001.2022174012746.uwssec.nc" + ) # @pytest.fixture @@ -29,104 +33,125 @@ def geo_file(fixturepath): @pytest.fixture def sst_file(): - return 'oisst.20220622' + return "oisst.20220622" @pytest.fixture def ndvi_file(): - return 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf' + return "NDVI.FM.c004.v2.0.WS.00-04.177.hdf" @pytest.fixture def geos_file_1(): - return 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4' + return "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4" @pytest.fixture def geos_file_2(): - return 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4' + return "GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1500.V01.nc4" @pytest.fixture def geos_land(): - return 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4' + return "GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4" @pytest.fixture def geos_ocean(): - return 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4' + return "GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1330.V01.nc4" @pytest.fixture def geos_constants(): - return 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4' + return "GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4" @pytest.fixture def ref_file(fixturepath): - return os.path.join(fixturepath, 'ref_ancillary.nc') + return os.path.join(fixturepath, "ref_ancillary.nc") # this is only temporary. At some point I'll write a better test def test_l1b(fixturepath, l1b_file, geo_file): - viirs = rd.ReadData(satellite='snpp', - sensor='viirs', - file_name_l1b=l1b_file, - file_name_geo=geo_file) + viirs = rd.ReadData( + satellite="snpp", sensor="viirs", file_name_l1b=l1b_file, file_name_geo=geo_file + ) geo = viirs.read_viirs_geo() l1b = viirs.read_viirs_l1b(geo.solar_zenith.values) - geo_vars = ['latitude', 'longitude', 'height', 'range', 'solar_zenith', 'solar_azimuth', - 'sensor_zenith', 'sensor_azimuth', 'land_water_mask', 'quality_flag', - 'relative_azimuth', 'sunglint_angle', 'scattering_angle'] + geo_vars = [ + "latitude", + "longitude", + "height", + "range", + "solar_zenith", + "solar_azimuth", + "sensor_zenith", + "sensor_azimuth", + "land_water_mask", + "quality_flag", + "relative_azimuth", + "sunglint_angle", + "scattering_angle", + ] assert list(geo).sort() == geo_vars.sort() - assert 'M11' in l1b.data_vars + assert "M11" in l1b.data_vars def test_sst(fixturepath, sst_file, ref_file): - viirs = rd.ReadData(satellite='snpp', - sensor='viirs') + viirs = rd.ReadData(satellite="snpp", sensor="viirs") geo = viirs.read_viirs_geo() - ancillary = rd.ReadAncillary(latitude=geo.latitude.values, - longitude=geo.longitude.values, - resolution=1, - ancillary_dir=fixturepath, - sst_file=sst_file) + ancillary = rd.ReadAncillary( + latitude=geo.latitude.values, + longitude=geo.longitude.values, + resolution=1, + ancillary_dir=fixturepath, + sst_file=sst_file, + ) sst = ancillary.get_sst() - check_differences(ref_file, sst, 'sst') + check_differences(ref_file, sst, "sst") def test_ndvi(fixturepath, ndvi_file, ref_file): - viirs = rd.ReadData(satellite='snpp', - sensor='viirs') + viirs = rd.ReadData(satellite="snpp", sensor="viirs") geo = viirs.read_viirs_geo() - ancillary = rd.ReadAncillary(latitude=geo.latitude.values, - longitude=geo.longitude.values, - resolution=1, - ancillary_dir=fixturepath, - ndvi_file=ndvi_file) + ancillary = rd.ReadAncillary( + latitude=geo.latitude.values, + longitude=geo.longitude.values, + resolution=1, + ancillary_dir=fixturepath, + ndvi_file=ndvi_file, + ) ndvi = ancillary.get_ndvi() - check_differences(ref_file, ndvi, 'ndvi') - - -def test_geos(fixturepath, geos_file_1, geos_file_2, geos_land, geos_ocean, - geos_constants, ref_file): - viirs = rd.ReadData(satellite='snpp', - sensor='viirs') + check_differences(ref_file, ndvi, "ndvi") + + +def test_geos( + fixturepath, + geos_file_1, + geos_file_2, + geos_land, + geos_ocean, + geos_constants, + ref_file, +): + viirs = rd.ReadData(satellite="snpp", sensor="viirs") geo = viirs.read_viirs_geo() - ancillary = rd.ReadAncillary(latitude=geo.latitude.values, - longitude=geo.longitude.values, - resolution=1, - ancillary_dir=fixturepath, - geos_file_1=geos_file_1, - geos_file_2=geos_file_2, - geos_land=geos_land, - geos_ocean=geos_ocean, - geos_constants=geos_constants) + ancillary = rd.ReadAncillary( + latitude=geo.latitude.values, + longitude=geo.longitude.values, + resolution=1, + ancillary_dir=fixturepath, + geos_file_1=geos_file_1, + geos_file_2=geos_file_2, + geos_land=geos_land, + geos_ocean=geos_ocean, + geos_constants=geos_constants, + ) geos_data = ancillary.get_geos() for var in list(geos_data.keys()): - check_differences(ref_file, geos_data[var], f'geos_{var}') + check_differences(ref_file, geos_data[var], f"geos_{var}") def check_differences(ref_file, test_data, var_name): diff --git a/tests/test_restoral.py b/tests/test_restoral.py new file mode 100644 index 0000000..a88aacb --- /dev/null +++ b/tests/test_restoral.py @@ -0,0 +1,82 @@ +import os + +import numpy as np +import pytest +import xarray as xr +import yaml + +import mvcm.restoral as restoral + +# data +# thresholds +# scene_flag + + +@pytest.fixture +def fixturepath(): + return os.path.join(os.path.dirname(__file__), "fixtures") + + +@pytest.fixture +def data_path(): + return "/ships19/hercules/pveglio/mvcm_cleanup" + + +@pytest.fixture +def thresholds_file(fixturepath): + return os.path.join(fixturepath, "thresholds.mvcm.snpp.v0.0.1.yaml") + + +@pytest.fixture +def data_file(data_path): + return os.path.join(data_path, "viirs_data_A2022173.1312.nc") + + +@pytest.fixture +def thresholds(thresholds_file): + return yaml.safe_load(open(thresholds_file)) + + +@pytest.fixture +def scene_flags(data_path): + return os.path.join(data_path, "scene_flags.nc") + + +@pytest.fixture +def confidence(data_path): + return os.path.join(data_path, "ref_confidence.nc") + + +@pytest.fixture +def bits(data_path): + return os.path.join(data_path, "bits.nc") + + +# chk_spatial_var +def test_spatial_variability(): + pass + + +# chk_sunglint +def test_sunglint(): + pass + + +# chk_shallow_water +def test_shallow_water(): + pass + + +# chk_land +def test_land(): + pass + + +# chk_coast +def test_coast(): + pass + + +# chck_land_night +def test_land_night(): + pass diff --git a/tests/test_spectral_tests.py b/tests/test_spectral_tests.py index b2c18f4..37f279d 100644 --- a/tests/test_spectral_tests.py +++ b/tests/test_spectral_tests.py @@ -1,9 +1,9 @@ import os -import yaml -import xarray as xr -import numpy as np +import numpy as np import pytest +import xarray as xr +import yaml # from typing import Dict # from attrs import define, field, Factory @@ -13,28 +13,28 @@ import mvcm.spectral_tests as tst # CREATE PATHS @pytest.fixture def fixturepath(): - return os.path.join(os.path.dirname(__file__), 'fixtures') + return os.path.join(os.path.dirname(__file__), "fixtures") @pytest.fixture def data_path(): - return '/ships19/hercules/pveglio/mvcm_cleanup' + return "/ships19/hercules/pveglio/mvcm_cleanup" # SET FILENAME FIXTURES @pytest.fixture def thresholds_file(fixturepath): - return os.path.join(fixturepath, 'thresholds.mvcm.snpp.v0.0.1.yaml') + return os.path.join(fixturepath, "thresholds.mvcm.snpp.v0.0.1.yaml") @pytest.fixture def ref_confidence_file(data_path): - return os.path.join(data_path, 'ref_confidence.nc') + return os.path.join(data_path, "ref_confidence.nc") @pytest.fixture def data_file(data_path): - return os.path.join(data_path, 'viirs_data_A2022173.1312.nc') + return os.path.join(data_path, "viirs_data_A2022173.1312.nc") # SET DATA FIXTURES @@ -56,12 +56,17 @@ def ref_confidence(ref_confidence_file): def test_11um_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) + for scene_name in [ + "Ocean_Day", + "Ocean_Night", + "Polar_Day_Ocean", + "Polar_Night_Ocean", + ]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) - cmin, bit = SceneType.test_11um('M15', cmin) + cmin, bit = SceneType.test_11um("M15", cmin) assert np.allclose(cmin, ref_confidence.bt11um_confidence.values) @@ -69,52 +74,77 @@ def test_11um_test(data, thresholds, ref_confidence): def test_surface_temperature_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Land_Night', 'Polar_Night_Land']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) + for scene_name in ["Land_Night", "Polar_Night_Land"]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) - cmin, bit = SceneType.surface_temperature_test('M15', data, cmin) + cmin, bit = SceneType.surface_temperature_test("M15", data, cmin) assert np.allclose(cmin, ref_confidence.surface_temperature_confidence.values) def test_sst_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) - - cmin, bit = SceneType.sst_test('M15', 'M16', cmin) + for scene_name in [ + "Ocean_Day", + "Ocean_Night", + "Polar_Day_Ocean", + "Polar_Night_Ocean", + ]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) + + cmin, bit = SceneType.sst_test("M15", "M16", cmin) assert np.allclose(cmin, ref_confidence.sst_confidence.values) def test_bt_diff_86_11um(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) - - cmin, bit = SceneType.bt_diff_86_11um('M14-M15', cmin) + for scene_name in [ + "Ocean_Day", + "Ocean_Night", + "Polar_Day_Ocean", + "Polar_Night_Ocean", + ]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) + + cmin, bit = SceneType.bt_diff_86_11um("M14-M15", cmin) assert np.allclose(cmin, ref_confidence.diff86_11um_confidence.values) def test_11_12um_diff(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast', - 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean', - 'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert', - 'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Land_Night', 'Polar_Night_Land', - 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) - - cmin, bit = SceneType.test_11_12um_diff('M15-M16', cmin) + for scene_name in [ + "Land_Day", + "Land_Day_Coast", + "Land_Day_Desert", + "Land_Day_Desert_Coast", + "Ocean_Day", + "Ocean_Night", + "Polar_Day_Ocean", + "Polar_Night_Ocean", + "Polar_Day_Land", + "Polar_Day_Coast", + "Polar_Day_Desert", + "Polar_Day_Desert_Coast", + "Polar_Day_Snow", + "Land_Night", + "Polar_Night_Land", + "Polar_Night_Snow", + "Day_Snow", + "Night_Snow", + ]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) + + cmin, bit = SceneType.test_11_12um_diff("M15-M16", cmin) assert np.allclose(cmin, ref_confidence.diff11_12um_confidence.values) @@ -129,114 +159,146 @@ def test_bt_difference_11_4um_test_land(): def test_variability_11um_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Polar_Day_Ocean', 'Polar_Night_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) + for scene_name in ["Polar_Day_Ocean", "Polar_Night_Ocean"]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) - cmin, bit = SceneType.variability_11um_test('M15', cmin) + cmin, bit = SceneType.variability_11um_test("M15", cmin) assert np.allclose(cmin, ref_confidence.var11um_confidence.values) def test_oceanic_stratus_11_4um_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast', - 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Land', 'Polar_Day_Coast', - 'Polar_Day_Desert', 'Polar_Day_Desert_Coast', 'Polar_Day_Ocean', - 'Polar_Night_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) - - cmin, bit = SceneType.oceanic_stratus_11_4um_test('M15-M12', cmin) + for scene_name in [ + "Land_Day", + "Land_Day_Coast", + "Land_Day_Desert", + "Land_Day_Desert_Coast", + "Ocean_Day", + "Ocean_Night", + "Polar_Day_Land", + "Polar_Day_Coast", + "Polar_Day_Desert", + "Polar_Day_Desert_Coast", + "Polar_Day_Ocean", + "Polar_Night_Ocean", + ]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) + + cmin, bit = SceneType.oceanic_stratus_11_4um_test("M15-M12", cmin) assert np.allclose(cmin, ref_confidence.oceanic_stratus_11_4um_confidence.values) def test_nir_reflectance_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Ocean_Day', 'Polar_Day_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) + for scene_name in ["Ocean_Day", "Polar_Day_Ocean"]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) - cmin, bit = SceneType.nir_reflectance_test('M07', cmin) + cmin, bit = SceneType.nir_reflectance_test("M07", cmin) assert np.allclose(cmin, ref_confidence.NIR_reflectance_confidence.values) def test_vis_nir_ratio_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Ocean_Day', 'Polar_Day_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) + for scene_name in ["Ocean_Day", "Polar_Day_Ocean"]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) - cmin, bit = SceneType.vis_nir_ratio_test('M07-M05ratio', cmin) + cmin, bit = SceneType.vis_nir_ratio_test("M07-M05ratio", cmin) assert np.allclose(cmin, ref_confidence.Vis_NIR_ratio_confidence.values) def test_16_21um_reflectance_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Ocean_Day', 'Polar_Day_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) + for scene_name in ["Ocean_Day", "Polar_Day_Ocean"]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) - cmin, bit = SceneType.test_16_21um_reflectance('M10', cmin) + cmin, bit = SceneType.test_16_21um_reflectance("M10", cmin) assert np.allclose(cmin, ref_confidence.refl_16_21um_confidence.values) def test_visible_reflectance_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', - 'Land_Day_Desert_Coast', 'Polar_Day_Land', 'Polar_Day_Coast', - 'Polar_Day_Desert', 'Polar_Day_Desert_Coast']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) - - cmin, bit = SceneType.visible_reflectance_test('M128', cmin) + for scene_name in [ + "Land_Day", + "Land_Day_Coast", + "Land_Day_Desert", + "Land_Day_Desert_Coast", + "Polar_Day_Land", + "Polar_Day_Coast", + "Polar_Day_Desert", + "Polar_Day_Desert_Coast", + ]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) + + cmin, bit = SceneType.visible_reflectance_test("M128", cmin) assert np.allclose(cmin, ref_confidence.vis_reflectance_confidence.values) def test_GEMI_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Land_Day_Desert', 'Polar_Day_Desert']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) + for scene_name in ["Land_Day_Desert", "Polar_Day_Desert"]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) - cmin, bit = SceneType.gemi_test('GEMI', cmin) + cmin, bit = SceneType.gemi_test("GEMI", cmin) assert np.allclose(cmin, ref_confidence.gemi_confidence.values) def test_1_38um_high_clouds_test(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert', 'Land_Day_Desert_Coast', - 'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert', - 'Polar_Day_Desert_Coast', 'Day_Snow', 'Ocean_Day', 'Polar_Day_Ocean']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) - - cmin, bit = SceneType.test_1_38um_high_clouds('M09', cmin) + for scene_name in [ + "Land_Day", + "Land_Day_Coast", + "Land_Day_Desert", + "Land_Day_Desert_Coast", + "Polar_Day_Land", + "Polar_Day_Coast", + "Polar_Day_Desert", + "Polar_Day_Desert_Coast", + "Day_Snow", + "Ocean_Day", + "Polar_Day_Ocean", + ]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) + + cmin, bit = SceneType.test_1_38um_high_clouds("M09", cmin) assert np.allclose(cmin, ref_confidence.hi_clouds_confidence.values) def test_4_12um_thin_cirrus_btd(data, thresholds, ref_confidence): cmin = np.ones(data.latitude.shape) - for scene_name in ['Land_Night', 'Polar_Night_Land', 'Polar_Night_Snow', 'Night_Snow']: - SceneType = tst.CloudTests(data=data, - scene_name=scene_name, - thresholds=thresholds) - - cmin, bit = SceneType.thin_cirrus_4_12um_BTD_test('M13-M16', cmin) + for scene_name in [ + "Land_Night", + "Polar_Night_Land", + "Polar_Night_Snow", + "Night_Snow", + ]: + SceneType = tst.CloudTests( + data=data, scene_name=scene_name, thresholds=thresholds + ) + + cmin, bit = SceneType.thin_cirrus_4_12um_BTD_test("M13-M16", cmin) assert np.allclose(cmin, ref_confidence.thin_cirrus_confidence.values) -- GitLab