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

Refactor telemetry to use xarray #134

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

Refactor telemetry to use xarray #134

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8861320
Refactor telemetry
mraspaud Oct 18, 2024
91c42ff
Style refactor
mraspaud Oct 18, 2024
a421bae
Fix tests
mraspaud Oct 18, 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
10 changes: 5 additions & 5 deletions pygac/calibration/noaa.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,9 +79,9 @@ def calibrate(ds, custom_coeffs=None, coeffs_file=None):
corr
)

prt = ds["prt_counts"].data
ict = ds["ict_counts"].data
space = ds["space_counts"].data
prt = ds["PRT"].mean(dim="PRT_measurement").data
ict = ds["ICT"].mean(dim="back_scan").data[:, -3:]
space = ds["space_counts"].mean(dim="back_scan").data[:, -3:]
Comment on lines +83 to +84
Copy link
Member

Choose a reason for hiding this comment

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

Could

.data[:, -3:]

be replaced with

.sel(channel_name=slice(-3, None)).data

?


ir_channels_to_calibrate = [3, 4, 5]

Expand Down Expand Up @@ -525,8 +525,8 @@ def calibrate_thermal(counts, prt, ict, space, line_numbers, channel, cal):
nonzeros = np.logical_not(zeros)

space[zeros] = np.interp((zeros).nonzero()[0],
(nonzeros).nonzero()[0],
space[nonzeros])
(nonzeros).nonzero()[0],
space[nonzeros])

# convolving and smoothing PRT, ICT and SPACE values
if lines > 51:
Expand Down
4 changes: 2 additions & 2 deletions pygac/gac_klm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,8 +138,8 @@
("PRT", ">u2", (3, )),
("ch3_patch_temp", ">u2"),
("spare", ">u2"), ]),
("back_scan", ">u2", (30, )),
("space_data", ">u2", (50, )),
("back_scan", ">u2", (10, 3)),
("space_data", ">u2", (10, 5)),
("sync_delta", ">u2"),
("zero_fill5", ">i2"),
# AVHRR SENSOR DATA
Expand Down
34 changes: 14 additions & 20 deletions pygac/klm_reader.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,14 +36,10 @@

import datetime
import logging

try:
from enum import IntFlag
except ImportError:
# python version < 3.6, use a simple object without nice representation
IntFlag = object
from enum import IntFlag

import numpy as np
import xarray as xr

