Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Isccpg ng #82

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
Open

Isccpg ng #82

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
e2a047a
First try
May 25, 2022
167ab89
Apply sunz correction
Aug 16, 2022
04aebd9
Treat isccp as Meteosat-11
Aug 22, 2023
46326cb
Merge branch 'main' into isccpg_ng
Aug 22, 2023
e01802d
Merge branch 'main' into isccpg_ng
Jul 5, 2024
35ba9ef
Merge branch 'main' into isccpg_ng
Nov 25, 2024
08348cb
Update to read L1G data
Nov 27, 2024
4c0d51e
fix reader that broke in recent commits
ninahakansson Nov 27, 2024
918af9b
Add the reclibration for VIS channels.
ninahakansson Nov 28, 2024
10da412
Add homogenization
ninahakansson Nov 29, 2024
5a5a1c2
Checking in version one of correction.
ninahakansson Dec 6, 2024
737bd75
Got to version 2
ninahakansson Dec 6, 2024
852aa41
Merge branch 'main' into isccpg_ng
ninahakansson Dec 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 55 additions & 0 deletions bin/isccpng2pps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2019 level1c4pps developers
#
# This file is part of level1c4pps
#
# level1c4pps is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# level1c4pps is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with level1c4pps. If not, see <http://www.gnu.org/licenses/>.
# Author(s):

# Martin Raspaud <[email protected]>
# Nina Hakansson <[email protected]>
# Adam.Dybbroe <[email protected]>


"""Script to make ISCCP-NG level1c in PPS-format with pytroll."""

import argparse
from level1c4pps.isccpng2pps_lib import process_one_scene


