# ESA Climate Toolbox imports for accessing and plotting ESA CCI data
from esa_climate_toolbox.core import get_op              # Get predefined operations (e.g., time series, averages)
from esa_climate_toolbox.core import list_ecv_datasets   # List available datasets per ECV (Essential Climate Variable)
from esa_climate_toolbox.core import get_store           # Connect to an ESA data store
from esa_climate_toolbox.core import list_datasets       # List all datasets in a store
from esa_climate_toolbox.core import search
from esa_climate_toolbox.ops import plot                 # High-level plotting functions for CCI data

# Used for connecting to remote data sources (e.g., ESA CCI ODP)
from xcube.core.store import new_data_store

# For displaying geographic regions interactively (if desired)
from IPython.display import GeoJSON
import shapely.geometry   # Handling geometric objects like bounding boxes

# Core data science libraries
import pandas as pd       # For tabular data handling and time manipulation
import numpy as np        # For numerical operations
import xarray as xr       # For multi-dimensional climate data structures

# Mapping and plotting
import matplotlib.pyplot as plt       # Plotting library
import cartopy.crs as ccrs            # Cartographic projections for spatial data
import IPython.display 

# Notebook settings
import warnings
warnings.filterwarnings("ignore") # Suppressing warnings to keep notebook output clean
%matplotlib inline

#list_datasets()

list_ecv_datasets("LST")

# Open the ESA CCI zarr store
cci_zarr_store = new_data_store("esa-cci-zarr")

# Which data do we have in the ESA CCI zarr store ?
list_datasets("esa-cci-zarr")

data_id ='ESACCI-LST-L3S-LST-IRCDR_-0.01deg_1MONTHLY_DAY-1995-2020-fv3.00.zarr'

# Open the ESA CCI zarr store
cci_zarr_store.describe_data(data_id)

#open the LST dataset
LST_ds = cci_zarr_store.open_data(
    data_id=data_id)
display(LST_ds)

from esa_climate_toolbox.core import list_operations
list_operations()

# Set bounding box for London
lon_min, lon_max = -0.7, 0.5
lat_min, lat_max = 51.0, 51.9

#It is located between latitudes 51°40′ and 51°1′ N and longitudes 0°30′ W and 0°20′ E;
bbox = (lon_min, lat_min, lon_max, lat_max)

##Display bounding box for reference
IPython.display.GeoJSON(shapely.geometry.box(*bbox).__geo_interface__)

# Getting the subset operator
subset_spatial_op = get_op('subset_spatial')
lst_ds_LDN = subset_spatial_op(ds=LST_ds, region=bbox)

lst_ds_LDN

# We need to call our operator first, in this case: 'tseries_mean'
# Getting the time-series mean operator
ts_mean_op = get_op('tseries_mean')

# Apply the operation to our dataset , specifying the variable of interest 
lst_LDN_mean = ts_mean_op(
    ds=lst_ds_LDN,
    var='lst'    
)

#this operator gives us a dataset with same coordinates, dimension but with a new variable: the time series of the mean of LST in London
lst_LDN_mean

lst_LDN_mean_plot = lst_LDN_mean.lst_mean.compute()  

# Let's make our first plot of the computed yearly mean variable
# Plotting the global mean LST
fig = plt.figure(figsize=(12,3))
# Daily mean LST
(lst_LDN_mean_plot-273.15).plot(c='k',linewidth=1,label='monthly')
# Yearly mean LST
# Add the legend
plt.legend()
# Add a grid 
plt.grid(True, which='both',linestyle='--')
# Add a title and axis labels
plt.title('Monthly LST averaged over London')
plt.ylabel('LST (°C)')
plt.xlabel('Year')
#how to save the figure
#plt.savefig('Mean-LST-1995-2020.png')
plt.show()

#Now let's take a look at yearly mean LST , first we use the Toolbox again to compute the yearly mean:
resample_time = get_op('temporal_aggregation') #fetching the operator to temporally aggregate our data

lst_LDN_yr_mean  =resample_time(
    ds=lst_LDN_mean,
    method='mean',
    period='1YE'
) #note: it is also possible to use this operator for different periods, as well as computing maximum and minimum for a given period.

lst_LDN_yr_mean

lst_LDN_yr_mean_plot = lst_LDN_yr_mean.lst_mean.compute()  

