Skip to content

Commit 7e09e73

Browse files
iosefaiosefa
committed
Add tiled EPT scripts for mean intensity and CHM
1 parent 983d593 commit 7e09e73

2 files changed

Lines changed: 281 additions & 0 deletions

File tree

‎scripts/tiled_chm.py‎

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
1+
#!/usr/bin/env python3
2+
"""Create tiled CHM rasters from an EPT dataset using PyForestScan."""
3+
4+
from __future__ import annotations
5+
6+
import argparse
7+
import os
8+
from typing import Tuple
9+
10+
import geopandas as gpd
11+
import numpy as np
12+
from numpy.lib import recfunctions as rfn
13+
14+
from pyforestscan.calculate import calculate_chm
15+
from pyforestscan.handlers import create_geotiff, read_lidar
16+
from pyforestscan.utils import get_bounds_from_ept, get_srs_from_ept
17+
18+
DEFAULT_EPT = "/mnt/x/PROJECTS_2/Big_Island/ChangeHI_Trees/Dry_Forest/Data/Lidar/ept-full/ept.json"
19+
DEFAULT_POLYGON = "/mnt/x/PROJECTS_2/Big_Island/ChangeHI_Trees/Puuwaawaa/Extents/Ahupuaa_Laupahoehoe/laupahoehoe.gpkg"
20+
21+
22+
def ensure_height_above_ground(arr: np.ndarray) -> np.ndarray:
23+
"""Ensure HeightAboveGround exists (required by calculate_chm)."""
24+
if "HeightAboveGround" in arr.dtype.names:
25+
return arr
26+
if "Z" not in arr.dtype.names:
27+
raise ValueError("Point cloud must contain either 'HeightAboveGround' or 'Z'.")
28+
29+
# Fallback for datasets without precomputed HAG.
30+
hag = arr["Z"] - np.nanmin(arr["Z"])
31+
return rfn.append_fields(arr, "HeightAboveGround", hag, usemask=False)
32+
33+
34+
def load_polygon_wkt_and_bounds(polygon_path: str) -> Tuple[str, Tuple[float, float, float, float]]:
35+
"""Load polygon file, dissolve all features, and return WKT + bounds."""
36+
gdf = gpd.read_file(polygon_path)
37+
if gdf.empty:
38+
raise ValueError(f"Polygon file has no features: {polygon_path}")
39+
40+
# Prefer union_all when available (GeoPandas/Shapely newer versions).
41+
if hasattr(gdf.geometry, "union_all"):
42+
geom = gdf.geometry.union_all()
43+
else:
44+
geom = gdf.geometry.unary_union
45+
46+
if geom.is_empty:
47+
raise ValueError(f"Polygon geometry is empty after dissolve: {polygon_path}")
48+
49+
min_x, min_y, max_x, max_y = geom.bounds
50+
return geom.wkt, (min_x, max_x, min_y, max_y)
51+
52+
53+
def parse_args() -> argparse.Namespace:
54+
parser = argparse.ArgumentParser(description=__doc__)
55+
parser.add_argument("--ept", default=DEFAULT_EPT, help="Path/URL to ept.json")
56+
parser.add_argument(
57+
"--polygon",
58+
default=DEFAULT_POLYGON,
59+
help="Polygon file (e.g., .gpkg) used to clip lidar reads",
60+
)
61+
parser.add_argument("--output-dir", default="chm_tiles", help="Output directory for CHM tile GeoTIFFs")
62+
parser.add_argument("--tile-size", type=float, default=1000.0, help="Tile size in map units (x and y)")
63+
parser.add_argument("--xy-res", type=float, default=0.5, help="XY voxel/grid resolution")
64+
parser.add_argument(
65+
"--interpolation",
66+
default="linear",
67+
choices=["nearest", "linear", "cubic", "none"],
68+
help="Interpolation method for CHM gaps",
69+
)
70+
return parser.parse_args()
71+
72+
73+
def main() -> None:
74+
args = parse_args()
75+
os.makedirs(args.output_dir, exist_ok=True)
76+
77+
ept_min_x, ept_max_x, ept_min_y, ept_max_y, min_z, max_z = get_bounds_from_ept(args.ept)
78+
srs = get_srs_from_ept(args.ept)
79+
if not srs:
80+
raise ValueError("Could not determine SRS from EPT metadata.")
81+
82+
polygon_wkt, (poly_min_x, poly_max_x, poly_min_y, poly_max_y) = load_polygon_wkt_and_bounds(args.polygon)
83+
min_x = max(ept_min_x, poly_min_x)
84+
max_x = min(ept_max_x, poly_max_x)
85+
min_y = max(ept_min_y, poly_min_y)
86+
max_y = min(ept_max_y, poly_max_y)
87+
88+
if max_x <= min_x or max_y <= min_y:
89+
raise ValueError("Polygon bounds do not overlap EPT bounds.")
90+
91+
tile_w = tile_h = args.tile_size
92+
num_tiles_x = int(np.ceil((max_x - min_x) / tile_w))
93+
num_tiles_y = int(np.ceil((max_y - min_y) / tile_h))
94+
95+
processed = 0
96+
skipped = 0
97+
interpolation = None if args.interpolation == "none" else args.interpolation
98+
99+
for i in range(num_tiles_x):
100+
for j in range(num_tiles_y):
101+
tile_min_x = min_x + i * tile_w
102+
tile_max_x = min(max_x, tile_min_x + tile_w)
103+
tile_min_y = min_y + j * tile_h
104+
tile_max_y = min(max_y, tile_min_y + tile_h)
105+
106+
bounds = ([tile_min_x, tile_max_x], [tile_min_y, tile_max_y], [min_z, max_z])
107+
arrays = read_lidar(args.ept, srs=srs, bounds=bounds, crop_poly=True, poly=polygon_wkt)
108+
if not arrays or arrays[0].size == 0:
109+
skipped += 1
110+
continue
111+
112+
points = ensure_height_above_ground(arrays[0])
113+
chm, extent = calculate_chm(
114+
points,
115+
voxel_resolution=(args.xy_res, args.xy_res, 1.0),
116+
interpolation=interpolation,
117+
)
118+
119+
out_tif = os.path.join(args.output_dir, f"tile_{i}_{j}_chm.tif")
120+
create_geotiff(chm, out_tif, srs, extent)
121+
processed += 1
122+
123+
print(f"Finished. Wrote {processed} tiles, skipped {skipped} empty tiles.")
124+
125+
126+
if __name__ == "__main__":
127+
main()

