def extract_patch_data_and_extent(raster_path, center_point, patch_size, pad_value=0):
with rasterio.open(raster_path) as ds:
transform = ds.transform
width, height = ds.width, ds.height
# Convert center point to pixel coordinates
target_r, target_c = ds.index(center_point.x, center_point.y)
# Calculate window
window_r = target_r - patch_size // 2
window_c = target_c - patch_size // 2
window = rasterio.windows.Window(window_c, window_r, patch_size, patch_size)
# Compute actual intersection window
full_window = rasterio.windows.Window(0, 0, width, height)
actual_window = window.intersection(full_window)
# Read data
data = ds.read(window=actual_window)
# Pad data if needed
padded_patch = np.full((ds.count, patch_size, patch_size), pad_value, dtype=ds.dtypes[0])
row_off = max(0, -window_r)
col_off = max(0, -window_c)
padded_patch[:, row_off:row_off + actual_window.height, col_off:col_off + actual_window.width] = data
# Calculate real-world extent (lon/lat)
min_lon, max_lat = transform * (window_c, window_r)
max_lon, min_lat = transform * (window_c + patch_size, window_r + patch_size)
extent = (min_lon, max_lon, min_lat, max_lat)
return padded_patch, extent
def visualize_real_world_patch_data(raster_dict, patch_centers_gdf, boundary_gdf, patch_size, pad_value=0):
"""Visualizes the data of the first patch from each raster, using a dictionary for inputs."""
num_rasters = len(raster_dict)
fig, axes = plt.subplots(1, num_rasters, figsize=(6 * num_rasters, 7), constrained_layout=True)
if num_rasters == 1:
axes = [axes]
# Iterate through the dictionary's items (label, path)
for ax, (label, raster_path) in zip(axes, raster_dict.items()):
resolution = extract_resolution_from_filename(raster_path)
if resolution is None:
resolution = 10 # Default resolution
# Use the first valid center point for demonstration
center_point = patch_centers_gdf.geometry.iloc[0]
lon, lat = center_point.x, center_point.y
# Extract patch data and its geographic extent
patch_data, extent = extract_patch_data_and_extent(
raster_path, center_point, patch_size, pad_value
)
# Plot the patch data using its real-world extent
if patch_data.shape[0] == 1: # Grayscale
ax.imshow(patch_data[0], cmap='viridis', extent=extent, origin='upper', alpha=0.8)
else: # Assumes RGB for multi-band
ax.imshow(np.transpose(patch_data[:3], (1, 2, 0)), extent=extent, origin='upper', alpha=0.8)
# Overlay the boundary for context
boundary_gdf.plot(ax=ax, facecolor='none', edgecolor='lime', linewidth=2, alpha=1.0)
# Plot the center point
ax.scatter(lon, lat, color='red', s=30, zorder=10, label='Patch Center')
# Use the label from the dictionary for a cleaner title
ax.set_title(f"Modality: {label}\nResolution: {resolution}m")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.tick_params(axis='x', rotation=45)
ax.legend()
plt.suptitle("Real-World Spatial Extent of Patch 0 Data", fontsize=18)
plt.show()
def extract_patch_data(raster_path, center_point, patch_size, pad_value=0):
"""Extracts patch data around a geographic center point."""
with rasterio.open(raster_path) as ds:
# Convert center point to pixel coordinates
target_r, target_c = ds.index(center_point.x, center_point.y)
# Calculate window
window_r = target_r - patch_size // 2
window_c = target_c - patch_size // 2
window = rasterio.windows.Window(window_c, window_r, patch_size, patch_size)
# Read data using boundless=True to handle edges
patch_data = ds.read(window=window, boundless=True, fill_value=pad_value)
return patch_data
def visualize_patch_data(raster_files_dict, patch_centers_gdf, patch_size, pad_value=0):
"""Visualizes the data of the first patch from each raster in pixel coordinates (Row/Column)."""
num_rasters = len(raster_files_dict)
fig, axes = plt.subplots(1, num_rasters, figsize=(6 * num_rasters, 7), constrained_layout=True)
if num_rasters == 1:
axes = [axes]
for ax, (label, raster_path) in zip(axes, raster_files_dict.items()):
resolution = extract_resolution_from_filename(raster_path)
if resolution is None:
resolution = 10 # Default resolution
# Use the first valid center point for demonstration
center_point = patch_centers_gdf.geometry.iloc[0]
# Extract patch data
patch_data = extract_patch_data(
raster_path, center_point, patch_size, pad_value
)
# Plot the patch data using its pixel dimensions
if patch_data.shape[0] == 1: # Grayscale
ax.imshow(patch_data[0], cmap='viridis', origin='upper', alpha=0.8)
else: # RGB
ax.imshow(np.transpose(patch_data[:3], (1, 2, 0)), origin='upper', alpha=0.8)
# Plot the center pixel
center_pixel = patch_size // 2
ax.scatter(center_pixel, center_pixel, color='red', s=40, marker='+', zorder=10, label='Patch Center')
ax.set_title(f"{label}\nResolution: {resolution}m")
ax.set_xlabel("Column")
ax.set_ylabel("Row")
ax.legend()
plt.suptitle(f"Patch Data (Pixel Coordinates: {patch_size}x{patch_size})", fontsize=18)
plt.show()