from pygac.correct_tsm_issue import TSM_AFFECTED_INTERVALS_KLM, get_tsm_idx
from pygac.reader import Reader, ReaderError
Expand Down Expand Up @@ -662,9 +658,9 @@ def read(self, filename, fileobj=None):
self.filename = filename
LOG.info("Reading %s", self.filename)
with file_opener(fileobj or filename) as fd_:
self.ars_head, self.head = self.read_header(
self.archive_header, self.head = self.read_header(
filename, fileobj=fd_)
if self.ars_head:
if self.archive_header:
ars_offset = ars_header.itemsize
else:
ars_offset = 0
Expand Down Expand Up @@ -742,23 +738,21 @@ def get_telemetry(self):
space_counts: np.array

"""
prt_counts = np.mean(self.scans["telemetry"]["PRT"], axis=1)
hmf_array = self.scans
prt_counts = xr.DataArray(hmf_array["telemetry"]["PRT"],
dims=["scan_line_index", "PRT_measurement"])

# getting ICT counts

ict_counts = np.zeros((len(self.scans), 3))
ict_counts[:, 0] = np.mean(self.scans["back_scan"][:, 0::3], axis=1)
ict_counts[:, 1] = np.mean(self.scans["back_scan"][:, 1::3], axis=1)
ict_counts[:, 2] = np.mean(self.scans["back_scan"][:, 2::3], axis=1)
ict_counts = xr.DataArray(hmf_array["back_scan"],
dims=["scan_line_index", "back_scan", "channel_name"],
coords=dict(channel_name=["3b", "4", "5"]))

# getting space counts
space_counts = xr.DataArray(self.split_array_along_channel_3(hmf_array["space_data"]),
dims=["scan_line_index", "back_scan", "channel_name"],
coords=dict(channel_name=["1", "2", "3a", "3b", "4", "5"]))

space_counts = np.zeros((len(self.scans), 3))
space_counts[:, 0] = np.mean(self.scans["space_data"][:, 2::5], axis=1)
space_counts[:, 1] = np.mean(self.scans["space_data"][:, 3::5], axis=1)
space_counts[:, 2] = np.mean(self.scans["space_data"][:, 4::5], axis=1)

return prt_counts, ict_counts, space_counts
return xr.Dataset(dict(PRT=prt_counts, ICT=ict_counts, space_counts=space_counts))

def _get_lonlat_from_file(self):
"""Get the longitudes and latitudes."""
Expand Down
4 changes: 2 additions & 2 deletions pygac/lac_klm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,8 +137,8 @@
("PRT", ">u2", (3, )),
("ch3_patch_temp", ">u2"),
("spare", ">u2"), ]),
("back_scan", ">u2", (30, )),
("space_data", ">u2", (50, )),
("back_scan", ">u2", (10, 3)),
("space_data", ">u2", (10, 5)),
("sync_delta", ">u2"),
("zero_fill4", ">i2"),
# EARTH OBSERVATIONS
Expand Down
69 changes: 46 additions & 23 deletions pygac/pod_reader.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
IntFlag = object

import numpy as np
import xarray as xr
from pyorbital.geoloc import compute_pixels, get_lonlatalt
from pyorbital.geoloc_instrument_definitions import avhrr_gac

Expand Down Expand Up @@ -278,9 +279,9 @@ def read(self, filename, fileobj=None):
LOG.info("Reading %s", self.filename)
# choose the right header depending on the date
with file_opener(fileobj or filename) as fd_:
self.tbm_head, self.head = self.read_header(
self.archive_header, self.head = self.read_header(
filename, fileobj=fd_, header_date=self.header_date)
if self.tbm_head:
if self.archive_header:
tbm_offset = tbm_header.itemsize
else:
tbm_offset = 0
Expand Down Expand Up @@ -347,8 +348,10 @@ def _validate_tbm_header(cls, potential_tbm_header):
return potential_tbm_header.copy()

# This will raise a DecodingError if the data_set_name is not valid.
cls._decode_data_set_name(data_set_name)
return potential_tbm_header.copy()
data_set_name = cls._decode_data_set_name(data_set_name)
tbm_header = potential_tbm_header.copy()
tbm_header["data_set_name"] = data_set_name.decode().strip("\x80").encode()
Copy link
Member

Choose a reason for hiding this comment

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

Could the .decode().strip().encode() part be moved to _decode_data_set_name?

return tbm_header


@classmethod
Expand Down Expand Up @@ -558,29 +561,49 @@ def get_telemetry(self):
space_counts: np.array

"""
number_of_scans = self.scans["telemetry"].shape[0]
decode_tele = np.zeros((int(number_of_scans), 105))
decode_tele[:, ::3] = (self.scans["telemetry"] >> 20) & 1023
decode_tele[:, 1::3] = (self.scans["telemetry"] >> 10) & 1023
decode_tele[:, 2::3] = self.scans["telemetry"] & 1023

prt_counts = np.mean(decode_tele[:, 17:20], axis=1)
telemetry = self.scans["telemetry"]
decode_tele = self.unpack_telemetry(telemetry)

hmf_dtype = np.dtype([("frame_sync", "u2", (6, )),
("id", [("id", "u2"),
("spare", "u2")]),
("timecode", "u2", (4, )),
("telemetry", [("ramp_calibration", "u2", (5, )),
("PRT", "u2", (3, )),
("ch3_patch_temp", "u2"),
("spare", "u2"), ]),
("back_scan", "u2", (10, 3)),
("space_data", "u2", (10, 5)),
("sync", "u2")])

hmf_array = decode_tele[:, :103].view(dtype=hmf_dtype).squeeze()
prt_counts = xr.DataArray(hmf_array["telemetry"]["PRT"],
dims=["scan_line_index", "PRT_measurement"])

# getting ICT counts
Copy link
Member

Choose a reason for hiding this comment

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

I think you could also remove the comments


ict_counts = np.zeros((int(number_of_scans), 3))
ict_counts[:, 0] = np.mean(decode_tele[:, 22:50:3], axis=1)
ict_counts[:, 1] = np.mean(decode_tele[:, 23:51:3], axis=1)
ict_counts[:, 2] = np.mean(decode_tele[:, 24:52:3], axis=1)
ict_counts = xr.DataArray(hmf_array["back_scan"],
dims=["scan_line_index", "back_scan", "channel_name"],
coords=dict(channel_name=["3", "4", "5"]))

# getting space counts

space_counts = np.zeros((int(number_of_scans), 3))
space_counts[:, 0] = np.mean(decode_tele[:, 54:100:5], axis=1)
space_counts[:, 1] = np.mean(decode_tele[:, 55:101:5], axis=1)
space_counts[:, 2] = np.mean(decode_tele[:, 56:102:5], axis=1)

return prt_counts, ict_counts, space_counts
space_counts = xr.DataArray(hmf_array["space_data"],
dims=["scan_line_index", "back_scan", "channel_name"],
coords=dict(channel_name=["1", "2", "3", "4", "5"]))

return xr.Dataset(dict(PRT=prt_counts, ICT=ict_counts, space_counts=space_counts))

def unpack_telemetry(self, telemetry):
"""Unpack the telemetry from 10 to 16 bits."""
# Documentation (POD guide, section 3) says:
# "The telemetry data are stored in 140 bytes. The first 103 (10 bit) words are packed three (10 bit) words in
# four bytes, right justified. The last four byte group contains one (10 bit) word with 20 trailing bits. All
# unused bits are set to zero."
number_of_scans = telemetry.shape[0]
decode_tele = np.zeros((int(number_of_scans), 105), dtype=np.uint16)
decode_tele[:, ::3] = (telemetry >> 20) & 1023
decode_tele[:, 1::3] = (telemetry >> 10) & 1023
decode_tele[:, 2::3] = telemetry & 1023
return decode_tele

@staticmethod
def _get_ir_channels_to_calibrate():
Expand Down
Loading
Loading