‎scripts/tiled_mean_intensity.py‎

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
#!/usr/bin/env python3
2+
"""Create tiled mean-intensity rasters from an EPT dataset using PyForestScan."""
3+
4+
from __future__ import annotations
5+
6+
import argparse
7+
import os
8+
from typing import Tuple
9+
10+
import geopandas as gpd
11+
import numpy as np
12+
from numpy.lib import recfunctions as rfn
13+
from scipy.interpolate import griddata
14+
15+
from pyforestscan.calculate import calculate_voxel_stat
16+
from pyforestscan.handlers import create_geotiff, read_lidar
17+
from pyforestscan.utils import get_bounds_from_ept, get_srs_from_ept
18+
19+
DEFAULT_EPT = "/mnt/x/PROJECTS_2/Big_Island/ChangeHI_Trees/Dry_Forest/Data/Lidar/ept-full/ept.json"
20+
DEFAULT_POLYGON = "/mnt/x/PROJECTS_2/Big_Island/ChangeHI_Trees/Puuwaawaa/Extents/Ahupuaa_Laupahoehoe/laupahoehoe.gpkg"
21+
22+
23+
def ensure_height_above_ground(arr: np.ndarray) -> np.ndarray:
24+
"""Ensure HeightAboveGround exists (required by calculate_voxel_stat)."""
25+
if "HeightAboveGround" in arr.dtype.names:
26+
return arr
27+
if "Z" not in arr.dtype.names:
28+
raise ValueError("Point cloud must contain either 'HeightAboveGround' or 'Z'.")
29+
30+
# Fallback for datasets without precomputed HAG.
31+
hag = arr["Z"] - np.nanmin(arr["Z"])
32+
return rfn.append_fields(arr, "HeightAboveGround", hag, usemask=False)
33+
34+
35+
def interpolate_nans(grid: np.ndarray, method: str = "linear") -> np.ndarray:
36+
"""Fill NaN cells using scipy interpolation over valid neighboring cells."""
37+
out = np.array(grid, dtype=float, copy=True)
38+
nan_mask = np.isnan(out)
39+
if not np.any(nan_mask):
40+
return out
41+
42+
valid_mask = ~nan_mask
43+
if np.count_nonzero(valid_mask) < 3:
44+
return out
45+
46+
x_idx, y_idx = np.indices(out.shape)
47+
points = np.column_stack((x_idx[valid_mask], y_idx[valid_mask]))
48+
values = out[valid_mask]
49+
targets = np.column_stack((x_idx[nan_mask], y_idx[nan_mask]))
50+
51+
interp_vals = griddata(points, values, targets, method=method)
52+
if method != "nearest" and np.any(np.isnan(interp_vals)):
53+
# Fill extrapolation gaps with nearest values.
54+
nearest_vals = griddata(points, values, targets[np.isnan(interp_vals)], method="nearest")
55+
interp_vals[np.isnan(interp_vals)] = nearest_vals
56+
57+
out[nan_mask] = interp_vals
58+
return out
59+
60+
61+
def load_polygon_wkt_and_bounds(polygon_path: str) -> Tuple[str, Tuple[float, float, float, float]]:
62+
"""Load polygon file, dissolve all features, and return WKT + bounds."""
63+
gdf = gpd.read_file(polygon_path)
64+
if gdf.empty:
65+
raise ValueError(f"Polygon file has no features: {polygon_path}")
66+
67+
geom = gdf.geometry.unary_union
68+
if geom.is_empty:
69+
raise ValueError(f"Polygon geometry is empty after dissolve: {polygon_path}")
70+
71+
min_x, min_y, max_x, max_y = geom.bounds
72+
return geom.wkt, (min_x, max_x, min_y, max_y)
73+
74+
75+
def parse_args() -> argparse.Namespace:
76+
parser = argparse.ArgumentParser(description=__doc__)
77+
parser.add_argument("--ept", default=DEFAULT_EPT, help="Path/URL to ept.json")
78+
parser.add_argument(
79+
"--polygon",
80+
default=DEFAULT_POLYGON,
81+
help="Polygon file (e.g., .gpkg) used to clip lidar reads",
82+
)
83+
parser.add_argument("--output-dir", default="mean_intensity_tiles", help="Output directory for GeoTIFF tiles")
84+
parser.add_argument("--tile-size", type=float, default=1000.0, help="Tile size in map units (x and y)")
85+
parser.add_argument("--xy-res", type=float, default=0.5, help="XY voxel/grid resolution")
86+
parser.add_argument("--z-res", type=float, default=1.0, help="Z voxel resolution")
87+
parser.add_argument(
88+
"--interpolation",
89+
default="linear",
90+
choices=["nearest", "linear", "cubic", "none"],
91+
help="Interpolation method for filling NaN cells in output rasters",
92+
)
93+
return parser.parse_args()
94+
95+
96+
def main() -> None:
97+
args = parse_args()
98+
os.makedirs(args.output_dir, exist_ok=True)
99+
100+
ept_min_x, ept_max_x, ept_min_y, ept_max_y, min_z, max_z = get_bounds_from_ept(args.ept)
101+
srs = get_srs_from_ept(args.ept)
102+
if not srs:
103+
raise ValueError("Could not determine SRS from EPT metadata.")
104+
105+
polygon_wkt, (poly_min_x, poly_max_x, poly_min_y, poly_max_y) = load_polygon_wkt_and_bounds(args.polygon)
106+
min_x = max(ept_min_x, poly_min_x)
107+
max_x = min(ept_max_x, poly_max_x)
108+
min_y = max(ept_min_y, poly_min_y)
109+
max_y = min(ept_max_y, poly_max_y)
110+
111+
if max_x <= min_x or max_y <= min_y:
112+
raise ValueError("Polygon bounds do not overlap EPT bounds.")
113+
114+
tile_w = tile_h = args.tile_size
115+
num_tiles_x = int(np.ceil((max_x - min_x) / tile_w))
116+
num_tiles_y = int(np.ceil((max_y - min_y) / tile_h))
117+
118+
processed = 0
119+
skipped = 0
120+
121+
for i in range(num_tiles_x):
122+
for j in range(num_tiles_y):
123+
tile_min_x = min_x + i * tile_w
124+
tile_max_x = min(max_x, tile_min_x + tile_w)
125+
tile_min_y = min_y + j * tile_h
126+
tile_max_y = min(max_y, tile_min_y + tile_h)
127+
128+
bounds = ([tile_min_x, tile_max_x], [tile_min_y, tile_max_y], [min_z, max_z])
129+
arrays = read_lidar(args.ept, srs=srs, bounds=bounds, crop_poly=True, poly=polygon_wkt)
130+
if not arrays or arrays[0].size == 0:
131+
skipped += 1
132+
continue
133+
134+
points = ensure_height_above_ground(arrays[0])
135+
mean_intensity, extent = calculate_voxel_stat(
136+
points,
137+
voxel_resolution=(args.xy_res, args.xy_res, args.z_res),
138+
dimension="Intensity",
139+
stat="mean",
140+
z_index_range=None,
141+
)
142+
143+
if args.interpolation != "none":
144+
mean_intensity = interpolate_nans(mean_intensity, method=args.interpolation)
145+
146+
out_tif = os.path.join(args.output_dir, f"tile_{i}_{j}_mean_intensity.tif")
147+
create_geotiff(mean_intensity, out_tif, srs, extent)
148+
processed += 1
149+
150+
print(f"Finished. Wrote {processed} tiles, skipped {skipped} empty tiles.")
151+
152+
153+
if __name__ == "__main__":
154+
main()

0 commit comments

Comments
 (0)