Skip to content

Commit

Permalink
Merge pull request #97 from salomoneliassonSMHI/SBAF_NN
Browse files Browse the repository at this point in the history 
Sbaf nn
  • Loading branch information
@ninahakansson
ninahakansson authored Nov 25, 2024
2 parents deed643 + 3db4362 commit 68482cd
Show file tree
Hide file tree
Showing 3 changed files with 210 additions and 15 deletions.
3 changes: 1 addition & 2 deletions bin/vgac2pps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,10 +46,9 @@
parser.add_argument('-all_ch', '--all_channels', action='store_true',
help="Save all 21 channels to level1c4pps file.")
parser.add_argument('-n19', '--as_noaa19',
choices=[version for version in SBAF],
default=None,
help=("Save only the AVHRR (1,2, 3B, 4, 5) channels to level1c4pps file. "
"And apply SBAFs to the channels."))
"And apply SBAFs to the channels. "))
parser.add_argument('-pps_ch', '--pps_channels', action='store_true',
help="Save only the necessary (for PPS) channels to level1c4pps file.")
parser.add_argument('-avhrr_ch', '--avhrr_channels', action='store_true',
Expand Down
4 changes: 4 additions & 0 deletions continuous_integration/environment.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,13 @@ name: test-environment
channels:
- conda-forge
dependencies:
- tensorflow
- scikit-learn
- sphinx
- scipy
- h5py
- python-geotiepoints
- matplotlib
- mock
- numpy<2.0.0
- satpy>0.41.1
Expand All @@ -18,3 +21,4 @@ dependencies:
- pip:
- trollsift
- pygac
- git+https://github.com/foua-pps/sbafs_ann@main
218 changes: 205 additions & 13 deletions level1c4pps/vgac2pps_lib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
import os
import time
from satpy.scene import Scene
from sbafs_ann.convert_vgac import convert_to_vgac_with_nn
from level1c4pps import (get_encoding, compose_filename,
set_header_and_band_attrs_defaults,
rename_latitude_longitude,
Expand All @@ -38,8 +39,6 @@
import logging
import numpy as np

# Example:

logger = logging.getLogger("vgac2pps")

# Order of BANDNAMES decides order of channels in file. Not important
Expand All @@ -66,7 +65,8 @@

MBAND_PPS = ["M05", "M07", "M09", "M10", "M11", "M12", "M14", "M15", "M16"]

MBAND_AVHRR = ["M05", "M07", "M12", "M15", "M16"]
# "M10", "M14" are not AVHRR channels, but needed for NN SABAF
MBAND_AVHRR = ["M05", "M07", "M12", "M15", "M16", "M10", "M14"]

MBAND_DEFAULT = ["M05", "M07", "M09", "M10", "M11", "M12", "M14", "M15", "M16"]

Expand Down Expand Up @@ -323,17 +323,143 @@
"slope": 1.002,
"offset": -0.69,
"comment": "Based on collocation data for VZA < 3 and SZA 0-180",
},
},
"KNMI": {
"r06": {
"viirs_channel": "M05",
"slope": 0.9401,
"offset": 0.629,
"comment": "VZA<60, SZA<60, Delta(VZA) < 5, Delta(Scat-angle) < 5",
},
"r09": {
"viirs_channel": "M07",
"slope": 0.9345,
"offset": -1.304,
"comment": "VZA<60, SZA<60, Delta(VZA) < 5, Delta(Scat-angle) < 5",
},
"tb37_night": {
"viirs_channel": "M12",
"slope": 0.9934,
"offset": 1.659,
"min_sunzenith": 89,
"comment": " VZA<60, SZA>95, Delta(VZA) < 5",
},
"tb37_day": {
"viirs_channel": "M12",
"slope": 0.9572,
"offset": 9.572,
"max_sunzenith": 80,
"comment": "VZA<60, SZA<60, Delta(VZA) < 5, Delta(Scat-angle) < 5",
},
"tb37_twilight": {
"viirs_channel": "M12",
"slope": 0.9753,
"offset": 5.6155,
"min_sunzenith": 80,
"max_sunzenith": 89,
"comment": "The linear average of the SBAFs for t37_day and t37_night. 80<= SZA <89",
},
"tb11": {
"viirs_channel": "M15",
"slope": 0.9986,
"offset": 0.600,
"comment": "",
},
"tb12": {
"viirs_channel": "M16",
"slope": 0.9943,
"offset": 1.328,
"comment": "",
},
},
"KNMI_v2": {
"r06": {
"viirs_channel": "M05",
"slope": 0.9401,
"offset": 0.629,
"comment": "VZA<60, SZA<60, Delta(VZA) < 5, Delta(Scat-angle) < 5",
},
"r09": {
"viirs_channel": "M07",
"slope": 0.9345,
"offset": -1.304,
"comment": "VZA<60, SZA<60, Delta(VZA) < 5, Delta(Scat-angle) < 5",
},
"tb37_day": {
"viirs_channel": "M12",
"slope": 0.9572,
"offset": 9.572,
"max_sunzenith": 70,
"comment": "VZA<60, SZA<60, Delta(VZA) < 5, Delta(Scat-angle) < 5",
},
"tb37_night": {
"viirs_channel": "M12",
"slope": 0.9934,
"offset": 1.659,
"min_sunzenith": 85,
"comment": "VZA<60, SZA>95, Delta(VZA) < 5",
},
"tb37_twilight": {
"viirs_channel": "M12",
"slope": None,
"offset": None,
"min_sunzenith": 70,
"max_sunzenith": 85,
"comment": "70 < SZA < 85: BT = (1-f)*BT(day) + f*BT(night), f=(SZA-70)/15",
},
"tb11": {
"viirs_channel": "M15",
"slope": 0.9986,
"offset": 0.600,
"comment": "VZA<60, SZA>95, Delta(VZA) < 5",
},
"tb12": {
"viirs_channel": "M16",
"slope": None,
"offset": None,
"comment": "VZA<60, SZA>95, Delta(VZA) < 5. Applies SBAF(tb11)-1.1646*(tb11-tb12)-0.235",
},
},
"NN_v1": {
"cfg_file_day": "ch7_SATZ_less_15_SUNZ_0_89_TD_1_min.yaml",
"cfg_file_night": "ch4_SATZ_less_25_SUNZ_90_180_TD_5_min.yaml",
"cfg_file_twilight": None,
"comment": "NN based on AVHRR and VGAC matchups using all AVHRR heritage channels"
},
"NN_v2": {
"cfg_file_day": "ch7_satz_max_25_SUNZ_0_80_tdiff_300_sec_20241031.yaml",
"cfg_file_night": "ch4_satz_max_25_SUNZ_90_180_tdiff_300_sec_20241031.yaml",
"cfg_file_twilight": "ch7_satz_max_25_SUNZ_80_89_tdiff_300_sec_20241031.yaml",
"comment": "NN based on AVHRR and VGAC matchups using all AVHRR heritage channels"
},
"NN_v3": {
"cfg_file_day": "ch7_satz_max_15_SUNZ_0_80_tdiff_120_sec_20241120.yaml",
"cfg_file_night": "ch4_satz_max_15_SUNZ_90_180_tdiff_120_sec_20241120.yaml",
"cfg_file_twilight": "ch7_satz_max_15_SUNZ_80_89_tdiff_120_sec_20241120.yaml",
"comment": "NN based on AVHRR and VGAC matchups using all AVHRR heritage channels"
}
}
}


