Agricultural land from satellite data

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

# Year range for data download
start_year = 1992
end_year = 2022  # Adjust based on latest available data

import cdsapi
import zipfile
import xarray as xr
import numpy as np
import pandas as pd
import geopandas as gpd
import regionmask
from pathlib import Path
import os

# Set up data directories
data_dir = Path("../data")
satellite_dir = data_dir / "satellite_landcover"
satellite_dir.mkdir(exist_ok=True)

print(f"\nSatellite data directory: {satellite_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')

# Download land cover data for each year
for yyyy in range(start_year, end_year + 1):
    print(f"Downloading land cover data for {yyyy}")

    # Check if file already exists
    zip_file = satellite_dir / f'land_cover_{admin_id}_{yyyy}.zip'
    if zip_file.exists():
        # print(f"  - File {zip_file.name} already exists. Skipping download.")
        continue
    
    # Select version based on year
    if yyyy <= 2015:
        lc_version = "v2_0_7cds"
    else:
        lc_version = "v2_1_1"
    
    dataset = "satellite-land-cover"
    request = {
        "variable": "all",
        "year": [str(yyyy)],
        "version": [lc_version],
        "area": [lat_max, lon_min, lat_min, lon_max]  # [North, West, South, East]
    }

    client = cdsapi.Client()
    client.retrieve(dataset, request).download(str(zip_file))
    print(f"  - Downloaded {zip_file.name}")

print("\nLand cover download complete!")

# Inspect one zip file to create mask
first_year = start_year
zip_file = satellite_dir / f'land_cover_{admin_id}_{first_year}.zip'

with zipfile.ZipFile(zip_file, 'r') as zip_ref:
    file_list = zip_ref.namelist()
    nc_files_in_zip = [f for f in file_list if f.endswith('.nc')]
    
    if nc_files_in_zip:
        lc_filename = nc_files_in_zip[0]
        zip_ref.extract(lc_filename, satellite_dir)
    else:
        print("No NetCDF file found in zip!")

# Load land cover data to get grid
lc_file = satellite_dir / lc_filename
lc_ds = xr.open_dataset(lc_file)

# Subset to our region (initial bounding box)
lc_region = lc_ds.sel(
    lon=slice(lon_min, lon_max),
    lat=slice(lat_max, lat_min)
)

# Create admin region mask: 1 where inside region, 0 elsewhere
lc_admin_mask_raw = admin_mask.mask(lc_region.lon, lc_region.lat)
lc_admin_mask = xr.DataArray(
    (~np.isnan(lc_admin_mask_raw.values)).astype(float),
    coords={'lat': lc_region.lat, 'lon': lc_region.lon},
    dims=['lat', 'lon']
)

print(f"Admin region area (grid cells): {lc_admin_mask.sum().values:.0f}")
print(f"Grid shape: {lc_admin_mask.shape}")

# Clean up extracted file
lc_file.unlink()

# Initialize list to store results
ag_data = []
latest_ag_mask = None  # Store latest year for saving as raster
latest_lc_region = None  # Store latest region data for coordinate info

for yyyy in range(start_year, end_year + 1):
    # Extract and load NetCDF file
    zip_file = satellite_dir / f'land_cover_{admin_id}_{yyyy}.zip'
    
    with zipfile.ZipFile(zip_file, 'r') as zip_ref:
        file_list = zip_ref.namelist()
        nc_files_in_zip = [f for f in file_list if f.endswith('.nc')]
        
        if nc_files_in_zip:
            lc_filename = nc_files_in_zip[0]
            zip_ref.extract(lc_filename, satellite_dir)
        else:
            print(f"  - No NetCDF file found in {zip_file.name}. Skipping.")
            continue

    # Load land cover data
    lc_file = satellite_dir / lc_filename
    lc_ds = xr.open_dataset(lc_file)

    # Subset to our region
    lc_region = lc_ds.sel(
        lon=slice(lon_min, lon_max),
        lat=slice(lat_max, lat_min)
    )

    # Create agricultural land mask (LCCS classes 10-40)
    ag_mask = ((lc_region['lccs_class'] >= 10) & (lc_region['lccs_class'] <= 40)).astype(float)

    # Calculate agricultural land fraction within admin region
    ag_fraction = (ag_mask * lc_admin_mask).sum().values / lc_admin_mask.sum().values
    
    # Store result
    ag_data.append({'year': yyyy, 'ag_fraction': ag_fraction})
    
    # Save latest year for raster output
    if yyyy == end_year:
        latest_ag_mask = ag_mask.copy()
        latest_lc_region = lc_region.copy()
    
    # Clean up extracted file
    lc_file.unlink()

# Create pandas DataFrame
ag_timeseries = pd.DataFrame(ag_data)
ag_timeseries

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

# Save to CSV
csv_file = output_dir / f'agricultural_land_fraction_{admin_id}.csv'
ag_timeseries.to_csv(csv_file, index=False)
print(f"Saved agricultural land fraction timeseries to: {csv_file}")
print(f"\nSummary statistics:")
print(ag_timeseries['ag_fraction'].describe())

# Save the most recent year as GeoTIFF for use in tutorial
if latest_ag_mask is not None and latest_lc_region is not None:
    # Set spatial dimensions and add CRS information
    latest_ag_mask = latest_ag_mask.rio.set_spatial_dims(x_dim='lon', y_dim='lat')
    latest_ag_mask.rio.write_crs("EPSG:4326", inplace=True)
    
    # Save as GeoTIFF
    tif_file = output_dir / f"satellite_landcover_{admin_id}.tif"
    latest_ag_mask.rio.to_raster(tif_file, driver="GTiff")
    print(f"Saved {end_year} satellite land cover raster to: {tif_file}")
else:
    print("No latest year data available to save as raster.")