Add historic IASI L2 format (#10)

* - Update python example
- Update MetopVariable to work with any AbstractArray
- Add IASI L2 V10 format
- Expose IASI L2 GIARD variables

* - Add CHANGELOG
- update IASI information
- update documentation page

* Update python docs.

Author: lupemba

CHANGELOG.md

@@ -0,0 +1,7 @@
+# Changelog
+
+## Unreleased
+- No items yet.
+
+## v0.1
+- (**BREAKING**) Initial release with support for IASI and ASCAT formats.

Project.toml

@@ -1,7 +1,7 @@
 name = "MetopDatasets"
 uuid = "0c26954c-4046-4b98-a13c-f9377ca4b9b7"
 authors = ["lupemba <simon.koklupemba@eumetsat.int> and contributors"]
-version = "0.0.7"
+version = "0.1.0"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

docs/src/IASI.md

@@ -246,17 +246,8 @@ longitude, latitude = let
 end
 
 # Read the pressure levels 
-temp_pressure_levels, hum_pressure_levels = let 
-    giard = MetopDatasets.read_first_record(ds, MetopDatasets.GIADR_IASI_SND_02_V11)
-    
-    scale_factor_temp = MetopDatasets.get_scale_factor(MetopDatasets.GIADR_IASI_SND_02_V11, :pressure_levels_temp)
-    temp_level = giard.pressure_levels_temp/10^scale_factor_temp
-
-    scale_factor_humidity = MetopDatasets.get_scale_factor(MetopDatasets.GIADR_IASI_SND_02_V11, :pressure_levels_humidity)
-    humidity_level = giard.pressure_levels_humidity/10^scale_factor_humidity
-
-    temp_level, humidity_level
-end
+temp_pressure_levels = cfvariable(ds, "pressure_levels_temp", maskingvalue = NaN)[:]
+hum_pressure_levels  = cfvariable(ds, "pressure_levels_humidity", maskingvalue = NaN)[:]
 
 # Plot figure
 fig = let

docs/src/index.md

@@ -16,11 +16,11 @@ It is also possible to use MetopDatasets.jl from Python. See [Use with Python](@
 This section gives a very short overview of the core functionalities. The [MetopDataset](@ref) page is recommend for more information and more specific examples are given in the Example section e.g. [ASCAT](@ref). The [NCDatasets documentation](https://alexander-barth.github.io/NCDatasets.jl/stable/) is also a great resource for information on how to use the datasets. 
 
 ### Installation
-MetopDatasets.jl can be installed via Pkg and the url to the GitHub repository.
+MetopDatasets.jl can be installed via Pkg.
 
 ```julia
 import Pkg
-Pkg.add(url="https://github.com/eumetsat/MetopDatasets.jl#main") 
+Pkg.add("MetopDatasets")
 ```
 
 ### Read data from a Metop Native binary file
@@ -361,6 +361,4 @@ end
 ## Development status and versioning
 The package was previously named **MetopNative.jl** and was hosted on the [EUMETSAT GitLab](https://gitlab.eumetsat.int/eumetlab/cross-cutting-tools/MetopNative.jl).
 
-This package is still in the early development phase and is not yet registered. The plan is to register the package once it is complete for ASCAT and IASI products. The package will not follow any specific version system before it is registered.
-
-After registration, the aim will be to follow the [semantic versioning](https://semver.org/) system. Note that the package will start as 0.x.y to signal that breaking changes are to be expected. This is done to allow for more rapid development and because major dependencies like the [CommonDataModel.jl](https://github.com/JuliaGeo/CommonDataModel.jl) are not at the version 1.0 milestone yet. It is therefore recommended to use [Pkg environments](https://pkgdocs.julialang.org/v1/compatibility/) for projects with MetopDatasets.jl to handle comparability and ensure reproducibility. Please note that many of the "breaking changes" will be small and only affect specific use cases. This could for example be the correction of a single variable name. Any breaking change will be listed in the [release description](https://github.com/eumetsat/MetopDatasets.jl/releases).
+The aim is to follow the [semantic versioning](https://semver.org/) system. Note that the package is still 0.x.y to signal that breaking changes are to be expected. This is done to allow for more rapid development and because major dependencies like the [CommonDataModel.jl](https://github.com/JuliaGeo/CommonDataModel.jl) are not at the version 1.0 milestone yet. It is therefore recommended to use [Pkg environments](https://pkgdocs.julialang.org/v1/compatibility/) for projects with MetopDatasets.jl to handle comparability and ensure reproducibility. Please note that many of the "breaking changes" will be small and only affect specific use cases. This could for example be the correction of a single variable name. All breaking changes are marked in the [changelog](https://github.com/eumetsat/MetopDatasets.jl/blob/main/CHANGELOG.md).

docs/src/python.md

@@ -14,9 +14,9 @@ julia --version
 pip --version
 ```
 This guide is tested with the following versions 
-- Python 3.9.12
+- Python 3.12.8
 - julia version 1.11.1
-- pip 21.2.4
+- pip 24.2
 ### Installing python packages
 We use pip to install `juliacall` and `numpy`. We need `juliacall` to interface with julia and `numpy` is just needed to demonstrate compatibility with numpy arrays.
 
@@ -32,7 +32,7 @@ import juliacall
 # make separate module
 jl = juliacall.newmodule("MetopDatasetsPy") 
 jl.seval("import Pkg")
-jl.Pkg.add(url="https://github.com/eumetsat/MetopDatasets.jl")
+jl.Pkg.add("MetopDatasets")
 ```
 
 ### Example
@@ -48,7 +48,7 @@ jl.seval("using MetopDatasets")
 The dataset is simply read with `MetopDataset`. Only the metadata is read straight away. The variables can be read on demand.
 ```python
 test_file = "/tcenas/home/lupemba/Documents/data/IASI_xxx_1C_M01_20240819103856Z_20240819104152Z_N_C_20240819112911Z"
-ds = jl.MetopDataset(test_file)
+ds = jl.MetopDataset(test_file, maskingvalue=float('nan'))
 ```
 The dataset has a method equivalent to `__repr__` so the structure of the dataset can be shown easily. The julia `keys` function can be used to only list variable names.
 ```python

src/Instruments/IASI/DiskArrayExtra.jl

@@ -93,8 +93,3 @@ function DiskArrays.readblock!(disk_array::IasiWaveNumberDiskArray, aout,
     aout .= sample_width .* (number_of_first_sample .+ i_channel .- 2)
     return nothing
 end
-
-function CDM.dimnames(disk_array::IasiWaveNumberDiskArray)
-    spectrum_dim = get_field_dimensions(disk_array.record_type, :gs1cspect)[1]
-    return [spectrum_dim, RECORD_DIM_NAME]
-end

src/Instruments/IASI/IASI_dimensions.jl

@@ -75,6 +75,29 @@ end
 
 ############  IASI_SND_02  ################
 
+function _get_fixed_dimensions(T::Type{IASI_SND_02_V11})
+    dimensions_dict = Dict(
+        "lat_lon" => 2,
+        "cloud_formations" => 3,
+        "solar_sat_zenith_azimuth" => 4,
+        "xtrack_sounder_pixels" => 120
+    )
+    return dimensions_dict
+end
+
+function _get_fixed_dimensions(T::Type{IASI_SND_02_V10})
+    dimensions_dict = Dict(
+        "lat_lon" => 2,
+        "cloud_formations" => 3,
+        "solar_sat_zenith_azimuth" => 4,
+        "xtrack_sounder_pixels" => 120,
+        IASI_L2_V10_03_PRESSURE_DIM => 2,
+        "n_surface_temperatures" => 2,
+        "error_data_dims" => 2
+    )
+    return dimensions_dict
+end
+
 function get_dimensions(T::Type{<:IASI_SND_02},
         data_record_layouts::Vector{<:RecordLayout})::Dict{String, <:Integer}
     dimensions_dict = Dict{String, Integer}()
@@ -99,29 +122,35 @@ function get_dimensions(T::Type{<:IASI_SND_02},
         end
     end
 
-    # default others
-    dimensions_dict["lat_lon"] = 2 #The coordinates have opposite order of L1C!
-    dimensions_dict["cloud_formations"] = 3
-    dimensions_dict["solar_sat_zenith_azimuth"] = 4
-    dimensions_dict["xtrack_sounder_pixels"] = 120
+    merge!(dimensions_dict, _get_fixed_dimensions(T))
 
     return dimensions_dict
 end
 
 function get_field_dimensions(
-        T::Type{<:IASI_SND_02}, layout::FlexibleRecordLayout, field_name::Symbol)
+        T::Type{<:IASI_SND_02}, field_name::Symbol)
     res = String[]
 
     if !(fieldtype(T, field_name) <: Array)
         return res
     end
 
-    dimension_dict = get_dimensions(T, [layout])
+    dimension_dict = _get_fixed_dimensions(T)
     array_size = _get_array_size(T, field_name)
 
     for d in array_size
         if d isa Symbol
             push!(res, string(d))
+        elseif T <: IASI_SND_02_V10 && d == 2
+            if field_name == :pressure_levels_ozone
+                push!(res, IASI_L2_V10_03_PRESSURE_DIM)
+            elseif field_name == :surface_temperature
+                push!(res, "n_surface_temperatures")
+            elseif field_name == :data_sizes
+                push!(res, "error_data_dims")
+            else
+                push!(res, "lat_lon")
+            end
         else
             names = [dim_name for (dim_name, dim_val) in dimension_dict if dim_val == d]
             push!(res, string(first(names)))

src/Instruments/IASI/IASI_giadr_02.jl

@@ -2,15 +2,39 @@
 # License: MIT
 
 # Add IASI Level 2 meta data record.
-GIADR_IASI_SND_02_V11_format = joinpath(@__DIR__, "csv_formats/GIADR_IASI_SND_02_V11.csv")
+const IASI_L2_V10_03_PRESSURE_DIM = "n_o3_profiles"
 
 abstract type GIADR_IASI_SND_02 <: GlobalInternalAuxillary end
 
+GIADR_IASI_SND_02_V11_format = joinpath(@__DIR__, "csv_formats/GIADR_IASI_SND_02_V11.csv")
+GIADR_IASI_SND_02_V10_format = joinpath(@__DIR__, "csv_formats/GIADR_IASI_SND_02_V10.csv")
+
 eval(record_struct_expression(GIADR_IASI_SND_02_V11_format, GIADR_IASI_SND_02))
+eval(record_struct_expression(GIADR_IASI_SND_02_V10_format, GIADR_IASI_SND_02))
 
 get_instrument_subclass(::Type{<:GIADR_IASI_SND_02}) = 1
 
-function get_flexible_dim_fields(::Type{<:GIADR_IASI_SND_02})
+function get_giard_varnames(::Type{GIADR_IASI_SND_02_V11})
+    return (
+        :pressure_levels_temp,
+        :pressure_levels_humidity,
+        :pressure_levels_ozone,
+        :surface_emissivity_wavelengths,
+        :forli_layer_heights_co,
+        :forli_layer_heights_hno3,
+        :forli_layer_heights_o3,
+        :brescia_altitudes_so2)
+end
+
+function get_giard_varnames(::Type{GIADR_IASI_SND_02_V10})
+    return (
+        :pressure_levels_temp,
+        :pressure_levels_ozone,
+        :pressure_levels_humidity,
+        :surface_emissivity_wavelengths)
+end
+
+function get_flexible_dim_fields(::Type{GIADR_IASI_SND_02_V11})
     return Dict(
         :num_temperature_pcs => :NPCT,
         :num_pressure_levels_temp => :NLT,
@@ -25,7 +49,15 @@ function get_flexible_dim_fields(::Type{<:GIADR_IASI_SND_02})
         :forli_num_layers_o3 => :NL_O3)
 end
 
-function get_iasi_l2_flex_size(giard::GIADR_IASI_SND_02)
+function get_flexible_dim_fields(::Type{GIADR_IASI_SND_02_V10})
+    return Dict(
+        :num_pressure_levels_temp => :NLT,
+        :num_pressure_levels_ozone => :NLO,
+        :num_pressure_levels_humidity => :NLQ,
+        :num_surface_emissivity_wavelengths => :NEW)
+end
+
+function get_iasi_l2_flex_size(giard::T) where {T <: GIADR_IASI_SND_02}
     flex_size_prod = Dict{Symbol, Int64}()
     giard_size_fields = get_flexible_dim_fields(typeof(giard))
 
@@ -36,6 +68,11 @@ function get_iasi_l2_flex_size(giard::GIADR_IASI_SND_02)
     end
 
     # computed dimensions
+    _add_computed_dimension!(flex_size_prod, T)
+    return flex_size_prod
+end
+
+function _add_computed_dimension!(flex_size_prod::Dict, ::Type{GIADR_IASI_SND_02_V11})
     flex_size_prod[:NEVA_CO] = ceil(Int64, flex_size_prod[:NL_CO] / 2)
     flex_size_prod[:NEVE_CO] = flex_size_prod[:NEVA_CO] * flex_size_prod[:NL_CO]
     flex_size_prod[:NEVA_HNO3] = ceil(Int64, flex_size_prod[:NL_HNO3] / 2)
@@ -45,6 +82,21 @@ function get_iasi_l2_flex_size(giard::GIADR_IASI_SND_02)
     flex_size_prod[:NERRT] = Int64(flex_size_prod[:NPCT] * (flex_size_prod[:NPCT] + 1) / 2)
     flex_size_prod[:NERRW] = Int64(flex_size_prod[:NPCW] * (flex_size_prod[:NPCW] + 1) / 2)
     flex_size_prod[:NERRO] = Int64(flex_size_prod[:NPCO] * (flex_size_prod[:NPCO] + 1) / 2)
+    return nothing
+end
 
-    return flex_size_prod
+function _add_computed_dimension!(flex_size_prod::Dict, ::Type{GIADR_IASI_SND_02_V10})
+    return nothing
+end
+
+function get_field_dimensions(T::Type{GIADR_IASI_SND_02_V11}, field::Symbol)
+    return [string(get_raw_format_dim(T)[field][1])]
+end
+
+function get_field_dimensions(T::Type{GIADR_IASI_SND_02_V10}, field::Symbol)
+    if field == :pressure_levels_ozone
+        return [IASI_L2_V10_03_PRESSURE_DIM, string(get_raw_format_dim(T)[field][2])]
+    end
+
+    return [string(get_raw_format_dim(T)[field][1])]
 end

src/Instruments/IASI/IASI_records.jl

@@ -15,21 +15,28 @@ function data_record_type(header::MainProductHeader, product_type::Val{:IASI_xxx
     end
 end
 
-IASI_SND_02_V11_format = joinpath(@__DIR__, "csv_formats/IASI_SND_02_V11.csv")
+######### level 2###########
 
 abstract type IASI_SND_02 <: DataRecord end
+
+IASI_SND_02_V11_format = joinpath(@__DIR__, "csv_formats/IASI_SND_02_V11.csv")
+IASI_SND_02_V10_format = joinpath(@__DIR__, "csv_formats/IASI_SND_02_V10.csv")
+
 # create data structure, add description and scale factors
 eval(record_struct_expression(IASI_SND_02_V11_format, IASI_SND_02))
+eval(record_struct_expression(IASI_SND_02_V10_format, IASI_SND_02))
 
 function data_record_type(header::MainProductHeader, product_type::Val{:IASI_SND_02})::Type
     if header.format_major_version == 11
         return IASI_SND_02_V11
+    elseif header.format_major_version == 10
+        return IASI_SND_02_V10
     else
         error("No format found for format major version :$(header.format_major_version)")
     end
 end
 
-function get_flexible_dim_fields(::Type{<:IASI_SND_02})
+function get_flexible_dim_fields(::Type{IASI_SND_02_V11})
     return Dict(
         :nerr => :NERR,
         :co_nbr => :CO_NBR,
@@ -38,18 +45,26 @@ function get_flexible_dim_fields(::Type{<:IASI_SND_02})
     )
 end
 
+function get_flexible_dim_fields(::Type{IASI_SND_02_V10})
+    return Dict{Symbol, Symbol}(
+    )
+end
+
+_get_giard_type(T::Type{IASI_SND_02_V11}) = GIADR_IASI_SND_02_V11
+_get_giard_type(T::Type{IASI_SND_02_V10}) = GIADR_IASI_SND_02_V10
+
 function _get_flexible_dims_file(file_pointer::IO, T::Type{<:IASI_SND_02})
     pos = position(file_pointer)
     seekstart(file_pointer)
-    giard = read_first_record(file_pointer, GIADR_IASI_SND_02_V11)
+    giard = read_first_record(file_pointer, _get_giard_type(T))
     flexible_dims_file = get_iasi_l2_flex_size(giard)
     seek(file_pointer, pos)
     return flexible_dims_file
 end
 
-# IASI_SND_02 needs to have a fill value for VInteger and BitString
+# IASI_SND_02_V11  needs to have a fill value for VInteger and BitString
 function get_missing_value(
-        record_type::Type{<:IASI_SND_02}, ::Type{VInteger{T}}, field::Symbol) where {T}
+        record_type::Type{IASI_SND_02_V11}, ::Type{VInteger{T}}, field::Symbol) where {T}
 
     # Only set missing value for variable fields to avoid masking real data
     if !fixed_size(record_type, field)
@@ -60,7 +75,7 @@ function get_missing_value(
 end
 
 function get_missing_value(
-        record_type::Type{<:IASI_SND_02}, ::Type{BitString{N}}, field::Symbol) where {N}
+        record_type::Type{IASI_SND_02_V11}, ::Type{BitString{N}}, field::Symbol) where {N}
 
     # Only set missing value for variable fields to avoid masking real data
     if !fixed_size(record_type, field)
@@ -70,3 +85,7 @@ function get_missing_value(
         return nothing
     end
 end
+
+## 
+const IASI_L2_V10_ERROR_DATA_NAME = :error_data
+const IASI_L2_V10_ERROR_DATA_DESCRIPTION = "Contents depend on MDR.FLG_STER field (Data is not parsed)"