Skip to content
Snippets Groups Projects
Commit d4650c0a authored by Paolo Veglio's avatar Paolo Veglio
Browse files

started testing the code. need to correct the way the final cmin before restoral is computed

parent eb1f7455
No related branches found
No related tags found
No related merge requests found
Pipeline #43453 failed
Hide whitespace changes
Inline Side-by-side
mvcm/main.py
+ 127
− 64
View file @ d4650c0a
...@@ -7,28 +7,28 @@ import numpy as np ...@@ -7,28 +7,28 @@ import numpy as np
from glob import glob from glob import glob
import read_data as rd import mvcm.read_data as rd
import spectral_tests as tst import mvcm.spectral_tests as tst
import scene as scn import mvcm.scene as scn
import restoral from mvcm.restoral import Restoral
# #################################################################### # # #################################################################### #
# TEST CASE # TEST CASE
# data: # data:
_datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires' _datapath = '/ships19/hercules/pveglio/mvcm_viirs_hires'
_fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0] _fname_mod02 = glob(f'{_datapath}/VNP02MOD.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0]
_fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1454.001.*.uwssec.nc')[0] _fname_mod03 = glob(f'{_datapath}/VNP03MOD.A2022173.1312.001.*.uwssec.nc')[0]
_fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1454.001.*.uwssec_bowtie_restored.nc')[0] _fname_img02 = glob(f'{_datapath}/VNP02IMG.A2022173.1312.001.*.uwssec_bowtie_corrected.nc')[0]
_fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1454.001.*.uwssec.nc')[0] _fname_img03 = glob(f'{_datapath}/VNP03IMG.A2022173.1312.001.*.uwssec.nc')[0]
# thresholds: # 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: # ancillary files:
_geos_atm_1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20220622_1200.V01.nc4' _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_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_land = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20220622_1330.V01.nc4'
_geos_ocean = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20220622_1430.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' _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' _ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.177.hdf'
_sst_file = 'oisst.20220622' _sst_file = 'oisst.20220622'
...@@ -51,7 +51,8 @@ def main(satellite: str = 'snpp', ...@@ -51,7 +51,8 @@ def main(satellite: str = 'snpp',
geos_ocean: str = _geos_ocean, geos_ocean: str = _geos_ocean,
geos_constants: str = _geos_constants, geos_constants: str = _geos_constants,
ndvi_file: str = _ndvi_file, 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. """Main function for the MVCM.
Parameters Parameters
...@@ -84,8 +85,10 @@ def main(satellite: str = 'snpp', ...@@ -84,8 +85,10 @@ def main(satellite: str = 'snpp',
file name of GEOS5 constants file name of GEOS5 constants
ndvi_file: str ndvi_file: str
NDVI file name NDVI file name
sst_file: sst_file: str
Sea surface temperature file name Sea surface temperature file name
eco_file: str
Ecosystem File
Returns Returns
------- -------
...@@ -103,6 +106,7 @@ def main(satellite: str = 'snpp', ...@@ -103,6 +106,7 @@ def main(satellite: str = 'snpp',
'GEOS_constants': f'{geos_constants}', 'GEOS_constants': f'{geos_constants}',
'NDVI': f'{ndvi_file}', 'NDVI': f'{ndvi_file}',
'SST': f'{sst_file}', 'SST': f'{sst_file}',
'ECO': f'{eco_file}',
'ANC_DIR': f'{data_path}/ancillary' 'ANC_DIR': f'{data_path}/ancillary'
} }
...@@ -112,67 +116,126 @@ def main(satellite: str = 'snpp', ...@@ -112,67 +116,126 @@ def main(satellite: str = 'snpp',
sunglint_angle = thresholds['Sun_Glint']['bounds'][3] 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 # Initialize the confidence arrays for the various test groups
cmin_G1 = np.ones(viirs_data.M01.shape) cmin_G1 = np.ones(viirs_data.M11.shape)
cmin_G2 = np.ones(viirs_data.M01.shape) cmin_G2 = np.ones(viirs_data.M11.shape)
cmin_G3 = np.ones(viirs_data.M01.shape) cmin_G3 = np.ones(viirs_data.M11.shape)
cmin_G4 = np.ones(viirs_data.M01.shape) cmin_G4 = np.ones(viirs_data.M11.shape)
cmin_G5 = np.ones(viirs_data.M01.shape) cmin_G5 = np.ones(viirs_data.M11.shape)
########################################################## ##########################################################
"""
for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert_Coast', for scene_name in ['Land_Day', 'Land_Day_Coast', 'Land_Day_Desert_Coast',
'Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean', 'Ocean_Day', 'Ocean_Night', 'Polar_Day_Ocean', 'Polar_Night_Ocean',
'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert', 'Polar_Day_Land', 'Polar_Day_Coast', 'Polar_Day_Desert',
'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Land_Night', 'Polar_Day_Desert_Coast', 'Polar_Day_Snow', 'Land_Night',
'Polar_Night_Land', 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']: 'Polar_Night_Land', 'Polar_Night_Snow', 'Day_Snow', 'Night_Snow']:
""" MyScene = tst.CloudTests(data=viirs_data,
for scene_name in ['Ocean_Day']: scene_name=scene_name,
MyScene = tst.CloudTests(viirs_data, scene_name, thresholds) thresholds=thresholds)
cmin_G1, bit1 = MyScene.test_11um('M15', cmin_G1, test_name='11um_Test') cmin_G1, bit01 = MyScene.test_11um('M15', cmin_G1)
cmin_G1, bit2 = MyScene.sst_test('M15', 'M16', cmin_G1, test_name='SST_Test') cmin_G1, bit02 = MyScene.surface_temperature_test('M15', viirs_data, cmin_G1)
cmin_G1, bit03 = MyScene.sst_test('M15', 'M16', cmin_G1)
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, bit04 = MyScene.bt_diff_86_11um('M14-M15', cmin_G2)
cmin_G2, bit5 = MyScene.oceanic_stratus_11_4um_test('M15-M12', cmin_G2, cmin_G2, bit05 = MyScene.test_11_12um_diff('M15-M16', cmin_G2)
test_name='11-4um_Oceanic_Stratus_Test') 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_G3, bit6 = MyScene.nir_reflectance_test('M07', cmin_G3, test_name='NIR_Reflectance_Test') # cmin_G2, bit08 = MyScene.bt_difference_11_4um_test_land('M15-M12', cmin_G2)
cmin_G3, bit7 = MyScene.vis_nir_ratio_test('M07-M05ratio', cmin_G3, test_name='Vis/NIR_Ratio_Test') cmin_G2, bit09 = MyScene.oceanic_stratus_11_4um_test('M15-M12', cmin_G2)
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.nir_reflectance_test('M07', cmin_G3)
cmin_G3, bit10 = MyScene.gemi_test('GEMI', cmin_G3, test_name='GEMI_Test') cmin_G3, bit11 = MyScene.vis_nir_ratio_test('M07-M05ratio', cmin_G3)
cmin_G3, bit12 = MyScene.test_16_21um_reflectance('M10', cmin_G3)
cmin_G4, bit11 = MyScene.test_1_38um_high_clouds('M09', cmin_G4, test_name='1.38um_High_Cloud_Test') # cmin_G3, bit13 = MyScene.visible_reflectance_test('M128', cmin_G3)
# cmin_G3, bit14 = MyScene.gemi_test('GEMI', cmin_G3)
cmin_G5, bit12 = MyScene.thin_cirrus_4_12um_BTD_test('M13-M16', cmin_G5,
test_name='4-12um_BTD_Thin_Cirrus_Test') cmin_G4, bit15 = MyScene.test_1_38um_high_clouds('M09', cmin_G4)
cmin = cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5 cmin_G5, bit16 = MyScene.thin_cirrus_4_12um_BTD_test('M13-M16', cmin_G5)
if scene_name in ['Ocean_Day']: # I need to modify this so that it's computed for each scene and that the exponent
# total_bit = np.full(viirs_data.M01.shape, 4) # is changed according to the number of groups in which tests are run
total_bit = bit1 + bit2 + bit4 + bit11 cmin = np.power(cmin_G1 * cmin_G2 * cmin_G3 * cmin_G4 * cmin_G5, 1/5)
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)
Restore = Restoral(data=viirs_data, thresholds=thresholds, scene_flag=scene_flags)
# idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['ice'] == 0) &
# (scene_flags['uniform'] == 1) & (cmin <= 0.99) & (cmin >= 0.05)) sunglint_bits = bit01 * bit03 * bit05 * bit15
# cmin[idx] = restoral.spatial(viirs_data, thresholds['Sun_Glint'], scene_flags, cmin)[idx] sh_water_bits = {'ir': bit01 * bit05,
'nir_1': bit10,
idx = np.nonzero((scene_flags['water'] == 1) & (scene_flags['sunglint'] == 1) & 'nir_2': bit15,
(scene_flags['uniform'] == 1) & (cmin <= 0.95)) 'sst': bit03
cmin[idx] = restoral.sunglint(viirs_data, thresholds['Sun_Glint'], total_bit, cmin)[idx] }
land_bits = bit15 * bit05
# MVCM = tst.ComputeTests(viirs_data, scene_name, thresholds) # make sure that the land_bits are calculated properly
# cmin = MVCM.clear_sky_restoral(cmin) idx = np.nonzero(scene_flags['desert'] == 0)
land_bits[idx] = land_bits[idx] * bit10[idx] * bit07[idx] * bit11[idx]
# return cmin coast_bits = bit05
np.savez('test_ams', confidence=cmin, lat=viirs_data.latitude.values, lon=viirs_data.longitude.values) 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)
########################################################## ##########################################################
......
mvcm/read_data.py
+ 59
− 0
View file @ d4650c0a
...@@ -3,6 +3,7 @@ import numpy as np ...@@ -3,6 +3,7 @@ import numpy as np
import numpy.typing as npt import numpy.typing as npt
import ancillary_data as anc import ancillary_data as anc
import mvcm.scene as scn
import os import os
import logging import logging
...@@ -454,7 +455,65 @@ class ReadAncillary(CollectInputs): ...@@ -454,7 +455,65 @@ class ReadAncillary(CollectInputs):
return ancillary_data 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_flags = scn.find_scene(viirs, sunglint_angle)
scene = scn.scene_id(scene_flags) scene = scn.scene_id(scene_flags)
......
mvcm/restoral.py
+ 17
− 14
View file @ d4650c0a
...@@ -97,16 +97,16 @@ class Restoral(object): ...@@ -97,16 +97,16 @@ class Restoral(object):
var_bt = spatial_var(bt, threshold['var_11um']) var_bt = spatial_var(bt, threshold['var_11um'])
var_reflectance = spatial_var(reflectance, threshold['var_0_86um']) 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)) (var_bt == 1) & (var_reflectance == 1))
confidence[idx] = 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['night'] == 1) & (var_bt == 1)) |
((self.scene_flag['day'] == 1) & (var_bt == 1) & (var_reflectance == 1)))) ((self.scene_flag['day'] == 1) & (var_bt == 1) & (var_reflectance == 1))))
confidence[idx] = 0.96 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 confidence[idx] = 0.67
logger.info('spatial_variability restoral test performed successfully') logger.info('spatial_variability restoral test performed successfully')
...@@ -170,18 +170,19 @@ class Restoral(object): ...@@ -170,18 +170,19 @@ class Restoral(object):
rad = self.data.M15.values 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.self.values m20m22diff = self.data.M12.values - self.data.M13.values
m22m31diff = self.data.M13.values - self.data.M15.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 = np.full((rad.ravel().shape[0], 3), np.array(threshold['bt11']))
thr11um[self.data.ecosystem.values.ravel() == 8, :] = np.array(threshold['bt11_desert']) thr11um[self.data.ecosystem.values.ravel() == 8, :] = np.array(threshold['bt11_desert'])
thr11um = thr11um - utils.bt11_elevation_adjustment(self.data.height.values).ravel() for i in range(3):
thr11um = thr11um.reshape(rad) 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 = 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[self.data.ecosystem.values.ravel() == 8] = np.array(threshold['m05m04_ratio_desert'])
thr_ratio = thr_ratio.reshape(rad) thr_ratio = thr_ratio.reshape(rad.shape)
idx = np.nonzero((bit == 1) & (self.scene_flag['day'] == 1) & (self.scene_flag['land'] == 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) & (self.scene_flag['snow'] == 0) & (self.scene_flag['ice'] == 0) &
...@@ -213,7 +214,7 @@ class Restoral(object): ...@@ -213,7 +214,7 @@ class Restoral(object):
confidence: np.ndarray) -> np.ndarray: confidence: np.ndarray) -> np.ndarray:
ndvi = utils.compute_ndvi(self.data.M05.values, self.data.M07.values) 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) & idx = np.nonzero((self.scene_flag['day'] == 1) & (self.scene_flag['snow'] == 0) &
(self.scene_flag['ice'] == 0) & (self.scene_flag['coast'] == 1) & (self.scene_flag['ice'] == 0) & (self.scene_flag['coast'] == 1) &
...@@ -232,19 +233,21 @@ class Restoral(object): ...@@ -232,19 +233,21 @@ class Restoral(object):
rad = self.data.M15.values rad = self.data.M15.values
threshold = self.thresholds['Land_Restoral'] threshold = self.thresholds['Land_Restoral']
thr11um = np.full((self.data.height.values.ravel().shape[0], 3), np.array(threshold['bt11'])) 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) & idx = np.nonzero((self.scene_flag['night'] == 1) & (self.scene_flag['snow'] == 0) &
(self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) & (self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) &
(rad.ravel() > thr11um[0])) (rad > thr11um[:, :, 0]))
confidence[idx] = 0.95 confidence[idx] = 0.95
idx = np.nonzero((self.scene_flag['night'] == 1) & (self.scene_flag['snow'] == 0) & idx = np.nonzero((self.scene_flag['night'] == 1) & (self.scene_flag['snow'] == 0) &
(self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) & (self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) &
(rad.ravel() > thr11um[1])) (rad > thr11um[:, :, 1]))
confidence[idx] = 0.96 confidence[idx] = 0.96
idx = np.nonzero((self.scene_flag['night'] == 1) & (self.scene_flag['snow'] == 0) & idx = np.nonzero((self.scene_flag['night'] == 1) & (self.scene_flag['snow'] == 0) &
(self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) & (self.scene_flag['ice'] == 0) & (bit == 1) & (confidence <= 0.95) &
(rad.ravel() > thr11um[2])) (rad > thr11um[:, :, 2]))
confidence[idx] = 1 confidence[idx] = 1
logger.info('land night restoral test performed successfully') logger.info('land night restoral test performed successfully')
......
mvcm/scene.py
+ 9
− 9
View file @ d4650c0a
...@@ -92,11 +92,11 @@ def test_scene(): ...@@ -92,11 +92,11 @@ def test_scene():
def find_scene(data, sunglint_angle): 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') # 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) # eco = np.array(data['eco'].values, dtype=np.uint8)
snowfr = data['geos_snowfr'].values snowfr = data['geos_snow_fraction'].values
icefr = data['geos_icefr'].values icefr = data['geos_ice_fraction'].values
lsf = np.array(data['land_water_mask'].values, dtype=np.uint8) lsf = np.array(data['land_water_mask'].values, dtype=np.uint8)
lat = data['latitude'].values lat = data['latitude'].values
lon = data['longitude'].values lon = data['longitude'].values
...@@ -107,7 +107,7 @@ def find_scene(data, sunglint_angle): ...@@ -107,7 +107,7 @@ def find_scene(data, sunglint_angle):
b086 = data['M07'].values b086 = data['M07'].values
elev = data['height'].values elev = data['height'].values
ndvibk = data['ndvi'].values ndvibk = data['ndvi'].values
sfctmp = data['geos_sfct'].values sfctmp = data['geos_surface_temperature'].values
dim1 = data.latitude.shape[0] dim1 = data.latitude.shape[0]
dim2 = data.latitude.shape[1] dim2 = data.latitude.shape[1]
...@@ -240,9 +240,9 @@ def find_scene(data, sunglint_angle): ...@@ -240,9 +240,9 @@ def find_scene(data, sunglint_angle):
scene_flag['vrused'] = np.ones((dim1, dim2)) scene_flag['vrused'] = np.ones((dim1, dim2))
scene_flag['vrused'][idx] = 0 scene_flag['vrused'][idx] = 0
snow_fraction = data['geos_snowfr'] snow_fraction = data['geos_snow_fraction']
perm_ice_fraction = data['geos_landicefr'] perm_ice_fraction = data['geos_land_ice_fraction']
ice_fraction = data['geos_icefr'] ice_fraction = data['geos_ice_fraction']
idx = tuple(np.nonzero((snow_fraction > 0.10) & (snow_fraction <= 1.0))) idx = tuple(np.nonzero((snow_fraction > 0.10) & (snow_fraction <= 1.0)))
scene_flag['map_snow'][idx] = 1 scene_flag['map_snow'][idx] = 1
......
mvcm/utility_functions.py 0 → 100644
+ 47
− 0
View file @ d4650c0a
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment