Understanding Agricultural Land Cover Changes

admin_id = "EL64"  # Central Greece

import xarray as xr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd
import regionmask
from pathlib import Path
import rioxarray as rxr
import os

# Set up data directories
data_dir = Path("../data")
corine_output_dir = data_dir / admin_id / "corine_land_cover"
satellite_output_dir = data_dir / admin_id / "satellite_land_cover"

print(f"\nData directories:")
print(f"  CORINE: {corine_output_dir}")
print(f"  Satellite: {satellite_output_dir}")

# Read NUTS shapefiles
regions_dir = data_dir / 'regions'
nuts_shp = regions_dir / 'NUTS_RG_20M_2024_4326' / 'NUTS_RG_20M_2024_4326.shp'
nuts_gdf = gpd.read_file(nuts_shp)

# Select the region of interest
sel_gdf = nuts_gdf[nuts_gdf['NUTS_ID'] == admin_id]
print(f"Region: {sel_gdf['NUTS_NAME'].values[0]} ({admin_id})")
print(f"Country: {sel_gdf['CNTR_CODE'].values[0]}")
lon_min, lat_min, lon_max, lat_max = sel_gdf.geometry.total_bounds

# Create a regionmask from the admin region geometry
admin_mask = regionmask.from_geopandas(sel_gdf, names='NUTS_ID')

# Load CORINE CSV files
csv_corine_categories = corine_output_dir / "land_cover_categories.csv"
csv_corine_agri = corine_output_dir / "agricultural_land_breakdown.csv"

# Read the data
df_corine_categories = pd.read_csv(csv_corine_categories)
df_corine_agri = pd.read_csv(csv_corine_agri)

print("CORINE land cover categories:")
print(f"  Years: {df_corine_categories['year'].min():.0f} - {df_corine_categories['year'].max():.0f}")
print(f"  Records: {len(df_corine_categories)}")
print(f"\nMain categories sample:")
print(df_corine_categories.head())
print(f"\nAgricultural breakdown sample:")
print(df_corine_agri.head())

# Load satellite CSV file
csv_satellite = satellite_output_dir / f"agricultural_land_fraction_{admin_id}.csv"

# Read the data
df_satellite = pd.read_csv(csv_satellite)
# Convert fraction to percentage for easier comparison
df_satellite['agricultural'] = df_satellite['ag_fraction'] * 100

print("Satellite land cover data:")
print(f"  Years: {df_satellite['year'].min():.0f} - {df_satellite['year'].max():.0f}")
print(f"  Records: {len(df_satellite)}")
print(f"\nSample:")
print(df_satellite.head())

fig, ax = plt.subplots(figsize=(14, 6))

# Plot CORINE agricultural land
ax.plot(df_corine_categories['year'], df_corine_categories['agricultural'], 
        marker='o', linewidth=2.5, markersize=10, label='CORINE Land Cover', 
        color='#8B4513', linestyle='-', alpha=0.8)

# Plot satellite agricultural land
ax.plot(df_satellite['year'], df_satellite['agricultural'], 
        marker='s', linewidth=2, markersize=6, label='Satellite Land Cover', 
        color='#228B22', linestyle='--', alpha=0.8)

ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Agricultural Land Coverage (%)', fontsize=12)
ax.set_title(f'Agricultural Land Change in {admin_id} - Dataset Comparison', fontsize=14, fontweight='bold')
ax.legend(fontsize=11, loc='upper left')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Find overlapping years (CORINE years within satellite range)
corine_years = df_corine_categories['year'].values
satellite_years = df_satellite['year'].values

# Get closest satellite values for each CORINE year
comparison_data = []
for corine_year in corine_years:
    # Find exact match or closest year in satellite data
    if corine_year in satellite_years:
        sat_val = df_satellite[df_satellite['year'] == corine_year]['agricultural'].values[0]
        cor_val = df_corine_categories[df_corine_categories['year'] == corine_year]['agricultural'].values[0]
        
        comparison_data.append({
            'year': corine_year,
            'corine': cor_val,
            'satellite': sat_val,
            'difference': sat_val - cor_val,
            'percent_diff': ((sat_val - cor_val) / cor_val) * 100
        })

df_comparison = pd.DataFrame(comparison_data)

print("Dataset Comparison (overlapping years):")
print("="*90)
print(f"{'Year':<8} {'CORINE (%)':>12} {'Satellite (%)':>14} {'Diff (pp)':>12} {'Diff (%)':>12}")
print("-"*90)