# Let's make our first plot of the computed yearly mean variable
# Plotting the global mean LST
fig = plt.figure(figsize=(12,3))
# Daily mean LST
(lst_LDN_yr_mean_plot-273.15).plot(c='k',linewidth=1,label='yearly')
# Yearly mean LST
# Add the legend
plt.legend()
# Add a grid 
plt.grid(True, which='both',linestyle='--')
# Add a title and axis labels
plt.title('Yearly LST averaged over London')
plt.ylabel('LST (°C)')
plt.xlabel('Year')
#how to save the figure
#plt.savefig('Mean-LST-1995-2020.png')
plt.show()

#Type your code here
london_lat, london_lon = 51.511448, -0.116414 # Kings College 
ald_lat, ald_lon =  51.8, -0.6  #Aldbury

#...

### Solution
time_slice="2019-07-01", "2019-07-01"
subset_temporal_op=get_op('subset_temporal')
lst_ds_sub_july=subset_temporal_op(ds=lst_ds_LDN,time_range=time_slice )
plot_op = get_op('plot')
plot(lst_ds_sub_july,
     var='lst', 
    indexers={'time': '2011-06-01'},  # Specify the date
    properties="cmap='viridis'") # Pass other properties here, such as color
london_lat, london_lon = 51.511448, -0.116414 # KC is 51.511448, and the longitude is -0.116414
stalb_lat, stalb_lon =  51.8, -0.6  #51.755001, and the longitude is -0.336000
plt.plot(london_lon, london_lat, linestyle='--', marker='o', color='k', label='London')
plt.plot(stalb_lon, stalb_lat, linestyle='--', marker='o', color='r', label='Aldbury')
#Aldbury ; Latitude · 51.800 ; Longitude · -0.600
plt.legend()

#Type your code here

#Solution
climatology_op = get_op('climatology')
lst_ds_clim = climatology_op(ds=lst_ds_LDN, var='lst')

#Type your code here

#Solution 1
from xcube.core.store import new_data_store
cci_zarr_store = new_data_store("esa-cci-zarr")
lc_ds = cci_zarr_store.open_data(
    'ESACCI-LC-L4-LCCS-Map-300m-P1Y-1992-2015-v2.0.7b.zarr'
)
lc_ds_sub = subset_temporal_op(ds=lc_ds, time_range=['2015-07-03', '2015-07-03'])

#Type your code here

#Solution 2
select_var_op = get_op('select_var')
lc_ds_sub = select_var_op(ds=lc_ds_sub, var="lccs_class")
normalize_op = get_op('normalize')
lc_ds_sub = normalize_op(ds=lc_ds_sub)

coregister_op = get_op('coregister')

downsampled_lc = coregister_op(
    ds_primary=lst_ds_sub_july,
    ds_replica=lc_ds_sub,
    method_ds="mode"
)
downsampled_lc

#Type your code here

# Solution
plot_op(
    ds=downsampled_lc,           # Dataset variable (should be xarray.Dataset or DataArray)
    var='lccs_class',            # Variable name
    title="Land Cover Map",
    properties="cmap='tab20'"    # Choose a colormap suitable for categorical data
)

# Colour map:
from matplotlib.colors import ListedColormap, BoundaryNorm

#Define some simplified example LCCS classes (you can expand as needed)
class_ids = [10, 50, 100, 130, 150, 160, 190, 200, 210, 220]

class_names = [
    "Cropland",
    "Forest", 
    "Shrubland",
    "Grassland & Herbaceous",
    "Sparse Vegetation",
    "Wetlands",
    "Urban",
    "Bare Areas",
    "Water",
    "Snow/Ice"
]
colors = [
    "#ffff64", "#006400", "#8ca000", "#ffb432", "#ffebaf",
    "#00785a", "#c31400", "#fff5d7", "#0046c8", "#ffffff"
]

# Create colormap and normalization
cmap = ListedColormap(colors)

#Type your code here

# Solution 
norm = BoundaryNorm(class_ids + [max(class_ids)+1], ncolors=cmap.N)
#Plot
plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
img = downsampled_lc.lccs_class[0,:,:].plot(
    ax=ax,
    cmap=cmap,
    norm=norm,
    transform=ccrs.PlateCarree(),
    add_colorbar=False
)
ax.coastlines()
london_lat, london_lon = 51.511448, -0.116414 # KC is 51.511448, and the longitude is -0.116414
ald_lat, ald_lon =  51.815 , -0.58 #, and the longitude is -0.336000
plt.plot(london_lon, london_lat, linestyle='--', marker='x', color='white', label='Kings College')
plt.plot(ald_lon, ald_lat, linestyle='--', marker='x', color='k', label='Aldbury')
#Aldbury ; Latitude · 51.800 ; Longitude · -0.600

