Skip to content

Commit

Permalink
Reduce allocations and a bit of memory. (#57)
Browse files Browse the repository at this point in the history
  • Loading branch information
@evetion
evetion authored Aug 14, 2023
1 parent 439c2fa commit 97d0d7f
Show file tree
Hide file tree
Showing 8 changed files with 236 additions and 191 deletions.
86 changes: 45 additions & 41 deletions src/GEDI/L2A.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ function points(
nts = Vector{NamedTuple}()
HDF5.h5open(granule.url, "r") do file
for track in tracks
if in(track, keys(file))
if haskey(file, track)
for track_nt ∈ points(granule, file, track, step, bbox, ground, canopy, filtered)
if !isempty(track_nt.height)
track_nt.height[track_nt.height.==fill_value] .= NaN
Expand All @@ -75,20 +75,22 @@ function points(
canopy = false,
filtered = true,
)
group = open_group(file, track)

if !isnothing(bbox)
# find data that falls withing bbox
if ground
x_grd = file["$track/lon_lowestmode"][:]::Vector{Float64}
y_grd = file["$track/lat_lowestmode"][:]::Vector{Float64}
x_grd = read_dataset(group, "lon_lowestmode")::Vector{Float64}
y_grd = read_dataset(group, "lat_lowestmode")::Vector{Float64}

ind = (x_grd .> bbox.X[1]) .& (y_grd .> bbox.Y[1]) .& (x_grd .< bbox.X[2]) .& (y_grd .< bbox.Y[2])
start_grd = findfirst(ind)
stop_grd = findlast(ind)
end

if canopy
x_can = file["$track/lon_highestreturn"][:]::Vector{Float64}
y_can = file["$track/lat_highestreturn"][:]::Vector{Float64}
x_can = read_dataset(group, "lon_highestreturn")::Vector{Float64}
y_can = read_dataset(group, "lat_highestreturn")::Vector{Float64}

ind = (x_can .> bbox.X[1]) .& (y_can .> bbox.Y[1]) .& (x_can .< bbox.X[2]) .& (y_can .< bbox.Y[2])
start_can = findfirst(ind)
Expand All @@ -101,8 +103,8 @@ function points(
start = start_grd
stop = stop_grd
else
start = minimum([start_grd, start_can])
stop = maximum([stop_grd, stop_can])
start = min(start_grd, start_can)
stop = max(stop_grd, stop_can)
end

elseif ground
Expand All @@ -117,7 +119,7 @@ function points(
# no data found
@warn "no data found within bbox of track $track in $(file.filename)"

power = occursin("Full power", read_attribute(file["$track"], "description")::String)
power = occursin("Full power", read_attribute(group, "description")::String)

if canopy
nt_canopy = (
Expand Down Expand Up @@ -188,63 +190,64 @@ function points(
end
else
start = 1
stop = length(file["$track/lon_highestreturn"])
stop = length(open_dataset(group, "lon_highestreturn"))

if ground
x_grd = file["$track/lon_lowestmode"][start:step:stop]::Vector{Float64}
y_grd = file["$track/lat_lowestmode"][start:step:stop]::Vector{Float64}
x_grd = open_dataset(group, "lon_lowestmode")[start:step:stop]::Vector{Float64}
y_grd = open_dataset(group, "lat_lowestmode")[start:step:stop]::Vector{Float64}
end

if canopy
x_can = file["$track/lon_highestreturn"][start:step:stop]::Vector{Float64}
y_can = file["$track/lat_highestreturn"][start:step:stop]::Vector{Float64}
x_can = open_dataset(group, "lon_highestreturn")[start:step:stop]::Vector{Float64}
y_can = open_dataset(group, "lat_highestreturn")[start:step:stop]::Vector{Float64}
end
end

# now that we have the start and stop extents
height_error = file["$track/elevation_bin0_error"][start:step:stop]::Vector{Float32}
height_reference = file["$track/digital_elevation_model"][start:step:stop]::Vector{Float32}
height_error = open_dataset(group, "elevation_bin0_error")[start:step:stop]::Vector{Float32}
height_reference = open_dataset(group, "digital_elevation_model")[start:step:stop]::Vector{Float32}
height_reference[height_reference.==-999999.0] .= NaN
intensity = file["$track/energy_total"][start:step:stop]::Vector{Float32}
aid = file["$track/selected_algorithm"][start:step:stop]::Vector{UInt8}
intensity = open_dataset(group, "energy_total")[start:step:stop]::Vector{Float32}
aid = open_dataset(group, "selected_algorithm")[start:step:stop]::Vector{UInt8}

if canopy
height_can = file["$track/elev_highestreturn"][start:step:stop]::Vector{Float32}
height_can = open_dataset(group, "elev_highestreturn")[start:step:stop]::Vector{Float32}
height_can[(height_can.<-1000.0).&(height_can.>25000.0)] .= NaN
end

if ground
height_grd = file["$track/elev_lowestmode"][start:step:stop]::Vector{Float32}
height_grd = open_dataset(group, "elev_lowestmode")[start:step:stop]::Vector{Float32}
height_grd[(height_grd.<-1000.0).&(height_grd.>25000.0)] .= NaN
end
datetime = file["$track/delta_time"][start:step:stop]::Vector{Float64}
datetime = open_dataset(group, "delta_time")[start:step:stop]::Vector{Float64}

# Quality
quality = file["$track/quality_flag"][start:step:stop]::Vector{UInt8}
rx_assess_quality_flag = file["$track/rx_assess/quality_flag"][start:step:stop]::Vector{UInt8}
degrade_flag = file["$track/degrade_flag"][start:step:stop]::Vector{UInt8}
stale_return_flag = file["$track/geolocation/stale_return_flag"][start:step:stop]::Vector{UInt8}
surface_flag = file["$track/surface_flag"][start:step:stop]::Vector{UInt8}
nmodes = file["$track/num_detectedmodes"][start:step:stop]::Vector{UInt8}

rx_maxamp = file["$track/rx_assess/rx_maxamp"][start:step:stop]::Vector{Float32}
sd_corrected = file["$track/rx_assess/sd_corrected"][start:step:stop]::Vector{Float32}
quality = open_dataset(group, "quality_flag")[start:step:stop]::Vector{UInt8}
rx_assess_quality_flag = open_dataset(group, "rx_assess/quality_flag")[start:step:stop]::Vector{UInt8}
degrade_flag = open_dataset(group, "degrade_flag")[start:step:stop]::Vector{UInt8}
stale_return_flag = open_dataset(group, "geolocation/stale_return_flag")[start:step:stop]::Vector{UInt8}
surface_flag = open_dataset(group, "surface_flag")[start:step:stop]::Vector{UInt8}
nmodes = open_dataset(group, "num_detectedmodes")[start:step:stop]::Vector{UInt8}

rx_maxamp = open_dataset(group, "rx_assess/rx_maxamp")[start:step:stop]::Vector{Float32}
sd_corrected = open_dataset(group, "rx_assess/sd_corrected")[start:step:stop]::Vector{Float32}
rx_maxamp_f = rx_maxamp ./ sd_corrected

# Algorithm
zcross = similar(aid, Float32)
toploc = similar(aid, Float32)
algrun = similar(aid, Bool)
for algorithm = 1:6
zcross[aid.==algorithm] = file["$track/rx_processing_a$algorithm/zcross"][start:step:stop][aid.==algorithm]
toploc[aid.==algorithm] = file["$track/rx_processing_a$algorithm/toploc"][start:step:stop][aid.==algorithm]
algrun[aid.==algorithm] =
file["$track/rx_processing_a$algorithm/rx_algrunflag"][start:step:stop][aid.==algorithm]
am = aid .== algorithm
zcross[am] = open_dataset(group, "rx_processing_a$algorithm/zcross")[start:step:stop][am]
toploc[am] = open_dataset(group, "rx_processing_a$algorithm/toploc")[start:step:stop][am]
algrun[am] =
open_dataset(group, "rx_processing_a$algorithm/rx_algrunflag")[start:step:stop][am]
end

sensitivity = file["$track/sensitivity"][start:step:stop]::Vector{Float32}
sun_angle = file["$track/solar_elevation"][start:step:stop]::Vector{Float32}
power = occursin("Full power", read_attribute(file["$track"], "description")::String)
sensitivity = open_dataset(group, "sensitivity")[start:step:stop]::Vector{Float32}
sun_angle = open_dataset(group, "solar_elevation")[start:step:stop]::Vector{Float32}
power = occursin("Full power", read_attribute(group, "description")::String)

# Quality flags as defined by [^l3]
m = trues(length(quality))
Expand Down Expand Up @@ -336,7 +339,7 @@ function lines(
nts = Vector{NamedTuple}()
HDF5.h5open(granule.url, "r") do file
for track in tracks
if in(track, keys(file))
if haskey(file, track)
for track_df ∈ points(granule, file, track, step, bbox, ground, canopy, filtered)
line = Line(track_df.longitude, track_df.latitude, Float64.(track_df.height))
nt = (geom = line, track = track, strong_beam = track_df.strong_beam[1], granule = granule.id)
Expand All @@ -354,7 +357,7 @@ end
Return the bounds of the GEDI granule.
!!! warning
This opens the .h5 file to read all tracks, so it is very slow.
"""
function bounds(granule::GEDI_Granule)
Expand All @@ -364,9 +367,10 @@ function bounds(granule::GEDI_Granule)
max_ys = -Inf
HDF5.h5open(granule.url, "r") do file
for track ∈ gedi_tracks
if in(track, keys(file))
min_x, max_x = extrema(file["$track/lon_lowestmode"][:])
min_y, max_y = extrema(file["$track/lat_lowestmode"][:])
if haskey(file, track)
group = open_group(file, track)
min_x, max_x = extrema(read_dataset(group, "lon_lowestmode"))
min_y, max_y = extrema(read_dataset(group, "lat_lowestmode"))
min_xs = min(min_xs, min_x)
min_ys = min(min_ys, min_y)
max_xs = max(max_xs, max_x)
Expand Down
60 changes: 32 additions & 28 deletions src/ICESat-2/ATL03.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,10 +39,10 @@ function points(
end
nts = Vector{NamedTuple}()
HDF5.h5open(granule.url, "r") do file
t_offset = file["ancillary_data/atlas_sdp_gps_epoch"][1]::Float64 + gps_offset
t_offset = open_dataset(file, "ancillary_data/atlas_sdp_gps_epoch")[1]::Float64 + gps_offset

for track ∈ tracks
if in(track, keys(file)) && in("heights", keys(file[track]))
if haskey(file, track) && haskey(open_group(file, track), "heights")
track_nt = points(granule, file, track, t_offset, step, bbox)
if !isempty(track_nt.height)
track_nt.height[track_nt.height.==fill_value] .= NaN
Expand Down Expand Up @@ -70,10 +70,11 @@ function lines(
end
nts = Vector{NamedTuple}()
HDF5.h5open(granule.url, "r") do file
t_offset = read(file, "ancillary_data/atlas_sdp_gps_epoch")[1]::Float64 + gps_offset
t_offset = open_dataset(file, "ancillary_data/atlas_sdp_gps_epoch")[1]::Float64 + gps_offset

for track ∈ tracks
if in(track, keys(file)) && in("heights", keys(file[track]))
if haskey(file, track) && haskey(open_group(file, track), "heights")

track_df = points(granule, file, track, t_offset, step, bbox)
line = Line(track_df.longitude, track_df.latitude, Float64.(track_df.height))
i = div(length(track_df.datetime), 2) + 1
Expand Down Expand Up @@ -101,9 +102,10 @@ function points(
bbox::Union{Nothing,Extent} = nothing,
)

group = open_group(file, track)
if !isnothing(bbox)
x = file["$track/heights/lon_ph"][:]::Vector{Float64}
y = file["$track/heights/lat_ph"][:]::Vector{Float64}
x = read_dataset(group, "heights/lon_ph")::Vector{Float64}
y = read_dataset(group, "heights/lat_ph")::Vector{Float64}

# find index of points inside of bbox
ind = (x .> bbox.X[1]) .& (y .> bbox.Y[1]) .& (x .< bbox.X[2]) .& (y .< bbox.Y[2])
Expand All @@ -112,8 +114,8 @@ function points(

if isnothing(start)
@warn "no data found within bbox of track $track in $(file.filename)"
spot_number = attrs(file["$track"])["atlas_spot_number"]::String
atlas_beam_type = attrs(file["$track"])["atlas_beam_type"]::String
spot_number = read_attribute(group, "atlas_spot_number")::String
atlas_beam_type = read_attribute(group, "atlas_beam_type")::String

nt = (
longitude = Float64[],
Expand All @@ -138,39 +140,39 @@ function points(
y = y[start:step:stop]
else
start = 1
stop = length(file["$track/heights/lon_ph"])
x = file["$track/heights/lon_ph"][start:step:stop]::Vector{Float64}
y = file["$track/heights/lat_ph"][start:step:stop]::Vector{Float64}
stop = length(open_dataset(group, "heights/lon_ph"))
x = open_dataset(group, "heights/lon_ph")[start:step:stop]::Vector{Float64}
y = open_dataset(group, "heights/lat_ph")[start:step:stop]::Vector{Float64}
end

height = file["$track/heights/h_ph"][start:step:stop]::Vector{Float32}
datetime = file["$track/heights/delta_time"][start:step:stop]::Vector{Float64}
height = open_dataset(group, "heights/h_ph")[start:step:stop]::Vector{Float32}
datetime = open_dataset(group, "heights/delta_time")[start:step:stop]::Vector{Float64}

# NOT SURE WHY ONLY THE FIRST CONFIDENCE FLAG WAS CHOSEN.. MIGHT NEED TO REVISIT
signal_confidence = file["$track/heights/signal_conf_ph"][1, start:step:stop]::Vector{Int8}
quality = file["$track/heights/quality_ph"][start:step:stop]::Vector{Int8}
signal_confidence = open_dataset(group, "heights/signal_conf_ph")[1, start:step:stop]::Vector{Int8}
quality = open_dataset(group, "heights/quality_ph")[start:step:stop]::Vector{Int8}

# Mapping between segment and photon
seg_cnt = file["$track/geolocation/segment_ph_cnt"][:]::Vector{Int32}
seg_cnt = read(open_dataset(group, "geolocation/segment_ph_cnt"))::Vector{Int32}
ph_ind = count2index(seg_cnt)
ph_ind = ph_ind[start:step:stop]

# extract data posted at segment frequency and map to photon frequency
segment = file["$track/geolocation/segment_id"][:]::Vector{Int32}
segment = segment[ph_ind]
segment = read(open_dataset(group, "geolocation/segment_id"))[ph_ind]::Vector{Int32}
# segment = segment[ph_ind]

sun_angle = file["$track/geolocation/solar_elevation"][:]::Vector{Float32}
sun_angle = sun_angle[ph_ind]
sun_angle = read(open_dataset(group, "geolocation/solar_elevation"))[ph_ind]::Vector{Float32}
# sun_angle = sun_angle[ph_ind]

uncertainty = file["$track/geolocation/sigma_h"][:]::Vector{Float32}
uncertainty = uncertainty[ph_ind]
uncertainty = read(open_dataset(group, "geolocation/sigma_h"))[ph_ind]::Vector{Float32}
# uncertainty = uncertainty[ph_ind]

height_ref = file["$track/geophys_corr/dem_h"][:]::Vector{Float32}
height_ref = height_ref[ph_ind]
height_ref = read(open_dataset(group, "geophys_corr/dem_h"))[ph_ind]::Vector{Float32}
# height_ref = height_ref[ph_ind]

# extract attributes
spot_number = attrs(file["$track"])["atlas_spot_number"]::String
atlas_beam_type = attrs(file["$track"])["atlas_beam_type"]::String
spot_number = read_attribute(group, "atlas_spot_number")::String
atlas_beam_type = read_attribute(group, "atlas_beam_type")::String

# convert from unix time to julia date time
datetime = unix2datetime.(datetime .+ t_offset)
Expand All @@ -193,6 +195,7 @@ function points(
return nt
end


"""
classify(granule::ICESat2_Granule{:ATL03}, atl08::Union{ICESat2_Granule{:ATL08},Nothing} = nothing, tracks = icesat2_tracks)
Expand All @@ -210,10 +213,10 @@ function classify(

dfs = Vector{NamedTuple}()
HDF5.h5open(granule.url, "r") do file
t_offset = read(file, "ancillary_data/atlas_sdp_gps_epoch")[1]::Float64 + gps_offset
t_offset = open_dataset(file, "ancillary_data/atlas_sdp_gps_epoch")[1]::Float64 + gps_offset

for track ∈ tracks
if in(track, keys(file)) && in("heights", keys(file[track]))
if haskey(file, track) && haskey(open_group(file, track), "heights")
track_df = points(granule, file, track, t_offset)

mapping = atl03_mapping(atl08, track)
Expand All @@ -236,6 +239,7 @@ function classify(
dfs
end


function create_mapping(dfsegment, unique_segments)
index_map = Dict{Int64,Int64}()
for unique_segment in unique_segments
Expand Down
36 changes: 19 additions & 17 deletions src/ICESat-2/ATL06.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,9 +37,10 @@ function points(
end
nts = Vector{NamedTuple}()
HDF5.h5open(granule.url, "r") do file
t_offset = file["ancillary_data/atlas_sdp_gps_epoch"][1]::Float64 + gps_offset
t_offset = open_dataset(file, "ancillary_data/atlas_sdp_gps_epoch")[1]::Float64 + gps_offset
for track ∈ tracks
if in(track, keys(file)) && in("land_ice_segments", keys(file[track]))
if haskey(file, track) && haskey(open_group(file, track), "land_ice_segments")

track_nt = points(granule, file, track, t_offset, step, bbox)
if !isempty(track_nt.height)
track_nt.height[track_nt.height.==fill_value] .= NaN
Expand All @@ -59,11 +60,12 @@ function points(
step = 1,
bbox = bbox::Union{Nothing,Extent} = nothing,
)
group = open_group(file, track)

# subset by bbox ?
if !isnothing(bbox)
x = file["$track/land_ice_segments/longitude"][:]::Vector{Float64}
y = file["$track/land_ice_segments/latitude"][:]::Vector{Float64}
x = read_dataset(group, "land_ice_segments/longitude")::Vector{Float64}
y = read_dataset(group, "land_ice_segments/latitude")::Vector{Float64}

# find index of points inside of bbox
ind = (x .> bbox.X[1]) .& (y .> bbox.Y[1]) .& (x .< bbox.X[2]) .& (y .< bbox.Y[2])
Expand All @@ -73,8 +75,8 @@ function points(
if isnothing(start)
@warn "no data found within bbox of track $track in $(file.filename)"

spot_number = attrs(file["$track"])["atlas_spot_number"]::String
atlas_beam_type = attrs(file["$track"])["atlas_beam_type"]::String
spot_number = read_attribute(group, "atlas_spot_number")::String
atlas_beam_type = read_attribute(group, "atlas_beam_type")::String

nt = (;
longitude = Float64[],
Expand All @@ -96,19 +98,19 @@ function points(
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}
stop = length(open_dataset(group, "land_ice_segments/longitude"))
x = open_dataset(group, "land_ice_segments/longitude")[start:step:stop]::Vector{Float64}
y = open_dataset(group, "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
z = open_dataset(group, "land_ice_segments/h_li")[start:step:stop]::Vector{Float32}
sigma_geo_h = open_dataset(group, "land_ice_segments/sigma_geo_h")[start:step:stop]::Vector{Float32}
h_li_sigma = open_dataset(group, "land_ice_segments/h_li_sigma")[start:step:stop]::Vector{Float32}
t = open_dataset(group, "land_ice_segments/delta_time")[start:step:stop]::Vector{Float64}
q = open_dataset(group, "land_ice_segments/atl06_quality_summary")[start:step:stop]::Vector{Int8}
dem = open_dataset(group, "land_ice_segments/dem/dem_h")[start:step:stop]::Vector{Float32}
spot_number = read_attribute(group, "atlas_spot_number")::String
atlas_beam_type = read_attribute(group, "atlas_beam_type")::String
times = unix2datetime.(t .+ t_offset)

sigma_geo_h[sigma_geo_h.==fill_value] .= NaN
Expand Down
Loading

0 comments on commit 97d0d7f

Please sign in to comment.