for _, row in df_comparison.iterrows():
    print(f"{row['year']:<8.0f} {row['corine']:>12.2f} {row['satellite']:>14.2f} "
          f"{row['difference']:>12.2f} {row['percent_diff']:>11.2f}%")

print("\nSummary Statistics:")
print(f"  Mean absolute difference: {df_comparison['difference'].abs().mean():.2f} percentage points")
print(f"  Mean relative difference: {df_comparison['percent_diff'].abs().mean():.2f}%")
print(f"  Max difference: {df_comparison['difference'].abs().max():.2f} pp in {df_comparison.loc[df_comparison['difference'].abs().idxmax(), 'year']:.0f}")

fig, ax = plt.subplots(figsize=(14, 6))

ax.plot(df_corine_categories['year'], df_corine_categories['agricultural'], 
        marker='o', linewidth=2.5, markersize=8, label='Agricultural', color='#8B4513')
ax.plot(df_corine_categories['year'], df_corine_categories['forest'], 
        marker='s', linewidth=2.5, markersize=8, label='Forest', color='#228B22')
ax.plot(df_corine_categories['year'], df_corine_categories['urban'], 
        marker='^', linewidth=2.5, markersize=8, label='Urban', color='#808080')
ax.plot(df_corine_categories['year'], df_corine_categories['water'], 
        marker='d', linewidth=2.5, markersize=8, label='Water', color='#4169E1')

ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Land Cover (%)', fontsize=12)
ax.set_title(f'Land Cover Change Over Time in {admin_id}', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print trends
print("Land Cover Trends (1990 to 2018):")
for category in ['agricultural', 'forest', 'urban', 'water']:
    start_val = df_corine_categories[category].iloc[0]
    end_val = df_corine_categories[category].iloc[-1]
    change = end_val - start_val
    print(f"  {category.capitalize()}: {start_val:.2f}% → {end_val:.2f}% (change: {change:+.2f} pp)")

# Stacked area plot of agricultural land types
fig, ax = plt.subplots(figsize=(14, 7))

# Prepare data for stacking
agri_cols = [col for col in df_corine_agri.columns if col != 'year']
years = df_corine_agri['year'].values

# Create stacked area plot
ax.stackplot(years, 
             *[df_corine_agri[col].values for col in agri_cols],
             labels=agri_cols,
             alpha=0.8)

ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Agricultural Land Coverage (%)', fontsize=12)
ax.set_title(f'Agricultural Land Types Over Time in {admin_id}', fontsize=14, fontweight='bold')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize=9)
ax.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

# Print summary of agricultural types
print("\nAgricultural land types (2018):")
latest_agri = df_corine_agri.iloc[-1]
for col in agri_cols:
    val = latest_agri[col]
    if val > 0.5:  # Only show significant types
        print(f"  {col}: {val:.2f}%")

# Install folium if not available
try:
    import folium
    from folium import plugins
except ImportError:
    !pip install folium
    import folium
    from folium import plugins

# Load CORINE land cover raster (latest year)
corine_tif = corine_output_dir / f"corine_landcover_{admin_id}.tif"
corine_raster = rxr.open_rasterio(corine_tif)
corine_raster = corine_raster.squeeze()

print(f"Loaded CORINE raster: {corine_tif.name}")
print(f"Shape: {corine_raster.shape}")
print(f"Coordinate ranges: lon=[{float(corine_raster.x.min()):.2f}, {float(corine_raster.x.max()):.2f}], lat=[{float(corine_raster.y.min()):.2f}, {float(corine_raster.y.max()):.2f}]")
print(f"Unique classes present: {np.unique(corine_raster.values[~np.isnan(corine_raster.values)])}")