def convert_to_noaa19(scene, sbaf_version):
""" Applies AVHRR SBAF to VGAC channels"""
def convert_to_noaa19_neural_network(scene, sbaf_version):
"""Applies AVHRR SBAF to VGAC channels using NN approach"""

if sbaf_version in ["NN_v1", "NN_v2", "NN_v3"]:
day_cfg_file = SBAF[sbaf_version]['cfg_file_day']
night_cfg_file = SBAF[sbaf_version]['cfg_file_night']
twilight_cfg_file = SBAF[sbaf_version]['cfg_file_twilight']
else:
logger.exception(f"Unrecognized NN version, {sbaf_version}")
scene = convert_to_vgac_with_nn(scene, day_cfg_file, night_cfg_file, twilight_cfg_file)

logger.info(f'Created NN version {sbaf_version}')

logger.info(f"Using SBAF_{sbaf_version}")

for avhhr_chan, scaling in SBAF[sbaf_version].items():
def convert_to_noaa19_linear(scene, SBAF):
""" Apply linear regression"""

for avhhr_chan, scaling in SBAF.items():
viirs_channel = scaling["viirs_channel"]
offset = scaling["offset"]
comment = scaling["comment"]
Expand All @@ -343,17 +469,83 @@ def convert_to_noaa19(scene, sbaf_version):
filt = filt & (scene["sunzenith"].values >= scaling["min_sunzenith"])
if "max_sunzenith" in scaling:
filt = filt & (scene["sunzenith"].values < scaling["max_sunzenith"])
scene[viirs_channel].values = np.where(
filt,
slope * scene[viirs_channel].values + offset,
scene[viirs_channel].values
)
scene[viirs_channel].values = np.where(filt,
slope * scene[viirs_channel].values + offset,
scene[viirs_channel].values
)
logger.info(f"{avhhr_chan:<13} = {slope:<6}*{viirs_channel:<3}+{offset:<5} ({comment})")


