Add bbox subsetting for ATL06 data

src/ICESat-2/ATL06.jl

@@ -1,5 +1,5 @@
 """
-    points(g::ICESat2_Granule{:ATL06}, tracks=icesat2_tracks, step=1)
+    points(g::ICESat2_Granule{:ATL06}, tracks=icesat2_tracks, step=1, bbox::Union{Nothing,NamedTuple{}} = nothing)
 
 Retrieve the points for a given ICESat-2 ATL06 (Land Ice) granule as a list of namedtuples, one for each beam.
 The names of the tuples are based on the following fields:
@@ -25,37 +25,79 @@ function points(
     granule::ICESat2_Granule{:ATL06};
     tracks = icesat2_tracks,
     step = 1,
+    bbox::Union{Nothing,NamedTuple{}} = nothing,
 )
     dfs = Vector{NamedTuple}()
     HDF5.h5open(granule.url, "r") do file
         t_offset = read(file, "ancillary_data/atlas_sdp_gps_epoch")[1]::Float64 + gps_offset
         for (i, track) in enumerate(tracks)
             if in(track, keys(file)) && in("land_ice_segments", keys(file[track]))
-                track_nt = points(granule, file, track, t_offset, step)
-                track_nt.height[track_nt.height.==fill_value] .= NaN
+                track_nt = points(granule, file, track, t_offset, step, bbox)
+                if !isempty(track_nt.height)
+                    track_nt.height[track_nt.height.==fill_value] .= NaN
+                end
                 push!(dfs, track_nt)
             end
         end
     end
     return dfs
 end
 
-
 function points(
     ::ICESat2_Granule{:ATL06},
     file::HDF5.H5DataStore,
     track::AbstractString,
     t_offset::Float64,
     step = 1,
+    bbox = bbox::Union{Nothing,NamedTuple{}} = nothing,
 )
-    z = file["$track/land_ice_segments/h_li"][1:step:end]::Vector{Float32}
-    sigma_geo_h = file["$track/land_ice_segments/sigma_geo_h"][1:step:end]::Vector{Float32}
-    h_li_sigma = file["$track/land_ice_segments/h_li_sigma"][1:step:end]::Vector{Float32}
-    x = file["$track/land_ice_segments/longitude"][1:step:end]::Vector{Float64}
-    y = file["$track/land_ice_segments/latitude"][1:step:end]::Vector{Float64}
-    t = file["$track/land_ice_segments/delta_time"][1:step:end]::Vector{Float64}
-    q = file["$track/land_ice_segments/atl06_quality_summary"][1:step:end]::Vector{Int8}
-    dem = file["$track/land_ice_segments/dem/dem_h"][1:step:end]::Vector{Float32}
+    # subset by bbox ?
+    if !isnothing(bbox)
+        x = file["$track/land_ice_segments/longitude"][:]::Vector{Float64}
+        y = file["$track/land_ice_segments/latitude"][:]::Vector{Float64}
+
+        # find index of points inside of bbox
+        ind = (x .> bbox.min_x) .& (y .> bbox.min_y) .& (x .< bbox.max_x) .& (y .< bbox.max_y)
+        start = findfirst(ind)
+        stop = findlast(ind)
+
+        if isnothing(start)
+            @warn "no data found within bbox: $(file.filename)"
+
+            spot_number = attrs(file["$track"])["atlas_spot_number"]::String
+            atlas_beam_type = attrs(file["$track"])["atlas_beam_type"]::String
+
+            nt = (;
+                longitude = Vector{Float64}[],
+                latitude = Vector{Float64}[],
+                height = Vector{Float32}[],
+                height_error = Vector{Float64}[],
+                datetime =Vector{Dates.DateTime}[],
+                quality = BitVector[],
+                track = Fill(track, 0),
+                strong_beam = Fill(atlas_beam_type == "strong", 0),
+                detector_id = Fill(parse(Int8, spot_number), 0),
+                height_reference = Vector{Float32}[],
+            )
+            return nt
+        end
+
+        # only include x and y data within bbox
+        x = x[start:step:stop]
+        y = y[start:step:stop]
+    else
+        start = 1
+        stop = length(file["$track/land_ice_segments/longitude"])
+        x = file["$track/land_ice_segments/longitude"][start:step:stop]::Vector{Float64}
+        y = file["$track/land_ice_segments/latitude"][start:step:stop]::Vector{Float64}
+    end
+
+    z = file["$track/land_ice_segments/h_li"][start:step:stop]::Vector{Float32}
+    sigma_geo_h = file["$track/land_ice_segments/sigma_geo_h"][start:step:stop]::Vector{Float32}
+    h_li_sigma = file["$track/land_ice_segments/h_li_sigma"][start:step:stop]::Vector{Float32}
+    t = file["$track/land_ice_segments/delta_time"][start:step:stop]::Vector{Float64}
+    q = file["$track/land_ice_segments/atl06_quality_summary"][start:step:stop]::Vector{Int8}
+    dem = file["$track/land_ice_segments/dem/dem_h"][start:step:stop]::Vector{Float32}
     spot_number = attrs(file["$track"])["atlas_spot_number"]::String
     atlas_beam_type = attrs(file["$track"])["atlas_beam_type"]::String
     times = unix2datetime.(t .+ t_offset)