# Define CORINE class names and colors
corine_classes = {
    1: ('Continuous urban fabric', '#E6004D'),
    2: ('Discontinuous urban fabric', '#FF0000'),
    3: ('Industrial or commercial units', '#CC4DF2'),
    4: ('Road and rail networks', '#CC0000'),
    5: ('Port areas', '#E6CCCC'),
    6: ('Airports', '#E6CCE6'),
    7: ('Mineral extraction sites', '#A600CC'),
    8: ('Dump sites', '#A64DCC'),
    9: ('Construction sites', '#FF4DFF'),
    10: ('Green urban areas', '#FFA6FF'),
    11: ('Sport and leisure facilities', '#FFE6FF'),
    12: ('Non-irrigated arable land', '#FFFFA8'),
    13: ('Permanently irrigated land', '#FFFF00'),
    14: ('Rice fields', '#E6E600'),
    15: ('Vineyards', '#E68000'),
    16: ('Fruit trees and berry plantations', '#F2A64D'),
    17: ('Olive groves', '#E6A600'),
    18: ('Pastures', '#E6E64D'),
    19: ('Annual crops associated with permanent crops', '#FFE6A6'),
    20: ('Complex cultivation patterns', '#FFE6CC'),
    21: ('Land principally occupied by agriculture', '#F2CC80'),
    22: ('Agro-forestry areas', '#80FF00'),
    23: ('Broad-leaved forest', '#00A600'),
    24: ('Coniferous forest', '#4DFF00'),
    25: ('Mixed forest', '#CCF24D'),
    26: ('Natural grasslands', '#A6FF80'),
    27: ('Moors and heathland', '#A6E64D'),
    28: ('Sclerophyllous vegetation', '#A6F200'),
    29: ('Transitional woodland-shrub', '#00CCF2'),
    30: ('Beaches, dunes, sands', '#CCE6FF'),
    31: ('Bare rocks', '#4D4DFF'),
    32: ('Sparsely vegetated areas', '#CCCCCC'),
    33: ('Burnt areas', '#000000'),
    34: ('Glaciers and perpetual snow', '#4DFF4D'),
    35: ('Inland marshes', '#00CCF2'),
    36: ('Peat bogs', '#80F2E6'),
    37: ('Salt marshes', '#00FFA6'),
    38: ('Salines', '#A6FFE6'),
    39: ('Intertidal flats', '#E6F2FF'),
    40: ('Water courses', '#00CCF2'),
    41: ('Water bodies', '#80F2E6'),
    42: ('Coastal lagoons', '#00FFA6'),
    43: ('Estuaries', '#A6FFE6'),
    44: ('Sea and ocean', '#E6F2FF')
}

# Calculate center of region
center_lat = (lat_min + lat_max) / 2
center_lon = (lon_min + lon_max) / 2

# Create base map
m = folium.Map(
    location=[center_lat, center_lon],
    zoom_start=9,
    tiles=None
)

# Add basemap options
folium.TileLayer('OpenStreetMap', name='OpenStreetMap').add_to(m)
folium.TileLayer(
    tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
    attr='Esri World Imagery',
    name='Satellite',
    overlay=False,
    control=True
).add_to(m)
folium.TileLayer(
    tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{z}/{y}/{x}',
    attr='Esri World Topo',
    name='Topographic',
    overlay=False,
    control=True
).add_to(m)

# Add region boundary
folium.GeoJson(
    sel_gdf.geometry,
    name=f'Region {admin_id}',
    style_function=lambda x: {
        'fillColor': 'none',
        'color': 'red',
        'weight': 3,
        'fillOpacity': 0
    }
).add_to(m)

# Create a simplified overlay showing agricultural vs non-agricultural land
# Agricultural classes are 12-22
corine_array = corine_raster.values
y_coords = corine_raster.y.values
x_coords = corine_raster.x.values

# Create points for agricultural land (green) and other land (various colors)
agri_points = []
other_points = []

# Sample the data for performance (every 5th point)
step = 5
for i in range(0, len(y_coords), step):
    for j in range(0, len(x_coords), step):
        val = corine_array[i, j]
        if not np.isnan(val):
            val_int = int(val)
            if 12 <= val_int <= 22:  # Agricultural
                agri_points.append([y_coords[i], x_coords[j]])
            else:
                other_points.append([y_coords[i], x_coords[j]])

# Add heatmaps for agricultural land
if len(agri_points) > 0:
    plugins.HeatMap(
        agri_points,
        name='Agricultural Land',
        min_opacity=0.4,
        max_zoom=13,
        radius=12,
        blur=15,
        gradient={0.4: 'yellow', 0.6: 'orange', 0.8: '#8B4513', 1.0: '#654321'}
    ).add_to(m)

# Add layer control
folium.LayerControl().add_to(m)

# Add title
title_html = f'''
<div style="position: fixed; 
     top: 10px; left: 50px; width: 450px; height: 70px; 
     background-color: white; border:2px solid grey; z-index:9999; 
     font-size:16px; padding: 10px; box-shadow: 2px 2px 6px rgba(0,0,0,0.3);
     border-radius: 5px;">
     <b>CORINE Land Cover - Agricultural Land Distribution</b><br>
     <span style="font-size:14px;">{admin_id} - Latest available (2018)</span><br>
     <i style="font-size:12px; color:#666;">Toggle layers and zoom to explore</i>
</div>
'''
m.get_root().html.add_child(folium.Element(title_html))

# Display map
m