From 31a0241384f312cdf1358f5df4958a857c3bd44e Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 5 Mar 2026 10:28:25 +0100 Subject: [PATCH 1/3] vizsurfdata: Python scripts for visualizing eCLM surfdata files Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.6 --- vizsurfdata/README.md | 40 + vizsurfdata/compare_surfdata.py | 1182 +++++++++++++++++++++++ vizsurfdata/requirements.txt | 4 + vizsurfdata/vizualize_surfdata.py | 1496 +++++++++++++++++++++++++++++ 4 files changed, 2722 insertions(+) create mode 100644 vizsurfdata/README.md create mode 100644 vizsurfdata/compare_surfdata.py create mode 100644 vizsurfdata/requirements.txt create mode 100644 vizsurfdata/vizualize_surfdata.py diff --git a/vizsurfdata/README.md b/vizsurfdata/README.md new file mode 100644 index 0000000..4a75861 --- /dev/null +++ b/vizsurfdata/README.md @@ -0,0 +1,40 @@ +# Python scripts for visualizing and comparing eCLM surface data files + +## Scripts + +### visualize_surfdata.py + +Visualizes all variables in a CLM5 surface data NetCDF file, generating PDF figures and an HTML report. + +```bash +python visualize_surfdata.py +``` + +**Output:** +- `_figures/` - Directory containing PDF figures +- `.html` - Interactive HTML report with embedded figures + +**Sections:** Site location map, basic parameters, land cover, natural PFTs, crop types, soil profiles, monthly LAI/SAI, urban parameters, emission factors, and summary. + +### compare_surfdata.py + +Compares two surface data files side-by-side, highlighting differences with color-coded visualizations. + +```bash +python compare_surfdata.py [-o output_name] +``` + +**Arguments:** +- `file1.nc`, `file2.nc` - NetCDF files to compare +- `label1`, `label2` - Labels for legends and naming (e.g., "2005" "2006" or "control" "experiment") +- `-o output_name` - Optional custom output name (default: `comparison__vs_`) + +**Output:** +- `_figures/` - Directory containing PDF figures +- `.html` - Interactive HTML comparison report + +**Color coding:** Green = no change, Red = increased, Blue = decreased + +## Requirements + +See [`../requirements.txt`](../requirements.txt). `cartopy` is optional — site location maps are skipped gracefully if it is not installed. diff --git a/vizsurfdata/compare_surfdata.py b/vizsurfdata/compare_surfdata.py new file mode 100644 index 0000000..9646407 --- /dev/null +++ b/vizsurfdata/compare_surfdata.py @@ -0,0 +1,1182 @@ +#!/usr/bin/env python3 +""" +Surface Data Comparison Script for CLM5 Surface Data Files +=========================================================== +This script compares two CLM5 surface data NetCDF files and creates +visualizations highlighting the differences, saving them as PDFs and +generating an HTML report. + +Author: Generated for TSMP-PDAF preprocessing +Usage: python compare_surfdata.py [-o output_name] +""" + +import os +import sys +import argparse +import base64 +from datetime import datetime +import numpy as np +import matplotlib +matplotlib.use('Agg') # Non-interactive backend +import matplotlib.pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages +import matplotlib.patches as mpatches +from netCDF4 import Dataset + +# Import shared constants and utilities from visualize_surfdata +from visualize_surfdata import ( + NATPFT_NAMES, CFT_NAMES, URBAN_TYPES, MONTH_NAMES, SOIL_DEPTHS, + get_scalar_value, get_1d_array, fig_to_base64, + create_site_location_figure +) + + +def extract_differing_parts(name1, name2): + """Extract the parts of two filenames that differ.""" + # Remove common prefix + min_len = min(len(name1), len(name2)) + prefix_end = 0 + for i in range(min_len): + if name1[i] != name2[i]: + break + prefix_end = i + 1 + + # Remove common suffix + suffix_start1 = len(name1) + suffix_start2 = len(name2) + for i in range(1, min_len - prefix_end + 1): + if name1[-i] != name2[-i]: + break + suffix_start1 = len(name1) - i + 1 + suffix_start2 = len(name2) - i + 1 + + diff1 = name1[prefix_end:suffix_start1].strip('_') + diff2 = name2[prefix_end:suffix_start2].strip('_') + + if diff1 and diff2: + return f"{diff1}_vs_{diff2}" + else: + # Fallback to timestamp + return datetime.now().strftime('%Y%m%d_%H%M%S') + + +def compute_difference(val1, val2): + """Compute difference and percent change between two values.""" + diff = val2 - val1 + if val1 != 0: + pct_change = (diff / abs(val1)) * 100 + elif val2 != 0: + pct_change = np.inf if diff > 0 else -np.inf + else: + pct_change = 0 + return diff, pct_change + + +def get_change_color(diff, threshold=0.01): + """Get color based on difference: green=same, red=increased, blue=decreased.""" + if abs(diff) < threshold: + return '#48bb78' # Green - no change + elif diff > 0: + return '#e53e3e' # Red - increased + else: + return '#3182ce' # Blue - decreased + + +def plot_comparison_bars(ax, labels, values1, values2, title, units, label1, label2): + """Create side-by-side bar chart comparing two datasets.""" + x = np.arange(len(labels)) + width = 0.35 + + bars1 = ax.bar(x - width/2, values1, width, label=label1, color='#4299e1', edgecolor='#2d3748', alpha=0.8) + bars2 = ax.bar(x + width/2, values2, width, label=label2, color='#ed8936', edgecolor='#2d3748', alpha=0.8) + + ax.set_xticks(x) + ax.set_xticklabels(labels, rotation=45, ha='right', fontsize=8) + ax.set_ylabel(f'[{units}]') + ax.set_title(title, fontsize=11, fontweight='bold') + ax.legend(fontsize=8) + ax.grid(axis='y', alpha=0.3) + + return bars1, bars2 + + +def plot_difference_profile(ax, depths, values1, values2, title, units, label1, label2): + """Create a soil profile comparison with difference shading.""" + diff = values2 - values1 + + # Plot both profiles + ax.plot(values1, depths, 'o-', color='#4299e1', linewidth=2, markersize=8, label=label1) + ax.plot(values2, depths, 's-', color='#ed8936', linewidth=2, markersize=8, label=label2) + + # Shade the difference + ax.fill_betweenx(depths, values1, values2, alpha=0.3, + color='#48bb78' if np.mean(diff) >= 0 else '#e53e3e') + + ax.invert_yaxis() + ax.set_xlabel(f'{title} [{units}]') + ax.set_ylabel('Depth (m)') + ax.set_title(title, fontsize=11, fontweight='bold') + ax.legend(fontsize=8) + ax.grid(alpha=0.3) + + +def plot_monthly_comparison(ax, data1, data2, pft_idx, title, units, label1, label2, pft_name): + """Plot monthly time series comparison for a specific PFT.""" + months = range(12) + + vals1 = data1[:, pft_idx, 0, 0] + vals2 = data2[:, pft_idx, 0, 0] + + ax.plot(months, vals1, 'o-', color='#4299e1', linewidth=2, markersize=6, label=label1) + ax.plot(months, vals2, 's-', color='#ed8936', linewidth=2, markersize=6, label=label2) + ax.fill_between(months, vals1, vals2, alpha=0.3, color='#a0aec0') + + ax.set_xticks(months) + ax.set_xticklabels(MONTH_NAMES) + ax.set_xlabel('Month') + ax.set_ylabel(f'[{units}]') + ax.set_title(f'{title} - {pft_name}', fontsize=10, fontweight='bold') + ax.legend(fontsize=8) + ax.grid(alpha=0.3) + + +def create_comparison_table(ax, params, title, label1='File 1', label2='File 2'): + """Create a comparison table with color-coded differences.""" + ax.axis('off') + + # Table data: [Parameter, Label1, Label2, Difference, Change%] + col_labels = ['Parameter', label1, label2, 'Difference', 'Change %'] + + table_data = [] + cell_colors = [] + + for param in params: + name, val1, val2, units = param + diff, pct = compute_difference(val1, val2) + + # Format values + if isinstance(val1, float): + if abs(val1) < 0.01 and val1 != 0: + val1_str = f'{val1:.2e}' + val2_str = f'{val2:.2e}' + diff_str = f'{diff:.2e}' + else: + val1_str = f'{val1:.4f}' + val2_str = f'{val2:.4f}' + diff_str = f'{diff:.4f}' + else: + val1_str = str(val1) + val2_str = str(val2) + diff_str = str(diff) + + if np.isinf(pct): + pct_str = 'N/A' + else: + pct_str = f'{pct:+.1f}%' + + table_data.append([f'{name} [{units}]', val1_str, val2_str, diff_str, pct_str]) + + # Color based on difference + color = get_change_color(diff, threshold=0.001) + cell_colors.append(['white', 'white', 'white', color, color]) + + table = ax.table(cellText=table_data, colLabels=col_labels, + loc='center', cellLoc='center', colWidths=[0.35, 0.15, 0.15, 0.15, 0.12]) + + table.auto_set_font_size(False) + table.set_fontsize(9) + table.scale(1.2, 1.5) + + # Style header + for j in range(len(col_labels)): + table[(0, j)].set_facecolor('#667eea') + table[(0, j)].set_text_props(color='white', fontweight='bold') + + # Apply cell colors + for i, row_colors in enumerate(cell_colors): + for j, color in enumerate(row_colors): + if color != 'white': + table[(i + 1, j)].set_facecolor(color) + table[(i + 1, j)].set_text_props(color='white', fontweight='bold') + + ax.set_title(title, fontsize=12, fontweight='bold', pad=20) + + +def create_html_report(figures_data, file1, file2, label1, label2, output_dir, output_name): + """Generate an HTML comparison report with all figures embedded.""" + + html_template = """ + + + + + Surface Data Comparison Report + + + +
+
+

CLM5 Surface Data Comparison Report

+

Side-by-side comparison of land surface parameters

+
+
+ {label1} + {file1} +
+
+ {label2} + {file2} +
+
+
+
+
+ No significant change +
+
+
+ Increased in {label2} +
+
+
+ Decreased in {label2} +
+
+
+ Generated: {timestamp} +
+
+ + + + {sections} + +
+

Generated by CLM5 Surface Data Comparison Tool

+

For discussion of model input parameter changes

