File:Earth like planet rainfall if distance 1p03 au interpolated downscaled 1 1.png
Original file (1,456 × 650 pixels, file size: 329 KB, MIME type: image/png)
Captions
Summary
[edit]DescriptionEarth like planet rainfall if distance 1p03 au interpolated downscaled 1 1.png |
English: Rainfall of earth-like planet, if distance is 1.03 AU |
Date | |
Source | Own work |
Author | Merikanto |
This image is based on eduplanet simulation and diku map. Downscaling with Python3 script. In script many ChatGPT generated functions.
Some codes to produce simulation is in
File:Earth_like_planet_temperature_if_distance_1p03_au_1_1.png
Python script to downscale and visualize rainfall.
-
- extract and process Eduplanet climate data test
-
- python3
-
- 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
- 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
- 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
- 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
- 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
- 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
- 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
- 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)
- 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
- 3
- main program .....
- iname1="./simu_103/resultat.nc"
iname1="./simu_103/diagfi.nc"
- 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)
- plt.imshow(map1, origin="upper")
- plt.show()
- quit(-1)
shape2=np.shape(map1)
from scipy.ndimage import generic_filter
dem=map1
- demkm=np.flipud(dem)/1000
demkm=dem/1000
sea_mask=np.copy(dem)
sea_mask[sea_mask<0]
sea_mask[sea_mask>0]=1
- sea_mask2=np.array(sea_mask)
slope, aspect = calculate_slope_aspect(dem)
- 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
- 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,:,:,:]
- print (np.shape(tass))
tas=np.mean(tass, axis=0)
rain=np.mean(rainss, axis=0)*365*24*3600
- plt.imshow(rain, origin="upper")
- plt.show()
- print (np.shape(rain))
- 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))
- plt.imshow(tas)
- plt.show()
- 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")
- wind_u=interpolate_raster(usl1, shape3, "cubic")
- wind_v=interpolate_raster(vsl1, shape3, "cubic")
- wind_w=interpolate_raster(wsl1, shape3, "cubic")
- plt.imshow(wind_angledeg1)
- plt.imshow(sea_mask)
- plt.imshow(distance_to_sea)
- plt.imshow(aspect)
- plt.show()
- quit(-1)
- downscaled wind
- maybe bok
- 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)
- 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
- distancewind0 = distance_against_wind_with_randomness(dem,wind_u_big, wind_v_big, sea_mask, angle_range_degrees)
- distancewind=distance_against_wind(dem,wind_u_big, wind_v_big, sea_mask)
- quit(-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
- 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)
- shadow_map = calculate_shadow_map(dem,wind_u_big, wind_v_big, light_angle_degrees, angle_range_degrees)
- 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)
- 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)
- plt.imshow(distancewind0)
- plt.imshow(shadow_map_smooth)
- plt.show()
- quit(-1)
wind_ds_angle1=np.arctan2(dsu, dsv)
wind_ds_angledeg1=np.rad2deg(wind_ds_angle1)
- 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))
- 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))
- print (np.shape(coslatsb))
- quit(-1)
- inputy=[ coslatsb, dem, latrasb, lonrasb]
- 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)
- accurain=get_sized_raster(dem, lonsb, latsb, rain, lons, lats)
- 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")
- plt.imshow(rain, cmap="viridis_r")
- plt.show()
- plt.imshow(accurain, cmap="viridis_r")
- plt.show()
- plotvar=rain
plotvar=(saccurain+accurain)/2
- plotvar=big_coarserain ## no interpolation
- plotvar=saccurain ## 2x interpolated
- 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"
- cmap1="RdYlBu"
- cmap1="viridis_r"
- cmap1="my_cmap"
plt.imshow(plotvar, origin="upper", cmap=cmap1, interpolation="none",vmin=0, vmax=2000, extent=[-180,180,-90,90])
- plt.imshow(map1)
- 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)
- plt.legend(fontsize=18)
- plt.contour(map1, levels=[5000], color="brown")
plt.show()
Licensing
[edit]- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 07:00, 7 July 2024 | 1,456 × 650 (329 KB) | Merikanto (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
There are no pages that use this file.
Metadata
This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Software used | |
---|---|
Horizontal resolution | 39.37 dpc |
Vertical resolution | 39.37 dpc |