# Library
import os
import geopandas as gpd
import rasterio
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import rasterio.mask
from rasterio.plot import show
from rasterio.warp import transform_geom

# Define file path
shapefile_path = r'\GIS\MartapuraBasin_WGS.shp'
lulc_folder = r'\Modelling\LULC\RAW'
output_excel = r'\Modelling\LULC\LULC_Analysis.xlsx'
output_image = r'\Modelling\LULC'

# Define land use categories and runoff coefficients
land_use_classes = {
    1: ('Evergreen Needleleaf Forests', 0.30, '#05450a'),
    2: ('Evergreen Broadleaf Forests', 0.25, '#086a10'),
    3: ('Deciduous Needleleaf Forests', 0.35, '#54a708'),
    4: ('Deciduous Broadleaf Forests', 0.30, '#78d203'),
    5: ('Mixed Forests', 0.30, '#009900'),
    6: ('Closed Shrublands', 0.45, '#c6b044'),
    7: ('Open Shrublands', 0.40, '#dcd159'),
    8: ('Woody Savannas', 0.20, '#dade48'),
    9: ('Savannas', 0.20, '#fbff13'),
    10: ('Grasslands', 0.40, '#b6ff05'),
    11: ('Permanent Wetlands', 0.60, '#27ff87'),
    12: ('Croplands', 0.75, '#c24f44'),
    13: ('Urban and Built-Up', 0.90, '#a5a5a5'),
    14: ('Cropland/Natural Vegetation Mosaics', 0.60, '#ff6d4c'),
    15: ('Snow and Ice', 0.05, '#69fff8'),
    16: ('Barren', 0.05, '#f9ffa4'),
    17: ('Water Bodies', 0.05, '#1c0dff'),
    0: ('NoData', 0.00, '#ffffff')  # NoData (to be excluded)
}

# Load shapefile
shapefile = gpd.read_file(shapefile_path)

# Function to reproject shapefile to match raster CRS
def reproject_shapefile(shapefile, raster_crs):
    return shapefile.to_crs(raster_crs)

# Function to clip raster
def clip_raster(raster_path, shapes):
    with rasterio.open(raster_path) as src:
        raster_crs = src.crs  # Get raster CRS
        
        # Reproject shapefile if necessary
        reprojected_shapefile = reproject_shapefile(shapefile, raster_crs)
        
        # Transform shapes to match raster CRS
        shapes_transformed = [transform_geom(reprojected_shapefile.crs, raster_crs, geom.__geo_interface__) for geom in reprojected_shapefile.geometry]
        
        # Clip raster
        clipped, transform = rasterio.mask.mask(src, shapes_transformed, crop=True)
        meta = src.meta.copy()
        meta.update({"height": clipped.shape[1], "width": clipped.shape[2], "transform": transform})
    return clipped[0], meta

# Initialize data storage
results = []

# Process each LULC file
for file in sorted(os.listdir(lulc_folder)):
    if file.endswith('.tif'):
        year = int(''.join(filter(str.isdigit, file)))  # Extract numeric year
        file_path = os.path.join(lulc_folder, file)
        
        # Clip raster
        lulc_data, meta = clip_raster(file_path, shapefile.geometry)
        
        # Exclude NoData (0) values
        valid_mask = lulc_data != 0  # Create a mask for valid values
        lulc_data = lulc_data[valid_mask]  # Apply mask
        
        if lulc_data.size == 0:  # Skip if all values are NoData
            continue
        
        # Compute area and percentage
        unique, counts = np.unique(lulc_data, return_counts=True)
        total_area = sum(counts)
        
        for code, count in zip(unique, counts):
            if code in land_use_classes:
                land_use_name, rc, color = land_use_classes[code]
                area_km2 = (count * meta['transform'][0] * -meta['transform'][4]) / 1e6  # Convert to sq km
                percentage = (count / total_area) * 100
                weighted_rc = area_km2 * rc
                results.append([year, code, land_use_name, area_km2, percentage, rc, weighted_rc, color])

# Convert to DataFrame
results_df = pd.DataFrame(results, columns=['Year', 'Category', 'Land Use Name', 'Area (sq km)', 'Percentage (%)', 'Runoff Coefficient', 'Weighted RC', 'Color'])

# Corrected way to compute Total Average RC per Year
total_avg_rc = results_df.groupby('Year')[['Weighted RC', 'Area (sq km)']].apply(lambda x: x['Weighted RC'].sum() / x['Area (sq km)'].sum()).rename('Total Average RC')

# Merge Total Average RC into DataFrame
results_df = results_df.merge(total_avg_rc, on='Year')

