Skip to content

Commit

Permalink
Partial ICESat cleanup and add ability to run without internet (#81)
Browse files Browse the repository at this point in the history 
This commit allows asp_plot to run offline by setting `--add_basemap False --plot_icesat False`.

The SlideRule parameters are also refactored and updated to better select "all", "ground", "canopy", and "top_of_canopy" measurements.

Finally, the SlideRule request parameters are now stored with the parquet files, instead of in the filename. If a request is made with different parameters, the existing file is checked and regenerated if the request has different parameters than the original.

A number of TODOs remain in `altimetry.py`, but these should be handled along with remaining updates on issues #62 and #40 on dedicated, smaller pull requests.

---------

Co-authored-by: Ben Purinton <[email protected]>
  • Loading branch information
@dshean @bpurinton
dshean and bpurinton authored Mar 1, 2025
1 parent ac1c131 commit 6610a8f
Show file tree
Hide file tree
Showing 18 changed files with 312 additions and 810 deletions.
16 changes: 13 additions & 3 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,9 +109,9 @@ Options:
logs will be examined to find it. If not
found, no difference plots will be
generated.
--add_basemap BOOLEAN If True, add a contextily basemap to the
figure, which requires internet connection.
Default: True.
--add_basemap BOOLEAN If True, add a basemaps to the figures,
which requires internet connection. Default:
True.
--plot_icesat BOOLEAN If True, plot an ICESat-2 difference plot
with the DEM result. This requires internet
connection to pull ICESat data. Default:
Expand All @@ -124,6 +124,16 @@ Options:
name of ASP processing.
```

### Running without internet connection

If you add these two flags as `False` to the `asp_plot` command, you can run it without internet connection:

```
--add_basemap False --plot_icesat False
```

Otherwise, basemaps will be fetched using contextly and ICESat-2 data will be fetched by SlideRule.

## CLI usage: `csm_camera_plot`

The `csm_camera_plot` command-line tool is a wrapper for outputting a summary plot after running tools like `bundle_adjust` and `jitter_solve`. The inputs must be [CSM camera files](https://stereopipeline.readthedocs.io/en/stable/examples/csm.html). Currently, this tool only supports CSM linescan cameras, such as those from WorldView satellites.
Expand Down
161 changes: 112 additions & 49 deletions asp_plot/altimetry.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
import json
import logging
import os

Expand All @@ -16,8 +17,6 @@
from asp_plot.stereopair_metadata_parser import StereopairMetadataParser
from asp_plot.utils import ColorBar, Raster, glob_file, save_figure

icesat2.init("slideruleearth.io", verbose=True)

logging.basicConfig(level=logging.WARNING)
logger = logging.getLogger(__name__)

Expand Down Expand Up @@ -46,6 +45,9 @@ def __init__(
self.atl06sr_processing_levels = atl06sr_processing_levels
self.atl06sr_processing_levels_filtered = atl06sr_processing_levels_filtered

# Initialize the SlideRule session (requires active connection)
icesat2.init("slideruleearth.io", verbose=True)

# TODO: Implement alongside request_atl03sr below
# if atl03sr is not None and not isinstance(atl03sr, gpd.GeoDataFrame):
# raise ValueError("ATL03 must be a GeoDataFrame if provided.")
Expand Down Expand Up @@ -81,13 +83,13 @@ def __init__(

def request_atl06sr_multi_processing(
self,
processing_levels=["ground", "canopy", "top_of_canopy"],
processing_levels=["all", "ground", "canopy", "top_of_canopy"],
res=20,
len=40,
ats=20,
cnt=10,
maxi=5,
h_sigma_quantile=0.95,
maxi=6,
h_sigma_quantile=1.0,
save_to_parquet=False,
filename="atl06sr",
region=None,
Expand All @@ -98,81 +100,131 @@ def request_atl06sr_multi_processing(
# See parameter discussion on: https://github.com/SlideRuleEarth/sliderule/issues/448
# "srt": -1 tells the server side code to look at the ATL03 confidence array for each photon
# and choose the confidence level that is highest across all five surface type entries.
parms_dict = {
# cnf options: {"atl03_tep", "atl03_not_considered", "atl03_background", "atl03_within_10m", \
# "atl03_low", "atl03_medium", "atl03_high"}
# Note reduced count for limited number of ground photons

# TODO: use the WorldCover values to determine if we should report canopy or top of canopy
# This is tricky, and not clear yet how to go about this. For now, just request all processing levels.

# TODO: Use more generic variable names and strings for functions that are not just limited to atl06
# This can be done when we implement additional requests for other data types

# Shared parameters for all processing levels
shared_parms = {
"poly": region,
"res": res,
"len": len,
"ats": ats,
"maxi": maxi,
"samples": {
"esa_worldcover": {
"asset": "esa-worldcover-10meter",
}
},
}

# Custom parameters for each processing level
custom_parms = {
"all": {
"cnf": "atl03_high",
"srt": -1,
"cnt": cnt,
},
"ground": {
"cnf": 0,
"cnf": "atl03_low",
"srt": -1,
"cnt": 5,
"atl08_class": "atl08_ground",
},
"canopy": {
"cnf": 0,
"cnf": "atl03_medium",
"srt": -1,
"cnt": 5,
"atl08_class": "atl08_canopy",
},
"top_of_canopy": {
"cnf": 0,
"cnf": "atl03_medium",
"srt": -1,
"cnt": 5,
"atl08_class": "atl08_top_of_canopy",
},
}

parms_dict = {
key: parms for key, parms in parms_dict.items() if key in processing_levels
# Filter custom_parms to only include requested processing levels
custom_parms = {
key: parms
for key, parms in custom_parms.items()
if key in processing_levels
}

for key, parms in parms_dict.items():
parms["poly"] = region
parms["res"] = res
parms["len"] = len
parms["ats"] = ats
parms["cnt"] = cnt
parms["maxi"] = maxi
parms["samples"] = {
"esa_worldcover": {
"asset": "esa-worldcover-10meter",
}
}
for key, custom_parm in custom_parms.items():
parms = {**shared_parms, **custom_parm}

fn_base = f"{filename}_res{res}_len{len}_cnt{cnt}_ats{ats}_maxi{maxi}_{key}"
fn_base = f"{filename}_{key}"

print(f"\nICESat-2 ATL06 request processing for: {key}")
fn = f"{fn_base}.parquet"

print(parms)

# Check for existing file with matching parameters
if os.path.exists(fn):
print(f"Existing file found, reading in: {fn}")
atl06sr = gpd.read_parquet(fn)

# Check for parameters in the column
if "sliderule_parameters" in atl06sr.columns:
try:
file_parms = json.loads(atl06sr["sliderule_parameters"].iloc[0])
parms_copy = parms.copy()
parms_copy["poly"] = str(parms_copy["poly"])

if str(parms_copy) != str(file_parms):
print("Parameters don't match request. Regenerating...")
atl06sr = icesat2.atl06p(parms)
if save_to_parquet:
self._save_to_parquet(fn, atl06sr, parms)
except Exception as e:
print(f"Error checking sliderule_parameters column: {e}")
else:
print("No parameters column found, regenerating...")
atl06sr = icesat2.atl06p(parms)
if save_to_parquet:
self._save_to_parquet(fn, atl06sr, parms)
else:
atl06sr = icesat2.atl06p(parms)
if save_to_parquet:
atl06sr.to_parquet(fn)
self._save_to_parquet(fn, atl06sr, parms)

self.atl06sr_processing_levels[key] = atl06sr

print(f"Filtering ATL06-SR {key}")
fn = f"{fn_base}_filtered.parquet"

if os.path.exists(fn):
print(f"Existing file found, reading in: {fn}")
atl06sr_filtered = gpd.read_parquet(fn)
else:
# From Aimee Gibbons:
# I'd recommend anything cycle 03 and later, due to pointing issues before cycle 03.
atl06sr_filtered = atl06sr[atl06sr["cycle"] >= 3]

# Remove bad fits using high percentile of `h_sigma`, the error estimate for the least squares fit model.
# TODO: not sure about h_sigma quantile...might throw out too much. Maybe just remove 0 values?
atl06sr_filtered = atl06sr_filtered[
atl06sr_filtered["h_sigma"]
< atl06sr_filtered["h_sigma"].quantile(h_sigma_quantile)
]
# Also need to filter out 0 values, not sure what these are caused by, but also very bad points.
atl06sr_filtered = atl06sr_filtered[atl06sr_filtered["h_sigma"] != 0]
# From Aimee Gibbons:
# I'd recommend anything cycle 03 and later, due to pointing issues before cycle 03.
atl06sr_filtered = atl06sr[atl06sr["cycle"] >= 3]

if save_to_parquet:
atl06sr_filtered.to_parquet(fn)
# Remove bad fits using high percentile of `h_sigma`, the error estimate for the least squares fit model.
# TODO: not sure about h_sigma quantile...might throw out too much. Maybe just remove 0 values?
atl06sr_filtered = atl06sr_filtered[
atl06sr_filtered["h_sigma"]
< atl06sr_filtered["h_sigma"].quantile(h_sigma_quantile)
]
# Also need to filter out 0 values, not sure what these are caused by, but also very bad points.
atl06sr_filtered = atl06sr_filtered[atl06sr_filtered["h_sigma"] != 0]

self.atl06sr_processing_levels_filtered[key] = atl06sr_filtered

def _save_to_parquet(self, fn, df, parms):
"""Save SlideRule dataframe to parquet including SlideRule parameters"""
# We could save the parameters to the parquet metadata, but this
# was proving rather difficult.
parms_copy = parms.copy()
parms_copy["poly"] = str(parms_copy["poly"])
df["sliderule_parameters"] = json.dumps(parms_copy)
df.to_parquet(fn)

def filter_esa_worldcover(self, filter_out="water", retain_only=None):
# Value Description
# 10 Tree cover
Expand Down Expand Up @@ -225,12 +277,13 @@ def predefined_temporal_filter_atl06sr(self, date=None):

for key in original_keys:
print(
f"\nFiltering ATL06 with 15 day pad, 90 day pad, and seasonal pad around {date} for: {key}"
f"\nFiltering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around {date} for: {key}"
)
atl06sr = self.atl06sr_processing_levels_filtered[key]

fifteen_day = atl06sr[abs(atl06sr.index - date) <= pd.Timedelta(days=15)]
fortyfive_day = atl06sr[abs(atl06sr.index - date) <= pd.Timedelta(days=45)]
ninetyone_day = atl06sr[abs(atl06sr.index - date) <= pd.Timedelta(days=91)]

image_season = date.strftime("%b")
if image_season in ["Dec", "Jan", "Feb"]:
Expand Down Expand Up @@ -258,6 +311,10 @@ def predefined_temporal_filter_atl06sr(self, date=None):
self.atl06sr_processing_levels_filtered[f"{key}_45_day_pad"] = (
fortyfive_day
)
if not ninetyone_day.empty:
self.atl06sr_processing_levels_filtered[f"{key}_91_day_pad"] = (
ninetyone_day
)
if not season_filter.empty:
self.atl06sr_processing_levels_filtered[f"{key}_seasonal"] = (
season_filter
Expand Down Expand Up @@ -391,7 +448,7 @@ def alignment_report(

def plot_atl06sr_time_stamps(
self,
key="ground",
key="all",
title="ICESat-2 ATL06-SR Time Stamps",
cmap="inferno",
map_crs="4326",
Expand Down Expand Up @@ -477,7 +534,7 @@ def plot_atl06sr_time_stamps(

def plot_atl06sr(
self,
key="ground",
key="all",
plot_beams=False,
plot_dem=False,
column_name="h_mean",
Expand All @@ -499,6 +556,7 @@ def plot_atl06sr(

if plot_dem:
ctx_kwargs = {}
# We downsample to speed plotting. This is not carried over into any analysis.
dem_downsampled = gu.Raster(self.dem_fn, downsample=10)
cb = ColorBar(perc_range=(2, 98))
cb.get_clim(dem_downsampled.data)
Expand All @@ -512,6 +570,8 @@ def plot_atl06sr(
)
ax.set_title(None)

# TODO: Implement optional hillshade plotting

if plot_beams:
color_dict = {
1: "red",
Expand Down Expand Up @@ -555,6 +615,9 @@ def plot_atl06sr(
if ctx_kwargs:
ctx.add_basemap(ax=ax, **ctx_kwargs)

ax.set_xticks([])
ax.set_yticks([])

fig.suptitle(f"{title}\n{key} (n={atl06sr.shape[0]})", size=10)
fig.tight_layout()
if save_dir and fig_fn:
Expand Down Expand Up @@ -594,7 +657,7 @@ def atl06sr_to_dem_dh(self):

def mapview_plot_atl06sr_to_dem(
self,
key="ground",
key="all",
clim=None,
plot_aligned=False,
save_dir=None,
Expand Down Expand Up @@ -640,7 +703,7 @@ def mapview_plot_atl06sr_to_dem(

def histogram(
self,
key="ground",
key="all",
title="Histogram",
plot_aligned=False,
save_dir=None,
Expand Down
20 changes: 17 additions & 3 deletions asp_plot/cli/asp_plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@
"--add_basemap",
prompt=False,
default=True,
help="If True, add a contextily basemap to the figure, which requires internet connection. Default: True.",
help="If True, add a basemaps to the figures, which requires internet connection. Default: True.",
)
@click.option(
"--plot_icesat",
Expand Down Expand Up @@ -177,7 +177,7 @@ def main(
)

# Geometry plot
plotter = SceneGeometryPlotter(directory)
plotter = SceneGeometryPlotter(directory, add_basemap=add_basemap)
plotter.dg_geom_plot(
save_dir=plots_directory, fig_fn=f"{next(figure_counter):02}.png"
)
Expand All @@ -187,14 +187,28 @@ def main(
icesat = Altimetry(directory=directory, dem_fn=asp_dem)

icesat.request_atl06sr_multi_processing(
processing_levels=["ground"],
processing_levels=["all", "ground"],
save_to_parquet=False,
)

icesat.filter_esa_worldcover(filter_out="water")

icesat.predefined_temporal_filter_atl06sr()

icesat.mapview_plot_atl06sr_to_dem(
key="all",
save_dir=plots_directory,
fig_fn=f"{next(figure_counter):02}.png",
**ctx_kwargs,
)

icesat.histogram(
key="all",
plot_aligned=False,
save_dir=plots_directory,
fig_fn=f"{next(figure_counter):02}.png",
)

icesat.mapview_plot_atl06sr_to_dem(
key="ground_seasonal",
save_dir=plots_directory,
Expand Down
Loading

0 comments on commit 6610a8f

Please sign in to comment.