# -----------------------------------------------------------------------------
# Main:
if __name__ == "__main__":
""" Create PPS-format level1c data
From a list of hirt files hrit create a level1c file for pps.
"""
parser = argparse.ArgumentParser(
description=('Script to produce a PPS-level1c file for a '
'ISCCP-NG scene.'))
parser.add_argument('files', metavar='fileN', type=str, nargs='+',
help='List of ISCCP-NG files to process')
parser.add_argument('-o', '--out_dir', type=str, nargs='?',
required=False, default='.',
help="Output directory where to store level1c file.")
parser.add_argument('-ne', '--nc_engine', type=str, nargs='?',
required=False, default='h5netcdf',
help="Engine for saving netcdf files netcdf4 or h5netcdf (default).")
parser.add_argument('-on', '--orbit_number', type=int, nargs='?',
required=False, default=99999,
help="Orbit number (default is 99999).")
options = parser.parse_args()
process_one_scene(options.files, options.out_dir,
engine=options.nc_engine,
orbit_n=options.orbit_number)
15 changes: 10 additions & 5 deletions level1c4pps/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -314,7 +314,7 @@ def get_band_encoding(dataset, bandnames, pps_tagnames, chunks=None):
# pygac qual flags
enc = {'dtype': 'int16', 'zlib': True,
'complevel': 4, '_FillValue': -32001.0}
elif name in ['scanline_timestamps']:
elif name in ['scanline_timestamps', 'pixel_time']:
# pygac scanline_timestamps
enc = {'dtype': 'int64',
'zlib': True,
Expand Down Expand Up @@ -435,7 +435,8 @@ def set_header_and_band_attrs_defaults(scene, BANDNAMES, PPS_TAGNAMES, REFL_BAND
# Assume factor applied if available as attribute.
scene[band].attrs['sun_earth_distance_correction_applied'] = 'True'
fix_sun_earth_distance_correction_factor(scene, band, irch.attrs['start_time'])
scene[band].attrs['wavelength'] = scene[band].attrs['wavelength'][0:3]
if 'wavelength' in scene[band].attrs:
scene[band].attrs['wavelength'] = scene[band].attrs['wavelength'][0:3]
scene[band].attrs['sun_zenith_angle_correction_applied'] = 'False'
if idtag in PPS_TAGNAMES_TO_IMAGE_NR:
scene[band].attrs['name'] = PPS_TAGNAMES_TO_IMAGE_NR[idtag]
Expand All @@ -452,7 +453,7 @@ def set_header_and_band_attrs_defaults(scene, BANDNAMES, PPS_TAGNAMES, REFL_BAND

# Add time coordinate. To make cfwriter aware that we want 3D data.
scene[band].coords['time'] = irch.attrs['start_time']

# Remove some attributes and coordinates
for attr in RENAME_VARS:
if attr in scene[band].attrs:
Expand All @@ -466,11 +467,15 @@ def set_header_and_band_attrs_defaults(scene, BANDNAMES, PPS_TAGNAMES, REFL_BAND
attr_value = scene[band].attrs.pop(attr, None)
if attr not in scene.attrs:
scene.attrs[attr] = attr_value
for coord_name in ['acq_time', 'latitude', 'longitude']:
for coord_name in ['acq_time', 'latitude', 'longitude', 'start_time', 'end_time', "crs"]:
try:
del scene[band].coords[coord_name]
except KeyError:
pass
try:
del scene[band].encoding['coordinates']
except (AttributeError, KeyError):
pass
return nimg


Expand Down Expand Up @@ -568,7 +573,7 @@ def compose_filename(scene, out_path, instrument, band=None):
end_time = band.attrs['end_time']
platform_name = scene.attrs['platform']
orbit_number = int(scene.attrs['orbit_number'])
out_path_with_dates = start_time.strftime(out_path)
out_path_with_dates = datetime.strftime(dt64_to_datetime(start_time), out_path)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe worth mentioning in the command line interface help that the out_path can be reformatted in this way.

filename = os.path.join(
out_path_with_dates,
"S_NWC_{:s}_{:s}_{:05d}_{:s}Z_{:s}Z.nc".format(
Expand Down
1 change: 0 additions & 1 deletion level1c4pps/eumgacfdr2pps_lib.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,6 @@ def set_header_and_band_attrs(scene, orbit_n=99999):
continue
if band in REFL_BANDS:
scene[band].attrs['sun_earth_distance_correction_applied'] = 'True'
del scene[band].encoding['coordinates']
return nimg


Expand Down
268 changes: 268 additions & 0 deletions level1c4pps/isccpng2pps_lib.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,268 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2019 level1c4pps developers
#
# This file is part of level1c4pps
#
# level1c4pps is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# level1c4pps is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with level1c4pps. If not, see <http://www.gnu.org/licenses/>.
# Author(s):

# Nina Hakansson <[email protected]>

"""Functions to convert ISCCP Next Generation level-1-G data to a NWCSAF/PPS level-1c formatet netCDF/CF file."""

import os
import time
import xarray as xr
import numpy as np
from satpy.scene import Scene
from level1c4pps import (get_encoding,
dt64_to_datetime,
compose_filename,
apply_sunz_correction,
rename_latitude_longitude,
update_angle_attributes, get_header_attrs,
set_header_and_band_attrs_defaults,
convert_angles,
adjust_lons_to_valid_range)

import logging


logger = logging.getLogger('isccpng2pps')

BANDNAMES = ['refl_01_60um',
'refl_00_65um',
'refl_00_86um',
'temp_03_80um',
'temp_08_60um',
'temp_09_70um',
'temp_11_00um',
'temp_12_00um',
'temp_13_30um',
'temp_06_20um',
'temp_07_30um']


REFL_BANDS = ['refl_01_60um',
'refl_00_65um',
'refl_00_86um']

ANGLE_NAMES = ['solar_zenith_angle', 'satellite_zenith_angle',
'solar_azimuth_angle', 'satellite_azimuth_angle']

PPS_TAGNAMES = {'refl_01_60um': "ch_r16",
'refl_00_65um': "ch_r06",
'refl_00_86um': "ch_r09",
'temp_03_80um': "ch_tb37",
'temp_08_60um': "ch_tb85",
'temp_09_70um': "ch_tbxx",
'temp_11_00um': "ch_tb11",
'temp_12_00um': "ch_tb12",
'temp_13_30um': "ch_tb133",
'temp_06_20um': "ch_tbxx",
'temp_07_30um': "ch_tb73"}

BANDNAMES = list(PPS_TAGNAMES.keys())

coef_slope_chan = ['refl_00_65um', 'refl_00_86um', 'refl_01_60um',
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe nicer to have put all coefficients related stuff in a yaml file with a format something like this:

coeffcients:
  channels: 
    - refl_00_65um 
    - refl_00_86um
    - refl_01_60um  
    - temp_03_80um
    - temp_06_20um
    - temp_07_30um
    - temp_08_60um  
    - temp_09_70um
    - temp_11_00um
    - temp_12_00um
    - temp_13_30um'
  satellite: 
    GOES-16: 
      day:
        slope: [1.00152, 0.946758, 1.00669, 0.926764, 0.991329, 0.999300, 1.01977, 1.03453, 0.996781, 1.02343, 0.863748]
        offset: [0.000369844, 0.193868, 0.708022, 18.4795, 1.58729, -0.0312020, -4.93224, -7.83331, 1.13395, -5.51940, 30.5574]
      night: 
        ...

'temp_03_80um', 'temp_06_20um', 'temp_07_30um', 'temp_08_60um',
'temp_09_70um', 'temp_11_00um', 'temp_12_00um', 'temp_13_30um']
coef_slope_list = {
'270_day': [1.00152, 0.946758, 1.00669, 0.926764, 0.991329, 0.999300, 1.01977, 1.03453, 0.996781, 1.02343, 0.863748],
'270_nig': [1.00000, 1.00000, 1.00000, 0.950274, 0.989769, 1.00713, 1.01844, 1.01421, 0.995356, 1.02101, 0.905866],
'271_day': [1.00298, 0.946479, 1.00762, 0.928358, 0.991436, 0.999576, 1.01939, 1.03306, 0.996656, 1.02405, 0.894260],
'271_nig': [1.00000, 1.00000, 1.00000, 0.950015, 0.989834, 1.00698, 1.01804, 1.01361, 0.995194, 1.02155, 0.930947],
'173_day': [1.00164, 0.944191, 1.00680, 0.929371, 1.00905, 1.00090, 1.01065, 1.02516, 0.996015, 1.02885, 0.871186],
'173_nig': [1.00000, 1.00000, 1.00000, 0.950006, 1.00397, 1.00545, 1.00990, 1.01022, 0.994367, 1.02563, 0.912099],
'55_day': [0.999586, 0.999657, 0.997100, 1.01192, 1.00193, 1.00093, 1.00202, 0.998663, 0.999412, 0.999663, 0.945356],
'55_nig': [1.00000, 1.00000, 1.00000, 1.00772, 1.00129, 1.00127, 1.00180, 0.999506, 0.999569, 0.999708, 0.961137],
}
coef_offset_list = {
'270_day': [0.000369844, 0.193868, 0.708022, 18.4795, 1.58729, -0.0312020, -4.93224, -7.83331, 1.13395, -5.51940, 30.5574],
'270_nig': [0.00000, 0.00000, 0.00000, 11.8884, 1.97083, -2.13396, -4.75410, -2.48929, 1.46463, -4.95287, 19.1043],
'271_day': [-0.00818343, 0.194117, 0.718540, 17.7929, 1.57229, -0.0877037, -4.83326, -7.49875, 1.17318, -5.66102, 23.6271],
'271_nig': [0.00000, 0.00000, 0.00000, 11.9196, 1.96629, -2.07708, -4.64315, -2.38190, 1.51330, -5.07539, 13.6827],
'173_day': [-4.45259e-05, 0.214025, 0.703233, 17.3952, -2.21718, -0.347357, -2.62520, -5.68860, 1.37321, -6.76384, 28.8614],
'173_nig': [0.00000, 0.00000, 0.00000, 11.9012, -0.963594, -1.56746, -2.52102, -1.75993, 1.76101, -5.99385, 17.7467],
'55_day': [0.00410180, 0.0160834, 0.101325, -2.83943, -0.417072, -0.181784, -0.486767, 0.304194, 0.134653, 0.0706959, 12.3763],
'55_nig': [0.00000, 0.00000, 0.00000, -1.61703, -0.258592, -0.278584, -0.442302, 0.0806762, 0.0877965, 0.0597462, 8.17879]
}

coef_slope_list_v1 = {
'270_day': [1.00166, 0.948021, 1.00000, 0.926764, 0.991329, 0.999300,
1.01977, 1.03453, 0.996781, 1.02343, 0.863748],
'270_all': [1.00000, 1.00000, 1.00000, 0.940848, 0.990431, 1.00483,
1.01947, 1.01835, 0.996272, 1.02239, 0.890045],
'271_day': [1.00313, 0.947739, 1.00000, 0.928358, 0.991436, 0.999576,
1.01939, 1.03306, 0.996656, 1.02405, 0.894260],
'271_all': [1.00000, 1.00000, 1.00000, 0.939285, 0.990504, 1.00481,
1.01907, 1.01759, 0.996129, 1.02298, 0.917607],
'173_day': [1.00169, 0.945509, 1.00000, 0.929371, 1.00905, 1.00090,
1.01065, 1.02516, 0.996015, 1.02885, 0.871186],
'173_all': [1.00000, 1.00000, 1.00000, 0.938706, 1.00526, 1.00413, 1.01050,
1.01338, 0.995399, 1.02742, 0.896835],
'55_day': [0.999577, 0.999643, 1.00000, 1.01192, 1.00193, 1.00093, 1.00202,
0.998663, 0.999412, 0.999663, 0.945356],
'55_all': [1.00000, 1.00000, 1.00000, 1.00882, 1.00146, 1.00121, 1.00194,
0.999321, 0.999513, 0.999684, 0.955400]
}
coef_offset_list_v1 = {
'270_day': [-0.00207069, 0.161403, 0.00000, 18.4795, 1.58729, -0.0312020,
-4.93224, -7.83331, 1.13395, -5.51940, 30.5574],
'270_all': [0.00000, 0.00000, 0.00000, 14.4615, 1.80946, -1.50406,
-4.94939, -3.61878, 1.23287, -5.28908, 23.4361],
'271_day': [-0.0107727, 0.161054, 0.00000, 17.7929, 1.57229, -0.0877037,
-4.83326, -7.49875, 1.17318, -5.66102, 23.6271],
'271_all': [0.00000, 0.00000, 0.00000, 14.7561, 1.80263, -1.48189,
-4.84116, -3.46801, 1.27653, -5.42165, 17.3323],
'173_day': [-0.00100716, 0.162379, 0.00000, 17.3952, -2.21718, -0.347357,
-2.62520, -5.68860, 1.37321, -6.76384, 28.8614],
'173_all': [0.00000, 0.00000, 0.00000, 14.8503, -1.29354, -1.20664,
-2.63734, -2.61762, 1.49858, -6.43452, 21.9239],
'55_day': [0.00444712, 0.0162967, 0.00000, -2.83943, -0.417072, -0.181784,
-0.486767, 0.304194, 0.134653, 0.0706959, 12.3763],
'55_all': [0.00000, 0.00000, 0.00000, -1.94343, -0.302762, -0.259928,
-0.471424, 0.131567, 0.104622, 0.0658078, 9.71430]
}
satellite_names = {270: "GOES-16", # ABI
271: "GOES-17", # ABI
173: "Himawari-8", # AHI
55: "Meteosat-8", # MSG1 SEVIRI
70: "Meteosat-11"}


def get_encoding_isccpng(scene):
"""Get netcdf encoding for all datasets."""
return get_encoding(scene,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should fit on one line.

BANDNAMES,
PPS_TAGNAMES,
chunks=None)


def set_header_and_band_attrs(scene, orbit_n=00000):
"""Set and delete some attributes."""
nimg = 0 # name of first dataset is image0
# Set some header attributes:
irch = scene['temp_11_00um']
irch.attrs['instrument'] = "seviri"
scene.attrs['source'] = "isccpng2pps.py"
scene.attrs['platform_name'] = "meteosat11"
nimg = set_header_and_band_attrs_defaults(scene, BANDNAMES, PPS_TAGNAMES, REFL_BANDS, irch, orbit_n=orbit_n)
for band in REFL_BANDS:
if band not in scene:
continue
scene[band].attrs['sun_zenith_angle_correction_applied'] = 'False'
return nimg


def homogenize_channel(scene, wmo_id, illum, band, sol_zen):
"""Homogenize one satellite for one band."""
data = scene[band]
dwmo_id = scene["wmo_id"].values
ch_index = coef_slope_chan.index(band)
k = coef_slope_list[str(wmo_id) + illum][ch_index]
c = coef_offset_list[str(wmo_id) + illum][ch_index]
logger.info(
f"Homogenizing channel {band} for {satellite_names[wmo_id]} for illum:{illum}"
f" with respect to SEVIRI on MSG4 using y={k} * x + {c}")
dn_discr = 90
if (illum == '_day'):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no need for parenthesis

update = ((dwmo_id == wmo_id) & (data > 0) & (sol_zen < dn_discr))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no need for the outermost parenthesis

scene[band].values = np.where(update, scene[band].values * k + c, scene[band].values)

if (illum == '_nig'):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same as above

update = ((dwmo_id == wmo_id) & (data > 0) & (sol_zen >= dn_discr))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same as above

scene[band].values = np.where(update, scene[band].values * k + c, scene[band].values)


def homogenize(scene):
"""Homogenize data to Meteosat-11."""
sol_zen = scene["sunzenith"]
for band in BANDNAMES:
for wmo_id in [270, 271, 173, 55]:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

for wmo_id in satellite_names:
     for illumination in ["_day", "_night"]:
         homogenize_channel(scene, wmo_id, illumination, band, sol_zen)

homogenize_channel(scene, wmo_id, '_day', band, sol_zen)
homogenize_channel(scene, wmo_id, '_nig', band, sol_zen)

def recalibrate_meteosat(scene):
"""Nominal calibration is applied, redo with meirnik calibration."""
from satpy.readers.seviri_base import MeirinkCalibrationHandler
start_time = dt64_to_datetime(scene["refl_00_65um"].attrs["start_time"])
channel_name = {"refl_00_65um": "VIS006",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I prefer to define constants like this and the ones below in the top of the module.

"refl_00_86um": "VIS008",
"refl_01_60um": "IR_016"}
platform_id = {55: 321,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do you actually need the 321 and 324 ids?

70: 324}
# Normally read from file (HRIT). We do not have access to nominal calibration
calibration_nominal = {2021: {321: {"VIS006": 24.974,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess you could use 55 and 70 here directly-

"VIS008": 32.377,
"IR_016": 23.710},
324: {"VIS006": 21.241,
"VIS008": 27.921,
"IR_016": 23.112}
}}
for wmo_id in platform_id:
for band in channel_name:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you could do: for wmo_id, platform_id in PLATFORM_ID.items(): and for band, channel_name in CHANNEL_NAME.items(): to increase readability if you follow the advise above.

dwmo_id = scene["wmo_id"]
old_gain = calibration_nominal[start_time.year][platform_id[wmo_id]][channel_name[band]]
meirink = MeirinkCalibrationHandler(calib_mode="MEIRINK-2023")
new_gain = 1000.0 * meirink.get_slope(platform_id[wmo_id], channel_name[band], start_time)
update = scene["wmo_id"].values == wmo_id
logger.info(
f"Recalibrating channel {band} for {satellite_names[wmo_id]} using y *= {new_gain} /{old_gain}")
scene[band].values = np.where(update, scene[band].values * new_gain / old_gain, scene[band].values)


def fix_pixel_time(scene):
"""Fix the time pixel variable, original file does not contain units."""
del scene["pixel_time"].coords["crs"]
scene["pixel_time"].encoding['coordinates'] = "lon lat"
scene["pixel_time"].data = scene["pixel_time"].data * np.timedelta64(1, 's') + scene['temp_11_00um'].attrs["start_time"]


def process_one_scene(scene_files, out_path,
engine='h5netcdf',
orbit_n=0):
"""Make level 1c files in PPS-format."""
tic = time.time()
scn_ = Scene(reader='multiple_sensors_isccpng_l1g_nc', filenames=scene_files)
scn_.load(BANDNAMES + ANGLE_NAMES + ["wmo_id", "pixel_time", "lon", "lat"])

# one ir channel
irch = scn_['temp_11_00um']

set_header_and_band_attrs(scn_, orbit_n=orbit_n)
fix_pixel_time(scn_)
# rename_latitude_longitude(scn_)
adjust_lons_to_valid_range(scn_)
convert_angles(scn_, delete_azimuth=True)
update_angle_attributes(scn_, irch)
recalibrate_meteosat(scn_)
homogenize(scn_)
apply_sunz_correction(scn_, REFL_BANDS)

filename = compose_filename(scn_, out_path, instrument='seviri', band=irch)
encoding = get_encoding_isccpng(scn_)

scn_.save_datasets(writer='cf',
filename=filename,
header_attrs=get_header_attrs(scn_, band=irch, sensor='seviri'),
engine=engine,
include_lonlats=False,
flatten_attrs=True,
encoding=encoding)
logging.info("Saved file {:s} after {:3.1f} seconds".format(
os.path.basename(filename),
time.time()-tic))
return filename
1 change: 1 addition & 0 deletions setup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@
url="https://github.com/foua-pps/level1c4pps",
packages=find_packages(),
scripts=['bin/seviri2pps.py',
'bin/isccpng2pps.py',
'bin/gac2pps.py',
'bin/mersi2pps.py',
'bin/mersi22pps.py',
Expand Down
Loading