diff --git a/mvcm/ancillary.pyx b/mvcm/ancillary.pyx
index 235e3be75bc53dfa2461929fb2e3070be5b719a5..7957107e620e4785b6586de306fc97ff3c0a0e0a 100644
--- a/mvcm/ancillary.pyx
+++ b/mvcm/ancillary.pyx
@@ -1,6 +1,42 @@
# cython: language_level=3
# cython: c_string_Type=unicode, c_string_encoding=utf8
+# Define structures required for snow_mask
+cdef extern from "c_tools/include/snow_mask.h":
+ struct Radiances:
+ float m01
+ float m02
+ float m04
+ float m20
+ float m29
+ float m31
+ float nir
+
+ struct Ancillary:
+ int greenland
+ int antarctic
+ int land
+ int water
+ int sunglint
+ int high_elevation
+ float elevation
+ float surface_temperature
+
+ struct Thresholds:
+ float bt11
+ float ndsi
+ float greenland_ndsi
+ float ref2
+ float diff85_11
+ float diff85_11hel
+ float diff37_11
+ float diff37_11hel
+ float mnir
+ float lsfcdif
+ float wsfcdif
+ float bt1
+ float prd_ndvi_const
+
cdef extern void get_Reynolds_SST(float *, float *, int, int, int, char *, float *)
cdef extern void get_NDVI_background(float *, float *, int, int, int, char *, float *)
cdef extern void get_Olson_eco(float *, float *, int, int, int, char *, unsigned char *)
@@ -12,6 +48,8 @@ cdef extern void check_reg_uniformity(int, int, int, int, int, unsigned char, un
float *, float *, unsigned char *, int *, int *, int *, int *)
cdef extern int get_b1_thresholds(float, float, float *, float *, float *, float *)
cdef extern float conf_test(float, float, float, float, float)
+cdef extern float snow_mask(float, unsigned char, Radiances, Ancillary, Thresholds)
+
import cython
from cython.view cimport array as cvarray
@@ -24,6 +62,59 @@ np.import_array()
ctypedef np.float_t DTYPE_t
DTYPE = np.float32
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.initializedcheck(False)
+def py_snow_mask(np.ndarray[float, ndim=2] latitude, np.ndarray[unsigned char, ndim=2] lsf, radiances, ancillary, thresholds):
+
+ cdef Radiances rad = Radiances(0., 0., 0., 0., 0., 0., 0.)
+
+ cdef Ancillary anc = Ancillary(0, 0, 0, 0, 0, 0, 0., 0., 0.)
+
+ cdef Thresholds thr = Thresholds(
+ thresholds["bt11"],
+ thresholds["ndsi"],
+ thresholds["greenland_ndsi"],
+ thresholds["ref2"],
+ thresholds["diff85_11"],
+ thresholds["diff85_11hel"],
+ thresholds["diff37_11"],
+ thresholds["diff37_11hel"],
+ thresholds["mnir"],
+ thresholds["lsfcdif"],
+ thresholds["wsfcdif"],
+ thresholds["bt1"],
+ thresholds["prd_ndvi_const"]
+ )
+
+ cdef np.ndarray ndsi = np.zeros((latitude.shape[0], latitude.shape[1]), order='C', dtype=np.float32)
+
+ for i in range(np.shape(latitude)[0]):
+ for j in range(np.shape(latitude)[1]):
+
+ rad.m01 = radiances["m01"][i][j]
+ rad.m02 = radiances["m02"][i][j]
+ rad.m04 = radiances["m04"][i][j]
+ rad.m20 = radiances["m20"][i][j]
+ rad.m29 = radiances["m29"][i][j]
+ rad.m31 = radiances["m31"][i][j]
+ rad.nir = radiances["nir"][i][j]
+
+ anc.greenland = ancillary["greenland"][i][j]
+ anc.antarctic = ancillary["antarctic"][i][j]
+ anc.land = ancillary["land"][i][j]
+ anc.water = ancillary["water"][i][j]
+ anc.sunglint = ancillary["sunglint"][i][j]
+ anc.high_elevation = ancillary["high_elevation"][i][j]
+ anc.elevation = ancillary["elevation"][i][j]
+ anc.surface_temperature = ancillary["surface_temperature"][i][j]
+
+ ndsi[i][j] = snow_mask(latitude[i][j], lsf[i][j], rad, anc, thr)
+
+ return ndsi
+
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.initializedcheck(False)
@@ -130,15 +221,6 @@ def py_get_GEOS(np.ndarray[float, ndim=1] lat, np.ndarray[float, ndim=1] lon, in
return geos_dict
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.initializedcheck(False)
-def py_snow_mask(char *satname, unsigned char lsf):
- # need to have for loop here to compute all the pixels since the function, as with everything else,
- # is run per pixel.
- pass
-
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.initializedcheck(False)
diff --git a/mvcm/c_tools/include/snow_mask.h b/mvcm/c_tools/include/snow_mask.h
new file mode 100644
index 0000000000000000000000000000000000000000..c01b14621358e7bd66cb2d64805f3cedbb701dab
--- /dev/null
+++ b/mvcm/c_tools/include/snow_mask.h
@@ -0,0 +1,41 @@
+struct Radiances {
+ float m01;
+ float m02;
+ float m04;
+ float m20;
+ float m29;
+ float m31;
+ float m35;
+ float nir;
+};
+
+struct Ancillary {
+ int greenland;
+ int antarctic;
+ int land;
+ int water;
+ int sunglint;
+ int high_elevation;
+ float elevation;
+ float surface_temperature;
+};
+
+struct Thresholds {
+ float bt11;
+ float ndsi;
+ float greenland_ndsi;
+ float ref2;
+ float diff85_11;
+ float diff85_11hel;
+ float diff37_11;
+ float diff37_11hel;
+ float mnir;
+ float lsfcdif;
+ float wsfcdif;
+ float bt1;
+ float prd_ndvi_const;
+};
+
+typedef struct Radiances Radiances;
+typedef struct Ancillary Ancillary;
+typedef struct Thresholds Thresholds;
diff --git a/mvcm/c_tools/snow_mask.c b/mvcm/c_tools/snow_mask.c
index f7c1c212915f3a36e1d07be448f79e54eaa39b96..dfb394aae0ffa7fbabc01c1551d5db12866a90c2 100644
--- a/mvcm/c_tools/snow_mask.c
+++ b/mvcm/c_tools/snow_mask.c
@@ -46,14 +46,28 @@ Revision History:
#include <stdio.h>
#include <math.h>
-//#include "granule.h"
-//#include "pixel.h"
-//#include "mask_processing_constants.h"
-//#include "thresholds.h"
+#include "snow_mask.h"
-void snow_mask(char *satname, unsigned char lsf)
+float snow_mask(float latitude, unsigned char lsf, Radiances rad, Ancillary anc, Thresholds thresholds)
{
+ // printf("m01: %f\n", rad.m01);
+ // printf("m02: %f\n", rad.m02);
+ // printf("m04: %f\n", rad.m04);
+ // printf("m20: %f\n", rad.m20);
+ // printf("m29: %f\n", rad.m29);
+ // printf("m31: %f\n", rad.m31);
+ // printf("nir: %f\n", rad.nir);
+ //
+ // printf("greenland: %d\n", anc.greenland);
+ // printf("antarctic: %d\n", anc.antarctic);
+ // printf("land: %d\n", anc.land);
+ // printf("water: %d\n", anc.water);
+ // printf("sunglint: %d\n", anc.sunglint);
+ // printf("high_elevation: %d\n", anc.high_elevation);
+ // printf("elevation: %f\n", anc.elevation);
+ // printf("surface_temperature: %f\n", anc.surface_temperature);
+ //
// Declarations
@@ -63,15 +77,8 @@ void snow_mask(char *satname, unsigned char lsf)
int I_bad;
int hi_elev_snow;
- float m04;
- float m01;
- float m02;
- float m20;
- float m29;
- float m31;
- float m35;
- float nirband;
float ndsi;
+ float snow_from_ndsi = 0;
float ndvi;
float pred_ndvi;
float diff;
@@ -81,74 +88,47 @@ void snow_mask(char *satname, unsigned char lsf)
float sfcdif;
float bad_data = -999.0;
-// Initializations
- m04 = pxin.rad[pxin.bnd555];
- m01 = pxin.rad[pxin.bnd065];
- m02 = pxin.rad[pxin.bnd086];
- m20 = pxin.rad[pxin.bnd375];
- m29 = pxin.rad[pxin.bnd855];
- m31 = pxin.rad[pxin.bnd11];
-
diff = 0.0;
diff2 = 0.0;
hi_elev_snow = 0;
-
-// Select NIR channel for NDSI based on platform name.
-// Band 6 for Terra, band 7 for Aqua.
-// if(strcmp(satname, "MODIS_Terra") == 0) {
-// nirband = pxin.rad[pxin.bnd164];
-// }
-// else if(strcmp(satname, "MODIS_Aqua") == 0) {
-// nirband = pxin.rad[pxin.bnd213];
-// }
- if(strcmp(satname, "VIIRS_NPP") == 0) {
- nirband = pxin.rad[pxin.bnd164];
- }
- else if(strcmp(satname, "VIIRS_N20") == 0) {
- nirband = pxin.rad[pxin.bnd164];
- }
- else {
- printf("Platform name not recognized in snow_mask \n");
- }
-
I_bad = rintf(bad_data);
// Begin snow/ice detection.
// First, check the 11 micron brightness temperature.
- if (rintf(m31) != I_bad && m31 <= sm_bt11[0]) {
+ if (rintf(rad.m31) != I_bad && rad.m31 <= thresholds.bt11) {
// Get NDSI index
- if (rintf(m04) != I_bad && rintf(nirband) != I_bad) {
+ if (rintf(rad.m04) != I_bad && rintf(rad.nir) != I_bad) {
- ndsi = (m04 - nirband) / (m04 + nirband);
- ndvi = (m02 - m01) / (m02 + m01);
- pred_ndvi = prd_ndvi_const[0] + (3.0 * ndsi);
+ ndsi = (rad.m04 - rad.nir) / (rad.m04 + rad.nir);
+ ndvi = (rad.m02 - rad.m01) / (rad.m02 + rad.m01);
+ pred_ndvi = thresholds.prd_ndvi_const + (3.0 * ndsi);
- if(pxin.land && (!pxin.Greenland) && (!pxin.Antarctic) ) {
- if(pxin.elev > 1000.0) {
+ if(anc.land && (!anc.greenland) && (!anc.antarctic) ) {
+ if(anc.elevation > 1000.0) {
hi_elev_snow = 1;
- if(ndsi > 0.0) pxin.ndsi_snow = 1;
+ if(ndsi > 0.0) snow_from_ndsi = 1;
}
- else if ( (ndsi > sm_ndsi[0] && m02 > sm_ref2[0]) ||
- (ndsi > 0.1 && ndsi <= sm_ndsi[0] && ndvi <= pred_ndvi) ) {
- pxin.ndsi_snow = 1;
+ else if ( (ndsi > thresholds.ndsi && rad.m02 > thresholds.ref2) ||
+ (ndsi > 0.1 && ndsi <= thresholds.ndsi && ndvi <= pred_ndvi) ) {
+ snow_from_ndsi = 1;
}
}
- else if( pxin.land && pxin.Greenland) {
- if( ndsi > Grnlnd_ndsi[0] ) pxin.ndsi_snow = 1;
+ else if(anc.land && anc.greenland) {
+ if( ndsi > thresholds.greenland_ndsi ) snow_from_ndsi = 1;
}
else {
- if (ndsi > sm_ndsi[0] && m02 > sm_ref2[0]) pxin.ndsi_snow = 1;
+ if (ndsi > thresholds.ndsi && rad.m02 > thresholds.ref2) snow_from_ndsi = 1;
}
// printf("ndsi: %f %f %f %f %f %f %f %f %d\n", m01, m02, m04, nirband, ndsi, ndvi, prd_ndvi_const[0], pred_ndvi, pxin.ndsi_snow);
- if(pxin.ndsi_snow) {
+ if(snow_from_ndsi) {
// If in sunglint region, disregard if between -60 and 50 degrees latitude.
- if(pxin.water && pxin.sunglint && pxin.lat <= 50.0 && pxin.lat >= -60.0) {
- pxin.ndsi_snow = 0;
+ if(anc.water && anc.sunglint && latitude <= 50.0 && latitude >= -60.0) {
+ snow_from_ndsi = 0;
}
// Check for ice clouds mis-identified as snow. Note modified BTD thresholds for Greenland.
@@ -156,25 +136,25 @@ void snow_mask(char *satname, unsigned char lsf)
diff2 = -99.0;
sth85_11 = -99.0;
sth37_11 = -99.0;
- if(rintf(m29) != I_bad && rintf(m31) != I_bad) {
- diff = m29 - m31;
- sth85_11 = sm85_11[0];
- if(rintf(m20) != I_bad && !(pxin.Greenland && pxin.hi_elev)) {
- diff2 = m20 - m31;
- sth37_11 = sm37_11[0];
+ if(rintf(rad.m29) != I_bad && rintf(rad.m31) != I_bad) {
+ diff = rad.m29 - rad.m31;
+ sth85_11 = thresholds.diff85_11;
+ if(rintf(rad.m20) != I_bad && !(anc.greenland && anc.high_elevation)) {
+ diff2 = rad.m20 - rad.m31;
+ sth37_11 = thresholds.diff37_11;
if(diff >= sth85_11 && diff2 >= sth37_11) {
- pxin.ndsi_snow = 0;
+ snow_from_ndsi = 0;
}
}
else {
- if(pxin.Greenland && pxin.hi_elev) {
- sth85_11 = sm85_11hel[0];
+ if(anc.greenland && anc.high_elevation) {
+ sth85_11 = thresholds.diff85_11hel;
}
else {
- sth85_11 = sm85_11[0];
+ sth85_11 = thresholds.diff85_11;
}
if(diff >= sth85_11) {
- pxin.ndsi_snow = 0;
+ snow_from_ndsi = 0;
}
}
}
@@ -183,36 +163,36 @@ void snow_mask(char *satname, unsigned char lsf)
// Check for water clouds mis-identified as snow.
diff = -99.0;
sth37_11 = -99.0;
- if(m20 != I_bad && m31 != I_bad) {
- if( !(pxin.hi_elev && pxin.Greenland) ) {
- diff = m20 - m31;
- if(pxin.hi_elev) {
- sth37_11 = sm37_11hel[0];
+ if(rad.m20 != I_bad && rad.m31 != I_bad) {
+ if( !(anc.high_elevation && anc.greenland) ) {
+ diff = rad.m20 - rad.m31;
+ if(anc.high_elevation) {
+ sth37_11 = thresholds.diff37_11hel;
}
else {
- sth37_11 = sm37_11[0];
+ sth37_11 = thresholds.diff37_11;
}
if(diff >= sth37_11) {
- pxin.ndsi_snow = 0;
+ snow_from_ndsi = 0;
}
}
}
// Check value of NIR reflectance.
- if( !pxin.hi_elev && !hi_elev_snow) {
- if(nirband > sm_mnir[0]) {
- pxin.ndsi_snow = 0;
+ if( !anc.high_elevation && !hi_elev_snow) {
+ if(rad.nir > thresholds.mnir) {
+ snow_from_ndsi = 0;
}
// Check 11 um vs. surface temperature.
- sfcdif = pxin.sfctmp - m31;
+ sfcdif = anc.surface_temperature - rad.m31;
if(lsf == 1 || lsf == 2) {
// Land or coast.
- if(sfcdif > sm_lsfcdif[0]) pxin.ndsi_snow = 0;
+ if(sfcdif > thresholds.lsfcdif) snow_from_ndsi = 0;
}
else if(lsf == 0 || lsf > 2) {
// Water.
- if(sfcdif > sm_wsfcdif[0] && pxin.sfctmp > sm_bt1[0]) pxin.ndsi_snow = 0;
+ if(sfcdif > thresholds.wsfcdif && anc.surface_temperature > thresholds.bt1) snow_from_ndsi = 0;
}
}
@@ -222,5 +202,5 @@ void snow_mask(char *satname, unsigned char lsf)
}
}
-
+ return snow_from_ndsi;
}
diff --git a/mvcm/main.py b/mvcm/main.py
index f565bc8eeec8c2ee35ca458e942b1b23925afcdf..7cc41d4a0faefaf8a0d59f0a3f851fcaba5ee288 100644
--- a/mvcm/main.py
+++ b/mvcm/main.py
@@ -235,10 +235,10 @@ def main(
# orbit_number = f.getncattr("OrbitNumber")
attrs = {attr: f.getncattr(attr) for attr in f.ncattrs()}
- viirs_data = rd.get_data(satellite, sensor, file_names, hires=use_hires)
+ viirs_data = rd.get_data(satellite, sensor, file_names, thresholds, hires=use_hires)
sunglint_angle = thresholds["Sun_Glint"]["bounds"][3]
- scene_flags = scn.find_scene(viirs_data, sunglint_angle)
+ scene_flags = scn.find_scene(viirs_data, sunglint_angle, thresholds["Snow_Mask"])
restore = Restoral(
data=viirs_data,
@@ -331,12 +331,12 @@ def main(
cmin_g2, bits["04"] = my_scene.bt_diff_86_11um("M14-M15", cmin_g2, bits["04"])
cmin_g2, bits["05"] = my_scene.test_11_12um_diff("M15-M16", cmin_g2, bits["05"])
cmin_g2, bits["06"] = my_scene.variability_11um_test(m15_name, cmin_g2, bits["06"])
- cmin_g2, bits["07"] = my_scene.bt_difference_11_4um_test_ocean(
- "M15-M12", cmin_g2, bits["07"]
- )
- cmin_g2, bits["08"] = my_scene.bt_difference_11_4um_test_land(
- "M15-M16", "M15-M12", cmin_g2, bits["08"]
- )
+ # cmin_g2, bits["07"] = my_scene.bt_difference_11_4um_test_ocean(
+ # "M15-M12", cmin_g2, bits["07"]
+ # )
+ # cmin_g2, bits["08"] = my_scene.bt_difference_11_4um_test_land(
+ # "M15-M16", "M15-M12", cmin_g2, bits["08"]
+ # )
cmin_g2, bits["09"] = my_scene.oceanic_stratus_11_4um_test("M15-M12", cmin_g2, bits["09"])
# Group 3
diff --git a/mvcm/preprocess_thresholds.py b/mvcm/preprocess_thresholds.py
index 16066f12c27f6ea70e051bba3b8855da50e03f95..5d8f7d4a9bcdb0988e75d6608de7f4549c0ecec3 100644
--- a/mvcm/preprocess_thresholds.py
+++ b/mvcm/preprocess_thresholds.py
@@ -472,7 +472,7 @@ def polar_night_thresholds(data, thresholds, scene, test_name, scene_idx):
# get_nl_thresholds
-def land_night_thresholds(data, threshold, coast=True):
+def land_night_thresholds(m15_m16, geos_tpw, threshold, coast=True):
"""Compyte land_night_thresholds."""
if coast is False:
lo_val = threshold["bt_diff_bounds"][0]
@@ -481,20 +481,20 @@ def land_night_thresholds(data, threshold, coast=True):
hi_val_thr = threshold["thr_mid"][1]
conf_range = threshold["thr_mid"][2]
- a = (data["M15-M16"].values - lo_val) / (hi_val - lo_val)
+ a = (m15_m16 - 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])
+ idx = np.nonzero(m15_m16 > 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])
+ idx = np.nonzero(m15_m16 < threshold["bt_diff_bounds"][1])
locut[idx] = threshold["thr_hi"][0]
midpt[idx] = threshold["thr_hi"][1]
hicut[idx] = threshold["thr_hi"][2]
@@ -515,7 +515,7 @@ def land_night_thresholds(data, threshold, coast=True):
b1 = threshold["coeffs"][1]
b2 = threshold["coeffs"][2]
- thr = b0 + threshold["int_adj"] + b1 * data.geos_tpw.values + b2 * data.geos_tpw.values**2
+ thr = b0 + threshold["int_adj"] + b1 * geos_tpw + b2 * geos_tpw**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]))
diff --git a/mvcm/read_data.py b/mvcm/read_data.py
index 3459f6c0445523796634f8902eb49fe1cad965e5..80b634b9b9c86eb607d776f934797c346b9e6b8f 100644
--- a/mvcm/read_data.py
+++ b/mvcm/read_data.py
@@ -701,7 +701,9 @@ class ReadAncillary(CollectInputs):
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, thresholds: dict, hires: bool = False
+) -> xr.Dataset:
"""Store every variable in one xarray dataset.
Parameters
@@ -765,7 +767,9 @@ def get_data(satellite: str, sensor: str, file_names: dict, hires: bool = False)
viirs_data.update(Ancillary.pack_data())
- scene_flags = scn.find_scene(viirs_data, geo_data.sunglint_angle.values)
+ scene_flags = scn.find_scene(
+ viirs_data, geo_data.sunglint_angle.values, thresholds["Snow_Mask"]
+ )
scene = scn.scene_id(scene_flags)
scene_xr = xr.Dataset()
diff --git a/mvcm/scene.py b/mvcm/scene.py
index 18ce837deb1d6efea66d190ebcb6e6dec5c0cdf8..90b764835ad912762f9e1a8e02e0841684b0fc72 100644
--- a/mvcm/scene.py
+++ b/mvcm/scene.py
@@ -1,13 +1,15 @@
"""Functions that define scene type."""
-import logging
-from typing import Dict
+import logging # noqa
import numpy as np
import xarray as xr
import ancillary_data as anc
+# from typing import Dict
+
+
# lsf: land sea flag
_scene_list = [
"Ocean_Day",
@@ -121,7 +123,9 @@ def test_scene():
"""
-def find_scene(data: xr.Dataset, sunglint_angle: float) -> Dict:
+def find_scene(
+ data: xr.Dataset, sunglint_angle: float, snow_mask_thresholds: dict | None = None
+) -> dict:
"""Define basic pixel information for scene classification."""
eco = np.array(data["ecosystem"].values, dtype=np.uint8)
# logger.warning('"eco" has been renamed "ecosystem" in the most recent version of
@@ -296,6 +300,31 @@ def find_scene(data: xr.Dataset, sunglint_angle: float) -> Dict:
scene_flag["vrused"] = np.ones((dim1, dim2), dtype=np.int8)
scene_flag["vrused"][idx] = 0
+ # Define parameters to compute Snow and Ice flags
+ if snow_mask_thresholds is not None:
+ radiances = {
+ "m01": data.M05.values,
+ "m02": data.M07.values,
+ "m04": data.M04.values,
+ "m20": data.M12.values,
+ "m29": data.M14.values,
+ "m31": data.M15.values,
+ "nir": data.M10.values,
+ }
+ ancillary = {
+ "greenland": scene_flag["greenland"],
+ "antarctic": scene_flag["antarctica"],
+ "land": scene_flag["land"],
+ "water": scene_flag["water"],
+ "sunglint": scene_flag["sunglint"],
+ "high_elevation": scene_flag["high_elevation"],
+ "elevation": elev,
+ "surface_temperature": sfctmp,
+ }
+ scene_flag["ndsi_snow"] = anc.py_snow_mask(
+ lat, lsf, radiances, ancillary, snow_mask_thresholds
+ )
+
snow_fraction = data["geos_snow_fraction"]
perm_ice_fraction = data["geos_land_ice_fraction"]
ice_fraction = data["geos_ice_fraction"]
@@ -340,7 +369,6 @@ def find_scene(data: xr.Dataset, sunglint_angle: float) -> Dict:
| (lsf == 5) & (scene_flag["ndsi_snow"] == 1)
)
scene_flag["ice"][idx] = 1
-
idx = np.nonzero(
(scene_flag["day"] == 1)
& (scene_flag["water"] == 1)
@@ -353,7 +381,7 @@ def find_scene(data: xr.Dataset, sunglint_angle: float) -> Dict:
idx = np.nonzero(
(scene_flag["day"] == 1)
& (scene_flag["land"] == 1)
- & (lat >= 48.0)
+ & (lat >= -48.0)
& (lat <= -34.0)
& (lon >= 165.0)
& (lon <= 180.0)
@@ -365,9 +393,8 @@ def find_scene(data: xr.Dataset, sunglint_angle: float) -> Dict:
& (scene_flag["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["map_snow"] == 1) | (scene_flag["new_zealand"] == 1))
)
scene_flag["snow"][idx] = 1
@@ -385,7 +412,7 @@ def find_scene(data: xr.Dataset, sunglint_angle: float) -> Dict:
scene_flag["snow"][idx] = 0
idx = np.nonzero(
- (scene_flag["day"] == 0) & (scene_flag["map_snow"] == 1) & (sfctmp > 280.0) & (elev < 500.0)
+ (scene_flag["day"] == 0) & (scene_flag["map_ice"] == 1) & (sfctmp > 280.0) & (elev < 500.0)
)
scene_flag["ice"][idx] = 0
@@ -394,7 +421,6 @@ def find_scene(data: xr.Dataset, sunglint_angle: float) -> Dict:
idx = np.nonzero((lat <= -11) & (lat > -40) & (lon <= 155) & (lon > 110))
scene_flag["australia"][idx] = 1
-
# Check regional uniformity
# Check for data border pixels
# eco_out = np.array(np.empty(eco.shape), dtype=np.uint8)
diff --git a/mvcm/spectral_tests.py b/mvcm/spectral_tests.py
index 11c98959d12321d7246756d2d6983a9c3cf4b75f..999f1c69ac3b4266d4d770f4d0263011207769cb 100644
--- a/mvcm/spectral_tests.py
+++ b/mvcm/spectral_tests.py
@@ -335,41 +335,84 @@ class CloudTests:
# 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"]
+ australia = self.data.Australia.values[self.scene_idx]
+ ndvi = self.data.ndvi.values[self.scene_idx]
+ coast = scene_flag["coast"][self.scene_idx]
+ sh_lake = scene_flag["sh_lake"][self.scene_idx]
+ rad11_12 = self.data[band1].values[self.scene_idx]
+ rad11_4 = self.data[band2].values[self.scene_idx]
+ tpw = self.data.geos_tpw.values[self.scene_idx]
+
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
-
- rad11_12 = self.data[band1].values
- rad11_4 = self.data[band2].values
- thr = preproc.land_night_thresholds(self.data, threshold, coast=False)
- idx_lake = np.where((sh_lake == 1) | (rad11_12 < 0))
- # test_bit[self.scene_idx] = 1 # THIS NEEDS TO BE CALCULATED PROPERLY
+ idx_lake = np.where(
+ (
+ ((australia == 0) & (ndvi < threshold["ndvi"]) & (coast == 0))
+ | ((australia == 1) & (ndvi < threshold["ndvi_australia"]) & (coast == 0))
+ )
+ & ((sh_lake == 1) | (rad11_12 < 0))
+ )
+ print(f"idx shape: {np.shape(idx)}")
+ thr = preproc.land_night_thresholds(rad11_12[idx], tpw[idx], threshold, coast=False)
thr[:3, idx_lake] = thr[:3, idx_lake] + 2
- temp_thr = np.squeeze(thr[:, idx[0]])
- temp_rad = rad11_4[idx]
- # if np.ndim(thr) == 1:
- # thr = thr.reshape(thr.shape[0], 1)
- # bits["test"][idx] = self.compute_test_bits(
- # temp_rad[idx], thr[1, idx], "<="
- # )
- kwargs["confidence"][idx] = conf.conf_test_new(temp_rad, temp_thr)
+ print(f"thr shape: {thr.shape}")
+ if np.shape(idx[0])[0] > 0:
+ kwargs["confidence"][idx] = conf.conf_test_new(rad11_4[idx], 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
+ # tmp_idx = np.ones(australia.shape)
+ # tmp_idx[idx] = 0
+ # idx = np.where(tmp_idx == 1)
+ # idx_lake = np.where((tmp_idx == 1) & ((sh_lake == 1) | (rad11_12 < 0)))
- kwargs["confidence"][not_idx] = conf.conf_test_dble(
- rad11_4[not_idx], thr[:, not_idx]
+ idx = np.where(
+ ~(
+ ((australia == 0) & (ndvi < threshold["ndvi"]) & (coast == 0))
+ | ((australia == 1) & (ndvi < threshold["ndvi_australia"]) & (coast == 0))
+ )
)
+ idx_lake = np.where(
+ ~(
+ ((australia == 0) & (ndvi < threshold["ndvi"]) & (coast == 0))
+ | ((australia == 1) & (ndvi < threshold["ndvi_australia"]) & (coast == 0))
+ & ((sh_lake == 1) | (rad11_12 < 0))
+ )
+ )
+ print(f"idx shape: {np.shape(idx)}")
+ print(f"rad shape: {np.shape(rad11_12)}")
+ if np.shape(idx[0])[0] > 0:
+ thr = preproc.land_night_thresholds(
+ rad11_12[idx], tpw[idx], threshold, coast=True
+ )
+ thr[0:3, idx_lake] = thr[0:3, idx_lake] - 1
+ thr[3:6, idx_lake] = thr[3:6, idx_lake] + 1
+ print(f"thr shape: {thr.shape}")
+
+ kwargs["confidence"][idx] = conf.conf_test_dble(rad11_4[idx], thr[:, idx])
+
+ # not_idx = np.ones(australia.shape, dtype=bool)
+ # not_idx[idx] = False
+ # idx_lake = np.where((sh_lake == 1) | (rad11_12 < 0))
+
+ # thr = preproc.land_night_thresholds(rad11_12, tpw, threshold, coast=False)
+
+ # test_bit[self.scene_idx] = 1 # THIS NEEDS TO BE CALCULATED PROPERLY
+ # thr[:3, idx_lake] = thr[:3, idx_lake] + 2
+ # temp_thr = np.squeeze(thr[:, idx[0]])
+ # temp_rad = rad11_4[idx]
+
+ # kwargs["confidence"][idx] = conf.conf_test_new(temp_rad, temp_thr)
+
+ # thr = preproc.land_night_thresholds(rad11_12, tpw, 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"]:
# rad11_12 = self.data[band1].values[self.scene_idx]
@@ -484,9 +527,18 @@ class CloudTests:
bits["qa"][qaidx] = 1
else:
bits["qa"][self.scene_idx] = 1
+ if self.scene_name in ["Day_Snow"]:
+ thr = preproc.polar_night_thresholds(
+ self.data,
+ threshold,
+ self.scene_name,
+ "11-4um_Oceanic_Stratus_Test",
+ self.scene_idx,
+ )
+ else:
+ thr = threshold["thr"]
logger.info(f'Running test 11-4um_Oceanic_Stratus_Test on "{self.scene_name}"\n')
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])
diff --git a/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml b/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml
index 99cc058820454e00b280e11c51a5e159d636eaa6..a77fab45ee35ef3ff0b9715b984bd848cbd52a0d 100644
--- a/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml
+++ b/thresholds/thresholds.mvcm.snpp.v0.0.1.yaml
@@ -270,10 +270,10 @@ Polar_Day_Land:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 0.3
- perform: True
+ perform: False
11-4um_Oceanic_Stratus_Test:
thr: [-16.0, -14.0, -12.0, 1.0, 1.0]
- perform: True
+ perform: False
pdlh20: [215.0, 220.0, 225.0, 1.0, 1.0]
Visible_Reflectance_Test:
thr: [0.207, 0.169, 0.132, 1.0, 1.0]
@@ -339,10 +339,10 @@ Polar_Day_Coast:
coeffs: [3.0, 1.0]
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
- perform: True
+ perform: False
11-4um_Oceanic_Stratus_Test:
thr: [-16.0, -14.0, -12.0, 1.0, 1.0]
- perform: True
+ perform: False
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
Visible_Reflectance_Test:
thr: [0.207, 0.169, 0.132, 1.0, 1.0]
@@ -358,11 +358,11 @@ Polar_Day_Desert:
coeffs: [3.5, 1.0]
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
- perform: True
+ perform: False
11-4um_Oceanic_Stratus_Test:
# thr: [-22.0, -20.0, -18.0, -2.0, 0.0, 2.0, 1.0, 1.0]
thr: [2.0, 0.0, -2.0, -18.0, -20.0, -22.0, 1.0, 1.0]
- perform: True # I disable this test while I figure out if I need
+ perform: False # I disable this test while I figure out if I need
bt_cutoff: 320.0 # the bt_cutoff parameter
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
Visible_Reflectance_Test:
@@ -384,11 +384,11 @@ Polar_Day_Desert_Coast:
coeffs: [3.5, 1.0]
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
- perform: True
+ perform: False
11-4um_Oceanic_Stratus_Test:
# thr: [-23.0, -21.0, -19.0, -2.0, 0.0, 2.0, 1.0, 1.0]
thr: [2.0, 0.0, -2.0, -19.0, -21.0, -23.0, 1.0, 1.0]
- perform: True # see Polar_Day_Desert
+ perform: False # see Polar_Day_Desert
bt_cutoff: 320.0
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
Visible_Reflectance_Test:
@@ -405,21 +405,22 @@ Polar_Day_Snow:
coeffs: [3.0, 1.0]
cmult: 0 # I NEED TO WORK ON THIS
adj: 0 # I NEED TO WORK ON THIS
- perform: True
+ perform: False
dpsbt1: 230.0
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
1.38um_High_Cloud_Test:
thr: [0.060, 0.0525, 0.045, 1.0, 1.0]
tpw: 0.75
perform: True
- 11-4um_BT_Difference_Test_Land:
+ # 11-4um_BT_Difference_Test_Land:
+ 11-4um_Oceanic_Stratus_Test:
low: [20.0, 18.0, 16.0, 1.0]
mid1: [18.0, 16.0, 2.0, 1.0]
mid2: [0.0, 0.0, 0.0, 0.0]
mid3: [0.0, 0.0, 0.0, 0.0]
high: [18.0, 16.0, 14.0, 1.0]
- bt_11_bounds: [230.0, 0.0, 0.0, 245.0]
- perform: True
+ bt11_bounds: [230.0, 0.0, 0.0, 245.0]
+ perform: False
dps11_12adj: 0.8
Polar_Night_Snow:
@@ -435,7 +436,7 @@ Polar_Night_Snow:
mid2: [0.00, 0.00, 0.00, 0.0] # pn_4_12m2
mid3: [0.00, 0.00, 0.00, 0.0] # pn_4_12m3
high: [2.50, 2.00, 1.50, 1.0] # pn_4_12h
- bt11_bounds: [235.0, 0.0, 0.0, 265.0] # bt_11_bounds
+ bt11_bounds: [235.0, 0.0, 0.0, 265.0] # bt11_bounds
perform: True
7.3-11um_BTD_Mid_Level_Cloud_Test:
low: [-1.00, 0.00, 1.00, 1.0, 1.0] # pn_7_11l
@@ -479,12 +480,12 @@ Polar_Day_Ocean:
perform: True
8.6-11um_Test:
thr: [-0.50, -1.00, -1.50, 1.0, 1.0]
- perform: True
+ perform: False
11-12um_Cirrus_Test:
coeffs: [3.0, 1.0]
cmult: 0.3
adj: 1.25
- perform: True
+ perform: False
NIR_Reflectance_Test:
thr: [0.062, 0.043, 0.029, 1.0, 1.0]
coeffs: [1.7291, 0.0715, -0.0026, 0.000025889]
@@ -514,7 +515,7 @@ Polar_Day_Ocean:
perform: True
11-4um_Oceanic_Stratus_Test:
thr: [-11.0, -9.0, -7.0, 1.0, 1.0]
- perform: True
+ perform: False
Water_Vapor_High_Clouds_Test: [215.0, 220.0, 225.0, 1.0, 1.0]
pdo_b26pfm: 1.0
pdo_b2bias_adj: 0.96
@@ -574,13 +575,21 @@ Day_Snow:
thr: [0.060, 0.0525, 0.045, 1.0, 1.0]
tpw: 0.75
perform: True
- 11-4um_BT_Difference_Test_Land:
+ # 11-4um_BT_Difference_Test_Land:
+ # low: [20.0, 18.0, 16.0, 1.0]
+ # mid1: [18.0, 16.0, 2.0, 1.0]
+ # mid2: [0.0, 0.0, 0.0, 0.0]
+ # mid3: [0.0, 0.0, 0.0, 0.0]
+ # high: [18.0, 16.0, 14.0, 1.0]
+ # bt11_bounds: [230.0, 0.0, 0.0, 245.0]
+ # perform: True
+ 11-4um_Oceanic_Stratus_Test:
low: [20.0, 18.0, 16.0, 1.0]
mid1: [18.0, 16.0, 2.0, 1.0]
mid2: [0.0, 0.0, 0.0, 0.0]
mid3: [0.0, 0.0, 0.0, 0.0]
high: [18.0, 16.0, 14.0, 1.0]
- bt_11_bounds: [230.0, 0.0, 0.0, 245.0]
+ bt11_bounds: [230.0, 0.0, 0.0, 245.0]
perform: True
ds11_12adj: 0.8
ds11_12lcmult: 0.3
@@ -605,7 +614,7 @@ Night_Snow:
mid2: [0.0, 0.0, 0.0, 0.0]
mid3: [0.0, 0.0, 0.0, 0.0]
high: [18.0, 16.0, 14.0, 1.0]
- bt_11_bounds: [230.0, 0.0, 0.0, 245.0]
+ bt11_bounds: [230.0, 0.0, 0.0, 245.0]
bt_thr: 230.0
tpw_thr: 0.2
perform: True
@@ -679,18 +688,18 @@ Night_Snow_Inversion:
n65_11: [10.0, 1.0]
Snow_Mask:
- sm_bt11: 280.0
- sm_ndsi: 0.31
- Grnlnd_ndsi: 0.66
- sm_ref2: 0.1056
- sm85_11: 0.0
- sm85_11hel: 1.5
- sm37_11: 11.0
- sm37_11hel: 16.0
- sm_mnir: 0.196
- sm_lsfcdif: 20.0
- sm_wsfcdif: 20.0
- sm_bt1: 273.0
+ bt11: 280.0
+ ndsi: 0.31
+ greenland_ndsi: 0.66
+ ref2: 0.1056
+ diff85_11: 0.0
+ diff85_11hel: 1.5
+ diff37_11: 11.0
+ diff37_11hel: 16.0
+ mnir: 0.196
+ lsfcdif: 20.0
+ wsfcdif: 20.0
+ bt1: 273.0
prd_ndvi_const: -0.2015
Coastal_NDVI_Thresholds: