2. 1. extract and process Eduplanet climate data test
3. 1. python3
4. 1. 07.07.2024 0000.0003

import xarray as xr

import matplotlib.pyplot as plt import numpy as np import math

import skimage

import matplotlib as mpl from matplotlib.colors import LinearSegmentedColormap, ListedColormap

import scipy from scipy.ndimage import distance_transform_edt from scipy.ndimage import gaussian_filter

from scipy.interpolate import RectBivariateSpline from scipy.interpolate import griddata from scipy.interpolate import interpn

from scipy.ndimage import generic_filter

from sklearn.linear_model import LinearRegression from sklearn import svm from sklearn.ensemble import RandomForestRegressor from sklearn.preprocessing import LabelEncoder from sklearn.preprocessing import StandardScaler from sklearn.svm import SVR

from pygam import LinearGAM, s from pygam import GammaGAM, InvGaussGAM, LogisticGAM, PoissonGAM , ExpectileGAM

def downscale_wind(dem0, windu0, windv0, windw0):

   shapeb=(180,180)
   shapeb=(360,360)
   shape2=(180,360)
   dem=interpola_tosize(dem0, shapeb)
   windu=interpola_tosize(windu0, shapeb)
   windv=interpola_tosize(windv0, shapeb)
   windw=interpola_tosize(windw0, shapeb)
   #plt.imshow(windw0)
   #plt.show()
   #quit(-1)

   wind_strength=1 ## near reality
   dem_shape=np.shape(dem)
   dem_size=dem_shape[0]
   # Lasketaan gradientti (x- ja y-suunnassa)
   grad_y, grad_x = np.gradient(dem)
   # Tuulen suunta ja voimakkuus
   #wind_direction = np.array([1, 0, 0])  # länsituuli, aluksi ilman z-komponenttia
   
   #wind_strength = 1.0
   # Tuulivektorit (16x16)
   wind_size = 360
   wind_field = np.zeros((wind_size, wind_size, 3))
   #wind_field[:, :, 0] = wind_direction[0] * wind_strength
   #wind_field[:, :, 1] = wind_direction[1] * wind_strength
   wind_field[:, :, 0] = windu
   wind_field[:, :, 1] = windv
   # Skaalataan tuulivektorit
   influence = 0.1 * 300  # Poikkeutusvaikutuksen kerroin
   scale_factor = dem.shape[0] // wind_field.shape[0]
   for i in range(wind_field.shape[0]):
       for j in range(wind_field.shape[1]):
           dem_i = i * scale_factor
           dem_j = j * scale_factor
           wind_field[i, j] = adjust_wind_vector(wind_field[i, j], grad_x[dem_i, dem_j], grad_y[dem_i, dem_j], wind_strength, influence)
   # Visualisoidaan DEM ja tuulikenttä
   #plt.figure(figsize=(12, 12))
   # Visualisoidaan DEM
   #plt.subplot(2, 2, 1)
   #plt.imshow(dem, cmap='terrain')
   #plt.colorbar(label='Elevation')
   #plt.title('Simple DEM with Gaussian Peaks')
   # Visualisoidaan tuulikenttä DEM:n päällä
   #plt.subplot(2, 2, 2)
   #plt.imshow(dem, cmap='terrain', alpha=0.5)
   #for i in range(wind_field.shape[0]):
   #    for j in range(wind_field.shape[1]):
   #        plt.arrow(j * scale_factor, i * scale_factor, wind_field[i, j, 0], wind_field[i, j, 1],
   #              head_width=1, head_length=1, fc='blue', ec='blue')
   #plt.title('Wind Field Adjusted by DEM Gradient')
   # Visualisoidaan tuulivirtausviivat
   #plt.subplot(2, 2, 3)
   x = np.linspace(0, dem_size - 1, wind_size)
   y = np.linspace(0, dem_size - 1, wind_size)
   X, Y = np.meshgrid(x, y)
   U = wind_field[:, :, 0]
   V = wind_field[:, :, 1]
   #plt.imshow(dem, cmap='terrain', alpha=0.5, extent=(0, dem_size, 0, dem_size))
   #plt.streamplot(X, Y, U, V, color='blue')
   #plt.title('Wind Streamlines Adjusted by DEM Gradient')
   #Visualisoidaan tuulen z-komponentti
   #plt.subplot(2, 2, 4)
   Z = np.zeros((dem_size, dem_size))
   for i in range(wind_field.shape[0]):
       for j in range(wind_field.shape[1]):
           Z[i*scale_factor:(i+1)*scale_factor, j*scale_factor:(j+1)*scale_factor] = wind_field[i, j, 2]

   #plt.imshow(Z, cmap='viridis', extent=(0, dem_size, 0, dem_size))
   #plt.imshow(U, cmap='viridis', extent=(0, dem_size, 0, dem_size))    
   #plt.streamplot(X, Y, U, V, color='blue')
   #plt.colorbar(label='Wind Z-Component')
   #plt.title('Wind Z-Component')
   #plt.tight_layout()
   #plt.show()
   uusu0=U
   uusv0=V
   uusz0=Z
   uusu=interpolate_raster(uusu0, shape2, "cubic")
   uusv=interpolate_raster(uusv0, shape2, "cubic")
   uusz=interpolate_raster(uusz0, shape2, "cubic")
   #print (np.shape(U))
   #plt.imshow(uusz)
   #plt.show()

   return(uusu, uusv, uusz)

