First support for ICESat GLAH06

.JuliaFormatter.toml

@@ -1,5 +1,3 @@
-always_for_in = true
-always_use_return = true
-margin = 80
-remove_extra_newlines = true
-short_to_long_function_def = true
+format_docstrings=true
+margin=120
+join_lines_based_on_source=true

src/ICESat-2/ATL06.jl

@@ -4,18 +4,18 @@
 Retrieve the points for a given ATL06 (Land Ice) granule as a list of namedtuples, one for each beam.
 The names of the tuples are based on the following fields:
 
-| Column           | Field                                   | Description                                       |
-|------------------|-----------------------------------------|---------------------------------------------------|
-| `longitude`        | `land_ice_segments/longitude`             | Longitude of segment center, WGS84, East=+        |
-| `latitude`         | `land_ice_segments/latitude`              | Latitude of segment center, WGS84, North=+        |
-| `height`           | `land_ice_segments/h_li`                  | Standard land-ice segment height                  |
-| `height_error`     | `land_ice_segments/sigma_geo_h`           | Total vertical geolocation error                  |
-| `datetime`         | `land_ice_segments/delta_time`            | + `ancillary_data/atlas_sdp_gps_epoch` + `gps_offset` |
-| `quality`          | `land_ice_segments/atl06_quality_summary` | Boolean flag indicating the best-quality subset   |
-| `track`            | `gt1l` - `gt3r` groups                    | -                                                 |
-| `power`            | `-`                                       | "strong" or "weak" laser power                    |
-| `detector_id`      | `atlas_spot_number attribute`             | -                                                 |
-| `height_reference` | `land_ice_segments/dem/dem_h`             | Height of the (best available) DEM                |
+| Column           | Field                                   | Description                                       | Units                         |
+|------------------|-----------------------------------------|---------------------------------------------------|-------------------------------|
+| `longitude`        | `land_ice_segments/longitude`             | Longitude of segment center, WGS84, East=+        | decimal degrees              |
+| `latitude`         | `land_ice_segments/latitude`              | Latitude of segment center, WGS84, North=+        | decimal degrees              |
+| `height`           | `land_ice_segments/h_li`                  | Standard land-ice segment height                  | m above the WGS 84 ellipsoid |
+| `height_error`     | `land_ice_segments/sigma_geo_h`           | Total vertical geolocation error                  | m above the WGS 84 ellipsoid |
+| `datetime`         | `land_ice_segments/delta_time`            | + `ancillary_data/atlas_sdp_gps_epoch` + `gps_offset` | date-time                  |
+| `quality`          | `land_ice_segments/atl06_quality_summary` | Boolean flag indicating the best-quality subset   | 1 = high quality             |
+| `track`            | `gt1l` - `gt3r` groups                    | -                                                 | -                            |
+| `power`            | `-`                                       | "strong" or "weak" laser power                    | -                            |
+| `detector_id`      | `atlas_spot_number attribute`             | -                                                 | -                            |
+| `height_reference` | `land_ice_segments/dem/dem_h`             | Height of the (best available) DEM                | -                            |
 
 You can combine the output in a `DataFrame` with `reduce(vcat, DataFrame.(points(g)))`.
 """

src/ICESat/GLAH06.jl

@@ -0,0 +1,76 @@
+"""
+    points(g::ICESat_Granule{:GLAH06})
+
+Retrieve the points for a given GLAH06 (Land Ice) granule as a list of namedtuples
+The names of the tuples are based on the following fields:
+
+| Variable           | Original Field                        | Description                                           | Units                   |
+|:------------------ |:------------------------------------- |:----------------------------------------------------- |:----------------------- |
+| `longitude`        | `Data_40HZ/Geolocation/d_lon`         | Longitude of segment center, WGS84, East=+            | decimal degrees         |
+| `latitude`         | `Data_40HZ/Geolocation/d_lat`         | Latitude of segment center, WGS84, North=+            | decimal degrees         |
+| `height`           | `Data_40HZ/Elevation_Surfaces/d_elev` | + `Data_40HZ/Elevation_Corrections/d_satElevCorr`     | m above WGS84 ellipsoid |
+| `datetime`         | `Data_40HZ/DS_UTCTime_40`             | Precise time of aquisiton                             | date-time               |
+| `quality` [^1]     | `Data_40HZ/Quality/elev_use_flg`      | & `Data_40HZ/Quality/sigma_att_flg` = 0               |                         |
+|                    | & `Data_40HZ/Waveform/i_numPk` = 1    | & `Data_40HZ.Elevation_Corrections/d_satElevCorr` < 3 | 1 = high quality        |
+| `height_reference` | `land_ice_segments/dem/dem_h`         | Height of the (best available) DEM                    | height above WGS84      |
+
+You can get the output in a `DataFrame` with `DataFrame(points(g)`.
+
+[^1]: Smith, B., Fricker, H. A., Gardner, A. S., Medley, B., Nilsson, J., Paolo, F. S., ... & Zwally, H. J. (2020). Pervasive ice sheet mass loss reflects competing ocean and atmosphere processes. Science, 368(6496), 1239-1242.
+"""
+function points(granule::ICESat_Granule{:GLAH06}; step = 1)
+    HDF5.h5open(granule.url, "r") do file
+        height = file["Data_40HZ/Elevation_Surfaces/d_elev"][1:step:end]::Vector{Float64}
+        valid = height .!= icesat_fill
+        height = height[valid]
+
+        saturation_correction =
+            file["Data_40HZ/Elevation_Corrections/d_satElevCorr"][1:step:end][valid]::Vector{Float64}
+        saturation_correction[(saturation_correction.==icesat_fill)] .= 0.0
+        height .+= saturation_correction
+
+        longitude = file["Data_40HZ/Geolocation/d_lon"][1:step:end][valid]::Vector{Float64}
+        longitude[longitude.>180] .= longitude[longitude.>180] .- 360.0  # translate from 0 - 360
+        latitude = file["Data_40HZ/Geolocation/d_lat"][1:step:end][valid]::Vector{Float64}
+
+        datetime = file["Data_40HZ/DS_UTCTime_40"][1:step:end][valid]::Vector{Float64}
+
+        quality = file["Data_40HZ/Quality/elev_use_flg"][1:step:end][valid]::Vector{Int8}
+        sigma_att_flg = file["Data_40HZ/Quality/sigma_att_flg"][1:step:end][valid]::Vector{Int8}
+        i_numPk = file["Data_40HZ/Waveform/i_numPk"][1:step:end][valid]::Vector{Int32}
+
+        height_ref = file["Data_40HZ/Geophysical/d_DEM_elv"][1:step:end][valid]::Vector{Float64}
+        height_ref[height_ref.==icesat_fill] .= NaN
+
+        datetime = unix2datetime.(datetime .+ j2000_offset)
+
+        # convert from TOPEX/POSEIDON to WGS84 ellipsoid using Proj.jl
+        # This pipeline was validated against MATLAB's geodetic2ecef -> ecef2geodetic
+        pipe = Proj.proj_create(
+            "+proj=pipeline +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +inv +proj=longlat +a=6378136.3 +rf=298.257 +e=0.08181922146 +step +proj=cart +a=6378136.3 +rf=298.257 +step +inv +proj=cart +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1",
+        )
+
+        pts = Proj.proj_trans.(pipe, Proj.PJ_FWD, zip(longitude, latitude, height_ref))
+        height_ref = [x[3] for x in pts]::Vector{Float64}
+
+        pts = Proj.proj_trans.(pipe, Proj.PJ_FWD, zip(longitude, latitude, height))
+        height = [x[3] for x in pts]::Vector{Float64}
+
+        # no need to update latitude as differences are well below the precision of the instrument (~1.0e-06 degrees)
+        # latitude = [x[1] for x in pts]::Vector{Float64}
+
+        gt = (
+            longitude = longitude,
+            latitude = latitude,
+            height = height,
+            datetime = datetime,
+            # quality defined according [^1]
+            quality = (quality .== 0) .&
+                      (sigma_att_flg .== 0) .&
+                      (i_numPk .== 1) .&
+                      (saturation_correction .< 3),
+            height_ref = height_ref,
+        )
+        return gt
+    end
+end

src/ICESat/GLAH14.jl

@@ -1,6 +1,4 @@
-const icesat_fill = 1.7976931348623157E308
-
-function points(granule::ICESat_Granule{:GLAH14}; step=1)
+function points(granule::ICESat_Granule{:GLAH14}; step = 1)
     HDF5.h5open(granule.url, "r") do file
 
         zt = file["Data_40HZ/Elevation_Surfaces/d_elev"][1:step:end]::Vector{Float64}
@@ -13,7 +11,7 @@ function points(granule::ICESat_Granule{:GLAH14}; step=1)
 
         x = file["Data_40HZ/Geolocation/d_lon"][1:step:end]::Vector{Float64}
         m .&= (x .!= icesat_fill)
-        x[x .> 180] .= x[x .> 180] .- 360.0  # translate from 0 - 360
+        x[x.>180] .= x[x.>180] .- 360.0  # translate from 0 - 360
         y = file["Data_40HZ/Geolocation/d_lat"][1:step:end]::Vector{Float64}
         m .&= (y .!= icesat_fill)
 
@@ -28,7 +26,7 @@ function points(granule::ICESat_Granule{:GLAH14}; step=1)
         gain_value = file["Data_40HZ/Waveform/i_gval_rcv"][1:step:end][m]::Vector{Int32}
 
         dem = file["Data_40HZ/Geophysical/d_DEM_elv"][1:step:end][m]::Vector{Float64}
-        dem[dem .== icesat_fill] .= NaN
+        dem[dem.==icesat_fill] .= NaN
 
         times = unix2datetime.(t .+ j2000_offset)
 
@@ -48,8 +46,8 @@ function points(granule::ICESat_Granule{:GLAH14}; step=1)
             gain = gain_value,
             reflectivity = ref_flag,
             attitude = att_flag,
-            saturation = sat_corr_flag
-            )
+            saturation = sat_corr_flag,
+        )
         gt
     end
 end

src/ICESat/ICESat.jl

@@ -1,9 +1,8 @@
 using Dates
 
 const j2000_offset = datetime2unix(DateTime(2000, 1, 1, 12, 0, 0))
-const fill_value = 3.4028235f38
 const icesat_inclination = 86.0  # actually 92, so this is 180. - 92.
-
+const icesat_fill = 1.7976931348623157E308
 
 mutable struct ICESat_Granule{product} <: Granule
     id::String
@@ -12,8 +11,8 @@ mutable struct ICESat_Granule{product} <: Granule
 end
 ICESat_Granule(product, args...) = ICESat_Granule{product}(args...)
 
-function Base.copy(g::ICESat_Granule{product}) where product
-    ICESat_Granule(product, g.id, g.url, g.info)
+function Base.copy(g::ICESat_Granule{product}) where {product}
+    return ICESat_Granule(product, g.id, g.url, g.info)
 end
 
 function bounds(granule::ICESat_Granule)
@@ -25,14 +24,12 @@ function bounds(granule::ICESat_Granule)
             min_y = parse(Float64, read(nt["geospatial_lat_min"])),
             max_y = parse(Float64, read(nt["geospatial_lat_max"])),
             min_z = -1000,
-            max_z = 8000
-         )
-        @info ntb
-        ntb
+            max_z = 8000,
+        )
+        return ntb
     end
 end
 
-
 Base.isfile(g::ICESat_Granule) = Base.isfile(g.url)
 
 """Derive info based on file id.
