corentin kuster

import requests
import pandas as pd
import io

url = "https://ec.europa.eu/eurostat/api/dissemination/sdmx/2.1/data/T2020_RD330/?format=CSV"
urlData = requests.get(url).content
renew_df = pd.read_csv(io.StringIO(urlData.decode('utf-8')))
renew_df.iloc[:,0] = renew_df.iloc[:,0].apply(lambda x: x.split(';')[-1])
renew_df.rename(columns={'freq;nrg_bal;unit;geoTIME_PERIOD': 'CODE'},inplace=True)
renew_df

import geopandas as gpd

all_data = pd.concat([renew_df,energy_df,gas_df,elc_df,typo_df,stock_dfeu_df[['geometry']]],join='inner',axis=1).replace(': ',np.nan) #eurostat set nan values as ':'
all_data_gdf  = gpd.GeoDataFrame(all_data)

import xarray as xr
import rioxarray
import rasterio

with xr.open_dataset('tg_0.50deg_reg_1995-2016_v14.0.nc') as file:
    dsCDF = file.rio.write_crs("epsg:4326", inplace=True)
    
dsCDF

#Create typical year from the 20 years'set with daily data point (grid X 366 days)
typical_year = dsCDF.groupby('time.dayofyear').median('time') 

#Split typical year into seasonal bins (grid X 2 seasons)
day_bins = [0,106,289,366] #Limits of the 2 seasons (day number over the year)
typical_Seasonal_year = typical_year.groupby_bins('dayofyear',day_bins,labels=['cold1','hot','cold2']).median('dayofyear') #Split the year into 3 bins
Cold = typical_Seasonal_year.sel(dayofyear_bins = ['cold2','cold1']).mean('dayofyear_bins') #Average of the 2 "Cold" season bins
Hot = typical_Seasonal_year.sel(dayofyear_bins = ['hot'])

#Note that this could have been done before defining Hot and Cold ¯_ツ_/¯
Hot.rio.write_crs("epsg:4326", inplace=True) #set crs to our raster
Hot = Hot.rio.reproject(dst_crs='EPSG:3857') #reproject because geojson in EPSG:3857
Hot = Hot.rio.clip(eu)
Hot = Hot.where(Hot['tg'] != -9999.)  #conserve only necesaary values

Cold.rio.write_crs("epsg:4326", inplace=True)
... #idem

hot_bins = Hot.quantile([0,0.25,0.5,0.75,1]).tg.values
Hot_discr = np.digitize(Hot.tg.values, hot_bins)
Hot_discr = Hot_discr.astype(np.int32)

#idem for Cold raster

from rasterio.features import shapes

west, east = min(Hot.x.values), max(Hot.x.values) #get corners coordinates
south, north =  min(Hot.y.values), max(Hot.y.values) #get corners coordinates
width = len(Hot.x.values)
height = len(Hot.y.values)

res_x = np.diff(Hot.x.values).mean() #get the x gradient
res_y = np.diff(Hot.y.values).mean() #get the y gradient

transform = rasterio.transform.from_origin(west, north, res_x, -res_y) #Return an Affine transformation for a georeferenced raster 

data = []
mask = Cold_discr != 5
for shp, val in shapes(Cold_discr,mask=mask, transform=transform):
    record = {'properties':{'bin': val}, 'geometry' : shp}
    data.append(record)

#let's make a geodataframe from it
Hot_df = gpd.GeoDataFrame.from_features(data,crs=3857) 
Hot_df = Hot_df.dissolve(by='bin')

for index, row in Hot_df.iterrows(): #for each bin
    temp_poly = row.geometry.buffer(100*1000) #buffer the current geometry
    temp_poly = temp_poly.intersection(eu)    #intersect it with EU geometry (outline)
    other_polys = Hot_df.loc[Hot_df.index != index].geometry.values.unary_union() #union of all the other polygon
    to_remove = temp_poly.intersection(other_polys) #intersection to retrieve what to remove (because of the previous buffer)
    temp_poly = temp_poly.difference(to_remove) #removal of teh extra bits

    Hot_df.loc[index,['geometry']]= temp_poly #overwrite geometry in dataframe

all_data_dis = all_data.apply(lambda col: pd.qcut(col, 4, labels=False, duplicates = 'drop')) #qcut is a Quantile-based discretization function.
all_data_dis = all_data_dis+1 #we add 1 so the value range from 1 to 4 (and not 0 to 3)

final_data.to_postgis(name, con, schema=None, if_exists='fail', index=False, index_label=None, chunksize=None, dtype=None)

final_data.to_file("data.shp")

final_data.to_file("data.geojson", driver='GeoJSON')
...