# Calculate Impervious Value (1 if RC >= 0.8, else 0)
results_df['Impervious Value'] = results_df['Runoff Coefficient'].apply(lambda x: 1 if x >= 0.8 else 0)

# Remove Color column before saving
results_df = results_df.drop(columns=['Color'])

# Save to Excel
results_df.to_excel(output_excel, index=False)

print(f"LULC analysis completed. Results saved to: {output_excel}")

# Plot LULC changes over time
fig, axes = plt.subplots(6, 4, figsize=(20, 20))  # Adjust grid size as needed
axes = axes.ravel()

for idx, year in enumerate(range(2001, 2024)):
    file_path = os.path.join(lulc_folder, f'LULC_{year}.tif')
    if os.path.exists(file_path):
        lulc_data, _ = clip_raster(file_path, shapefile.geometry)
        color_map = np.zeros((lulc_data.shape[0], lulc_data.shape[1], 3), dtype=np.uint8)
        for code, (_, _, color) in land_use_classes.items():
            mask = lulc_data == code
            rgb = tuple(int(color[i:i+2], 16) for i in (1, 3, 5))
            color_map[mask] = rgb
        
        axes[idx].imshow(color_map)
        axes[idx].set_title(f'LULC {year}')
        axes[idx].axis('on')
        
plt.tight_layout()
plt.savefig(output_image, dpi=300)
plt.show()

print("LULC analysis completed. Visualization generated.")

years_to_plot = [2014, 2020, 2021, 2022]
fig, axes_lulc = plt.subplots(1, 4, figsize=(20, 5))

for i, year_plot in enumerate(years_to_plot):
    file_path = os.path.join(lulc_folder, f'LULC_{year_plot}.tif')
    if os.path.exists(file_path):
        lulc_data, _ = clip_raster(file_path, shapefile.geometry)
        color_map = np.zeros((*lulc_data.shape, 3), dtype=np.uint8)
        for code, (_, _, color) in land_use_classes.items():
            mask = lulc_data == code
            rgb = tuple(int(color[i:i+2], 16) for i in (1, 3, 5))
            color_map[mask] = rgb
        axes_lulc[i].imshow(color_map)
        axes_lulc[i].set_title(f'LULC {year_plot}')
        axes_lulc[i].axis('off')
    else:
        axes_lulc[i].set_visible(False)

plt.tight_layout()
plt.show()

import pandas as pd
import matplotlib.pyplot as plt

# Convert results to DataFrame
df = pd.DataFrame(results, columns=['Year', 'Category', 'Land Use Name', 'Area (sq km)', 'Percentage', 'Runoff Coefficient', 'Weighted RC', 'Color'])

# Pivot the DataFrame for trend analysis
trend_df = df.pivot(index='Year', columns='Land Use Name', values='Area (sq km)')

# Create the figure and axes
fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=True)

# --- Plot 1: Overall Land Use Change Trends ---
axes[0].set_title('Overall Land Use Change Trends on Martapura Watershed (2001-2023)')
for category in trend_df.columns:
    axes[0].plot(trend_df.index, trend_df[category], label=category, marker='o')

axes[0].set_ylabel('Area (sq km)')
axes[0].legend(title='Land Use Categories', bbox_to_anchor=(1.05, 1), loc='upper left')
axes[0].grid(True)

# --- Plot 2: Urban Areas ---
special_categories = ['Urban and Built-Up']
trend_special = trend_df[special_categories]  # Filter only these categories

axes[1].set_title('Urban Area Trends on Martapura Watershed (2001-2023)')
for category in special_categories:
    axes[1].plot(trend_special.index, trend_special[category], label=category, marker='o')

axes[1].set_xlabel('Year')
axes[1].set_ylabel('Area (sq km)')
axes[1].legend(title='Land Use Categories', bbox_to_anchor=(1.05, 1), loc='upper left')
axes[1].grid(True)

# --- Plot 3: Croplands and Wetlands ---
special_categories = ['Croplands','Permanent Wetlands']
trend_special = trend_df[special_categories]  # Filter only these categories

axes[2].set_title('Croplands and Wetlands Trends on Martapura Watershed (2001-2023)')
for category in special_categories:
    axes[2].plot(trend_special.index, trend_special[category], label=category, marker='o')

axes[2].set_xlabel('Year')
axes[2].set_ylabel('Area (sq km)')
axes[2].legend(title='Land Use Categories', bbox_to_anchor=(1.05, 1), loc='upper left')
axes[2].grid(True)

# Adjust layout and show the plot
plt.tight_layout()
plt.show()