diff --git a/Project.toml b/Project.toml index b21b089..44cae07 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MetopDatasets" uuid = "0c26954c-4046-4b98-a13c-f9377ca4b9b7" authors = ["lupemba and contributors"] -version = "0.0.6" +version = "0.0.7" [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" diff --git a/README.md b/README.md index a120934..fa28bf1 100644 --- a/README.md +++ b/README.md @@ -1,17 +1,16 @@ -# Under development -This package is still in the early development phase. Note that the package was previously named **MetopNative.jl** and was hosted on the [EUMETSAT GitLab](https://gitlab.eumetsat.int/eumetlab/cross-cutting-tools/MetopNative.jl) - # MetopDatasets.jl - +[![documentation stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://eumetsat.github.io/MetopDatasets.jl/stable/) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://eumetsat.github.io/MetopDatasets.jl/dev/) [![Build Status](https://github.com/eumetsat/MetopDatasets.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/eumetsat/MetopDatasets.jl/actions/workflows/CI.yml?query=branch%3Amain) [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) -MetopDatasets.jl is a package for reading products from the [METOP satellites](https://www.eumetsat.int/our-satellites/metop-series) using the native binary format specified for each product. The METOP satellites are part of the EUMETSAT-POLAR-SYSTEM (EPS) and produce long series of near real-time weather and climate observation. Learn more and access the products on [EUMETSATs user-portal](https://user.eumetsat.int/dashboard). +MetopDatasets.jl is a package for reading products from the [METOP satellites](https://www.eumetsat.int/our-satellites/metop-series) using the native binary format specified for each product. The METOP satellites are part of the EUMETSAT-POLAR-SYSTEM (EPS) and have produced near real-time, global weather and climate observation since 2007. Learn more METOP and the data access on [EUMETSATs user-portal](https://user.eumetsat.int/dashboard). MetopDatasets.jl exports the `MetopDataset` API which is an implementation of the [CommonDataModel.jl](https://github.com/JuliaGeo/CommonDataModel.jl) interface and thus provides data access similar to e.g. [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl) and [GRIBDatasets.jl](https://github.com/JuliaGeo/GRIBDatasets.jl). -Only a subset of the METOP native formats are supported currently but we are contiguously adding formats. The goal is to support all publicly available native METOP products. +Only a subset of the METOP native formats are supported currently but we are continuously adding formats. The goal is to support all publicly available [native METOP products](https://data.eumetsat.int/extended?query=&filter=satellite__Metop&filter=availableFormats__Native). See [supported formats](https://eumetsat.github.io/MetopDatasets.jl/dev/#Supported-formats) for more information + +It is also possible to use MetopDatasets.jl from Python. See [section in documentation](https://eumetsat.github.io/MetopDatasets.jl/dev/python) for more information. ## Copyright and License This code is licensed under MIT license. See file LICENSE for details on the usage and distribution terms. @@ -21,357 +20,53 @@ This code is licensed under MIT license. See file LICENSE for details on the usa * [Jonas Wilzewski](mailto://jonas.wilzewski@eumetsat.int) - *Contributor* - [EUMETSAT](http://www.eumetsat.int) ## Installation -MetopDatasets.jl can be installed via Pkg and the url to the github.com repo. +MetopDatasets.jl can be installed via Pkg and the url to the GitHub repository. ```julia import Pkg Pkg.add(url="https://github.com/eumetsat/MetopDatasets.jl") ``` -## Uninstall -```julia -import Pkg -Pkg.rm("MetopDatasets") -``` - - -## Reference documents -The package is implemented based on -- [EPS Generic Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/pdf_gen_pfs_13e3f0feb7.pdf) -### Supported formats -- [ASCAT Level 1: Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/ASCAT_Level_1_Product_Format_V12_Annexe_50fe72d349.pdf) -- [ASCAT Level 2 Soil Moisture: Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/pdf_ten_0343_eps_ascatl2_pfs_f509981295.pdf) -- [IASI Level 1 Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/pdf_iasi_level_1_pfs_2105bc9ccf.pdf) - -### Formats under development -- [IASI Level 2 Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/pdf_ten_980760_eps_iasi_l2_f9511c26d2.pdf) - - - ## Examples - -### Read data from a Metop Native binary file -To read a Metop Native binary file: - +Open data set and list variables ```julia using MetopDatasets -ds = MetopDataset("ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nat") +ds = MetopDataset("ASCA_SZO_1B_M03_20230329063300Z_20230329063556Z_N_C_20230329081417Z"); +keys(ds) ``` REPL output: ``` -Dataset: -Group: / - -Dimensions - num_band = 3 - xtrack = 82 - atrack = 3264 - -Variables - record_start_time (3264) - Datatype: Dates.DateTime (Float64) - Dimensions: atrack - Attributes: - description = Record header start time - units = seconds since 2000-1-1 0:0:0 - - record_stop_time (3264) - Datatype: Dates.DateTime (Float64) - Dimensions: atrack - Attributes: - description = Record header stop time - units = seconds since 2000-1-1 0:0:0 - - degraded_inst_mdr (3264) - Datatype: UInt8 (UInt8) - Dimensions: atrack - Attributes: - description = Quality of MDR has been degraded from nominal due to an instrument degradation. - - degraded_proc_mdr (3264) - Datatype: UInt8 (UInt8) - Dimensions: atrack - Attributes: - description = Quality of MDR has been degraded from nominal due to a processing degradation. - - utc_line_nodes (3264) - Datatype: Dates.DateTime (Float64) - Dimensions: atrack - Attributes: - description = UTC time of line of nodes - units = seconds since 2000-1-1 0:0:0 - - abs_line_number (3264) - Datatype: Int32 (Int32) - Dimensions: atrack - Attributes: - description = Absolute (unique) counter for the line of nodes (from format version 12.0 onwards only) - - sat_track_azi (3264) - Datatype: Float64 (UInt16) - Dimensions: atrack - Attributes: - description = Azimuth angle bearing (range: 0 to 360) of nadir track velocity - scale_factor = 0.010000000000000002 - - as_des_pass (3264) - Datatype: UInt8 (UInt8) - Dimensions: atrack - Attributes: - description = Ascending/descending pass indicator - - swath_indicator (82 × 3264) - Datatype: UInt8 (UInt8) - Dimensions: xtrack × atrack - Attributes: - description = Swath (0=LEFT, 1=RIGHT) - - latitude (82 × 3264) - Datatype: Float64 (Int32) - Dimensions: xtrack × atrack - Attributes: - description = Latitude (-90 to 90 deg) - scale_factor = 1.0e-6 - - longitude (82 × 3264) - Datatype: Float64 (Int32) - Dimensions: xtrack × atrack - Attributes: - description = Longitude (0 to 360 deg) - scale_factor = 1.0e-6 - - sigma0_trip (3 × 82 × 3264) - Datatype: Float64 (Int32) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Sigma0 triplet, re-sampled to swath grid, for 3 beams (fore, mid, aft) - scale_factor = 1.0e-6 - - kp (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Kp for re-sampled sigma0 tripplet. Values between 0 and 1 - scale_factor = 0.0001 - - inc_angle_trip (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Incidence angle for re-sampled sigma0 tripplet. - scale_factor = 0.010000000000000002 - - azi_angle_trip (3 × 82 × 3264) - Datatype: Float64 (Int16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Incidence angle for re-sampled sigma0 tripplet. Values range from -180 to +180, where minus is west and plus is east. - scale_factor = 0.010000000000000002 - - num_val_trip (3 × 82 × 3264) - Datatype: UInt32 (UInt32) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Number of full resolution sigma0 values contributing to the re-sampled sigma0 tripplet. - - f_kp (3 × 82 × 3264) - Datatype: UInt8 (UInt8) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to the quality of the Kp estimate (0=NOMINAL, 1=NON-NOMINAL) - - f_usable (3 × 82 × 3264) - Datatype: UInt8 (UInt8) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to the usability of the sigma0 tripplet (0=GOOD, 1=USABLE, 2=NOT USABLE) - - f_f (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to non-nominal amount of input raw data to calculate echo corrections (value between 0 and 1 shows the fraction of original samples -affected) - scale_factor = 0.001 - - f_v (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to non enough amount of input raw data to calculate echo corrections (value between 0 and 1 shows the fraction of original samples affected) - scale_factor = 0.001 - - f_oa (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to lack of accuracy of orbit/atticute knowledge (value between 0 and 1 shows the fraction of original samples affected) - scale_factor = 0.001 - - f_sa (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to solar array reflection contamination (value between 0 and 1 shows the fraction of original samples affected) - scale_factor = 0.001 - - f_tel (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to non-nominal telemetry check results (value between 0 and 1 shows the fraction of original samples affected) - scale_factor = 0.001 - - f_ref (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to non-nominal raw echo correction reference functions (value between 0 and 1 shows the fraction of original samples affected) - scale_factor = 0.001 - - f_land (3 × 82 × 3264) - Datatype: Float64 (UInt16) - Dimensions: num_band × xtrack × atrack - Attributes: - description = Flag related to presence of land in the re-sampled sigma0 tripllet (value between 0 and 1 shows the fraction of original samples affected) - scale_factor = 0.001 - -Global attributes - product_name = ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z - parent_product_name_1 = ASCA_xxx_1A_M01_20190109125700Z_20190109143859Z_N_O_20190109134742Z - parent_product_name_2 = xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx - parent_product_name_3 = xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx - parent_product_name_4 = xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx - instrument_id = ASCA - instrument_model = 1 - product_type = SZR - processing_level = 1B - spacecraft_id = M01 - sensing_start = 2019-01-09T12:57:00 - sensing_end = 2019-01-09T14:38:58 - sensing_start_theoretical = 2019-01-09T12:57:00 - sensing_end_theoretical = 2019-01-09T14:39:00 - processing_centre = CGS1 - processor_major_version = 10 - processor_minor_version = 0 - format_major_version = 12 - format_minor_version = 0 - processing_time_start = 2019-01-09T13:48:16 - processing_time_end = 2019-01-09T15:27:30 - processing_mode = N - disposition_mode = O - receiving_ground_station = SVL - receive_time_start = 2019-01-09T12:55:25 - receive_time_end = 2019-01-09T15:24:35 - orbit_start = 32745 - orbit_end = 32746 - actual_product_size = 26618899 - state_vector_time = 2019-01-09T12:27:10 - semi_major_axis = 7204713107 - eccentricity = 1312 - inclination = 98722 - perigee_argument = 67970 - right_ascension = 70895 - mean_anomaly = 292186 - x_position = -5122760992 - y_position = 5061020490 - z_position = 2084232 - x_velocity = 1170132 - y_velocity = 1168504 - z_velocity = 7355678 - earth_sun_distance_ratio = 983395 - location_tolerance_radial = 0 - location_tolerance_crosstrack = 0 - location_tolerance_alongtrack = 0 - yaw_error = 0 - roll_error = 0 - pitch_error = 0 - subsat_latitude_start = 71792 - subsat_longitude_start = -24448 - subsat_latitude_end = 69881 - subsat_longitude_end = -52970 - leap_second = 0 - leap_second_utc = - total_records = 3283 - total_mphr = 1 - total_sphr = 1 - total_ipr = 9 - total_geadr = 1 - total_giadr = 0 - total_veadr = 5 - total_viadr = 2 - total_mdr = 3264 - count_degraded_inst_mdr = 0 - count_degraded_proc_mdr = 0 - count_degraded_inst_mdr_blocks = 0 - count_degraded_proc_mdr_blocks = 0 - duration_of_product = 6118000 - milliseconds_of_data_present = 6056250 - milliseconds_of_data_missing = 0 - subsetted_product = F +("record_start_time", "record_stop_time", "degraded_inst_mdr", "degraded_proc_mdr", "utc_line_nodes", "abs_line_number", "sat_track_azi", "as_des_pass", "swath_indicator", "latitude", "longitude", "sigma0_trip", "kp", "inc_angle_trip", "azi_angle_trip", "num_val_trip", "f_kp", "f_usable", "f_land", +"lcr", "flagfield") ``` - -The variables can be loaded from MetopDataset by indexing the dataset. The variable then works as a lazy array loading the data on indexing: +Display variable information ```julia -ds["latitude"][2,4] +ds["latitude"] ``` REPL output: ``` -66.707944 +latitude (42 × 48) + Datatype: Float64 (Int32) + Dimensions: xtrack × atrack + Attributes: + description = Latitude (-90 to 90 deg) + scale_factor = 1.0e-6 ``` -It is also possible to load the complete array +Load the complete variable ```julia ds["latitude"][:,:] ``` REPL output: ``` -82×3264 Matrix{Float64}: - 66.862 66.7841 66.706 66.6276 66.5489 … 65.6304 65.5487 65.4667 65.3844 65.3019 - 66.9431 66.865 66.7866 66.7079 66.629 65.7077 65.6257 65.5434 65.4609 65.3782 - ⋮ ⋱ ⋮ - 74.2237 74.1153 74.007 73.8986 73.7903 … 72.5493 72.4409 72.3325 72.2241 72.1157 - 74.2342 74.1259 74.0175 73.9092 73.8008 72.5597 72.4513 72.3429 72.2345 72.1261 -``` -Data from the main product header is accessed as attributes. -```julia -ds.attrib["instrument_id"] -``` -REPL output: -``` -"ASCA" -``` -### Convert a Metop Native binary file to netCDF - -A Metop Native binary file can be converted to netCDF using the [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl) package. This -is possible because both MetopDatasets.jl and NCDatasets.jl implement the [CommonDataModel.jl](https://github.com/JuliaGeo/CommonDataModel.jl) interface. - -```julia -using MetopDatasets -using NCDatasets - -input_file = "ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nat" -output_file = "ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nc" - -ds_native = MetopDataset(input_file) -ds_nc = NCDataset(output_file, "c") -NCDatasets.write(ds_nc, ds_native) - -close(ds_native) -close(ds_nc) -``` -It is also possible to use the safe `do` syntax that ensures the files are closed correctly even in the case of exceptions. - -```julia -using MetopDatasets -using NCDatasets - -input_file = "ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nat" -output_file = "ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nc" - -MetopDataset(input_file) do ds_native - NCDataset(output_file, "c") do ds_nc - NCDatasets.write(ds_nc, ds_native) - end -end -``` +42×48 Matrix{Float64}: + -33.7308 -33.949 … -43.7545 -43.9721 + -33.6969 -33.9152 -43.7252 -43.9429 + -33.6624 -33.8808 -43.695 -43.9127 + -33.6274 -33.8458 -43.6639 -43.8818 + ⋮ ⋱ + -30.1606 -30.3748 -39.9343 -40.1446 + -30.0909 -30.3049 … -39.8538 -40.0638 + -30.0206 -30.2344 -39.7726 -39.9823 +``` +See documentation page for more information. \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index d897df3..89d69c5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -23,8 +23,12 @@ makedocs(; ), pages = [ "Introduction" => "index.md", - "Public API" => "public_api.md", - "Internal API" => "internal_api.md" + "Use with Python" => "python.md", + "Examples" => [ + "ASCAT" => "ASCAT.md", + "IASI" => "IASI.md" + ], + "Full API" => "full_api.md" ] ) diff --git a/docs/src/ASCAT.md b/docs/src/ASCAT.md new file mode 100644 index 0000000..d11de07 --- /dev/null +++ b/docs/src/ASCAT.md @@ -0,0 +1,273 @@ +## ASCAT + +The Advanced Scatterometer (ASCAT) is an instrument on the METOP satellites. It is a fixed beam C-band scatterometer with 6 antenna beams used to retrieve wind vectors over the ocean and soil moisture over the land. The level 1 products are normalised backscatter coefficients (sigma0) that the instrument measures directly. The level 2 products are soil moisture and wind vectors derived from the backscatter. The wind vector product does not use the native binary format so they are not described further here. It should be noted that the winds are available as NetCDF. +For more information see [ASCAT User Guide](https://user.eumetsat.int/s3/eup-strapi-media/ASCAT_User_Guide_1618a3d741.pdf) + +All data used in examples are available from [EUMETSAT data store](https://data.eumetsat.int/extended?query=&filter=instrument__ASCAT&filter=availableFormats__Native) + +## SZF beam statistics +The level 1B SZF files contain full resolution backscatter from the 6 antenna beams. See the [ASCAT level 1B SZF product page](https://data.eumetsat.int/product/EO:EUM:DAT:METOP:ASCSZF1B) for more information. +This example shows how we can get the 5% and 95% quantile of the backscatter from each beam only considering the observations with no flags. It should be noted that observations over land are are flagged with a land-flag so theses statistics only consider oceans. + +Read and inspect the variables: +```julia +using MetopDatasets, Statistics +ds = MetopDataset("ASCA_SZF_1B_M03_20241217091500Z_20241217105659Z_N_O_20241217105652Z.nat") + +# Read arrays immediately +beam_numbers = Array(ds["beam_number"]) +sigma0 = Array(ds["sigma0_full"]) +flags = Array(ds["flagfield"]) + +@show size(beam_numbers) +@show size(sigma0) +@show size(flags) + +no_flags = flags .== 0 +println("Number of observations with no flag: ", sum(no_flags)) +``` +REPL output: +``` +size(beam_numbers) = (43236,) +size(sigma0) = (192, 43236) +size(flags) = (192, 43236) +Number of observations with no flag: 5948489 +``` +We can see that `beam_numbers` only have 1 dimension but `sigma0` and `flags` have two. This is because there is only one beam number per data record but each data record contain a line with 192 observations across the swath. +Printing the quantiles for each beam. +```julia +for i in 1:6 # same as `i in sort(unique(beam_numbers))` + in_beam = beam_numbers .== i + sigma0_beam = sigma0[:, in_beam] + no_flags_beam = no_flags[:, in_beam] + + sigma_q05 = quantile(sigma0_beam[no_flags_beam],0.05) + sigma_q95 = quantile(sigma0_beam[no_flags_beam],0.95) + + rounded_range = round.((sigma_q05,sigma_q95); digits=2) + + println("Beam: ", i, ", Sigma0 (q05,q95) range: ", rounded_range ) +end +``` +REPL output: +``` +Beam: 1, Sigma0 (q05,q95) range: (-29.84, -12.31) +Beam: 2, Sigma0 (q05,q95) range: (-27.22, -5.63) +Beam: 3, Sigma0 (q05,q95) range: (-29.94, -12.76) +Beam: 4, Sigma0 (q05,q95) range: (-28.75, -12.96) +Beam: 5, Sigma0 (q05,q95) range: (-25.5, -5.2) +Beam: 6, Sigma0 (q05,q95) range: (-28.4, -12.02) +``` +We can now see the sigma0 quintile ranges for each beam. + +## SZR and SZO plot +The level 1B SZR files contain backscatter values from the 6 ASCAT beams, resampled to a common 12.5 x 12.5 km grid. There are three beams imaging each side of the satellite so this results in 3 backscatter measurements for each grid point. The level 1B SZO files are identical but resampled to a 25 x 25 km grid instead. See [ASCAT level 1B SZR product page](https://data.eumetsat.int/product/EO:EUM:DAT:METOP:ASCSZR1B) and [ASCAT level 1B SZO product page](https://data.eumetsat.int/product/EO:EUM:DAT:METOP:ASCSZO1B) for more information. + +These examples will work to both plot SZR and SZO data. The examples are made using the following packages. +``` +[13f3f980] CairoMakie v0.12.16 +[db073c08] GeoMakie v0.7.9 +``` +Note that the plot can become interactive by replacing CarioMakie with GLMakie. See [Makie documentation](https://docs.makie.org/) for more information on plotting. +### Global plot of SZR or SZO +This example plots the backscatter from and entire orbit. The backscatter is just plotted for one beam direction. +```julia +using MetopDatasets +using CairoMakie, GeoMakie, Statistics + +# Open the dataset. +ds = MetopDataset("ASCA_SZR_1B_M01_20241217081500Z_20241217095658Z_N_O_20241217090832Z.nat"); + +# Select which beam to plot +beam_names = ["fore-beam", "mid-beam", "aft-beam"] +beam_index = 3 + +# read variables for beam +latitude = Array(ds["latitude"]) +longitude = Array(ds["longitude"]) +sigma0 = ds["sigma0_trip"][beam_index,:,:] +usable_flag = ds["f_usable"][beam_index,:,:] + +# change longitude to be -180 to 180 +longitude[longitude .> 180] .-= 360 + +# check flags +usable = usable_flag .!= 2 + +# get sigma0 ranges +sigma_q05 = quantile(sigma0[usable],0.05) +sigma_q95 = quantile(sigma0[usable],0.95) + +# plot the data +fig = let + + # Create figure and axis + fig = Figure() + ax = GeoAxis(fig[1, 1], + title = "Backscatter for " * beam_names[beam_index], + xlabel = "longitude", + ylabel = "latitude") + + # plot "good" and "usable" data with color map + scatter!(ax, longitude[usable], latitude[usable], + color = sigma0[usable], colorrange = (sigma_q05,sigma_q95), markersize = 1) + + # plot "bad" data as red + scatter!(ax, longitude[.!usable], latitude[.!usable], + color = :red, markersize = 1) + + # Add colorbar + Colorbar(fig[1,2], colorrange = (sigma_q05,sigma_q95), label="dB") + + # Add coastlines + lines!(ax, GeoMakie.coastlines()) + fig +end +``` +![Global SZR plot](global_szr_plot.png) + +The plot shows the left and right swath ASCAT mapped for an orbit. + + +### Plotting a subset of SZR or SZO +This example shows how to plot the backscatter for a small subset of the SZR product. +```julia +using MetopDatasets +using CairoMakie, GeoMakie, Statistics + +ds = MetopDataset("ASCA_SZR_1B_M01_20241217081500Z_20241217095658Z_N_O_20241217090832Z.nat"); + +# Select subset record subset to plot +record_subset = 150:350 + +# define the two swaths. +grid_nodes = ds.dim["xtrack"] +swath_width = div(grid_nodes,2) +left_swath = 1:swath_width +right_swath = (swath_width+1):grid_nodes + +# read variables for right_swath of the subset. +latitude = ds["latitude"][right_swath, record_subset] +longitude = ds["longitude"][right_swath, record_subset] +sigma0 = ds["sigma0_trip"][:, right_swath, record_subset] +usable_flag = ds["f_usable"][:, right_swath, record_subset] + +# change longitude to be -180 to 180 +longitude[longitude .> 180] .-= 360 + +# check flags +usable = usable_flag .!= 2 + +# get sigma0 ranges +sigma_q05 = quantile(sigma0[usable],0.05) +sigma_q95 = quantile(sigma0[usable],0.95) + +# plot the data +beam_names = ["fore-beam", "mid-beam", "aft-beam"] +fig = let + fig = Figure() + + # plot each beam + for i in 1:3 + sigma_i = sigma0[i,:,:] + usable_i = usable[i,:,:] + + # Create new axis + ax = GeoAxis(fig[1, i], + title = beam_names[i], + xlabel = "longitude", + ylabel = "latitude", + limits = (extrema(longitude), extrema(latitude)) + ) + + # plot "good" and "usable" data with color map + scatter!(ax, longitude[usable_i], latitude[usable_i], + color = sigma_i[usable_i], colorrange = (sigma_q05,sigma_q95), markersize = 6) + + # plot "bad" data as red + scatter!(ax, longitude[.!usable_i], latitude[.!usable_i], + color = :red, markersize = 6) + + # add coastlines + lines!(ax, GeoMakie.coastlines(), color=:black) + end + + # add colorbar + Colorbar(fig[1,4], colorrange = (sigma_q05,sigma_q95), label="dB") + + fig +end +``` +![Subset SZR plot](subset_szr_plot.png) + +The subset shows the measured backscatter from the right swath over the mediterranean and across Italy. The red spots shows the measurements flagged as "bad". This is probably caused by radio frequency interference from ground equipment resulting in a low signal-to-noise ratio. + +## Soil moisture example +The ASCAT Level 2 SMR product contains relative soil moisture on a 12.5 x 12.5 km grid. The SMO product is identical but resampled to a 25 x 25 km grid instead. +See [ASCAT level 2 SMR product page](https://data.eumetsat.int/product/EO:EUM:DAT:METOP:SOMO12) and [ASCAT level 2 SMO product page](https://data.eumetsat.int/product/EO:EUM:DAT:METOP:SOMO25) for more information. + +These examples will work to plot both SMR and SMO data. The examples are made using the following packages. +``` +[13f3f980] CairoMakie v0.12.16 +[db073c08] GeoMakie v0.7.9 +``` +### Plotting a subset of SMR or SMO +This example shows how to plot the soil moisture for from a SMR or SMO product. +```julia +using MetopDatasets +using CairoMakie, GeoMakie, Dates + +ds = MetopDataset("ASCA_SMR_02_M01_20241217081500Z_20241217095658Z_N_O_20241217090845Z.nat"); + +# Select subset record subset to plot +record_subset = 150:350 + +# define the two swaths. +grid_nodes = ds.dim["xtrack"] +swath_width = div(grid_nodes,2) +left_swath = 1:swath_width +right_swath = (swath_width+1):grid_nodes + +# get sensing date +sensing_date = Date(DateTime(ds.attrib["sensing_start"])) + +# read variables for right_swath of the subset. +latitude = ds["latitude"][right_swath, record_subset] +longitude = ds["longitude"][right_swath, record_subset] +soil_moisture = ds["soil_moisture"][right_swath, record_subset] +flags = ds["aggregated_quality_flag"][ right_swath, record_subset] + +# change longitude to be -180 to 180 +longitude[longitude .> 180] .-= 360 + +# The out of range soil moisture is flagged with 255 +usable = flags .!= 255 + +# Set the unusable soil moisture to NaN +soil_moisture[.!usable] .= NaN + +fig = let + fig = Figure() + + # Create new axis + ax = GeoAxis(fig[1, 1], + title = "Soil moisture $(sensing_date)", + xlabel = "longitude", + ylabel = "latitude", + limits = (extrema(longitude), extrema(latitude)) + ) + + # plot with scatter + scatter!(ax, longitude[:], latitude[:], color=soil_moisture[:], colorrange=(0.0,100.0)) + + # add coastlines + lines!(ax, GeoMakie.coastlines(), color=:black) + + # add colorbar + Colorbar(fig[1,2], colorrange = (0.0,100.0), label="soil moisture %") + + fig +end +``` +![Subset SZR plot](subset_smr_plot.png) +The plot shows the retrieved soil moisture. The filtering on flags in this example is very rough. It might be good to have a close look at the individual flag variables depending on the application. \ No newline at end of file diff --git a/docs/src/IASI.md b/docs/src/IASI.md new file mode 100644 index 0000000..12f4817 --- /dev/null +++ b/docs/src/IASI.md @@ -0,0 +1,208 @@ +## IASI + +The Infrared Atmospheric Sounding Interferometer (IASI) is an instrument on the METOP satellites. It is a hyper-spectral sensor measuring the radiation from the atmosphere and earth in 8461 spectral channels. These measurements are used to derive used to derive a plethora of geophysical variables (e.g. temperature and humidity profiles). This makes IASI a key data source for numerical weather prediction (NWP) and applications in atmospheric chemistry and monitoring of essential climate variables. +See [ASI Level 1: Product Guide](https://user.eumetsat.int/s3/eup-strapi-media/pdf_iasi_pg_487c765315.pdf) and [IASI Level 2: Product Guide](https://user.eumetsat.int/s3/eup-strapi-media/IASI_Level_2_Product_Guide_8f61a2369f.pdf) for more information. + + +## Static plot of L1C spectra +This example is made using the following packages. +``` +[13f3f980] CairoMakie v0.12.16 +[db073c08] GeoMakie v0.7.9 +``` + +The key variable is the "gs1cspect" which contains the radiance spectrum measured by the IASI instrument. The spectrum from a full orbit is almost 2 GB of data. In this example we will just load one data record with observation locations and plot the spectrum of two observations. The spectra are converted from radiances to brightness temperature since this is often most convenient when interpreting the IASI spectra. The cloud cover is also read from the file and included as legends on the plot of the spectra. Cloud cover is essential to understanding the spectra and normally only cloud free observations are assimilated in NWP. + +```julia +using MetopDatasets +using CairoMakie, GeoMakie + +ds = MetopDataset("IASI_xxx_1C_M01_20240925202059Z_20240925220258Z_N_O_20240925211316Z.nat"); + +# Select a single data record +data_record_index = 670 + +# Select 2 points to plot the full spectrum +selected_points = CartesianIndex.([(2,15), (1,27)]) +selected_colors = [:red, :black] + +# Read the geolocation of the data record +longitude, latitude = let + lon_lat = ds["ggeosondloc"][:,:,:,data_record_index] + lon_lat[1,:,:], lon_lat[2,:,:] +end + +# read spectrum as brightness temperature +wavenumber = ds["spectra_wavenumber"][:, data_record_index] +selected_T_b = let + spectra = ds["gs1cspect"][:,:,:,data_record_index] + # convert select spectra to brightness temperature + [brightness_temperature.(spectra[:,i] , wavenumber) + for i in selected_points] +end + +# Read cloud cover for the selected points +selected_cloud_cover = ds["geumavhrr1bcldfrac"][:,:,data_record_index][selected_points] + +# plot the data +fig = let + + fig = Figure() + + # axis to plot geolocation + ax = GeoAxis(fig[1, 1], + title = "Observations", + xlabel = "longitude", + ylabel = "latitude", + limits = (extrema(longitude), extrema(latitude))) + + # plot all observations from data record in gray + scatter!(ax, longitude[:], latitude[:], color=:gray) + # plot selected observations in color + scatter!(ax, longitude[selected_points], latitude[selected_points], + color=selected_colors, marker=:xcross, markersize = 15) + + # Add coastlines + lines!(ax, GeoMakie.coastlines()) + + # Plot the selected spectra + ax2 = Axis(fig[2, 1], + title = "Selected spectra", + xlabel = "Wavenumber (cm-1)", + ylabel = "Brightness Temperature (K)", + limits=((500,3500), (150,300))) + + for i in eachindex(selected_points) + cloud_cover_i = Int(selected_cloud_cover[i]) + lines!(ax2, wavenumber./100, selected_T_b[i], color = selected_colors[i], + label = "$(cloud_cover_i)% cloud cover") + end + + # Add legends + axislegend(ax2, position = :rb) + + fig +end +``` +![Static IASI spectrum](static_IASI.png) +The top plot shows a row of observations west of the French coast. Two observation have been marked with a colored X. The spectrum of these observations are shown below. + + + +## Advanced interactive plots of L1C spectra +This example shows how the IASI spectra can be shown interactively on a responsive map. The example relies on [Tyler.jl](https://makieorg.github.io/Tyler.jl/v0.2.0/) and [GLMakie.jl](https://docs.makie.org/stable/explanations/backends/glmakie#glmakie). It is recommend to first try some examples from the Tyler documentation if you are new to this package. +``` +[e9467ef8] GLMakie v0.10.16 +[e170d443] Tyler v0.2.0 +``` +This example will use Tyler to create an interactive background map similar to google maps in the top panel of the figure. Then all the locations of the observation will be plotted on the map colored based on the cloud cover. The bottom panel will display the spectrum of a single observation. Clicking an observation in the top panel will load that spectrum and plot it in the bottom panel. The loading of the spectrum will be done lazily to conserve memory. +```julia +using Tyler, GLMakie, MetopDatasets + +ds = MetopDataset("IASI_xxx_1C_M01_20240925202059Z_20240925220258Z_N_O_20240925211316Z.nat"); + +# read geolocation points and data shape +pts, pts_original_size = let + lon_lat = ds["ggeosondloc"][:,:,:,:] + + # make sure longitude is from -180 to 180 + lon_lat[1, :,:,:][ lon_lat[1,:,:,:] .>180] .-= 360 + + # convert the points to tuples + lon_lat = tuple.(lon_lat[1, :,:,:],lon_lat[2, :,:,:]) + + # store the original shape of the points + pts_original_size = size(lon_lat) + + # Flatten the points and convert them to web_mercator (the coordinate system used by Tyler) + MT = Tyler.MapTiles + pts = [Point2f(MT.project(lon_lat[i], MT.wgs84, MT.web_mercator)) + for i in eachindex(lon_lat)] + + pts, pts_original_size +end; + +# read cloud fraction +cloud_fraction = Float32.(ds["geumavhrr1bcldfrac"][:]); + +# helper function to read the spectrum for a single point +function read_spectrum_pts(ds, index::CartesianIndex) + # read spectrum and wavenumber + spectrum = ds["gs1cspect"][:,Tuple(index)...] + wavenumber = ds["spectra_wavenumber"][:, Tuple(index)[end]] + + # covert to brightness temperature + T_B = brightness_temperature.(spectrum, wavenumber) + wavenumber_cm = wavenumber./100 + + # join brightness temperature and wavenumber to points + spectrum_pts = Point2f.(tuple.(wavenumber_cm, T_B)) + return spectrum_pts +end + +# read an initial spectrum +spectrum_pts = read_spectrum_pts(ds, CartesianIndex(1,1,1)); + +# create the inter active plot. +fig = let + fig = Figure() + + # select background map and initial zoom + provider = Tyler.TileProviders.Esri(:WorldImagery); + extent = Tyler.Extent(X = (-10, 10), Y = (-10, 10)); + + # create background map + ax1 = Axis(fig[1,1]) + m = Tyler.Map(extent; provider, figure=fig, axis=ax1); + wait(m); + + # Plot observation points with cloud cover + objscatter = scatter!(ax1, pts, color = cloud_fraction, + colorrange = (0,100), colormap=:grays, markersize=15) + # hack from https://github.com/MakieOrg/Tyler.jl/issues/109 + translate!(objscatter, 0, 0, 10) + + # Plot a red cross on top of the selected point + selected_point = Observable(pts[1:1]) + selected_scatter = scatter!(ax1, selected_point, + color = :red, markersize=10, marker =:xcross) + # hack from https://github.com/MakieOrg/Tyler.jl/issues/109 + translate!(selected_scatter, 0, 0, 11) + + # Add colorbar + Colorbar(fig[1, 2], limits = (0.0,100.0), colormap =:grays, label = "Cloud Fraction") + hidedecorations!(ax1) + + # Create the second plot for the spectrum + ax2 = Axis(fig[2, 1], + title = "Observation spectrum", + xlabel = "Wavenumber (cm-1)", + ylabel = "Brightness Temperature (K)") + + # plot the spectrum + spectrum_observable = Observable(spectrum_pts) + lines!(ax2, spectrum_observable) + + # Add event handler to update the plot when the user click on a new observation + obs_func = on(events(m.axis).mousebutton) do event + if event.button == Mouse.left && event.action == Mouse.press + plt, i = pick(m.axis) + if plt == objscatter # check if an observation was clicked + # get the CartesianIndex + cartesian_i = CartesianIndices(pts_original_size)[i] + # load the selected spectrum and update the spectrum plot + spectrum_observable[] = read_spectrum_pts(ds, cartesian_i) + # update the red x + selected_point[] = pts[i:i] + end + end + end + + fig +end +``` +![Interactive IASI spectrum](interactive_IASI.png) +It is now possible to interactively explore the nearly 100 000 observations from an obit of IASI with the background map giving important context. + +## Level 2 profiles +MetopDatasets.jl does not support IASI level 2 products yet. diff --git a/docs/src/full_api.md b/docs/src/full_api.md new file mode 100644 index 0000000..59c9b2b --- /dev/null +++ b/docs/src/full_api.md @@ -0,0 +1,16 @@ + +## Public functions and types + +```@autodocs +Modules = [MetopDatasets] +Private = false +Order = [:function, :type, :macro] +``` + + +## Internal functions and types +```@autodocs +Modules = [MetopDatasets] +Public = false +Order = [:function, :type, :macro] +``` \ No newline at end of file diff --git a/docs/src/global_szr_plot.png b/docs/src/global_szr_plot.png new file mode 100644 index 0000000..a108001 Binary files /dev/null and b/docs/src/global_szr_plot.png differ diff --git a/docs/src/index.md b/docs/src/index.md index 446dacd..c1d7d0b 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,14 +1,349 @@ ```@meta CurrentModule = MetopDatasets ``` +# MetopDatasets.jl -# MetopDatasets +MetopDatasets.jl is a package for reading products from the [METOP satellites](https://www.eumetsat.int/our-satellites/metop-series) using the native binary format specified for each product. The METOP satellites are part of the EUMETSAT-POLAR-SYSTEM (EPS) and have produced near real-time, global weather and climate observation since 2007. Learn more and access the products on [EUMETSATs user-portal](https://user.eumetsat.int/dashboard). -Documentation for [MetopDatasets](https://github.com/eumetsat/MetopDatasets.jl). +MetopDatasets.jl exports the `MetopDataset` API which is an implementation of the [CommonDataModel.jl](https://github.com/JuliaGeo/CommonDataModel.jl) interface and thus provides data access similar to e.g. [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl) and [GRIBDatasets.jl](https://github.com/JuliaGeo/GRIBDatasets.jl). -## Copyright and License +Only a subset of the METOP native formats are supported currently, but we are continuously adding formats. The goal is to support all publicly available [native METOP products](https://data.eumetsat.int/extended?query=&filter=satellite__Metop&filter=availableFormats__Native). See [Supported formats](@ref) for more information. -This code is licensed under MIT license. See file LICENSE for details on the usage and distribution terms. +It is also possible to use MetopDatasets.jl from Python. See [Use with Python](@ref) for more information. -## More documentation is under way +## Quick start +This section gives a very short overview of the core functionalities. 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. + +```julia +import Pkg +Pkg.add(url="https://github.com/eumetsat/MetopDatasets.jl#main") +``` + +### Read data from a Metop Native binary file +Read a Metop Native binary file and display meta data: + +```julia +using MetopDatasets +ds = MetopDataset("ASCA_SZO_1B_M03_20230329063300Z_20230329063556Z_N_C_20230329081417Z") +``` +REPL output: +``` +Dataset: +Group: / + +Dimensions + num_band = 3 + xtrack = 42 + atrack = 48 + +Variables + record_start_time (48) + Datatype: Dates.DateTime (Float64) + Dimensions: atrack + Attributes: + description = Record header start time + units = seconds since 2000-1-1 0:0:0 + + record_stop_time (48) + Datatype: Dates.DateTime (Float64) + Dimensions: atrack + Attributes: + description = Record header stop time + units = seconds since 2000-1-1 0:0:0 + + degraded_inst_mdr (48) + Datatype: UInt8 (UInt8) + Dimensions: atrack + Attributes: + description = Quality of MDR has been degraded from nominal due to an instrument degradation. + + degraded_proc_mdr (48) + Datatype: UInt8 (UInt8) + Dimensions: atrack + Attributes: + description = Quality of MDR has been degraded from nominal due to a processing degradation. + + utc_line_nodes (48) + Datatype: Dates.DateTime (Float64) + Dimensions: atrack + Attributes: + description = UTC time of line of nodes + units = seconds since 2000-1-1 0:0:0 + + abs_line_number (48) + Datatype: Int32 (Int32) + Dimensions: atrack + Attributes: + description = Absolute (unique) counter for the line of nodes (from format version 12.0 onwards only) + + sat_track_azi (48) + Datatype: Float64 (UInt16) + Dimensions: atrack + Attributes: + description = Azimuth angle bearing (range: 0 to 360) of nadir track velocity + scale_factor = 0.010000000000000002 + + as_des_pass (48) + Datatype: UInt8 (UInt8) + Dimensions: atrack + Attributes: + description = Ascending/descending pass indicator (0=DESC, 1=ASC) + + swath_indicator (42 × 48) + Datatype: UInt8 (UInt8) + Dimensions: xtrack × atrack + Attributes: + description = Swath (0=LEFT, 1=RIGHT) + + latitude (42 × 48) + Datatype: Float64 (Int32) + Dimensions: xtrack × atrack + Attributes: + description = Latitude (-90 to 90 deg) + scale_factor = 1.0e-6 + + longitude (42 × 48) + Datatype: Float64 (Int32) + Dimensions: xtrack × atrack + Attributes: + description = Longitude (0 to 360 deg) + scale_factor = 1.0e-6 + + sigma0_trip (3 × 42 × 48) + Datatype: Float64 (Int32) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Sigma0 triplet, re-sampled to swath grid, for 3 beams (fore, mid, aft) + scale_factor = 1.0e-6 + + kp (3 × 42 × 48) + Datatype: Float64 (UInt16) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Kp for re-sampled sigma0 triplet. Values between 0 and 1 + scale_factor = 0.0001 + + inc_angle_trip (3 × 42 × 48) + Datatype: Float64 (UInt16) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Incidence angle for re-sampled sigma0 triplet. + scale_factor = 0.010000000000000002 + + azi_angle_trip (3 × 42 × 48) + Datatype: Float64 (Int16) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Azimuth angle of the up-wind direction for a given measurement triplet (range: -180 to +180, where minus is west and plus is east with respect to North) + scale_factor = 0.010000000000000002 + + num_val_trip (3 × 42 × 48) + Datatype: UInt32 (UInt32) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Number of full resolution sigma0 values contributing to the re-sampled sigma0 triplet. + + f_kp (3 × 42 × 48) + Datatype: UInt8 (UInt8) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Flag related to the quality of the Kp estimate (0=NOMINAL, 1=NON-NOMINAL) + + f_usable (3 × 42 × 48) + Datatype: UInt8 (UInt8) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Flag related to the usability of the sigma0 triplet (0=GOOD, 1=USABLE, 2=NOT USABLE) + + f_land (3 × 42 × 48) + Datatype: Float64 (UInt16) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Flag related to presence of land in the re-sampled sigma0 triplet (based on land mask; value between 0 and 1 +shows the fraction of original samples affected) + scale_factor = 0.001 + + lcr (3 × 42 × 48) + Datatype: Float64 (UInt16) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Land Contamination Ratio estimate (based on SRF) + scale_factor = 0.0001 + + flagfield (3 × 42 × 48) + Datatype: UInt32 (UInt32) + Dimensions: num_band × xtrack × atrack + Attributes: + description = Flag field containing quality information + +Global attributes + product_name = ASCA_SZO_1B_M03_20230329063300Z_20230329063556Z_N_C_20230329081417Z + parent_product_name_1 = ASCA_xxx_1A_M03_20230329063300Z_20230329063559Z_N_C_20230329081221Z + parent_product_name_2 = xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx + parent_product_name_3 = xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx + parent_product_name_4 = xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx + instrument_id = ASCA + instrument_model = 1 + product_type = SZO + processing_level = 1B + spacecraft_id = M03 + sensing_start = 2023-03-29T06:33:00 + sensing_end = 2023-03-29T06:35:56 + sensing_start_theoretical = 2023-03-29T06:06:00 + sensing_end_theoretical = 2023-03-29T07:48:00 + processing_centre = CGS2 + processor_major_version = 11 + processor_minor_version = 3 + format_major_version = 13 + format_minor_version = 1 + processing_time_start = 2023-03-29T08:14:17 + processing_time_end = 2023-03-29T08:14:25 + processing_mode = N + disposition_mode = C + receiving_ground_station = SVL + receive_time_start = 2023-03-29T07:41:22 + receive_time_end = 2023-03-29T07:42:22 + orbit_start = 22777 + orbit_end = 22778 + actual_product_size = 171868 + state_vector_time = 2023-03-29T05:33:10.999 + semi_major_axis = 7204607302 + eccentricity = 1053 + inclination = 98685 + perigee_argument = 57931 + right_ascension = 149006 + mean_anomaly = 302003 + x_position = -3668934509 + y_position = -6195727325 + z_position = -20970818 + x_velocity = -1426267 + y_velocity = 827418 + z_velocity = 7356929 + earth_sun_distance_ratio = 998227 + location_tolerance_radial = 0 + location_tolerance_crosstrack = 0 + location_tolerance_alongtrack = 0 + yaw_error = 0 + roll_error = 0 + pitch_error = 0 + subsat_latitude_start = -32200 + subsat_longitude_start = 38892 + subsat_latitude_end = -42446 + subsat_longitude_end = 35658 + leap_second = 0 + leap_second_utc = + total_records = 67 + total_mphr = 1 + total_sphr = 1 + total_ipr = 9 + total_geadr = 1 + total_giadr = 0 + total_veadr = 5 + total_viadr = 2 + total_mdr = 48 + count_degraded_inst_mdr = 0 + count_degraded_proc_mdr = 0 + count_degraded_inst_mdr_blocks = 0 + count_degraded_proc_mdr_blocks = 0 + duration_of_product = 176250 + milliseconds_of_data_present = 176250 + milliseconds_of_data_missing = 0 + subsetted_product = F +``` +The variables can be loaded from MetopDataset by indexing the dataset. The variable then works as a lazy array loading the data on indexing: +```julia +ds["latitude"][2,4] +``` +REPL output: +``` +-34.351729 +``` +It is also possible to load the complete array +```julia +ds["latitude"][:,:] +``` +REPL output: +``` +42×48 Matrix{Float64}: + -33.7308 -33.949 … -43.7545 -43.9721 + -33.6969 -33.9152 -43.7252 -43.9429 + -33.6624 -33.8808 -43.695 -43.9127 + -33.6274 -33.8458 -43.6639 -43.8818 + ⋮ ⋱ + -30.1606 -30.3748 -39.9343 -40.1446 + -30.0909 -30.3049 … -39.8538 -40.0638 + -30.0206 -30.2344 -39.7726 -39.9823 +``` +Data from the main product header is accessed as attributes. +```julia +ds.attrib["instrument_id"] +``` +REPL output: +``` +"ASCA" +``` +### Convert a Metop Native binary file to netCDF + +A Metop Native binary file can be converted to netCDF using the [NCDatasets.jl](https://github.com/Alexander-Barth/NCDatasets.jl) package. This +is possible because both MetopDatasets.jl and NCDatasets.jl implement the [CommonDataModel.jl](https://github.com/JuliaGeo/CommonDataModel.jl) interface. + +```julia +using MetopDatasets +using NCDatasets + +input_file = "ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nat" +output_file = "ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nc" + +ds_native = MetopDataset(input_file) +ds_nc = NCDataset(output_file, "c") +NCDatasets.write(ds_nc, ds_native) + +close(ds_native) +close(ds_nc) +``` +It is also possible to use the safe `do` syntax that ensures the files are closed correctly even in the case of exceptions. + +```julia +using MetopDatasets +using NCDatasets + +input_file = "ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nat" +output_file = "ASCA_SZR_1B_M01_20190109125700Z_20190109143858Z_N_O_20190109134816Z.nc" + +MetopDataset(input_file) do ds_native + NCDataset(output_file, "c") do ds_nc + NCDatasets.write(ds_nc, ds_native) + end +end +``` + +## Supported formats +- ASCAT Level 1B +- ASCAT Level 2 Soil Moisture +- IASI Level 1C + +### Formats not yet supported +- AMSU-A Level 1B +- AVHRR Level 1B +- GOME-2 Level 1B +- HIRS Level 1B +- IASI Level 2 Combined Sounding +- MHS Level 1B + +### Reference documents +- [EPS Generic Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/pdf_gen_pfs_13e3f0feb7.pdf) +- [ASCAT Level 1: Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/ASCAT_Level_1_Product_Format_V12_Annexe_50fe72d349.pdf) +- [ASCAT Level 2 Soil Moisture: Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/pdf_ten_0343_eps_ascatl2_pfs_f509981295.pdf) +- [IASI Level 1 Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/pdf_iasi_level_1_pfs_2105bc9ccf.pdf) +- [IASI Level 2 Product Format Specification](https://user.eumetsat.int/s3/eup-strapi-media/pdf_ten_980760_eps_iasi_l2_f9511c26d2.pdf) + +## 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). \ No newline at end of file diff --git a/docs/src/interactive_IASI.png b/docs/src/interactive_IASI.png new file mode 100644 index 0000000..57f0b08 Binary files /dev/null and b/docs/src/interactive_IASI.png differ diff --git a/docs/src/internal_api.md b/docs/src/internal_api.md deleted file mode 100644 index 5ab7a3f..0000000 --- a/docs/src/internal_api.md +++ /dev/null @@ -1,7 +0,0 @@ - -## Internal functions and types -```@autodocs -Modules = [MetopDatasets] -Public = false -Order = [:function, :type, :macro] -``` \ No newline at end of file diff --git a/docs/src/public_api.md b/docs/src/public_api.md deleted file mode 100644 index 5a4d6e2..0000000 --- a/docs/src/public_api.md +++ /dev/null @@ -1,9 +0,0 @@ - - -# MetopDataset - -This is the main interface. - -```@docs -MetopDataset -``` diff --git a/docs/src/python.md b/docs/src/python.md new file mode 100644 index 0000000..fe31bdd --- /dev/null +++ b/docs/src/python.md @@ -0,0 +1,77 @@ +## Use with Python +It should be note that the current focus of MetopDatasets is julia and that python is considered secondary. For python alternatives see +[Eugene](https://anaconda.org/Eumetsat/eugene) ([documentation](https://www-cdn.eumetsat.int/files/2020-04/pdf_ten_02030_ug_eugene.pdf)) or [Satpy](https://satpy.readthedocs.io/en/stable/index.html) (supports many EUMETSAT formats but only limited support for Metop). + +This guide gives a basic example of using MetopDatasets in python via juliacall. For more information see [juliacall documentation](https://juliapy.github.io/PythonCall.jl/stable/juliacall/) for more information. In the future we might make a python wrapper package for MetopDatasets to ease installation and usage. + +## Installation +The installation part just needs to be run once. +### Prerequisites +Julia, Python and Pip all needs to be installed on the machine. This can be checked with the following bash commands. +```bash +python --version +julia --version +pip --version +``` +This guide is tested with the following versions +- Python 3.9.12 +- julia version 1.11.1 +- pip 21.2.4 +### 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. + +```bash +pip install juliacall +pip install numpy +``` + +### Installing MetopDatasets.jl +Install `MetopDatasets.jl` via `juliacall` by running the following Python code. +```python +import juliacall +# make separate module +jl = juliacall.newmodule("MetopDatasetsPy") +jl.seval("import Pkg") +jl.Pkg.add(url="https://github.com/eumetsat/MetopDatasets.jl") +``` + +### Example +You are now ready to use `MetopDatasets.jl` in python. Below are snippets of python code showing a simple example. +#### Loading MetopDatasets in the python session. +```python +import juliacall +import numpy as np +jl = juliacall.newmodule("MetopDatasetsPy") +jl.seval("using MetopDatasets") +``` +#### Reading a dataset +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) +``` +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 +jl.keys(ds) +``` +The individual variables can also be inspected. +```python +ds["gs1cspect"] +``` +The individual variables can be loaded and used like `np.arrays`. The record time is a small variable so we can load it all into memory. +```python +record_start_time = jl.Array(ds["record_start_time"]) +print("record_start_time") +print(np.shape(record_start_time)) +print(np.min(record_start_time)) +print(np.max(record_start_time)) +``` + +It is also possible to just load a slice of a variable. The size of the IASI spectra of an entire orbit is around 2 GB but we can easily load a subset into memory. +```python +spectra_index = 2300 +single_channel_slice = ds["gs1cspect"][spectra_index,:,:,0:10] +print("single_channel_slice") +print(np.shape(single_channel_slice)) +print(np.mean(single_channel_slice)) +``` diff --git a/docs/src/static_IASI.png b/docs/src/static_IASI.png new file mode 100644 index 0000000..165fafe Binary files /dev/null and b/docs/src/static_IASI.png differ diff --git a/docs/src/subset_smr_plot.png b/docs/src/subset_smr_plot.png new file mode 100644 index 0000000..d135bc4 Binary files /dev/null and b/docs/src/subset_smr_plot.png differ diff --git a/docs/src/subset_szr_plot.png b/docs/src/subset_szr_plot.png new file mode 100644 index 0000000..3147a57 Binary files /dev/null and b/docs/src/subset_szr_plot.png differ diff --git a/src/Instruments/IASI/DiskArrayExtra.jl b/src/Instruments/IASI/DiskArrayExtra.jl index 5b443ef..e7ea712 100644 --- a/src/Instruments/IASI/DiskArrayExtra.jl +++ b/src/Instruments/IASI/DiskArrayExtra.jl @@ -54,8 +54,8 @@ end """ IasiWaveNumberDiskArray <: AbstractMetopDiskArray{Float64, 2} -The `IasiWaveNumberDiskArray` is a disk array that computes the wave number of the IASI spectrum. The -wave number is computed using `:idefnsfirst1b` and `:idefspectdwn1b` from each data record. +The `IasiWaveNumberDiskArray` is a disk array that computes the wavenumber of the IASI spectrum. The +wavenumber is computed using `:idefnsfirst1b` and `:idefspectdwn1b` from each data record. """ struct IasiWaveNumberDiskArray <: AbstractMetopDiskArray{Float64, 2} record_type::Type{<:IASI_XXX_1C} @@ -65,12 +65,12 @@ struct IasiWaveNumberDiskArray <: AbstractMetopDiskArray{Float64, 2} size::Tuple{Int64, Int64} end -const IASI_WAVE_NUMBER_NAME = :spectra_wave_number -const IASI_WAVE_NUMBER_DESCRIPTION = "Wave number of IASI 1C spectra samples" +const IASI_WAVENUMBER_NAME = :spectra_wavenumber +const IASI_WAVENUMBER_DESCRIPTION = "Wavenumber of IASI 1C spectra samples" function IasiWaveNumberDiskArray( ds::MetopDataset{R}, field_name::Symbol) where {R <: IASI_XXX_1C} - @assert field_name == IASI_WAVE_NUMBER_NAME + @assert field_name == IASI_WAVENUMBER_NAME number_of_first_sample = MetopDiskArray( ds.file_pointer, ds.data_record_chunks, :idefnsfirst1b) diff --git a/src/Instruments/IASI/IASI.jl b/src/Instruments/IASI/IASI.jl index c9d214d..d3f98f8 100644 --- a/src/Instruments/IASI/IASI.jl +++ b/src/Instruments/IASI/IASI.jl @@ -6,3 +6,4 @@ include("IASI_dimensions.jl") include("IASI_giadr.jl") include("DiskArrayExtra.jl") include("VariableExtra.jl") +include("utils.jl") diff --git a/src/Instruments/IASI/VariableExtra.jl b/src/Instruments/IASI/VariableExtra.jl index db74a29..f1977d0 100644 --- a/src/Instruments/IASI/VariableExtra.jl +++ b/src/Instruments/IASI/VariableExtra.jl @@ -3,11 +3,11 @@ function CDM.varnames(ds::MetopDataset{R}) where {R <: IASI_XXX_1C} if ds.auto_convert - public_names = (string(IASI_WAVE_NUMBER_NAME), _standard_varnames(ds)...) + public_names = (string(IASI_WAVENUMBER_NAME), default_varnames(ds)...) return public_names end - return _standard_varnames(ds) + return default_varnames(ds) end # Overload to enable special scaling for IASI spectrum. @@ -20,27 +20,27 @@ function CDM.variable( N = ndims(sepctrum_disk_array) return MetopVariable{T, N, R}(ds, sepctrum_disk_array) - elseif ds.auto_convert && Symbol(varname) == IASI_WAVE_NUMBER_NAME - wave_number_disk_array = IasiWaveNumberDiskArray(ds, Symbol(varname)) - T = eltype(wave_number_disk_array) - N = ndims(wave_number_disk_array) + elseif ds.auto_convert && Symbol(varname) == IASI_WAVENUMBER_NAME + wavenumber_disk_array = IasiWaveNumberDiskArray(ds, Symbol(varname)) + T = eltype(wavenumber_disk_array) + N = ndims(wavenumber_disk_array) - return MetopVariable{T, N, R}(ds, wave_number_disk_array) + return MetopVariable{T, N, R}(ds, wavenumber_disk_array) end - return _standard_variable(ds, varname) + return default_variable(ds, varname) end # Overload to add :gircimage scale_factor from giadr. Not given in the format specs like the rest. function get_cf_attributes(ds::MetopDataset{R}, field::Symbol, auto_convert::Bool)::Dict{Symbol, Any} where {R <: IASI_XXX_1C} - if field == IASI_WAVE_NUMBER_NAME + if field == IASI_WAVENUMBER_NAME return Dict{Symbol, Any}( :units => "m-1" ) end - cf_attributes = _standard_cf_attributes(R, field, auto_convert) + cf_attributes = default_cf_attributes(R, field, auto_convert) if field == :gircimage giadr = read_first_record(ds, GIADR_IASI_XXX_1C_V11) @@ -53,9 +53,9 @@ end function CDM.attrib( v::MetopVariable{T, N, R}, name::CDM.SymbolOrString) where {T, N, R <: IASI_XXX_1C} - if (v.disk_array.field_name == IASI_WAVE_NUMBER_NAME) && (string(name) == "description") - return IASI_WAVE_NUMBER_DESCRIPTION + if (v.disk_array.field_name == IASI_WAVENUMBER_NAME) && (string(name) == "description") + return IASI_WAVENUMBER_DESCRIPTION end - return _standard_attrib(v, name) + return default_attrib(v, name) end diff --git a/src/Instruments/IASI/utils.jl b/src/Instruments/IASI/utils.jl new file mode 100644 index 0000000..36acf25 --- /dev/null +++ b/src/Instruments/IASI/utils.jl @@ -0,0 +1,32 @@ +# Copyright (c) 2024 EUMETSAT +# License: MIT + +""" + brightness_temperature(I::T, wavenumber::Real, default=T(NaN)) where T <: Real + +Converting the IASI L1C spectrum from radiances to brightness temperature. +Note that the wavenumber must in meters^-1 + +# Example +```julia-repl +julia> file_path = "test/testData/IASI_xxx_1C_M01_20240819103856Z_20240819104152Z_N_C_20240819112911Z" +julia> ds = MetopDataset(file_path); +julia> spectrum = ds["gs1cspect"][:,1,1,1] +julia> wavenumber = ds["spectra_wavenumber"][:, 1] +julia> # convert from radiances to brightness temperature +julia> T_B = brightness_temperature.(spectrum, wavenumber) +``` +""" +function brightness_temperature(I::T, wavenumber::Real, default = T(NaN)) where {T <: Real} + # IASI Level 2: Product Generation Specification + # Equation 26 + c1 = 1.1910427e-16 + c2 = 1.4387752e-2 + + if I > 0 + T_brightness = c2 * wavenumber / log(1 + c1 * wavenumber^3 / I) + return T_brightness + else + return default + end +end diff --git a/src/InterfaceDataModel/Dataset.jl b/src/InterfaceDataModel/Dataset.jl index 4ed948b..9938e51 100644 --- a/src/InterfaceDataModel/Dataset.jl +++ b/src/InterfaceDataModel/Dataset.jl @@ -91,7 +91,7 @@ function MetopDataset( end ## Extend CommonDataModel.AbstractDataset interface -function _standard_varnames(ds::MetopDataset{R}) where {R} +function default_varnames(ds::MetopDataset{R}) where {R} filed_names_no_record_header = string.(fieldnames(R))[2:end] public_fields = ( "record_start_time", "record_stop_time", filed_names_no_record_header...) @@ -99,7 +99,7 @@ function _standard_varnames(ds::MetopDataset{R}) where {R} end function CDM.varnames(ds::MetopDataset{R}) where {R} - return _standard_varnames(ds) + return default_varnames(ds) end # needed for CommonDataModel 0.3.6 and older diff --git a/src/InterfaceDataModel/variable.jl b/src/InterfaceDataModel/variable.jl index deb70ca..751cae6 100644 --- a/src/InterfaceDataModel/variable.jl +++ b/src/InterfaceDataModel/variable.jl @@ -14,10 +14,10 @@ end ### helper functions to get_cf_attributes. function get_cf_attributes(ds::MetopDataset{R}, field::Symbol, auto_convert::Bool)::Dict{Symbol, Any} where {R <: DataRecord} - return _standard_cf_attributes(R, field, auto_convert) # logic is factored out for reusability + return default_cf_attributes(R, field, auto_convert) # logic is factored out for reusability end -function _standard_cf_attributes( +function default_cf_attributes( R::Type{<:DataRecord}, field::Symbol, auto_convert::Bool)::Dict{Symbol, Any} cf_attributes = Dict{Symbol, Any}() @@ -55,10 +55,10 @@ end ## Extend CommonDataModel.AbstractVariable interface function CDM.variable(ds::MetopDataset, varname::CDM.SymbolOrString) - return _standard_variable(ds, varname) # logic is factored out for reusability + return default_variable(ds, varname) # logic is factored out for reusability end -function _standard_variable(ds::MetopDataset{R}, varname::CDM.SymbolOrString) where {R} +function default_variable(ds::MetopDataset{R}, varname::CDM.SymbolOrString) where {R} disk_array = MetopDiskArray(ds.file_pointer, ds.data_record_chunks, Symbol(varname); auto_convert = ds.auto_convert) T = eltype(disk_array) @@ -97,7 +97,7 @@ function CDM.attribnames(v::MetopVariable) return ("description", string.(cf_attributes)...) end -function _standard_attrib(v::MetopVariable, name::CDM.SymbolOrString) +function default_attrib(v::MetopVariable, name::CDM.SymbolOrString) if !(string(name) in CDM.attribnames(v)) error("$name not found") end @@ -112,7 +112,7 @@ function _standard_attrib(v::MetopVariable, name::CDM.SymbolOrString) end function CDM.attrib(v::MetopVariable, name::CDM.SymbolOrString) - return _standard_attrib(v, name) + return default_attrib(v, name) end Base.size(v::MetopVariable) = size(v.disk_array) diff --git a/src/MetopDatasets.jl b/src/MetopDatasets.jl index d84f2b9..c24d17c 100644 --- a/src/MetopDatasets.jl +++ b/src/MetopDatasets.jl @@ -26,7 +26,13 @@ const RECORD_DIM_NAME = "atrack" export MetopDataset -# public functions -@compat public read_first_record, scale_iasi_spectrum, max_giadr_channel +# helper functions +export read_first_record, scale_iasi_spectrum, max_giadr_channel, brightness_temperature + +# Function and types needed to extend the interface +@compat public record_struct_expression, data_record_type +@compat public get_cf_attributes, default_cf_attributes, default_variable +@compat public AbstractMetopDiskArray, MetopDiskArray, MetopVariable +@compat public MainProductHeader, RecordChunk, DataRecord end diff --git a/src/genericTypes/main_product_header.jl b/src/genericTypes/main_product_header.jl index d3f5cd4..267d42b 100644 --- a/src/genericTypes/main_product_header.jl +++ b/src/genericTypes/main_product_header.jl @@ -113,7 +113,7 @@ function native_read(io::IO, T::Type{MainProductHeader})::MainProductHeader main_header_content = String(ntoh.(main_header_content)) # Check for corrupted CRLF header - if occursin('\r',main_header_content) + if occursin('\r', main_header_content) msg = raw"Corrupted file with \r in main product header." msg *= "\n This issue often occurs when transferring native binary files via FTP using ASCII mode." msg *= "\n The native binary files should always be transferred in binary mode." diff --git a/test/ASCAT.jl b/test/ASCAT.jl index b096ba8..611f64c 100644 --- a/test/ASCAT.jl +++ b/test/ASCAT.jl @@ -106,11 +106,11 @@ if isdir("testData") #TODO test data should be handle as an artifact or somethin ds = MetopDataset(SZF_with_dummy_in_mid) # check the number of records - total_count = parse(Int,ds.attrib["total_mdr"]) + total_count = parse(Int, ds.attrib["total_mdr"]) data_count = ds.dim[MetopDatasets.RECORD_DIM_NAME] dummy_count = 3 # the product have 3 dummy records @test total_count == - (data_count + dummy_count) + (data_count + dummy_count) # check that longitude and latitude are in the correct range longitude = Array(ds["longitude_full"]) @@ -120,15 +120,17 @@ if isdir("testData") #TODO test data should be handle as an artifact or somethin @test all((-90 .<= latitude) .& (latitude .<= 90)) # test time stamps - @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_localisation"][1]) < Second(2) - @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_localisation"][end]) < Second(2) + @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_localisation"][1]) < + Second(2) + @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_localisation"][end]) < + Second(2) end @testset "SZR no dummy record" begin ds = MetopDataset(SZR_V13_test_file) # check the number of records - total_count = parse(Int,ds.attrib["total_mdr"]) + total_count = parse(Int, ds.attrib["total_mdr"]) data_count = ds.dim[MetopDatasets.RECORD_DIM_NAME] @test total_count == data_count @@ -140,39 +142,49 @@ if isdir("testData") #TODO test data should be handle as an artifact or somethin @test all((-90 .<= latitude) .& (latitude .<= 90)) # test time stamps - @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_line_nodes"][1]) < Second(2) - @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_line_nodes"][end]) < Second(2) + @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_line_nodes"][1]) < + Second(2) + @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_line_nodes"][end]) < + Second(2) end @testset "SZO " begin - ds= MetopDataset(SZO_V13_test_file) + ds = MetopDataset(SZO_V13_test_file) # check times to sample start and end of file. - @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_line_nodes"][1]) < Second(2) - @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_line_nodes"][end]) < Second(2) + @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_line_nodes"][1]) < + Second(2) + @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_line_nodes"][end]) < + Second(2) end @testset "SZF v11" begin - ds= MetopDataset(SZF_V11_test_file) + ds = MetopDataset(SZF_V11_test_file) # check times to sample start and end of file. - @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_localisation"][1]) < Second(2) - @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_localisation"][end]) < Second(2) + @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_localisation"][1]) < + Second(2) + @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_localisation"][end]) < + Second(2) end @testset "SMR v13" begin - ds= MetopDataset(SMR_V12_test_file) + ds = MetopDataset(SMR_V12_test_file) # check times to sample start and end of file. - @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_line_nodes"][1]) < Second(2) - @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_line_nodes"][end]) < Second(2) + @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_line_nodes"][1]) < + Second(2) + @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_line_nodes"][end]) < + Second(2) end @testset "SMO v13" begin - ds= MetopDataset(SMO_V12_test_file) + ds = MetopDataset(SMO_V12_test_file) # check times to sample start and end of file. - @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_line_nodes"][1]) < Second(2) - @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_line_nodes"][end]) < Second(2) + @test abs(DateTime(ds.attrib["sensing_start"]) - ds["utc_line_nodes"][1]) < + Second(2) + @test abs(DateTime(ds.attrib["sensing_end"]) - ds["utc_line_nodes"][end]) < + Second(2) end end diff --git a/test/IASI.jl b/test/IASI.jl index 3771938..7a91883 100644 --- a/test/IASI.jl +++ b/test/IASI.jl @@ -56,11 +56,11 @@ import CommonDataModel as CDM @test l1c_spectra[92, 1, 1, 1]≈0.0006165 atol=2e-5 @test l1c_spectra[92, 1, 1, 1] isa Float32 - l1c_wave_number = ds["spectra_wave_number"] + l1c_wavenumber = ds["spectra_wavenumber"] # test wavelength of CO2 Q (spectral index 92) should be 15 microns - @test all((1 ./ l1c_wave_number[92, :]) .≈ 14.97566454511419e-6) - @test l1c_wave_number[1, 1] isa Float64 - @test CDM.dimnames(l1c_wave_number) == ["spectral", "atrack"] + @test all((1 ./ l1c_wavenumber[92, :]) .≈ 14.97566454511419e-6) + @test l1c_wavenumber[1, 1] isa Float64 + @test CDM.dimnames(l1c_wavenumber) == ["spectral", "atrack"] # test that channels with no scaling is 0 @test all(l1c_spectra[8462:end, 1:3, 1, 4] .== 0) @@ -98,7 +98,6 @@ import CommonDataModel as CDM @test l1c_spectra[92, 1, 1, 1]≈0.0006165 atol=2e-5 @test l1c_spectra[92, 1, 1, 1] isa Float64 close(ds_high) - end end @@ -122,8 +121,8 @@ end scaled_spectra = MetopDatasets.scale_iasi_spectrum(selected_spectra, giadr) @test scaled_spectra[92, 1, 1]≈0.0006165 atol=2e-5 - # spectra_wave_number is only computed for auto_convert = true - @test !("spectra_wave_number" in keys(ds)) + # spectra_wavenumber is only computed for auto_convert = true + @test !("spectra_wavenumber" in keys(ds)) close(ds) end