Agricultural land from CORINE land cover

admin_id = "EL64"  # Example admin ID for Central Greece

import xarray as xr
import numpy as np
import pandas as pd
import geopandas as gpd
import regionmask
from pathlib import Path
import rioxarray as rxr
import os
import re

# Set up data directories
data_dir = Path("../data")
corine_dir = data_dir / "corine_landcover"

print(f"\nCORINE data directory: {corine_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"Found {admin_id} region: {sel_gdf['NUTS_NAME'].values[0]}")
print(f"Bounding box: {sel_gdf.geometry.total_bounds}")
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')

# CORINE Land Cover class names (official nomenclature)
corine_classes = {
    # Artificial surfaces (1-11)
    1: 'Continuous urban fabric', 2: 'Discontinuous urban fabric',
    3: 'Industrial or commercial units', 4: 'Road and rail networks',
    5: 'Port areas', 6: 'Airports', 7: 'Mineral extraction sites',
    8: 'Dump sites', 9: 'Construction sites', 10: 'Green urban areas',
    11: 'Sport and leisure facilities',
    # Agricultural areas (12-22)
    12: 'Non-irrigated arable land', 13: 'Permanently irrigated land',
    14: 'Rice fields', 15: 'Vineyards', 16: 'Fruit trees and berry plantations',
    17: 'Olive groves', 18: 'Pastures', 19: 'Annual crops associated with permanent crops',
    20: 'Complex cultivation patterns', 21: 'Land principally occupied by agriculture',
    22: 'Agro-forestry areas',
    # Forest and semi-natural areas (23-34)
    23: 'Broad-leaved forest', 24: 'Coniferous forest', 25: 'Mixed forest',
    26: 'Natural grasslands', 27: 'Moors and heathland', 28: 'Sclerophyllous vegetation',
    29: 'Transitional woodland-shrub', 30: 'Beaches, dunes, sands',
    31: 'Bare rocks', 32: 'Sparsely vegetated areas', 33: 'Burnt areas',
    34: 'Glaciers and perpetual snow',
    # Wetlands (35-39)
    35: 'Inland marshes', 36: 'Peat bogs', 37: 'Salt marshes',
    38: 'Salines', 39: 'Intertidal flats',
    # Water bodies (40-44)
    40: 'Water courses', 41: 'Water bodies', 42: 'Coastal lagoons',
    43: 'Estuaries', 44: 'Sea and ocean'
}

# Define land cover class groups
agricultural_classes = list(range(12, 23))  # Classes 12-22
forest_classes = list(range(23, 30))  # Classes 23-29
urban_classes = list(range(1, 12))  # Classes 1-11
water_classes = list(range(40, 45))  # Classes 40-44

print(f"Total CORINE classes: {len(corine_classes)}")
print(f"Agricultural classes (12-22): {len(agricultural_classes)} types")

# List all CORINE tif files
tif_files = list(corine_dir.glob("*/*/*.tif"))
print(f"Found {len(tif_files)} CORINE tif files:")
for f in tif_files:
    print(f"  {f.name}")

# Read all tif files into a single xarray dataset
datasets = []
years = []

for tif_file in tif_files:
    # Read using rioxarray
    ds = rxr.open_rasterio(tif_file)
    
    # Extract year from filename - look for pattern like CLC1990, CLC2000, etc.
    year_match = re.search(r'CLC(\d{4})', tif_file.name)
    if year_match:
        year = int(year_match.group(1))
        years.append(year)
        print(f"Processing year: {year}")
    else:
        print(f"Warning: Could not extract year from {tif_file.name}")
        continue
    
    datasets.append(ds)

# Sort datasets and years together by year
sorted_pairs = sorted(zip(years, datasets), key=lambda x: x[0])
years, datasets = zip(*sorted_pairs)

# Concatenate along new time dimension
corine_ds = xr.concat(datasets, dim="time")
corine_ds = corine_ds.rename({"band": "land_cover", "x": "lon", "y": "lat"})

# Assign time based on extracted years
corine_ds["time"] = pd.to_datetime([f"{year}-01-01" for year in years])
corine_ds = corine_ds.assign_coords(time=corine_ds.time)
corine_ds = corine_ds.squeeze("land_cover")  # remove land_cover dimension

# Replace fill value (-128) with NaN
corine_ds = corine_ds.where(corine_ds != -128)

print(f"\nCORINE dataset shape: {corine_ds.shape}")
print(f"Years available: {years}")
print(f"Coordinate ranges: lon=[{float(corine_ds.lon.min()):.2f}, {float(corine_ds.lon.max()):.2f}], lat=[{float(corine_ds.lat.min()):.2f}, {float(corine_ds.lat.max()):.2f}]")

# Calculate grid cell areas (accounting for latitude)
# Earth radius in meters
R = 6371000

# Get lat/lon spacing
lat = corine_ds.lat.values
lon = corine_ds.lon.values
dlat = np.abs(np.diff(lat).mean())
dlon = np.abs(np.diff(lon).mean())

# Calculate area for each latitude (in km²)
lat_rad = np.deg2rad(lat)
# Area = R² * dλ * (sin(φ₂) - sin(φ₁))
areas = R**2 * np.deg2rad(dlon) * np.abs(
    np.sin(lat_rad + np.deg2rad(dlat/2)) - 
    np.sin(lat_rad - np.deg2rad(dlat/2))
) / 1e6  # Convert to km²

# Create 2D area array
area_2d = np.tile(areas[:, np.newaxis], (1, len(lon)))

print(f"Grid cell areas range from {areas.min():.2f} to {areas.max():.2f} km²")
print(f"Total study area: {area_2d.sum():.2f} km²")

# Calculate land cover percentages for each year
def calculate_coverage(data_array, class_list, area_2d):
    """Calculate areal coverage for specified classes"""
    mask = np.isin(data_array.values, class_list)
    total_area = area_2d[~np.isnan(data_array.values)].sum()
    class_area = area_2d[mask & ~np.isnan(data_array.values)].sum()
    return (class_area / total_area) * 100

results = []
for i, year in enumerate(corine_ds.time.dt.year.values):
    data = corine_ds.isel(time=i)
    
    # Calculate percentages for major categories
    agri_pct = calculate_coverage(data, agricultural_classes, area_2d)
    forest_pct = calculate_coverage(data, forest_classes, area_2d)
    urban_pct = calculate_coverage(data, urban_classes, area_2d)
    water_pct = calculate_coverage(data, water_classes, area_2d)
    
    results.append({
        'year': year,
        'agricultural': agri_pct,
        'forest': forest_pct,
        'urban': urban_pct,
        'water': water_pct
    })
    
    print(f"{year}: Agricultural={agri_pct:.1f}%, Forest={forest_pct:.1f}%, Urban={urban_pct:.1f}%, Water={water_pct:.1f}%")

df_coverage = pd.DataFrame(results)
print("\nLand cover categories:")
print(df_coverage)

# Detailed breakdown of agricultural land by crop type
agri_class_names = {
    12: 'Non-irrigated arable', 13: 'Permanently irrigated',
    14: 'Rice fields', 15: 'Vineyards', 16: 'Fruit trees',
    17: 'Olive groves', 18: 'Pastures', 19: 'Annual+permanent crops',
    20: 'Complex cultivation', 21: 'Agriculture+natural',
    22: 'Agro-forestry'
}

agri_breakdown = []
for i, year in enumerate(corine_ds.time.dt.year.values):
    data = corine_ds.isel(time=i)
    row = {'year': year}
    
    for cls in agricultural_classes:
        pct = calculate_coverage(data, [cls], area_2d)
        row[agri_class_names[cls]] = pct
    
    agri_breakdown.append(row)

df_agri = pd.DataFrame(agri_breakdown)
print("Agricultural land breakdown:")
print(df_agri)

# Create output directory
output_dir = data_dir / admin_id / "corine_land_cover"
output_dir.mkdir(parents=True, exist_ok=True)

# Save main land cover categories
csv_file_main = output_dir / "land_cover_categories.csv"
df_coverage.to_csv(csv_file_main, index=False)
print(f"Saved main land cover categories to: {csv_file_main}")

# Save agricultural land breakdown
csv_file_agri = output_dir / "agricultural_land_breakdown.csv"
df_agri.to_csv(csv_file_agri, index=False)
print(f"Saved agricultural land breakdown to: {csv_file_agri}")

# Save the most recent year as GeoTIFF for use in tutorial
latest_year_idx = -1
latest_data = corine_ds.isel(time=latest_year_idx)
latest_year = int(corine_ds.time.dt.year.values[latest_year_idx])

# Set spatial dimensions and CRS for rioxarray
latest_data = latest_data.rio.set_spatial_dims(x_dim='lon', y_dim='lat')
latest_data.rio.write_crs("EPSG:4326", inplace=True)

# Save as GeoTIFF
tif_file = output_dir / f"corine_landcover_{admin_id}.tif"
latest_data.rio.to_raster(tif_file, driver="GTiff")
print(f"Saved processed {latest_year} CORINE land cover raster to: {tif_file}")