def convert_to_noaa19_KNMI_v2(scene, sbaf_version):
""" Apply 1 channel linear regression SBAF for KNMI version 2"""

# I need to save the t11 values before the SBAF adjustment as they are needed for the tb12 SBAF
tb11_original = scene["M15"].values.copy()
for avhrr_chan, scaling in SBAF[sbaf_version].items():
viirs_channel = scaling["viirs_channel"]
offset = scaling["offset"]
comment = scaling["comment"]
slope = scaling["slope"]
filt = np.ones_like(scene[viirs_channel].values, dtype=bool)
if "min_sunzenith" in scaling:
filt = filt & (scene["sunzenith"].values >= scaling["min_sunzenith"])
if "max_sunzenith" in scaling:
filt = filt & (scene["sunzenith"].values < scaling["max_sunzenith"])

if slope is not None:
# Then simple linear regression
scene[viirs_channel].values = np.where(
filt,
slope * scene[viirs_channel].values + offset,
scene[viirs_channel].values
)
logger.info(f"{avhrr_chan:<13} = {slope:<6}*{viirs_channel:<3}+{offset:<5} ({comment})")
else:
if avhrr_chan == "tb37_twilight":
# 70 < SZA < 85: BT = (1-f)*BT(day) + f*BT(night), f=(SZA-70)/15

f = (scene["sunzenith"].values - scaling["min_sunzenith"])/15
tb37_day_slope = scaling["tb37_day"]["slope"]
tb37_day_offset = scaling["tb37_day"]["offset"]
tb37_night_slope = scaling["tb37_night"]["slope"]
tb37_night_offset = scaling["tb37_night"]["offset"]
tb37_day = tb37_day_slope * scene[viirs_channel].values + tb37_day_offset
tb37_night = tb37_night_slope * scene[viirs_channel].values + tb37_night_offset

scene[viirs_channel].values = np.where(
filt,
(1-f)*tb37_day + f*tb37_night,
scene[viirs_channel].values
)

elif avhrr_chan == "tb12":
# BT(5) = BT(ch4)-1.1646*(BT(M15)-BT(M16))-0.235
scene[viirs_channel].values = scene["M15"].values-1.1646*(tb11_original-scene["M16"].values)-0.235
else:
logger.exception(f'Unknown channel, {avhrr_chan}, or missing slope parameter')
logger.info(f"{avhrr_chan:<13}: ({comment})")


def convert_to_noaa19(scene, sbaf_version):
""" Applies AVHRR SBAF to VGAC channels"""

logger.info(f"Using SBAF_{sbaf_version}")

if "NN" in sbaf_version:
convert_to_noaa19_neural_network(scene, sbaf_version)
elif sbaf_version == "KNMI_v2":
convert_to_noaa19_KNMI_v2(scene, sbaf_version)
else:
convert_to_noaa19_linear(scene, SBAF[sbaf_version])

if "npp" in scene.attrs["platform"].lower():
scene.attrs["platform"] = "vgacsnpp"
scene.attrs["platform"] = scene.attrs["platform"].replace("noaa", "vgac")

# These are only needed for the NN, but they must be removed before saving
del scene["M10"]
del scene["M14"]


def get_encoding_viirs(scene):
"""Get netcdf encoding for all datasets."""
Expand Down

0 comments on commit 68482cd

Please sign in to comment.