Add support for 3D visualization

giswqs commented 5 months ago

NEON data cube

https://github.com/opengeos/HyperCoast/assets/5016453/9ce75341-cc73-4ed7-8bb2-47221878876e

https://github.com/opengeos/HyperCoast/assets/5016453/b9730de9-c4f0-4cd4-ac77-27cadaa12dfe

giswqs commented 5 months ago

Display an RGB image on top of an image cube.

import hypercoast
import numpy as np
import pyvista as pv
import matplotlib.pyplot as plt

url = "https://github.com/opengeos/datasets/releases/download/hypercoast/NEON_D02_SERC_DP3_368000_4306000_reflectance.h5"
filepath = "data/neon.h5"
hypercoast.download_file(url, filepath)

dataset = hypercoast.read_neon(filepath) # Load the dataset as an xarray dataset

da = dataset["reflectance"]  # xarray data array
values = da.to_numpy()
print(values.shape)  # (1000, 1000, 426)

# Select three bands for RGB
band_red = 100  # Replace with the index of your red band
band_green = 60  # Replace with the index of your green band
band_blue = 30  # Replace with the index of your blue band

# Extract the selected bands and replace NaNs with zero
red = np.nan_to_num(values[:, :, band_red])
green = np.nan_to_num(values[:, :, band_green])
blue = np.nan_to_num(values[:, :, band_blue])

# Normalize the bands to [0, 1] for display
def normalize_band(band):
    band_min, band_max = np.min(band), np.max(band)
    return (band - band_min) / (band_max - band_min)

red_normalized = normalize_band(red)
green_normalized = normalize_band(green)
blue_normalized = normalize_band(blue)

# Stack the bands to create an RGB image
rgb_image = np.dstack((red_normalized, green_normalized, blue_normalized))

# Display the RGB image using Matplotlib
plt.imshow(rgb_image)
plt.axis('off')
plt.show()

# Create the spatial reference for the image cube
grid = pv.ImageData()

# Set the grid dimensions: shape because we want to inject our values on the POINT data
grid.dimensions = values.shape

# Edit the spatial reference
grid.origin = (0, 0, 0)  # The bottom left corner of the data set
grid.spacing = (1, 1, 1)  # These are the cell sizes along each axis

# Add the data values to the cell data
grid.point_data["values"] = np.nan_to_num(values).flatten(order="F")  # Flatten the array

# Plot the image cube with the RGB image overlay
p = pv.Plotter()
p.add_mesh(grid, cmap="jet", show_edges=False, clim=(0, 0.5))

x_dim, y_dim = rgb_image.shape[0], rgb_image.shape[1]
z_dim = 1
im = pv.ImageData(dimensions=(x_dim, y_dim, z_dim))

# Add scalar data, you may also need to flatten this
im.point_data["rgb_image"] = rgb_image.reshape(-1, 3)*2

grid_z_max = grid.bounds[5]
im.origin = (0, 0, grid_z_max)
p.add_mesh(im, show_edges=False, show_scalar_bar=False, clim=(0, 0.2), rgb=True)

# Show the plot
p.show()

https://github.com/opengeos/HyperCoast/assets/5016453/760cd110-d24c-4c1f-926f-6c923adfc759

opengeos / HyperCoast

Add support for 3D visualization #45