def downscale_machinelearn_log(temp, inputy, slons, slats, blons, blats):

   ## NOK
   #temp_min=np.min(temp)
   #temp_max=np.max(temp)
   #height_min=np.min(dem_data)
   #height_max=np.max(dem_data)
   temp=np.log10(temp)
   temp=np.nan_to_num(temp,0.01, neginf=0.01)

   avar1=inputy[0]
   shape2=np.shape(avar1) 
   shape1=np.shape(temp)

   #temptar=interpola_tosize(temp, shape2)
   temptar0=interpolate_raster(temp, (40,40), "linear")
   temptar=interpolate_raster(temptar0, shape2, "cubic")

   #np.inf_to_num(temp,4)
   avar1=inputy[0]
   bvar1=inputy[1]
   cvar1=inputy[2]
   dvar1=inputy[3]
   evar1=inputy[4]  
   fvar1=temptar 
   #gvar1=inputy[5]   

   outnx=shape2[1]
   outny=shape2[0]

   #plt.imshow(y1)
   #temptar=get_interpolated_raster(bvar1, blons, blats, temp, slons, slats)

   #plt.imshow(temptar)
   #plt.show()
   #quit(-1)
   #inputy=[  dem, latrasb, lonrasb, relative_wind_angle_deg,  distance_to_sea  ]
   avar0=get_sized_raster(temp, slons, slats, avar1, blons, blats)
   bvar0=get_sized_raster(temp, slons, slats, bvar1, blons, blats)
   cvar0=get_sized_raster(temp, slons, slats, cvar1, blons, blats)
   dvar0=get_sized_raster(temp, slons, slats, dvar1, blons, blats)
   evar0=get_sized_raster(temp, slons, slats, evar1, blons, blats)
   fvar0=get_sized_raster(temp, slons, slats, fvar1, blons, blats)
   #gvar0=get_sized_raster(temp, slons, slats, gvar1, blons, blats)
   #temptar=skimage.transform.resize(temp, np.shape(inputy[0]) ) 
   #bvar1=inputy[0]
   #bvar0=skimage.transform.resize(inputy[0], (17,17))
   #avar1=inputy[1]
   #avar0=skimage.transform.resize(inputy[1], (17,17))
   #cvar1=inputy[2]
   #cvar0=skimage.transform.resize(inputy[2], (17,17))
   #dvar1=inputy[3]
   #dvar0=skimage.transform.resize(inputy[3], (17,17))
   #ivar1=np.dstack([ avar1])
   #ivar0=np.dstack([ avar0])
   #ivar1=np.dstack([bvar1, avar1])
   #ivar0=np.dstack([ bvar0, avar0])
   #ivar1=np.dstack([dvar1, cvar1, bvar1, avar1])
   #ivar0=np.dstack([dvar0, cvar0, bvar0, avar0])
   #ivar1=np.dstack([avar1, dvar1, evar1])
   #ivar0=np.dstack([avar0, dvar0, evar0])  
   ivar1=np.dstack([avar1, dvar1, bvar1, evar1])
   ivar0=np.dstack([avar0, dvar0, bvar0, evar0])  
   #X2 = bvar1.flatten().reshape(-1, 1)
   #X = bvar0.flatten().reshape(-1, 1)
   #X2 = ivar1.flatten().reshape(-1, 1)
   #X = ivar0.flatten().reshape(-1, 1)
   X2 = ivar1.flatten().reshape(-1, 4)
   X = ivar0.flatten().reshape(-1, 4)

   #y = temp_interpolated.flatten()
   y = temp.flatten()
   scaler = StandardScaler().fit(X)
   X_train_scaled = scaler.transform(X)
   scaler2 = StandardScaler().fit(X2)
   X2_test_scaled = scaler.transform(X2)

   #model=SVR(kernel = 'linear') ##bestis, but slow
   #model=SVR(kernel = 'poly') # nok
   #model=SVR(kernel = 'rbf') #3nok
   #y_test_pred = svr_rbf.predict(X2_test_scaled)
   #y_test_pred = svr_poly.predict(X2)
   #Y2=y_test_pred
   
   #model = LinearRegression()
 
   #model = RandomForestRegressor(n_estimators=300, random_state=0, oob_score=True)
   #model = LinearGAM(s(0))
   #model = LinearGAM() #3 ok

   # from pygam import GammaGAM, InvGaussGAM, LogisticGAM, PoissonGAM , ExpectileGAM
   #model=GammaGAM()
   #model=InvGaussGAM()
   #model=LogisticGAM()
   #model= PoissonGAM () #3 nok
   model= ExpectileGAM () ## OK
 
   model=model.fit(X, y)
   Y2 = model.predict(X2)
   tempds = Y2.reshape(outny, outnx)
   
   #tout=tempds
   tout=np.power(10,tempds)
   #tout=temptar
   #plt.imshow(tout)
   #plt.show()
   #quit(-1)
   #tout=(temptar-tempds)/2
   return(tout )