+
+
+ +""" + + sections_html = "" + nav_links_html = "" + + for section in figures_data: + section_id = section['id'] + section_title = section['title'] + section_desc = section.get('description', '') + + nav_links_html += f'
  • {section_title}
    • \n' + + figures_html = "" + for fig_data in section['figures']: + pdf_filename = fig_data['pdf_name'] + caption = fig_data['caption'] + img_base64 = fig_data['base64'] + + figures_html += f""" +
    + {caption} +

    {caption}

    + Download PDF +
    + """ + + sections_html += f""" +
    +

    {section_title}

    +

    {section_desc}

    + {figures_html} +
    + """ + + html_content = html_template.format( + file1=os.path.basename(file1), + file2=os.path.basename(file2), + label1=label1, + label2=label2, + timestamp=datetime.now().strftime('%Y-%m-%d %H:%M:%S'), + nav_links=nav_links_html, + sections=sections_html + ) + + output_file = os.path.join(output_dir, f'{output_name}.html') + with open(output_file, 'w') as f: + f.write(html_content) + + return output_file + + +def main(): + """Main function to create comparison visualizations.""" + + parser = argparse.ArgumentParser( + description='Compare two CLM5 surface data NetCDF files', + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=""" +Examples: + python compare_surfdata.py file1.nc file2.nc 2005 2006 + python compare_surfdata.py file1.nc file2.nc control experiment -o my_comparison + """ + ) + parser.add_argument('file1', help='First (reference) NetCDF file') + parser.add_argument('file2', help='Second (comparison) NetCDF file') + parser.add_argument('label1', help='Label for the first file (used in legends and naming)') + parser.add_argument('label2', help='Label for the second file (used in legends and naming)') + parser.add_argument('-o', '--output', dest='output_name', default=None, + help='Output name for HTML and figures directory (default: comparison__vs_)') + + args = parser.parse_args() + + # Validate input files + if not os.path.exists(args.file1): + print(f"Error: File not found: {args.file1}") + sys.exit(1) + if not os.path.exists(args.file2): + print(f"Error: File not found: {args.file2}") + sys.exit(1) + + # Use labels for legends + label1 = args.label1 + label2 = args.label2 + + # Determine output name + if args.output_name: + output_name = args.output_name + else: + # Use labels for default output name + output_name = f'comparison_{label1}_vs_{label2}' + + # Create output directory + output_dir = os.path.dirname(args.file1) + if not output_dir: + output_dir = '.' + pdf_dir = os.path.join(output_dir, f'{output_name}_figures') + os.makedirs(pdf_dir, exist_ok=True) + + print(f"Reading NetCDF files...") + print(f" {label1}: {args.file1}") + print(f" {label2}: {args.file2}") + + nc1 = Dataset(args.file1, 'r') + nc2 = Dataset(args.file2, 'r') + + # Get coordinates for both sites + lon1 = get_scalar_value(nc1.variables['LONGXY']) + lat1 = get_scalar_value(nc1.variables['LATIXY']) + lon2 = get_scalar_value(nc2.variables['LONGXY']) + lat2 = get_scalar_value(nc2.variables['LATIXY']) + + figures_data = [] + + # ========================================================================= + # SECTION 0: Site Locations Map + # ========================================================================= + print("Creating Site Locations Map...") + site_map_figures = [] + fig_map = create_site_location_figure( + [lon1, lon2], [lat1, lat2], labels=[label1, label2]) + if fig_map is not None: + pdf_path = os.path.join(pdf_dir, '00_site_locations.pdf') + fig_map.savefig(pdf_path, bbox_inches='tight') + site_map_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': (f'Political map showing the locations of both sites: ' + f'{label1} ({lon1:.3f}°E, {lat1:.3f}°N) and ' + f'{label2} ({lon2:.3f}°E, {lat2:.3f}°N).'), + 'base64': fig_to_base64(fig_map) + }) + plt.close(fig_map) + + if site_map_figures: + figures_data.append({ + 'id': 'site-locations', + 'title': 'Site Locations', + 'description': 'Political map showing the geographic locations of both sites.', + 'figures': site_map_figures + }) + + # ========================================================================= + # SECTION 1: Scalar Parameters Comparison + # ========================================================================= + print("Creating Section 1: Scalar Parameters Comparison...") + section1_figures = [] + + fig, axes = plt.subplots(2, 2, figsize=(16, 12)) + + # Location and basic parameters + basic_params = [ + ('LONGXY', 'degrees E'), + ('LATIXY', 'degrees N'), + ('AREA', 'km^2'), + ('FMAX', 'unitless'), + ('zbedrock', 'm'), + ('SLOPE', 'degrees'), + ('STD_ELEV', 'm'), + ('gdp', 'unitless'), + ] + + basic_data = [] + for var_name, units in basic_params: + val1 = get_scalar_value(nc1.variables[var_name]) + val2 = get_scalar_value(nc2.variables[var_name]) + basic_data.append((var_name, val1, val2, units)) + + create_comparison_table(axes[0, 0], basic_data, 'Basic Site Parameters', label1, label2) + + # Land cover fractions + land_params = [ + ('PCT_NATVEG', '%'), + ('PCT_CROP', '%'), + ('PCT_LAKE', '%'), + ('PCT_WETLAND', '%'), + ('PCT_GLACIER', '%'), + ('peatf', 'fraction'), + ('LANDFRAC_PFT', 'fraction'), + ] + + land_data = [] + for var_name, units in land_params: + val1 = get_scalar_value(nc1.variables[var_name]) + val2 = get_scalar_value(nc2.variables[var_name]) + land_data.append((var_name, val1, val2, units)) + + # Add total urban + urban1 = np.sum(get_1d_array(nc1.variables['PCT_URBAN'])) + urban2 = np.sum(get_1d_array(nc2.variables['PCT_URBAN'])) + land_data.append(('PCT_URBAN (total)', urban1, urban2, '%')) + + create_comparison_table(axes[0, 1], land_data, 'Land Cover Fractions', label1, label2) + + # Soil parameters + soil_params = [ + ('SOIL_COLOR', 'index'), + ('mxsoil_color', 'count'), + ] + + soil_data = [] + for var_name, units in soil_params: + val1 = get_scalar_value(nc1.variables[var_name]) + val2 = get_scalar_value(nc2.variables[var_name]) + soil_data.append((var_name, val1, val2, units)) + + create_comparison_table(axes[1, 0], soil_data, 'Soil Parameters', label1, label2) + + # Emission factors + ef_params = [ + ('EF1_BTR', 'unitless'), + ('EF1_FET', 'unitless'), + ('EF1_FDT', 'unitless'), + ('EF1_SHR', 'unitless'), + ('EF1_GRS', 'unitless'), + ('EF1_CRP', 'unitless'), + ] + + ef_data = [] + for var_name, units in ef_params: + val1 = get_scalar_value(nc1.variables[var_name]) + val2 = get_scalar_value(nc2.variables[var_name]) + ef_data.append((var_name, val1, val2, units)) + + create_comparison_table(axes[1, 1], ef_data, 'Emission Factors', label1, label2) + + fig.suptitle('Scalar Parameters Comparison', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '01_scalar_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Comparison of scalar parameters: basic site info, land cover fractions, soil parameters, and emission factors.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'scalar-comparison', + 'title': 'Scalar Parameters', + 'description': 'Comparison of single-value parameters. Green indicates no change, red indicates increase in File 2, blue indicates decrease.', + 'figures': section1_figures + }) + + # ========================================================================= + # SECTION 2: Land Cover Comparison + # ========================================================================= + print("Creating Section 2: Land Cover Comparison...") + section2_figures = [] + + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + + # Major land cover comparison + land_labels = ['Natural Veg', 'Cropland', 'Urban', 'Lake', 'Wetland', 'Glacier'] + land_vals1 = [ + get_scalar_value(nc1.variables['PCT_NATVEG']), + get_scalar_value(nc1.variables['PCT_CROP']), + np.sum(get_1d_array(nc1.variables['PCT_URBAN'])), + get_scalar_value(nc1.variables['PCT_LAKE']), + get_scalar_value(nc1.variables['PCT_WETLAND']), + get_scalar_value(nc1.variables['PCT_GLACIER']) + ] + land_vals2 = [ + get_scalar_value(nc2.variables['PCT_NATVEG']), + get_scalar_value(nc2.variables['PCT_CROP']), + np.sum(get_1d_array(nc2.variables['PCT_URBAN'])), + get_scalar_value(nc2.variables['PCT_LAKE']), + get_scalar_value(nc2.variables['PCT_WETLAND']), + get_scalar_value(nc2.variables['PCT_GLACIER']) + ] + + plot_comparison_bars(axes[0], land_labels, land_vals1, land_vals2, + 'Major Land Cover Types', '%', label1, label2) + + # Difference bar chart + diffs = np.array(land_vals2) - np.array(land_vals1) + colors = [get_change_color(d, 0.01) for d in diffs] + axes[1].bar(land_labels, diffs, color=colors, edgecolor='#2d3748') + axes[1].axhline(y=0, color='black', linestyle='-', linewidth=0.5) + axes[1].set_ylabel('Difference (File2 - File1) [%]') + axes[1].set_title('Land Cover Change', fontsize=11, fontweight='bold') + axes[1].tick_params(axis='x', rotation=45) + axes[1].grid(axis='y', alpha=0.3) + + for i, (d, label) in enumerate(zip(diffs, land_labels)): + if abs(d) > 0.01: + axes[1].text(i, d + 0.1 * np.sign(d), f'{d:+.2f}', ha='center', fontsize=9) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '02_land_cover_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section2_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Side-by-side comparison of major land cover types and their differences.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'land-cover-comparison', + 'title': 'Land Cover Comparison', + 'description': 'Comparison of major land cover type fractions between the two files.', + 'figures': section2_figures + }) + + # ========================================================================= + # SECTION 3: Natural PFT Comparison + # ========================================================================= + print("Creating Section 3: Natural PFT Comparison...") + section3_figures = [] + + pct_nat_pft1 = get_1d_array(nc1.variables['PCT_NAT_PFT']) + pct_nat_pft2 = get_1d_array(nc2.variables['PCT_NAT_PFT']) + natpft_labels = [NATPFT_NAMES.get(i, f'PFT {i}')[:25] for i in range(15)] + + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + + # Side-by-side comparison + plot_comparison_bars(axes[0], natpft_labels, pct_nat_pft1, pct_nat_pft2, + 'Natural PFT Distribution', '%', label1, label2) + + # Difference plot + diffs = pct_nat_pft2 - pct_nat_pft1 + colors = [get_change_color(d, 0.1) for d in diffs] + y_pos = range(15) + axes[1].barh(y_pos, diffs, color=colors, edgecolor='#2d3748') + axes[1].axvline(x=0, color='black', linestyle='-', linewidth=0.5) + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(natpft_labels, fontsize=8) + axes[1].set_xlabel('Difference (File2 - File1) [%]') + axes[1].set_title('Natural PFT Change', fontsize=11, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '03_natural_pft_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section3_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Comparison of natural Plant Functional Type distributions.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'natural-pft-comparison', + 'title': 'Natural PFT Comparison', + 'description': 'Comparison of natural Plant Functional Type (PFT) distributions within the natural vegetation landunit.', + 'figures': section3_figures + }) + + # ========================================================================= + # SECTION 4: Crop Functional Types Comparison + # ========================================================================= + print("Creating Section 4: CFT Comparison...") + section4_figures = [] + + pct_cft1 = get_1d_array(nc1.variables['PCT_CFT']) + pct_cft2 = get_1d_array(nc2.variables['PCT_CFT']) + cft_indices = nc1.variables['cft'][:] + + # Find CFTs that are non-zero in either file + nonzero_mask = (pct_cft1 > 0.1) | (pct_cft2 > 0.1) + + if np.any(nonzero_mask): + nonzero_cft1 = pct_cft1[nonzero_mask] + nonzero_cft2 = pct_cft2[nonzero_mask] + nonzero_indices = cft_indices[nonzero_mask] + nonzero_labels = [CFT_NAMES.get(int(i), f'CFT {i}')[:25] for i in nonzero_indices] + + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + + plot_comparison_bars(axes[0], nonzero_labels, nonzero_cft1, nonzero_cft2, + 'Crop Type Distribution', '%', label1, label2) + + # Difference plot + diffs = nonzero_cft2 - nonzero_cft1 + colors = [get_change_color(d, 0.1) for d in diffs] + y_pos = range(len(nonzero_labels)) + axes[1].barh(y_pos, diffs, color=colors, edgecolor='#2d3748') + axes[1].axvline(x=0, color='black', linestyle='-', linewidth=0.5) + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(nonzero_labels, fontsize=8) + axes[1].set_xlabel('Difference (File2 - File1) [%]') + axes[1].set_title('Crop Type Change', fontsize=11, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04_cft_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Comparison of Crop Functional Type distributions (showing only non-zero CFTs).', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'cft-comparison', + 'title': 'Crop Functional Types Comparison', + 'description': 'Comparison of Crop Functional Type (CFT) distributions within the crop landunit.', + 'figures': section4_figures + }) + + # ========================================================================= + # SECTION 5: Soil Properties Comparison + # ========================================================================= + print("Creating Section 5: Soil Properties Comparison...") + section5_figures = [] + + fig, axes = plt.subplots(1, 3, figsize=(15, 6)) + + # Sand + sand1 = get_1d_array(nc1.variables['PCT_SAND']) + sand2 = get_1d_array(nc2.variables['PCT_SAND']) + plot_difference_profile(axes[0], SOIL_DEPTHS, sand1, sand2, 'Sand Content', '%', label1, label2) + + # Clay + clay1 = get_1d_array(nc1.variables['PCT_CLAY']) + clay2 = get_1d_array(nc2.variables['PCT_CLAY']) + plot_difference_profile(axes[1], SOIL_DEPTHS, clay1, clay2, 'Clay Content', '%', label1, label2) + + # Organic + org1 = get_1d_array(nc1.variables['ORGANIC']) + org2 = get_1d_array(nc2.variables['ORGANIC']) + plot_difference_profile(axes[2], SOIL_DEPTHS, org1, org2, 'Organic Matter', 'kg/m3', label1, label2) + + fig.suptitle('Soil Properties Comparison by Depth', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05_soil_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Vertical profiles of soil properties with shaded areas showing differences between files.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Soil difference table + fig, ax = plt.subplots(figsize=(14, 8)) + + soil_layer_data = [] + for i, depth in enumerate(SOIL_DEPTHS): + soil_layer_data.append((f'Sand @ {depth:.2f}m', sand1[i], sand2[i], '%')) + soil_layer_data.append((f'Clay @ {depth:.2f}m', clay1[i], clay2[i], '%')) + soil_layer_data.append((f'Organic @ {depth:.2f}m', org1[i], org2[i], 'kg/m3')) + + # Only show layers with differences + diff_layers = [(n, v1, v2, u) for n, v1, v2, u in soil_layer_data if abs(v2 - v1) > 0.01] + + if diff_layers: + create_comparison_table(ax, diff_layers[:15], 'Soil Properties with Differences', label1, label2) + else: + ax.text(0.5, 0.5, 'No significant differences in soil properties', + ha='center', va='center', transform=ax.transAxes, fontsize=14) + ax.axis('off') + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05b_soil_diff_table.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Table of soil properties showing layers with significant differences.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'soil-comparison', + 'title': 'Soil Properties Comparison', + 'description': 'Comparison of soil texture (sand, clay) and organic matter content across depth layers.', + 'figures': section5_figures + }) + + # ========================================================================= + # SECTION 6: Monthly LAI Comparison + # ========================================================================= + print("Creating Section 6: Monthly LAI Comparison...") + section6_figures = [] + + lai1 = nc1.variables['MONTHLY_LAI'][:] + lai2 = nc2.variables['MONTHLY_LAI'][:] + + # Combine PFT names + all_pft_names = {**NATPFT_NAMES, **CFT_NAMES} + + # Find PFTs with non-zero LAI in either file + active_pfts = [] + for pft in range(79): + if np.any(lai1[:, pft, 0, 0] > 0) or np.any(lai2[:, pft, 0, 0] > 0): + active_pfts.append(pft) + + # Plot comparison for top active PFTs + n_plots = min(6, len(active_pfts)) + if n_plots > 0: + fig, axes = plt.subplots(2, 3, figsize=(15, 10)) + axes = axes.flatten() + + for i, pft_idx in enumerate(active_pfts[:n_plots]): + pft_name = all_pft_names.get(pft_idx, f'PFT {pft_idx}')[:30] + plot_monthly_comparison(axes[i], lai1, lai2, pft_idx, 'LAI', 'm2/m2', label1, label2, pft_name) + + # Hide unused axes + for i in range(n_plots, 6): + axes[i].axis('off') + + fig.suptitle('Monthly LAI Comparison by PFT', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06_lai_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Monthly Leaf Area Index (LAI) comparison for active PFTs. Shaded area shows the difference.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # LAI difference heatmap + fig, ax = plt.subplots(figsize=(14, 8)) + + if len(active_pfts) > 0: + lai_diff_matrix = np.zeros((len(active_pfts[:15]), 12)) + for i, pft in enumerate(active_pfts[:15]): + lai_diff_matrix[i, :] = lai2[:, pft, 0, 0] - lai1[:, pft, 0, 0] + + im = ax.imshow(lai_diff_matrix, cmap='RdBu_r', aspect='auto', vmin=-np.max(np.abs(lai_diff_matrix)), vmax=np.max(np.abs(lai_diff_matrix))) + ax.set_xticks(range(12)) + ax.set_xticklabels(MONTH_NAMES) + ax.set_yticks(range(len(active_pfts[:15]))) + ax.set_yticklabels([all_pft_names.get(p, f'PFT {p}')[:25] for p in active_pfts[:15]], fontsize=8) + ax.set_title('LAI Difference (File2 - File1)', fontsize=12, fontweight='bold') + plt.colorbar(im, ax=ax, label='LAI difference (m2/m2)') + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06b_lai_diff_heatmap.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Heatmap of LAI differences by PFT and month. Blue = decrease, Red = increase in File 2.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'lai-comparison', + 'title': 'Monthly LAI Comparison', + 'description': 'Comparison of monthly Leaf Area Index values for active Plant Functional Types.', + 'figures': section6_figures + }) + + # ========================================================================= + # SECTION 7: Urban Parameters Comparison + # ========================================================================= + print("Creating Section 7: Urban Parameters Comparison...") + section7_figures = [] + + fig, axes = plt.subplots(2, 3, figsize=(15, 10)) + + # Urban percent + pct_urban1 = get_1d_array(nc1.variables['PCT_URBAN']) + pct_urban2 = get_1d_array(nc2.variables['PCT_URBAN']) + plot_comparison_bars(axes[0, 0], URBAN_TYPES, pct_urban1, pct_urban2, + 'Urban Fraction', '%', label1, label2) + + # Canyon HWR + hwr1 = get_1d_array(nc1.variables['CANYON_HWR']) + hwr2 = get_1d_array(nc2.variables['CANYON_HWR']) + plot_comparison_bars(axes[0, 1], URBAN_TYPES, hwr1, hwr2, + 'Canyon H/W Ratio', 'ratio', label1, label2) + + # Roof height + ht1 = get_1d_array(nc1.variables['HT_ROOF']) + ht2 = get_1d_array(nc2.variables['HT_ROOF']) + plot_comparison_bars(axes[0, 2], URBAN_TYPES, ht1, ht2, + 'Roof Height', 'm', label1, label2) + + # Building temperature + tb1 = get_1d_array(nc1.variables['T_BUILDING_MIN']) - 273.15 + tb2 = get_1d_array(nc2.variables['T_BUILDING_MIN']) - 273.15 + plot_comparison_bars(axes[1, 0], URBAN_TYPES, tb1, tb2, + 'Min Building Temp', 'C', label1, label2) + + # Roof fraction + rf1 = get_1d_array(nc1.variables['WTLUNIT_ROOF']) + rf2 = get_1d_array(nc2.variables['WTLUNIT_ROOF']) + plot_comparison_bars(axes[1, 1], URBAN_TYPES, rf1, rf2, + 'Roof Fraction', 'fraction', label1, label2) + + # Pervious road fraction + pr1 = get_1d_array(nc1.variables['WTROAD_PERV']) + pr2 = get_1d_array(nc2.variables['WTROAD_PERV']) + plot_comparison_bars(axes[1, 2], URBAN_TYPES, pr1, pr2, + 'Pervious Road Fraction', 'fraction', label1, label2) + + fig.suptitle('Urban Parameters Comparison', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '07_urban_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section7_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Comparison of urban parameters across three density classes.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'urban-comparison', + 'title': 'Urban Parameters Comparison', + 'description': 'Comparison of urban canyon and building parameters for three density classes.', + 'figures': section7_figures + }) + + # ========================================================================= + # SECTION 8: Summary of All Differences + # ========================================================================= + print("Creating Section 8: Summary of Differences...") + section8_figures = [] + + fig = plt.figure(figsize=(16, 12)) + + # Count differences + all_diffs = [] + + # Scalar differences + scalar_vars = ['LONGXY', 'LATIXY', 'AREA', 'FMAX', 'zbedrock', 'SLOPE', 'STD_ELEV', + 'PCT_NATVEG', 'PCT_CROP', 'PCT_LAKE', 'PCT_WETLAND', 'PCT_GLACIER', + 'SOIL_COLOR', 'peatf', 'gdp', 'LAKEDEPTH'] + + for var in scalar_vars: + try: + v1 = get_scalar_value(nc1.variables[var]) + v2 = get_scalar_value(nc2.variables[var]) + diff = abs(v2 - v1) + if diff > 0.001: + all_diffs.append((var, v1, v2, diff)) + except: + pass + + # Create summary visualization + ax1 = fig.add_subplot(2, 2, 1) + + n_total = len(scalar_vars) + 15 + 64 + 30 + 20 # approximate total variables + n_changed = len(all_diffs) + n_same = n_total - n_changed + + ax1.pie([n_same, n_changed], labels=['Unchanged', 'Changed'], + colors=['#48bb78', '#e53e3e'], autopct='%1.1f%%', startangle=90) + ax1.set_title('Overall Parameter Changes', fontsize=12, fontweight='bold') + + # Top differences table + ax2 = fig.add_subplot(2, 2, 2) + ax2.axis('off') + + if all_diffs: + # Sort by absolute difference + sorted_diffs = sorted(all_diffs, key=lambda x: x[3], reverse=True)[:10] + top_diff_data = [(name, v1, v2, 'various') for name, v1, v2, _ in sorted_diffs] + create_comparison_table(ax2, top_diff_data, 'Top 10 Scalar Differences', label1, label2) + else: + ax2.text(0.5, 0.5, 'No significant scalar differences found', + ha='center', va='center', fontsize=14) + + # Land cover change summary + ax3 = fig.add_subplot(2, 2, 3) + land_changes = np.array(land_vals2) - np.array(land_vals1) + colors = [get_change_color(d, 0.01) for d in land_changes] + bars = ax3.barh(land_labels, land_changes, color=colors, edgecolor='#2d3748') + ax3.axvline(x=0, color='black', linestyle='-', linewidth=0.5) + ax3.set_xlabel('Change in Coverage (%)') + ax3.set_title('Land Cover Changes Summary', fontsize=12, fontweight='bold') + ax3.grid(axis='x', alpha=0.3) + + # PFT change summary + ax4 = fig.add_subplot(2, 2, 4) + pft_diffs = pct_nat_pft2 - pct_nat_pft1 + significant_pft_changes = [(NATPFT_NAMES.get(i, f'PFT {i}')[:20], pft_diffs[i]) + for i in range(15) if abs(pft_diffs[i]) > 0.1] + + if significant_pft_changes: + labels, values = zip(*significant_pft_changes) + colors = [get_change_color(v, 0.1) for v in values] + ax4.barh(range(len(labels)), values, color=colors, edgecolor='#2d3748') + ax4.set_yticks(range(len(labels))) + ax4.set_yticklabels(labels, fontsize=9) + ax4.axvline(x=0, color='black', linestyle='-', linewidth=0.5) + ax4.set_xlabel('Change in PFT Coverage (%)') + ax4.set_title('Significant PFT Changes', fontsize=12, fontweight='bold') + ax4.grid(axis='x', alpha=0.3) + else: + ax4.text(0.5, 0.5, 'No significant PFT changes', + ha='center', va='center', fontsize=14) + ax4.axis('off') + + fig.suptitle('Summary of All Differences', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '08_summary.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section8_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Summary overview of all differences between the two surface data files.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'summary', + 'title': 'Summary of Differences', + 'description': 'A comprehensive summary of all parameter differences for quick reference.', + 'figures': section8_figures + }) + + # Close NetCDF files + nc1.close() + nc2.close() + + # Generate HTML report + print("Generating HTML report...") + html_file = create_html_report(figures_data, args.file1, args.file2, label1, label2, output_dir, output_name) + + print(f"\nComparison complete!") + print(f"PDF figures saved in: {pdf_dir}") + print(f"HTML report saved as: {html_file}") + + return html_file + + +if __name__ == '__main__': + main() diff --git a/vizsurfdata/requirements.txt b/vizsurfdata/requirements.txt new file mode 100644 index 0000000..573ce58 --- /dev/null +++ b/vizsurfdata/requirements.txt @@ -0,0 +1,4 @@ +numpy +matplotlib +netCDF4 +cartopy diff --git a/vizsurfdata/vizualize_surfdata.py b/vizsurfdata/vizualize_surfdata.py new file mode 100644 index 0000000..ff31d28 --- /dev/null +++ b/vizsurfdata/vizualize_surfdata.py @@ -0,0 +1,1496 @@ +#!/usr/bin/env python3 +""" +Surface Data Visualization Script for CLM5 Surface Data Files +============================================================== +This script reads a CLM5 surface data NetCDF file and creates comprehensive +visualizations for all variables, saving them as PDFs and generating an +HTML report for easy sharing and discussion. + +Author: Generated for TSMP-PDAF preprocessing +Usage: python visualize_surfdata.py +""" + +import os +import sys +import base64 +from datetime import datetime +import numpy as np +import matplotlib +matplotlib.use('Agg') # Non-interactive backend +import matplotlib.pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages +import matplotlib.patches as mpatches +from netCDF4 import Dataset + +try: + import cartopy.crs as ccrs + import cartopy.feature as cfeature + HAS_CARTOPY = True +except ImportError: + HAS_CARTOPY = False + print("Warning: cartopy not available. Site location map will be skipped.") + +# PFT/CFT names sourced from the official CLM5 documentation: +# NATPFT (Table 2.2.1): https://escomp.github.io/CTSM/release-clm5.0/tech_note/Ecosystem/CLM50_Tech_Note_Ecosystem.html#id15 +# CFT (Table 2.26.1): https://escomp.github.io/CTSM/release-clm5.0/tech_note/Crop_Irrigation/CLM50_Tech_Note_Crop_Irrigation.html#id20 +# +# At import time we try to parse those tables from the live HTML (stdlib only, 3-second +# timeout). If the network is unavailable the hardcoded fallbacks below are used instead. + +def _fetch_pft_names_from_docs(): + """Parse NATPFT and CFT name tables from the official CLM5 HTML documentation. + + Returns (natpft_dict, cft_dict), where each value is either a populated + {int: str} dict or None when fetching/parsing failed. + """ + import urllib.request + from html.parser import HTMLParser + + _NATPFT_URL = ( + "https://escomp.github.io/CTSM/release-clm5.0/tech_note/" + "Ecosystem/CLM50_Tech_Note_Ecosystem.html" + ) + _CFT_URL = ( + "https://escomp.github.io/CTSM/release-clm5.0/tech_note/" + "Crop_Irrigation/CLM50_Tech_Note_Crop_Irrigation.html" + ) + + class _TableParser(HTMLParser): + """Collect cell text for every row inside a .""" + def __init__(self, table_id): + super().__init__() + self._target = table_id + self._active = False + self._in_cell = False + self._cur_row = [] + self._cur_text = "" + self.rows = [] + + def handle_starttag(self, tag, attrs): + d = dict(attrs) + if tag == "table" and d.get("id") == self._target: + self._active = True + if self._active: + if tag == "tr": + self._cur_row = [] + elif tag in ("td", "th"): + self._in_cell = True + self._cur_text = "" + + def handle_endtag(self, tag): + if not self._active: + return + if tag == "table": + self._active = False + elif tag == "tr": + if self._cur_row: + self.rows.append(self._cur_row) + self._cur_row = [] + elif tag in ("td", "th"): + self._cur_row.append(self._cur_text.strip()) + self._in_cell = False + + def handle_data(self, data): + if self._in_cell: + self._cur_text += data + + def _parse(url, table_id): + try: + with urllib.request.urlopen(url, timeout=3) as resp: + html = resp.read().decode("utf-8", errors="replace") + except Exception: + return None + parser = _TableParser(table_id) + parser.feed(html) + result = {} + for row in parser.rows: + if len(row) < 2: + continue + try: + result[int(row[0])] = row[1] + except ValueError: + pass # skip header rows whose first cell is not an integer + return result if result else None + + natpft = _parse(_NATPFT_URL, table_id="id15") + cft = _parse(_CFT_URL, table_id="id20") + return natpft, cft + + +# Hardcoded fallback — names transcribed from the official docs (see URLs above). +# PFT 9 is "Broadleaf evergreen shrub – temperate" (the "temperate" qualifier matters; +# a tropical broadleaf evergreen shrub is not represented in CLM5). +_NATPFT_FALLBACK = { + 0: "Bare Ground", + 1: "Needleleaf evergreen tree – temperate", + 2: "Needleleaf evergreen tree – boreal", + 3: "Needleleaf deciduous tree – boreal", + 4: "Broadleaf evergreen tree – tropical", + 5: "Broadleaf evergreen tree – temperate", + 6: "Broadleaf deciduous tree – tropical", + 7: "Broadleaf deciduous tree – temperate", + 8: "Broadleaf deciduous tree – boreal", + 9: "Broadleaf evergreen shrub – temperate", + 10: "Broadleaf deciduous shrub – temperate", + 11: "Broadleaf deciduous shrub – boreal", + 12: "C3 arctic grass", + 13: "C3 grass", + 14: "C4 grass", +} + +_CFT_FALLBACK = { + 15: "C3 unmanaged rainfed crop", 16: "C3 unmanaged irrigated crop", + 17: "Temperate corn rainfed", 18: "Temperate corn irrigated", + 19: "Spring wheat rainfed", 20: "Spring wheat irrigated", + 21: "Winter wheat rainfed", 22: "Winter wheat irrigated", + 23: "Temperate soybean rainfed", 24: "Temperate soybean irrigated", + 25: "Barley rainfed", 26: "Barley irrigated", + 27: "Winter barley rainfed", 28: "Winter barley irrigated", + 29: "Rye rainfed", 30: "Rye irrigated", + 31: "Winter rye rainfed", 32: "Winter rye irrigated", + 33: "Cassava rainfed", 34: "Cassava irrigated", + 35: "Citrus rainfed", 36: "Citrus irrigated", + 37: "Cocoa rainfed", 38: "Cocoa irrigated", + 39: "Coffee rainfed", 40: "Coffee irrigated", + 41: "Cotton rainfed", 42: "Cotton irrigated", + 43: "Datepalm rainfed", 44: "Datepalm irrigated", + 45: "Foddergrass rainfed", 46: "Foddergrass irrigated", + 47: "Grapes rainfed", 48: "Grapes irrigated", + 49: "Groundnuts rainfed", 50: "Groundnuts irrigated", + 51: "Millet rainfed", 52: "Millet irrigated", + 53: "Oilpalm rainfed", 54: "Oilpalm irrigated", + 55: "Potatoes rainfed", 56: "Potatoes irrigated", + 57: "Pulses rainfed", 58: "Pulses irrigated", + 59: "Rapeseed rainfed", 60: "Rapeseed irrigated", + 61: "Rice rainfed", 62: "Rice irrigated", + 63: "Sorghum rainfed", 64: "Sorghum irrigated", + 65: "Sugarbeet rainfed", 66: "Sugarbeet irrigated", + 67: "Sugarcane rainfed", 68: "Sugarcane irrigated", + 69: "Sunflower rainfed", 70: "Sunflower irrigated", + 71: "Miscanthus rainfed", 72: "Miscanthus irrigated", + 73: "Switchgrass rainfed", 74: "Switchgrass irrigated", + 75: "Tropical corn rainfed", 76: "Tropical corn irrigated", + 77: "Tropical soybean rainfed", 78: "Tropical soybean irrigated", +} + +_fetched_natpft, _fetched_cft = _fetch_pft_names_from_docs() +NATPFT_NAMES = _fetched_natpft if _fetched_natpft is not None else _NATPFT_FALLBACK +CFT_NAMES = _fetched_cft if _fetched_cft is not None else _CFT_FALLBACK + +URBAN_TYPES = ["Tall Building District (TBD)", "High Density (HD)", "Medium Density (MD)"] + +MONTH_NAMES = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', + 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] + +# Soil layer depths (approximate centers in meters) +SOIL_DEPTHS = [0.01, 0.04, 0.09, 0.16, 0.26, 0.40, 0.59, 0.83, 1.14, 1.56] + + +def get_scalar_value(var): + """Extract a scalar value from a potentially multi-dimensional variable.""" + data = var[:] + if hasattr(data, 'mask'): + data = np.ma.filled(data, np.nan) + return float(np.squeeze(data)) + + +def get_1d_array(var): + """Extract a 1D array from a variable.""" + data = var[:] + if hasattr(data, 'mask'): + data = np.ma.filled(data, np.nan) + return np.squeeze(data) + + +def create_scalar_card(ax, name, value, units, long_name): + """Create a card-style visualization for a scalar variable.""" + ax.set_xlim(0, 10) + ax.set_ylim(0, 10) + ax.axis('off') + + # Background + rect = mpatches.FancyBboxPatch((0.5, 0.5), 9, 9, boxstyle="round,pad=0.1", + facecolor='#f0f4f8', edgecolor='#2c5282', linewidth=2) + ax.add_patch(rect) + + # Variable name + ax.text(5, 7.5, name, ha='center', va='center', fontsize=14, fontweight='bold', color='#2c5282') + + # Value + if isinstance(value, float): + if abs(value) < 0.01 or abs(value) > 10000: + value_str = f'{value:.4e}' + else: + value_str = f'{value:.4f}' + else: + value_str = str(value) + + ax.text(5, 5, value_str, ha='center', va='center', fontsize=20, fontweight='bold', color='#1a365d') + + # Units + ax.text(5, 3, f'[{units}]', ha='center', va='center', fontsize=10, color='#4a5568') + + # Long name (wrapped) + wrapped_name = '\n'.join([long_name[i:i+30] for i in range(0, len(long_name), 30)]) + ax.text(5, 1.5, wrapped_name, ha='center', va='center', fontsize=8, color='#718096', style='italic') + + +def plot_soil_profile(ax, depths, values, title, units, color='#3182ce'): + """Create a soil profile plot (vertical bar chart).""" + ax.barh(range(len(depths)), values, color=color, edgecolor='#2c5282', alpha=0.8) + ax.set_yticks(range(len(depths))) + ax.set_yticklabels([f'{d:.2f}m' for d in depths]) + ax.invert_yaxis() + ax.set_xlabel(f'{title} [{units}]') + ax.set_ylabel('Depth') + ax.set_title(title, fontsize=12, fontweight='bold') + ax.grid(axis='x', alpha=0.3) + + # Add value labels + for i, v in enumerate(values): + ax.text(v + max(values)*0.02, i, f'{v:.1f}', va='center', fontsize=8) + + +def plot_pie_chart(ax, labels, values, title, min_pct=0.5): + """Create a pie chart with smart label handling.""" + # Filter out zero or very small values + mask = values > min_pct + filtered_labels = [l for l, m in zip(labels, mask) if m] + filtered_values = values[mask] + + if len(filtered_values) == 0: + ax.text(0.5, 0.5, 'No significant values', ha='center', va='center', transform=ax.transAxes) + ax.set_title(title, fontsize=12, fontweight='bold') + return + + colors = plt.cm.tab20(np.linspace(0, 1, len(filtered_values))) + + wedges, texts, autotexts = ax.pie(filtered_values, labels=None, autopct='%1.1f%%', + colors=colors, pctdistance=0.75) + + # Add legend with labels + ax.legend(wedges, [f'{l}: {v:.1f}%' for l, v in zip(filtered_labels, filtered_values)], + loc='center left', bbox_to_anchor=(1, 0.5), fontsize=8) + + ax.set_title(title, fontsize=12, fontweight='bold') + + +def plot_bar_chart(ax, labels, values, title, units, color='#3182ce', horizontal=False): + """Create a bar chart.""" + x = range(len(labels)) + + if horizontal: + ax.barh(x, values, color=color, edgecolor='#2c5282', alpha=0.8) + ax.set_yticks(x) + ax.set_yticklabels(labels, fontsize=8) + ax.set_xlabel(f'{title} [{units}]') + else: + ax.bar(x, values, color=color, edgecolor='#2c5282', alpha=0.8) + ax.set_xticks(x) + ax.set_xticklabels(labels, rotation=45, ha='right', fontsize=8) + ax.set_ylabel(f'[{units}]') + + ax.set_title(title, fontsize=12, fontweight='bold') + ax.grid(axis='y' if not horizontal else 'x', alpha=0.3) + + +def plot_monthly_timeseries(ax, data, pft_indices, title, units, pft_names): + """Plot monthly time series for selected PFTs.""" + colors = plt.cm.tab10(np.linspace(0, 1, len(pft_indices))) + + for i, (pft_idx, color) in enumerate(zip(pft_indices, colors)): + values = data[:, pft_idx, 0, 0] # time, pft, lat, lon + if np.any(values > 0): + label = pft_names.get(pft_idx, f'PFT {pft_idx}') + # Truncate label if too long + if len(label) > 25: + label = label[:22] + '...' + ax.plot(range(12), values, marker='o', label=label, color=color, linewidth=2) + + ax.set_xticks(range(12)) + ax.set_xticklabels(MONTH_NAMES) + ax.set_xlabel('Month') + ax.set_ylabel(f'[{units}]') + ax.set_title(title, fontsize=12, fontweight='bold') + ax.legend(loc='upper left', bbox_to_anchor=(1, 1), fontsize=7) + ax.grid(alpha=0.3) + + +def plot_grouped_bar(ax, data, group_labels, bar_labels, title, units): + """Create a grouped bar chart.""" + x = np.arange(len(group_labels)) + width = 0.25 + n_bars = len(bar_labels) + + colors = plt.cm.Set2(np.linspace(0, 1, n_bars)) + + for i, (bar_label, color) in enumerate(zip(bar_labels, colors)): + offset = (i - n_bars/2 + 0.5) * width + ax.bar(x + offset, data[i], width, label=bar_label, color=color, edgecolor='#2c5282') + + ax.set_xticks(x) + ax.set_xticklabels(group_labels, fontsize=9) + ax.set_ylabel(f'[{units}]') + ax.set_title(title, fontsize=12, fontweight='bold') + ax.legend(fontsize=8) + ax.grid(axis='y', alpha=0.3) + + +def plot_heatmap(ax, data, x_labels, y_labels, title, units, cmap='viridis'): + """Create a heatmap.""" + im = ax.imshow(data, cmap=cmap, aspect='auto') + + ax.set_xticks(range(len(x_labels))) + ax.set_xticklabels(x_labels, fontsize=9) + ax.set_yticks(range(len(y_labels))) + ax.set_yticklabels(y_labels, fontsize=9) + + # Add value annotations + for i in range(len(y_labels)): + for j in range(len(x_labels)): + val = data[i, j] + text_color = 'white' if val > np.mean(data) else 'black' + ax.text(j, i, f'{val:.2f}', ha='center', va='center', color=text_color, fontsize=8) + + ax.set_title(f'{title} [{units}]', fontsize=12, fontweight='bold') + plt.colorbar(im, ax=ax, label=units) + + +def fig_to_base64(fig): + """Convert a matplotlib figure to base64 string for HTML embedding.""" + import io + buf = io.BytesIO() + fig.savefig(buf, format='png', dpi=100, bbox_inches='tight') + buf.seek(0) + return base64.b64encode(buf.read()).decode('utf-8') + + +def create_site_location_figure(lons, lats, labels=None, zoom_margin=10, figsize=(10, 7)): + """Create a figure with a political map showing one or more site locations. + + Parameters + ---------- + lons, lats : lists of float + Longitude and latitude values of the site(s). + labels : list of str, optional + Labels shown next to each site marker. + zoom_margin : float + Degrees of padding around the site(s) on the map. + figsize : tuple + Figure size in inches. + + Returns None if cartopy is not available. + """ + if not HAS_CARTOPY: + return None + + import cartopy.crs as ccrs + import cartopy.feature as cfeature + + fig = plt.figure(figsize=figsize) + ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) + + # Compute map extent + if len(lons) == 1: + lon_c, lat_c = lons[0], lats[0] + extent = [lon_c - zoom_margin, lon_c + zoom_margin, + lat_c - zoom_margin * 0.7, lat_c + zoom_margin * 0.7] + else: + lon_min, lon_max = min(lons), max(lons) + lat_min, lat_max = min(lats), max(lats) + span = max(lon_max - lon_min, lat_max - lat_min) + margin = max(zoom_margin, span + 5) + lon_c = (lon_min + lon_max) / 2 + lat_c = (lat_min + lat_max) / 2 + extent = [lon_c - margin, lon_c + margin, + lat_c - margin * 0.7, lat_c + margin * 0.7] + + ax.set_extent(extent, crs=ccrs.PlateCarree()) + + # Add map features + ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='#c8e6f5') + ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='#f2efe8') + ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8, edgecolor='#555555') + ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.7, edgecolor='#888888', linestyle='-') + ax.add_feature(cfeature.LAKES.with_scale('50m'), facecolor='#c8e6f5', alpha=0.8) + ax.add_feature(cfeature.RIVERS.with_scale('50m'), linewidth=0.4, edgecolor='#6baed6', alpha=0.6) + + # Plot site markers + site_colors = ['#e53e3e', '#3182ce', '#38a169', '#d69e2e', '#805ad5'] + for i, (lon, lat) in enumerate(zip(lons, lats)): + color = site_colors[i % len(site_colors)] + ax.plot(lon, lat, marker='*', color=color, markersize=20, + transform=ccrs.PlateCarree(), zorder=5, + markeredgecolor='white', markeredgewidth=0.5) + label_text = labels[i] if labels else f'({lon:.3f}°E, {lat:.3f}°N)' + ax.text(lon + 0.3, lat + 0.3, label_text, + transform=ccrs.PlateCarree(), fontsize=9, fontweight='bold', + color=color, + bbox=dict(boxstyle='round,pad=0.3', facecolor='white', + edgecolor=color, linewidth=1.2, alpha=0.9), + zorder=6) + + # Gridlines + gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, + linewidth=0.5, color='gray', alpha=0.5, linestyle='--') + gl.top_labels = False + gl.right_labels = False + + if len(lons) == 1: + ax.set_title(f'Site Location ({lons[0]:.3f}°E, {lats[0]:.3f}°N)', + fontsize=13, fontweight='bold') + else: + ax.set_title('Site Locations', fontsize=13, fontweight='bold') + + return fig + + +def create_html_report(figures_data, nc_file, output_dir): + """Generate an HTML report with all figures embedded.""" + + html_template = """ + + + + + Surface Data Visualization Report + + + +
    +
    +

    CLM5 Surface Data Visualization Report

    +

    Comprehensive analysis of land surface parameters for model input review

    +
    + Source File: {nc_file}
    + Generated: {timestamp}
    + Site: BE-Lon (Belgium - Lonzee)
    + Coordinates: {lon:.3f}E, {lat:.3f}N +
    +
    + + + + {sections} + +
    +

    Generated by CLM5 Surface Data Visualization Tool

    +

    For discussion of model input parameters

    +
    +
    + +""" + + sections_html = "" + nav_links_html = "" + + for section in figures_data: + section_id = section['id'] + section_title = section['title'] + section_desc = section.get('description', '') + + nav_links_html += f'
  • {section_title}
    • \n' + + figures_html = "" + for fig_data in section['figures']: + pdf_filename = fig_data['pdf_name'] + caption = fig_data['caption'] + img_base64 = fig_data['base64'] + + figures_html += f""" +
    + {caption} +

    {caption}

    + Download PDF +
    + """ + + sections_html += f""" +
    +

    {section_title}

    +

    {section_desc}

    + {figures_html} +
    + """ + + # Get coordinates from the first figure data + lon = figures_data[0].get('lon', 0) + lat = figures_data[0].get('lat', 0) + + html_content = html_template.format( + nc_file=os.path.basename(nc_file), + timestamp=datetime.now().strftime('%Y-%m-%d %H:%M:%S'), + lon=lon, + lat=lat, + nav_links=nav_links_html, + sections=sections_html + ) + + # Create HTML filename based on NC filename (replace .nc with .html) + nc_basename = os.path.basename(nc_file) + if nc_basename.endswith('.nc'): + html_basename = nc_basename[:-3] + '.html' + else: + html_basename = nc_basename + '.html' + + output_file = os.path.join(output_dir, html_basename) + with open(output_file, 'w') as f: + f.write(html_content) + + return output_file + + +def main(nc_file): + """Main function to create all visualizations.""" + + if not os.path.exists(nc_file): + print(f"Error: File not found: {nc_file}") + sys.exit(1) + + # Create output directory based on NC filename (replace .nc with _figures) + output_dir = os.path.dirname(nc_file) + if not output_dir: + output_dir = '.' + + nc_basename = os.path.basename(nc_file) + if nc_basename.endswith('.nc'): + figures_dirname = nc_basename[:-3] + '_figures' + else: + figures_dirname = nc_basename + '_figures' + + pdf_dir = os.path.join(output_dir, figures_dirname) + os.makedirs(pdf_dir, exist_ok=True) + + print(f"Reading NetCDF file: {nc_file}") + nc = Dataset(nc_file, 'r') + + # Get coordinates + lon = get_scalar_value(nc.variables['LONGXY']) + lat = get_scalar_value(nc.variables['LATIXY']) + + figures_data = [] + + # ========================================================================= + # SECTION 1: Location and Basic Parameters + # ========================================================================= + print("Creating Section 1: Location and Basic Parameters...") + section1_figures = [] + + # Figure 1.0: Site location map + fig_map = create_site_location_figure([lon], [lat]) + if fig_map is not None: + pdf_path = os.path.join(pdf_dir, '00_site_location.pdf') + fig_map.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': f'Political map showing the site location at {lon:.3f}°E, {lat:.3f}°N.', + 'base64': fig_to_base64(fig_map) + }) + plt.close(fig_map) + + # Figure 1.1: Location and basic info + fig, axes = plt.subplots(3, 4, figsize=(16, 12)) + fig.suptitle('Basic Site Parameters', fontsize=16, fontweight='bold') + + scalar_vars = [ + ('LONGXY', 'degrees east'), + ('LATIXY', 'degrees north'), + ('AREA', 'km^2'), + ('FMAX', 'unitless'), + ('SOIL_COLOR', 'index'), + ('mxsoil_color', 'unitless'), + ('zbedrock', 'm'), + ('SLOPE', 'degrees'), + ('STD_ELEV', 'm'), + ('LAKEDEPTH', 'm'), + ('gdp', 'unitless'), + ('abm', 'month') + ] + + for ax, (var_name, units) in zip(axes.flatten(), scalar_vars): + var = nc.variables[var_name] + value = get_scalar_value(var) + long_name = var.long_name if hasattr(var, 'long_name') else var_name + create_scalar_card(ax, var_name, value, units, long_name) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '01_basic_parameters.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Basic site parameters including location, area, soil color, bedrock depth, slope, and economic indicators.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'basic-params', + 'title': 'Basic Site Parameters', + 'description': 'Fundamental parameters describing the site location, topography, and basic soil/lake properties.', + 'figures': section1_figures, + 'lon': lon, + 'lat': lat + }) + + # ========================================================================= + # SECTION 2: Land Cover Fractions + # ========================================================================= + print("Creating Section 2: Land Cover Fractions...") + section2_figures = [] + + # Figure 2.1: Major land cover pie chart + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + + # Major land units + pct_natveg = get_scalar_value(nc.variables['PCT_NATVEG']) + pct_crop = get_scalar_value(nc.variables['PCT_CROP']) + pct_urban = np.sum(get_1d_array(nc.variables['PCT_URBAN'])) + pct_lake = get_scalar_value(nc.variables['PCT_LAKE']) + pct_wetland = get_scalar_value(nc.variables['PCT_WETLAND']) + pct_glacier = get_scalar_value(nc.variables['PCT_GLACIER']) + + land_labels = ['Natural Vegetation', 'Cropland', 'Urban', 'Lake', 'Wetland', 'Glacier'] + land_values = np.array([pct_natveg, pct_crop, pct_urban, pct_lake, pct_wetland, pct_glacier]) + + plot_pie_chart(axes[0], land_labels, land_values, 'Major Land Cover Types', min_pct=0.1) + + # Land fractions bar chart + axes[1].bar(land_labels, land_values, color=['#48bb78', '#ecc94b', '#a0aec0', '#4299e1', '#9f7aea', '#63b3ed'], + edgecolor='#2d3748') + axes[1].set_ylabel('Percent (%)') + axes[1].set_title('Land Cover Distribution', fontsize=12, fontweight='bold') + axes[1].tick_params(axis='x', rotation=45) + for i, v in enumerate(land_values): + if v > 0: + axes[1].text(i, v + 1, f'{v:.1f}%', ha='center', fontsize=9) + axes[1].grid(axis='y', alpha=0.3) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '02_land_cover_major.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section2_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Distribution of major land cover types: natural vegetation, cropland, urban, lake, wetland, and glacier.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'land-cover', + 'title': 'Land Cover Fractions', + 'description': 'Overview of land surface cover type distribution. Key variables for understanding the dominant land use at the site.', + 'figures': section2_figures + }) + + # ========================================================================= + # SECTION 3: Natural PFT Distribution + # ========================================================================= + print("Creating Section 3: Natural PFT Distribution...") + section3_figures = [] + + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + + # Get PCT_NAT_PFT data + pct_nat_pft = get_1d_array(nc.variables['PCT_NAT_PFT']) + natpft_labels = [NATPFT_NAMES.get(i, f'PFT {i}') for i in range(15)] + + # Pie chart for non-zero PFTs + plot_pie_chart(axes[0], natpft_labels, pct_nat_pft, + f'Natural PFT Distribution\n(within {pct_natveg:.1f}% natural veg)', min_pct=0.5) + + # Horizontal bar chart + colors = plt.cm.Greens(np.linspace(0.3, 0.9, 15)) + y_pos = range(15) + axes[1].barh(y_pos, pct_nat_pft, color=colors, edgecolor='#2d3748') + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(natpft_labels, fontsize=9) + axes[1].set_xlabel('Percent of Natural Vegetation Landunit (%)') + axes[1].set_title('Natural Plant Functional Types', fontsize=12, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + + # Add value labels for non-zero + for i, v in enumerate(pct_nat_pft): + if v > 0: + axes[1].text(v + 1, i, f'{v:.1f}%', va='center', fontsize=8) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '03_natural_pft.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section3_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Distribution of natural Plant Functional Types (PFTs). Values represent percentages within the natural vegetation landunit.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'natural-pft', + 'title': 'Natural PFT Distribution', + 'description': 'Natural Plant Functional Type distribution within the natural vegetation landunit. Important for understanding vegetation dynamics.', + 'figures': section3_figures + }) + + # ========================================================================= + # SECTION 4: Crop Functional Types + # ========================================================================= + print("Creating Section 4: Crop Functional Types...") + section4_figures = [] + + # Get PCT_CFT data + pct_cft = get_1d_array(nc.variables['PCT_CFT']) + cft_indices = nc.variables['cft'][:] + + # Only show non-zero CFTs + nonzero_mask = pct_cft > 0.1 + nonzero_cft_values = pct_cft[nonzero_mask] + nonzero_cft_indices = cft_indices[nonzero_mask] + nonzero_cft_labels = [CFT_NAMES.get(int(i), f'CFT {i}') for i in nonzero_cft_indices] + + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + + if len(nonzero_cft_values) > 0: + # Pie chart + plot_pie_chart(axes[0], nonzero_cft_labels, nonzero_cft_values, + f'Crop Type Distribution\n(within {pct_crop:.1f}% cropland)', min_pct=0.1) + + # Bar chart + colors = plt.cm.YlOrBr(np.linspace(0.3, 0.9, len(nonzero_cft_values))) + y_pos = range(len(nonzero_cft_values)) + axes[1].barh(y_pos, nonzero_cft_values, color=colors, edgecolor='#2d3748') + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(nonzero_cft_labels, fontsize=9) + axes[1].set_xlabel('Percent of Crop Landunit (%)') + axes[1].set_title('Crop Functional Types (non-zero only)', fontsize=12, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + + for i, v in enumerate(nonzero_cft_values): + axes[1].text(v + 0.5, i, f'{v:.1f}%', va='center', fontsize=9) + else: + axes[0].text(0.5, 0.5, 'No crops present', ha='center', va='center', transform=axes[0].transAxes) + axes[1].text(0.5, 0.5, 'No crops present', ha='center', va='center', transform=axes[1].transAxes) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04_crop_functional_types.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Distribution of Crop Functional Types (CFTs). Values represent percentages within the crop landunit.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Nitrogen fertilizer for crops + fig, ax = plt.subplots(figsize=(14, 6)) + fertnitro = get_1d_array(nc.variables['CONST_FERTNITRO_CFT']) + nonzero_fert_mask = fertnitro > 0 + + if np.any(nonzero_fert_mask): + nonzero_fert = fertnitro[nonzero_fert_mask] + nonzero_fert_indices = cft_indices[nonzero_fert_mask] + nonzero_fert_labels = [CFT_NAMES.get(int(i), f'CFT {i}') for i in nonzero_fert_indices] + + colors = plt.cm.Blues(np.linspace(0.4, 0.9, len(nonzero_fert))) + x_pos = range(len(nonzero_fert)) + ax.bar(x_pos, nonzero_fert, color=colors, edgecolor='#2d3748') + ax.set_xticks(x_pos) + ax.set_xticklabels(nonzero_fert_labels, rotation=45, ha='right', fontsize=9) + ax.set_ylabel('Nitrogen Fertilizer (gN/m2/yr)') + ax.set_title('Nitrogen Fertilizer Application by Crop Type', fontsize=12, fontweight='bold') + ax.grid(axis='y', alpha=0.3) + + for i, v in enumerate(nonzero_fert): + ax.text(i, v + 0.3, f'{v:.1f}', ha='center', fontsize=9) + else: + ax.text(0.5, 0.5, 'No fertilizer data', ha='center', va='center', transform=ax.transAxes) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04b_nitrogen_fertilizer.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Nitrogen fertilizer application rates for different crop types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'crop-types', + 'title': 'Crop Functional Types', + 'description': 'Crop Functional Type (CFT) distribution and associated nitrogen fertilizer application rates.', + 'figures': section4_figures + }) + + # ========================================================================= + # SECTION 5: Soil Properties + # ========================================================================= + print("Creating Section 5: Soil Properties...") + section5_figures = [] + + fig, axes = plt.subplots(1, 3, figsize=(15, 6)) + + # PCT_SAND + pct_sand = get_1d_array(nc.variables['PCT_SAND']) + plot_soil_profile(axes[0], SOIL_DEPTHS, pct_sand, 'Sand Content', '%', '#f6ad55') + + # PCT_CLAY + pct_clay = get_1d_array(nc.variables['PCT_CLAY']) + plot_soil_profile(axes[1], SOIL_DEPTHS, pct_clay, 'Clay Content', '%', '#fc8181') + + # ORGANIC + organic = get_1d_array(nc.variables['ORGANIC']) + plot_soil_profile(axes[2], SOIL_DEPTHS, organic, 'Organic Matter', 'kg/m3', '#68d391') + + fig.suptitle('Soil Properties by Depth', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05_soil_properties.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Vertical profiles of soil texture (sand and clay content) and organic matter density across 10 soil layers.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Combined soil texture visualization + fig, ax = plt.subplots(figsize=(10, 8)) + pct_silt = 100 - pct_sand - pct_clay # Calculate silt + + # Stacked bar chart for soil texture + x = range(len(SOIL_DEPTHS)) + width = 0.6 + + ax.barh(x, pct_sand, width, label='Sand', color='#f6ad55', edgecolor='#2d3748') + ax.barh(x, pct_silt, width, left=pct_sand, label='Silt', color='#a0aec0', edgecolor='#2d3748') + ax.barh(x, pct_clay, width, left=pct_sand+pct_silt, label='Clay', color='#fc8181', edgecolor='#2d3748') + + ax.set_yticks(x) + ax.set_yticklabels([f'{d:.2f}m' for d in SOIL_DEPTHS]) + ax.invert_yaxis() + ax.set_xlabel('Percent (%)') + ax.set_ylabel('Depth') + ax.set_title('Soil Texture Composition by Depth', fontsize=12, fontweight='bold') + ax.legend(loc='upper right') + ax.set_xlim(0, 100) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05b_soil_texture.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Stacked bar chart showing soil texture composition (sand, silt, clay) at each soil layer.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'soil', + 'title': 'Soil Properties', + 'description': 'Soil texture and organic matter content profiles. Critical inputs for soil hydrology and biogeochemistry.', + 'figures': section5_figures + }) + + # ========================================================================= + # SECTION 6: Monthly LAI and Vegetation Parameters + # ========================================================================= + print("Creating Section 6: Monthly Vegetation Parameters...") + section6_figures = [] + + # Find PFTs with non-zero LAI + lai_data = nc.variables['MONTHLY_LAI'][:] + sai_data = nc.variables['MONTHLY_SAI'][:] + + # Identify active PFTs (those with non-zero LAI at any month) + active_pfts = [] + for pft in range(79): + if np.any(lai_data[:, pft, 0, 0] > 0): + active_pfts.append(pft) + + # Combine PFT names + all_pft_names = {**NATPFT_NAMES, **CFT_NAMES} + + # Plot LAI + fig, ax = plt.subplots(figsize=(14, 7)) + plot_monthly_timeseries(ax, lai_data, active_pfts[:10], + 'Monthly Leaf Area Index (LAI)', 'm2/m2', all_pft_names) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06_monthly_lai.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Monthly Leaf Area Index (LAI) time series for active Plant Functional Types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Plot SAI + fig, ax = plt.subplots(figsize=(14, 7)) + plot_monthly_timeseries(ax, sai_data, active_pfts[:10], + 'Monthly Stem Area Index (SAI)', 'm2/m2', all_pft_names) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06b_monthly_sai.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Monthly Stem Area Index (SAI) time series for active Plant Functional Types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Plot canopy heights + height_top = nc.variables['MONTHLY_HEIGHT_TOP'][:] + height_bot = nc.variables['MONTHLY_HEIGHT_BOT'][:] + + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + plot_monthly_timeseries(axes[0], height_top, active_pfts[:10], + 'Monthly Canopy Top Height', 'm', all_pft_names) + plot_monthly_timeseries(axes[1], height_bot, active_pfts[:10], + 'Monthly Canopy Bottom Height', 'm', all_pft_names) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06c_canopy_heights.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Monthly canopy top and bottom heights for active Plant Functional Types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'vegetation', + 'title': 'Monthly Vegetation Parameters', + 'description': 'Seasonal cycle of vegetation properties including LAI, SAI, and canopy heights.', + 'figures': section6_figures + }) + + # ========================================================================= + # SECTION 7: Urban Parameters + # ========================================================================= + print("Creating Section 7: Urban Parameters...") + section7_figures = [] + + # Urban percent by type + pct_urban = get_1d_array(nc.variables['PCT_URBAN']) + + fig, axes = plt.subplots(2, 3, figsize=(15, 10)) + + # Urban fraction + axes[0, 0].bar(URBAN_TYPES, pct_urban, color=['#e53e3e', '#dd6b20', '#d69e2e'], edgecolor='#2d3748') + axes[0, 0].set_ylabel('Percent (%)') + axes[0, 0].set_title('Urban Fraction by Density Type', fontsize=11, fontweight='bold') + axes[0, 0].tick_params(axis='x', rotation=15) + for i, v in enumerate(pct_urban): + axes[0, 0].text(i, v + 0.1, f'{v:.2f}%', ha='center', fontsize=9) + + # Canyon height-to-width ratio + canyon_hwr = get_1d_array(nc.variables['CANYON_HWR']) + axes[0, 1].bar(URBAN_TYPES, canyon_hwr, color='#805ad5', edgecolor='#2d3748') + axes[0, 1].set_ylabel('Height/Width Ratio') + axes[0, 1].set_title('Canyon Height-to-Width Ratio', fontsize=11, fontweight='bold') + axes[0, 1].tick_params(axis='x', rotation=15) + + # Roof height + ht_roof = get_1d_array(nc.variables['HT_ROOF']) + axes[0, 2].bar(URBAN_TYPES, ht_roof, color='#3182ce', edgecolor='#2d3748') + axes[0, 2].set_ylabel('Height (m)') + axes[0, 2].set_title('Roof Height', fontsize=11, fontweight='bold') + axes[0, 2].tick_params(axis='x', rotation=15) + + # Building temperature minimum + t_building = get_1d_array(nc.variables['T_BUILDING_MIN']) + axes[1, 0].bar(URBAN_TYPES, t_building - 273.15, color='#e53e3e', edgecolor='#2d3748') # Convert to Celsius + axes[1, 0].set_ylabel('Temperature (C)') + axes[1, 0].set_title('Min. Building Interior Temp.', fontsize=11, fontweight='bold') + axes[1, 0].tick_params(axis='x', rotation=15) + + # Roof and pervious road fractions + wtlunit_roof = get_1d_array(nc.variables['WTLUNIT_ROOF']) + wtroad_perv = get_1d_array(nc.variables['WTROAD_PERV']) + + x = np.arange(len(URBAN_TYPES)) + width = 0.35 + axes[1, 1].bar(x - width/2, wtlunit_roof, width, label='Roof Fraction', color='#805ad5', edgecolor='#2d3748') + axes[1, 1].bar(x + width/2, wtroad_perv, width, label='Pervious Road Fraction', color='#38a169', edgecolor='#2d3748') + axes[1, 1].set_xticks(x) + axes[1, 1].set_xticklabels(URBAN_TYPES, rotation=15) + axes[1, 1].set_ylabel('Fraction') + axes[1, 1].set_title('Urban Surface Fractions', fontsize=11, fontweight='bold') + axes[1, 1].legend(fontsize=8) + + # Wall and roof thickness + thick_roof = get_1d_array(nc.variables['THICK_ROOF']) + thick_wall = get_1d_array(nc.variables['THICK_WALL']) + + axes[1, 2].bar(x - width/2, thick_roof, width, label='Roof Thickness', color='#4299e1', edgecolor='#2d3748') + axes[1, 2].bar(x + width/2, thick_wall, width, label='Wall Thickness', color='#ed8936', edgecolor='#2d3748') + axes[1, 2].set_xticks(x) + axes[1, 2].set_xticklabels(URBAN_TYPES, rotation=15) + axes[1, 2].set_ylabel('Thickness (m)') + axes[1, 2].set_title('Building Element Thickness', fontsize=11, fontweight='bold') + axes[1, 2].legend(fontsize=8) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '07_urban_basic.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section7_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Basic urban parameters for three density classes: Tall Building District, High Density, and Medium Density.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Urban radiative properties + fig, axes = plt.subplots(2, 2, figsize=(14, 10)) + + # Emissivities + em_vars = ['EM_IMPROAD', 'EM_PERROAD', 'EM_ROOF', 'EM_WALL'] + em_data = [get_1d_array(nc.variables[v]) for v in em_vars] + em_labels = ['Impervious Road', 'Pervious Road', 'Roof', 'Wall'] + + x = np.arange(len(URBAN_TYPES)) + width = 0.2 + colors = plt.cm.Set2(np.linspace(0, 1, 4)) + + for i, (data, label, color) in enumerate(zip(em_data, em_labels, colors)): + axes[0, 0].bar(x + (i - 1.5) * width, data, width, label=label, color=color, edgecolor='#2d3748') + + axes[0, 0].set_xticks(x) + axes[0, 0].set_xticklabels(URBAN_TYPES, rotation=15) + axes[0, 0].set_ylabel('Emissivity') + axes[0, 0].set_title('Surface Emissivities', fontsize=11, fontweight='bold') + axes[0, 0].legend(fontsize=8) + axes[0, 0].set_ylim(0, 1) + + # Direct albedos (VIS band) + alb_vars_dir = ['ALB_IMPROAD_DIR', 'ALB_PERROAD_DIR', 'ALB_ROOF_DIR', 'ALB_WALL_DIR'] + alb_data_dir = [] + for v in alb_vars_dir: + data = nc.variables[v][:] + alb_data_dir.append(data[0, :, 0, 0]) # VIS band + + for i, (data, label, color) in enumerate(zip(alb_data_dir, em_labels, colors)): + axes[0, 1].bar(x + (i - 1.5) * width, data, width, label=label, color=color, edgecolor='#2d3748') + + axes[0, 1].set_xticks(x) + axes[0, 1].set_xticklabels(URBAN_TYPES, rotation=15) + axes[0, 1].set_ylabel('Albedo (VIS)') + axes[0, 1].set_title('Direct Albedo (Visible Band)', fontsize=11, fontweight='bold') + axes[0, 1].legend(fontsize=8) + + # Thermal conductivity (first level) + tk_vars = ['TK_ROOF', 'TK_WALL', 'TK_IMPROAD'] + tk_labels = ['Roof', 'Wall', 'Impervious Road'] + tk_data = [] + for v in tk_vars: + data = nc.variables[v][:] + tk_data.append(data[0, :, 0, 0]) # First level + + for i, (data, label) in enumerate(zip(tk_data, tk_labels)): + axes[1, 0].bar(x + (i - 1) * 0.25, data, 0.25, label=label, + color=plt.cm.coolwarm(i/2), edgecolor='#2d3748') + + axes[1, 0].set_xticks(x) + axes[1, 0].set_xticklabels(URBAN_TYPES, rotation=15) + axes[1, 0].set_ylabel('Thermal Conductivity (W/m*K)') + axes[1, 0].set_title('Thermal Conductivity (Layer 1)', fontsize=11, fontweight='bold') + axes[1, 0].legend(fontsize=8) + + # Heat capacity (first level) + cv_vars = ['CV_ROOF', 'CV_WALL', 'CV_IMPROAD'] + cv_labels = ['Roof', 'Wall', 'Impervious Road'] + cv_data = [] + for v in cv_vars: + data = nc.variables[v][:] + cv_data.append(data[0, :, 0, 0] / 1e6) # Convert to MJ for readability + + for i, (data, label) in enumerate(zip(cv_data, cv_labels)): + axes[1, 1].bar(x + (i - 1) * 0.25, data, 0.25, label=label, + color=plt.cm.coolwarm(i/2), edgecolor='#2d3748') + + axes[1, 1].set_xticks(x) + axes[1, 1].set_xticklabels(URBAN_TYPES, rotation=15) + axes[1, 1].set_ylabel('Heat Capacity (MJ/m3*K)') + axes[1, 1].set_title('Volumetric Heat Capacity (Layer 1)', fontsize=11, fontweight='bold') + axes[1, 1].legend(fontsize=8) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '07b_urban_thermal.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section7_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Urban radiative and thermal properties: emissivities, albedos, thermal conductivity, and heat capacity.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'urban', + 'title': 'Urban Parameters', + 'description': 'Urban canyon and building parameters for three density classes. Important for urban energy balance simulations.', + 'figures': section7_figures + }) + + # ========================================================================= + # SECTION 8: Emission Factors and Other Parameters + # ========================================================================= + print("Creating Section 8: Emission Factors and Other Parameters...") + section8_figures = [] + + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + + # Emission factors (isoprene) + ef_vars = ['EF1_BTR', 'EF1_FET', 'EF1_FDT', 'EF1_SHR', 'EF1_GRS', 'EF1_CRP'] + ef_labels = ['Broadleaf Trees', 'Fineleaf Evergreen Trees', 'Fineleaf Deciduous Trees', + 'Shrubs', 'Grasses', 'Crops'] + ef_values = [get_scalar_value(nc.variables[v]) for v in ef_vars] + + colors = plt.cm.Greens(np.linspace(0.3, 0.9, len(ef_vars))) + axes[0].bar(ef_labels, ef_values, color=colors, edgecolor='#2d3748') + axes[0].set_ylabel('Emission Factor') + axes[0].set_title('Isoprene Emission Factors by Vegetation Type', fontsize=11, fontweight='bold') + axes[0].tick_params(axis='x', rotation=45) + axes[0].set_yscale('log') + for i, v in enumerate(ef_values): + axes[0].text(i, v * 1.1, f'{v:.0f}', ha='center', fontsize=9) + + # Glacier elevation classes + glc_mec = get_1d_array(nc.variables['GLC_MEC']) + pct_glc_mec = get_1d_array(nc.variables['PCT_GLC_MEC']) + topo_glc_mec = get_1d_array(nc.variables['TOPO_GLC_MEC']) + + ax2 = axes[1].twinx() + x = range(len(pct_glc_mec)) + axes[1].bar(x, pct_glc_mec, color='#63b3ed', alpha=0.7, label='Glacier %', edgecolor='#2d3748') + ax2.plot(x, topo_glc_mec, 'ro-', label='Mean Elevation', linewidth=2, markersize=8) + + axes[1].set_xticks(x) + axes[1].set_xticklabels([f'{glc_mec[i]:.0f}-{glc_mec[i+1]:.0f}m' for i in range(len(glc_mec)-1)], rotation=45) + axes[1].set_xlabel('Elevation Class') + axes[1].set_ylabel('Glacier Percent (%)', color='#3182ce') + ax2.set_ylabel('Mean Elevation (m)', color='red') + axes[1].set_title('Glacier Elevation Classes', fontsize=11, fontweight='bold') + axes[1].legend(loc='upper left') + ax2.legend(loc='upper right') + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '08_emission_glacier.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section8_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Isoprene emission factors by vegetation type and glacier distribution across elevation classes.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Harvest parameters + fig, ax = plt.subplots(figsize=(10, 6)) + + harvest_vars = ['CONST_HARVEST_VH1', 'CONST_HARVEST_VH2', 'CONST_HARVEST_SH1', + 'CONST_HARVEST_SH2', 'CONST_HARVEST_SH3'] + harvest_labels = ['Primary Forest', 'Primary Non-Forest', 'Secondary Mature Forest', + 'Secondary Young Forest', 'Secondary Non-Forest'] + harvest_values = [get_scalar_value(nc.variables[v]) for v in harvest_vars] + + colors = plt.cm.YlOrBr(np.linspace(0.3, 0.9, len(harvest_vars))) + ax.bar(harvest_labels, harvest_values, color=colors, edgecolor='#2d3748') + ax.set_ylabel('Harvest Rate (gC/m2/yr)') + ax.set_title('Constant Harvest Rates by Land Type', fontsize=12, fontweight='bold') + ax.tick_params(axis='x', rotation=30) + ax.grid(axis='y', alpha=0.3) + + for i, v in enumerate(harvest_values): + if v > 0: + ax.text(i, v + 10, f'{v:.1f}', ha='center', fontsize=9) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '08b_harvest.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section8_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Constant harvest rates for different land cover types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'other', + 'title': 'Emission Factors & Other Parameters', + 'description': 'Biogenic emission factors, glacier elevation classes, and harvest parameters.', + 'figures': section8_figures + }) + + # ========================================================================= + # SECTION 9: Summary Overview + # ========================================================================= + print("Creating Section 9: Summary Overview...") + section9_figures = [] + + # Create a comprehensive summary figure + fig = plt.figure(figsize=(16, 12)) + + # Land cover pie (small) + ax1 = fig.add_subplot(2, 3, 1) + nonzero_land = land_values > 0.1 + ax1.pie(land_values[nonzero_land], labels=[l for l, m in zip(land_labels, nonzero_land) if m], + autopct='%1.1f%%', colors=plt.cm.Set3(np.linspace(0, 1, sum(nonzero_land)))) + ax1.set_title('Land Cover Types', fontsize=11, fontweight='bold') + + # Soil texture summary + ax2 = fig.add_subplot(2, 3, 2) + mean_sand = np.mean(pct_sand) + mean_clay = np.mean(pct_clay) + mean_silt = 100 - mean_sand - mean_clay + ax2.pie([mean_sand, mean_silt, mean_clay], labels=['Sand', 'Silt', 'Clay'], + autopct='%1.1f%%', colors=['#f6ad55', '#a0aec0', '#fc8181']) + ax2.set_title('Mean Soil Texture', fontsize=11, fontweight='bold') + + # Active PFTs summary + ax3 = fig.add_subplot(2, 3, 3) + active_labels = [all_pft_names.get(p, f'PFT {p}')[:20] for p in active_pfts[:8]] + active_lai_mean = [np.mean(lai_data[:, p, 0, 0]) for p in active_pfts[:8]] + ax3.barh(range(len(active_labels)), active_lai_mean, color=plt.cm.Greens(np.linspace(0.3, 0.9, len(active_labels)))) + ax3.set_yticks(range(len(active_labels))) + ax3.set_yticklabels(active_labels, fontsize=8) + ax3.set_xlabel('Mean LAI') + ax3.set_title('Top Active PFTs by LAI', fontsize=11, fontweight='bold') + + # Key scalars table + ax4 = fig.add_subplot(2, 3, (4, 6)) + ax4.axis('off') + + key_params = [ + ('Longitude', f'{lon:.3f} E'), + ('Latitude', f'{lat:.3f} N'), + ('Natural Vegetation', f'{pct_natveg:.1f}%'), + ('Cropland', f'{pct_crop:.1f}%'), + ('Urban', f'{np.sum(pct_urban):.2f}%'), + ('Bedrock Depth', f'{get_scalar_value(nc.variables["zbedrock"]):.2f} m'), + ('Slope', f'{get_scalar_value(nc.variables["SLOPE"]):.2f} deg'), + ('Soil Color Index', f'{int(get_scalar_value(nc.variables["SOIL_COLOR"]))}'), + ('Max Saturated Fraction', f'{get_scalar_value(nc.variables["FMAX"]):.3f}'), + ('Peatland Fraction', f'{get_scalar_value(nc.variables["peatf"]):.3f}'), + ] + + table_data = [[p[0], p[1]] for p in key_params] + table = ax4.table(cellText=table_data, colLabels=['Parameter', 'Value'], + loc='center', cellLoc='left', colWidths=[0.4, 0.3]) + table.auto_set_font_size(False) + table.set_fontsize(11) + table.scale(1.2, 1.8) + + # Style the table + for i in range(len(table_data) + 1): + for j in range(2): + cell = table[(i, j)] + if i == 0: + cell.set_facecolor('#667eea') + cell.set_text_props(color='white', fontweight='bold') + elif i % 2 == 0: + cell.set_facecolor('#f7fafc') + else: + cell.set_facecolor('#edf2f7') + + ax4.set_title('Key Site Parameters Summary', fontsize=12, fontweight='bold', pad=20) + + fig.suptitle('Surface Data Overview - BE-Lon Site', fontsize=16, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '09_summary.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section9_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Comprehensive summary of key surface data parameters for the BE-Lon site.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + figures_data.append({ + 'id': 'summary', + 'title': 'Summary Overview', + 'description': 'A comprehensive overview combining key parameters for quick reference during discussions.', + 'figures': section9_figures + }) + + # Close NetCDF file + nc.close() + + # Generate HTML report + print("Generating HTML report...") + html_file = create_html_report(figures_data, nc_file, output_dir) + + print(f"\nVisualization complete!") + print(f"PDF figures saved in: {pdf_dir}") + print(f"HTML report saved as: {html_file}") + + return html_file + + +if __name__ == '__main__': + if len(sys.argv) < 2: + # Default to the NC file in the same directory + script_dir = os.path.dirname(os.path.abspath(__file__)) + nc_file = os.path.join(script_dir, 'surfdata_BE-Lon_hist_78pfts_CMIP6_simyr2005_c260203.nc') + if not os.path.exists(nc_file): + print("Usage: python visualize_surfdata.py ") + print("No default file found.") + sys.exit(1) + else: + nc_file = sys.argv[1] + + main(nc_file) From 736b7fdc6879490281d64dfac62f904d3642595f Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 5 Mar 2026 15:27:09 +0100 Subject: [PATCH 2/3] vizsurfdata: Scripts adapted for distributed surface data --- vizsurfdata/README.md | 68 +- vizsurfdata/compare_surfdata.py | 1262 +++++++++++++++-------- vizsurfdata/requirements.txt | 2 +- vizsurfdata/vizualize_surfdata.py | 1571 ++++++++++++++++++++--------- 4 files changed, 2007 insertions(+), 896 deletions(-) diff --git a/vizsurfdata/README.md b/vizsurfdata/README.md index 4a75861..f6fe3b1 100644 --- a/vizsurfdata/README.md +++ b/vizsurfdata/README.md @@ -1,40 +1,84 @@ # Python scripts for visualizing and comparing eCLM surface data files +Both scripts auto-detect whether the input file(s) are **single-site** +(`lsmlat = lsmlon = 1`) or a **regional grid** and activate the appropriate +visualisation mode automatically. + ## Scripts ### visualize_surfdata.py -Visualizes all variables in a CLM5 surface data NetCDF file, generating PDF figures and an HTML report. +Visualizes all variables in a CLM5 surface data NetCDF file, generating PDF +figures and an HTML report. ```bash python visualize_surfdata.py ``` **Output:** -- `_figures/` - Directory containing PDF figures -- `.html` - Interactive HTML report with embedded figures +- `_figures/` – directory containing one PDF per section +- `.html` – self-contained HTML report with all figures embedded + +**Sections and mode-specific behaviour:** -**Sections:** Site location map, basic parameters, land cover, natural PFTs, crop types, soil profiles, monthly LAI/SAI, urban parameters, emission factors, and summary. +| # | Section | Single-site | Regional grid | +|---|--------------------------|---------------------------------------------------------------|-----------------------------------------------------------------------------------------------------------------------| +| 1 | Domain / site overview | Cartopy point map + scalar cards | Cartopy bounding-box map + spatial maps of AREA, FMAX, SOIL\_COLOR, zbedrock, SLOPE, STD\_ELEV, LAKEDEPTH, peatf, gdp | +| 2 | Land cover fractions | Pie chart + bar chart | Spatial maps of PCT\_NATVEG, PCT\_CROP, PCT\_URBAN (total), PCT\_LAKE, PCT\_WETLAND, PCT\_GLACIER | +| 3 | Natural PFTs | Pie chart + horizontal bar chart | Spatial maps for each active PFT (max > 1 %) + domain-mean bar chart | +| 4 | Crop functional types | Pie chart + bar chart + fertiliser bar | Spatial maps for active CFTs + domain-mean fertiliser bar | +| 5 | Soil properties | Vertical profiles (sand, clay, organic) + stacked texture bar | Depth-mean spatial maps + domain-mean profiles with spatial std | +| 6 | Monthly vegetation | LAI / SAI / height time series per PFT | Domain-mean LAI / SAI time series + annual-mean LAI spatial maps per PFT | +| 7 | Urban parameters | Bar charts per density class | PCT\_URBAN spatial maps + domain-mean bar charts for building parameters | +| 8 | Emission factors & other | Bar charts + glacier elevation chart + harvest bars | Spatial maps of EF and harvest variables + domain-mean glacier chart | +| 9 | Summary | Pie charts + scalar table | Grid of key spatial maps | ### compare_surfdata.py -Compares two surface data files side-by-side, highlighting differences with color-coded visualizations. +Compares two CLM5 surface data files side-by-side, highlighting differences +with colour-coded visualisations. ```bash python compare_surfdata.py [-o output_name] ``` **Arguments:** -- `file1.nc`, `file2.nc` - NetCDF files to compare -- `label1`, `label2` - Labels for legends and naming (e.g., "2005" "2006" or "control" "experiment") -- `-o output_name` - Optional custom output name (default: `comparison__vs_`) +- `file1.nc`, `file2.nc` – NetCDF files to compare (must be the same grid type) +- `label1`, `label2` – labels used in legends and output naming (e.g. `"2005" "2006"` or `"control" "experiment"`) +- `-o output_name` – optional custom output name (default: `comparison__vs_`) **Output:** -- `_figures/` - Directory containing PDF figures -- `.html` - Interactive HTML comparison report +- `_figures/` – directory containing one PDF per section +- `.html` – self-contained HTML comparison report + +**Colour coding:** Green = no change, Red = increased in file 2, Blue = decreased in file 2 -**Color coding:** Green = no change, Red = increased, Blue = decreased +**Mode-specific behaviour:** + +| Mode | Comparison approach | +|------|---------------------| +| Single-site | Scalar comparison tables, side-by-side bar charts, overlaid soil profiles, LAI time-series overlays | +| Regional grid | Difference maps (file 2 − file 1), domain-mean comparison bar charts, diverging-colourmap spatial plots | ## Requirements -See [`../requirements.txt`](../requirements.txt). `cartopy` is optional — site location maps are skipped gracefully if it is not installed. +Install dependencies from the local `requirements.txt`: + +```bash +pip install -r vizsurfdata/requirements.txt +``` + +| Package | Required | Purpose | +|---------|----------|---------| +| `numpy` | yes | array operations | +| `matplotlib` | yes | all plots | +| `netCDF4` | yes | reading `.nc` files | +| `cartopy` | optional | site / domain overview maps (skipped gracefully if absent) | + +## Notes + +- Both scripts write output next to the input NetCDF file. +- For regional grids with many active PFTs or CFTs the spatial-map figures can + be large; PDF rendering may take a minute or two. +- Domain-mean values shown in regional mode are averaged over land cells only + (determined from `LANDFRAC_PFT` or `PFTDATA_MASK`). diff --git a/vizsurfdata/compare_surfdata.py b/vizsurfdata/compare_surfdata.py index 9646407..3531d48 100644 --- a/vizsurfdata/compare_surfdata.py +++ b/vizsurfdata/compare_surfdata.py @@ -27,7 +27,9 @@ from visualize_surfdata import ( NATPFT_NAMES, CFT_NAMES, URBAN_TYPES, MONTH_NAMES, SOIL_DEPTHS, get_scalar_value, get_1d_array, fig_to_base64, - create_site_location_figure + create_site_location_figure, + detect_grid_type, get_field_2d, plot_map, plot_diff_map, plot_map_grid, + create_domain_map_figure, plot_monthly_timeseries_regional, ) @@ -560,140 +562,221 @@ def main(): nc1 = Dataset(args.file1, 'r') nc2 = Dataset(args.file2, 'r') - # Get coordinates for both sites - lon1 = get_scalar_value(nc1.variables['LONGXY']) - lat1 = get_scalar_value(nc1.variables['LATIXY']) - lon2 = get_scalar_value(nc2.variables['LONGXY']) - lat2 = get_scalar_value(nc2.variables['LATIXY']) + # Detect grid type (both files must match) + is_regional = detect_grid_type(nc1) + if detect_grid_type(nc2) != is_regional: + print("Warning: files have different grid types. Proceeding with " + f"{'regional' if is_regional else 'single-site'} mode (from file 1).") + print(f"Grid mode: {'regional' if is_regional else 'single-site'}") + + if is_regional: + lons1 = get_field_2d(nc1.variables['LONGXY']) + lats1 = get_field_2d(nc1.variables['LATIXY']) + lons2 = get_field_2d(nc2.variables['LONGXY']) + lats2 = get_field_2d(nc2.variables['LATIXY']) + # Use file-1 grid for plotting (assumed identical for same-domain comparisons) + lons = lons1 + lats = lats1 + # Land mask from file 1 + if 'LANDFRAC_PFT' in nc1.variables: + land_mask = get_field_2d(nc1.variables['LANDFRAC_PFT']) > 0 + elif 'PFTDATA_MASK' in nc1.variables: + land_mask = get_field_2d(nc1.variables['PFTDATA_MASK']).astype(bool) + else: + land_mask = np.ones(lons.shape, dtype=bool) + else: + lon1 = get_scalar_value(nc1.variables['LONGXY']) + lat1 = get_scalar_value(nc1.variables['LATIXY']) + lon2 = get_scalar_value(nc2.variables['LONGXY']) + lat2 = get_scalar_value(nc2.variables['LATIXY']) figures_data = [] # ========================================================================= - # SECTION 0: Site Locations Map + # SECTION 0: Domain / Site Locations Map # ========================================================================= - print("Creating Site Locations Map...") - site_map_figures = [] - fig_map = create_site_location_figure( - [lon1, lon2], [lat1, lat2], labels=[label1, label2]) - if fig_map is not None: - pdf_path = os.path.join(pdf_dir, '00_site_locations.pdf') - fig_map.savefig(pdf_path, bbox_inches='tight') - site_map_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': (f'Political map showing the locations of both sites: ' - f'{label1} ({lon1:.3f}°E, {lat1:.3f}°N) and ' - f'{label2} ({lon2:.3f}°E, {lat2:.3f}°N).'), - 'base64': fig_to_base64(fig_map) - }) - plt.close(fig_map) - - if site_map_figures: + print("Creating Section 0: Domain / Site Locations Map...") + section0_figures = [] + + if is_regional: + fig_map = create_domain_map_figure(lons, lats) + if fig_map is not None: + lon_min = float(np.nanmin(lons)) + lon_max = float(np.nanmax(lons)) + lat_min = float(np.nanmin(lats)) + lat_max = float(np.nanmax(lats)) + pdf_path = os.path.join(pdf_dir, '00_domain_overview.pdf') + fig_map.savefig(pdf_path, bbox_inches='tight') + section0_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': (f'Regional domain: {lon_min:.2f}\u00b0\u2013{lon_max:.2f}\u00b0E, ' + f'{lat_min:.2f}\u00b0\u2013{lat_max:.2f}\u00b0N ' + f'({lons.shape[1]}\u00d7{lons.shape[0]} cells). ' + f'Comparing {label1} vs {label2}.'), + 'base64': fig_to_base64(fig_map) + }) + plt.close(fig_map) + else: + fig_map = create_site_location_figure( + [lon1, lon2], [lat1, lat2], labels=[label1, label2]) + if fig_map is not None: + pdf_path = os.path.join(pdf_dir, '00_site_locations.pdf') + fig_map.savefig(pdf_path, bbox_inches='tight') + section0_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': (f'Site locations: {label1} ({lon1:.3f}\u00b0E, {lat1:.3f}\u00b0N) ' + f'and {label2} ({lon2:.3f}\u00b0E, {lat2:.3f}\u00b0N).'), + 'base64': fig_to_base64(fig_map) + }) + plt.close(fig_map) + + if section0_figures: figures_data.append({ - 'id': 'site-locations', - 'title': 'Site Locations', - 'description': 'Political map showing the geographic locations of both sites.', - 'figures': site_map_figures + 'id': 'locations', + 'title': 'Domain / Site Locations', + 'description': 'Geographic overview of the compared domain or sites.', + 'figures': section0_figures }) # ========================================================================= - # SECTION 1: Scalar Parameters Comparison + # SECTION 1: Basic Parameters Comparison # ========================================================================= - print("Creating Section 1: Scalar Parameters Comparison...") + print("Creating Section 1: Basic Parameters Comparison...") section1_figures = [] - fig, axes = plt.subplots(2, 2, figsize=(16, 12)) - - # Location and basic parameters - basic_params = [ - ('LONGXY', 'degrees E'), - ('LATIXY', 'degrees N'), - ('AREA', 'km^2'), - ('FMAX', 'unitless'), - ('zbedrock', 'm'), - ('SLOPE', 'degrees'), - ('STD_ELEV', 'm'), - ('gdp', 'unitless'), - ] - - basic_data = [] - for var_name, units in basic_params: - val1 = get_scalar_value(nc1.variables[var_name]) - val2 = get_scalar_value(nc2.variables[var_name]) - basic_data.append((var_name, val1, val2, units)) - - create_comparison_table(axes[0, 0], basic_data, 'Basic Site Parameters', label1, label2) - - # Land cover fractions - land_params = [ - ('PCT_NATVEG', '%'), - ('PCT_CROP', '%'), - ('PCT_LAKE', '%'), - ('PCT_WETLAND', '%'), - ('PCT_GLACIER', '%'), - ('peatf', 'fraction'), - ('LANDFRAC_PFT', 'fraction'), - ] - - land_data = [] - for var_name, units in land_params: - val1 = get_scalar_value(nc1.variables[var_name]) - val2 = get_scalar_value(nc2.variables[var_name]) - land_data.append((var_name, val1, val2, units)) - - # Add total urban - urban1 = np.sum(get_1d_array(nc1.variables['PCT_URBAN'])) - urban2 = np.sum(get_1d_array(nc2.variables['PCT_URBAN'])) - land_data.append(('PCT_URBAN (total)', urban1, urban2, '%')) - - create_comparison_table(axes[0, 1], land_data, 'Land Cover Fractions', label1, label2) - - # Soil parameters - soil_params = [ - ('SOIL_COLOR', 'index'), - ('mxsoil_color', 'count'), - ] - - soil_data = [] - for var_name, units in soil_params: - val1 = get_scalar_value(nc1.variables[var_name]) - val2 = get_scalar_value(nc2.variables[var_name]) - soil_data.append((var_name, val1, val2, units)) - - create_comparison_table(axes[1, 0], soil_data, 'Soil Parameters', label1, label2) - - # Emission factors - ef_params = [ - ('EF1_BTR', 'unitless'), - ('EF1_FET', 'unitless'), - ('EF1_FDT', 'unitless'), - ('EF1_SHR', 'unitless'), - ('EF1_GRS', 'unitless'), - ('EF1_CRP', 'unitless'), - ] - - ef_data = [] - for var_name, units in ef_params: - val1 = get_scalar_value(nc1.variables[var_name]) - val2 = get_scalar_value(nc2.variables[var_name]) - ef_data.append((var_name, val1, val2, units)) - - create_comparison_table(axes[1, 1], ef_data, 'Emission Factors', label1, label2) - - fig.suptitle('Scalar Parameters Comparison', fontsize=14, fontweight='bold') - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '01_scalar_comparison.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section1_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Comparison of scalar parameters: basic site info, land cover fractions, soil parameters, and emission factors.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + # Helper: domain mean of a 2-D field, masked to land cells + def _dmean(nc, var_name): + d = get_field_2d(nc.variables[var_name]) + return float(np.nanmean(np.where(land_mask, d, np.nan))) + + if is_regional: + # Difference maps for key spatial fields + basic_diff_vars = [ + ('AREA', 'km\u00b2 diff', 'RdBu_r'), + ('FMAX', 'unitless diff', 'RdBu_r'), + ('SOIL_COLOR','index diff', 'RdBu_r'), + ('zbedrock', 'm diff', 'RdBu_r'), + ('SLOPE', 'deg diff', 'RdBu_r'), + ('STD_ELEV', 'm diff', 'RdBu_r'), + ('peatf', 'fraction diff', 'RdBu_r'), + ('gdp', 'unitless diff', 'RdBu_r'), + ] + diff_data, diff_titles, diff_units = [], [], [] + for var_name, units_str, _ in basic_diff_vars: + if var_name in nc1.variables and var_name in nc2.variables: + d1 = get_field_2d(nc1.variables[var_name]) + d2 = get_field_2d(nc2.variables[var_name]) + if d1.ndim == 2 and d2.ndim == 2: + diff_data.append(d2 - d1) + diff_titles.append(f'\u0394 {var_name}') + diff_units.append(units_str) + + n = len(diff_data) + ncols = 3 + nrows = max(1, (n + ncols - 1) // ncols) + fig, axes = plt.subplots(nrows, ncols, + figsize=(ncols * 4.2, nrows * 3.2), squeeze=False) + axes_flat = axes.flatten() + for i in range(n): + plot_diff_map(axes_flat[i], diff_data[i], lons, lats, + diff_titles[i], diff_units[i]) + for i in range(n, len(axes_flat)): + axes_flat[i].axis('off') + fig.suptitle(f'Basic Parameters \u2013 Difference ({label2} \u2212 {label1})', + fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '01_basic_diff_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': f'Difference maps of basic parameters ({label2} \u2212 {label1}). Red = increase, Blue = decrease.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Domain-mean comparison table + scalar_table_params = [ + ('AREA', 'km\u00b2'), ('FMAX', 'unitless'), ('zbedrock', 'm'), + ('SLOPE', 'deg'), ('STD_ELEV', 'm'), ('peatf', 'fraction'), + ('PCT_NATVEG', '%'), ('PCT_CROP', '%'), ('PCT_LAKE', '%'), + ('PCT_WETLAND', '%'), ('PCT_GLACIER', '%'), ('gdp', 'unitless'), + ] + tbl_data = [] + for var_name, units_str in scalar_table_params: + if var_name in nc1.variables and var_name in nc2.variables: + v1 = _dmean(nc1, var_name) + v2 = _dmean(nc2, var_name) + tbl_data.append((f'{var_name} (domain mean)', v1, v2, units_str)) + urban1_mean = float(np.nanmean(np.where(land_mask, + np.array(nc1.variables['PCT_URBAN'][:], dtype=float).sum(axis=0), np.nan))) + urban2_mean = float(np.nanmean(np.where(land_mask, + np.array(nc2.variables['PCT_URBAN'][:], dtype=float).sum(axis=0), np.nan))) + tbl_data.append(('PCT_URBAN total (domain mean)', urban1_mean, urban2_mean, '%')) + + fig, ax = plt.subplots(figsize=(16, 8)) + create_comparison_table(ax, tbl_data, 'Domain-Mean Parameter Comparison', label1, label2) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '01b_domain_mean_table.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean comparison of key parameters.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + else: + fig, axes = plt.subplots(2, 2, figsize=(16, 12)) + + basic_params = [ + ('LONGXY', 'degrees E'), ('LATIXY', 'degrees N'), + ('AREA', 'km^2'), ('FMAX', 'unitless'), + ('zbedrock', 'm'), ('SLOPE', 'degrees'), + ('STD_ELEV', 'm'), ('gdp', 'unitless'), + ] + basic_data = [(v, get_scalar_value(nc1.variables[v]), + get_scalar_value(nc2.variables[v]), u) for v, u in basic_params] + create_comparison_table(axes[0, 0], basic_data, 'Basic Site Parameters', label1, label2) + + land_params = [ + ('PCT_NATVEG', '%'), ('PCT_CROP', '%'), ('PCT_LAKE', '%'), + ('PCT_WETLAND', '%'), ('PCT_GLACIER', '%'), + ('peatf', 'fraction'), ('LANDFRAC_PFT', 'fraction'), + ] + land_data = [(v, get_scalar_value(nc1.variables[v]), + get_scalar_value(nc2.variables[v]), u) for v, u in land_params] + land_data.append(('PCT_URBAN (total)', + float(np.sum(get_1d_array(nc1.variables['PCT_URBAN']))), + float(np.sum(get_1d_array(nc2.variables['PCT_URBAN']))), '%')) + create_comparison_table(axes[0, 1], land_data, 'Land Cover Fractions', label1, label2) + + soil_data = [(v, get_scalar_value(nc1.variables[v]), + get_scalar_value(nc2.variables[v]), u) + for v, u in [('SOIL_COLOR', 'index'), ('mxsoil_color', 'count')]] + create_comparison_table(axes[1, 0], soil_data, 'Soil Parameters', label1, label2) + + ef_params = [('EF1_BTR', 'unitless'), ('EF1_FET', 'unitless'), ('EF1_FDT', 'unitless'), + ('EF1_SHR', 'unitless'), ('EF1_GRS', 'unitless'), ('EF1_CRP', 'unitless')] + ef_data = [(v, get_scalar_value(nc1.variables[v]), + get_scalar_value(nc2.variables[v]), u) for v, u in ef_params] + create_comparison_table(axes[1, 1], ef_data, 'Emission Factors', label1, label2) + + fig.suptitle('Scalar Parameters Comparison', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '01_scalar_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Comparison of scalar parameters: basic site info, land cover fractions, soil parameters, and emission factors.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'scalar-comparison', - 'title': 'Scalar Parameters', - 'description': 'Comparison of single-value parameters. Green indicates no change, red indicates increase in File 2, blue indicates decrease.', + 'title': 'Basic Parameters Comparison', + 'description': ('Difference maps and domain-mean comparison table.' if is_regional else + 'Comparison of single-value parameters. Green = no change, red = increase, blue = decrease.'), 'figures': section1_figures }) @@ -703,58 +786,136 @@ def main(): print("Creating Section 2: Land Cover Comparison...") section2_figures = [] - fig, axes = plt.subplots(1, 2, figsize=(14, 6)) - - # Major land cover comparison land_labels = ['Natural Veg', 'Cropland', 'Urban', 'Lake', 'Wetland', 'Glacier'] - land_vals1 = [ - get_scalar_value(nc1.variables['PCT_NATVEG']), - get_scalar_value(nc1.variables['PCT_CROP']), - np.sum(get_1d_array(nc1.variables['PCT_URBAN'])), - get_scalar_value(nc1.variables['PCT_LAKE']), - get_scalar_value(nc1.variables['PCT_WETLAND']), - get_scalar_value(nc1.variables['PCT_GLACIER']) - ] - land_vals2 = [ - get_scalar_value(nc2.variables['PCT_NATVEG']), - get_scalar_value(nc2.variables['PCT_CROP']), - np.sum(get_1d_array(nc2.variables['PCT_URBAN'])), - get_scalar_value(nc2.variables['PCT_LAKE']), - get_scalar_value(nc2.variables['PCT_WETLAND']), - get_scalar_value(nc2.variables['PCT_GLACIER']) - ] - - plot_comparison_bars(axes[0], land_labels, land_vals1, land_vals2, - 'Major Land Cover Types', '%', label1, label2) - - # Difference bar chart - diffs = np.array(land_vals2) - np.array(land_vals1) - colors = [get_change_color(d, 0.01) for d in diffs] - axes[1].bar(land_labels, diffs, color=colors, edgecolor='#2d3748') - axes[1].axhline(y=0, color='black', linestyle='-', linewidth=0.5) - axes[1].set_ylabel('Difference (File2 - File1) [%]') - axes[1].set_title('Land Cover Change', fontsize=11, fontweight='bold') - axes[1].tick_params(axis='x', rotation=45) - axes[1].grid(axis='y', alpha=0.3) - - for i, (d, label) in enumerate(zip(diffs, land_labels)): - if abs(d) > 0.01: - axes[1].text(i, d + 0.1 * np.sign(d), f'{d:+.2f}', ha='center', fontsize=9) - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '02_land_cover_comparison.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section2_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Side-by-side comparison of major land cover types and their differences.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + if is_regional: + # Domain-mean land cover values (needed also in Section 8) + urban_sum1 = np.array(nc1.variables['PCT_URBAN'][:], dtype=float).sum(axis=0) + urban_sum2 = np.array(nc2.variables['PCT_URBAN'][:], dtype=float).sum(axis=0) + land_vals1 = [ + _dmean(nc1, 'PCT_NATVEG'), _dmean(nc1, 'PCT_CROP'), + float(np.nanmean(np.where(land_mask, urban_sum1, np.nan))), + _dmean(nc1, 'PCT_LAKE'), _dmean(nc1, 'PCT_WETLAND'), _dmean(nc1, 'PCT_GLACIER'), + ] + land_vals2 = [ + _dmean(nc2, 'PCT_NATVEG'), _dmean(nc2, 'PCT_CROP'), + float(np.nanmean(np.where(land_mask, urban_sum2, np.nan))), + _dmean(nc2, 'PCT_LAKE'), _dmean(nc2, 'PCT_WETLAND'), _dmean(nc2, 'PCT_GLACIER'), + ] + + # Difference maps + lc_diff_specs = [ + ('PCT_NATVEG', '%'), ('PCT_CROP', '%'), ('PCT_LAKE', '%'), + ('PCT_WETLAND', '%'), ('PCT_GLACIER', '%'), + ] + diff_data, diff_titles, diff_units_list = [], [], [] + for var_name, units_str in lc_diff_specs: + if var_name in nc1.variables and var_name in nc2.variables: + d1 = get_field_2d(nc1.variables[var_name]) + d2 = get_field_2d(nc2.variables[var_name]) + diff_data.append(d2 - d1) + diff_titles.append(f'\u0394 {var_name}') + diff_units_list.append(units_str) + diff_data.append(urban_sum2 - urban_sum1) + diff_titles.append('\u0394 PCT_URBAN (total)') + diff_units_list.append('%') + + ncols = 3 + n = len(diff_data) + nrows = max(1, (n + ncols - 1) // ncols) + fig, axes = plt.subplots(nrows, ncols, + figsize=(ncols * 4.2, nrows * 3.2), squeeze=False) + axes_flat = axes.flatten() + for i in range(n): + plot_diff_map(axes_flat[i], diff_data[i], lons, lats, + diff_titles[i], diff_units_list[i]) + for i in range(n, len(axes_flat)): + axes_flat[i].axis('off') + fig.suptitle(f'Land Cover \u2013 Difference ({label2} \u2212 {label1})', + fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '02_landcover_diff_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section2_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': (f'Difference maps of land cover fractions ({label2} \u2212 {label1}). ' + 'Red = increase, Blue = decrease.'), + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Domain-mean comparison bar chart + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + plot_comparison_bars(axes[0], land_labels, land_vals1, land_vals2, + 'Land Cover Domain Means', '%', label1, label2) + diffs = np.array(land_vals2) - np.array(land_vals1) + colors = [get_change_color(d, 0.01) for d in diffs] + axes[1].bar(land_labels, diffs, color=colors, edgecolor='#2d3748') + axes[1].axhline(y=0, color='black', linestyle='-', linewidth=0.5) + axes[1].set_ylabel('Difference (domain-mean) [%]') + axes[1].set_title('Land Cover Change (domain mean)', fontsize=11, fontweight='bold') + axes[1].tick_params(axis='x', rotation=45) + axes[1].grid(axis='y', alpha=0.3) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '02b_landcover_mean_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section2_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean land cover comparison and differences.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + else: + land_vals1 = [ + get_scalar_value(nc1.variables['PCT_NATVEG']), + get_scalar_value(nc1.variables['PCT_CROP']), + np.sum(get_1d_array(nc1.variables['PCT_URBAN'])), + get_scalar_value(nc1.variables['PCT_LAKE']), + get_scalar_value(nc1.variables['PCT_WETLAND']), + get_scalar_value(nc1.variables['PCT_GLACIER']) + ] + land_vals2 = [ + get_scalar_value(nc2.variables['PCT_NATVEG']), + get_scalar_value(nc2.variables['PCT_CROP']), + np.sum(get_1d_array(nc2.variables['PCT_URBAN'])), + get_scalar_value(nc2.variables['PCT_LAKE']), + get_scalar_value(nc2.variables['PCT_WETLAND']), + get_scalar_value(nc2.variables['PCT_GLACIER']) + ] + + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + plot_comparison_bars(axes[0], land_labels, land_vals1, land_vals2, + 'Major Land Cover Types', '%', label1, label2) + + diffs = np.array(land_vals2) - np.array(land_vals1) + colors = [get_change_color(d, 0.01) for d in diffs] + axes[1].bar(land_labels, diffs, color=colors, edgecolor='#2d3748') + axes[1].axhline(y=0, color='black', linestyle='-', linewidth=0.5) + axes[1].set_ylabel('Difference (File2 - File1) [%]') + axes[1].set_title('Land Cover Change', fontsize=11, fontweight='bold') + axes[1].tick_params(axis='x', rotation=45) + axes[1].grid(axis='y', alpha=0.3) + + for i, (d, lbl) in enumerate(zip(diffs, land_labels)): + if abs(d) > 0.01: + axes[1].text(i, d + 0.1 * np.sign(d), f'{d:+.2f}', ha='center', fontsize=9) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '02_land_cover_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section2_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Side-by-side comparison of major land cover types and their differences.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'land-cover-comparison', 'title': 'Land Cover Comparison', - 'description': 'Comparison of major land cover type fractions between the two files.', + 'description': ('Difference maps and domain-mean land cover comparison.' if is_regional else + 'Comparison of major land cover type fractions between the two files.'), 'figures': section2_figures }) @@ -764,95 +925,236 @@ def main(): print("Creating Section 3: Natural PFT Comparison...") section3_figures = [] - pct_nat_pft1 = get_1d_array(nc1.variables['PCT_NAT_PFT']) - pct_nat_pft2 = get_1d_array(nc2.variables['PCT_NAT_PFT']) natpft_labels = [NATPFT_NAMES.get(i, f'PFT {i}')[:25] for i in range(15)] - fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + if is_regional: + # Domain-mean PCT_NAT_PFT per PFT index (shape: natpft, lat, lon) + raw1 = np.array(nc1.variables['PCT_NAT_PFT'][:], dtype=float) + raw2 = np.array(nc2.variables['PCT_NAT_PFT'][:], dtype=float) + n_pft = raw1.shape[0] + pct_nat_pft1 = np.array([ + float(np.nanmean(np.where(land_mask, raw1[p], np.nan))) for p in range(n_pft) + ]) + pct_nat_pft2 = np.array([ + float(np.nanmean(np.where(land_mask, raw2[p], np.nan))) for p in range(n_pft) + ]) - # Side-by-side comparison - plot_comparison_bars(axes[0], natpft_labels, pct_nat_pft1, pct_nat_pft2, - 'Natural PFT Distribution', '%', label1, label2) - - # Difference plot - diffs = pct_nat_pft2 - pct_nat_pft1 - colors = [get_change_color(d, 0.1) for d in diffs] - y_pos = range(15) - axes[1].barh(y_pos, diffs, color=colors, edgecolor='#2d3748') - axes[1].axvline(x=0, color='black', linestyle='-', linewidth=0.5) - axes[1].set_yticks(y_pos) - axes[1].set_yticklabels(natpft_labels, fontsize=8) - axes[1].set_xlabel('Difference (File2 - File1) [%]') - axes[1].set_title('Natural PFT Change', fontsize=11, fontweight='bold') - axes[1].grid(axis='x', alpha=0.3) - - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '03_natural_pft_comparison.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section3_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Comparison of natural Plant Functional Type distributions.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) - - figures_data.append({ - 'id': 'natural-pft-comparison', - 'title': 'Natural PFT Comparison', - 'description': 'Comparison of natural Plant Functional Type (PFT) distributions within the natural vegetation landunit.', - 'figures': section3_figures - }) - - # ========================================================================= - # SECTION 4: Crop Functional Types Comparison - # ========================================================================= - print("Creating Section 4: CFT Comparison...") - section4_figures = [] - - pct_cft1 = get_1d_array(nc1.variables['PCT_CFT']) - pct_cft2 = get_1d_array(nc2.variables['PCT_CFT']) - cft_indices = nc1.variables['cft'][:] + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + plot_comparison_bars(axes[0], natpft_labels, pct_nat_pft1, pct_nat_pft2, + 'Natural PFT Domain Means', '%', label1, label2) + diffs = pct_nat_pft2 - pct_nat_pft1 + colors = [get_change_color(d, 0.1) for d in diffs] + y_pos = range(n_pft) + axes[1].barh(y_pos, diffs, color=colors, edgecolor='#2d3748') + axes[1].axvline(x=0, color='black', linestyle='-', linewidth=0.5) + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(natpft_labels, fontsize=8) + axes[1].set_xlabel('Difference domain mean (File2 \u2212 File1) [%]') + axes[1].set_title('Natural PFT Change (domain mean)', fontsize=11, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '03_natural_pft_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section3_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': ('Domain-mean natural PFT distribution comparison ' + f'({label2} \u2212 {label1}).'), + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - # Find CFTs that are non-zero in either file - nonzero_mask = (pct_cft1 > 0.1) | (pct_cft2 > 0.1) + # Spatial diff maps for the most changed PFTs + abs_diffs = np.abs(diffs) + top_pft_indices = np.argsort(abs_diffs)[::-1][:6] + top_pft_indices = [p for p in top_pft_indices if abs_diffs[p] > 0.01] + if top_pft_indices: + ncols = 3 + nrows = max(1, (len(top_pft_indices) + ncols - 1) // ncols) + fig, axes = plt.subplots(nrows, ncols, + figsize=(ncols * 4.2, nrows * 3.2), squeeze=False) + axes_flat = axes.flatten() + for idx, p in enumerate(top_pft_indices): + diff2d = raw2[p] - raw1[p] + plot_diff_map(axes_flat[idx], diff2d, lons, lats, + f'\u0394 {natpft_labels[p]}', '%') + for idx in range(len(top_pft_indices), len(axes_flat)): + axes_flat[idx].axis('off') + fig.suptitle(f'Natural PFT Difference Maps ({label2} \u2212 {label1})', + fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '03b_natural_pft_diff_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section3_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': ('Spatial difference maps for the most changed natural PFTs. ' + 'Red = increase, Blue = decrease.'), + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - if np.any(nonzero_mask): - nonzero_cft1 = pct_cft1[nonzero_mask] - nonzero_cft2 = pct_cft2[nonzero_mask] - nonzero_indices = cft_indices[nonzero_mask] - nonzero_labels = [CFT_NAMES.get(int(i), f'CFT {i}')[:25] for i in nonzero_indices] + else: + pct_nat_pft1 = get_1d_array(nc1.variables['PCT_NAT_PFT']) + pct_nat_pft2 = get_1d_array(nc2.variables['PCT_NAT_PFT']) fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + plot_comparison_bars(axes[0], natpft_labels, pct_nat_pft1, pct_nat_pft2, + 'Natural PFT Distribution', '%', label1, label2) - plot_comparison_bars(axes[0], nonzero_labels, nonzero_cft1, nonzero_cft2, - 'Crop Type Distribution', '%', label1, label2) - - # Difference plot - diffs = nonzero_cft2 - nonzero_cft1 + diffs = pct_nat_pft2 - pct_nat_pft1 colors = [get_change_color(d, 0.1) for d in diffs] - y_pos = range(len(nonzero_labels)) + y_pos = range(15) axes[1].barh(y_pos, diffs, color=colors, edgecolor='#2d3748') axes[1].axvline(x=0, color='black', linestyle='-', linewidth=0.5) axes[1].set_yticks(y_pos) - axes[1].set_yticklabels(nonzero_labels, fontsize=8) + axes[1].set_yticklabels(natpft_labels, fontsize=8) axes[1].set_xlabel('Difference (File2 - File1) [%]') - axes[1].set_title('Crop Type Change', fontsize=11, fontweight='bold') + axes[1].set_title('Natural PFT Change', fontsize=11, fontweight='bold') axes[1].grid(axis='x', alpha=0.3) plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '04_cft_comparison.pdf') + pdf_path = os.path.join(pdf_dir, '03_natural_pft_comparison.pdf') fig.savefig(pdf_path, bbox_inches='tight') - section4_figures.append({ + section3_figures.append({ 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Comparison of Crop Functional Type distributions (showing only non-zero CFTs).', + 'caption': 'Comparison of natural Plant Functional Type distributions.', 'base64': fig_to_base64(fig) }) plt.close(fig) + figures_data.append({ + 'id': 'natural-pft-comparison', + 'title': 'Natural PFT Comparison', + 'description': ('Domain-mean comparison and spatial difference maps of natural PFTs.' + if is_regional else + 'Comparison of natural Plant Functional Type (PFT) distributions ' + 'within the natural vegetation landunit.'), + 'figures': section3_figures + }) + + # ========================================================================= + # SECTION 4: Crop Functional Types Comparison + # ========================================================================= + print("Creating Section 4: CFT Comparison...") + section4_figures = [] + + cft_indices = nc1.variables['cft'][:] + + if is_regional: + raw_cft1 = np.array(nc1.variables['PCT_CFT'][:], dtype=float) + raw_cft2 = np.array(nc2.variables['PCT_CFT'][:], dtype=float) + n_cft = raw_cft1.shape[0] + pct_cft1 = np.array([ + float(np.nanmean(np.where(land_mask, raw_cft1[c], np.nan))) for c in range(n_cft) + ]) + pct_cft2 = np.array([ + float(np.nanmean(np.where(land_mask, raw_cft2[c], np.nan))) for c in range(n_cft) + ]) + nonzero_mask = (pct_cft1 > 0.1) | (pct_cft2 > 0.1) + + if np.any(nonzero_mask): + nonzero_cft1 = pct_cft1[nonzero_mask] + nonzero_cft2 = pct_cft2[nonzero_mask] + nonzero_indices = cft_indices[nonzero_mask] + nonzero_labels = [CFT_NAMES.get(int(i), f'CFT {i}')[:25] for i in nonzero_indices] + + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + plot_comparison_bars(axes[0], nonzero_labels, nonzero_cft1, nonzero_cft2, + 'Crop Type Domain Means', '%', label1, label2) + diffs_cft = nonzero_cft2 - nonzero_cft1 + colors = [get_change_color(d, 0.1) for d in diffs_cft] + y_pos = range(len(nonzero_labels)) + axes[1].barh(y_pos, diffs_cft, color=colors, edgecolor='#2d3748') + axes[1].axvline(x=0, color='black', linestyle='-', linewidth=0.5) + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(nonzero_labels, fontsize=8) + axes[1].set_xlabel('Difference domain mean (File2 \u2212 File1) [%]') + axes[1].set_title('Crop Type Change (domain mean)', fontsize=11, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04_cft_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean CFT distribution comparison (active CFTs only).', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Spatial diff maps for the most changed CFTs + nonzero_all = np.where(nonzero_mask)[0] + abs_diffs_cft = np.abs(diffs_cft) + order = np.argsort(abs_diffs_cft)[::-1] + top_cft_local = [i for i in order if abs_diffs_cft[i] > 0.01][:6] + top_cft_global = [nonzero_all[i] for i in top_cft_local] + top_cft_labels = [nonzero_labels[i] for i in top_cft_local] + + if top_cft_global: + ncols = 3 + nrows = max(1, (len(top_cft_global) + ncols - 1) // ncols) + fig, axes = plt.subplots(nrows, ncols, + figsize=(ncols * 4.2, nrows * 3.2), squeeze=False) + axes_flat = axes.flatten() + for idx, (g, lbl) in enumerate(zip(top_cft_global, top_cft_labels)): + diff2d = raw_cft2[g] - raw_cft1[g] + plot_diff_map(axes_flat[idx], diff2d, lons, lats, + f'\u0394 {lbl}', '%') + for idx in range(len(top_cft_global), len(axes_flat)): + axes_flat[idx].axis('off') + fig.suptitle(f'CFT Difference Maps ({label2} \u2212 {label1})', + fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04b_cft_diff_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': ('Spatial difference maps for the most changed CFTs. ' + 'Red = increase, Blue = decrease.'), + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + else: + pct_cft1 = get_1d_array(nc1.variables['PCT_CFT']) + pct_cft2 = get_1d_array(nc2.variables['PCT_CFT']) + nonzero_mask = (pct_cft1 > 0.1) | (pct_cft2 > 0.1) + + if np.any(nonzero_mask): + nonzero_cft1 = pct_cft1[nonzero_mask] + nonzero_cft2 = pct_cft2[nonzero_mask] + nonzero_indices = cft_indices[nonzero_mask] + nonzero_labels = [CFT_NAMES.get(int(i), f'CFT {i}')[:25] for i in nonzero_indices] + + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + plot_comparison_bars(axes[0], nonzero_labels, nonzero_cft1, nonzero_cft2, + 'Crop Type Distribution', '%', label1, label2) + + diffs_cft = nonzero_cft2 - nonzero_cft1 + colors = [get_change_color(d, 0.1) for d in diffs_cft] + y_pos = range(len(nonzero_labels)) + axes[1].barh(y_pos, diffs_cft, color=colors, edgecolor='#2d3748') + axes[1].axvline(x=0, color='black', linestyle='-', linewidth=0.5) + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(nonzero_labels, fontsize=8) + axes[1].set_xlabel('Difference (File2 - File1) [%]') + axes[1].set_title('Crop Type Change', fontsize=11, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04_cft_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Comparison of Crop Functional Type distributions (showing only non-zero CFTs).', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + figures_data.append({ 'id': 'cft-comparison', 'title': 'Crop Functional Types Comparison', - 'description': 'Comparison of Crop Functional Type (CFT) distributions within the crop landunit.', + 'description': ('Domain-mean CFT comparison and spatial difference maps.' + if is_regional else + 'Comparison of Crop Functional Type (CFT) distributions within the crop landunit.'), 'figures': section4_figures }) @@ -862,67 +1164,132 @@ def main(): print("Creating Section 5: Soil Properties Comparison...") section5_figures = [] - fig, axes = plt.subplots(1, 3, figsize=(15, 6)) - - # Sand - sand1 = get_1d_array(nc1.variables['PCT_SAND']) - sand2 = get_1d_array(nc2.variables['PCT_SAND']) - plot_difference_profile(axes[0], SOIL_DEPTHS, sand1, sand2, 'Sand Content', '%', label1, label2) - - # Clay - clay1 = get_1d_array(nc1.variables['PCT_CLAY']) - clay2 = get_1d_array(nc2.variables['PCT_CLAY']) - plot_difference_profile(axes[1], SOIL_DEPTHS, clay1, clay2, 'Clay Content', '%', label1, label2) + if is_regional: + # PCT_SAND, PCT_CLAY, ORGANIC shape: (nlevsoi, lsmlat, lsmlon) + raw_sand1 = np.array(nc1.variables['PCT_SAND'][:], dtype=float) + raw_sand2 = np.array(nc2.variables['PCT_SAND'][:], dtype=float) + raw_clay1 = np.array(nc1.variables['PCT_CLAY'][:], dtype=float) + raw_clay2 = np.array(nc2.variables['PCT_CLAY'][:], dtype=float) + raw_org1 = np.array(nc1.variables['ORGANIC'][:], dtype=float) + raw_org2 = np.array(nc2.variables['ORGANIC'][:], dtype=float) + n_lev = raw_sand1.shape[0] + + # Domain-mean profiles + sand1 = np.array([float(np.nanmean(np.where(land_mask, raw_sand1[l], np.nan))) + for l in range(n_lev)]) + sand2 = np.array([float(np.nanmean(np.where(land_mask, raw_sand2[l], np.nan))) + for l in range(n_lev)]) + clay1 = np.array([float(np.nanmean(np.where(land_mask, raw_clay1[l], np.nan))) + for l in range(n_lev)]) + clay2 = np.array([float(np.nanmean(np.where(land_mask, raw_clay2[l], np.nan))) + for l in range(n_lev)]) + org1 = np.array([float(np.nanmean(np.where(land_mask, raw_org1[l], np.nan))) + for l in range(n_lev)]) + org2 = np.array([float(np.nanmean(np.where(land_mask, raw_org2[l], np.nan))) + for l in range(n_lev)]) + depths = np.array(SOIL_DEPTHS[:n_lev]) + + # Domain-mean profile comparison + fig, axes = plt.subplots(1, 3, figsize=(15, 6)) + plot_difference_profile(axes[0], depths, sand1, sand2, 'Sand Content (domain mean)', '%', label1, label2) + plot_difference_profile(axes[1], depths, clay1, clay2, 'Clay Content (domain mean)', '%', label1, label2) + plot_difference_profile(axes[2], depths, org1, org2, 'Organic Matter (domain mean)', 'kg/m\u00b3', label1, label2) + fig.suptitle('Soil Properties \u2013 Domain-Mean Profiles', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05_soil_profiles.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean soil profiles with shaded difference areas.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - # Organic - org1 = get_1d_array(nc1.variables['ORGANIC']) - org2 = get_1d_array(nc2.variables['ORGANIC']) - plot_difference_profile(axes[2], SOIL_DEPTHS, org1, org2, 'Organic Matter', 'kg/m3', label1, label2) + # Depth-mean difference maps + sand_dm1 = np.nanmean(raw_sand1, axis=0) + sand_dm2 = np.nanmean(raw_sand2, axis=0) + clay_dm1 = np.nanmean(raw_clay1, axis=0) + clay_dm2 = np.nanmean(raw_clay2, axis=0) + org_dm1 = np.nanmean(raw_org1, axis=0) + org_dm2 = np.nanmean(raw_org2, axis=0) + + fig, axes = plt.subplots(1, 3, figsize=(14, 4), squeeze=False) + axes_flat = axes.flatten() + plot_diff_map(axes_flat[0], sand_dm2 - sand_dm1, lons, lats, '\u0394 Sand (depth mean)', '%') + plot_diff_map(axes_flat[1], clay_dm2 - clay_dm1, lons, lats, '\u0394 Clay (depth mean)', '%') + plot_diff_map(axes_flat[2], org_dm2 - org_dm1, lons, lats, '\u0394 Organic (depth mean)', 'kg/m\u00b3') + fig.suptitle(f'Soil Difference Maps ({label2} \u2212 {label1})', + fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05b_soil_diff_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': ('Depth-mean difference maps for soil properties. ' + 'Red = increase, Blue = decrease.'), + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - fig.suptitle('Soil Properties Comparison by Depth', fontsize=14, fontweight='bold') - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '05_soil_comparison.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section5_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Vertical profiles of soil properties with shaded areas showing differences between files.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + else: + sand1 = get_1d_array(nc1.variables['PCT_SAND']) + sand2 = get_1d_array(nc2.variables['PCT_SAND']) + clay1 = get_1d_array(nc1.variables['PCT_CLAY']) + clay2 = get_1d_array(nc2.variables['PCT_CLAY']) + org1 = get_1d_array(nc1.variables['ORGANIC']) + org2 = get_1d_array(nc2.variables['ORGANIC']) + depths = np.array(SOIL_DEPTHS) + + fig, axes = plt.subplots(1, 3, figsize=(15, 6)) + plot_difference_profile(axes[0], depths, sand1, sand2, 'Sand Content', '%', label1, label2) + plot_difference_profile(axes[1], depths, clay1, clay2, 'Clay Content', '%', label1, label2) + plot_difference_profile(axes[2], depths, org1, org2, 'Organic Matter', 'kg/m3', label1, label2) + + fig.suptitle('Soil Properties Comparison by Depth', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05_soil_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Vertical profiles of soil properties with shaded areas showing differences between files.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - # Soil difference table - fig, ax = plt.subplots(figsize=(14, 8)) + # Soil difference table + fig, ax = plt.subplots(figsize=(14, 8)) - soil_layer_data = [] - for i, depth in enumerate(SOIL_DEPTHS): - soil_layer_data.append((f'Sand @ {depth:.2f}m', sand1[i], sand2[i], '%')) - soil_layer_data.append((f'Clay @ {depth:.2f}m', clay1[i], clay2[i], '%')) - soil_layer_data.append((f'Organic @ {depth:.2f}m', org1[i], org2[i], 'kg/m3')) + soil_layer_data = [] + for i, depth in enumerate(SOIL_DEPTHS): + soil_layer_data.append((f'Sand @ {depth:.2f}m', sand1[i], sand2[i], '%')) + soil_layer_data.append((f'Clay @ {depth:.2f}m', clay1[i], clay2[i], '%')) + soil_layer_data.append((f'Organic @ {depth:.2f}m', org1[i], org2[i], 'kg/m3')) - # Only show layers with differences - diff_layers = [(n, v1, v2, u) for n, v1, v2, u in soil_layer_data if abs(v2 - v1) > 0.01] + diff_layers = [(n, v1, v2, u) for n, v1, v2, u in soil_layer_data if abs(v2 - v1) > 0.01] - if diff_layers: - create_comparison_table(ax, diff_layers[:15], 'Soil Properties with Differences', label1, label2) - else: - ax.text(0.5, 0.5, 'No significant differences in soil properties', - ha='center', va='center', transform=ax.transAxes, fontsize=14) - ax.axis('off') + if diff_layers: + create_comparison_table(ax, diff_layers[:15], 'Soil Properties with Differences', label1, label2) + else: + ax.text(0.5, 0.5, 'No significant differences in soil properties', + ha='center', va='center', transform=ax.transAxes, fontsize=14) + ax.axis('off') - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '05b_soil_diff_table.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section5_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Table of soil properties showing layers with significant differences.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05b_soil_diff_table.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Table of soil properties showing layers with significant differences.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'soil-comparison', 'title': 'Soil Properties Comparison', - 'description': 'Comparison of soil texture (sand, clay) and organic matter content across depth layers.', + 'description': ('Domain-mean soil profiles and depth-mean difference maps.' + if is_regional else + 'Comparison of soil texture (sand, clay) and organic matter content across depth layers.'), 'figures': section5_figures }) @@ -932,73 +1299,156 @@ def main(): print("Creating Section 6: Monthly LAI Comparison...") section6_figures = [] - lai1 = nc1.variables['MONTHLY_LAI'][:] - lai2 = nc2.variables['MONTHLY_LAI'][:] + lai1 = np.array(nc1.variables['MONTHLY_LAI'][:], dtype=float) + lai2 = np.array(nc2.variables['MONTHLY_LAI'][:], dtype=float) - # Combine PFT names all_pft_names = {**NATPFT_NAMES, **CFT_NAMES} - - # Find PFTs with non-zero LAI in either file - active_pfts = [] - for pft in range(79): - if np.any(lai1[:, pft, 0, 0] > 0) or np.any(lai2[:, pft, 0, 0] > 0): - active_pfts.append(pft) - - # Plot comparison for top active PFTs - n_plots = min(6, len(active_pfts)) - if n_plots > 0: - fig, axes = plt.subplots(2, 3, figsize=(15, 10)) - axes = axes.flatten() - - for i, pft_idx in enumerate(active_pfts[:n_plots]): - pft_name = all_pft_names.get(pft_idx, f'PFT {pft_idx}')[:30] - plot_monthly_comparison(axes[i], lai1, lai2, pft_idx, 'LAI', 'm2/m2', label1, label2, pft_name) - - # Hide unused axes - for i in range(n_plots, 6): - axes[i].axis('off') - - fig.suptitle('Monthly LAI Comparison by PFT', fontsize=14, fontweight='bold') + n_pft_lai = lai1.shape[1] + + if is_regional: + # Domain-mean LAI per PFT per month: (12, n_pft) + dm_lai1 = np.array([ + [float(np.nanmean(np.where(land_mask, lai1[m, p], np.nan))) + for p in range(n_pft_lai)] + for m in range(12) + ]) + dm_lai2 = np.array([ + [float(np.nanmean(np.where(land_mask, lai2[m, p], np.nan))) + for p in range(n_pft_lai)] + for m in range(12) + ]) + + # Active PFTs (non-zero domain-mean in either file) + active_pfts = [p for p in range(n_pft_lai) + if np.any(dm_lai1[:, p] > 0) or np.any(dm_lai2[:, p] > 0)] + + # Time-series comparison for up to 6 active PFTs + n_plots = min(6, len(active_pfts)) + if n_plots > 0: + fig, axes_grid = plt.subplots(2, 3, figsize=(15, 10)) + axes_grid = axes_grid.flatten() + months = range(12) + for i, pft_idx in enumerate(active_pfts[:n_plots]): + pft_name = all_pft_names.get(pft_idx, f'PFT {pft_idx}')[:30] + ax = axes_grid[i] + ax.plot(months, dm_lai1[:, pft_idx], 'o-', color='#4299e1', + linewidth=2, markersize=6, label=label1) + ax.plot(months, dm_lai2[:, pft_idx], 's-', color='#ed8936', + linewidth=2, markersize=6, label=label2) + ax.fill_between(months, dm_lai1[:, pft_idx], dm_lai2[:, pft_idx], + alpha=0.3, color='#a0aec0') + ax.set_xticks(months) + ax.set_xticklabels(MONTH_NAMES, rotation=45, ha='right', fontsize=7) + ax.set_ylabel('[m\u00b2/m\u00b2]') + ax.set_title(f'LAI \u2013 {pft_name}', fontsize=9, fontweight='bold') + ax.legend(fontsize=7) + ax.grid(alpha=0.3) + for i in range(n_plots, 6): + axes_grid[i].axis('off') + fig.suptitle('Monthly LAI Comparison (domain mean) by PFT', + fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06_lai_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': ('Domain-mean monthly LAI comparison for active PFTs. ' + 'Shaded area shows the difference.'), + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Domain-mean LAI difference heatmap + fig, ax = plt.subplots(figsize=(14, 8)) + if active_pfts: + shown = active_pfts[:15] + lai_diff_matrix = dm_lai2[:, shown].T - dm_lai1[:, shown].T # (n_shown, 12) + vmax = float(np.nanmax(np.abs(lai_diff_matrix))) or 1.0 + im = ax.imshow(lai_diff_matrix, cmap='RdBu_r', aspect='auto', + vmin=-vmax, vmax=vmax) + ax.set_xticks(range(12)) + ax.set_xticklabels(MONTH_NAMES) + ax.set_yticks(range(len(shown))) + ax.set_yticklabels([all_pft_names.get(p, f'PFT {p}')[:25] for p in shown], fontsize=8) + ax.set_title(f'LAI Difference (domain mean, {label2} \u2212 {label1})', + fontsize=12, fontweight='bold') + plt.colorbar(im, ax=ax, label='LAI diff (m\u00b2/m\u00b2)') plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '06_lai_comparison.pdf') + pdf_path = os.path.join(pdf_dir, '06b_lai_diff_heatmap.pdf') fig.savefig(pdf_path, bbox_inches='tight') section6_figures.append({ 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Monthly Leaf Area Index (LAI) comparison for active PFTs. Shaded area shows the difference.', + 'caption': ('Heatmap of domain-mean LAI differences by PFT and month. ' + 'Blue = decrease, Red = increase in File 2.'), 'base64': fig_to_base64(fig) }) plt.close(fig) - # LAI difference heatmap - fig, ax = plt.subplots(figsize=(14, 8)) - - if len(active_pfts) > 0: - lai_diff_matrix = np.zeros((len(active_pfts[:15]), 12)) - for i, pft in enumerate(active_pfts[:15]): - lai_diff_matrix[i, :] = lai2[:, pft, 0, 0] - lai1[:, pft, 0, 0] - - im = ax.imshow(lai_diff_matrix, cmap='RdBu_r', aspect='auto', vmin=-np.max(np.abs(lai_diff_matrix)), vmax=np.max(np.abs(lai_diff_matrix))) - ax.set_xticks(range(12)) - ax.set_xticklabels(MONTH_NAMES) - ax.set_yticks(range(len(active_pfts[:15]))) - ax.set_yticklabels([all_pft_names.get(p, f'PFT {p}')[:25] for p in active_pfts[:15]], fontsize=8) - ax.set_title('LAI Difference (File2 - File1)', fontsize=12, fontweight='bold') - plt.colorbar(im, ax=ax, label='LAI difference (m2/m2)') + else: + # Active PFTs (single-site: index [0, 0]) + active_pfts = [] + for pft in range(n_pft_lai): + if np.any(lai1[:, pft, 0, 0] > 0) or np.any(lai2[:, pft, 0, 0] > 0): + active_pfts.append(pft) + + n_plots = min(6, len(active_pfts)) + if n_plots > 0: + fig, axes_grid = plt.subplots(2, 3, figsize=(15, 10)) + axes_grid = axes_grid.flatten() + + for i, pft_idx in enumerate(active_pfts[:n_plots]): + pft_name = all_pft_names.get(pft_idx, f'PFT {pft_idx}')[:30] + plot_monthly_comparison(axes_grid[i], lai1, lai2, pft_idx, + 'LAI', 'm2/m2', label1, label2, pft_name) + + for i in range(n_plots, 6): + axes_grid[i].axis('off') + + fig.suptitle('Monthly LAI Comparison by PFT', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06_lai_comparison.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Monthly Leaf Area Index (LAI) comparison for active PFTs. Shaded area shows the difference.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # LAI difference heatmap + fig, ax = plt.subplots(figsize=(14, 8)) + if active_pfts: + shown = active_pfts[:15] + lai_diff_matrix = np.zeros((len(shown), 12)) + for i, pft in enumerate(shown): + lai_diff_matrix[i, :] = lai2[:, pft, 0, 0] - lai1[:, pft, 0, 0] + + vmax = float(np.nanmax(np.abs(lai_diff_matrix))) or 1.0 + im = ax.imshow(lai_diff_matrix, cmap='RdBu_r', aspect='auto', + vmin=-vmax, vmax=vmax) + ax.set_xticks(range(12)) + ax.set_xticklabels(MONTH_NAMES) + ax.set_yticks(range(len(shown))) + ax.set_yticklabels([all_pft_names.get(p, f'PFT {p}')[:25] for p in shown], fontsize=8) + ax.set_title('LAI Difference (File2 - File1)', fontsize=12, fontweight='bold') + plt.colorbar(im, ax=ax, label='LAI difference (m2/m2)') - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '06b_lai_diff_heatmap.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section6_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Heatmap of LAI differences by PFT and month. Blue = decrease, Red = increase in File 2.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06b_lai_diff_heatmap.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Heatmap of LAI differences by PFT and month. Blue = decrease, Red = increase in File 2.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'lai-comparison', 'title': 'Monthly LAI Comparison', - 'description': 'Comparison of monthly Leaf Area Index values for active Plant Functional Types.', + 'description': ('Domain-mean monthly LAI timeseries comparison and difference heatmap.' + if is_regional else + 'Comparison of monthly Leaf Area Index values for active Plant Functional Types.'), 'figures': section6_figures }) @@ -1008,43 +1458,59 @@ def main(): print("Creating Section 7: Urban Parameters Comparison...") section7_figures = [] - fig, axes = plt.subplots(2, 3, figsize=(15, 10)) + def _urban_dmean(nc, var_name, urban_idx): + """Domain-mean of a (numurbl, lat, lon) urban variable for one urban class.""" + d = np.array(nc.variables[var_name][:], dtype=float) + return float(np.nanmean(np.where(land_mask, d[urban_idx], np.nan))) + + n_urban = len(URBAN_TYPES) + + if is_regional: + pct_urban1 = np.array([_urban_dmean(nc1, 'PCT_URBAN', i) for i in range(n_urban)]) + pct_urban2 = np.array([_urban_dmean(nc2, 'PCT_URBAN', i) for i in range(n_urban)]) + hwr1 = np.array([_urban_dmean(nc1, 'CANYON_HWR', i) for i in range(n_urban)]) + hwr2 = np.array([_urban_dmean(nc2, 'CANYON_HWR', i) for i in range(n_urban)]) + ht1 = np.array([_urban_dmean(nc1, 'HT_ROOF', i) for i in range(n_urban)]) + ht2 = np.array([_urban_dmean(nc2, 'HT_ROOF', i) for i in range(n_urban)]) + tb1 = np.array([_urban_dmean(nc1, 'T_BUILDING_MIN', i) for i in range(n_urban)]) - 273.15 + tb2 = np.array([_urban_dmean(nc2, 'T_BUILDING_MIN', i) for i in range(n_urban)]) - 273.15 + rf1 = np.array([_urban_dmean(nc1, 'WTLUNIT_ROOF', i) for i in range(n_urban)]) + rf2 = np.array([_urban_dmean(nc2, 'WTLUNIT_ROOF', i) for i in range(n_urban)]) + pr1 = np.array([_urban_dmean(nc1, 'WTROAD_PERV', i) for i in range(n_urban)]) + pr2 = np.array([_urban_dmean(nc2, 'WTROAD_PERV', i) for i in range(n_urban)]) + else: + pct_urban1 = get_1d_array(nc1.variables['PCT_URBAN']) + pct_urban2 = get_1d_array(nc2.variables['PCT_URBAN']) + hwr1 = get_1d_array(nc1.variables['CANYON_HWR']) + hwr2 = get_1d_array(nc2.variables['CANYON_HWR']) + ht1 = get_1d_array(nc1.variables['HT_ROOF']) + ht2 = get_1d_array(nc2.variables['HT_ROOF']) + tb1 = get_1d_array(nc1.variables['T_BUILDING_MIN']) - 273.15 + tb2 = get_1d_array(nc2.variables['T_BUILDING_MIN']) - 273.15 + rf1 = get_1d_array(nc1.variables['WTLUNIT_ROOF']) + rf2 = get_1d_array(nc2.variables['WTLUNIT_ROOF']) + pr1 = get_1d_array(nc1.variables['WTROAD_PERV']) + pr2 = get_1d_array(nc2.variables['WTROAD_PERV']) - # Urban percent - pct_urban1 = get_1d_array(nc1.variables['PCT_URBAN']) - pct_urban2 = get_1d_array(nc2.variables['PCT_URBAN']) + fig, axes = plt.subplots(2, 3, figsize=(15, 10)) plot_comparison_bars(axes[0, 0], URBAN_TYPES, pct_urban1, pct_urban2, - 'Urban Fraction', '%', label1, label2) - - # Canyon HWR - hwr1 = get_1d_array(nc1.variables['CANYON_HWR']) - hwr2 = get_1d_array(nc2.variables['CANYON_HWR']) + 'Urban Fraction' + (' (domain mean)' if is_regional else ''), + '%', label1, label2) plot_comparison_bars(axes[0, 1], URBAN_TYPES, hwr1, hwr2, - 'Canyon H/W Ratio', 'ratio', label1, label2) - - # Roof height - ht1 = get_1d_array(nc1.variables['HT_ROOF']) - ht2 = get_1d_array(nc2.variables['HT_ROOF']) + 'Canyon H/W Ratio' + (' (domain mean)' if is_regional else ''), + 'ratio', label1, label2) plot_comparison_bars(axes[0, 2], URBAN_TYPES, ht1, ht2, - 'Roof Height', 'm', label1, label2) - - # Building temperature - tb1 = get_1d_array(nc1.variables['T_BUILDING_MIN']) - 273.15 - tb2 = get_1d_array(nc2.variables['T_BUILDING_MIN']) - 273.15 + 'Roof Height' + (' (domain mean)' if is_regional else ''), + 'm', label1, label2) plot_comparison_bars(axes[1, 0], URBAN_TYPES, tb1, tb2, - 'Min Building Temp', 'C', label1, label2) - - # Roof fraction - rf1 = get_1d_array(nc1.variables['WTLUNIT_ROOF']) - rf2 = get_1d_array(nc2.variables['WTLUNIT_ROOF']) + 'Min Building Temp' + (' (domain mean)' if is_regional else ''), + '\u00b0C', label1, label2) plot_comparison_bars(axes[1, 1], URBAN_TYPES, rf1, rf2, - 'Roof Fraction', 'fraction', label1, label2) - - # Pervious road fraction - pr1 = get_1d_array(nc1.variables['WTROAD_PERV']) - pr2 = get_1d_array(nc2.variables['WTROAD_PERV']) + 'Roof Fraction' + (' (domain mean)' if is_regional else ''), + 'fraction', label1, label2) plot_comparison_bars(axes[1, 2], URBAN_TYPES, pr1, pr2, - 'Pervious Road Fraction', 'fraction', label1, label2) + 'Pervious Road Fraction' + (' (domain mean)' if is_regional else ''), + 'fraction', label1, label2) fig.suptitle('Urban Parameters Comparison', fontsize=14, fontweight='bold') plt.tight_layout() @@ -1052,7 +1518,9 @@ def main(): fig.savefig(pdf_path, bbox_inches='tight') section7_figures.append({ 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Comparison of urban parameters across three density classes.', + 'caption': ('Domain-mean urban parameter comparison across three density classes.' + if is_regional else + 'Comparison of urban parameters across three density classes.'), 'base64': fig_to_base64(fig) }) plt.close(fig) @@ -1060,7 +1528,9 @@ def main(): figures_data.append({ 'id': 'urban-comparison', 'title': 'Urban Parameters Comparison', - 'description': 'Comparison of urban canyon and building parameters for three density classes.', + 'description': ('Domain-mean urban canyon and building parameter comparison.' + if is_regional else + 'Comparison of urban canyon and building parameters for three density classes.'), 'figures': section7_figures }) @@ -1070,33 +1540,32 @@ def main(): print("Creating Section 8: Summary of Differences...") section8_figures = [] - fig = plt.figure(figsize=(16, 12)) - - # Count differences + # Key variables compared: use domain means for regional, scalar values for single-site + summary_vars = ['AREA', 'FMAX', 'zbedrock', 'SLOPE', 'STD_ELEV', + 'PCT_NATVEG', 'PCT_CROP', 'PCT_LAKE', 'PCT_WETLAND', 'PCT_GLACIER', + 'SOIL_COLOR', 'peatf', 'gdp', 'LAKEDEPTH'] all_diffs = [] - - # Scalar differences - scalar_vars = ['LONGXY', 'LATIXY', 'AREA', 'FMAX', 'zbedrock', 'SLOPE', 'STD_ELEV', - 'PCT_NATVEG', 'PCT_CROP', 'PCT_LAKE', 'PCT_WETLAND', 'PCT_GLACIER', - 'SOIL_COLOR', 'peatf', 'gdp', 'LAKEDEPTH'] - - for var in scalar_vars: + for var in summary_vars: try: - v1 = get_scalar_value(nc1.variables[var]) - v2 = get_scalar_value(nc2.variables[var]) + if is_regional: + v1 = _dmean(nc1, var) + v2 = _dmean(nc2, var) + else: + v1 = get_scalar_value(nc1.variables[var]) + v2 = get_scalar_value(nc2.variables[var]) diff = abs(v2 - v1) if diff > 0.001: all_diffs.append((var, v1, v2, diff)) - except: + except Exception: pass - # Create summary visualization - ax1 = fig.add_subplot(2, 2, 1) + fig = plt.figure(figsize=(16, 12)) - n_total = len(scalar_vars) + 15 + 64 + 30 + 20 # approximate total variables + # Pie: changed vs unchanged + ax1 = fig.add_subplot(2, 2, 1) + n_total = len(summary_vars) + len(pct_nat_pft1) + 30 # approximate n_changed = len(all_diffs) - n_same = n_total - n_changed - + n_same = max(0, n_total - n_changed) ax1.pie([n_same, n_changed], labels=['Unchanged', 'Changed'], colors=['#48bb78', '#e53e3e'], autopct='%1.1f%%', startangle=90) ax1.set_title('Overall Parameter Changes', fontsize=12, fontweight='bold') @@ -1104,23 +1573,23 @@ def main(): # Top differences table ax2 = fig.add_subplot(2, 2, 2) ax2.axis('off') - if all_diffs: - # Sort by absolute difference sorted_diffs = sorted(all_diffs, key=lambda x: x[3], reverse=True)[:10] top_diff_data = [(name, v1, v2, 'various') for name, v1, v2, _ in sorted_diffs] - create_comparison_table(ax2, top_diff_data, 'Top 10 Scalar Differences', label1, label2) + create_comparison_table(ax2, top_diff_data, + 'Top Differences (domain mean)' if is_regional else 'Top 10 Scalar Differences', + label1, label2) else: - ax2.text(0.5, 0.5, 'No significant scalar differences found', - ha='center', va='center', fontsize=14) + ax2.text(0.5, 0.5, 'No significant differences found', + ha='center', va='center', fontsize=14) # Land cover change summary ax3 = fig.add_subplot(2, 2, 3) land_changes = np.array(land_vals2) - np.array(land_vals1) colors = [get_change_color(d, 0.01) for d in land_changes] - bars = ax3.barh(land_labels, land_changes, color=colors, edgecolor='#2d3748') + ax3.barh(land_labels, land_changes, color=colors, edgecolor='#2d3748') ax3.axvline(x=0, color='black', linestyle='-', linewidth=0.5) - ax3.set_xlabel('Change in Coverage (%)') + ax3.set_xlabel('Change in Coverage (%)' + (' [domain mean]' if is_regional else '')) ax3.set_title('Land Cover Changes Summary', fontsize=12, fontweight='bold') ax3.grid(axis='x', alpha=0.3) @@ -1128,21 +1597,20 @@ def main(): ax4 = fig.add_subplot(2, 2, 4) pft_diffs = pct_nat_pft2 - pct_nat_pft1 significant_pft_changes = [(NATPFT_NAMES.get(i, f'PFT {i}')[:20], pft_diffs[i]) - for i in range(15) if abs(pft_diffs[i]) > 0.1] - + for i in range(len(pft_diffs)) if abs(pft_diffs[i]) > 0.1] if significant_pft_changes: - labels, values = zip(*significant_pft_changes) + lbls, values = zip(*significant_pft_changes) colors = [get_change_color(v, 0.1) for v in values] - ax4.barh(range(len(labels)), values, color=colors, edgecolor='#2d3748') - ax4.set_yticks(range(len(labels))) - ax4.set_yticklabels(labels, fontsize=9) + ax4.barh(range(len(lbls)), values, color=colors, edgecolor='#2d3748') + ax4.set_yticks(range(len(lbls))) + ax4.set_yticklabels(lbls, fontsize=9) ax4.axvline(x=0, color='black', linestyle='-', linewidth=0.5) - ax4.set_xlabel('Change in PFT Coverage (%)') + ax4.set_xlabel('Change in PFT Coverage (%)' + (' [domain mean]' if is_regional else '')) ax4.set_title('Significant PFT Changes', fontsize=12, fontweight='bold') ax4.grid(axis='x', alpha=0.3) else: ax4.text(0.5, 0.5, 'No significant PFT changes', - ha='center', va='center', fontsize=14) + ha='center', va='center', fontsize=14) ax4.axis('off') fig.suptitle('Summary of All Differences', fontsize=14, fontweight='bold') diff --git a/vizsurfdata/requirements.txt b/vizsurfdata/requirements.txt index 573ce58..69dfe88 100644 --- a/vizsurfdata/requirements.txt +++ b/vizsurfdata/requirements.txt @@ -1,4 +1,4 @@ numpy matplotlib netCDF4 -cartopy +cartopy # optional – required for domain/site location maps; skipped gracefully if absent diff --git a/vizsurfdata/vizualize_surfdata.py b/vizsurfdata/vizualize_surfdata.py index ff31d28..4673231 100644 --- a/vizsurfdata/vizualize_surfdata.py +++ b/vizsurfdata/vizualize_surfdata.py @@ -202,6 +202,160 @@ def get_1d_array(var): return np.squeeze(data) +# ============================================================================= +# Grid-type detection +# ============================================================================= + +def detect_grid_type(nc): + """Return True if the file contains a regional (multi-cell) grid.""" + return nc.dimensions['lsmlat'].size > 1 or nc.dimensions['lsmlon'].size > 1 + + +# ============================================================================= +# Regional data helpers +# ============================================================================= + +def get_field_2d(var): + """Return a 2-D (nlat, nlon) float array for a (lsmlat, lsmlon) variable. + Any leading singleton dimensions are squeezed away.""" + data = var[:] + if hasattr(data, 'mask'): + data = np.ma.filled(data, np.nan) + return np.array(np.squeeze(data), dtype=float) + + +def _robust_clim(data): + """Return (vmin, vmax) clipped to the 2nd-98th percentile of finite values.""" + finite = data[np.isfinite(data)] + if finite.size == 0: + return 0.0, 1.0 + vmin = float(np.percentile(finite, 2)) + vmax = float(np.percentile(finite, 98)) + if vmin == vmax: + vmin -= 0.5 + vmax += 0.5 + return vmin, vmax + + +def plot_map(ax, data2d, lons2d, lats2d, title, units, cmap='viridis'): + """Plot a 2-D spatial field using pcolormesh with a colourbar.""" + vmin, vmax = _robust_clim(data2d) + pcm = ax.pcolormesh(lons2d, lats2d, data2d, cmap=cmap, + vmin=vmin, vmax=vmax, shading='auto') + plt.colorbar(pcm, ax=ax, label=units, fraction=0.046, pad=0.04) + ax.set_title(title, fontsize=10, fontweight='bold') + ax.set_xlabel('Lon (\u00b0E)', fontsize=7) + ax.set_ylabel('Lat (\u00b0N)', fontsize=7) + ax.tick_params(labelsize=7) + + +def plot_diff_map(ax, diff2d, lons2d, lats2d, title, units): + """Plot a difference field (file2 - file1) with a symmetric diverging colourmap.""" + finite = diff2d[np.isfinite(diff2d)] + vmax = float(np.percentile(np.abs(finite), 98)) if finite.size else 1.0 + if vmax == 0: + vmax = 1.0 + pcm = ax.pcolormesh(lons2d, lats2d, diff2d, cmap='RdBu_r', + vmin=-vmax, vmax=vmax, shading='auto') + plt.colorbar(pcm, ax=ax, label=units, fraction=0.046, pad=0.04) + ax.set_title(title, fontsize=10, fontweight='bold') + ax.set_xlabel('Lon (\u00b0E)', fontsize=7) + ax.set_ylabel('Lat (\u00b0N)', fontsize=7) + ax.tick_params(labelsize=7) + + +def plot_map_grid(data_list, titles, lons2d, lats2d, suptitle, + units=None, cmap='viridis', ncols=3, cell_size=(4.2, 3.2)): + """Create a grid of 2-D spatial maps. + + units : str or list of str (one per map) + cmap : str or list of str (one per map) + """ + n = len(data_list) + if n == 0: + fig, ax = plt.subplots(figsize=(6, 4)) + ax.text(0.5, 0.5, 'No data', ha='center', va='center', transform=ax.transAxes) + ax.axis('off') + return fig + units_list = ([units] * n) if (units is None or isinstance(units, str)) else list(units) + if units is None: + units_list = [''] * n + cmap_list = [cmap] * n if isinstance(cmap, str) else list(cmap) + nrows = max(1, (n + ncols - 1) // ncols) + fig, axes = plt.subplots(nrows, ncols, + figsize=(ncols * cell_size[0], nrows * cell_size[1]), + squeeze=False) + axes_flat = axes.flatten() + for i in range(n): + plot_map(axes_flat[i], data_list[i], lons2d, lats2d, + titles[i], units_list[i], cmap=cmap_list[i]) + for i in range(n, len(axes_flat)): + axes_flat[i].axis('off') + fig.suptitle(suptitle, fontsize=14, fontweight='bold') + plt.tight_layout() + return fig + + +def create_domain_map_figure(lons2d, lats2d, figsize=(10, 7)): + """Create a domain overview map with the bounding box highlighted. + Returns None if cartopy is not available.""" + if not HAS_CARTOPY: + return None + lon_min = float(np.nanmin(lons2d)) + lon_max = float(np.nanmax(lons2d)) + lat_min = float(np.nanmin(lats2d)) + lat_max = float(np.nanmax(lats2d)) + margin_lon = max(2.0, (lon_max - lon_min) * 0.15) + margin_lat = max(2.0, (lat_max - lat_min) * 0.15) + title = (f'Domain Extent ' + f'{lon_min:.2f}\u00b0\u2013{lon_max:.2f}\u00b0E, ' + f'{lat_min:.2f}\u00b0\u2013{lat_max:.2f}\u00b0N ' + f'({lons2d.shape[1]}\u00d7{lons2d.shape[0]} cells)') + fig = plt.figure(figsize=figsize) + ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) + extent = [lon_min - margin_lon, lon_max + margin_lon, + lat_min - margin_lat, lat_max + margin_lat] + ax.set_extent(extent, crs=ccrs.PlateCarree()) + ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='#c8e6f5') + ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='#f2efe8') + ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8, edgecolor='#555555') + ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.7, edgecolor='#888888') + ax.add_feature(cfeature.LAKES.with_scale('50m'), facecolor='#c8e6f5', alpha=0.8) + ax.add_feature(cfeature.RIVERS.with_scale('50m'), linewidth=0.4, edgecolor='#6baed6', alpha=0.6) + ax.plot([lon_min, lon_max, lon_max, lon_min, lon_min], + [lat_min, lat_min, lat_max, lat_max, lat_min], + color='#e53e3e', linewidth=2.5, transform=ccrs.PlateCarree(), zorder=5) + gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, + linewidth=0.5, color='gray', alpha=0.5, linestyle='--') + gl.top_labels = False + gl.right_labels = False + ax.set_title(title, fontsize=12, fontweight='bold') + return fig + + +def plot_monthly_timeseries_regional(ax, data, pft_indices, title, units, + pft_names, land_mask): + """Plot monthly time series for selected PFTs as spatial means over land cells.""" + colors = plt.cm.tab10(np.linspace(0, 1, max(len(pft_indices), 1))) + for pft_idx, color in zip(pft_indices, colors): + values = np.array([ + float(np.nanmean(np.where(land_mask, data[m, pft_idx, :, :], np.nan))) + for m in range(12) + ]) + if np.any(values > 0): + label = pft_names.get(pft_idx, f'PFT {pft_idx}') + if len(label) > 25: + label = label[:22] + '...' + ax.plot(range(12), values, marker='o', label=label, color=color, linewidth=2) + ax.set_xticks(range(12)) + ax.set_xticklabels(MONTH_NAMES) + ax.set_xlabel('Month') + ax.set_ylabel(f'[{units}] (spatial mean)') + ax.set_title(title, fontsize=12, fontweight='bold') + ax.legend(loc='upper left', bbox_to_anchor=(1, 1), fontsize=7) + ax.grid(alpha=0.3) + + def create_scalar_card(ax, name, value, units, long_name): """Create a card-style visualization for a scalar variable.""" ax.set_xlim(0, 10) @@ -445,7 +599,7 @@ def create_site_location_figure(lons, lats, labels=None, zoom_margin=10, figsize return fig -def create_html_report(figures_data, nc_file, output_dir): +def create_html_report(figures_data, nc_file, output_dir, metadata_str=''): """Generate an HTML report with all figures embedded.""" html_template = """ @@ -617,8 +771,7 @@ def create_html_report(figures_data, nc_file, output_dir):
    Source File: {nc_file}
    Generated: {timestamp}
    - Site: BE-Lon (Belgium - Lonzee)
    - Coordinates: {lon:.3f}E, {lat:.3f}N + {metadata_str}
    @@ -671,15 +824,10 @@ def create_html_report(figures_data, nc_file, output_dir): """ - # Get coordinates from the first figure data - lon = figures_data[0].get('lon', 0) - lat = figures_data[0].get('lat', 0) - html_content = html_template.format( nc_file=os.path.basename(nc_file), timestamp=datetime.now().strftime('%Y-%m-%d %H:%M:%S'), - lon=lon, - lat=lat, + metadata_str=metadata_str, nav_links=nav_links_html, sections=sections_html ) @@ -722,9 +870,34 @@ def main(nc_file): print(f"Reading NetCDF file: {nc_file}") nc = Dataset(nc_file, 'r') - # Get coordinates - lon = get_scalar_value(nc.variables['LONGXY']) - lat = get_scalar_value(nc.variables['LATIXY']) + # Detect grid type + is_regional = detect_grid_type(nc) + print(f"Grid mode: {'regional' if is_regional else 'single-site'}") + + # Get coordinates / domain info + if is_regional: + lons = get_field_2d(nc.variables['LONGXY']) + lats = get_field_2d(nc.variables['LATIXY']) + lon_min = float(np.nanmin(lons)) + lon_max = float(np.nanmax(lons)) + lat_min = float(np.nanmin(lats)) + lat_max = float(np.nanmax(lats)) + # Land mask for spatial averages + if 'LANDFRAC_PFT' in nc.variables: + land_mask = get_field_2d(nc.variables['LANDFRAC_PFT']) > 0 + elif 'PFTDATA_MASK' in nc.variables: + land_mask = get_field_2d(nc.variables['PFTDATA_MASK']).astype(bool) + else: + land_mask = np.ones(lons.shape, dtype=bool) + metadata_str = (f'Mode: Regional grid ' + f'({lons.shape[1]}\u00d7{lons.shape[0]} cells)
    ' + f'Domain: ' + f'{lon_min:.2f}\u00b0\u2013{lon_max:.2f}\u00b0E, ' + f'{lat_min:.2f}\u00b0\u2013{lat_max:.2f}\u00b0N') + else: + lon = get_scalar_value(nc.variables['LONGXY']) + lat = get_scalar_value(nc.variables['LATIXY']) + metadata_str = f'Coordinates: {lon:.3f}\u00b0E, {lat:.3f}\u00b0N' figures_data = [] @@ -734,60 +907,108 @@ def main(nc_file): print("Creating Section 1: Location and Basic Parameters...") section1_figures = [] - # Figure 1.0: Site location map - fig_map = create_site_location_figure([lon], [lat]) - if fig_map is not None: - pdf_path = os.path.join(pdf_dir, '00_site_location.pdf') - fig_map.savefig(pdf_path, bbox_inches='tight') + if is_regional: + # Domain overview map + fig_map = create_domain_map_figure(lons, lats) + if fig_map is not None: + pdf_path = os.path.join(pdf_dir, '00_domain_overview.pdf') + fig_map.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': (f'Domain extent: {lon_min:.2f}\u00b0\u2013{lon_max:.2f}\u00b0E, ' + f'{lat_min:.2f}\u00b0\u2013{lat_max:.2f}\u00b0N ' + f'({lons.shape[1]}\u00d7{lons.shape[0]} cells).'), + 'base64': fig_to_base64(fig_map) + }) + plt.close(fig_map) + + # Grid of basic parameter maps + basic_map_vars = [ + ('AREA', 'km\u00b2', 'viridis'), + ('FMAX', 'unitless', 'Blues'), + ('SOIL_COLOR','index', 'tab20b'), + ('zbedrock', 'm', 'copper'), + ('SLOPE', 'degrees', 'YlOrRd'), + ('STD_ELEV', 'm', 'terrain'), + ('LAKEDEPTH', 'm', 'Blues'), + ('peatf', 'fraction', 'Greens'), + ('gdp', 'unitless', 'YlOrBr'), + ] + data_list, titles_list, units_list, cmaps_list = [], [], [], [] + for var_name, units_str, cmap_str in basic_map_vars: + if var_name in nc.variables: + d = get_field_2d(nc.variables[var_name]) + if d.ndim == 2: + data_list.append(d) + titles_list.append(var_name) + units_list.append(units_str) + cmaps_list.append(cmap_str) + fig = plot_map_grid(data_list, titles_list, lons, lats, + 'Basic Parameters \u2013 Spatial Distribution', + units=units_list, cmap=cmaps_list, ncols=3) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '01_basic_parameters.pdf') + fig.savefig(pdf_path, bbox_inches='tight') section1_figures.append({ 'pdf_name': os.path.basename(pdf_path), - 'caption': f'Political map showing the site location at {lon:.3f}°E, {lat:.3f}°N.', - 'base64': fig_to_base64(fig_map) + 'caption': 'Spatial maps of basic parameters: area, FMAX, soil colour, bedrock depth, slope, elevation std, lake depth, peatland fraction, GDP.', + 'base64': fig_to_base64(fig) }) - plt.close(fig_map) - - # Figure 1.1: Location and basic info - fig, axes = plt.subplots(3, 4, figsize=(16, 12)) - fig.suptitle('Basic Site Parameters', fontsize=16, fontweight='bold') - - scalar_vars = [ - ('LONGXY', 'degrees east'), - ('LATIXY', 'degrees north'), - ('AREA', 'km^2'), - ('FMAX', 'unitless'), - ('SOIL_COLOR', 'index'), - ('mxsoil_color', 'unitless'), - ('zbedrock', 'm'), - ('SLOPE', 'degrees'), - ('STD_ELEV', 'm'), - ('LAKEDEPTH', 'm'), - ('gdp', 'unitless'), - ('abm', 'month') - ] - - for ax, (var_name, units) in zip(axes.flatten(), scalar_vars): - var = nc.variables[var_name] - value = get_scalar_value(var) - long_name = var.long_name if hasattr(var, 'long_name') else var_name - create_scalar_card(ax, var_name, value, units, long_name) + plt.close(fig) - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '01_basic_parameters.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section1_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Basic site parameters including location, area, soil color, bedrock depth, slope, and economic indicators.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + else: + # Figure 1.0: Site location map + fig_map = create_site_location_figure([lon], [lat]) + if fig_map is not None: + pdf_path = os.path.join(pdf_dir, '00_site_location.pdf') + fig_map.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': f'Political map showing the site location at {lon:.3f}\u00b0E, {lat:.3f}\u00b0N.', + 'base64': fig_to_base64(fig_map) + }) + plt.close(fig_map) + + # Figure 1.1: Location and basic info + fig, axes = plt.subplots(3, 4, figsize=(16, 12)) + fig.suptitle('Basic Site Parameters', fontsize=16, fontweight='bold') + + scalar_vars = [ + ('LONGXY', 'degrees east'), + ('LATIXY', 'degrees north'), + ('AREA', 'km^2'), + ('FMAX', 'unitless'), + ('SOIL_COLOR', 'index'), + ('mxsoil_color', 'unitless'), + ('zbedrock', 'm'), + ('SLOPE', 'degrees'), + ('STD_ELEV', 'm'), + ('LAKEDEPTH', 'm'), + ('gdp', 'unitless'), + ('abm', 'month') + ] + + for ax, (var_name, units) in zip(axes.flatten(), scalar_vars): + var = nc.variables[var_name] + value = get_scalar_value(var) + long_name = var.long_name if hasattr(var, 'long_name') else var_name + create_scalar_card(ax, var_name, value, units, long_name) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '01_basic_parameters.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section1_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Basic site parameters including location, area, soil color, bedrock depth, slope, and economic indicators.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'basic-params', - 'title': 'Basic Site Parameters', + 'title': 'Basic Site Parameters' if not is_regional else 'Domain Overview & Basic Parameters', 'description': 'Fundamental parameters describing the site location, topography, and basic soil/lake properties.', 'figures': section1_figures, - 'lon': lon, - 'lat': lat }) # ========================================================================= @@ -796,47 +1017,94 @@ def main(nc_file): print("Creating Section 2: Land Cover Fractions...") section2_figures = [] - # Figure 2.1: Major land cover pie chart - fig, axes = plt.subplots(1, 2, figsize=(14, 6)) - - # Major land units - pct_natveg = get_scalar_value(nc.variables['PCT_NATVEG']) - pct_crop = get_scalar_value(nc.variables['PCT_CROP']) - pct_urban = np.sum(get_1d_array(nc.variables['PCT_URBAN'])) - pct_lake = get_scalar_value(nc.variables['PCT_LAKE']) - pct_wetland = get_scalar_value(nc.variables['PCT_WETLAND']) - pct_glacier = get_scalar_value(nc.variables['PCT_GLACIER']) - - land_labels = ['Natural Vegetation', 'Cropland', 'Urban', 'Lake', 'Wetland', 'Glacier'] - land_values = np.array([pct_natveg, pct_crop, pct_urban, pct_lake, pct_wetland, pct_glacier]) - - plot_pie_chart(axes[0], land_labels, land_values, 'Major Land Cover Types', min_pct=0.1) - - # Land fractions bar chart - axes[1].bar(land_labels, land_values, color=['#48bb78', '#ecc94b', '#a0aec0', '#4299e1', '#9f7aea', '#63b3ed'], - edgecolor='#2d3748') - axes[1].set_ylabel('Percent (%)') - axes[1].set_title('Land Cover Distribution', fontsize=12, fontweight='bold') - axes[1].tick_params(axis='x', rotation=45) - for i, v in enumerate(land_values): - if v > 0: - axes[1].text(i, v + 1, f'{v:.1f}%', ha='center', fontsize=9) - axes[1].grid(axis='y', alpha=0.3) + if is_regional: + # Spatial maps of major land cover fractions + lc_vars = [ + ('PCT_NATVEG', 'Natural Vegetation (%)', 'Greens'), + ('PCT_CROP', 'Cropland (%)', 'YlOrBr'), + ('PCT_LAKE', 'Lake (%)', 'Blues'), + ('PCT_WETLAND', 'Wetland (%)', 'PuBuGn'), + ('PCT_GLACIER', 'Glacier (%)', 'cool'), + ] + data_list, titles_list, units_list, cmaps_list = [], [], [], [] + for var_name, title_str, cmap_str in lc_vars: + if var_name in nc.variables: + data_list.append(get_field_2d(nc.variables[var_name])) + titles_list.append(title_str) + units_list.append('%') + cmaps_list.append(cmap_str) + # Total urban (sum over density classes) + if 'PCT_URBAN' in nc.variables: + pct_urban_3d = np.array(nc.variables['PCT_URBAN'][:]) + if hasattr(pct_urban_3d, 'mask'): + pct_urban_3d = np.ma.filled(pct_urban_3d, 0.0) + data_list.append(pct_urban_3d.sum(axis=0)) + titles_list.append('Urban total (%)') + units_list.append('%') + cmaps_list.append('Reds') + fig = plot_map_grid(data_list, titles_list, lons, lats, + 'Land Cover Fractions \u2013 Spatial Distribution', + units=units_list, cmap=cmaps_list, ncols=3) + pdf_path = os.path.join(pdf_dir, '02_land_cover_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section2_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Spatial maps of land cover fractions: natural vegetation, cropland, lake, wetland, glacier, and total urban.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '02_land_cover_major.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section2_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Distribution of major land cover types: natural vegetation, cropland, urban, lake, wetland, and glacier.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + # Also store domain-mean land values for use in later sections + pct_natveg = float(np.nanmean(np.where(land_mask, get_field_2d(nc.variables['PCT_NATVEG']), np.nan))) + pct_crop = float(np.nanmean(np.where(land_mask, get_field_2d(nc.variables['PCT_CROP']), np.nan))) + pct_urban_arr = data_list[titles_list.index('Urban total (%)')] + pct_urban = np.array([ + float(np.nanmean(np.where(land_mask, get_field_2d(nc.variables['PCT_URBAN'])[i], np.nan))) + for i in range(nc.variables['PCT_URBAN'].shape[0]) + ]) + + else: + # Figure 2.1: Major land cover pie chart + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + + pct_natveg = get_scalar_value(nc.variables['PCT_NATVEG']) + pct_crop = get_scalar_value(nc.variables['PCT_CROP']) + pct_urban = get_1d_array(nc.variables['PCT_URBAN']) + pct_lake = get_scalar_value(nc.variables['PCT_LAKE']) + pct_wetland = get_scalar_value(nc.variables['PCT_WETLAND']) + pct_glacier = get_scalar_value(nc.variables['PCT_GLACIER']) + + land_labels = ['Natural Vegetation', 'Cropland', 'Urban', 'Lake', 'Wetland', 'Glacier'] + land_values = np.array([pct_natveg, pct_crop, np.sum(pct_urban), + pct_lake, pct_wetland, pct_glacier]) + + plot_pie_chart(axes[0], land_labels, land_values, 'Major Land Cover Types', min_pct=0.1) + + axes[1].bar(land_labels, land_values, + color=['#48bb78', '#ecc94b', '#a0aec0', '#4299e1', '#9f7aea', '#63b3ed'], + edgecolor='#2d3748') + axes[1].set_ylabel('Percent (%)') + axes[1].set_title('Land Cover Distribution', fontsize=12, fontweight='bold') + axes[1].tick_params(axis='x', rotation=45) + for i, v in enumerate(land_values): + if v > 0: + axes[1].text(i, v + 1, f'{v:.1f}%', ha='center', fontsize=9) + axes[1].grid(axis='y', alpha=0.3) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '02_land_cover_major.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section2_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Distribution of major land cover types: natural vegetation, cropland, urban, lake, wetland, and glacier.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'land-cover', 'title': 'Land Cover Fractions', - 'description': 'Overview of land surface cover type distribution. Key variables for understanding the dominant land use at the site.', + 'description': 'Overview of land surface cover type distribution.', 'figures': section2_figures }) @@ -846,45 +1114,94 @@ def main(nc_file): print("Creating Section 3: Natural PFT Distribution...") section3_figures = [] - fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + if is_regional: + pct_nat_pft_data = np.array(nc.variables['PCT_NAT_PFT'][:], dtype=float) + if hasattr(nc.variables['PCT_NAT_PFT'][:], 'mask'): + pct_nat_pft_data = np.ma.filled(nc.variables['PCT_NAT_PFT'][:], np.nan) + n_natpft = pct_nat_pft_data.shape[0] + + # Maps for PFTs present anywhere in the domain (max > 1%) + active_idx = [i for i in range(n_natpft) + if np.nanmax(pct_nat_pft_data[i]) > 1.0] + data_list = [pct_nat_pft_data[i] for i in active_idx] + titles_list = [NATPFT_NAMES.get(i, f'PFT {i}')[:22] for i in active_idx] + + if data_list: + ncols = min(4, len(data_list)) + fig = plot_map_grid(data_list, titles_list, lons, lats, + 'Natural PFT Distribution (% of natural veg landunit)', + units='%', cmap='YlGn', ncols=ncols) + pdf_path = os.path.join(pdf_dir, '03_natural_pft_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section3_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': f'Spatial maps of natural PFT fractions (showing {len(active_idx)} PFTs with max > 1%).', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Domain-mean bar chart + natpft_labels = [NATPFT_NAMES.get(i, f'PFT {i}') for i in range(n_natpft)] + pct_nat_pft_mean = np.array([ + float(np.nanmean(np.where(land_mask, pct_nat_pft_data[i], np.nan))) + for i in range(n_natpft) + ]) + fig, ax = plt.subplots(figsize=(12, 6)) + colors_bar = plt.cm.Greens(np.linspace(0.3, 0.9, n_natpft)) + ax.barh(range(n_natpft), pct_nat_pft_mean, color=colors_bar, edgecolor='#2d3748') + ax.set_yticks(range(n_natpft)) + ax.set_yticklabels(natpft_labels, fontsize=9) + ax.set_xlabel('Domain-mean percent of Natural Vegetation Landunit (%)') + ax.set_title('Natural PFTs \u2013 Domain Mean', fontsize=12, fontweight='bold') + ax.grid(axis='x', alpha=0.3) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '03b_natural_pft_mean.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section3_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean natural PFT fractions (land-cell average).', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - # Get PCT_NAT_PFT data - pct_nat_pft = get_1d_array(nc.variables['PCT_NAT_PFT']) - natpft_labels = [NATPFT_NAMES.get(i, f'PFT {i}') for i in range(15)] + # Keep domain-mean for summary use + pct_nat_pft = pct_nat_pft_mean - # Pie chart for non-zero PFTs - plot_pie_chart(axes[0], natpft_labels, pct_nat_pft, - f'Natural PFT Distribution\n(within {pct_natveg:.1f}% natural veg)', min_pct=0.5) + else: + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) - # Horizontal bar chart - colors = plt.cm.Greens(np.linspace(0.3, 0.9, 15)) - y_pos = range(15) - axes[1].barh(y_pos, pct_nat_pft, color=colors, edgecolor='#2d3748') - axes[1].set_yticks(y_pos) - axes[1].set_yticklabels(natpft_labels, fontsize=9) - axes[1].set_xlabel('Percent of Natural Vegetation Landunit (%)') - axes[1].set_title('Natural Plant Functional Types', fontsize=12, fontweight='bold') - axes[1].grid(axis='x', alpha=0.3) + pct_nat_pft = get_1d_array(nc.variables['PCT_NAT_PFT']) + natpft_labels = [NATPFT_NAMES.get(i, f'PFT {i}') for i in range(15)] - # Add value labels for non-zero - for i, v in enumerate(pct_nat_pft): - if v > 0: - axes[1].text(v + 1, i, f'{v:.1f}%', va='center', fontsize=8) + plot_pie_chart(axes[0], natpft_labels, pct_nat_pft, + f'Natural PFT Distribution\n(within {pct_natveg:.1f}% natural veg)', min_pct=0.5) - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '03_natural_pft.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section3_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Distribution of natural Plant Functional Types (PFTs). Values represent percentages within the natural vegetation landunit.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + colors = plt.cm.Greens(np.linspace(0.3, 0.9, 15)) + y_pos = range(15) + axes[1].barh(y_pos, pct_nat_pft, color=colors, edgecolor='#2d3748') + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(natpft_labels, fontsize=9) + axes[1].set_xlabel('Percent of Natural Vegetation Landunit (%)') + axes[1].set_title('Natural Plant Functional Types', fontsize=12, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + for i, v in enumerate(pct_nat_pft): + if v > 0: + axes[1].text(v + 1, i, f'{v:.1f}%', va='center', fontsize=8) + + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '03_natural_pft.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section3_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Distribution of natural Plant Functional Types (PFTs). Values represent percentages within the natural vegetation landunit.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'natural-pft', 'title': 'Natural PFT Distribution', - 'description': 'Natural Plant Functional Type distribution within the natural vegetation landunit. Important for understanding vegetation dynamics.', + 'description': 'Natural Plant Functional Type distribution within the natural vegetation landunit.', 'figures': section3_figures }) @@ -894,82 +1211,126 @@ def main(nc_file): print("Creating Section 4: Crop Functional Types...") section4_figures = [] - # Get PCT_CFT data - pct_cft = get_1d_array(nc.variables['PCT_CFT']) cft_indices = nc.variables['cft'][:] - # Only show non-zero CFTs - nonzero_mask = pct_cft > 0.1 - nonzero_cft_values = pct_cft[nonzero_mask] - nonzero_cft_indices = cft_indices[nonzero_mask] - nonzero_cft_labels = [CFT_NAMES.get(int(i), f'CFT {i}') for i in nonzero_cft_indices] + if is_regional: + pct_cft_data = np.array(nc.variables['PCT_CFT'][:], dtype=float) + if hasattr(nc.variables['PCT_CFT'][:], 'mask'): + pct_cft_data = np.ma.filled(nc.variables['PCT_CFT'][:], np.nan) + n_cft = pct_cft_data.shape[0] + + # Maps for CFTs present anywhere (max > 1%) + active_cft_idx = [i for i in range(n_cft) + if np.nanmax(pct_cft_data[i]) > 1.0] + if active_cft_idx: + data_list = [pct_cft_data[i] for i in active_cft_idx] + titles_list = [CFT_NAMES.get(int(cft_indices[i]), f'CFT {cft_indices[i]}')[:22] + for i in active_cft_idx] + ncols = min(4, len(data_list)) + fig = plot_map_grid(data_list, titles_list, lons, lats, + 'Crop Functional Types (% of crop landunit)', + units='%', cmap='YlOrBr', ncols=ncols) + pdf_path = os.path.join(pdf_dir, '04_cft_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': f'Spatial maps of crop functional types (showing {len(active_cft_idx)} CFTs with max > 1%).', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Domain-mean fertilizer bar chart + fertnitro_data = np.array(nc.variables['CONST_FERTNITRO_CFT'][:], dtype=float) + if hasattr(nc.variables['CONST_FERTNITRO_CFT'][:], 'mask'): + fertnitro_data = np.ma.filled(nc.variables['CONST_FERTNITRO_CFT'][:], np.nan) + fert_mean = np.array([ + float(np.nanmean(np.where(land_mask, fertnitro_data[i], np.nan))) + for i in range(n_cft) + ]) + nonzero_fert_mask = fert_mean > 0 + if np.any(nonzero_fert_mask): + nz_fert = fert_mean[nonzero_fert_mask] + nz_fert_labels = [CFT_NAMES.get(int(cft_indices[i]), f'CFT {cft_indices[i]}') + for i in range(n_cft) if nonzero_fert_mask[i]] + fig, ax = plt.subplots(figsize=(14, 6)) + colors = plt.cm.Blues(np.linspace(0.4, 0.9, len(nz_fert))) + ax.bar(range(len(nz_fert)), nz_fert, color=colors, edgecolor='#2d3748') + ax.set_xticks(range(len(nz_fert))) + ax.set_xticklabels(nz_fert_labels, rotation=45, ha='right', fontsize=9) + ax.set_ylabel('N Fertilizer \u2013 Domain Mean (gN/m\u00b2/yr)') + ax.set_title('Nitrogen Fertilizer \u2013 Domain Mean by Crop Type', fontsize=12, fontweight='bold') + ax.grid(axis='y', alpha=0.3) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04b_nitrogen_fertilizer.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean nitrogen fertilizer application rates for non-zero crop types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - fig, axes = plt.subplots(1, 2, figsize=(16, 7)) - - if len(nonzero_cft_values) > 0: - # Pie chart - plot_pie_chart(axes[0], nonzero_cft_labels, nonzero_cft_values, - f'Crop Type Distribution\n(within {pct_crop:.1f}% cropland)', min_pct=0.1) - - # Bar chart - colors = plt.cm.YlOrBr(np.linspace(0.3, 0.9, len(nonzero_cft_values))) - y_pos = range(len(nonzero_cft_values)) - axes[1].barh(y_pos, nonzero_cft_values, color=colors, edgecolor='#2d3748') - axes[1].set_yticks(y_pos) - axes[1].set_yticklabels(nonzero_cft_labels, fontsize=9) - axes[1].set_xlabel('Percent of Crop Landunit (%)') - axes[1].set_title('Crop Functional Types (non-zero only)', fontsize=12, fontweight='bold') - axes[1].grid(axis='x', alpha=0.3) - - for i, v in enumerate(nonzero_cft_values): - axes[1].text(v + 0.5, i, f'{v:.1f}%', va='center', fontsize=9) - else: - axes[0].text(0.5, 0.5, 'No crops present', ha='center', va='center', transform=axes[0].transAxes) - axes[1].text(0.5, 0.5, 'No crops present', ha='center', va='center', transform=axes[1].transAxes) - - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '04_crop_functional_types.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section4_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Distribution of Crop Functional Types (CFTs). Values represent percentages within the crop landunit.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) - - # Nitrogen fertilizer for crops - fig, ax = plt.subplots(figsize=(14, 6)) - fertnitro = get_1d_array(nc.variables['CONST_FERTNITRO_CFT']) - nonzero_fert_mask = fertnitro > 0 - - if np.any(nonzero_fert_mask): - nonzero_fert = fertnitro[nonzero_fert_mask] - nonzero_fert_indices = cft_indices[nonzero_fert_mask] - nonzero_fert_labels = [CFT_NAMES.get(int(i), f'CFT {i}') for i in nonzero_fert_indices] - - colors = plt.cm.Blues(np.linspace(0.4, 0.9, len(nonzero_fert))) - x_pos = range(len(nonzero_fert)) - ax.bar(x_pos, nonzero_fert, color=colors, edgecolor='#2d3748') - ax.set_xticks(x_pos) - ax.set_xticklabels(nonzero_fert_labels, rotation=45, ha='right', fontsize=9) - ax.set_ylabel('Nitrogen Fertilizer (gN/m2/yr)') - ax.set_title('Nitrogen Fertilizer Application by Crop Type', fontsize=12, fontweight='bold') - ax.grid(axis='y', alpha=0.3) - - for i, v in enumerate(nonzero_fert): - ax.text(i, v + 0.3, f'{v:.1f}', ha='center', fontsize=9) else: - ax.text(0.5, 0.5, 'No fertilizer data', ha='center', va='center', transform=ax.transAxes) - - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '04b_nitrogen_fertilizer.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section4_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Nitrogen fertilizer application rates for different crop types.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + pct_cft = get_1d_array(nc.variables['PCT_CFT']) + nonzero_mask = pct_cft > 0.1 + nonzero_cft_values = pct_cft[nonzero_mask] + nonzero_cft_indices = cft_indices[nonzero_mask] + nonzero_cft_labels = [CFT_NAMES.get(int(i), f'CFT {i}') for i in nonzero_cft_indices] + + fig, axes = plt.subplots(1, 2, figsize=(16, 7)) + if len(nonzero_cft_values) > 0: + plot_pie_chart(axes[0], nonzero_cft_labels, nonzero_cft_values, + f'Crop Type Distribution\n(within {pct_crop:.1f}% cropland)', min_pct=0.1) + colors = plt.cm.YlOrBr(np.linspace(0.3, 0.9, len(nonzero_cft_values))) + y_pos = range(len(nonzero_cft_values)) + axes[1].barh(y_pos, nonzero_cft_values, color=colors, edgecolor='#2d3748') + axes[1].set_yticks(y_pos) + axes[1].set_yticklabels(nonzero_cft_labels, fontsize=9) + axes[1].set_xlabel('Percent of Crop Landunit (%)') + axes[1].set_title('Crop Functional Types (non-zero only)', fontsize=12, fontweight='bold') + axes[1].grid(axis='x', alpha=0.3) + for i, v in enumerate(nonzero_cft_values): + axes[1].text(v + 0.5, i, f'{v:.1f}%', va='center', fontsize=9) + else: + axes[0].text(0.5, 0.5, 'No crops present', ha='center', va='center', transform=axes[0].transAxes) + axes[1].text(0.5, 0.5, 'No crops present', ha='center', va='center', transform=axes[1].transAxes) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04_crop_functional_types.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Distribution of Crop Functional Types (CFTs). Values represent percentages within the crop landunit.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(14, 6)) + fertnitro = get_1d_array(nc.variables['CONST_FERTNITRO_CFT']) + nonzero_fert_mask = fertnitro > 0 + if np.any(nonzero_fert_mask): + nonzero_fert = fertnitro[nonzero_fert_mask] + nonzero_fert_indices = cft_indices[nonzero_fert_mask] + nonzero_fert_labels = [CFT_NAMES.get(int(i), f'CFT {i}') for i in nonzero_fert_indices] + colors = plt.cm.Blues(np.linspace(0.4, 0.9, len(nonzero_fert))) + ax.bar(range(len(nonzero_fert)), nonzero_fert, color=colors, edgecolor='#2d3748') + ax.set_xticks(range(len(nonzero_fert))) + ax.set_xticklabels(nonzero_fert_labels, rotation=45, ha='right', fontsize=9) + ax.set_ylabel('Nitrogen Fertilizer (gN/m\u00b2/yr)') + ax.set_title('Nitrogen Fertilizer Application by Crop Type', fontsize=12, fontweight='bold') + ax.grid(axis='y', alpha=0.3) + for i, v in enumerate(nonzero_fert): + ax.text(i, v + 0.3, f'{v:.1f}', ha='center', fontsize=9) + else: + ax.text(0.5, 0.5, 'No fertilizer data', ha='center', va='center', transform=ax.transAxes) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '04b_nitrogen_fertilizer.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section4_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Nitrogen fertilizer application rates for different crop types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'crop-types', @@ -984,61 +1345,126 @@ def main(nc_file): print("Creating Section 5: Soil Properties...") section5_figures = [] - fig, axes = plt.subplots(1, 3, figsize=(15, 6)) - - # PCT_SAND - pct_sand = get_1d_array(nc.variables['PCT_SAND']) - plot_soil_profile(axes[0], SOIL_DEPTHS, pct_sand, 'Sand Content', '%', '#f6ad55') - - # PCT_CLAY - pct_clay = get_1d_array(nc.variables['PCT_CLAY']) - plot_soil_profile(axes[1], SOIL_DEPTHS, pct_clay, 'Clay Content', '%', '#fc8181') - - # ORGANIC - organic = get_1d_array(nc.variables['ORGANIC']) - plot_soil_profile(axes[2], SOIL_DEPTHS, organic, 'Organic Matter', 'kg/m3', '#68d391') - - fig.suptitle('Soil Properties by Depth', fontsize=14, fontweight='bold') - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '05_soil_properties.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section5_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Vertical profiles of soil texture (sand and clay content) and organic matter density across 10 soil layers.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) - - # Combined soil texture visualization - fig, ax = plt.subplots(figsize=(10, 8)) - pct_silt = 100 - pct_sand - pct_clay # Calculate silt - - # Stacked bar chart for soil texture - x = range(len(SOIL_DEPTHS)) - width = 0.6 - - ax.barh(x, pct_sand, width, label='Sand', color='#f6ad55', edgecolor='#2d3748') - ax.barh(x, pct_silt, width, left=pct_sand, label='Silt', color='#a0aec0', edgecolor='#2d3748') - ax.barh(x, pct_clay, width, left=pct_sand+pct_silt, label='Clay', color='#fc8181', edgecolor='#2d3748') + if is_regional: + # Load 3-D soil arrays: (nlevsoi, nlat, nlon) + sand_3d = np.array(nc.variables['PCT_SAND'][:], dtype=float) + clay_3d = np.array(nc.variables['PCT_CLAY'][:], dtype=float) + org_3d = np.array(nc.variables['ORGANIC'][:], dtype=float) + for arr in [sand_3d, clay_3d, org_3d]: + if hasattr(arr, 'mask'): + arr = np.ma.filled(arr, np.nan) + + # Maps of depth-mean values + sand_mean2d = np.nanmean(sand_3d, axis=0) + clay_mean2d = np.nanmean(clay_3d, axis=0) + org_mean2d = np.nanmean(org_3d, axis=0) + + fig = plot_map_grid( + [sand_mean2d, clay_mean2d, org_mean2d], + ['Sand \u2013 depth mean', 'Clay \u2013 depth mean', 'Organic matter \u2013 depth mean'], + lons, lats, + 'Soil Properties \u2013 Depth-Averaged Spatial Distribution', + units=['%', '%', 'kg/m\u00b3'], + cmap=['YlOrBr', 'Reds', 'Greens'], + ncols=3 + ) + pdf_path = os.path.join(pdf_dir, '05_soil_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Depth-averaged spatial maps of sand content, clay content, and organic matter density.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Domain-mean soil profiles (mean ± std over land cells) + def _masked_profile(arr3d): + out_mean, out_std = [], [] + for lev in range(arr3d.shape[0]): + d = arr3d[lev] + vals = d[land_mask & np.isfinite(d)] + out_mean.append(float(np.nanmean(vals)) if vals.size else np.nan) + out_std.append(float(np.nanstd(vals)) if vals.size else np.nan) + return np.array(out_mean), np.array(out_std) + + sand_m, sand_s = _masked_profile(sand_3d) + clay_m, clay_s = _masked_profile(clay_3d) + org_m, org_s = _masked_profile(org_3d) + + fig, axes = plt.subplots(1, 3, figsize=(15, 6)) + for ax, m, s, title, units_str, col in zip( + axes, + [sand_m, clay_m, org_m], + [sand_s, clay_s, org_s], + ['Sand Content', 'Clay Content', 'Organic Matter'], + ['%', '%', 'kg/m\u00b3'], + ['#f6ad55', '#fc8181', '#68d391']): + ax.barh(range(len(SOIL_DEPTHS)), m, color=col, edgecolor='#2d3748', alpha=0.8, + xerr=s, error_kw=dict(ecolor='#2d3748', capsize=3)) + ax.set_yticks(range(len(SOIL_DEPTHS))) + ax.set_yticklabels([f'{d:.2f}m' for d in SOIL_DEPTHS]) + ax.invert_yaxis() + ax.set_xlabel(f'{title} [{units_str}]') + ax.set_title(f'{title}\n(domain mean \u00b1 std)', fontsize=11, fontweight='bold') + ax.grid(axis='x', alpha=0.3) + fig.suptitle('Soil Properties by Depth \u2013 Domain Mean', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05b_soil_profiles.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean soil profiles (error bars show spatial std) for sand, clay, and organic matter.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - ax.set_yticks(x) - ax.set_yticklabels([f'{d:.2f}m' for d in SOIL_DEPTHS]) - ax.invert_yaxis() - ax.set_xlabel('Percent (%)') - ax.set_ylabel('Depth') - ax.set_title('Soil Texture Composition by Depth', fontsize=12, fontweight='bold') - ax.legend(loc='upper right') - ax.set_xlim(0, 100) + # Store domain-mean 1-D profiles for summary use + pct_sand = sand_m + pct_clay = clay_m - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '05b_soil_texture.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section5_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Stacked bar chart showing soil texture composition (sand, silt, clay) at each soil layer.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + else: + fig, axes = plt.subplots(1, 3, figsize=(15, 6)) + pct_sand = get_1d_array(nc.variables['PCT_SAND']) + plot_soil_profile(axes[0], SOIL_DEPTHS, pct_sand, 'Sand Content', '%', '#f6ad55') + pct_clay = get_1d_array(nc.variables['PCT_CLAY']) + plot_soil_profile(axes[1], SOIL_DEPTHS, pct_clay, 'Clay Content', '%', '#fc8181') + organic = get_1d_array(nc.variables['ORGANIC']) + plot_soil_profile(axes[2], SOIL_DEPTHS, organic, 'Organic Matter', 'kg/m3', '#68d391') + fig.suptitle('Soil Properties by Depth', fontsize=14, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05_soil_properties.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Vertical profiles of soil texture (sand and clay content) and organic matter density across 10 soil layers.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(10, 8)) + pct_silt = 100 - pct_sand - pct_clay + x = range(len(SOIL_DEPTHS)) + width = 0.6 + ax.barh(x, pct_sand, width, label='Sand', color='#f6ad55', edgecolor='#2d3748') + ax.barh(x, pct_silt, width, left=pct_sand, label='Silt', color='#a0aec0', edgecolor='#2d3748') + ax.barh(x, pct_clay, width, left=pct_sand + pct_silt, label='Clay', color='#fc8181', edgecolor='#2d3748') + ax.set_yticks(x) + ax.set_yticklabels([f'{d:.2f}m' for d in SOIL_DEPTHS]) + ax.invert_yaxis() + ax.set_xlabel('Percent (%)') + ax.set_ylabel('Depth') + ax.set_title('Soil Texture Composition by Depth', fontsize=12, fontweight='bold') + ax.legend(loc='upper right') + ax.set_xlim(0, 100) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '05b_soil_texture.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section5_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Stacked bar chart showing soil texture composition (sand, silt, clay) at each soil layer.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'soil', @@ -1053,65 +1479,109 @@ def main(nc_file): print("Creating Section 6: Monthly Vegetation Parameters...") section6_figures = [] - # Find PFTs with non-zero LAI - lai_data = nc.variables['MONTHLY_LAI'][:] - sai_data = nc.variables['MONTHLY_SAI'][:] - - # Identify active PFTs (those with non-zero LAI at any month) - active_pfts = [] - for pft in range(79): - if np.any(lai_data[:, pft, 0, 0] > 0): - active_pfts.append(pft) - - # Combine PFT names + lai_data = np.array(nc.variables['MONTHLY_LAI'][:], dtype=float) + sai_data = np.array(nc.variables['MONTHLY_SAI'][:], dtype=float) all_pft_names = {**NATPFT_NAMES, **CFT_NAMES} - # Plot LAI - fig, ax = plt.subplots(figsize=(14, 7)) - plot_monthly_timeseries(ax, lai_data, active_pfts[:10], - 'Monthly Leaf Area Index (LAI)', 'm2/m2', all_pft_names) - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '06_monthly_lai.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section6_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Monthly Leaf Area Index (LAI) time series for active Plant Functional Types.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) - - # Plot SAI - fig, ax = plt.subplots(figsize=(14, 7)) - plot_monthly_timeseries(ax, sai_data, active_pfts[:10], - 'Monthly Stem Area Index (SAI)', 'm2/m2', all_pft_names) - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '06b_monthly_sai.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section6_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Monthly Stem Area Index (SAI) time series for active Plant Functional Types.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) - - # Plot canopy heights - height_top = nc.variables['MONTHLY_HEIGHT_TOP'][:] - height_bot = nc.variables['MONTHLY_HEIGHT_BOT'][:] + if is_regional: + # Active PFTs: non-zero anywhere in domain + active_pfts = [pft for pft in range(lai_data.shape[1]) + if np.nanmax(lai_data[:, pft, :, :]) > 0] + + # Domain-mean timeseries + fig, ax = plt.subplots(figsize=(14, 7)) + plot_monthly_timeseries_regional(ax, lai_data, active_pfts[:10], + 'Monthly LAI \u2013 Domain Mean', 'm\u00b2/m\u00b2', + all_pft_names, land_mask) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06_monthly_lai.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean monthly LAI for active PFTs (spatial average over land cells).', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(14, 7)) + plot_monthly_timeseries_regional(ax, sai_data, active_pfts[:10], + 'Monthly SAI \u2013 Domain Mean', 'm\u00b2/m\u00b2', + all_pft_names, land_mask) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06b_monthly_sai.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Domain-mean monthly SAI for active PFTs.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Annual-mean LAI maps for top active PFTs + n_map_pfts = min(8, len(active_pfts)) + if n_map_pfts > 0: + data_list = [np.nanmean(lai_data[:, p, :, :], axis=0) for p in active_pfts[:n_map_pfts]] + titles_list = [all_pft_names.get(p, f'PFT {p}')[:22] for p in active_pfts[:n_map_pfts]] + ncols = min(4, n_map_pfts) + fig = plot_map_grid(data_list, titles_list, lons, lats, + 'Annual-Mean LAI by PFT', + units='m\u00b2/m\u00b2', cmap='YlGn', ncols=ncols) + pdf_path = os.path.join(pdf_dir, '06c_lai_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': f'Spatial maps of annual-mean LAI for the {n_map_pfts} most active PFTs.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) - fig, axes = plt.subplots(1, 2, figsize=(14, 6)) - plot_monthly_timeseries(axes[0], height_top, active_pfts[:10], - 'Monthly Canopy Top Height', 'm', all_pft_names) - plot_monthly_timeseries(axes[1], height_bot, active_pfts[:10], - 'Monthly Canopy Bottom Height', 'm', all_pft_names) - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '06c_canopy_heights.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section6_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Monthly canopy top and bottom heights for active Plant Functional Types.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + else: + # Identify active PFTs (non-zero LAI at any month at the single cell) + active_pfts = [pft for pft in range(79) + if np.any(lai_data[:, pft, 0, 0] > 0)] + + fig, ax = plt.subplots(figsize=(14, 7)) + plot_monthly_timeseries(ax, lai_data, active_pfts[:10], + 'Monthly Leaf Area Index (LAI)', 'm2/m2', all_pft_names) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06_monthly_lai.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Monthly Leaf Area Index (LAI) time series for active Plant Functional Types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + fig, ax = plt.subplots(figsize=(14, 7)) + plot_monthly_timeseries(ax, sai_data, active_pfts[:10], + 'Monthly Stem Area Index (SAI)', 'm2/m2', all_pft_names) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06b_monthly_sai.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Monthly Stem Area Index (SAI) time series for active Plant Functional Types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + height_top = nc.variables['MONTHLY_HEIGHT_TOP'][:] + height_bot = nc.variables['MONTHLY_HEIGHT_BOT'][:] + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + plot_monthly_timeseries(axes[0], height_top, active_pfts[:10], + 'Monthly Canopy Top Height', 'm', all_pft_names) + plot_monthly_timeseries(axes[1], height_bot, active_pfts[:10], + 'Monthly Canopy Bottom Height', 'm', all_pft_names) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '06c_canopy_heights.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section6_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Monthly canopy top and bottom heights for active Plant Functional Types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'vegetation', @@ -1126,12 +1596,62 @@ def main(nc_file): print("Creating Section 7: Urban Parameters...") section7_figures = [] - # Urban percent by type - pct_urban = get_1d_array(nc.variables['PCT_URBAN']) - + if is_regional: + # Maps of PCT_URBAN per density class + pct_urban_3d = np.array(nc.variables['PCT_URBAN'][:], dtype=float) + data_list = [pct_urban_3d[i] for i in range(len(URBAN_TYPES))] + titles_list = [f'PCT_URBAN\n{t}' for t in URBAN_TYPES] + fig = plot_map_grid(data_list, titles_list, lons, lats, + 'Urban Fraction by Density Class', + units='%', cmap='Reds', ncols=3) + pdf_path = os.path.join(pdf_dir, '07_urban_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section7_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Spatial maps of urban fraction for Tall Building District, High Density, and Medium Density classes.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Domain-mean bar charts for building parameters + def _urban_mean(var_name): + d = np.array(nc.variables[var_name][:], dtype=float) + return np.array([ + float(np.nanmean(np.where(land_mask, d[i], np.nan))) + for i in range(len(URBAN_TYPES)) + ]) + + pct_urban = _urban_mean('PCT_URBAN') + canyon_hwr = _urban_mean('CANYON_HWR') + ht_roof = _urban_mean('HT_ROOF') + t_building = _urban_mean('T_BUILDING_MIN') + wtlunit_roof = _urban_mean('WTLUNIT_ROOF') + wtroad_perv = _urban_mean('WTROAD_PERV') + thick_roof = _urban_mean('THICK_ROOF') + thick_wall = _urban_mean('THICK_WALL') + em_improad = _urban_mean('EM_IMPROAD') + em_perroad = _urban_mean('EM_PERROAD') + em_roof = _urban_mean('EM_ROOF') + em_wall = _urban_mean('EM_WALL') + + # For single-site: read all urban parameter vectors directly + if not is_regional: + pct_urban = get_1d_array(nc.variables['PCT_URBAN']) + canyon_hwr = get_1d_array(nc.variables['CANYON_HWR']) + ht_roof = get_1d_array(nc.variables['HT_ROOF']) + t_building = get_1d_array(nc.variables['T_BUILDING_MIN']) + wtlunit_roof = get_1d_array(nc.variables['WTLUNIT_ROOF']) + wtroad_perv = get_1d_array(nc.variables['WTROAD_PERV']) + thick_roof = get_1d_array(nc.variables['THICK_ROOF']) + thick_wall = get_1d_array(nc.variables['THICK_WALL']) + em_improad = get_1d_array(nc.variables['EM_IMPROAD']) + em_perroad = get_1d_array(nc.variables['EM_PERROAD']) + em_roof = get_1d_array(nc.variables['EM_ROOF']) + em_wall = get_1d_array(nc.variables['EM_WALL']) + + # Shared bar chart for building parameters (works for both modes) fig, axes = plt.subplots(2, 3, figsize=(15, 10)) - # Urban fraction axes[0, 0].bar(URBAN_TYPES, pct_urban, color=['#e53e3e', '#dd6b20', '#d69e2e'], edgecolor='#2d3748') axes[0, 0].set_ylabel('Percent (%)') axes[0, 0].set_title('Urban Fraction by Density Type', fontsize=11, fontweight='bold') @@ -1139,45 +1659,31 @@ def main(nc_file): for i, v in enumerate(pct_urban): axes[0, 0].text(i, v + 0.1, f'{v:.2f}%', ha='center', fontsize=9) - # Canyon height-to-width ratio - canyon_hwr = get_1d_array(nc.variables['CANYON_HWR']) axes[0, 1].bar(URBAN_TYPES, canyon_hwr, color='#805ad5', edgecolor='#2d3748') axes[0, 1].set_ylabel('Height/Width Ratio') axes[0, 1].set_title('Canyon Height-to-Width Ratio', fontsize=11, fontweight='bold') axes[0, 1].tick_params(axis='x', rotation=15) - # Roof height - ht_roof = get_1d_array(nc.variables['HT_ROOF']) axes[0, 2].bar(URBAN_TYPES, ht_roof, color='#3182ce', edgecolor='#2d3748') axes[0, 2].set_ylabel('Height (m)') axes[0, 2].set_title('Roof Height', fontsize=11, fontweight='bold') axes[0, 2].tick_params(axis='x', rotation=15) - # Building temperature minimum - t_building = get_1d_array(nc.variables['T_BUILDING_MIN']) - axes[1, 0].bar(URBAN_TYPES, t_building - 273.15, color='#e53e3e', edgecolor='#2d3748') # Convert to Celsius - axes[1, 0].set_ylabel('Temperature (C)') + axes[1, 0].bar(URBAN_TYPES, t_building - 273.15, color='#e53e3e', edgecolor='#2d3748') + axes[1, 0].set_ylabel('Temperature (\u00b0C)') axes[1, 0].set_title('Min. Building Interior Temp.', fontsize=11, fontweight='bold') axes[1, 0].tick_params(axis='x', rotation=15) - # Roof and pervious road fractions - wtlunit_roof = get_1d_array(nc.variables['WTLUNIT_ROOF']) - wtroad_perv = get_1d_array(nc.variables['WTROAD_PERV']) - x = np.arange(len(URBAN_TYPES)) width = 0.35 axes[1, 1].bar(x - width/2, wtlunit_roof, width, label='Roof Fraction', color='#805ad5', edgecolor='#2d3748') - axes[1, 1].bar(x + width/2, wtroad_perv, width, label='Pervious Road Fraction', color='#38a169', edgecolor='#2d3748') + axes[1, 1].bar(x + width/2, wtroad_perv, width, label='Pervious Road Fraction', color='#38a169', edgecolor='#2d3748') axes[1, 1].set_xticks(x) axes[1, 1].set_xticklabels(URBAN_TYPES, rotation=15) axes[1, 1].set_ylabel('Fraction') axes[1, 1].set_title('Urban Surface Fractions', fontsize=11, fontweight='bold') axes[1, 1].legend(fontsize=8) - # Wall and roof thickness - thick_roof = get_1d_array(nc.variables['THICK_ROOF']) - thick_wall = get_1d_array(nc.variables['THICK_WALL']) - axes[1, 2].bar(x - width/2, thick_roof, width, label='Roof Thickness', color='#4299e1', edgecolor='#2d3748') axes[1, 2].bar(x + width/2, thick_wall, width, label='Wall Thickness', color='#ed8936', edgecolor='#2d3748') axes[1, 2].set_xticks(x) @@ -1186,12 +1692,16 @@ def main(nc_file): axes[1, 2].set_title('Building Element Thickness', fontsize=11, fontweight='bold') axes[1, 2].legend(fontsize=8) + lbl_suffix = ' (domain mean)' if is_regional else '' + fig.suptitle(f'Urban Parameters{lbl_suffix}', fontsize=14, fontweight='bold') plt.tight_layout() pdf_path = os.path.join(pdf_dir, '07_urban_basic.pdf') fig.savefig(pdf_path, bbox_inches='tight') section7_figures.append({ 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Basic urban parameters for three density classes: Tall Building District, High Density, and Medium Density.', + 'caption': ('Basic urban parameters for three density classes (domain-mean values shown for regional grids).' + if is_regional else + 'Basic urban parameters for three density classes.'), 'base64': fig_to_base64(fig) }) plt.close(fig) @@ -1199,10 +1709,8 @@ def main(nc_file): # Urban radiative properties fig, axes = plt.subplots(2, 2, figsize=(14, 10)) - # Emissivities - em_vars = ['EM_IMPROAD', 'EM_PERROAD', 'EM_ROOF', 'EM_WALL'] - em_data = [get_1d_array(nc.variables[v]) for v in em_vars] em_labels = ['Impervious Road', 'Pervious Road', 'Roof', 'Wall'] + em_data = [em_improad, em_perroad, em_roof, em_wall] x = np.arange(len(URBAN_TYPES)) width = 0.2 @@ -1222,8 +1730,15 @@ def main(nc_file): alb_vars_dir = ['ALB_IMPROAD_DIR', 'ALB_PERROAD_DIR', 'ALB_ROOF_DIR', 'ALB_WALL_DIR'] alb_data_dir = [] for v in alb_vars_dir: - data = nc.variables[v][:] - alb_data_dir.append(data[0, :, 0, 0]) # VIS band + data = np.array(nc.variables[v][:], dtype=float) + if is_regional: + # domain mean over (numurbl, lat, lon): axis 0 = numrad, 1 = numurbl + alb_data_dir.append(np.array([ + float(np.nanmean(np.where(land_mask, data[0, i, :, :], np.nan))) + for i in range(len(URBAN_TYPES)) + ])) + else: + alb_data_dir.append(data[0, :, 0, 0]) # VIS band for i, (data, label, color) in enumerate(zip(alb_data_dir, em_labels, colors)): axes[0, 1].bar(x + (i - 1.5) * width, data, width, label=label, color=color, edgecolor='#2d3748') @@ -1239,8 +1754,14 @@ def main(nc_file): tk_labels = ['Roof', 'Wall', 'Impervious Road'] tk_data = [] for v in tk_vars: - data = nc.variables[v][:] - tk_data.append(data[0, :, 0, 0]) # First level + d = np.array(nc.variables[v][:], dtype=float) + if is_regional: + tk_data.append(np.array([ + float(np.nanmean(np.where(land_mask, d[0, i, :, :], np.nan))) + for i in range(len(URBAN_TYPES)) + ])) + else: + tk_data.append(d[0, :, 0, 0]) for i, (data, label) in enumerate(zip(tk_data, tk_labels)): axes[1, 0].bar(x + (i - 1) * 0.25, data, 0.25, label=label, @@ -1248,7 +1769,7 @@ def main(nc_file): axes[1, 0].set_xticks(x) axes[1, 0].set_xticklabels(URBAN_TYPES, rotation=15) - axes[1, 0].set_ylabel('Thermal Conductivity (W/m*K)') + axes[1, 0].set_ylabel('Thermal Conductivity (W/m\u00b7K)') axes[1, 0].set_title('Thermal Conductivity (Layer 1)', fontsize=11, fontweight='bold') axes[1, 0].legend(fontsize=8) @@ -1257,8 +1778,14 @@ def main(nc_file): cv_labels = ['Roof', 'Wall', 'Impervious Road'] cv_data = [] for v in cv_vars: - data = nc.variables[v][:] - cv_data.append(data[0, :, 0, 0] / 1e6) # Convert to MJ for readability + d = np.array(nc.variables[v][:], dtype=float) + if is_regional: + cv_data.append(np.array([ + float(np.nanmean(np.where(land_mask, d[0, i, :, :], np.nan))) + for i in range(len(URBAN_TYPES)) + ]) / 1e6) + else: + cv_data.append(d[0, :, 0, 0] / 1e6) # Convert to MJ for readability for i, (data, label) in enumerate(zip(cv_data, cv_labels)): axes[1, 1].bar(x + (i - 1) * 0.25, data, 0.25, label=label, @@ -1293,81 +1820,120 @@ def main(nc_file): print("Creating Section 8: Emission Factors and Other Parameters...") section8_figures = [] - fig, axes = plt.subplots(1, 2, figsize=(14, 6)) - - # Emission factors (isoprene) - ef_vars = ['EF1_BTR', 'EF1_FET', 'EF1_FDT', 'EF1_SHR', 'EF1_GRS', 'EF1_CRP'] + ef_vars = ['EF1_BTR', 'EF1_FET', 'EF1_FDT', 'EF1_SHR', 'EF1_GRS', 'EF1_CRP'] ef_labels = ['Broadleaf Trees', 'Fineleaf Evergreen Trees', 'Fineleaf Deciduous Trees', 'Shrubs', 'Grasses', 'Crops'] - ef_values = [get_scalar_value(nc.variables[v]) for v in ef_vars] - - colors = plt.cm.Greens(np.linspace(0.3, 0.9, len(ef_vars))) - axes[0].bar(ef_labels, ef_values, color=colors, edgecolor='#2d3748') - axes[0].set_ylabel('Emission Factor') - axes[0].set_title('Isoprene Emission Factors by Vegetation Type', fontsize=11, fontweight='bold') - axes[0].tick_params(axis='x', rotation=45) - axes[0].set_yscale('log') - for i, v in enumerate(ef_values): - axes[0].text(i, v * 1.1, f'{v:.0f}', ha='center', fontsize=9) - - # Glacier elevation classes - glc_mec = get_1d_array(nc.variables['GLC_MEC']) - pct_glc_mec = get_1d_array(nc.variables['PCT_GLC_MEC']) - topo_glc_mec = get_1d_array(nc.variables['TOPO_GLC_MEC']) - - ax2 = axes[1].twinx() - x = range(len(pct_glc_mec)) - axes[1].bar(x, pct_glc_mec, color='#63b3ed', alpha=0.7, label='Glacier %', edgecolor='#2d3748') - ax2.plot(x, topo_glc_mec, 'ro-', label='Mean Elevation', linewidth=2, markersize=8) - - axes[1].set_xticks(x) - axes[1].set_xticklabels([f'{glc_mec[i]:.0f}-{glc_mec[i+1]:.0f}m' for i in range(len(glc_mec)-1)], rotation=45) - axes[1].set_xlabel('Elevation Class') - axes[1].set_ylabel('Glacier Percent (%)', color='#3182ce') - ax2.set_ylabel('Mean Elevation (m)', color='red') - axes[1].set_title('Glacier Elevation Classes', fontsize=11, fontweight='bold') - axes[1].legend(loc='upper left') - ax2.legend(loc='upper right') - - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '08_emission_glacier.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section8_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Isoprene emission factors by vegetation type and glacier distribution across elevation classes.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) - - # Harvest parameters - fig, ax = plt.subplots(figsize=(10, 6)) - - harvest_vars = ['CONST_HARVEST_VH1', 'CONST_HARVEST_VH2', 'CONST_HARVEST_SH1', - 'CONST_HARVEST_SH2', 'CONST_HARVEST_SH3'] - harvest_labels = ['Primary Forest', 'Primary Non-Forest', 'Secondary Mature Forest', - 'Secondary Young Forest', 'Secondary Non-Forest'] - harvest_values = [get_scalar_value(nc.variables[v]) for v in harvest_vars] - colors = plt.cm.YlOrBr(np.linspace(0.3, 0.9, len(harvest_vars))) - ax.bar(harvest_labels, harvest_values, color=colors, edgecolor='#2d3748') - ax.set_ylabel('Harvest Rate (gC/m2/yr)') - ax.set_title('Constant Harvest Rates by Land Type', fontsize=12, fontweight='bold') - ax.tick_params(axis='x', rotation=30) - ax.grid(axis='y', alpha=0.3) - - for i, v in enumerate(harvest_values): - if v > 0: - ax.text(i, v + 10, f'{v:.1f}', ha='center', fontsize=9) + if is_regional: + # Spatial maps of emission factors + data_list = [get_field_2d(nc.variables[v]) for v in ef_vars] + fig = plot_map_grid(data_list, ef_labels, lons, lats, + 'Isoprene Emission Factors \u2013 Spatial Distribution', + units='unitless', cmap='YlGn', ncols=3) + pdf_path = os.path.join(pdf_dir, '08_emission_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section8_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Spatial maps of isoprene emission factors for the six vegetation functional types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Harvest maps + harvest_vars = ['CONST_HARVEST_VH1', 'CONST_HARVEST_VH2', 'CONST_HARVEST_SH1', + 'CONST_HARVEST_SH2', 'CONST_HARVEST_SH3'] + harvest_labels = ['Primary Forest', 'Primary Non-Forest', 'Sec. Mature Forest', + 'Sec. Young Forest', 'Sec. Non-Forest'] + harvest_data = [get_field_2d(nc.variables[v]) for v in harvest_vars] + fig = plot_map_grid(harvest_data, harvest_labels, lons, lats, + 'Constant Harvest Rates \u2013 Spatial Distribution', + units='gC/m\u00b2/yr', cmap='YlOrBr', ncols=3) + pdf_path = os.path.join(pdf_dir, '08b_harvest_maps.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section8_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Spatial maps of constant harvest rates for different land cover types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + # Glacier elevation (domain-mean of PCT_GLC_MEC, GLC_MEC is 1-D) + glc_mec = get_1d_array(nc.variables['GLC_MEC']) + pct_glc_data = np.array(nc.variables['PCT_GLC_MEC'][:], dtype=float) + topo_glc_data = np.array(nc.variables['TOPO_GLC_MEC'][:], dtype=float) + pct_glc_mec = np.array([ + float(np.nanmean(np.where(land_mask, pct_glc_data[i], np.nan))) + for i in range(pct_glc_data.shape[0]) + ]) + topo_glc_mec = np.array([ + float(np.nanmean(np.where(land_mask, topo_glc_data[i], np.nan))) + for i in range(topo_glc_data.shape[0]) + ]) - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '08b_harvest.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section8_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Constant harvest rates for different land cover types.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + else: + ef_values = [get_scalar_value(nc.variables[v]) for v in ef_vars] + + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + colors = plt.cm.Greens(np.linspace(0.3, 0.9, len(ef_vars))) + axes[0].bar(ef_labels, ef_values, color=colors, edgecolor='#2d3748') + axes[0].set_ylabel('Emission Factor') + axes[0].set_title('Isoprene Emission Factors by Vegetation Type', fontsize=11, fontweight='bold') + axes[0].tick_params(axis='x', rotation=45) + axes[0].set_yscale('log') + for i, v in enumerate(ef_values): + axes[0].text(i, v * 1.1, f'{v:.0f}', ha='center', fontsize=9) + + glc_mec = get_1d_array(nc.variables['GLC_MEC']) + pct_glc_mec = get_1d_array(nc.variables['PCT_GLC_MEC']) + topo_glc_mec = get_1d_array(nc.variables['TOPO_GLC_MEC']) + + ax2 = axes[1].twinx() + x_glc = range(len(pct_glc_mec)) + axes[1].bar(x_glc, pct_glc_mec, color='#63b3ed', alpha=0.7, label='Glacier %', edgecolor='#2d3748') + ax2.plot(x_glc, topo_glc_mec, 'ro-', label='Mean Elevation', linewidth=2, markersize=8) + axes[1].set_xticks(x_glc) + axes[1].set_xticklabels([f'{glc_mec[i]:.0f}-{glc_mec[i+1]:.0f}m' for i in range(len(glc_mec)-1)], rotation=45) + axes[1].set_xlabel('Elevation Class') + axes[1].set_ylabel('Glacier Percent (%)', color='#3182ce') + ax2.set_ylabel('Mean Elevation (m)', color='red') + axes[1].set_title('Glacier Elevation Classes', fontsize=11, fontweight='bold') + axes[1].legend(loc='upper left') + ax2.legend(loc='upper right') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '08_emission_glacier.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section8_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Isoprene emission factors by vegetation type and glacier distribution across elevation classes.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + harvest_vars = ['CONST_HARVEST_VH1', 'CONST_HARVEST_VH2', 'CONST_HARVEST_SH1', + 'CONST_HARVEST_SH2', 'CONST_HARVEST_SH3'] + harvest_labels = ['Primary Forest', 'Primary Non-Forest', 'Secondary Mature Forest', + 'Secondary Young Forest', 'Secondary Non-Forest'] + harvest_values = [get_scalar_value(nc.variables[v]) for v in harvest_vars] + + fig, ax = plt.subplots(figsize=(10, 6)) + colors = plt.cm.YlOrBr(np.linspace(0.3, 0.9, len(harvest_vars))) + ax.bar(harvest_labels, harvest_values, color=colors, edgecolor='#2d3748') + ax.set_ylabel('Harvest Rate (gC/m\u00b2/yr)') + ax.set_title('Constant Harvest Rates by Land Type', fontsize=12, fontweight='bold') + ax.tick_params(axis='x', rotation=30) + ax.grid(axis='y', alpha=0.3) + for i, v in enumerate(harvest_values): + if v > 0: + ax.text(i, v + 10, f'{v:.1f}', ha='center', fontsize=9) + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '08b_harvest.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section8_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Constant harvest rates for different land cover types.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'other', @@ -1382,88 +1948,121 @@ def main(nc_file): print("Creating Section 9: Summary Overview...") section9_figures = [] - # Create a comprehensive summary figure - fig = plt.figure(figsize=(16, 12)) - - # Land cover pie (small) - ax1 = fig.add_subplot(2, 3, 1) - nonzero_land = land_values > 0.1 - ax1.pie(land_values[nonzero_land], labels=[l for l, m in zip(land_labels, nonzero_land) if m], - autopct='%1.1f%%', colors=plt.cm.Set3(np.linspace(0, 1, sum(nonzero_land)))) - ax1.set_title('Land Cover Types', fontsize=11, fontweight='bold') - - # Soil texture summary - ax2 = fig.add_subplot(2, 3, 2) - mean_sand = np.mean(pct_sand) - mean_clay = np.mean(pct_clay) - mean_silt = 100 - mean_sand - mean_clay - ax2.pie([mean_sand, mean_silt, mean_clay], labels=['Sand', 'Silt', 'Clay'], - autopct='%1.1f%%', colors=['#f6ad55', '#a0aec0', '#fc8181']) - ax2.set_title('Mean Soil Texture', fontsize=11, fontweight='bold') - - # Active PFTs summary - ax3 = fig.add_subplot(2, 3, 3) - active_labels = [all_pft_names.get(p, f'PFT {p}')[:20] for p in active_pfts[:8]] - active_lai_mean = [np.mean(lai_data[:, p, 0, 0]) for p in active_pfts[:8]] - ax3.barh(range(len(active_labels)), active_lai_mean, color=plt.cm.Greens(np.linspace(0.3, 0.9, len(active_labels)))) - ax3.set_yticks(range(len(active_labels))) - ax3.set_yticklabels(active_labels, fontsize=8) - ax3.set_xlabel('Mean LAI') - ax3.set_title('Top Active PFTs by LAI', fontsize=11, fontweight='bold') - - # Key scalars table - ax4 = fig.add_subplot(2, 3, (4, 6)) - ax4.axis('off') - - key_params = [ - ('Longitude', f'{lon:.3f} E'), - ('Latitude', f'{lat:.3f} N'), - ('Natural Vegetation', f'{pct_natveg:.1f}%'), - ('Cropland', f'{pct_crop:.1f}%'), - ('Urban', f'{np.sum(pct_urban):.2f}%'), - ('Bedrock Depth', f'{get_scalar_value(nc.variables["zbedrock"]):.2f} m'), - ('Slope', f'{get_scalar_value(nc.variables["SLOPE"]):.2f} deg'), - ('Soil Color Index', f'{int(get_scalar_value(nc.variables["SOIL_COLOR"]))}'), - ('Max Saturated Fraction', f'{get_scalar_value(nc.variables["FMAX"]):.3f}'), - ('Peatland Fraction', f'{get_scalar_value(nc.variables["peatf"]):.3f}'), - ] - - table_data = [[p[0], p[1]] for p in key_params] - table = ax4.table(cellText=table_data, colLabels=['Parameter', 'Value'], - loc='center', cellLoc='left', colWidths=[0.4, 0.3]) - table.auto_set_font_size(False) - table.set_fontsize(11) - table.scale(1.2, 1.8) - - # Style the table - for i in range(len(table_data) + 1): - for j in range(2): - cell = table[(i, j)] - if i == 0: - cell.set_facecolor('#667eea') - cell.set_text_props(color='white', fontweight='bold') - elif i % 2 == 0: - cell.set_facecolor('#f7fafc') - else: - cell.set_facecolor('#edf2f7') - - ax4.set_title('Key Site Parameters Summary', fontsize=12, fontweight='bold', pad=20) - - fig.suptitle('Surface Data Overview - BE-Lon Site', fontsize=16, fontweight='bold') - plt.tight_layout() - pdf_path = os.path.join(pdf_dir, '09_summary.pdf') - fig.savefig(pdf_path, bbox_inches='tight') - section9_figures.append({ - 'pdf_name': os.path.basename(pdf_path), - 'caption': 'Comprehensive summary of key surface data parameters for the BE-Lon site.', - 'base64': fig_to_base64(fig) - }) - plt.close(fig) + nc_basename_stem = os.path.splitext(os.path.basename(nc_file))[0] + + if is_regional: + # Grid of key spatial maps + summary_vars = [ + ('PCT_NATVEG', 'Natural Vegetation (%)', 'Greens'), + ('PCT_CROP', 'Cropland (%)', 'YlOrBr'), + ('SOIL_COLOR', 'Soil Color Index', 'tab20b'), + ('FMAX', 'Max Saturated Fraction', 'Blues'), + ('zbedrock', 'Bedrock Depth (m)', 'copper'), + ('peatf', 'Peatland Fraction', 'Greens'), + ] + data_list, titles_list, cmaps_list = [], [], [] + for var_name, title_str, cmap_str in summary_vars: + if var_name in nc.variables: + d = get_field_2d(nc.variables[var_name]) + if d.ndim == 2: + data_list.append(d) + titles_list.append(title_str) + cmaps_list.append(cmap_str) + fig = plot_map_grid(data_list, titles_list, lons, lats, + f'Key Parameters Overview \u2013 {nc_basename_stem}', + units='', cmap=cmaps_list, ncols=3) + pdf_path = os.path.join(pdf_dir, '09_summary.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section9_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Spatial overview of key surface data parameters.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) + + else: + fig = plt.figure(figsize=(16, 12)) + + ax1 = fig.add_subplot(2, 3, 1) + land_labels = ['Natural Vegetation', 'Cropland', 'Urban', 'Lake', 'Wetland', 'Glacier'] + land_values = np.array([ + pct_natveg, pct_crop, np.sum(pct_urban), + get_scalar_value(nc.variables['PCT_LAKE']), + get_scalar_value(nc.variables['PCT_WETLAND']), + get_scalar_value(nc.variables['PCT_GLACIER']) + ]) + nonzero_land = land_values > 0.1 + ax1.pie(land_values[nonzero_land], + labels=[l for l, m in zip(land_labels, nonzero_land) if m], + autopct='%1.1f%%', colors=plt.cm.Set3(np.linspace(0, 1, sum(nonzero_land)))) + ax1.set_title('Land Cover Types', fontsize=11, fontweight='bold') + + ax2 = fig.add_subplot(2, 3, 2) + mean_sand = float(np.mean(pct_sand)) + mean_clay = float(np.mean(pct_clay)) + mean_silt = 100 - mean_sand - mean_clay + ax2.pie([mean_sand, mean_silt, mean_clay], labels=['Sand', 'Silt', 'Clay'], + autopct='%1.1f%%', colors=['#f6ad55', '#a0aec0', '#fc8181']) + ax2.set_title('Mean Soil Texture', fontsize=11, fontweight='bold') + + ax3 = fig.add_subplot(2, 3, 3) + active_labels = [all_pft_names.get(p, f'PFT {p}')[:20] for p in active_pfts[:8]] + active_lai_mean = [float(np.mean(lai_data[:, p, 0, 0])) for p in active_pfts[:8]] + ax3.barh(range(len(active_labels)), active_lai_mean, + color=plt.cm.Greens(np.linspace(0.3, 0.9, max(len(active_labels), 1)))) + ax3.set_yticks(range(len(active_labels))) + ax3.set_yticklabels(active_labels, fontsize=8) + ax3.set_xlabel('Mean LAI') + ax3.set_title('Top Active PFTs by LAI', fontsize=11, fontweight='bold') + + ax4 = fig.add_subplot(2, 3, (4, 6)) + ax4.axis('off') + key_params = [ + ('Longitude', f'{lon:.3f}\u00b0E'), + ('Latitude', f'{lat:.3f}\u00b0N'), + ('Natural Vegetation', f'{pct_natveg:.1f}%'), + ('Cropland', f'{pct_crop:.1f}%'), + ('Urban (total)', f'{float(np.sum(pct_urban)):.2f}%'), + ('Bedrock Depth', f'{get_scalar_value(nc.variables["zbedrock"]):.2f} m'), + ('Slope', f'{get_scalar_value(nc.variables["SLOPE"]):.2f}\u00b0'), + ('Soil Color Index', f'{int(get_scalar_value(nc.variables["SOIL_COLOR"]))}'), + ('Max Saturated Frac.', f'{get_scalar_value(nc.variables["FMAX"]):.3f}'), + ('Peatland Fraction', f'{get_scalar_value(nc.variables["peatf"]):.3f}'), + ] + table_data = [[p[0], p[1]] for p in key_params] + table = ax4.table(cellText=table_data, colLabels=['Parameter', 'Value'], + loc='center', cellLoc='left', colWidths=[0.4, 0.3]) + table.auto_set_font_size(False) + table.set_fontsize(11) + table.scale(1.2, 1.8) + for i in range(len(table_data) + 1): + for j in range(2): + cell = table[(i, j)] + if i == 0: + cell.set_facecolor('#667eea') + cell.set_text_props(color='white', fontweight='bold') + elif i % 2 == 0: + cell.set_facecolor('#f7fafc') + else: + cell.set_facecolor('#edf2f7') + ax4.set_title('Key Site Parameters Summary', fontsize=12, fontweight='bold', pad=20) + + fig.suptitle(f'Surface Data Overview \u2013 {nc_basename_stem}', + fontsize=16, fontweight='bold') + plt.tight_layout() + pdf_path = os.path.join(pdf_dir, '09_summary.pdf') + fig.savefig(pdf_path, bbox_inches='tight') + section9_figures.append({ + 'pdf_name': os.path.basename(pdf_path), + 'caption': 'Comprehensive summary of key surface data parameters.', + 'base64': fig_to_base64(fig) + }) + plt.close(fig) figures_data.append({ 'id': 'summary', 'title': 'Summary Overview', - 'description': 'A comprehensive overview combining key parameters for quick reference during discussions.', + 'description': 'A comprehensive overview combining key parameters for quick reference.', 'figures': section9_figures }) @@ -1472,7 +2071,7 @@ def main(nc_file): # Generate HTML report print("Generating HTML report...") - html_file = create_html_report(figures_data, nc_file, output_dir) + html_file = create_html_report(figures_data, nc_file, output_dir, metadata_str) print(f"\nVisualization complete!") print(f"PDF figures saved in: {pdf_dir}") From c2a3007e9b3f6da3eba79d5c3cca65cc06c33d7c Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 5 Mar 2026 17:09:40 +0100 Subject: [PATCH 3/3] vizsurfdata: typo fix in script name --- vizsurfdata/{vizualize_surfdata.py => visualize_surfdata.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename vizsurfdata/{vizualize_surfdata.py => visualize_surfdata.py} (100%) diff --git a/vizsurfdata/vizualize_surfdata.py b/vizsurfdata/visualize_surfdata.py similarity index 100% rename from vizsurfdata/vizualize_surfdata.py rename to vizsurfdata/visualize_surfdata.py