# Add lat/lon gridlines
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.5, color="gray", alpha=0.7, linestyle="--")

ax.set_title("Land Cover Map in London")

# Add custom legend
from matplotlib.patches import Patch
legend_handles = [Patch(color=colors[i], label=class_names[i]) for i in range(len(class_ids))]
plt.legend(handles=legend_handles, bbox_to_anchor=(1.05, 1), loc="upper left")

plt.tight_layout()

#Type your code here

#Solution
#class_ids = [10, 50, 100, 130, 150, 160, 190, 200, 210, 220]
#[10, 50, 100, 130, 150, 160, 200, 220]
urban_mask = downsampled_lc['lccs_class'] == 190 # Boolean mask for urban points
urban_mask_da = urban_mask.astype(int).to_dataset(name="urban_mask")  # Convert to dataset with int values

# Apply the threshold to create the non-urban mask
non_urban_mask = ~urban_mask  # Invert the urban mask
non_urban_mask_da = non_urban_mask.astype(int).to_dataset(name="non_urban_mask")  # Convert to dataset with int values

# Apply the threshold to create another mask, f.i vegetation
#veg_mask =  downsampled_lc['lccs_class'] <= 50
#veg_mask_da =  veg_mask.astype(int).to_dataset(name="veg_mask")  # Convert to dataset with int values
# Plot the urban mask
plot_op(
    ds=urban_mask_da,          # Dataset variable
    var='urban_mask',          # Variable name
    title="Urban Mask",
    properties="cmap='Greys'"  # Use a greyscale colormap
)


# Plot the non-urban mask
plot_op(
    ds=non_urban_mask_da,      # Dataset variable
    var='non_urban_mask',      # Variable name
    title="Non-Urban Mask",
    properties="cmap='Greys'"  # Use a greyscale colormap
)


# Plot the vegetation max
#plot_op( ds=veg_mask_da,  var='veg_mask',     title="Veg Mask",  properties="cmap='Greys'"  )

#Type your code here

#Solution
lst_data_array = lst_ds_clim['lst']  # Replace 'lst' with the correct variable name

# Flip the urban_mask along the latitude axis
urban_mask_flipped = np.flipud(urban_mask)  # Flip along the vertical axis

# Apply the flipped mask to the LST data
urban_lst = xr.where(urban_mask_flipped, lst_data_array, np.nan)
non_urban_lst = xr.where(~urban_mask_flipped, lst_data_array, np.nan)

# Plot an example time step
time_step = lst_data_array.time[5]  # Choose the sixth time step for June
plt.figure(figsize=(12, 6))

# Plot Urban LST
plt.subplot(1, 2, 1)
(urban_lst-273.15).sel(time=time_step).plot(cmap="plasma_r",vmin=25, vmax = 35)
plt.title(f"Urban LST at {time_step.values}")
plt.plot(london_lon, london_lat, linestyle='--', marker='x', color='red', label='Kings College')
plt.plot(stalb_lon, stalb_lat, linestyle='--', marker='x', color='k', label='Aldbury')
# Plot Non-Urban LST
plt.subplot(1, 2, 2)
(non_urban_lst-273.15).sel(time=time_step).plot(cmap="plasma_r",vmin=25, vmax = 35)
plt.plot(london_lon, london_lat, linestyle='--', marker='x', color='red', label='Kings College')
plt.plot(stalb_lon, stalb_lat, linestyle='--', marker='x', color='k', label='Aldbury')
plt.title(f"Non-Urban LST at {time_step.values}")

plt.tight_layout()
plt.show()

#Type your code here

#Solution
#Compute the mean LST for urban and non-urban pixels over time
urban_mean = urban_lst.mean(dim=["lat", "lon"])
non_urban_mean = non_urban_lst.mean(dim=["lat", "lon"])
UHI = urban_mean - non_urban_mean 
# Plot the time series
plt.figure(figsize=(10, 6))
(urban_mean- 273.15).plot(label="Urban")
UHI.plot(label="UHI")
plt.title("UHI over London")
plt.ylabel("UHI (K)")
plt.xlabel("Time")
plt.legend()
plt.show()