# Copyright Iris contributors
#
# This file is part of Iris and is released under the BSD license.
# See LICENSE in the root of the repository for full licensing details.
"""Provides NIMROD file format capabilities.
.. z_reference:: iris.fileformats.nimrod
:tags: topic_load_save
API reference
"""
from enum import Enum
import glob
import os
import struct
import sys
import numpy as np
import iris
from iris.exceptions import TranslationError
import iris.fileformats.nimrod_load_rules
from iris.fileformats.nimrod_load_rules import Table
# general header (int16) elements 1-31 (Fortran bytes 1-62)
general_header_int16s = (
"vt_year",
"vt_month",
"vt_day",
"vt_hour",
"vt_minute",
"vt_second",
"dt_year",
"dt_month",
"dt_day",
"dt_hour",
"dt_minute",
"datum_type",
"datum_len",
"experiment_num",
"horizontal_grid_type",
"num_rows",
"num_cols",
"nimrod_version",
"field_code",
"vertical_coord_type",
"reference_vertical_coord_type",
"data_specific_float32_len",
"data_specific_int16_len",
"origin_corner",
"int_mdi",
"period_minutes",
"num_model_levels",
"proj_biaxial_ellipsoid",
"ensemble_member",
"model_origin_id",
"averagingtype",
)
# general header (float32) elements 32-59 (Fortran bytes 63-174)
general_header_float32s = (
"vertical_coord",
"reference_vertical_coord",
"y_origin",
"row_step",
"x_origin",
"column_step",
"float32_mdi",
"MKS_data_scaling",
"data_offset",
"x_offset",
"y_offset",
"true_origin_latitude",
"true_origin_longitude",
"true_origin_easting",
"true_origin_northing",
"tm_meridian_scaling",
"threshold_value_alt",
"threshold_value",
)
# data specific header (float32) elements 60-104 (Fortran bytes 175-354)
data_header_float32s = (
"tl_y",
"tl_x",
"tr_y",
"tr_x",
"br_y",
"br_x",
"bl_y",
"bl_x",
"sat_calib",
"sat_space_count",
"ducting_index",
"elevation_angle",
"neighbourhood_radius",
"threshold_vicinity_radius",
"recursive_filter_alpha",
"threshold_fuzziness",
"threshold_duration_fuzziness",
)
# data specific header (char) elements 105-107 (bytes 355-410)
# units, source and title
# data specific header (int16) elements 108-159 (Fortran bytes 411-512)
table_1_data_header_int16s = (
(
"radar_number",
"radar_sites",
"additional_radar_sites",
"clutter_map_number",
"calibration_type",
"bright_band_height",
"bright_band_intensity",
"bright_band_test_param_1",
"bright_band_test_param_2",
"infill_flag",
"stop_elevation",
"copy_vertical_coord",
"copy_reference_vertical_coord",
"copy_y_origin",
"copy_row_step",
"copy_x_origin",
"copy_column_step",
"copy_float32_mdi",
"copy_MKS_data_scaling",
"copy_data_offset",
"copy_x_offset",
"copy_y_offset",
"copy_true_origin_latitude",
"copy_true_origin_longitude",
"copy_tl_y",
"copy_tl_x",
"copy_tr_y",
"copy_tr_x",
"copy_br_y",
"copy_br_x",
"copy_bl_y",
"copy_bl_x",
"sensor_identifier",
"meteosat_identifier",
"availability_of_synop_meteosat",
"software_identifier",
"software_major_version",
"software_minor_version",
"software_micro_version",
)
+ tuple(f"data_header_int16_{i}" for i in range(48, 59))
+ ("period_seconds",)
)
# data specific header (int16) elements 108-159 (Fortran bytes 411-512)
table_2_data_header_int16s = (
(
"threshold_type",
"probability_method",
"recursive_filter_iterations",
"member_count",
"probability_period_of_event",
)
+ tuple(f"data_header_int16_{i}" for i in range(5, 59))
+ ("period_seconds",)
)
# data specific header (int16) elements 108-159 (Fortran bytes 411-512)
table_3_data_header_int16s = (
"data_header_int16_00",
"data_header_int16_01",
"data_header_int16_02",
"data_header_int16_03",
"data_header_int16_04",
"data_header_int16_05",
"soil_type",
"data_header_int16_07",
"data_header_int16_08",
"data_header_int16_09",
"data_header_int16_10",
) + tuple(f"data_header_int16_{i}" for i in range(11, 60))
# data specific header (int16) elements 108-159 (Fortran bytes 411-512)
table_4_data_header_int16s = (
"data_header_int16_00",
"data_header_int16_01",
"data_header_int16_02",
"data_header_int16_03",
"data_header_int16_04",
"data_header_int16_05",
"data_header_int16_06",
"radiation_code",
"data_header_int16_08",
"data_header_int16_09",
"data_header_int16_10",
) + tuple(f"data_header_int16_{i}" for i in range(11, 60))
def _read_chars(infile, num):
"""Read characters from the (big-endian) file."""
instr = infile.read(num)
result = struct.unpack(">%ds" % num, instr)[0]
result = result.decode()
return result
[docs]
class NimrodField:
"""A data field from a NIMROD file.
Capable of converting itself into a :class:`~iris.cube.Cube`
References
----------
Met Office (2003): Met Office Rain Radar Data from the NIMROD System.
NCAS British Atmospheric Data Centre, date of citation.
https://catalogue.ceda.ac.uk/uuid/82adec1f896af6169112d09cc1174499
"""
def __init__(self, from_file=None):
"""Create a NimrodField object and optionally read from an open file.
Example::
with open("nimrod_file", "rb") as infile:
field = NimrodField(infile)
"""
if from_file is not None:
self.read(from_file)
[docs]
def read(self, infile):
"""Read the next field from the given file object."""
self._read_header(infile)
self._read_data(infile)
def _read_header_subset(self, infile, names, dtype):
# Read contiguous header items of the same data type.
values = np.fromfile(infile, dtype=dtype, count=len(names))
if sys.byteorder == "little":
values.byteswap(True)
for i, name in enumerate(names):
setattr(self, name, values[i])
def _read_header(self, infile):
"""Load the 512 byte header (surrounded by 4-byte length)."""
leading_length = struct.unpack(">L", infile.read(4))[0]
if leading_length != 512:
raise TranslationError("Expected header leading_length of 512")
# general header (int16) elements 1-31 (bytes 1-62)
self._read_header_subset(infile, general_header_int16s, np.int16)
# general header (float32) elements 32-59 (bytes 63-174)
self._read_header_subset(infile, general_header_float32s, np.float32)
# skip unnamed floats
infile.seek(4 * (28 - len(general_header_float32s)), os.SEEK_CUR)
# data specific header (float32) elements 60-104 (bytes 175-354)
self._read_header_subset(infile, data_header_float32s, np.float32)
# skip unnamed floats
infile.seek(4 * (45 - len(data_header_float32s)), os.SEEK_CUR)
# data specific header (char) elements 105-107 (bytes 355-410)
self.units = _read_chars(infile, 8)
self.source = _read_chars(infile, 24)
self.title = _read_chars(infile, 24)
# determine which of Table 1, 2, 3 or 4 is being used
Table_3_field_codes = [
18,
144,
190,
191,
192,
193,
194,
196,
197,
198,
199,
201,
202,
203,
204,
218,
219,
301,
302,
901,
8229,
8230,
]
Table_4_field_codes = [90, 91, 92, 93, 96, 303]
default_float_threshold = self.float32_mdi
threshold_set = (
self.threshold_value != default_float_threshold
or self.threshold_value_alt != default_float_threshold
)
# The `Table` enum is defined in iris.fileformats.nimrod_load_rules
if self.field_code in Table_3_field_codes:
table = Table.table_3
data_header_int16s = table_3_data_header_int16s
elif self.field_code in Table_4_field_codes:
table = Table.table_4
data_header_int16s = table_4_data_header_int16s
elif threshold_set:
table = Table.table_2
data_header_int16s = table_2_data_header_int16s
else:
table = Table.table_1
data_header_int16s = table_1_data_header_int16s
self.table = table
# data specific header (int16) elements 108- (bytes 411-512)
self._read_header_subset(infile, data_header_int16s, np.int16)
# skip unnamed int16s
infile.seek(2 * (51 - len(data_header_int16s)), os.SEEK_CUR)
trailing_length = struct.unpack(">L", infile.read(4))[0]
if trailing_length != leading_length:
raise TranslationError(
"Expected header trailing_length of {}, got {}.".format(
leading_length, trailing_length
)
)
def _read_data(self, infile):
"""Read the data array: int8, int16, int32 or float32.
(surrounded by 4-byte length, at start and end)
"""
# what are we expecting?
num_data = int(self.num_rows) * int(self.num_cols)
num_data_bytes = int(num_data) * int(self.datum_len)
# format string for unpacking the file.read()
# 0:real
if self.datum_type == 0:
numpy_dtype = np.float32
# 1:int
elif self.datum_type == 1:
if self.datum_len == 1:
numpy_dtype = np.int8
elif self.datum_len == 2:
numpy_dtype = np.int16
elif self.datum_len == 4:
numpy_dtype = np.int32
else:
raise TranslationError("Undefined datum length %d" % self.datum_type)
# 2:byte
elif self.datum_type == 2:
numpy_dtype = np.byte
else:
raise TranslationError("Undefined data type")
leading_length = struct.unpack(">L", infile.read(4))[0]
if leading_length != num_data_bytes:
raise TranslationError(
"Expected data leading_length of %d" % num_data_bytes
)
self.data = np.fromfile(infile, dtype=numpy_dtype, count=num_data)
if sys.byteorder == "little":
self.data.byteswap(True)
trailing_length = struct.unpack(">L", infile.read(4))[0]
if trailing_length != leading_length:
raise TranslationError(
"Expected data trailing_length of %d" % num_data_bytes
)
# Form the correct shape.
self.data = self.data.reshape(self.num_rows, self.num_cols)
[docs]
def load_cubes(filenames, callback=None):
"""Load cubes from a list of NIMROD filenames.
Parameters
----------
filenames :
List of NIMROD filenames to load.
callback : optional
A function which can be passed on to :func:`iris.io.run_callback`.
Notes
-----
The resultant cubes may not be in the same order as in the files.
"""
if isinstance(filenames, str):
filenames = [filenames]
for filename in filenames:
for path in glob.glob(filename):
with open(path, "rb") as infile:
while True:
try:
field = NimrodField(infile)
except struct.error:
# End of file. Move on to the next file.
break
cube = iris.fileformats.nimrod_load_rules.run(field)
# Were we given a callback?
if callback is not None:
cube = iris.io.run_callback(callback, cube, field, filename)
if cube is None:
continue
yield cube