def calculate_shadow_map_smooth_base(dem, wind_field_x, wind_field_y, light_angle_degrees, angle_range_degrees, shadow_intensity_decrease, shadow_intensity_range):

   height, width = dem.shape
   shadow_map = np.ones((height, width))  # Initialize with maximum intensity
   wind_field_x = -wind_field_x
   wind_field_y = -wind_field_y
   light_angle_radians = np.radians(light_angle_degrees)
   light_slope = np.tan(light_angle_radians)

   for y in range(height):
       for x in range(width):
           cx, cy = x, y
           current_height = dem[cy, cx]
           shadow_intensity = shadow_intensity_range[1]

           while 0 <= cx < width and 0 <= cy < height:
               wind_x, wind_y = add_wind_randomness(wind_field_x[int(cy), int(cx)], wind_field_y[int(cy), int(cx)], angle_range_degrees)
               cx += wind_x
               cy += wind_y

               if 0 <= cx < width and 0 <= cy < height:
                   distance = np.sqrt(wind_x**2 + wind_y**2)
                   height_along_light = current_height + light_slope * distance

                   if dem[int(cy), int(cx)] > height_along_light:
                       shadow_intensity *= shadow_intensity_decrease  # Decrease intensity in the shadow
                       shadow_intensity = max(shadow_intensity_range[0], shadow_intensity)
                       shadow_map[int(cy), int(cx)] = shadow_intensity
                   else:
                       shadow_intensity *= shadow_intensity_decrease  # Gradually decrease intensity
                       shadow_intensity = max(shadow_intensity_range[0], shadow_intensity)
                       shadow_map[int(cy), int(cx)] = shadow_intensity

   return shadow_map

def calculate_shadow_map_smooth(dem, wind_field_x, wind_field_y, light_angle_degrees, angle_range_degrees, shadow_intensity_decrease, shadow_intensity_range):

   height, width = dem.shape
   shadow_map = np.ones((height, width))  # Initialize with maximum intensity
   wind_field_x = -wind_field_x
   wind_field_y = -wind_field_y
   light_angle_radians = np.radians(light_angle_degrees)
   light_slope = np.tan(light_angle_radians)

   for y in range(height):
       for x in range(width):
           cx, cy = x, y
           current_height = dem[cy, cx]
           shadow_intensity = shadow_intensity_range[1]

           while 0 <= cx < width and 0 <= cy < height:
               wind_x, wind_y = add_wind_randomness(wind_field_x[int(cy), int(cx)], wind_field_y[int(cy), int(cx)], angle_range_degrees)
               cx += wind_x
               cy += wind_y

               if 0 <= cx < width and 0 <= cy < height:
                   distance = np.sqrt(wind_x**2 + wind_y**2)
                   height_along_light = current_height + light_slope * distance

                   if dem[int(cy), int(cx)] > height_along_light:
                       shadow_intensity *= shadow_intensity_decrease  # Decrease intensity in the shadow
                       shadow_intensity = max(shadow_intensity_range[0], shadow_intensity)
                       shadow_map[int(cy), int(cx)] = shadow_intensity
                   else:
                       shadow_intensity *= shadow_intensity_decrease  # Gradually decrease intensity
                       shadow_intensity = max(shadow_intensity_range[0], shadow_intensity)
                       shadow_map[int(cy), int(cx)] = shadow_intensity

   return shadow_map

1. Calculate the shadow map

def calculate_shadow_map(dem, wind_field_x, wind_field_y, light_angle_degrees, angle_range_degrees):

   height, width = dem.shape
   shadow_map = np.zeros((height, width))

   light_angle_radians = np.radians(light_angle_degrees)
   light_slope = np.tan(light_angle_radians)

   for y in range(height):
       for x in range(width):
           cx, cy = x, y
           current_height = dem[cy, cx]
           shadow = False

           while 0 <= cx < width and 0 <= cy < height:
               wind_x, wind_y = add_wind_randomness(wind_field_x[int(cy), int(cx)], wind_field_y[int(cy), int(cx)], angle_range_degrees)
               cx += wind_x
               cy += wind_y

               if 0 <= cx < width and 0 <= cy < height:
                   distance = np.sqrt(wind_x**2 + wind_y**2)
                   height_along_light = current_height + light_slope * distance

                   if dem[int(cy), int(cx)] > height_along_light:
                       shadow = True
                       break

           if shadow:
               shadow_map[y, x] = 1

   return shadow_map

1. Calculate TPI

