Derive individual ATL03 photons from ATL08.

Project.toml

@@ -9,6 +9,7 @@ Aria2_jll = "9ab3bdc3-1250-5043-8fac-ac7e82d2cbc9"
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

src/ICESat-2/ATL08.jl

@@ -1,3 +1,6 @@
+using DataInterpolations: LinearInterpolation
+using Statistics
+
 """
     points(g::ICESat2_Granule{:ATL08}; tracks=icesat2_tracks, step=1, canopy=false, ground=true, bbox::Union{Nothing,Extent,NamedTuple} = nothing)
 
@@ -298,3 +301,228 @@
     )
     nt
 end
+
+function photons(
+    granule::ICESat2_Granule{:ATL08};
+    tracks = icesat2_tracks,
+    step = 1,
+)
+    nts = HDF5.h5open(granule.url, "r") do file
+        t_offset = open_dataset(file, "ancillary_data/atlas_sdp_gps_epoch")[1]::Float64 + gps_offset
+
+        # Determine number of loops over tracks and ground and/or canopy
+        ftracks = filter(track -> haskey(file, track) && haskey(open_group(file, track), "land_segments"), tracks)
+
+        map(ftracks) do track
+            track_nt = _photons(granule, file, track, t_offset)
+            replace!(x -> x === fill_value ? NaN : x, track_nt.height)
+            track_nt
+        end
+    end
+    return PartitionedTable(nts, granule)
+end
+
+round_step(x, step) = round(Int, x / step, RoundDown) * step
+
+function _photons(
+    ::ICESat2_Granule{:ATL08},
+    file::HDF5.H5DataStore,
+    track::AbstractString,
+    t_offset::Float64,
+)
+    group = open_group(file, track)
+    step = 1
+    h_te_best_fit = open_dataset(group, "land_segments/terrain/h_te_best_fit")[:]::Vector{Float32}
+    h_te_best_fit_20m = vec(open_dataset(group, "land_segments/terrain/h_te_best_fit_20m")[:, :]::Array{Float32})
+    #     he = repeat(open_dataset(group, "land_segments/terrain/h_te_uncertainty")[:]::Vector{Float32}, inner = 5)
+    subset_te_flag = open_dataset(group, "land_segments/terrain/subset_te_flag")[:, :]::Matrix{Int8}
+    h_canopy = open_dataset(group, "land_segments/canopy/h_canopy")[:]::Vector{Float32}
+    h_canopy_20m = vec(open_dataset(group, "land_segments/canopy/h_canopy_20m")[:, :]::Array{Float32})
+    #     he = repeat(open_dataset(group, "land_segments/canopy/h_canopy_uncertainty")[1:step:end]::Vector{Float32}, inner = 5)
+    subset_can_flag = open_dataset(group, "land_segments/canopy/subset_can_flag")[:, :]::Matrix{Int8}
+    longitude_20m = vec(open_dataset(group, "land_segments/longitude_20m")[:, 1:step:end]::Array{Float32})
+    longitude = open_dataset(group, "land_segments/longitude")[:]::Vector{Float32}
+    latitude_20m = vec(open_dataset(group, "land_segments/latitude_20m")[:, 1:step:end]::Array{Float32})
+    latitude = open_dataset(group, "land_segments/latitude")[:]::Vector{Float32}
+
+    mask20 = (h_te_best_fit_20m .!= fill_value) .& (longitude_20m .!= fill_value) .& (latitude_20m .!= fill_value)
+    mask = (h_te_best_fit .!= fill_value) .& (longitude .!= fill_value) .| (latitude .!= fill_value)
+
+    delta_time = open_dataset(group, "land_segments/delta_time")[1:step:end]::Vector{Float64}
+    lon_interpol = LinearInterpolation(Float64.(longitude)[mask], delta_time[mask])
+    lat_interpol = LinearInterpolation(Float64.(latitude)[mask], delta_time[mask])
+    time_interpol = LinearInterpolation(delta_time, 1:1:length(delta_time))
+    time20m_interpol = time_interpol(range(1 - 2 // 5, length(delta_time) + 2 // 5, length = 5 * length(delta_time)))
+    height20m_interpol = LinearInterpolation(h_te_best_fit_20m[mask20], time20m_interpol[mask20])
+    # lon20m_interpolation = LinearInterpolation(Float64.(vec(longitude_20m)), time20m_interpol)
+    # lat20m_interpolation = LinearInterpolation(Float64.(vec(latitude_20m)), time20m_interpol)
+
+    # sensitivity = repeat(open_dataset(group, "land_segments/snr")[1:step:end]::Vector{Float32}, inner = 5)
+    # clouds = repeat(open_dataset(group, "land_segments/layer_flag")[1:step:end]::Vector{Int8}, inner = 5)
+    # scattered = repeat(open_dataset(group, "land_segments/msw_flag")[1:step:end]::Vector{Int8}, inner = 5)
+    # saturated = repeat(open_dataset(group, "land_segments/sat_flag")[1:step:end]::Vector{Int8}, inner = 5)
+    # phr = repeat(open_dataset(group, "land_segments/ph_removal_flag")[1:step:end]::Vector{Int8}, inner = 5)
+    # dem = repeat(open_dataset(group, "land_segments/dem_h")[1:step:end]::Vector{Float32}, inner = 5)
+    # times = unix2datetime.(t .+ t_offset)
+    atlas_beam_type = read_attribute(group, "atlas_beam_type")::String
+    spot_number = read_attribute(group, "atlas_spot_number")::String
+
+    # asr = repeat(open_dataset(group, "land_segments/asr")[1:step:end]::Vector{Float32}, inner = 5)
+    # nph = repeat(open_dataset(group, "land_segments/n_seg_ph")[1:step:end]::Vector{Int32}, inner = 5)
+
+    ph_ndx_beg = read(open_dataset(group, "land_segments/ph_ndx_beg"))::Vector{Int64}
+    segment_id_beg = read(open_dataset(group, "land_segments/segment_id_beg"))::Vector{Int32}
+    segment_id_end = read(open_dataset(group, "land_segments/segment_id_end"))::Vector{Int32}
+    # n_seg_ph = read(open_dataset(group, "land_segments/n_seg_ph"))::Vector{Int32}
+
+    @assert issorted(segment_id_beg) "segment_id_beg is not sorted"
+    @assert issorted(segment_id_end) "segment_id_end is not sorted"
+
+    classed_pc_flag = read(group, "signal_photons/classed_pc_flag")::Array{Int8,1}
+    ph_segment_id = read(group, "signal_photons/ph_segment_id")::Array{Int32,1}
+    @assert issorted(ph_segment_id) "ph_segment_id is not sorted"
+
+    # Find the 100m and 20m indices for each individual photon
+    # Nothing lines up in terms of segments. Each 100m segment can
+    # use 1 to 5 20m segments, including jumps between segment numbers.
+    # The photon segment_ids are wider than the segment ids mentioned in the 100m segments
+    # and the number of photons used is less than provided.
+
+    # Find index to 100m segments for attributes
+    index100m = searchsortedfirst.(Ref(segment_id_beg), ph_segment_id) .- 1
+    index100m[insorted.(ph_segment_id, Ref(segment_id_beg))] .+= 1
+    index100me = searchsortedfirst.(Ref(segment_id_end), ph_segment_id)
+    is100m = index100m .== index100me
+    index100m[.!is100m] .= 1  # set to valid index
+
+    @assert index100m[ph_ndx_beg] == 1:length(ph_ndx_beg)
+
+    # 100m can group less than 5 20m segments (!)
+    # @info unique(segment_id_end .- segment_id_beg .+ 1)
+
+    # Find valid segments
+    # Per the ATBD 2.1.18 and 2.2.7, each 20m segment is only valid
+    # if there are at least 10 photons, and 3 classified photons.
+    # Canopy segments therefore also require at least 3 terrain photons
+    # on top of a required 3 canopy photons.
+    canopy_valid = falses(length(ph_segment_id))
+    terrain_valid = falses(length(ph_segment_id))
+    segments = unique(ph_segment_id)
+    segments_nph = zeros(Int16, size(segments))
+    @views for (i, segment) in enumerate(segments)
+        I = searchsorted(ph_segment_id, segment)
+        nt = count(==(1), classed_pc_flag[I])
+        nc = count(>=(2), classed_pc_flag[I])
+        tv = length(I) >= 10 & nt >= 3  # 10 photons, 3 ground
+        cv = tv & nc >= 3  # valid ground and 3 canopy
+        terrain_valid[I] .= tv
+        canopy_valid[I] .= cv
+        segments_nph[i] = nt
+    end
+
+    valid_segment = falses(length(ph_segment_id))
+    @views for (i, class) in enumerate(classed_pc_flag)
+        if class == 1
+            valid_segment[i] = terrain_valid[i]
+        elseif class >= 2
+            valid_segment[i] = canopy_valid[i]
+        end
+    end
+
+    # Assign segment id to 20m photons
+    segments20m = zeros(Int32, size(subset_te_flag))
+    @views for I in eachindex(segment_id_beg)
+        validt = subset_te_flag[:, I]
+        begi = segment_id_beg[I]
+        endi = segment_id_end[I]
+
+        number_of_segments = endi - begi + 1
+        if number_of_segments == 5
+            segments20m[:, I] = begi:endi
+        else
+            i = findfirst(x -> x >= 0, validt)
+            if length(i:5) < number_of_segments
+                # TODO Match by time?
+                @error "Can't place $number_of_segments segments ($begi:$endi) to $validt at $i"
+            else
+                segments20m[i:i+number_of_segments-1, I] = begi:endi
+            end
+        end
+    end
+    segments20mv = vec(segments20m)
+    is20m = insorted.(ph_segment_id, Ref(sort(unique(segments20mv))))
+    prev = 1
+    for i in eachindex(segments20mv)
+        segment = segments20mv[i]
+        if segment == 0
+            segments20mv[i] = prev
+        else
+            prev = segment
+        end
+    end
+    @assert issorted(segments20mv)
+
+    index20m = searchsortedfirst.(Ref(segments20mv), ph_segment_id)
+    index20m[.!is20m] .= 1
+
+    d = read(group, "signal_photons/d_flag")::Array{Int8,1}
+
+    # relative height about what?
+    ph_h = read(group, "signal_photons/ph_h")::Array{Float32,1}
+    ph_h[ph_h.==fill_value] .= NaN
+    delta_time = read(group, "signal_photons/delta_time")::Vector{Float64}
+    # lon = lon20m_interpolation(delta_time)
+    lon = lon_interpol(delta_time)
+    # lat = lat20m_interpolation(delta_time)
+    lat = lat_interpol(delta_time)
+    height = height20m_interpol(delta_time) .+ ph_h
+    times = unix2datetime.(delta_time .+ t_offset)
+    atlas_beam_type = read_attribute(group, "atlas_beam_type")::String
+
+    terrain_height100m = h_te_best_fit[index100m]
+    terrain_height100m[.!is100m] .= NaN
+    terrain_height100m[findall(==(fill_value), terrain_height100m)] .= NaN
+
+    terrain_height20m = h_te_best_fit_20m[index20m]
+    terrain_height20m[.!is20m] .= NaN
+    terrain_height20m[findall(==(fill_value), terrain_height20m)] .= NaN
+
+    canopy_height100m = h_canopy[index100m]
+    canopy_height100m[.!is100m] .= NaN
+    canopy_height100m[findall(==(fill_value), canopy_height100m)] .= NaN
+
+    canopy_height20m = h_canopy_20m[index20m]
+    canopy_height20m[.!is20m] .= NaN
+    canopy_height20m[findall(==(fill_value), canopy_height20m)] .= NaN
+
+    nt = (
+        longitude = lon,
+        latitude = lat,
+        height = height,
+        terrain_height20m = terrain_height20m,
+        terrain_height100m = terrain_height100m,
+        canopy_height20m = canopy_height20m,
+        canopy_height100m = canopy_height100m,
+        # height_error = he,
+        datetime = times,
+        segment_id = ph_segment_id,
+        valid_segment = valid_segment,
+        is20m = is20m,
+        is100m = is100m,
+        # quality = q,
+        # phr = Bool.(phr),
+        # sensitivity = sensitivity,
+        # scattered = scattered,
+        # saturated = saturated,
+        # clouds = Bool.(clouds),
+        track = Fill(track, length(times)),
+        strong_beam = Fill(atlas_beam_type == "strong", length(times)),
+        classification = classed_pc_flag,
+        d_flag = Bool.(d),
+        # height_reference = dem,
+        detector_id = Fill(parse(Int8, spot_number), length(times)),
+        # reflectance = asr,
+        # nphotons = nph,
+    )
+    nt
+end