@@ -41,11 +38,21 @@ The id is built up as follows, see 1.2.5 in the user guide
 ATL03_[yyyymmdd][hhmmss]_[ttttccss]_[vvv_rr].h5
 """
 function info(g::ICESat_Granule)
-    icesat_info(g.id)
+    return icesat_info(g.id)
 end
 
 function icesat_info(filename)
     id, _ = splitext(filename)
-    type, revision, orbit, cycle, track, segment, version, filetype = split(id, "_")
-    (type = Symbol(type), phase = parse(Int, orbit[1]), rgt = parse(Int, track[2]), instance = parse(Int, track[3:4]), cycle = parse(Int, cycle), segment = parse(Int, segment), version = parse(Int, version), revision = parse(Int, revision))
+    type, revision, orbit, cycle, track, segment, version, filetype =
+        split(id, "_")
+    return (
+        type = Symbol(type),
+        phase = parse(Int, orbit[1]),
+        rgt = parse(Int, track[2]),
+        instance = parse(Int, track[3:4]),
+        cycle = parse(Int, cycle),
+        segment = parse(Int, segment),
+        version = parse(Int, version),
+        revision = parse(Int, revision),
+    )
 end

src/SpaceLiDAR.jl

@@ -22,6 +22,7 @@ include("ICESat-2/ATL06.jl")
 include("ICESat-2/ATL08.jl")
 include("ICESat-2/ATL12.jl")
 include("ICESat/ICESat.jl")
+include("ICESat/GLAH06.jl")
 include("ICESat/GLAH14.jl")
 include("laz.jl")
 include("interpolate.jl")
@@ -43,6 +44,7 @@ precompile(download!, (ICESat2_Granule,))
 precompile(download!, (GEDI_Granule,))
 precompile(points, (GEDI_Granule,))
 precompile(points, (ICESat2_Granule,))
+precompile(points, (ICESat_Granule,))
 precompile(lines, (GEDI_Granule,))
 precompile(lines, (ICESat2_Granule,))
 precompile(angle, (Vector{Float32}, Vector{Float32}))

src/geoid.jl

@@ -3,7 +3,7 @@ using Downloads
 
 function to_egm2008!(table)
     Proj.enable_network!()
-    trans = Proj.Transformation("EPSG:4326", "EPSG:4326+3855")
+    trans = Proj.Transformation("EPSG:4979", "EPSG:3855")
     data = Proj.Coord.(table.y, table.x, table.z)
     return table[:, :z] .= getindex.(trans.(data), 3)
 end

test/runtests.jl

@@ -23,11 +23,9 @@ end
 download_artifact(v"0.2", "ATL03_20201121151145_08920913_005_01.h5")
 download_artifact(v"0.2", "ATL06_20220404104324_01881512_005_01.h5")
 download_artifact(v"0.2", "ATL08_20201121151145_08920913_005_01.h5")
-download_artifact(
-    v"0.1",
-    "GEDI02_A_2019242104318_O04046_01_T02343_02_003_02_V002.h5",
-)
+download_artifact(v"0.1", "GEDI02_A_2019242104318_O04046_01_T02343_02_003_02_V002.h5")
 download_artifact(v"0.1", "GLAH14_634_1102_001_0071_0_01_0001.H5")
+download_artifact(v"0.1", "GLAH06_634_2131_002_0084_4_01_0001.H5")
 
 @testset "SpaceLiDAR.jl" begin
     @testset "utils" begin
@@ -40,42 +38,24 @@ download_artifact(v"0.1", "GLAH14_634_1102_001_0071_0_01_0001.H5")
     end
 
     @testset "search" begin
-        @test length(find(:ICESat, "GLAH14")) > 0
-        @test length(
-            find(
-                :ICESat2,
-                "ATL03",
-                (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0),
-            ),
-        ) > 0
-        @test length(
-            find(
-                :ICESat2,
-                "ATL08",
-                (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0),
-            ),
-        ) > 0
-        @test length(
-            find(
-                :ICESat2,
-                "ATL06",
-                (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0),
-            ),
-        ) > 0
-        @test length(
-            find(
-                :GEDI,
-                "GEDI02_A",
-                (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0),
-            ),
-        ) > 0
+        @test length(find(:ICESat, "GLAH14", (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0))) > 0
+        @test length(find(:ICESat, "GLAH06", (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0))) > 0
+        @test length(find(:ICESat2, "ATL03", (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0))) > 0
+        @test length(find(:ICESat2, "ATL08", (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0))) > 0
+        @test length(find(:ICESat2, "ATL06", (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0))) > 0
+        @test length(find(:GEDI, "GEDI02_A", (min_x = 4.0, min_y = 40.0, max_x = 5.0, max_y = 50.0))) > 0
     end
 
     @testset "GLAH14" begin
         fn = joinpath(@__DIR__, "data/GLAH14_634_1102_001_0071_0_01_0001.H5")
         g = SpaceLiDAR.granule_from_file(fn)
         points = SpaceLiDAR.points(g)
     end
+    @testset "GLAH06" begin
+        fn = joinpath(@__DIR__, "data/GLAH06_634_2131_002_0084_4_01_0001.H5")
+        g = SpaceLiDAR.granule_from_file(fn)
+        points = SpaceLiDAR.points(g)
+    end
     @testset "ATL03" begin
         fn3 = joinpath(@__DIR__, "data/ATL03_20201121151145_08920913_005_01.h5")
         g3 = SpaceLiDAR.granule_from_file(fn3)