def tpi(dem):

   def local_mean(data):
       center = data[len(data) // 2]
       return center - np.mean(data)
   
   tpi_map = generic_filter(dem, local_mean, size=3)
   return tpi_map

1. Calculate TRI

def tri(dem):

   def local_roughness(data):
       center = data[len(data) // 2]
       return np.max(np.abs(data - center))
   
   tri_map = generic_filter(dem, local_roughness, size=3)
   return tri_map

1. Function to add randomness to the wind direction

def add_wind_randomness(wind_x, wind_y, angle_range_degrees):

   angle_range_radians = np.radians(angle_range_degrees)
   random_angle = np.random.uniform(-angle_range_radians, angle_range_radians)
   cos_angle = np.cos(random_angle)
   sin_angle = np.sin(random_angle)
   new_wind_x = wind_x * cos_angle - wind_y * sin_angle
   new_wind_y = wind_x * sin_angle + wind_y * cos_angle
   return new_wind_x, new_wind_y

1. Calculate the distance to the sea along the reversed wind direction with randomness

def distance_against_wind_with_randomness(dem, wind_field_x, wind_field_y, sea_mask, angle_range_degrees):

   height, width = dem.shape
   distance_map = np.full((height, width), np.inf)
   wind_field_x = -wind_field_x
   wind_field_y = -wind_field_y
   for y in range(height):
       for x in range(width):
           if sea_mask[y, x]:
               distance_map[y, x] = 0
           else:
               distance = 0
               cx, cy = x, y
               while 0 <= cx < width and 0 <= cy < height:
                   if sea_mask[int(cy), int(cx)]:
                       distance_map[y, x] = distance
                       break
                   wind_x, wind_y = add_wind_randomness(wind_field_x[int(cy), int(cx)], wind_field_y[int(cy), int(cx)], angle_range_degrees)
                   cx += wind_x
                   cy += wind_y
                   distance += 1
                   # Ensure cx and cy stay within bounds
                   if not (0 <= cx < width and 0 <= cy < height):
                       break

   return distance_map

1. Calculate the distance to the sea along the reversed wind direction

def distance_against_wind(dem, wind_field_x, wind_field_y, sea_mask):

   height, width = dem.shape
   distance_map = np.full((height, width), np.inf)
   wind_field_x = -wind_field_x
   wind_field_y = -wind_field_y

   for y in range(height):
       for x in range(width):
           if sea_mask[y, x]:
               distance_map[y, x] = 0
           else:
               distance = 0
               cx, cy = x, y
               while 0 <= cx < width and 0 <= cy < height:
                   if sea_mask[int(cy), int(cx)]:
                       distance_map[y, x] = distance
                       break
                   cx += wind_field_x[int(cy), int(cx)]
                   cy += wind_field_y[int(cy), int(cx)]
                   distance += 1
                   # Ensure cx and cy stay within bounds
                   if not (0 <= cx < width and 0 <= cy < height):
                       break

   return distance_map

1. Funktio tuulen suuntavektorin muuttamiseen gradientin mukaan ja z-komponentin laskemiseen

def adjust_wind_vector(wind, grad_x, grad_y, wind_strength, influence=0.1):

   # Lasketaan maaston gradientin kulma
   terrain_angle = np.arctan2(grad_y, grad_x)

   # Lasketaan tuulen kulma
   wind_angle = np.arctan2(wind[1], wind[0])

   # Erotus tuulen ja maaston kulman välillä
   angle_diff = np.mod(terrain_angle - wind_angle + np.pi, 2 * np.pi) - np.pi

   # Päätetään tuulen uusi suunta
   perpendicular = np.zeros(3)
   if angle_diff > 0:
       # Käännetään oikealle
       perpendicular[:2] = [-grad_y, grad_x]
   else:
       # Käännetään vasemmalle
       perpendicular[:2] = [grad_y, -grad_x]

   adjusted_wind = wind + influence * perpendicular

   # Lasketaan z-komponentti DEM:n gradientin avulla
   # Z-komponentti nousee länsirinteellä ja laskee itärinteellä
   if grad_x > 0:
       adjusted_wind[2] = influence * np.sqrt(grad_x**2 + grad_y**2)
   else:
       adjusted_wind[2] = -influence * np.sqrt(grad_x**2 + grad_y**2)

   # Normalisoidaan tulos
   norm = np.linalg.norm(adjusted_wind[:2])  # normalisoidaan vain x- ja y-komponentit
   if norm != 0:
       adjusted_wind[:2] = adjusted_wind[:2] / norm * wind_strength

   return adjusted_wind

def interpolate_raster(inras, shape2, method="cubic"):

   ## seems nok
   # Create meshgrid for sabras
   #method="cubic"
   shape1=np.shape(inras)
   #print("shape1 ", shape1)
   #quit(-1)

   tanx=shape2[1]

   soy=shape1[0] 
   sox=shape1[1]    

   #sablon_grid, sablat_grid = np.meshgrid(, sablats)
   grid_x, grid_y = np.mgrid[0:1:360j, 0:1:360j]
   #grid_x, grid_y = np.mgrid[0:1:tanx, 0:1:tany]

   inx=np.linspace(0,1,sox)
   iny=np.linspace(0,1,soy)

   xs1=np.repeat(inx,soy).reshape(sox,soy)
   xs2=np.ravel(xs1)
   ys1=np.repeat(iny,sox)
   ys1=np.reshape(ys1, (soy,sox))
   ys1=np.rot90(ys1)

   ys2=np.ravel(ys1)

   points=np.dstack([xs2,ys2])
   points=points[0, :, :]
   
   boints2=(inx,iny)
   values2=inras

   values=inras.ravel()
   #print("kkk")
   #print(np.shape(values))
   #print(np.shape(points))
   
   #quit(-1)
   #values=np.random(17*17)

   # Flatten the grids for interpolation
   #points = np.array([(lon, lat) for lon in biglons for lat in biglats])

   #values = bigras.flatten()

   # Perform the interpolation
   #output_raster = griddata(points, values, (grid_x, grid_y), method=method)
   output_r0 = griddata(points, values, (grid_y, grid_x), method=method)
   #output_r0 = interpn(boints2, values2, (361,361), method=method)

   #print(np.max(grid_y))
   output_r1=np.rot90(output_r0)
   output_r2=skimage.transform.resize(output_r1, shape2, anti_aliasing=True)
   #output_raster=np.reshape(output_r2, shape2)
   output_raster=np.flipud(output_r2)    
   
   return (output_raster)

def downscale_machinelearn_base(temp, inputy):

   #from pygam import LinearGAM, s
   #temp_min=np.min(temp)
   #temp_max=np.max(temp)
   #height_min=np.min(dem_data)
   #height_max=np.max(dem_data)

   temptar=skimage.transform.resize(temp, np.shape(inputy[0]) ) 
  
   bvar1=inputy[0]
   bvar0=skimage.transform.resize(inputy[0], (17,17))
   avar1=inputy[1]
   avar0=skimage.transform.resize(inputy[1], (17,17))
   cvar1=inputy[2]
   cvar0=skimage.transform.resize(inputy[2], (17,17))
   dvar1=inputy[3]
   dvar0=skimage.transform.resize(inputy[3], (17,17))

   ivar1=np.dstack([bvar1, avar1])
   ivar0=np.dstack([ bvar0, avar0])
   #ivar1=np.dstack([dvar1, cvar1, bvar1, avar1])
   #ivar0=np.dstack([dvar0, cvar0, bvar0, avar0])
   ivar1=np.dstack([dvar1, bvar1, avar1])
   ivar0=np.dstack([dvar0, bvar0, avar0])   
   #X2 = bvar1.flatten().reshape(-1, 1)
   #X = bvar0.flatten().reshape(-1, 1)
   #X2 = ivar1.flatten().reshape(-1, 1)
   #X = ivar0.flatten().reshape(-1, 1)
   X2 = ivar1.flatten().reshape(-1, 3)
   X = ivar0.flatten().reshape(-1, 3)

   #y = temp_interpolated.flatten()
   y = temp.flatten()
   scaler = StandardScaler().fit(X)
   X_train_scaled = scaler.transform(X)
   scaler2 = StandardScaler().fit(X2)
   X2_test_scaled = scaler.transform(X2)
   svr_lin = SVR(kernel = 'linear')
   svr_rbf = SVR(kernel = 'rbf')
   svr_poly = SVR(kernel = 'poly')
   #svr_lin.fit(X_train_scaled, y)
   #svr_rbf.fit(X_train_scaled, y)
   #svr_poly.fit(X_train_scaled, y)
   #svr_poly.fit(X, y)
   #model=SVR(kernel = 'linear') ##bestis
   #model=SVR(kernel = 'poly')
   #model=SVR(kernel = 'rbf')
   #y_test_pred = svr_rbf.predict(X2_test_scaled)
   #y_test_pred = svr_poly.predict(X2)
   #Y2=y_test_pred
   
   #model = LinearRegression()
 
   #model = RandomForestRegressor(n_estimators=30, random_state=0, oob_score=True)
   #model = LinearGAM(s(0))
   #model = LinearGAM()
   # from pygam import GammaGAM, InvGaussGAM, LogisticGAM, PoissonGAM , ExpectileGAM
   #model=GammaGAM()
   #model=InvGaussGAM()
   #model=LogisticGAM()   

   model=model.fit(X, y)
   Y2 = model.predict(X2)
   tempds = Y2.reshape(180, 360)

   tout=tempds
   #tout=temptar

   return(tout )

def find_nearest(array, value):

   """Find the nearest value in an array."""
   idx = (np.abs(array - value)).argmin()
   return idx

def get_interpolated_raster0(sabras, sablons, sablats, bigras, biglons, biglats, method='cubic'):

   ## seems nok
   # Create meshgrid for sabras
   method="nearest"

   sablons=sablons+180
   sablats=sablats+90 
   biglons=biglons+180
   biglats=biglats+90 
   sablon_grid, sablat_grid = np.meshgrid(sablons, sablats)
   
   # Flatten the grids for interpolation
   points = np.array([(lon, lat) for lon in biglons for lat in biglats])
   values = bigras.flatten()
   grid_x, grid_y = sablon_grid, sablat_grid

   # Perform the interpolation
   #output_raster = griddata(points, values, (grid_x, grid_y), method=method)
   output_raster = griddata(points, values, (grid_y, grid_x), method=method)   
   return output_raster

def get_sized_raster(sabras, sablons, sablats, bigras, biglons, biglats):

   """
   Generate a resized raster sabras from bigras by mapping sablons and sablats
   values to the nearest biglons and biglats values in bigras.

   Parameters:
   sabras (np.ndarray): Template raster array to fill.
   sablons (np.ndarray): Longitudes corresponding to sabras.
   sablats (np.ndarray): Latitudes corresponding to sabras.
   bigras (np.ndarray): Source raster array.
   biglons (np.ndarray): Longitudes corresponding to bigras.
   biglats (np.ndarray): Latitudes corresponding to bigras.

   Returns:
   np.ndarray: Resized raster with values from bigras mapped according to sablons and sablats.
   """
   # Get the shape of the template raster
   sab_shape = sabras.shape
   
   # Initialize the output raster
   output_raster = np.zeros(sab_shape)
   
   for i in range(sab_shape[0]):
       for j in range(sab_shape[1]):
           # Get the longitude and latitude for the curfrom scipy.interpolate import griddatarent sabras cell
           lon = sablons[j]
           lat = sablats[i]
           
           # Find the nearest bigras indices
           lon_idx = find_nearest(biglons, lon)
           lat_idx = find_nearest(biglats, lat)
           
           # Assign the value from bigras to the output raster
           output_raster[i, j] = bigras[lat_idx, lon_idx]
   
   return output_raster

def dstest2(lons, lats, tees, altitudes, dem, lonsb, latsb):

   print(np.shape(tees))
   print(np.array(lons))
   print(np.array(lats))
   print(np.array(altitudes))
   print(lonsb)
   print(latsb)
   print(tees)

   accutemp = np.copy(dem) * 0
   shape2 = np.shape(dem)
   mx2 = shape2[1]
   my2 = shape2[0]
   
   # Varmistetaan, että lons ja lats ovat nousevassa järjestyksessä
   lons1 = np.array(lons)
   lats1 = np.array(lats)
   lons1_sorted_indices = np.argsort(lons1)
   lats1_sorted_indices = np.argsort(lats1)
   
   lons1 = lons1[lons1_sorted_indices]
   lats1 = lats1[lats1_sorted_indices]
   
   tees_sorted = tees[:, lats1_sorted_indices, :][:, :, lons1_sorted_indices]
   
   altitudes1 = altitudes
   
   # Luodaan bikubiset interpolointifunktiot jokaiselle korkeudelle
   interpolators = []
   for k in range(len(altitudes1)):
       interpolator = RectBivariateSpline(lats1, lons1, tees_sorted[k, :, :])
       interpolators.append(interpolator)
   
   for iy in range(my2):
       for ix in range(mx2):
           lat2 = latsb[iy]
           lon2 = lonsb[ix]
           alt2 = dem[iy, ix]
           
           # Etsi lähimmät korkeudet
           altitudedex2 = np.searchsorted(altitudes1, alt2, side='left')
           altitudedex1 = altitudedex2 - 1
           
           # Korjataan indeksit olemaan rajojen sisällä
           altitudedex1 = np.clip(altitudedex1, 0, len(altitudes1) - 2)
           altitudedex2 = np.clip(altitudedex2, 1, len(altitudes1) - 1)
           
           # Lasketaan interpolointifunktioilla lämpötilat jokaiselle korkeudelle
           ta = interpolators[altitudedex1](lat2, lon2)[0, 0]
           tb = interpolators[altitudedex2](lat2, lon2)[0, 0]
           
           # Lineaarinen interpolointi korkeuden perusteella
           baseadelta = altitudes1[altitudedex2] - altitudes1[altitudedex1]
           adelta = alt2 - altitudes1[altitudedex1]
           if baseadelta != 0:
               temp = ((baseadelta - adelta) / baseadelta) * ta + (adelta / baseadelta) * tb
           else:
               temp = ta  # Jos baseadelta on nolla (epätodennäköistä)
           
           accutemp[iy, ix] = temp
   
   return accutemp

def interpola_tosize(indata, shape2):

   kind1='cubic'
   x = np.linspace(0, 1, indata.shape[1])
   y = np.linspace(0, 1, indata.shape[0])
   xx, yy = np.meshgrid(x, y)
   x_new = np.linspace(0, 1, shape2[1])
   y_new = np.linspace(0, 1, shape2[0])
   xx_new, yy_new = np.meshgrid(x_new, y_new)
   outdata = scipy.interpolate.griddata((xx.ravel(), yy.ravel()), indata.ravel(),(xx_new, yy_new), method=kind1)
   return(outdata)

1. Calculate slope and aspect

def calculate_slope_aspect(dem):

   height, width = dem.shape
   dx, dy = np.gradient(dem)

   slope = np.sqrt(dx**2 + dy**2)
   aspect = np.arctan2(dy, -dx) * 180 / np.pi

   # Adjust aspect to be in the range [0, 360]
   aspect = np.where(aspect < 0, 360 + aspect, aspect)

   return slope, aspect

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 3
   2. main program .....

1. iname1="./simu_103/resultat.nc"

iname1="./simu_103/diagfi.nc"

1. imapname1="./bmap1.nc"

imapname1="./map1.nc"

loca1=2600 loca2=loca1+382

shape3=(360,360)

dsin1 = xr.open_dataset(iname1, decode_times=False )

dsmap1 = xr.open_dataset(imapname1, decode_times=False )

map00=dsmap1.z.values

map1=np.flipud(map00)

1. plt.imshow(map1, origin="upper")

1. plt.show()

1. quit(-1)

shape2=np.shape(map1)

from scipy.ndimage import generic_filter

dem=map1

1. demkm=np.flipud(dem)/1000

demkm=dem/1000

sea_mask=np.copy(dem)

sea_mask[sea_mask<0] sea_mask[sea_mask>0]=1

1. sea_mask2=np.array(sea_mask)

slope, aspect = calculate_slope_aspect(dem)

1. Convert slope to degrees

slope_degrees = np.arctan(slope) * 180 / np.pi

distance_to_sea = distance_transform_edt(sea_mask)

mapshape=np.shape(map1) mapwidth=mapshape[1] mapheight=mapshape[0]

dxmap2=180/mapwidth dymap2=90/mapheight

1. 1. h2o_ice_surf

aire=np.array(dsin1.aire.values[:,:])

tass=np.array(dsin1.tsurf.values[loca1:loca2,:,:]) rainss=np.array(dsin1.rain.values[loca1:loca2,:,:]) teess=dsin1.temp.values[loca1:loca2,:,:,:] relhumss=dsin1.RH.values[loca1:loca2,:,:,:] uss=dsin1.u.values[loca1:loca2,:,:,:] vss=dsin1.v.values[loca1:loca2,:,:,:] wss=dsin1.w.values[loca1:loca2,:,:,:]

1. print (np.shape(tass))

tas=np.mean(tass, axis=0) rain=np.mean(rainss, axis=0)*365*24*3600

1. plt.imshow(rain, origin="upper")

1. plt.show()

1. print (np.shape(rain))
2. quit(-1)

shape1=np.shape(tas)

tees=np.mean(teess, axis=0) relhums=np.mean(relhumss, axis=0) us=np.mean(uss, axis=0) vs=np.mean(vss, axis=0) ws=np.mean(wss, axis=0)

print (np.shape(tas))

1. plt.imshow(tas)
2. plt.show()
3. quit(-1)

lats=np.array(dsin1.latitude) lons=np.array(dsin1.longitude)

altitudes=np.array(dsin1.altitude)

print (altitudes)

wkn1=6 # 7 ca 1.8 km

usl1=us[wkn1] vsl1=vs[wkn1] wsl1=ws[wkn1]

wind_horiz1=np.sqrt(usl1*usl1+wsl1*wsl1) wind_angle1=np.arctan2(usl1,vsl1) wind_angledeg1=np.rad2deg(wind_angle1)

wind_angledeg_big_1=interpolate_raster(wind_angledeg1, shape2, "cubic") wind_anglerad_big_1=interpolate_raster(wind_angle1, shape2, "cubic") wind_horiz_big_1=interpolate_raster(wind_horiz1, shape2, "cubic") wind_vert_big_1=interpolate_raster(wsl1, shape2, "cubic")

wind_u_big=interpolate_raster(usl1, shape2, "cubic") wind_v_big=interpolate_raster(vsl1, shape2, "cubic") wind_w_big=interpolate_raster(wsl1, shape2, "cubic")

1. wind_u=interpolate_raster(usl1, shape3, "cubic")
2. wind_v=interpolate_raster(vsl1, shape3, "cubic")
3. wind_w=interpolate_raster(wsl1, shape3, "cubic")
4. plt.imshow(wind_angledeg1)

1. plt.imshow(sea_mask)

1. plt.imshow(distance_to_sea)
2. plt.imshow(aspect)

1. plt.show()
2. quit(-1)

1. 1. downscaled wind
   2. maybe bok
      1. 1. dsu, dsv, dsw=downscale_wind(dem, wind_horiz_big_1, wind_anglerad_big_1, wind_vert_big_1)

dsu, dsv, dsw=downscale_wind(dem, wind_u_big, wind_v_big, wind_w_big)

1. distancewind1=distance_against_wind(dem, dsu, dsv, sea_mask)

angle_range_degrees = 30 # The range within which the wind direction can vary light_angle_degrees = 30 # The angle of the light source above the horizon

1. distancewind0 = distance_against_wind_with_randomness(dem,wind_u_big, wind_v_big, sea_mask, angle_range_degrees)

1. distancewind=distance_against_wind(dem,wind_u_big, wind_v_big, sea_mask)

1. quit(-1)

1. Parameters for randomness, light angle, and shadow intensity

angle_range_degrees = 5 # The range within which the wind direction can vary light_angle_degrees = 10 # The angle of the light source above the horizon shadow_intensity_decrease = 0.99 # Rate at which shadow intensity decreases shadow_intensity_range = (0.3, 3.0) # Range of shadow intensity values

1. Apply Gaussian filter to smooth the terrain

sigma = 25 # Standard deviation for Gaussian kernel. Adjust to control the level of smoothing. smoothed_dem = gaussian_filter(dem, sigma=sigma)

1. shadow_map = calculate_shadow_map(dem,wind_u_big, wind_v_big, light_angle_degrees, angle_range_degrees)

1. shadow_map_smooth = calculate_shadow_map_smooth(smoothed_dem, wind_u_big, wind_v_big, light_angle_degrees, angle_range_degrees, shadow_intensity_decrease, shadow_intensity_range)
2. shadow_map_smooth = calculate_shadow_map_smooth_base(dem, wind_u_big, wind_v_big, light_angle_degrees, angle_range_degrees, shadow_intensity_decrease, shadow_intensity_range)

tpimap=tpi(dem) trimap=tri(dem)

1. plt.imshow(distancewind0)

1. plt.imshow(shadow_map_smooth)

1. plt.show()

1. quit(-1)

wind_ds_angle1=np.arctan2(dsu, dsv) wind_ds_angledeg1=np.rad2deg(wind_ds_angle1)

1. relative_wind_angle_deg=abs(wind_angledeg_big_1-aspect)

relative_wind_angle_deg=abs(wind_ds_angledeg1-aspect)

lonsb=np.linspace(-180+dxmap2,18+0-dxmap2,360) latsb=np.linspace(-90+dymap2,90-dymap2,180)

latras=np.repeat(lats, len(lons)) latras=np.reshape(latras, np.shape(tas)) lonras=np.repeat(lons, len(lats)) lonras=np.rot90(np.reshape(lonras, np.shape(tas)))

latrasb=np.repeat(latsb, len(lonsb)) latrasb=np.flipud(np.reshape(latrasb, np.shape(map1))) lonrasb=np.repeat(lonsb, len(latsb))

1. lonrasb=np.flipud(np.rot90(np.reshape(lonrasb, np.shape(map1))))

lonrasb=np.flipud((np.reshape(lonrasb, np.shape(map1)))) coslatsb=np.cos(np.radians(latrasb))

1. print (np.shape(coslatsb))

1. quit(-1)

1. inputy=[ coslatsb, dem, latrasb, lonrasb]
2. inputy=[ dem, latrasb, lonrasb, relative_wind_angle_deg, distance_to_sea, distancewind, tpimap ]

inputy=[ dem, latrasb, lonrasb, relative_wind_angle_deg, distance_to_sea, dsw ]

accurain=downscale_machinelearn_log(rain, inputy, lons, lats, lonsb, latsb)

1. accurain=get_sized_raster(dem, lonsb, latsb, rain, lons, lats)
2. accurain=get_interpolated_raster(dem, lonsb, latsb, rain, lons, lats)

big_coarserain=interpolate_raster(rain, (40,40), "nearest") saccurain0=interpolate_raster(rain, (40,40), "linear") saccurain=interpolate_raster(saccurain0, shape2, "cubic")

1. plt.imshow(rain, cmap="viridis_r")

1. plt.show()

1. plt.imshow(accurain, cmap="viridis_r")

1. plt.show()

1. plotvar=rain

plotvar=(saccurain+accurain)/2

1. plotvar=big_coarserain ## no interpolation
2. plotvar=saccurain ## 2x interpolated

1. plotvar=accurain

colors = ["#ffffcc", "#a1dab4", "#41b6c4", "#2c7fb8", "#253494"] my_cmap = ListedColormap(colors, name="my_cmap")

my_cmap_r = my_cmap.reversed() mpl.colormaps.register(cmap=my_cmap) mpl.colormaps.register(cmap=my_cmap_r)

cmap1="Blues"

1. cmap1="RdYlBu"
2. cmap1="viridis_r"
3. cmap1="my_cmap"

plt.imshow(plotvar, origin="upper", cmap=cmap1, interpolation="none",vmin=0, vmax=2000, extent=[-180,180,-90,90])

1. 1. plt.imshow(map1)

1. plt.imshow(np.flipud(np.exp(shade1)), cmap="gray",alpha=0.5,extent=[-180,180,-90,90] )

plt.imshow(plotvar, origin="upper", cmap=cmap1,alpha=0.6, interpolation="none",vmin=0, vmax=2000, extent=[-180,180,-90,90]) plt.colorbar()

plt.contour(map1,origin="upper", levels=[0], colors=["k"], lw=30, alpha=0.9, extent=[-180,180,-90,90]) cs=plt.contour(plotvar,origin="upper", levels=[150,300,600,1200,1800], colors=["#0f0000"], alpha=0.7, extent=[-180,180,-90,90]) plt.clabel(cs, fontsize=18)

plt.suptitle("Earth-like planet that has distance 1.03 AU", fontsize=22) plt.title(" Mean rainfall "+str()+" mm/yr", fontsize=18) plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.grid(linestyle=":", color="k", lw=0.2)

1. plt.legend(fontsize=18)
2. plt.contour(map1, levels=[5000], color="brown")

plt.show()