From d4650c0a433e0ffe462bad2691bbaf8c051aff2f Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Tue, 28 Mar 2023 20:59:41 +0000
Subject: [PATCH] started testing the code. need to correct the way the final
cmin before restoral is computed
---
mvcm/main.py | 191 +++++++++++++++++++++++++-------------
mvcm/read_data.py | 59 ++++++++++++
mvcm/restoral.py | 31 ++++---
mvcm/scene.py | 18 ++--
mvcm/utility_functions.py | 47 ++++++++++
5 files changed, 259 insertions(+), 87 deletions(-)
create mode 100644 mvcm/utility_functions.py
diff --git a/mvcm/main.py b/mvcm/main.py
index 27fa0ed..1fc73b8 100644
--- a/mvcm/main.py
+++ b/mvcm/main.py
@@ -7,28 +7,28 @@ import numpy as np
from glob import glob
-import read_data as rd
-import spectral_tests as tst
-import scene as scn
-import restoral
+import mvcm.read_data as rd
+import mvcm.spectral_tests as tst
+import mvcm.scene as scn
+from mvcm.restoral import Restoral
# #################################################################### #
# TEST CASE
# data:
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
-_fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0]
-_fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1454.001.*.uwssec.nc')[0]
-_fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0]
-_fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1454.001.*.uwssec.nc')[0]
+_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.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_1430.V01.nc4'
-_geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1430.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'
@@ -51,7 +51,8 @@ def main(satellite: str = 'snpp',
geos_ocean: str = _geos_ocean,
geos_constants: str = _geos_constants,
ndvi_file: str = _ndvi_file,
- sst_file: str = _sst_file) -> None:
+ sst_file: str = _sst_file,
+ eco_file: str = _eco_file) -> None:
"""Main function for the MVCM.
Parameters
@@ -84,8 +85,10 @@ def main(satellite: str = 'snpp',
file name of GEOS5 constants
ndvi_file: str
NDVI file name
- sst_file:
+ sst_file: str
Sea surface temperature file name
+ eco_file: str
+ Ecosystem File
Returns
-------
@@ -103,6 +106,7 @@ def main(satellite: str = 'snpp',
'GEOS_constants': f'{geos_constants}',
'NDVI': f'{ndvi_file}',
'SST': f'{sst_file}',
+ 'ECO': f'{eco_file}',
'ANC_DIR': f'{data_path}/ancillary'
}
@@ -112,67 +116,126 @@ def main(satellite: str = 'snpp',
sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
- viirs_data = rd.get_data(satellite, sensor, file_names, sunglint_angle)
+ viirs_data = rd.get_data(satellite, sensor, file_names)
# Initialize the confidence arrays for the various test groups
- cmin_G1 = np.ones(viirs_data.M01.shape)
- cmin_G2 = np.ones(viirs_data.M01.shape)
- cmin_G3 = np.ones(viirs_data.M01.shape)
- cmin_G4 = np.ones(viirs_data.M01.shape)
- cmin_G5 = np.ones(viirs_data.M01.shape)
+ cmin_G1 = np.ones(viirs_data.M11.shape)
+ cmin_G2 = np.ones(viirs_data.M11.shape)
+ cmin_G3 = np.ones(viirs_data.M11.shape)
+ cmin_G4 = np.ones(viirs_data.M11.shape)
+ cmin_G5 = np.ones(viirs_data.M11.shape)
##########################################################
- """
for scene_name in ['Land_Day', 'Land_Day_Coast', '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']:
- """
- for scene_name in ['Ocean_Day']:
- MyScene = tst.CloudTests(viirs_data, scene_name, thresholds)
-
- cmin_G1, bit1 = MyScene.test_11um('M15', cmin_G1, test_name='11um_Test')
- cmin_G1, bit2 = MyScene.sst_test('M15', 'M16', cmin_G1, test_name='SST_Test')
-
- cmin_G2, bit3 = MyScene.bt_diff_86_11um('M14-M15', cmin_G2, test_name='8.6-11um_Test')
- cmin_G2, bit4 = MyScene.test_11_12um_diff('M15-M16', cmin_G2, test_name='11-12um_Cirrus_Test')
- cmin_G2, bit5 = MyScene.oceanic_stratus_11_4um_test('M15-M12', cmin_G2,
- test_name='11-4um_Oceanic_Stratus_Test')
-
- cmin_G3, bit6 = MyScene.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test')
- cmin_G3, bit7 = MyScene.vis_nir_ratio_test('M07-M05ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test')
- cmin_G3, bit8 = MyScene.nir_reflectance_test('M10', cmin_G3, test_name='1.6_2.1um_NIR_Reflectance_Test')
- cmin_G3, bit9 = MyScene.visible_reflectance_test('M128', cmin_G3, test_name='Visible_Reflectance_Test')
- cmin_G3, bit10 = MyScene.gemi_test('GEMI', cmin_G3, test_name='GEMI_Test')
-
- cmin_G4, bit11 = MyScene.test_1_38um_high_clouds('M09', cmin_G4, test_name='1.38um_High_Cloud_Test')
-
- cmin_G5, bit12 = MyScene.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
-
- if scene_name in ['Ocean_Day']:
- # total_bit = np.full(viirs_data.M01.shape, 4)
- total_bit = bit1 + bit2 + bit4 + bit11
- sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
- scene_flags = scn.find_scene(viirs_data, sunglint_angle)
-
- # 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(viirs_data, thresholds['Sun_Glint'], scene_flags, cmin)[idx]
-
- idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['sunglint'] == 1) &
- (scene_flags['uniform'] == 1) & (cmin <= 0.95))
- cmin[idx] = restoral.sunglint(viirs_data, thresholds['Sun_Glint'], total_bit, cmin)[idx]
-
- # MVCM = tst.ComputeTests(viirs_data, scene_name, thresholds)
- # cmin = MVCM.clear_sky_restoral(cmin)
-
- # return cmin
- np.savez('test_ams', confidence=cmin, lat=viirs_data.latitude.values, lon=viirs_data.longitude.values)
+ MyScene = tst.CloudTests(data=viirs_data,
+ scene_name=scene_name,
+ thresholds=thresholds)
+
+ cmin_G1, bit01 = MyScene.test_11um('M15', cmin_G1)
+ cmin_G1, bit02 = MyScene.surface_temperature_test('M15', viirs_data, cmin_G1)
+ cmin_G1, bit03 = MyScene.sst_test('M15', 'M16', cmin_G1)
+
+ cmin_G2, bit04 = MyScene.bt_diff_86_11um('M14-M15', cmin_G2)
+ cmin_G2, bit05 = MyScene.test_11_12um_diff('M15-M16', cmin_G2)
+ cmin_G2, bit06 = MyScene.variability_11um_test('M15', cmin_G2)
+ cmin_G2, bit07 = MyScene.bt_difference_11_4um_test_ocean('M15-M12', cmin_G2)
+ # cmin_G2, bit08 = MyScene.bt_difference_11_4um_test_land('M15-M12', cmin_G2)
+ cmin_G2, bit09 = MyScene.oceanic_stratus_11_4um_test('M15-M12', cmin_G2)
+
+ cmin_G3, bit10 = MyScene.nir_reflectance_test('M07', cmin_G3)
+ cmin_G3, bit11 = MyScene.vis_nir_ratio_test('M07-M05ratio', cmin_G3)
+ cmin_G3, bit12 = MyScene.test_16_21um_reflectance('M10', cmin_G3)
+ # cmin_G3, bit13 = MyScene.visible_reflectance_test('M128', cmin_G3)
+ # cmin_G3, bit14 = MyScene.gemi_test('GEMI', cmin_G3)
+
+ cmin_G4, bit15 = MyScene.test_1_38um_high_clouds('M09', cmin_G4)
+
+ cmin_G5, bit16 = MyScene.thin_cirrus_4_12um_BTD_test('M13-M16', cmin_G5)
+
+ # I need to modify this so that it's computed for each scene and that the exponent
+ # is changed according to the number of groups in which tests are run
+ cmin = np.power(cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5, 1/5)
+
+ scene_flags = scn.find_scene(viirs_data, sunglint_angle)
+ Restore = Restoral(data=viirs_data, thresholds=thresholds, scene_flag=scene_flags)
+
+ sunglint_bits = bit01 * bit03 * bit05 * bit15
+ sh_water_bits = {'ir': bit01 * bit05,
+ 'nir_1': bit10,
+ 'nir_2': bit15,
+ 'sst': bit03
+ }
+ land_bits = bit15 * bit05
+ # make sure that the land_bits are calculated properly
+ idx = np.nonzero(scene_flags['desert'] == 0)
+ land_bits[idx] = land_bits[idx] * bit10[idx] * bit07[idx] * bit11[idx]
+ coast_bits = bit05
+ land_night_bits = bit16
+
+ 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(cmin)[idx]
+
+ idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['uniform'] == 1) &
+ (scene_flags['sunglint'] == 1) & (cmin <= 0.95))
+ cmin[idx] = Restore.sunglint(sunglint_bits, 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, 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, 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[idx] = Restore.coast(coast_bits, cmin)[idx]
+
+ 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, cmin)[idx]
+
+ # bits translation MVCM-python -> MVCM-C
+ # 01 13 test_11um
+ # 02 27 surface_temperature_test
+ # 03 27 sst_test
+ # 04 24 bt_diff_86_11um
+ # 05 18 test_11_12um_diff
+ # 07 19 oceanic_stratus_11_4um_test
+ # 10 20 nir_reflectance_test
+ # 11 21 vis_nir_ratio
+ # 12 23 test_16_21um_reflectance
+ # 15 16 test_1_38um_high_clouds
+ # 16 17 thin_cirrus_4_12um_BTD_test
+
+ # if scene_name in ['Ocean_Day']:
+ # # total_bit = np.full(viirs_data.M01.shape, 4)
+ # total_bit = bit1 + bit2 + bit4 + bit11
+ # sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
+ # scene_flags = scn.find_scene(viirs_data, sunglint_angle)
+ #
+ # # 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(viirs_data, thresholds['Sun_Glint'], scene_flags, cmin)[idx]
+ #
+ # idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['sunglint'] == 1) &
+ # (scene_flags['uniform'] == 1) & (cmin <= 0.95))
+ # cmin[idx] = restoral.sunglint(viirs_data, thresholds['Sun_Glint'], total_bit, cmin)[idx]
+ #
+ # # MVCM = tst.ComputeTests(viirs_data, scene_name, thresholds)
+ # # cmin = MVCM.clear_sky_restoral(cmin)
+ #
+ # # return cmin
+ np.savez('test_w_restorals', confidence=cmin, lat=viirs_data.latitude.values, lon=viirs_data.longitude.values)
##########################################################
diff --git a/mvcm/read_data.py b/mvcm/read_data.py
index 2f1234a..4efd9c8 100644
--- a/mvcm/read_data.py
+++ b/mvcm/read_data.py
@@ -3,6 +3,7 @@ import numpy as np
import numpy.typing as npt
import ancillary_data as anc
+import mvcm.scene as scn
import os
import logging
@@ -454,7 +455,65 @@ class ReadAncillary(CollectInputs):
return ancillary_data
+def get_data(satellite: str,
+ sensor: str,
+ file_names: Dict) -> xr.Dataset:
+ """Store every variable in one xarray dataset
+
+ Parameters
+ ----------
+ satellite: str
+ satellite name
+ sensor: str
+ sensor name
+ file_names: Dict
+ dictionary containing all input file names
+
+ Returns
+ -------
+ 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)
+
+ geo_data = Viirs.read_viirs_geo()
+ viirs_data = Viirs.read_viirs_l1b(geo_data.solar_zenith.values)
+ viirs_data = Viirs.preprocess_viirs(geo_data, viirs_data)
+
+ Ancillary = ReadAncillary(latitude=geo_data.latitude.values,
+ longitude=geo_data.longitude.values,
+ resolution=1,
+ 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())
+
+ scene_flags = scn.find_scene(viirs_data, geo_data.sunglint_angle.values)
+ scene = scn.scene_id(scene_flags)
+
+ scene_xr = xr.Dataset()
+ for s in scn._scene_list:
+ 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'])
+
+ return input_data
+
+
"""
+
scene_flags = scn.find_scene(viirs, sunglint_angle)
scene = scn.scene_id(scene_flags)
diff --git a/mvcm/restoral.py b/mvcm/restoral.py
index 89998eb..8aeba69 100644
--- a/mvcm/restoral.py
+++ b/mvcm/restoral.py
@@ -97,16 +97,16 @@ class Restoral(object):
var_bt = spatial_var(bt, threshold['var_11um'])
var_reflectance = spatial_var(reflectance, threshold['var_0_86um'])
- idx = np.nonzero((self.confidence > 0.95) & (self.scene_flag['day'] == 1) &
+ idx = np.nonzero((confidence > 0.95) & (self.scene_flag['day'] == 1) &
(var_bt == 1) & (var_reflectance == 1))
confidence[idx] = 1
- idx = np.nonzero((self.confidence > 0.66) & (
+ idx = np.nonzero((confidence > 0.66) & (
((self.scene_flag['night'] == 1) & (var_bt == 1)) |
((self.scene_flag['day'] == 1) & (var_bt == 1) & (var_reflectance == 1))))
confidence[idx] = 0.96
- idx = np.nonzero((self.confidence <= 0.66) & (var_bt == 1))
+ idx = np.nonzero((confidence <= 0.66) & (var_bt == 1))
confidence[idx] = 0.67
logger.info('spatial_variability restoral test performed successfully')
@@ -170,18 +170,19 @@ class Restoral(object):
rad = self.data.M15.values
m05m04ratio = self.data.M08.values/self.data.M04.values
- m20m22diff = self.data.M12.values - self.data.M13.self.values
+ m20m22diff = self.data.M12.values - self.data.M13.values
m22m31diff = self.data.M13.values - self.data.M15.values
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 = thr11um - utils.bt11_elevation_adjustment(self.data.height.values).ravel()
- thr11um = thr11um.reshape(rad)
+ for i in range(3):
+ 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], 3), 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)
+ 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) &
@@ -213,7 +214,7 @@ class Restoral(object):
confidence: np.ndarray) -> np.ndarray:
ndvi = utils.compute_ndvi(self.data.M05.values, self.data.M07.values)
- threshold = self.threshold['Coastal_NDVI_and_Shallow_Water']
+ threshold = self.thresholds['Coastal_NDVI_Thresholds']
idx = np.nonzero((self.scene_flag['day'] == 1) & (self.scene_flag['snow'] == 0) &
(self.scene_flag['ice'] == 0) & (self.scene_flag['coast'] == 1) &
@@ -232,19 +233,21 @@ class Restoral(object):
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']))
- thr11um = threshold['bt11_land_night'] - utils.bt11_elevation_adjustment(self.data.height.values).ravel()
+ for i in range(3):
+ 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.ravel() > thr11um[0]))
+ (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.ravel() > thr11um[1]))
+ (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.ravel() > thr11um[2]))
+ (rad > thr11um[:, :, 2]))
confidence[idx] = 1
logger.info('land night restoral test performed successfully')
diff --git a/mvcm/scene.py b/mvcm/scene.py
index eaf6ee4..b1f12e3 100644
--- a/mvcm/scene.py
+++ b/mvcm/scene.py
@@ -92,11 +92,11 @@ def test_scene():
def find_scene(data, sunglint_angle):
- # 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_snowfr'].values
- icefr = data['geos_icefr'].values
+ 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
@@ -107,7 +107,7 @@ def find_scene(data, sunglint_angle):
b086 = data['M07'].values
elev = data['height'].values
ndvibk = data['ndvi'].values
- sfctmp = data['geos_sfct'].values
+ sfctmp = data['geos_surface_temperature'].values
dim1 = data.latitude.shape[0]
dim2 = data.latitude.shape[1]
@@ -240,9 +240,9 @@ def find_scene(data, sunglint_angle):
scene_flag['vrused'] = np.ones((dim1, dim2))
scene_flag['vrused'][idx] = 0
- snow_fraction = data['geos_snowfr']
- perm_ice_fraction = data['geos_landicefr']
- ice_fraction = data['geos_icefr']
+ 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
diff --git a/mvcm/utility_functions.py b/mvcm/utility_functions.py
new file mode 100644
index 0000000..69c920f
--- /dev/null
+++ b/mvcm/utility_functions.py
@@ -0,0 +1,47 @@
+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:
+ """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)
+ return local_mean
+
+
+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)
+ return local_std
+
+
+def spatial_var(rad: np.ndarray,
+ threshold: Dict) -> np.ndarray:
+ """Compute 3x3 spatial variability and returns number of pixels within certain threshold"""
+ 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.ones(test.shape)
+ var[test < threshold] = 0
+
+ return var.sum(axis=2)
+
+
+def compute_ndvi(b01: np.ndarray,
+ b02: np.ndarray) -> np.ndarray:
+ """Compute Normalized Difference Vegetation Index (NDVI)"""
+ 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
--
GitLab