File:Europe chelsa trace tjuly 14500BP 1.png

Size of this preview: 800 × 530 pixels. Other resolutions: 320 × 212 pixels | 640 × 424 pixels | 1,024 × 679 pixels | 1,328 × 880 pixels.

Original file ‎(1,328 × 880 pixels, file size: 507 KB, MIME type: image/png)

This is a file from the Wikimedia Commons. Information from its description page there is shown below.
Commons is a freely licensed media file repository. You can help.

Summary

DescriptionEurope chelsa trace tjuly 14500BP 1.png	English: Temperature of July, Europe, during 14500 yeas ago. Source of data is CHELSA papeloclimate trace dataset.
Date	16 February 2021
Source	Own work
Author	Merikanto

Camera location	50° 00′ 00″ N, 1° 00′ 00″ E	View this and other nearby images on: OpenStreetMap	50.000000; 1.000000

The source of raster data to produce this image is

CHELSA trace dataset.

https://chelsa-climate.org/ https://envicloud.wsl.ch/#/?prefix=chelsa%2Fchelsa_V1

Data is downloader via script, processed further and visualized with Nasa Panoply netcdf Viewer.

Karger, D.N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R.W., Zimmermann, N.E., Linder, P., Kessler, M. (2017). Climatologies at high resolution for the Earth land surface areas. Scientific Data. 4 170122. https://doi.org/10.1038/sdata.2017.122

Karger D.N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R.W., Zimmermann, N.E,, Linder, H.P., Kessler, M.. Data from: Climatologies at high resolution for the earth’s land surface areas. Dryad Digital Repository.http://dx.doi.org/doi:10.5061/dryad.kd1d4

New data downloader and processer 18.2.2021

NEW data downloader 5.3.2021



## python 3/windows 10

chelsa climate data downloader and downscaling test program

attempt to downscale current gts5 to past gts5
## also externel raster and rgb raster downscaler test
# suppors only chelsa v 1.0, 1.21

WARNING many datasets included, not all
WARNING if chelsa develops, this code will deprecate fast


6.3.2021
0000.0028

import numpy as np
import matplotlib.pyplot as plt
import scipy
from osgeo import gdal
from osgeo import osr
from osgeo import gdalconst
import requests
import os
import glob
from numpy import dtype
import rasterio
from rasterio.plot import show
from rasterio.plot import show_hist
from rasterio.mask import mask
from rasterio.windows import Window
from mpl_toolkits.basemap import shiftgrid
from mpl_toolkits.basemap import Basemap
from shapely.geometry import box
from shapely import speedups
speedups.disable()
import geopandas as gpd
from fiona.crs import from_epsg
import pycrs
import json
import skimage.transform as st
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn import svm, metrics
from pygam import LogisticGAM
from pygam import LinearGAM
from scipy.ndimage import gaussian_filter
import netCDF4
from scipy import exp, mgrid, signal, row_stack, column_stack, tile, misc
from osgeo import gdal
import sys
from struct import *
import pysolar
import pyeto
from climlab.solar.orbital import OrbitalTable
from climlab.solar.insolation import daily_insolation
import calendar
import datetime
from datetime import datetime, timezone

import pyvips





def rm_create_dir(path1):
try:
os.rmdir(path1)
except OSError:
print(" ")
 try:
os.mkdir(path1)
except OSError:
print ("Create dir %s failed" % path1)
else:
print ("Created dir %s " % path1)
 return(0)
def create_dirs():
 rm_create_dir("chelsa")
rm_create_dir("chelsa2")
rm_create_dir("work1")
rm_create_dir("work2")
rm_create_dir("work3")
rm_create_dir("work4")
rm_create_dir("work5")
rm_create_dir("work6")
rm_create_dir("work7")
rm_create_dir("work8")
rm_create_dir("output1")
 return(0)


def mini_create_dirs():
 rm_create_dir("work1")
rm_create_dir("work2")
rm_create_dir("work3")
rm_create_dir("work4")
rm_create_dir("work5")
rm_create_dir("work6")
rm_create_dir("work7")
rm_create_dir("work8")
rm_create_dir("output1")
 return(0)


def save_raster_compress(outname1, rastype1, rasdata1,lon1,lat1,lon2,lat2,cols, rows):
## testef w/gtiff ok
xmin=lon1
ymin=lat1
xmax=lon2
ymax=lat2
nx=cols
ny=rows
xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymax, 0, -yres)
 dst_ds = gdal.GetDriverByName(rastype1).Create(outname1, nx, ny, 1, gdal.GDT_Float32, options=['COMPRESS=DEFLATE'])
dst_ds.SetGeoTransform(geotransform) # specify coords
srs = osr.SpatialReference() # establish encoding
srs.ImportFromEPSG(4326) # lat/long
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
dst_ds.GetRasterBand(1).WriteArray(rasdata1) # write r-band to the raster
dst_ds.FlushCache() # write to disk
dst_ds = None
return(0)
def gaussian_blur_sub(in_array, size):
ny, nx = in_array.shape
ary = row_stack((tile(in_array[0,:], size).reshape((size, nx)),in_array,tile(in_array[-1,:], size).reshape((size, nx))))
arx = column_stack((tile(ary[:,0], size).reshape((size, ny + size*2)).T,ary,tile(ary[:,-1], size).reshape((size, ny + size*2)).T))
x, y = mgrid[-size:size+1, -size:size+1]
g = exp(-(x**2/float(size) + y**2/float(size)))
g = (g / g.sum()) #.astype(in_array.dtype)
return signal.convolve(arx, g, mode='valid')
def gaussian_blur(src1, dst1, size1):
source = gdal.Open(src1)
gt = source.GetGeoTransform()
proj = source.GetProjection()
band_array = source.GetRasterBand(1).ReadAsArray()
blurredArray = gaussian_blur_sub(band_array, size1)
#driver = gdal.GetDriverByName('GTIFF')
driver = gdal.GetDriverByName('NetCDF')
driver.Register()
cols = source.RasterXSize
rows = source.RasterYSize
bands = source.RasterCount
band = source.GetRasterBand(1)
datatype = band.DataType
output = driver.Create(dst1, cols, rows, bands, datatype)
output.SetGeoTransform(gt)
output.SetProjection(proj)
outBand = output.GetRasterBand(1)
outBand.WriteArray(blurredArray, 0, 0)
source = None # close raster
output = None
band = None
outBand = None
return(0)


def rasterize_shapefile(shpath1, outfile1, lon1, lat1, lon2, lat2, cols, rows):
outfile2=outfile1.replace("tif","nc")
tempofile1="temp1.tif"
kom1="gdal_rasterize -burn 255 -ot Float32"
kom2=" -te "+str(lon1)+" "+str(lat1)+" "+str(lon2)+" "+str(lat2)
kom3=" -ts "+str(cols)+" "+str(rows)
kom4=" "+shpath1+" "+tempofile1
kommandoo1=kom1+kom2+kom3+kom4
koma2="gdal_calc.py -A "+tempofile1+ " --calc=\"(A>200)*99\" --outfile=temp2.tif"
koma3="gdal_calc.py -A temp2.tif --calc=\"(A<1)*1\" --outfile="+outfile1
 #koma4="gdal_calc.py -A temp3.tif --calc=\"(A>97)*0\" --outfile=temp4.tif"
 kommandoo2="gdal_translate -of netcdf -ot Float32 "+outfile1+" "+tempofile1
 print (kommandoo1)
#print (kommandoo2)
os.system(kommandoo1)
os.system(koma2)
os.system(koma3)
#os.system(koma4)
#river_proximity("temp3.tif", "trivo.tif")
 os.system(kommandoo2)
return(0)
def chelsa_file_name(variable1, year1,month1):
fana1="nonefile_error"
 if(month1==0):
# annual
base1="CHELSA"
base2="_V1.2.1.tif"
variable2=variable1
fana1=base1+"_"+variable2+"_"+str(year1)+base2
else:
if (variable1=="bio"):
## bio, not month
base1="CHELSA"
base2="_V1.2.1.tif"
variable2=variable1
fana1=base1+"_"+variable2+str(month1).zfill(2)+"_"+str(year1)+base2
else:
## monthly variable
variable2=variable1
base1="CHELSA"
base2="_V1.2.1.tif"
fana1=base1+"_"+variable2+"_"+str(year1)+"_"+str(month1).zfill(2)+base2
 return(fana1)


def chelsa_timeseries_name(variable1, year1,month1):
fana1="nonefile_error"
 if(month1==0):
# annual
base1="CHELSA"
base2="_V1.2.1.tif"
variable2=variable1
fana1=base1+"_"+variable2+"_"+str(year1)+base2
else:
if (variable1=="bio"):
## bio, not month
base1="CHELSA"
base2="_V1.2.1.tif"
variable2=variable1
fana1=base1+"_"+variable2+str(month1).zfill(2)+"_"+str(year1)+base2
else:
## monthly variable
variable2=variable1
base1="CHELSA"
base2="_V1.2.1.tif"
fana1=base1+"_"+variable2+"_"+str(year1)+"_"+str(month1).zfill(2)+base2
 return(fana1)
def chelsa_simu_name(simu1,variable1, year1,month1):
fana1="nonefile_error"
base1="CHELSA"+"_"+simu1
base2="_1.tif"
 simu2=simu1
 if(simu1=="NCAR_CCSM4"):
simu2="PMIP_CCSM4"
 if(simu1=="PMIP_CCSM4"):
simu1="NCAR_CCSM4"
simu2="PMIP_CCSM4"
 base3="CHELSA"+"_"+simu2


 if(month1==0):
# annual
variable2=variable1
fana1=base3+"_"+variable2+"_"+str(year1)+base2
else:
if (variable1=="bio"):
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/pmip3/NCAR_CCSM4/CHELSA_PMIP_CCSM4_BIO_01.tif
## bio, not month
variable2="BIO"
fana1=base3+"_"+variable2+"_"+str(month1).zfill(2)+".tif"
else:
## monthly variable
variable2=variable1
if (variable2=="prec"):
if(month1==1):
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/pmip3/NCAR_CCSM4/CHELSA_PMIP_CCSM4_prec_01_1.tif
#fana1=base3+"_"+variable2+"_"+str(month1).zfill(2)+base2
fana1=base3+"_"+variable2+"_"+str(month1).zfill(2)+base2
else:
fana1=base3+"_"+variable2+"_"+str(month1)+base2
else:
fana1=base3+"_"+variable2+"_"+str(month1)+base2


 return(fana1)


def chelsa_trace_name(variable1, year1,month1):
fana1="nonefile_error"
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/pr/CHELSA_TraCE21k_pr_10_-100_V1.0.tif
base1="CHELSA_TraCE21k"
base2="_V1.0.tif"


 year2=int(year1/100)*-1
 #https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/tasmax/CHELSAtrace_tasmax_7_-145_V1.0.tif
 if(month1==0):
if (variable1=="orog"):
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/orog/CHELSA_TraCE21k_dem_-101_V1.0.tif
 #https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/orog/glacier_plus_dem_0_V1.2.tif
fana1="CHELSA_TraCE21k_dem_"+str(year2*-1)+"_V1.2.tif"
return(fana1)
else:
# annual
variable2=variable1
fana1=base1+"_"+variable2+"_"+str(year2)+base2
 else:
if (variable1=="bio"):
## bio, not month
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/bio/CHELSA_bio01_-100_V1.2.1.tif
##print("BIIB MAYBE NOK")
base1="CHELSA"
 base2="_V1.2.1.tif"
variable2="bio"
fana1=base1+"_"+variable2+str(month1).zfill(2)+"_"+str(year2)+base2
else:
## monthly variable
variable2=variable1
fana1=base1+"_"+variable2+"_"+str(month1)+"_"+str(year2)+base2
 return(fana1)


def chelsa_cruts_name(variable1, year1,month1):
fana1="nonefile_error"
 #https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_cruts/tmax/CHELSAcruts_tmax_10_1901_V.1.0.tif
base1="CHELSAcruts"
base2="_V.1.0.tif"
if(month1==0):
# annual NA NOK
 variable2=variable1
fana1=base1+"_"+variable2+"_"+str(year1)+base2
else:
if (variable1=="bio"):
 ## bio, not month NA NOK
variable2=variable1
fana1=base1+"_"+variable2+str(month1)+"_"+str(year1)+base2
else:
## monthly variable
variable2=variable1
fana1=base1+"_"+variable2+"_"+str(month1)+"_"+str(year1)+base2
 return(fana1)


def chelsa_climatologies_name(variable1, year1,month1):
fana1="nonefile_error"
 base1="CHELSA"
base2="_1979-2013_V1.2_land.tif"
 variable2=variable1
if(variable1=="tmean"):
variable2="temp10"
if(variable1=="tmax"):
variable2="temp10"
if(variable1=="tmin"):
variable2="temp10"
 if(variable1=="prec"):
base2="_V1.2_land.tif"
 #https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/climatologies/prec/CHELSA_prec_01_V1.2_land.tif
 if(month1<1):
print("Month below 0")
return(0)
 if(month1==0):
# annual NOK
variable2=variable1
fana1=base1+"_"+variable2+"_"+base2
else:
if (variable1=="bio"):
 ## bio, not month
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/climatologies/bio/CHELSA_bio10_12.tif
variable2="bio10"
fana1=base1+"_"+variable2+"_"+str(month1).zfill(2)+".tif"
else:
#print("CC Month")
#print(month1)
## monthly variable
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/climatologies/tmean/CHELSA_temp10_01_1979-2013_V1.2_land.tif
fana1=base1+"_"+variable2+"_"+str(month1).zfill(2)+base2
#print(fana1)


return(fana1)
def save_raster(outname1, rastype1, rasdata1,lon1,lat1,lon2,lat2,cols, rows):
## testef w/gtiff ok
xmin=lon1
ymin=lat1
xmax=lon2
ymax=lat2
nx=cols
ny=rows
xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymax, 0, -yres)
 dst_ds = gdal.GetDriverByName(rastype1).Create(outname1, nx, ny, 1, gdal.GDT_Float32, options=['COMPRESS=DEFLATE'])
dst_ds.SetGeoTransform(geotransform) # specify coords
srs = osr.SpatialReference() # establish encoding
srs.ImportFromEPSG(4326) # lat/long
dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
dst_ds.GetRasterBand(1).WriteArray(rasdata1) # write r-band to the raster
dst_ds.FlushCache() # write to disk
dst_ds = None
return(0)
def chelsa_exchelsa_name(simu1, variable1, year1,month1):
fana1="nonefile_error"
 base1="CHELSA"
base2=".tif"
 ## note change of variables!!!
variable2=variable1
variable1=simu1
 if(variable1=="pet"):
base2="_1979-2013.tif"
if(variable1=="pet"):
base2="_1979-2013.tif"
 if(variable2 in("dfg", "gsl","gst", "lgd")):
base2="_V1.2.1.sdat.tif"




if(month1==0):
# annual NOK
fana1=base1+"_"+variable2+"_"+str(year1)+base2
else:
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/exchelsa/srad/CHELSA_stot_7.tif
fana1=base1+"_"+variable2+"_"+str(month1)+base2
 return(fana1)


def get_chelsa_variable(repo1, simu1, variable1, year1, month1):
headers1 = {
'User-Agent': 'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/46.0.2490.80 Safari/537.36',
'Content-Type': 'text/html',
}
 baseurl1="https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1"
url0=baseurl1+"/"+repo1+"/"
 if(repo1=="timeseries"):
 pathname1=variable1
filename1=chelsa_file_name(variable1, year1,month1)
url1=url0+pathname1+"/"+filename1
outfilename1="./chelsa/"+filename1
print(url1)
print(outfilename1)
 if(repo1=="pmip3"):
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/pmip3/NCAR_CCSM4/CHELSA_PMIP_CCSM4_tmean_7_1.tif
 if(simu1=="DEM"):
url1="https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/pmip3/DEM/high_longlat.tif"
outfilename1="./chelsa/high_longlat.tif"
else:
pathname1=simu1
basename0=chelsa_simu_name(simu1,variable1, year1,month1)
filename1=basename0
url1=url0+pathname1+"/"+filename1
outfilename1="./chelsa/"+filename1
print(url1)
print(outfilename1)
 if(repo1=="chelsa_trace"):
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/tasmax/CHELSAtrace_tasmax_7_-145_V1.0.tif
pathname1=variable1
filename1=chelsa_trace_name(variable1, year1,month1)
url1=url0+pathname1+"/"+filename1
outfilename1="./chelsa/"+filename1
print(url1)
print(outfilename1)
 if(repo1=="chelsa_cruts"):
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_cruts/tmax/CHELSAcruts_tmax_10_1901_V.1.0.tif
pathname1=variable1
filename1=chelsa_cruts_name(variable1, year1,month1)
url1=url0+pathname1+"/"+filename1
outfilename1="./chelsa/"+filename1
print(url1)
print(outfilename1)
 if(repo1=="climatologies"):
pathname1=variable1
filename1=chelsa_climatologies_name(variable1, year1,month1)
url1=url0+pathname1+"/"+filename1
outfilename1="./chelsa/"+filename1
print(url1)
print(outfilename1)
 if(repo1=="exchelsa"):
## !! note use simu1 as folder, variable1 as variable
pathname1=simu1
filename1=chelsa_exchelsa_name(simu1, variable1, year1,month1)
url1=url0+pathname1+"/"+filename1
outfilename1="./chelsa/"+filename1
print(url1)
print(outfilename1)
 r1 = requests.get(url1, headers=headers1)
if(r1==0):
printf("Chelsa loading error")
return(-1)
 with open(outfilename1, "wb") as code:
code.write(r1.content)
 return(0)


def get_chelsa_trace_data(year1):
 for n in range(1,13):
get_chelsa_variable("chelsa_trace","", "pr",year1, n)
 for n in range(1,20):
get_chelsa_variable("chelsa_trace","", "bio", year1, n)
 for n in range(1,13):
get_chelsa_variable("chelsa_trace","", "tasmax", year1, n)
get_chelsa_variable("chelsa_trace","", "tasmin", year1, n)
 get_chelsa_variable("chelsa_trace","", "orog", year1, 0)
 return(0);


def getFeatures(gdf):
return [json.loads(gdf.to_json())['features'][0]['geometry']]
def clip_raster_to_outfile(inf1, outf1, lon1, lat1, lon2, lat2):
print("clip rout")
try:
data = rasterio.open(inf1)
except:
print("Dscaler.py: Raster file error, maybe corrupt")
return(-1)
print("Data")
print(data)
print(".. data")
 bbox=box(lon1, lat1, lon2, lat2)
geo=gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=from_epsg(4326))
coords = getFeatures(geo)
#print(coords)
out_img, out_transform = mask(data, shapes=coords, crop=True)
out_meta = data.meta.copy()
#epsg_code = int(data.crs.data['init'][1:])
epsg_code=int(4326)
height=out_img.shape[1]
width=out_img.shape[2]
#print(height)
#print (width)
out_meta.update({"driver": "GTiff","height": out_img.shape[1],"width": out_img.shape[2],"transform": out_transform,"crs": pycrs.parse.from_epsg_code(epsg_code).to_proj4()})
try:
with rasterio.open(outf1, "w", **out_meta) as dest:
dest.write(out_img)
except:
print("TIF write error")
 #clipped = rasterio.open(outf1)
return(0)


def get_chelsa_timeseries_variable(year1, variable1):
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/timeseries/tmean/CHELSA_tmean_1979_01_V1.2.1.tif
variable2=variable1
print("Timeseries variable ", year1,variable1)
 urlbase1="https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/timeseries/"
base1="CHELSA"
base2="_V1.2.1.tif"
headers1 = {
'User-Agent': 'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/46.0.2490.80 Safari/537.36',
'Content-Type': 'text/html',
}
 for n in range(1,13):
geturl1=urlbase1+variable1+"/"+base1+"_"+variable2+"_"+str(year1)+"_"+str(n).zfill(2)+base2
filename1="./chelsa/"+base1+"_"+variable2+"_"+str(year1)+"_"+str(n).zfill(2)+base2
print(geturl1)
print(filename1)
r1 = requests.get(geturl1, headers=headers1)
if(r1==0):
printf("Chelsa loading error")
return(-1)
 with open(filename1, "wb") as code:
code.write(r1.content)
return(0)


def chelsa_file_name_b(variable1, year1,month1):
fana1="nonefile_error"
 if(month1==0):
# annual
base1="CHELSA"
base2="_V1.2.1.tif"
variable2=variable1
fana1=base1+"_"+variable2+"_"+str(year1)+base2
else:
if (variable1=="bio"):
## bio, not month
base1="CHELSA"
base2="_V1.2.1.tif"
variable2=variable1
fana1=base1+"_"+variable2+str(month1).zfill(2)+"_"+str(year1)+base2
else:
## monthly variable
variable2=variable1
base1="CHELSA"
base2="_V1.2.1.tif"
fana1=base1+"_"+variable2+"_"+str(year1)+"_"+str(month1).zfill(2)+base2


return(fana1)


def cut_and_scale_chelsa_variable(repo1, simu1, variable1, year1, month1, lon1, lat1, lon2, lat2, width1, height1, width2, height2):


chelsadir1="./chelsa/"
chelsadir2="./work1/"
chelsadir3="./work2/"
chelsadir4="./work3/"
chelsadir5="./work4/"
 chen1=chelsa_timeseries_name(variable1, year1,month1)
 if(repo1=='timeseries'):
chen1=chelsa_timeseries_name(variable1, year1,month1)
 if(repo1=='trace'):
print("trace")
chen1=chelsa_trace_name(variable1, year1,month1)


inf1=chelsadir1+chen1
outf1=chelsadir2+chen1
outf2=chelsadir3+chen1
outf3=chelsadir4+chen1
outf3=outf3.replace(".tif", ".nc")
outf4=chelsadir5+chen1
outf5=outf4.replace(".tif", ".nc")
print(inf1)
print("clip")
clip_raster_to_outfile(inf1, outf1, lon1, lat1, lon2, lat2)
ds=0
print("warp")
ds = gdal.Warp(outf2, outf1, width=width1, height=height1, resampleAlg='bilinear')
ds=0
ds = gdal.Warp(outf4, outf1, width=width2, height=height2, resampleAlg='bilinear')
ds = 0
print("translate")
ds = gdal.Translate(outf3, outf2, format='NetCDF')
ds = gdal.Translate(outf5, outf4, format='NetCDF')
 ds=0
return(0)
def random_forest_multiple_vars( x1, y1, x2):
print ("RF /2")
xa1=np.array(x1).astype(float)
ya1=np.array(y1).astype(float)
xa2=np.array(x2).astype(float)
 #print (xa1)

quit(-1)

 sha1=np.shape(x1)
 dim2=sha1[1]
 x=xa1.reshape((-1, dim2))
#y=ya1.reshape((-1, 1))
y=ya1.ravel()
xb=xa2.reshape((-1, dim2))


#model = LinearRegression().fit(x, y)
#degree=3
#polyreg=make_pipeline(PolynomialFeatures(degree),LinearRegression())
#model=polyreg.fit(x,y)
 sc = StandardScaler()
xa = sc.fit_transform(x)
xba = sc.transform(xb)


model = RandomForestRegressor(n_estimators=10, max_features=2).fit(xa,y)


y2= model.predict(xba)
return(y2)


def linear_regression_singe_var( x1, y1, x2):
#print (x1)
#print(y1)
#return(0)
 #xa1=np.array(x1)
#ya1=np.array(y1)
#xa2=np.array(x2)
xa1=np.array(x1)
ya1=np.array(y1)
xa2=np.array(x2)
 sha1=np.shape(x1)
 #dim2=sha1[1]
 #x=xa1.reshape((-1, dim2))
#y=ya1.reshape((-1, 1))
#xb=xa2.reshape((-1, dim2))
x=xa1.reshape((-1, 1))
y=ya1.reshape((-1, 1))
xb=xa2.reshape((-1, 1))
 #print (x)
#print (y)
 #model = LinearRegression().fit(x, y)
#model = LogisticGAM().fit(x, y)
#model = LinearGAM().fit(x, y)
#model = RandomForestRegressor(n_estimators=10, max_features=2).fit(x,y)
degree=2
polyreg=make_pipeline(PolynomialFeatures(degree),LinearRegression())
#polyreg=make_pipeline(PolynomialFeatures(degree), )
model=polyreg.fit(x,y)
 ## warning not xb
y2= model.predict(xb)
#print(y2)
#print("LR")
return(y2)


def linear_regression_multiple_vars( x1, y1, x2):
print ("MLR")
xa1=np.array(x1)
ya1=np.array(y1)
xa2=np.array(x2)
 sha1=np.shape(x1)
 dim2=sha1[1]
 x=xa1.reshape((-1, dim2))
y=ya1.reshape((-1, 1))
xb=xa2.reshape((-1, dim2))


#model = LinearRegression().fit(x, y)
degree=3
polyreg=make_pipeline(PolynomialFeatures(degree),LinearRegression())
model=polyreg.fit(x,y)


y2= model.predict(xb)
return(y2)
def gam_multiple_vars( x1, y1, x2):
print ("GAM")
xa1=np.array(x1)
ya1=np.array(y1)
xa2=np.array(x2)
 #print (xa1)

quit(-1)

 sha1=np.shape(x1)
 dim2=sha1[1]
 x=xa1.reshape((-1, dim2))
y=ya1.reshape((-1, 1))
xb=xa2.reshape((-1, dim2))
 #sc = StandardScaler()
#xa = sc.fit_transform(x)
#xba = sc.transform(xb)
 #model = LogisticGAM().fit(x, y)
model = LinearGAM().fit(x, y)
 y2= model.predict(xb)
 return(y2)


def array2raster(newRasterfn,rasterfn,array):

   raster = gdal.Open(rasterfn)
   geotransform = raster.GetGeoTransform()
   originX = geotransform[0]
   originY = geotransform[3]
   pixelWidth = geotransform[1]
   pixelHeight = geotransform[5]
   cols = array.shape[1]
   rows = array.shape[0]

   driver = gdal.GetDriverByName('GTiff')
   outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
   outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
   outband = outRaster.GetRasterBand(1)
   outband.WriteArray(array)
   outRasterSRS = osr.SpatialReference()
   outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
   outRaster.SetProjection(outRasterSRS.ExportToWkt())
   outband.FlushCache()

def array2nc(newRasterfn,rasterfn,array):

   raster = gdal.Open(rasterfn)
   geotransform = raster.GetGeoTransform()
   originX = geotransform[0]
   originY = geotransform[3]
   pixelWidth = geotransform[1]
   pixelHeight = geotransform[5]
   cols = array.shape[1]
   rows = array.shape[0]

   driver = gdal.GetDriverByName('NetCDF')
   #outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
   outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
   outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
   outband = outRaster.GetRasterBand(1)
   outband.WriteArray(array)
   outRasterSRS = osr.SpatialReference()
   outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
   outRaster.SetProjection(outRasterSRS.ExportToWkt())
   outband.FlushCache()
   return(0)

def get_chelsa_timeseries_bundle(year1):
 for n in range(1,13):
get_chelsa_variable("timeseries","", "tmean", year1, n)
get_chelsa_variable("timeseries","", "prec", year1, n)


for n in range(1,20):
get_chelsa_variable("timeseries","", "bio", year1, n)


get_chelsa_variable("timeseries","", "gts5", year1, 0)
 return(0)


def cut_chelsa_bundle(repo1, simu1, year1, month1, lon1, lat1, lon2, lat2, width1, height1, width2, height2):
#gdal.SetConfigOption('CPL_LOG', 'NUL')
global set_biovars_off
 if(repo1=='timeseries'):
variable1="gts5"
cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, 0, lon1, lat1, lon2, lat2, width1, height1, width2, height2)


if (repo1=='timeseries'):
for n in range(1,13):
print("Variables ", n)
variable1="tmean"
cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, n, lon1, lat1, lon2, lat2, width1, height1, width2, height2)
variable1="prec"
cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, n, lon1, lat1, lon2, lat2, width1, height1, width2, height2)
 if (set_biovars_off==0):
for n in range(1,20):
print("Bio ", n)
variable1="bio"
cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, n, lon1, lat1, lon2, lat2, width1, height1, width2, height2)
 if (repo1=='trace'):
 for n in range(1,13):
print("Trace variables ", n)
variable1="tasmax"
cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, n, lon1, lat1, lon2, lat2, width1, height1, width2, height2)
variable1="tasmin"
cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, n, lon1, lat1, lon2, lat2, width1, height1, width2, height2)
variable1="pr"
cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, n, lon1, lat1, lon2, lat2, width1, height1, width2, height2)




if (set_biovars_off==0):
#for n in range(8,9):
for n in range(1,20):
print("Trace bio ", n)
variable1="bio"
cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, n, lon1, lat1, lon2, lat2, width1, height1, width2, height2)


 #cut_and_scale_chelsa_variable(repo1, simu1,variable1, year1, 0, lon1, lat1, lon2, lat2, width1, height1)
 return(0)


def downscale_chelsa_bundle(repo1, simu1, year1):
 print("Downscale test ...")
 bigs=[]
smalls=[]
 #variable1="tmean"
 for n in range(1,13):
#for n in range(6,8):
#for n in range(1,13,3):
print("Variables ", n)
 month1=n
maxname1="./work2/"+chelsa_timeseries_name(variable1, year1,month1)
minname1="./work4/"+chelsa_timeseries_name(variable1, year1,month1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 #print (np.shape(maxarray1))
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls.append(minna1)
bigs.append(maxa1)


#print("Target rasters ...")
month1=0
mintargetname1="./work4/"+chelsa_timeseries_name("gts5", year1,month1)
month1=0
maxtargetname1="./work2/"+chelsa_timeseries_name("gts5", year1,month1)
 maxtarget1 = gdal.Open(maxtargetname1)
maxtargetarray1 = np.array(maxtarget1.GetRasterBand(1).ReadAsArray())
mintarget1 = gdal.Open(mintargetname1)
mintargetarray1 = np.array(mintarget1.GetRasterBand(1).ReadAsArray())
 minnatarget1=mintargetarray1.flatten()
maxatarget1=maxtargetarray1.flatten()
 apoints1=list(zip(*smalls))
bpoints1=minnatarget1
cpoints1=list(zip(*bigs))
 #y2=gam_multiple_vars(apoints1, bpoints1, cpoints1)
#y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints1)
y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints1)
 youtarray1=y2.reshape( bheight1, bwidth1).astype(float)


 array2nc("./output1/downscaled.nc",maxtargetname1,youtarray1)
 return(0)




def tiff2nc(fromraster, toraster):

   raster = gdal.Open(fromraster)
   geotransform = raster.GetGeoTransform()
   originX = geotransform[0]
   originY = geotransform[3]
   pixelWidth = geotransform[1]
   pixelHeight = geotransform[5]
   array=raster.ReadAsArray()
   cols = array.shape[1]
   rows = array.shape[0]

   driver = gdal.GetDriverByName('NetCDF')

   outRaster = driver.Create(toraster, cols, rows, 1, gdal.GDT_Float32)
   outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
   outband = outRaster.GetRasterBand(1)
   outband.WriteArray(array)
   outRasterSRS = osr.SpatialReference()
   outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
   outRaster.SetProjection(outRasterSRS.ExportToWkt())
   outband.FlushCache()
   return(0)





def downscale_neo(repo1, simu1, year1):
print("NEO downscale test ...")
bigs=[]
smalls=[]
 #variable1="tmean"
for variable1 in ("tmean", "prec"):
print (variable1)
for n in range(1,13):
print("Variables ", n)
 month1=n
maxname1="./work2/"+chelsa_timeseries_name(variable1, year1,month1)
minname1="./work4/"+chelsa_timeseries_name(variable1, year1,month1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 #print (np.shape(maxarray1))
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls.append(minna1)
bigs.append(maxa1)


#print("Target rasters ...")
month1=0
#mintargetname1="./work4/"+chelsa_timeseries_name("gts5", year1,month1)
mintargetname1="./work4/neo_ndvi.tif"
month1=0
maxtargetname1="./work2/neo_ndvi.tif"
#maxtargetname1="./work2/"+chelsa_timeseries_name("gts5", year1,month1)
 maxtarget1 = gdal.Open(maxtargetname1)
maxtargetarray1 = np.array(maxtarget1.GetRasterBand(1).ReadAsArray())
mintarget1 = gdal.Open(mintargetname1)
mintargetarray1 = np.array(mintarget1.GetRasterBand(1).ReadAsArray())
 minnatarget1=mintargetarray1.flatten()
maxatarget1=maxtargetarray1.flatten()
 ## note warning this is data filter, very specific to dataset
## clamps to max, min values
 minnatarget1 = np.where(minnatarget1 < 0, 0, minnatarget1)
maxatarget1 = np.where(maxatarget1 < 0, 0, maxatarget1)
minnatarget1 = np.where(minnatarget1 > 1, 0, minnatarget1)
maxatarget1 = np.where(maxatarget1 > 1, 0, maxatarget1)
 apoints1=list(zip(*smalls))
bpoints1=minnatarget1
cpoints1=list(zip(*bigs))
 print("Downscaling ....")
 #y2=gam_multiple_vars(apoints1, bpoints1, cpoints1)
##y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints1)
y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints1)
 youtarray1=y2.reshape( bheight1, bwidth1).astype(float)


 array2nc("./output1/downscaled_neo_current.nc",maxtargetname1,youtarray1)
 return(0)


def process_neo_raster(inputfile1, outputfile0, referencefile0):
 outputfile1="./work2/"+outputfile0
outputfile2="./work4/"+outputfile0
 referencefile1="./work2/"+referencefile0
referencefile2="./work4/"+referencefile0
 raster_reso_to_reference_raster(inputfile1, referencefile1, outputfile1)
 raster_reso_to_reference_raster(inputfile1, referencefile2, outputfile2)
 outputnc1=outputfile1.replace(".tif", ".nc")
outputnc2=outputfile2.replace(".tif", ".nc")
 tiff2nc(outputfile1, outputnc1)
tiff2nc(outputfile2, outputnc2)
 return(0);




def calculate_tmean_from_trace_tasmax_tasmin(repo1, simu1, year2, lon1, lat1, lon2, lat2, width1, height1, width2, height2):

nok

 #for n in range(1,13):
for n in range(1,13):
print(n)
name1="./work2/"+chelsa_trace_name("tasmax", year2, n)
name2="./work2/"+chelsa_trace_name("tasmin", year2, n)
name3="./work2/"+chelsa_trace_name("tmean", year2, n)
name4=name3
name4=name4.replace("tif", "nc")
name5="./work4/"+chelsa_trace_name("tasmax", year2, n)
name6="./work4/"+chelsa_trace_name("tasmin", year2, n)
name7="./work4/"+chelsa_trace_name("tmean", year2, n)
name8=name7
name8=name8.replace("tif", "nc")
inimage1, lats1, lons1=load_raster_float32(name1)
inimage2, lats1, lons1=load_raster_float32(name2)
outimage1=(inimage1+inimage2)/2.0
outimage1f=np.array(outimage1).astype(float)
print(name3)
saveraster_float32(name3,"GTiff", outimage1f, lon1, lon2, lat1, lat2)
#saveraster_float32(name4,"NetCDF4", outimage1f, lon1, lon2, lat1, lat2)
jnimage1, lats1, lons1=load_raster_float32(name5)
jnimage2, lats1, lons1=load_raster_float32(name6)
joutimage1=(jnimage1+jnimage2)/2.0
joutimage1f=np.array(joutimage1).astype(float)
print(name7)
saveraster_float32(name7,"GTiff", joutimage1f, lon1, lon2, lat1, lat2)
#saveraster_float32(name8,"NetCDF4", outimage1f, lon1, lon2, lat1, lat2)


print ("DEBUG BREAK")
#quit(-1)
 return(0)


def trace_tmean_from_tasmax_tasmin(slice1):
 for n in range(1,13):
 aoutname1="./work2/"+chelsa_trace_name("tmean", slice1,n)
ainame1="./work2/"+chelsa_trace_name("tasmax", slice1,n)
ainame2="./work2/"+chelsa_trace_name("tasmin", slice1,n)
 boutname1="./work4/"+chelsa_trace_name("tmean", slice1,n)
biname1="./work4/"+chelsa_trace_name("tasmax", slice1,n)
biname2="./work4/"+chelsa_trace_name("tasmin", slice1,n)


kommand1="gdal_calc.py --calc=\"(A+B)/2\" --outfile="+aoutname1+" -A "+ainame1+" -B "+ainame2
kommand2="gdal_calc.py --calc=\"(A+B)/2\" --outfile="+boutname1+" -A "+biname1+" -B "+biname2
 os.system(kommand1)
os.system(kommand2)


def downscale_timeseries_trace(outname1,repo1, repo2, simu1, year1, year2):
print("Timeseries vs trace downscale test ...")
 slice1=int(year2/100)
 bigs1=[]
smalls1=[]
 variable1='tmean'
print("Timeseries")
if (variable1=='tmean'):
#for variable1 in ("tmean", "prec"):
#for variable1 in ("tasmax", "tasmin"):
print (variable1)
#for n in range(1,13):
for n in range(4,10):
print(" TS variables ", n)
 month1=n
maxname1="./work2/"+chelsa_timeseries_name(variable1, year1,month1)
minname1="./work4/"+chelsa_timeseries_name(variable1, year1,month1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 #print (np.shape(maxarray1))
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls1.append(minna1)
bigs1.append(maxa1)


bigs2=[]
smalls2=[]
if(variable1=='tmean'):
#for variable1 in ("tmean", "pr"):
#for variable1 in ("tasmax", "tasmin"):
print (variable1)
#for n in range(1,13):
for n in range(4,10):
print("n TRACE variables ", n)
 month1=n
maxname1="./work2/"+chelsa_trace_name(variable1, year2,month1)
minname1="./work4/"+chelsa_trace_name(variable1, year2,month1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 #print (np.shape(maxarray1))
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls2.append(minna1)
bigs2.append(maxa1)


 #print("Target rasters ...")
month1=0
mintargetname1="./work4/"+chelsa_timeseries_name("gts5", year1,month1)
maxtargetname1="./work2/"+chelsa_timeseries_name("gts5", year1,month1)
 #mintargetname1="./work4/neo_ndvi.tif"
#maxtargetname1="./work2/neo_ndvi.tif"
 #mintargetname1="./work4/neo_rgb.tif"
#maxtargetname1="./work2/neo_rgb.tif"
 #loadrgb(filename)
#writergb(salbuname, filename,arrays)
 maxtarget1 = gdal.Open(maxtargetname1)
maxtargetarray1 = np.array(maxtarget1.GetRasterBand(1).ReadAsArray())
mintarget1 = gdal.Open(mintargetname1)
mintargetarray1 = np.array(mintarget1.GetRasterBand(1).ReadAsArray())
 minnatarget1=mintargetarray1.flatten()
maxatarget1=maxtargetarray1.flatten()
 ## note warning this is data filter, very specific to dataset
## clamps to max, min values
 #minnatarget1 = np.where(minnatarget1 < 0, 0, minnatarget1)
#maxatarget1 = np.where(maxatarget1 < 0, 0, maxatarget1)
#minnatarget1 = np.where(minnatarget1 > 1, 0, minnatarget1)
#maxatarget1 = np.where(maxatarget1 > 1, 0, maxatarget1)
 apoints1=list(zip(*smalls1))
bpoints1=minnatarget1
cpoints1=list(zip(*bigs1))
 apoints2=list(zip(*smalls2))
bpoints2=minnatarget1
cpoints2=list(zip(*bigs2))
 print("Downscaling ....")
 ## test only
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints1)
##y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints1) ## nogood
#y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints1) ##ok


#y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints2) ##NOK
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints2)
y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints2) ##ok


youtarray1=y2.reshape( bheight1, bwidth1).astype(float)


 array2nc(outname1,maxtargetname1,youtarray1)
#array2raster("./output1/downscaled.tiff",maxtargetname1,youtarray1)
 return(0)






def loadrgb(filename):

   raster = gdal.Open(filename)
   geotransform = raster.GetGeoTransform()
   originX = geotransform[0]
   originY = geotransform[3]
   pixelWidth = geotransform[1]
   pixelHeight = geotransform[5]
   cols = array.shape[1]
   rows = array.shape[0]
   r=raster.ReadAsArray(1)
   g=raster.ReadAsArray(2)
   b=raster.ReadAsArray(3)
   r1=r.astype(float)
   g2=g.astype(float)
   b2=b.astype(float)
   rasa1=[]
   rasa1.append(r1)
   rasa1.append(g1)
   rasa1.append(b1) 
   rgbpoints1=list(zip(*rasa1))
   return(rgbpoints1)

pazka

def process_rgb_raster():
inputfile1="./neo/BlueMarbleNG_2004-08-01_rgb_3600x1800.TIFF"
referencefile1="./work2/CHELSA_gts5_1993_V1.2.1.tif"
outputfile1="./work2/neo_rgb.tif"
#raster_reso_to_reference_raster(inputfile1, referencefile1, outputfile1)
referencefile2="./work4/CHELSA_gts5_1993_V1.2.1.tif"
outputfile2="./work4/neo_rgb.tif"
 #raster_reso_to_reference_raster(inputfile1, referencefile2, outputfile2)
 outputnc1=outputfile1.replace(".tif", ".nc")
outputnc2=outputfile2.replace(".tif", ".nc")
 tiff2nc(outputfile1, outputnc1)
tiff2nc(outputfile2, outputnc2)
 return(0);
def cut_and_scale_rgb_file(neofile1, neofile2,lon1, lat1, lon2, lat2, width1, height1, width2, height2):
 chelsadir2="./work1/"
chelsadir3="./work2/"
chelsadir4="./work3/"
chelsadir5="./work4/"
 inf1=neofile1
outf1=chelsadir2+neofile2
outf2=chelsadir3+neofile2
outf3=chelsadir4+neofile2
outf3=outf3.replace(".tif", ".nc")
outf4=chelsadir5+neofile2
outf5=outf4.replace(".tif", ".nc")
print(inf1)
clip_raster_to_outfile(inf1, outf1, lon1, lat1, lon2, lat2)
ds=0
ds = gdal.Warp(outf2, outf1, width=width1, height=height1, resampleAlg='bilinear')
ds=0
ds = gdal.Warp(outf4, outf1, width=width2, height=height2, resampleAlg='bilinear')
ds = 0
ds = gdal.Translate(outf3, outf2, format='NetCDF')
ds = gdal.Translate(outf5, outf4, format='NetCDF')
 ds=0
return(0)


def rgb_downscale_timeseries_trace(repo1, repo2, simu1, year1, year2, lon1, lat1, lon2, lat2, width1, height1):
print("RGB downscale test ...")
 slice1=int(year2/100)
 bigs1=[]
smalls1=[]
 variable1='tmean'
print("Timeseries")
#if (variable1=='tmean'):
for variable1 in ("tmean", "prec"):
print (variable1)
for n in range(1,13):
print(" TS variables ", n)
 month1=n
maxname1="./work2/"+chelsa_timeseries_name(variable1, year1,month1)
minname1="./work4/"+chelsa_timeseries_name(variable1, year1,month1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 #print (np.shape(maxarray1))
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls1.append(minna1)
bigs1.append(maxa1)


bigs2=[]
smalls2=[]
#if(variable1=='tmean'):
for variable1 in ("tmean", "pr"):
print (variable1)
for n in range(1,13):
print("n TRACE variables ", n)
 month1=n
maxname1="./work2/"+chelsa_trace_name(variable1, year2,month1)
minname1="./work4/"+chelsa_trace_name(variable1, year2,month1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 #print (np.shape(maxarray1))
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls2.append(minna1)
bigs2.append(maxa1)


#print("Target rasters ...")
month1=0
#mintargetname1="./work4/"+chelsa_timeseries_name("gts5", year1,month1)
#maxtargetname1="./work2/"+chelsa_timeseries_name("gts5", year1,month1)
 #mintargetname1="./work4/neo_ndvi.tif"
#maxtargetname1="./work2/neo_ndvi.tif"
 mintargetname1="./work4/neo_rgb.tif"
maxtargetname1="./work2/neo_rgb.tif"
 #loadrgb(filename)
#writergb(salbuname, filename,arrays)
 maxtarget1 = gdal.Open(maxtargetname1)
maxtargetarray1 = np.array(maxtarget1.GetRasterBand(1).ReadAsArray())
mintarget1 = gdal.Open(mintargetname1)
mintargetarray1 = np.array(mintarget1.GetRasterBand(1).ReadAsArray())
minnatarget1=mintargetarray1.flatten()
maxatarget1=maxtargetarray1.flatten()
 maxtarget2 = gdal.Open(maxtargetname1)
maxtargetarray2 = np.array(maxtarget1.GetRasterBand(2).ReadAsArray())
mintarget2 = gdal.Open(mintargetname1)
mintargetarray2 = np.array(mintarget2.GetRasterBand(2).ReadAsArray())
minnatarget2=mintargetarray2.flatten()
maxatarget2=maxtargetarray2.flatten()
 maxtarget3 = gdal.Open(maxtargetname1)
maxtargetarray3 = np.array(maxtarget3.GetRasterBand(3).ReadAsArray())
mintarget3 = gdal.Open(mintargetname1)
mintargetarray3 = np.array(mintarget3.GetRasterBand(3).ReadAsArray())
minnatarget3=mintargetarray3.flatten()
maxatarget3=maxtargetarray3.flatten()


apoints1=list(zip(*smalls1))
cpoints1=list(zip(*bigs1))
apoints2=list(zip(*smalls2))
cpoints2=list(zip(*bigs2))
 bpoints1=minnatarget1
bpoints2=minnatarget1
y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints2) ##ok
youtarray1=y2.reshape( bheight1, bwidth1).astype(float)
bpoints1=minnatarget2
bpoints2=minnatarget2
y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints2) ##ok
youtarray2=y2.reshape( bheight1, bwidth1).astype(float)
bpoints1=minnatarget3
bpoints2=minnatarget3
y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints2) ##ok
youtarray3=y2.reshape( bheight1, bwidth1).astype(float)
 arrays=[]
arrays.append(youtarray1)
arrays.append(youtarray2)
arrays.append(youtarray3)
 #image_size = (width1, height1)
image_size = (height1, width1)
 lat = [lat1,lat2]
lon = [lon1,lon2]




r_pixels = youtarray1.astype(np.uint8)
g_pixels = youtarray2.astype(np.uint8)
b_pixels = youtarray3.astype(np.uint8)
 nx = image_size[0]
ny = image_size[1]
xmin, ymin, xmax, ymax = [min(lon), min(lat), max(lon), max(lat)]
xres = (xmax - xmin) / float(nx)
yres = (ymax - ymin) / float(ny)
geotransform = (xmin, xres, 0, ymax, 0, -yres)
 dst_ds = gdal.GetDriverByName('GTiff').Create('./output1/rgb.tif', ny, nx, 3, gdal.GDT_Byte)
dst_ds.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(3857)
dst_ds.SetProjection(srs.ExportToWkt())
dst_ds.GetRasterBand(1).WriteArray(r_pixels)
dst_ds.GetRasterBand(2).WriteArray(g_pixels)
dst_ds.GetRasterBand(3).WriteArray(b_pixels)
dst_ds.FlushCache()
dst_ds = None


#writergb("./work2/neo_rgb.tif", "./output/rgbtest.tif",arrays)
 #array2nc("./output1/downscaled.nc",maxtargetname1,youtarray1)
#array2raster("./output1/downscaled.tiff",maxtargetname1,youtarray1)
 return(0)


def calcu_daylength(dayOfYear, lat):
latInRad = np.deg2rad(lat)
declinationOfEarth = 23.45*np.sin(np.deg2rad(360.0*(283.0+dayOfYear)/365.0))
if (-np.tan(latInRad) * np.tan(np.deg2rad(declinationOfEarth)) <= -1.0):
return (24.0)
elif (-np.tan(latInRad) * np.tan(np.deg2rad(declinationOfEarth)) >= 1.0):
return (0.0)
else:
hourAngle = np.rad2deg(np.arccos(-np.tan(latInRad) * np.tan(np.deg2rad(declinationOfEarth))))
 return (2.0*hourAngle/15.0)


def calcu_daylength_of_month(month, latgrida1):
dof1=month*30+15
nax=np.shape(latgrida1)[1]-1
nay=np.shape(latgrida1)[0]-1
print(nax,nay)
daylengthgrida1=latgrida1*0.0
for nny in range(0,nay):
for nnx in range(0,nax):
latt1=latgrida1[nny,nnx]
daylengthss1=calcu_daylength(dof1, latt1)
daylengthgrida1[nny,nnx]=daylengthss1
 return(daylengthgrida1)




def days_in_month(month,year=2013):
return calendar.monthrange(year,month)[1]


def radiation(year, month, day, lat, lon,minute_step=30):
dayrad=0.0
for hour in range(24):
for minute in range(0,60,minute_step):
thr = hour + minute/60.
second=0
d = datetime(year, month, day, hour, minute, second, tzinfo=timezone.utc)
altitude_deg = pysolar.solar.get_altitude(lat, lon, d)
solar_rad = pysolar.solar.radiation.get_radiation_direct(d, altitude_deg)
altitude_deg= 90. - altitude_deg
solar_rad2=solar_rad*60*minute_step
dayrad=dayrad+solar_rad2
 return (dayrad)
def monthly_solar_incoming_radiation(year,month,lat):
minute_step=30
days=calendar.monthrange(year,month)[1]+1
dayrad=abs(lat)*0.0
#print(dayrad)
 #quit(-1)
edayrada=dayrad
for n in range(1,days):
dayrada=radiation(year, month, n, lat, 0,minute_step=30)
#print(dayrada)
dayradasum=np.sum(dayrada)
if(np.isnan(dayradasum)):
dayrada=edayrada


dayrad=dayrad+dayrada
edayrada=dayrada
 return(dayrad)
def load_raster_float32(filename):
print(filename)
raster = gdal.Open(filename)
print(raster)
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
d1 = geotransform[2]
d2 = geotransform[4]
array = raster.ReadAsArray()
 cols = array.shape[1]
rows = array.shape[0]
 limitX=cols*pixelWidth+originX
limitY=rows*pixelHeight+originY
lons=np.array(np.arange(originX, limitX,pixelWidth))
lats=np.array(np.arange(originY, limitY,pixelHeight))


return(array, lats, lons)




def cut_and_scale_raster_file(neofile1, neofile2,lon1, lat1, lon2, lat2, width1, height1, width2, height2):
 chelsadir2="./work1/"
chelsadir3="./work2/"
chelsadir4="./work3/"
chelsadir5="./work4/"
 inf1=neofile1
outf1=chelsadir2+neofile2
outf2=chelsadir3+neofile2
outf3=chelsadir4+neofile2
outf3=outf3.replace(".tif", ".nc")
outf4=chelsadir5+neofile2
outf5=outf4.replace(".tif", ".nc")
print(inf1)
clip_raster_to_outfile(inf1, outf1, lon1, lat1, lon2, lat2)
ds=0
ds = gdal.Warp(outf2, outf1, width=width1, height=height1, resampleAlg='bilinear')
ds=0
ds = gdal.Warp(outf4, outf1, width=width2, height=height2, resampleAlg='bilinear')
ds = 0
ds = gdal.Translate(outf3, outf2, format='NetCDF')
ds = gdal.Translate(outf5, outf4, format='NetCDF')
 ds=0
return(0)
def getFeatures(gdf):
return [json.loads(gdf.to_json())['features'][0]['geometry']]
def clip_raster_to_outfile_1(inf1, outf1, lon1, lat1, lon2, lat2):
data = rasterio.open(inf1)
bbox=box(lon1, lat1, lon2, lat2)
geo=gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=from_epsg(4326))
#geo=geo.to_crs(crs=data.crs.data)
coords = getFeatures(geo)
#print(coords)
out_img, out_transform = mask(data, shapes=coords, crop=True)
out_meta = data.meta.copy()
#epsg_code = int(data.crs.data['init'][1:])
epsg_code=int(4326)
height=out_img.shape[1]
width=out_img.shape[2]
#print(height)
#print (width)
out_meta.update({"driver": "GTiff","height": out_img.shape[1],"width": out_img.shape[2],"transform": out_transform,"crs": pycrs.parse.from_epsg_code(epsg_code).to_proj4()})
 with rasterio.open(outf1, "w", **out_meta) as dest:
dest.write(out_img)
 clipped = rasterio.open(outf1)
return(0)


def create_landsea(neofile1, neofile2,lon1, lat1, lon2, lat2, cols,rows):
 array1, lats, lons=load_raster_float32(neofile1)
#array1=np.array(array1).astype(float)
#array2=np.maximum(array1, 0, array1)
array1[array1 < 0.0] = 0.0
array1[array1 > 0.0] = 1.0
array2=array1.astype(np.float32)
array3=np.flipud(array2)
 drivername1="GTiff"
saveraster_float32(neofile2,drivername1, array3, lon1, lon2, lat1, lat2)
drivername2="NetCDF"
neofile3=neofile2.replace(".tif", ".nc")
saveraster_float32(neofile3,drivername2, array3, lon1, lon2, lat1, lat2)
 return(0)
def saveraster_float32(newname,drivername, array, lon1, lon2, lat1, lat2):
#print (array)
print("Save")
#print (array.dtype)
#print(array)
cols=array.shape[1]
rows=array.shape[0]
 array2=array.astype(np.float32)
pixelWidth = (lon2-lon1)/cols
pixelHeight = (lat2-lat1)/rows


print(cols, rows)
print (array2.shape)


driver = gdal.GetDriverByName(drivername)
outRaster = driver.Create(newname, cols, rows, 1, gdal.GDT_Float32, options=['COMPRESS=DEFLATE'])
outRaster.SetGeoTransform((lon1, pixelWidth, 0, lat1, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array2)
spatialReference = osr.SpatialReference()
 #spatialReference.SetUTM(zonenum, zone >= 'N')
#spatialReference.SetWellKnownGeogCS('WGS84')
#
#retval = dataset.SetProjection(wkt)


outRasterSRS = osr.SpatialReference()
#outRasterSRS.ImportFromWkt("EPSG:4326")
spatialReference.SetWellKnownGeogCS('WGS84')
wkt = spatialReference.ExportToWkt()
 outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
driver=0
return(0)


def sea_proximity(srcname, dstname):
src_ds = gdal.Open(srcname)
srcband=src_ds.GetRasterBand(1)
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create(dstname,src_ds.RasterXSize, src_ds.RasterYSize, 1,gdal.GetDataTypeByName('Float32'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
gdal.ComputeProximity(srcband,dstband,["VALUES='0,-1'","DISTUNITS=GEO"])
srcband = None
dstband = None
src_ds = None
dst_ds = None
return(0)

def gts(meantemps, basetemp, tempscale = 10, tempofset=273.16):

def gts(meantemps, basetemp, tempscale = 1, tempofset=0):
#rem note temperatures can be kelvins*10 c, not celcius
 ndays =[31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
gdds=basetemp*0.0
for n in range(0,12):
print(n)
gddx=meantemps[n]
gddx=(gddx/tempscale)-tempofset
 gddx=gddx-basetemp
 gddx[gddx < 0.0] = 0.0
gddx=gddx * ndays[n]
gdds=gdds+gddx


gdds=np.flipud(gdds)
return(gdds)


def load_temp_stack_1(year2):
#namebase1=CHELSA_tmean_1993_10_V1.2.1
slice2=int(year2/100)
temps=[]
 for n in range(1,13):
print(n)
#name1="./work2/CHELSA_tmean_1993_"+str(n).zfill(2)+"_V1.2.1.tif"
#name1="./work2/CHELSAtrace_tmean_"+str(n)+"_-"+str(slice2)+"_V1.0.tif"
name1="./work2/CHELSA_TraCE21k_tmean_"+str(n)+"_-"+str(slice2)+"_V1.0.tif"
print(name1)
temp, lons, lats=load_raster_float32(name1)
temp2=(temp*0.1)-273.16
temps.append(temp2)
 return(temps, lons, lats)
def load_precip_stack_1(year2):
slice2=int(year2/100)
temps=[]
 for n in range(1,13):
print(n)
#name1="./work4/CHELSA_prec_1993_"+str(n).zfill(2)+"_V1.2.1.tif"
name1="./work2/CHELSA_TraCE21k_pr_"+str(n)+"_-"+str(slice2)+"_V1.0.tif"
print(name1)
temp, lons, lats=load_raster_float32(name1)
#temp2=(temp*0.1)-273.16
temps.append(temp)
 return(temps, lons, lats)


def load_evap_stack_1():
 temps=[]
 for n in range(1,13):
print(n)
name1="./work2/th_evapo1_"+str(n)+".nc"
#print(name1)
temp, lons, lats=load_raster_float32(name1)
#temp2=(temp*0.1)-273.16
temps.append(temp)
 return(temps, lons, lats)




def koe_evaporation_thorntwaite(temperasters, latraster, lonraster):

maybe nok

 dimx, dimy=np.shape(latraster)
 tmarray=np.array(temperasters)
 latdegreesarray=latraster.flatten()
londegreesarray=lonraster.flatten()
 latradiansarray=latdegreesarray*0.0
 len1=len(latdegreesarray)
print(len1)
 for n in range(0,len1):
latdegree=latdegreesarray[n]
latradiansarray[n]=pyeto.deg2rad(latdegree)


print("tmarray")
print (np.shape(tmarray))
 tempix0=[]
 for iy in range(0,dimy):
for ix in range(0,dimx):
tamo=[]
for mo in range(0,12):
tt=tmarray[mo,iy,ix]
tamo.append(tt)
 tamos=np.array(tamo).astype(float)
tempix0.append(tamos)


tempix1=np.array(tempix0).astype(float)
tempix=(tempix1/10.0)-273.16
print(np.shape(tempix1))




 #quit(-1)
 #print (latradiansarray)
#quit(-1)
 #daylightarray = pyeto.monthly_mean_daylight_hours(latradiansarray)
len1=len(latradiansarray)
print(len1)
 #daylightarray=latradiansarray*0
daylights0=[]
 for n in range(0,len1):
daylighthours=pyeto.monthly_mean_daylight_hours(latradiansarray[n])
daylights0.append(np.array(daylighthours).astype(float))
 daylights=np.array(daylights0).astype(float)
 evapos0=[]
 for n in range(0,len1):
tempo=tempix[n]
daylo=daylights[n]
evapo1=pyeto.thornthwaite(tempo, daylo)
evapos0.append(np.array(evapo1).astype(float))
 evapos=np.array(evapos0).T
 print("Shape evapos")
print (np.shape(evapos))
 #quit(-1)
monthevapos=[]
 for m in range(0,12):
#evapotab0=evapos[0:][m]
evapotab0=evapos[m]
print (len(evapotab0))
evapotab1=evapotab0.reshape(dimy, dimx)
monthevapos.append(evapotab1)


#plt.imshow(monthevapos[7])
#plt.show()
 #quit(-1)
 pointer1=62000
temp1=tempix[pointer1]
lat1=latdegreesarray[pointer1]
lon1=londegreesarray[pointer1]
evapotabu1=evapos[6]
nayte_evapo1=evapotabu1[pointer1]
 print ("Nayte")
print(lat1, lon1)
print (temp1)
print(nayte_evapo1)
latitude_radians1 = pyeto.deg2rad(lat1)
mean_monthly_temperature = temp1
monthly_mean_daylight = pyeto.monthly_mean_daylight_hours(latitude_radians1)
evapo1=pyeto.thornthwaite(mean_monthly_temperature, monthly_mean_daylight)
print("Algo")
print(evapo1)


 return(monthevapos)


def calcu_evaporation_thornthwaite_2(year2, lon1,lat1, lon2, lat2):
## maybe ok
ndays =[31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
#def calcu_daylength(dayOfYear, lat):
 temps, lats, lons=load_temp_stack_1(year2)
longrida, latgrida = np.meshgrid(lons, lats)


termal_index1=temps[0]*0.0
 for n in range(-1,12):
print(n+1)
tempi=temps[n]
temp5=tempi*0.2
indexmonth1=np.power(temp5, 1.514)
where_are_NaNs = np.isnan(indexmonth1)
indexmonth1[where_are_NaNs] = 0.0001
#indexmonth1[indexmonth1>0] = 0
termal_index1=termal_index1+indexmonth1
 #plt.imshow(termal_index1)
#plt.show()
 #quit(-1)


## coarse value!
#alphacoef1=0.49239+termal_index1*0.01792
alphacoef1=0.49239+termal_index1*0.01792
 ak2=np.power(termal_index1,2.0)*7.71e-5
where_are_NaNs = np.isnan(ak2)
ak2[where_are_NaNs] = 0.0
ak3=6.72e-7*np.power(termal_index1,3.0)
where_are_NaNs = np.isnan(ak3)
ak3[where_are_NaNs] = 0.0
 alphacoef1=alphacoef1+ak2+ak3


 #plt.imshow(alphacoef1)
#plt.show()


solar_evapo=[]


for n in range(-1,12):
print(n)
#inname2="./work2/solar_"+str(n+1)+".tif"
#solar1,latas, lonas=load_raster_float32(inname2)
#evapo_solaronly= solar1* 0.0864 * 0.408/1000000*ndays[n+1]
 daylengthgrida1=calcu_daylength_of_month(n, latgrida)
#print(daylengthgrida1)
alla12=daylengthgrida1/12.0
anna30=ndays[n]/30.0
anna30=1
tempa=temps[n]
tempa[tempa < 0.0] = 0.0


#print(tempa)
terma10a=(10*tempa)/termal_index1
terma10b=np.power(terma10a, alphacoef1)
#pet_th=16*alla12*anna30*terma10b
##pet_nc=16.0*terma10b
pet_nc=16*alla12*anna30*terma10b
#plt.imshow(pet_nc)
#plt.show()
#quit(-1)
 evapoo1=np.array(pet_nc).astype(float)
 outname1="./work2/th_evapo1_"+str(n+1)+".nc"
evapoo1=np.flipud(evapoo1)
saveraster_float32(outname1,"NetCDF", evapoo1, lon1, lon2, lat1, lat2)
 return(solar_evapo)
def river_proximity(srcname, dstname):
src_ds = gdal.Open(srcname)
srcband=src_ds.GetRasterBand(1)
## jn waruning test code
#srcarray=src_ds.ReadAsArray()
#srcarray=srcarray*0
#srcband.WriteArray(srcarray)
 drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create(dstname,src_ds.RasterXSize, src_ds.RasterYSize, 1,gdal.GetDataTypeByName('Float32'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)


gdal.ComputeProximity(srcband,dstband,["VALUES='0,-1'","DISTUNITS=GEO"])
srcband = None
dstband = None
src_ds = None
dst_ds = None
return(0)


def downscale_timeseries_trace_neo(outname1,targetname1, simu1, year1, year2):
print("Timeseries vs trace downscale test ...")
 slice1=int(year2/100)
 bigs1=[]
smalls1=[]
 variable1='tmean'
print("Timeseries")
#if (variable1=='tmean'):
for variable1 in ("tmean", "prec"):
#for variable1 in ("tasmax", "tasmin"):
print (variable1)
for n in range(1,13):
#for n in range(4,10):
print(" TS variables ", n)
 month1=n
maxname1="./work2/"+chelsa_timeseries_name(variable1, year1,month1)
minname1="./work4/"+chelsa_timeseries_name(variable1, year1,month1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 #print (np.shape(maxarray1))
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls1.append(minna1)
bigs1.append(maxa1)


bigs2=[]
smalls2=[]
#if(variable1=='tmean'):
for variable1 in ("tmean", "pr"):
#for variable1 in ("tasmax", "tasmin"):
print (variable1)
for n in range(1,13):
#for n in range(4,10):
print("n TRACE variables ", n)
 month1=n
maxname1="./work2/"+chelsa_trace_name(variable1, year2,month1)
minname1="./work4/"+chelsa_trace_name(variable1, year2,month1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 #print (np.shape(maxarray1))
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls2.append(minna1)
bigs2.append(maxa1)


 #print("Target rasters ...")
month1=0
#mintargetname1="./work4/"+chelsa_timeseries_name("gts5", year1,month1)
#maxtargetname1="./work2/"+chelsa_timeseries_name("gts5", year1,month1)
 mintargetname1="./work4/"+targetname1
maxtargetname1="./work2/"+targetname1
 #loadrgb(filename)
#writergb(salbuname, filename,arrays)
 print(maxtargetname1)
print(mintargetname1)
 maxtarget1 = gdal.Open(maxtargetname1)
maxtargetarray1 = np.array(maxtarget1.GetRasterBand(1).ReadAsArray())
mintarget1 = gdal.Open(mintargetname1)
mintargetarray1 = np.array(mintarget1.GetRasterBand(1).ReadAsArray())
 minnatarget1=mintargetarray1.flatten()
maxatarget1=maxtargetarray1.flatten()
 ## note warning this is data filter, very specific to dataset
## clamps to max, min values
 flatten_to=1
flatten_upper_limit=1.0
flatten_lower_limit=0.0
if(flatten_to==1):
minnatarget1 = np.where(minnatarget1 < flatten_lower_limit, flatten_lower_limit, minnatarget1)
maxatarget1 = np.where(maxatarget1 < flatten_lower_limit, flatten_lower_limit, maxatarget1)
minnatarget1 = np.where(minnatarget1 > flatten_upper_limit, flatten_upper_limit, minnatarget1)
maxatarget1 = np.where(maxatarget1 > flatten_upper_limit, flatten_upper_limit, maxatarget1)


apoints1=list(zip(*smalls1))
bpoints1=minnatarget1
cpoints1=list(zip(*bigs1))
 apoints2=list(zip(*smalls2))
bpoints2=minnatarget1
cpoints2=list(zip(*bigs2))
 print("Downscaling ....")
 ## test only
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints1)
##y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints1) ## nogood
#y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints1) ##ok


#y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints2) ##NOK
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints2)
y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints2) ##ok


youtarray1=y2.reshape( bheight1, bwidth1).astype(float)


 array2nc(outname1,maxtargetname1,youtarray1)
#array2raster("./output1/downscaled.tiff",maxtargetname1,youtarray1)
 return(0)




def raster_reso_to_reference_raster(inputfile, referencefile, outputfile):
print(inputfile, referencefile, outputfile)
input = gdal.Open(inputfile, gdalconst.GA_ReadOnly)
inputProj = input.GetProjection()
inputTrans = input.GetGeoTransform()
reference = gdal.Open(referencefile)
referenceProj = reference.GetProjection()
referenceTrans = reference.GetGeoTransform()
bandreference = reference.GetRasterBand(1)
x = reference.RasterXSize
y = reference.RasterYSize
driver= gdal.GetDriverByName('GTiff')
output = driver.Create(outputfile,x,y,1,bandreference.DataType)
output.SetGeoTransform(referenceTrans)
output.SetProjection(referenceProj)
gdal.ReprojectImage(input,output,inputProj,referenceProj,gdalconst.GRA_Bilinear)
del output
return(0)




def downscale_own_dirs(outname1,maxtargetname1, mintargetname1, dira1, dira2, dirb1, dirb2):
print(" File against own files test ...")
 bigs1=[]
smalls1=[]
 for ff in os.listdir(dira1):
maxname1=dira1+"/"+ff
minname1=dira2+"/"+ff
print(minname1)
 max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
 minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
 bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls1.append(minna1)
bigs1.append(maxa1)


bigs2=[]
smalls2=[]
 print("")
print("KKKKKKK")
for ff in os.listdir(dirb1):
print(ff)
maxname1=dirb1+"/"+ff
minname1=dirb2+"/"+ff
print(minname1)
max1 = gdal.Open(maxname1)
maxarray1 = np.array(max1.GetRasterBand(1).ReadAsArray())
min1 = gdal.Open(minname1)
minarray1 = np.array(min1.GetRasterBand(1).ReadAsArray())
minna1=minarray1.flatten()
maxa1=maxarray1.flatten()
bheight1=np.shape(maxarray1)[0]
bwidth1=np.shape(maxarray1)[1]
bheight2=np.shape(minarray1)[0]
bwidth2=np.shape(minarray1)[1]
 smalls2.append(minna1)
bigs2.append(maxa1)


 #print("Target rasters ...")
month1=0
#mintargetname1="./work4/"+chelsa_timeseries_name("gts5", year1,month1)
#maxtargetname1="./work2/"+chelsa_timeseries_name("gts5", year1,month1)


#loadrgb(filename)
#writergb(salbuname, filename,arrays)
 maxtarget1 = gdal.Open(maxtargetname1)
maxtargetarray1 = np.array(maxtarget1.GetRasterBand(1).ReadAsArray())
mintarget1 = gdal.Open(mintargetname1)
mintargetarray1 = np.array(mintarget1.GetRasterBand(1).ReadAsArray())
 minnatarget1=mintargetarray1.flatten()
maxatarget1=maxtargetarray1.flatten()
 ## note warning this is data filter, very specific to dataset
## clamps to max, min values
 flatten_to_01=1
 if(flatten_to_01==1):
minnatarget1 = np.where(minnatarget1 < 0, 0, minnatarget1)
maxatarget1 = np.where(maxatarget1 < 0, 0, maxatarget1)
minnatarget1 = np.where(minnatarget1 > 1, 0, minnatarget1)
maxatarget1 = np.where(maxatarget1 > 1, 0, maxatarget1)
 apoints1=list(zip(*smalls1))
bpoints1=minnatarget1
cpoints1=list(zip(*bigs1))
 apoints2=list(zip(*smalls2))
bpoints2=minnatarget1
cpoints2=list(zip(*bigs2))
 print("Downscaling ....")
 ## test only
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints1)
##y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints1) ## nogood
#y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints1) ##ok


#y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints2) ##NOK
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints2)
y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints2) ##ok


youtarray1=y2.reshape( rows, cols).astype(float)


 save_points_to_netcdf(outfilename1, "var", lons1, lats1, youtarray1)
#array2nc(outname1,maxtargetname1,youtarray1)
#array2raster("./output1/downscaled.tiff",maxtargetname1,youtarray1)
 return(0)


def process_my_own_rasters(rasdir1, rasdir2, referenceneo1):
for ff in os.listdir(rasdir1):
fina1=rasdir1+"/"+ff
fina2=rasdir2+"/"+ff
print (fina1)
print (fina2)
print (referenceneo1)
raster_reso_to_reference_raster(fina1, referenceneo1, fina2)
 return(0)


def gdaldem_test1(finktio1, outformat1, infile1, outfile1, par1, par2, par3=0):
 k1="gdaldem "
k2=finktio1+" "
k3="-of "+outformat1+" "
 kommando1=k1+k2
 if(finktio1=="hillshade"):
p1="-az "+str(par1)+ " "
p2="-alt "+str(par2)+" "


kommando1=kommando1+p1+p2+infile1+" "+outfile1
print(kommando1)
 #return(0)
 os.system(kommando1)
 return(0)
def get_chelsa_pmip3_data(simu1):
 print("Get Chalsa PMIP3 data ", simu1)


for n in range(1,13):
get_chelsa_variable("pmip3",simu1, "prec",0, n)
 for n in range(1,13):
get_chelsa_variable("pmip3",simu1, "tmean", 0, n)
 for n in range(1,20):
get_chelsa_variable("pmip3",simu1, "bio", 0, n)
 get_chelsa_variable("pmip3","DEM", 0,0, 0)


return(0);
def get_chelsa_climatologies_data():
 print("Get Chelsa climatologies data ")
 for n in range(1,13):
get_chelsa_variable("climatologies","", "tmean", 0, n)
 for n in range(1,20):
get_chelsa_variable("climatologies","", "bio", 0, n)
 for n in range(1,13):
get_chelsa_variable("climatologies","", "prec",0, n)
#https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/climatologies/prec/CHELSA_prec_01_V1.2_land.tif


 return(0);


def get_chelsa_cruts_data(year1):
 for n in range(1,13):
get_chelsa_variable("chelsa_cruts","", "prec", year1, n)


for n in range(1,13):
get_chelsa_variable("chelsa_cruts","", "tmax", year1, n)
 for n in range(1,13):
get_chelsa_variable("chelsa_cruts","", "tmin", year1, n)


for n in range(1,20):
get_chelsa_variable("chelsa_cruts","", "bio", year1, n)
 return(0)


def cut_chelsa_dir(chelsadira1, lon1, lat1, lon2, lat2, width1, height1, width2, height2):
global set_biovars_off
 cutdir1="./work2/"
cutdir2="./work4/"
chelsadir2="./work1/"
chelsadir3="./work2/"
chelsadir4="./work3/"
chelsadir5="./work4/"
 for ff in os.listdir(chelsadira1):
bobo=0
 if(set_biovars_off==1):
if 'bio' in ff:
bobo=1
 if(set_biovars_off==1):
if 'BIO' in ff:
bobo=1


if(bobo==0):
fina1=chelsadira1+"/"+ff
tana1=cutdir1+"/"+ff
tana1=cutdir1+"/"+ff
inf1=fina1
outf1=chelsadir2+ff
outf2=chelsadir3+ff
outf3=chelsadir4+ff
outf3=outf3.replace(".tif", ".nc")
outf4=chelsadir5+ff
outf5=outf4.replace(".tif", ".nc")
print(inf1)
clip_raster_to_outfile(inf1, outf1, lon1, lat1, lon2, lat2)
ds=0
ds = gdal.Warp(outf2, outf1, width=width1, height=height1, resampleAlg='bilinear')
ds=0
ds = gdal.Warp(outf4, outf1, width=width2, height=height2, resampleAlg='bilinear')
ds = 0
ds = gdal.Translate(outf3, outf2, format='NetCDF')
ds = gdal.Translate(outf5, outf4, format='NetCDF')
ds=0
 return(0)




def listauz(dire1, patterni1):
listaa1 = glob.glob(dire1+"/"+patterni1)
return(listaa1)


def loadraster_stak(direktory1, pattern1, pattern2):
## raster stack
rasta=[]
for ff in os.listdir(direktory1):
fina1=direktory1+"/"+ff
if pattern2 in ff:
if pattern1 in ff:
print(ff)
ras, lonk, latt=load_raster_float32(fina1)
#ras2=np.array(ras).flatten().astype(float)
rasta.append(ras)
 return(rasta, lonk, latt)
def loadraster_stak_arrays(direktory1, pattern1, pattern2):
## stack of raster 1d array for downscaling
rasta=[]
for ff in os.listdir(direktory1):
fina1=direktory1+"/"+ff
if pattern2 in ff:
if pattern1 in ff:
print(ff)
ras, lonk, latt=load_raster_float32(fina1)
ras2=np.array(ras).flatten().astype(float)
rasta.append(ras2)
 return(rasta, lonk, latt)


def downscale_vars_directly(inname1, inname2, outname1, varname1, varname2,lon1, lon2, lat1, lat2):
workdir1="./work2"
workdir2="./work4"
 bigs1, lo, la=loadraster_stak_arrays(workdir1, varname1,".tif")
bigs2, lo, la=loadraster_stak_arrays(workdir1, varname2, ".tif")
smalls1, lo2, la2=loadraster_stak_arrays(workdir2, varname1,".tif" )
smalls2, lo2, la2=loadraster_stak_arrays(workdir2, varname2,".tif" )
 #warinh kelvin to celcius
for n in range(0,12):
smalls1[n]=smalls2[n]*0.1-273.16
bigs2[n]=bigs2[n]*0.1-273.16


 print (len(bigs1))
print (len(bigs2))
print (len(smalls1))
print (len(smalls2))
 print(len(bigs1[1]))
print(len(bigs2[1]))
print(len(smalls1[1]))
print(len(smalls2[1]))
#quit(-1)
#q1=bigs1[0]
#print(q1.shape())


maxtarget1 = gdal.Open(inname1)
 maxtargetarray1 = np.array(maxtarget1.GetRasterBand(1).ReadAsArray()).astype(float)
cols = int(maxtargetarray1.shape[1])
rows = int(maxtargetarray1.shape[0])
 mintarget1 = gdal.Open(inname2)
mintargetarray1 = np.array(mintarget1.GetRasterBand(1).ReadAsArray()).astype(float)
minnatarget1=mintargetarray1.flatten()
maxatarget1=maxtargetarray1.flatten()




#plt.imshow(maxtargetarray1)
 #plt.show()
 print(cols, rows)
 #quit(-1)


#print (len(maxtargetarray1) )
#print(len(mintargetarray1) )
print(len(maxatarget1) )
print(len(minnatarget1) )
 #print(varname1, varname2)
#print(inname1)
#print(inname2)


 ## note warning this is data filter, very specific to dataset
## clamps to max, min values
 flatten_to=1
flatten_upper_limit=100.0
flatten_lower_limit=0.0
if(flatten_to==1):
minnatarget1 = np.where(minnatarget1 < flatten_lower_limit, flatten_lower_limit, minnatarget1)
maxatarget1 = np.where(maxatarget1 < flatten_lower_limit, flatten_lower_limit, maxatarget1)
#minnatarget1 = np.where(minnatarget1 > flatten_upper_limit, flatten_upper_limit, minnatarget1)
#maxatarget1 = np.where(maxatarget1 > flatten_upper_limit, flatten_upper_limit, maxatarget1)
 #midas1=minnatarget1.reshape(100,300)
#plt.imshow(midas1)
 #plt.show()
 #quit(-1)




apoints1=list(zip(*smalls1))
bpoints1=minnatarget1
cpoints1=list(zip(*bigs1))
 apoints2=list(zip(*smalls2))
bpoints2=minnatarget1
cpoints2=list(zip(*bigs2))
 print("Downscaling ....")
 ## test only
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints1)
##y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints1) ## nogood
#y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints1) ##ok


#y2=linear_regression_multiple_vars(apoints1, bpoints1, cpoints2) ##NOK
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints2)
#y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints2) ##ok
 #print(np.shape(apoints1[0]))
#print(np.shape(bpoints1))
#print(np.shape(cpoints1[0]))
print(bpoints1)




y2=random_forest_multiple_vars(apoints1, bpoints1, cpoints2) ##ok
#y2=gam_multiple_vars(apoints1, bpoints1, cpoints2)


 yo1=y2.reshape( rows, cols)
yo2=np.flipud(yo1)


print(yo1)
#plt.imshow(yo1)
 #plt.show()
 #youtarray2=youtarray1[1]
 saveraster_float32(outname1,"NetCDF", yo2, lon1, lon2, lat1, lat2)
#save_points_to_netcdf(outname1, "var", lo, la, y2)
return(0)


def calculate_tmean_from_chelsa_cruts_data(repo1, simu1, year1, lon1, lat1, lon2, lat2, width1, height1, width2, height2):

nok

 #for n in range(1,13):
for n in range(1,2):
print(n)
namea1="./work2/"+chelsa_cruts_name("tmax", year1, n)
namea2="./work2/"+chelsa_cruts_name("tmin", year1, n)
namea3="./work2/"+chelsa_cruts_name("tmean", year1, n)
namea3a="./work3/"+chelsa_cruts_name("tmean", year1, n)
nameb1="./work4/"+chelsa_cruts_name("tmax", year1, n)
nameb2="./work4/"+chelsa_cruts_name("tmin", year1, n)
nameb3="./work4/"+chelsa_cruts_name("tmean", year1, n)
namea4=namea3
namea4=namea4.replace("tif", "nc")
namea4=namea4.replace("work2", "work3")
nameb4=nameb3
nameb4=nameb4.replace("tif", "nc")
print(namea1)
print(nameb1)


array1, la, lo=load_raster_float32(namea1)
array2, la, lo=load_raster_float32(namea2)
print (array1)
 array3=(array1+array2)/2.0
array4=np.array(array3).astype(float)
 #plt.imshow(array4)
 #plt.show()
save_raster(namea3, "GTiff", array4,lon1,lat1,lon2,lat2, width1, height1)
print("Namea4")
print(namea4)
#save_points_to_netcdf(namea4, "var", lo, la, array4)
#nc_save_test_1("out.nc", array4, "Temp", la, lo)
nc_write(namea4, array4, lo, la)
#save_raster(namea4, "NetCDF", array4,lon1,lat1,lon2,lat2, width1, height1)
 brray1, la2, lo2=load_raster_float32(nameb1)
brray2, la2, lo2=load_raster_float32(nameb2)
 brray3=(brray1+brray2)/2.0
brray4=np.array(brray3).astype(float)
 save_raster(nameb3, "GTiff", brray4,lon1,lat1,lon2,lat2, width2, height2)
#save_raster(nameb4, "NetCDF", brray4,lon1,lat1,lon2,lat2, width2, height2)


 return(0)


def nc_write(fn, valus, lonis, latis):
ds = netCDF4.Dataset(fn, 'w', format='NETCDF4')
lat = ds.createDimension('lat', len(latis))
lon = ds.createDimension('lon', len(lonis))
lats = ds.createVariable('lat', 'f4', ('lat',))
lons = ds.createVariable('lon', 'f4', ('lon',))
value = ds.createVariable('value', 'f4', ('lat', 'lon',))
value.units = 'Unknown'
lats[:] = latis
lons[:] = lonis
value[:, :] = valus
 ds.title='Variaapeli'
#ds.subtitle='Hauska NC juttu'
#ds.Subtitle='Sopo 1'
#ds.koe='Raikaa'
#ds.setncattr_string('Subitile', ['a', 'b'])
#ds.setncattr_string('Subtitle', ['Tittuli'])
ds.setncattr_string('Subtitle', 'Tittuli')
value.units = 'Unknown'
value.long_name='Var1'
value.subtitle='Variable 1'
 print('var size after adding second data', value.shape)
ds.close()
return(0)
def rotate_rasters_180_to_360(inadir1, outadir1):
 for ina1 in os.listdir(inadir1):
bobo=0
 if(set_biovars_off==1):
if 'bio' in ina1:
bobo=1
 if(set_biovars_off==1):
if 'BIO' in ina1:
bobo=1
 if (bobo==0):
iname1=inadir1+ina1
outname1=outadir1+ina1
r1,lats1, lons1=load_raster_float32(iname1)
r2, lons2 = shiftgrid(0.0, r1, lons1, start=True)
rastype1="GTiff"
cols=len(lons2)
rows=len(lats1)
lon1=min(lons2)
lon2=max(lons2)
lat1=min(lats1)
lat2=max(lats1)
save_raster_compress(outname1, rastype1, r2,lon1,lat1,lon2,lat2,cols, rows)


return(0)
def climlab_calculate_insolation_milankovic(year2,month1, lat1, lat2, cols, rows):
 insolation_monthly_0=[]
 yr1=year2/1000.0*-1
 #years1 = np.linspace(yr1,yr1,1)
latsat1=np.linspace(lat1, lat2, rows)
#daisat1=np.linspace(1.0, 365.0, 365)
 day1=(month1-1)*30+15
 orb = OrbitalTable.interp(kyear=yr1)
 insol1 = daily_insolation(lat=latsat1, day=day1, orb=orb)
daysinmonth1=calendar.monthrange(1992,month1)[1]+1
insol2=insol1*daysinmonth1*3600*24
 colu1=np.array(insol2)
rowsu1=np.repeat(colu1, cols)
rowsu2=np.reshape(rowsu1, (cols,rows))
rowsu3=rowsu2
#plt.imshow(rowsu2)
#plt.show()
return(rowsu3)
def calcu_evaporation_thornthwaite_simu(year2, lon1,lat1, lon2, lat2):
## maybe ok
ndays =[31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
#def calcu_daylength(dayOfYear, lat):
 #temps, lats, lons=load_temp_stack_1(year2)
 slice2=int(year2/100)
temps=[]
 for n in range(1,13):
print(n)
name1="./work2/CHELSA_PMIP_CCSM4_tmean_"+str(n)+"_1.tif"
temp, lats, lons=load_raster_float32(name1)
temp2=(temp*0.1)-273.16
temps.append(temp2)
 longrida, latgrida = np.meshgrid(lons, lats)


termal_index1=temps[0]*0.0
 for n in range(-1,12):
print(n+1)
tempi=temps[n]
temp5=tempi*0.2
indexmonth1=np.power(temp5, 1.514)
where_are_NaNs = np.isnan(indexmonth1)
indexmonth1[where_are_NaNs] = 0.0001
#indexmonth1[indexmonth1>0] = 0
termal_index1=termal_index1+indexmonth1
 #plt.imshow(termal_index1)
#plt.show()
 #quit(-1)


## coarse value!
#alphacoef1=0.49239+termal_index1*0.01792
alphacoef1=0.49239+termal_index1*0.01792
 ak2=np.power(termal_index1,2.0)*7.71e-5
where_are_NaNs = np.isnan(ak2)
ak2[where_are_NaNs] = 0.0
ak3=6.72e-7*np.power(termal_index1,3.0)
where_are_NaNs = np.isnan(ak3)
ak3[where_are_NaNs] = 0.0
 alphacoef1=alphacoef1+ak2+ak3


 #plt.imshow(alphacoef1)
#plt.show()


solar_evapo=[]


for n in range(-1,12):
print(n)
#inname2="./work2/solar_"+str(n+1)+".tif"
#solar1,latas, lonas=load_raster_float32(inname2)
#evapo_solaronly= solar1* 0.0864 * 0.408/1000000*ndays[n+1]
 daylengthgrida1=calcu_daylength_of_month(n, latgrida)
#print(daylengthgrida1)
alla12=daylengthgrida1/12.0
anna30=ndays[n]/30.0
anna30=1
tempa=temps[n]
tempa[tempa < 0.0] = 0.0


#print(tempa)
terma10a=(10*tempa)/termal_index1
terma10b=np.power(terma10a, alphacoef1)
 print(np.shape(terma10b))
print(np.shape(alla12))
print(np.shape(anna30))
 #pet_th=16*alla12*anna30*terma10b
##pet_nc=16.0*terma10b
pet_nc=16*alla12*anna30*terma10b
#plt.imshow(pet_nc)
#plt.show()
#quit(-1)
 evapoo1=np.array(pet_nc).astype(float)
 outname1="./work2/th_evapo1_"+str(n+1)+".nc"
evapoo1=np.flipud(evapoo1)
saveraster_float32(outname1,"NetCDF", evapoo1, lon1, lon2, lat1, lat2)
 return(solar_evapo)


def change_data_values_scale(srcname1, dstname1, ofset1, coef1, lon1,lat1,lon2,lat2, cols, rows):
 array1, lats, lons=load_raster_float32(srcname1)
array2=array1*coef1+ofset1
dstname2=dstname1
dstname2=dstname2.replace('tif','nc')
save_raster(dstname1, "GTiff", array2,lon1,lat1,lon2,lat2,cols, rows)
save_raster(dstname2, "NetCDF", array2,lon1,lat1,lon2,lat2,cols, rows)
#saveraster_float32(dstname1,"GTiff", array2, lon1, lon2, lat1, lat2)
 return(0)
def generate_celcius_rasters(dataset1, dataset2,simu1,year1, year2, lon1,lat1,lon2,lat2, width1, height1, width2, height2):
srcpath1="./work2"
srcpath2="./work4"
 if(dataset2=="pmip"):
dstpath1="./work7"
dstpath2="./work8"
for n in range(1,13):
tempname1=chelsa_simu_name(simu1,"tmean", year1,n)
precname1=chelsa_simu_name(simu1,"prec", year1,n)
srctempname1=srcpath1+"/"+tempname1
dsttempname1=dstpath1+"/"+tempname1
srctempname2=srcpath2+"/"+tempname1
dsttempname2=dstpath2+"/"+tempname1
srcprecname1=srcpath1+"/"+precname1
dstprecname1=dstpath1+"/"+precname1
srcprecname2=srcpath2+"/"+precname1
dstprecname2=dstpath2+"/"+precname1
print(tempname1)
print(srctempname1)
ofset1=-273.16
coef1=0.1
change_data_values_scale(srctempname1, dsttempname1,ofset1, coef1, lon1,lat1,lon2,lat2,width1, height1)
change_data_values_scale(srctempname2, dsttempname2,ofset1, coef1, lon1,lat1,lon2,lat2,width2, height2)
ofset1=0
coef1=0.1
change_data_values_scale(srcprecname1, dstprecname1,ofset1, coef1, lon1,lat1,lon2,lat2, width1, height1)
change_data_values_scale(srcprecname2, dstprecname2,ofset1, coef1, lon1,lat1,lon2,lat2, width2, height2)
#print("KKKK")
#exit(-1)
 if(dataset1=="climatologies"):
dstpath1="./work5"
dstpath2="./work6"
for n in range(1,13):
tempname1=chelsa_climatologies_name("tmean", year1,n)
precname1=chelsa_climatologies_name("prec", year1,n)
srctempname1=srcpath1+"/"+tempname1
dsttempname1=dstpath1+"/"+tempname1
srctempname2=srcpath2+"/"+tempname1
dsttempname2=dstpath2+"/"+tempname1
srcprecname1=srcpath1+"/"+precname1
dstprecname1=dstpath1+"/"+precname1
srcprecname2=srcpath2+"/"+precname1
dstprecname2=dstpath2+"/"+precname1
print(tempname1)
print(srctempname1)
ofset1=0.0
coef1=0.1
change_data_values_scale(srctempname1, dsttempname1,ofset1, coef1, lon1,lat1,lon2,lat2,width1, height1)
change_data_values_scale(srctempname2, dsttempname2,ofset1, coef1, lon1,lat1,lon2,lat2, width2, height2)
ofset1=0
coef1=1.0
change_data_values_scale(srcprecname1, dstprecname1,ofset1, coef1, lon1,lat1,lon2,lat2, width1, height1)
change_data_values_scale(srcprecname2, dstprecname2,ofset1, coef1, lon1,lat1,lon2,lat2, width2, height2)
 ######################
###############################
 if(dataset2=="trace"):
 calculate_tmean_from_trace_tasmax_tasmin("trace", "", year2, lon1, lat1, lon2, lat2, width1, height1, width2, height2)


dstpath1="./work7"
dstpath2="./work8"
for n in range(1,13):
#def chelsa_trace_name(variable1, year1,month1):
tempname1=chelsa_trace_name("tmean", year2,n)
precname1=chelsa_trace_name("pr", year2,n)
srctempname1=srcpath1+"/"+tempname1
dsttempname1=dstpath1+"/"+tempname1
srctempname2=srcpath2+"/"+tempname1
dsttempname2=dstpath2+"/"+tempname1
srcprecname1=srcpath1+"/"+precname1
dstprecname1=dstpath1+"/"+precname1
srcprecname2=srcpath2+"/"+precname1
dstprecname2=dstpath2+"/"+precname1
print(tempname1)
print(srctempname1)
ofset1=-273.16
coef1=0.1
change_data_values_scale(srctempname1, dsttempname1,ofset1, coef1, lon1,lat1,lon2,lat2,width1, height1)
change_data_values_scale(srctempname2, dsttempname2,ofset1, coef1, lon1,lat1,lon2,lat2,width2, height2)
ofset1=0
coef1=1.0
change_data_values_scale(srcprecname1, dstprecname1,ofset1, coef1, lon1,lat1,lon2,lat2, width1, height1)
change_data_values_scale(srcprecname2, dstprecname2,ofset1, coef1, lon1,lat1,lon2,lat2, width2, height2)
#print("KKKK")
#exit(-1)
 if(dataset1=="timeseries"):
dstpath1="./work5"
dstpath2="./work6"
for n in range(1,13):
tempname1=chelsa_timeseries_name("tmean", year1,n)
precname1=chelsa_timeseries_name("prec", year1,n)
srctempname1=srcpath1+"/"+tempname1
dsttempname1=dstpath1+"/"+tempname1
srctempname2=srcpath2+"/"+tempname1
dsttempname2=dstpath2+"/"+tempname1
srcprecname1=srcpath1+"/"+precname1
dstprecname1=dstpath1+"/"+precname1
srcprecname2=srcpath2+"/"+precname1
dstprecname2=dstpath2+"/"+precname1
print(tempname1)
print(srctempname1)
ofset1=-273.16
coef1=0.1
change_data_values_scale(srctempname1, dsttempname1,ofset1, coef1, lon1,lat1,lon2,lat2,width1, height1)
change_data_values_scale(srctempname2, dsttempname2,ofset1, coef1, lon1,lat1,lon2,lat2, width2, height2)
ofset1=0
coef1=1.0
change_data_values_scale(srcprecname1, dstprecname1,ofset1, coef1, lon1,lat1,lon2,lat2, width1, height1)
change_data_values_scale(srcprecname2, dstprecname2,ofset1, coef1, lon1,lat1,lon2,lat2, width2, height2)


def calcu_evaporation_thornthwaite_trace_ofset(year2, tempofset1, lon1,lat1, lon2, lat2):
## maybe ok
ndays =[31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
#def calcu_daylength(dayOfYear, lat):
 #temps, lats, lons=load_temp_stack_1(year2)
 slice2=int(year2/100)
temps=[]
 for n in range(1,13):
print(n)
name1="./work7/"+chelsa_trace_name("tmean",year2, n)
temp, lats, lons=load_raster_float32(name1)
temp2=temp+tempofset1
temps.append(temp2)
 longrida, latgrida = np.meshgrid(lons, lats)


termal_index1=temps[0]*0.0
 for n in range(-1,12):
print(n+1)
tempi=temps[n]
temp5=tempi*0.2
indexmonth1=np.power(temp5, 1.514)
where_are_NaNs = np.isnan(indexmonth1)
indexmonth1[where_are_NaNs] = 0.0001
#indexmonth1[indexmonth1>0] = 0
termal_index1=termal_index1+indexmonth1
 #plt.imshow(termal_index1)
#plt.show()
 #quit(-1)


## coarse value!
#alphacoef1=0.49239+termal_index1*0.01792
alphacoef1=0.49239+termal_index1*0.01792
 ak2=np.power(termal_index1,2.0)*7.71e-5
where_are_NaNs = np.isnan(ak2)
ak2[where_are_NaNs] = 0.0
ak3=6.72e-7*np.power(termal_index1,3.0)
where_are_NaNs = np.isnan(ak3)
ak3[where_are_NaNs] = 0.0
 alphacoef1=alphacoef1+ak2+ak3


 #plt.imshow(alphacoef1)
#plt.show()


solar_evapo=[]


for n in range(-1,12):
print(n)
#inname2="./work2/solar_"+str(n+1)+".tif"
#solar1,latas, lonas=load_raster_float32(inname2)
#evapo_solaronly= solar1* 0.0864 * 0.408/1000000*ndays[n+1]
 daylengthgrida1=calcu_daylength_of_month(n, latgrida)
#print(daylengthgrida1)
alla12=daylengthgrida1/12.0
anna30=ndays[n]/30.0
anna30=1
tempa=temps[n]
tempa[tempa < 0.0] = 0.0


#print(tempa)
terma10a=(10*tempa)/termal_index1
terma10b=np.power(terma10a, alphacoef1)
 print(np.shape(terma10b))
print(np.shape(alla12))
print(np.shape(anna30))
 #pet_th=16*alla12*anna30*terma10b
##pet_nc=16.0*terma10b
pet_nc=16*alla12*anna30*terma10b
#plt.imshow(pet_nc)
#plt.show()
#quit(-1)
 evapoo1=np.array(pet_nc).astype(float)
 outname1="./work2/th_evapo1_"+str(n+1)+".nc"
evapoo1=np.flipud(evapoo1)
saveraster_float32(outname1,"NetCDF", evapoo1, lon1, lon2, lat1, lat2)
 return(solar_evapo)

get_chelsa_variable("timeseries","", "tmean", 1992, 7)
get_chelsa_variable("timeseries","", "prec", 1992, 7)
get_chelsa_variable("timeseries","", "bio", 1992, 18)
get_chelsa_variable("timeseries","", "gdd5", 1992, 0)
get_chelsa_variable("pmip3","NCAR_CCSM4", "tmean", 0, 7)
get_chelsa_variable("pmip3","NCAR_CCSM4", "bio", 0, 12)
get_chelsa_variable("pmip3","DEM", 0, 0, 0)
get_chelsa_variable("chelsa_trace","", "tasmax", 14500, 7)
get_chelsa_variable("chelsa_trace","", "bio", 14500, 12)
get_chelsa_variable("chelsa_trace","", "orog", 14500, 0)
get_chelsa_variable("chelsa_cruts","", "tmax", 1992, 7)
get_chelsa_variable("chelsa_cruts","", "prec", 1992, 7)
get_chelsa_variable("climatologies","", "prec", 0, 7)
get_chelsa_variable("climatologies","", "bio", 0, 12)
get_chelsa_variable("exchelsa","srad", "stot", 0, 7)
get_chelsa_variable("exchelsa","snow", "snow_days", 1992, 0)
get_chelsa_variable("exchelsa","pet", "pet", 0, 7)
get_chelsa_variable("exchelsa","gst", "gsl", 1992, 0)


main proggis #########################

main control wariablesdownload_chelsa_data=1

download_chelsa_data_timeseries_trace=1
download_chelsa_data_simu_climatologies=0

dataset1="climatologies"
dataset2="pmip"

dataset1="timeseries"
dataset2="trace"


simu1="NCAR_CCSM4"
set_biovars_off=1 ## set this, if cutting of bio vars stuck the program!! problem is in trace bio vars
generate_360_rasters=0

!!! Warning problems with chelsa trace bio variables !!

cut_chelsa_data=0
cut_chelsa_directory=1
rasters_to_celcius=1
downscale_chelsa_data=1
calculate_chelsa_gts5=1
calculate_chelsa_ptopet=1
calculate_simu_ptopet=0
calculate_geographical_rasters=0
chelsa_trace_geographical_rasters=0
calculate_solar=0

additional stuff ...

neo_experiment=0
gaussian_blur_test=0
own_rasters_past_present_downscaling=0

maybe not functional, developing, debug and testing code only

testbench_develop=0

https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/chelsa_trace/pr/CHELSA_TraCE21k_pr_10_-100_V1.0.tifquit(-1)settingscurrent year from downscale

year1=1993

past simu year to downscale
year2=12500

year2=14500

year2=21000
year2=18000

europe

lon1=-30.0
lon2=60.0
lat1=30.0
lat2=70.0

big europe

lon1=-60.0
lon2=90.0
lat1=30.0
lat2=80.0

very big europe, "euraasia"

lon1=-60.0
lon2=140.0
lat1=30.0
lat2=80.0

beringia

lon1=-180
lon2=-90
lat1=50
lat2=80

beringia super 360

lon1=90
lon2=270
lat1=50
lat2=80

width1=30*20
height1=10*20
width2=30*10
height2=10*10

width1=40*20
height1=10*20
width2=40*10
height2=10*10


chelsadir1="./chelsa/"
chelsadir2="./chelsa2/"
chelsadira=chelsadir1


nlon, nlat = (width1, height1)
nlon2, nlat2 = (width2, height2)
lons = np.linspace(lon1, lon2, nlon)
lats = np.linspace(lat1,lat2, nlat)
lons2 = np.linspace(lon1, lon2, nlon2)
lats2 = np.linspace(lat1,lat2, nlat2)
longrid, latgrid = np.meshgrid(lons, lats)
longrid2, latgrid2 = np.meshgrid(lons2, lats2)


if(download_chelsa_data_timeseries_trace==1):
create_dirs()
get_chelsa_trace_data(year2)
get_chelsa_timeseries_bundle(year1)
get_chelsa_variable("timeseries","", "gts5", year1, 0)
if(download_chelsa_data_simu_climatologies==1):
create_dirs()
get_chelsa_climatologies_data()
get_chelsa_pmip3_data(simu1)
#get_chelsa_cruts_data(year1)
#listaa1 = listauz("./chelsa1","*tmax*")
#listaa1 = glob.glob('./chelsa1/*tmax*.tif')
#print(listaa1)


if(generate_360_rasters==1):
rotate_rasters_180_to_360(chelsadir1, chelsadir2)


if(cut_chelsa_data==1):
mini_create_dirs()
#cut_chelsa_dir(chelsadira,lon1, lat1, lon2, lat2, width1, height1, width2, height2)
 cut_chelsa_bundle("timeseries", "", year1, 0, lon1, lat1, lon2, lat2, width1, height1, width2, height2)
#slice1=year2
cut_chelsa_bundle("trace", "", year2, 0, lon1, lat1, lon2, lat2, width1, height1, width2, height2)
if(cut_chelsa_directory==1):
mini_create_dirs()
cut_chelsa_dir(chelsadira,lon1, lat1, lon2, lat2, width1, height1, width2, height2)
if(rasters_to_celcius==1):
generate_celcius_rasters(dataset1, dataset2, simu1,year1, year2, lon1,lat1,lon2,lat2, width1, height1, width2, height2 )


if(downscale_chelsa_data==1):
# downscale gdd
trace_tmean_from_tasmax_tasmin(year2)
downscale_timeseries_trace("./output1/gts5_downscaled.nc","timeseries", "trace", "", year1, year2)


if(calculate_chelsa_gts5==1):
temps00=load_temp_stack_1(year2)
temps, lons, lats=np.array(temps00)
gts5=gts(temps, 5.0)
#plt.imshow(gts5)
#plt.show()
saveraster_float32("./output1/gts5_calculated.nc","NetCDF", gts5, lon1, lon2, lat1, lat2)
if(calculate_chelsa_ptopet==1):
#calcu_evaporation_thornthwaite_2(year2, lon1,lat1, lon2, lat2)
# temperature bias, ofset is -9
calcu_evaporation_thornthwaite_trace_ofset(year2, 0, lon1,lat1, lon2, lat2)
precips, lats1, lons1=load_precip_stack_1(year2)
evaps, lats2, lons2=load_evap_stack_1()
annprecip=precips[0]*0.0
annevap=evaps[0]*0.0
for n in range(-1,12):
annprecip=precips[n]+annprecip
annevap=evaps[n]+annevap
 ptopet00=annprecip/annevap
 ptopet01=ptopet00
ptopet01[ptopet01>10.0]=10.0
ptopet01[ptopet01<0.0]=0.0
ptopet1=np.flipud(ptopet01)
#saveraster_float32("./output1/ptopet_calculated.tif","GTiff", ptopet1, lon1, lon2, lat1, lat2)
saveraster_float32("./output1/ptopet_calculated.nc","NetCDF", ptopet1, lon1, lon2, lat1, lat2)
if(calculate_simu_ptopet==1):
calcu_evaporation_thornthwaite_simu(year2, lon1,lat1, lon2, lat2)
#precips, lats1, lons1=load_precip_stack_1(year2)
precips=[]
 for n in range(1,13):
print(n)
#name1="./work4/CHELSA_prec_1993_"+str(n).zfill(2)+"_V1.2.1.tif"
## note correct "01" in precip files to !
name1="./work2/CHELSA_PMIP_CCSM4_prec_"+str(n)+"_1.tif"
#print(name1)
precip, lons, lats=load_raster_float32(name1)
#temp2=(temp*0.1)-273.16
precips.append(precip)
 evaps, lats2, lons2=load_evap_stack_1()
annprecip=precips[0]*0.0
annevap=evaps[0]*0.0
for n in range(-1,12):
annprecip=precips[n]+annprecip
annevap=evaps[n]+annevap
 #ptopet00=annprecip/annevap
#propably p in files is Px10
ptopet00=annprecip*0.1/annevap
 ptopet01=ptopet00
ptopet01[ptopet01>10.0]=10.0
ptopet01[ptopet01<0.0]=0.0
ptopet1=np.flipud(ptopet01)
#saveraster_float32("./output1/ptopet_calculated.tif","GTiff", ptopet1, lon1, lon2, lat1, lat2)
saveraster_float32("./output1/ptopet_calculated.nc","NetCDF", ptopet1, lon1, lon2, lat1, lat2)




if(calculate_geographical_rasters==1):
cut_and_scale_raster_file("./etopo1/ETOPO1_Ice_c_geotiff.tif", "etopo1.tif",lon1, lat1, lon2, lat2, width1, height1, width2, height2)
create_landsea("./work1/etopo1.tif", "./work2/landsea_etopo_current.tif",lon1, lat1, lon2, lat2, width1, height1)
create_landsea("./work4/etopo1.tif", "./work4/landsea_etopo_current.tif",lon1, lat1, lon2, lat2, width2, height2)
saveraster_float32("./work2/longrid_current.tif","GTiff", longrid, lon1, lon2, lat1, lat2)
saveraster_float32("./work2/latgrid_current.tif","GTiff", latgrid, lon1, lon2, lat1, lat2)
saveraster_float32("./work2/longrid_current.nc","NetCDF", longrid, lon1, lon2, lat1, lat2)
saveraster_float32("./work2/latgrid_current.nc","NetCDF", latgrid, lon1, lon2, lat1, lat2)
saveraster_float32("./work4/longrid_current.tif","GTiff", longrid2, lon1, lon2, lat1, lat2)
saveraster_float32("./work4/latgrid_current.tif","GTiff", latgrid2, lon1, lon2, lat1, lat2)
saveraster_float32("./work4/longrid_current.nc","NetCDF", longrid2, lon1, lon2, lat1, lat2)
saveraster_float32("./work4/latgrid_current.nc","NetCDF", latgrid2, lon1, lon2, lat1, lat2)
sea_proximity("./work2/landsea_etopo_current.tif", "./work2/distance_current.tif")
sea_proximity("./work4/landsea_etopo_current.tif", "./work4/distance_current.tif")
outfile1="rivers.tif"
shpath1="./shape/natural_earth/ne_10m_rivers_lake_centerlines.shp"
#shpath1="./shape/GSHHS_shp/i/GSHHS_i_L3.shp"
#shpath1="./shape/unesco/rivers.shp"
rasterize_shapefile(shpath1, "./work2/rivers.tif", lon1, lat1, lon2, lat2, nlon, nlat)
rasterize_shapefile(shpath1, "./work4/rivers.tif", lon1, lat1, lon2, lat2, nlon2, nlat2)
river_proximity("./work2/rivers.tif", "./work2/distance_current_rivers.tif")
river_proximity("./work4/rivers.tif", "./work4/distance_current_rivers.tif")


if(chelsa_trace_geographical_rasters==1):
slice2=int(year2/100)
slice2s=str(slice2)
get_chelsa_variable("chelsa_trace","", "orog", year2, 0)
demname1="glacier_plus_dem_"+slice2s+"_V1.2.tif"
neofile1="./chelsa/"+demname1
#./chelsa/glacier_plus_dem_180_V1.2.tif
 #neofile1=""
neofile2=demname1
cut_and_scale_raster_file(neofile1, neofile2,lon1, lat1, lon2, lat2, width1, height1, width2, height2)
create_landsea("./work2/"+demname1, "./work2/landsea_glacier.tif",lon1, lat1, lon2, lat2, width1, height1)
create_landsea("./work4/"+demname1, "./work4/landsea_glacier.tif",lon1, lat1, lon2, lat2, width2, height2)
sea_proximity("./work2/landsea_glacier.tif", "./work2/distance_glacier.tif")
sea_proximity("./work4/landsea_glacier.tif", "./work4/distance_glacier.tif")
#gdaldem_test1("hillshade", "NetCDF", "./work2/etopo1.tif", "out.nc", 270,30)


if(gaussian_blur_test==1):
src1="./output1/ptopet_calculated.tif"
dst1="./output1/ptopet_calculated_blur20.tif"
dst2=str.replace(dst1,"tif", "nc")
inimage1, lats, lons=load_raster_float32(src1)
#inimage2=np.array(inimage1).astype(double)
result1 = gaussian_filter(inimage1, sigma=5)
result2=np.array(result1).astype(float)
saveraster_float32(dst1,"GTiff", result1, lon1, lon2, lat1, lat2)
saveraster_float32(dst2,"NetCdf", result1, lon1, lon2, lat1, lat2)
if(neo_experiment==1):
print("Stub")
## neo, rgb experiment
inputneo1="./neo/MOD_NDVI_M_2020-07-01_rgb_3600x1800.FLOAT.TIFF"
#inputneo1="./neo/MOD_NDVI_M_2020-07-01_rgb_3600x1800.FLOAT.TIFF"
#inputneo1="./neo/UiO_PEX_PERPROB_5.0_20181128_2000_2016_NH.tif"
#inputneo1="./neo/UiO_PEX_MAGTM_5.0_20181127_2000_2016_NH.tif"
inputneo1="./neo/MOD15A2_M_LAI_2015-07-01_rgb_3600x1800.FLOAT.TIFF"
#inputneo1="./neo/CHELSA_gts5_1979_V1.2.1.tif"
#inputneo1="./neo/CHELSA_gts5_2010_V1.2.1.tif"
#referenceneo1="CHELSA_temp10_01_1979-2013_V1.2_land.tif"
#inputneo1="./neo/PermafrostNSIDC_2002-02-01_rgb_3600x1800.TIFF"
#inputneo1="./neo/BlueMarbleNG_2004-08-01_rgb_3600x1800.TIFF"
 referenceneo1="CHELSA_tmean_1993_07_V1.2.1.tif"
outputneo1="neo_p.tif"
outputneo2="neo_rgb.tif"
 #cut_and_scale_rgb_file(inputneo1, outputneo2,lon1, lat1, lon2, lat2, width1, height1, width2, height2)


#process_neo_raster(inputneo1, outputneo1, referenceneo1)
 #downscale_vars_directly("./work2/neo_p.tif", "./work4/neo_p.tif", "out_p.nc", "temp","tmean", lon1, lon2, lat1, lat2)
 #downscale_neo("timeseries", "", year1)
downscale_timeseries_trace_neo("./output1/neo_ndvi.nc",outputneo1, "", year1, year2)
 #rgb_downscale_timeseries_trace("timeseries", "trace", "", year1, year2, lon1, lat1, lon2, lat2, width1, height1)
## neo experiment
#process_neo_raster()
#downscale_neo("timeseries", "", year1)
if(calculate_solar==1):
print("Calcu solar, maybe slooow ...")
for n in range(1,12):
print(n)
#monthrad=monthly_solar_incoming_radiation(year1, n,latgrid)
#monthrad2=monthly_solar_incoming_radiation(year1, n,latgrid2)
monthrad=climlab_calculate_insolation_milankovic(year2,n,lat1, lat2, width1, height1)
monthrad2=climlab_calculate_insolation_milankovic(year2,n,lat1, lat2, width2, height2)
name1="./work2/solar_"+str(n)+".tif"
name2="./work2/solar_"+str(n)+".nc"
bname1="./work4/solar_"+str(n)+".tif"
bname2="./work4/solar_"+str(n)+".nc"
#saveraster_float32(name1,"GTiff", monthrad, lon1, lon2, lat1, lat2)
#saveraster_float32(name2,"NetCDF", monthrad, lon1, lon2, lat1, lat2)
saveraster_float32(bname1,"GTiff", monthrad2, lon1, lon2, lat1, lat2)
saveraster_float32(bname2,"NetCDF", monthrad2, lon1, lon2, lat1, lat2)
#monthrads.append(monthrads)
#yearad=yearad+monthrad
if(own_rasters_past_present_downscaling==1):
process_my_own_rasters("./juma1", "./juma2", "./work4/CHELSA_gts5_1993_V1.2.1.tif")
#process_my_own_rasters("./jama1", "./jama2", "./work4/CHELSA_gts5_1993_V1.2.1.tif")
#save_points_to_netcdf("outs.nc", "Grape", lons, lats, longrid)
#downscale_own_dirs("./output1/own_ndvi.nc","./work2/neo_ndvi.tif", "./work4/neo_ndvi.tif", "./jama1", "./jama2", "./juma1", "./juma2")

testbench site



if(testbench_develop==1):
## develop code only, maube not functional!
calculate_tmean_from_chelsa_cruts_data("chelsa_cruts", "", year1, lon1, lat1, lon2, lat2, width1, height1, width2, height2)
quit(-1)
bigs1, lo, la=loadraster_stak("./work2/", "tmean",".tif")
heina1=bigs1[6]
heina1[heina1<-400.0]=-400
heina1[heina1>300.0]=0.0
plt.imshow(heina1,plt.cm.Reds,extent=[min(la),max(la),min(lo), max(lo)] )
#cmap=plt.cm.Reds, interpolation='none',
plt.show()
quit(-1)
##-- testbench


print(".")
quit(-1)

Licensing

I, the copyright holder of this work, hereby publish it under the following license:

This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.

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.
share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.

File history

Click on a date/time to view the file as it appeared at that time.

	Date/Time	Thumbnail	Dimensions	User	Comment
current	14:27, 16 February 2021		1,328 × 880 (507 KB)	Merikanto	Uploaded own work with UploadWizard

File usage

No pages on the English Wikipedia use this file (pages on other projects are not listed).

File:Europe chelsa trace tjuly 14500BP 1.png

Summary

Licensing

Captions

Items portrayed in this file

depicts

creator

some value

copyright status

copyrighted

copyright license

Creative Commons Attribution-ShareAlike 4.0 International

source of file

original creation by uploader

inception

16 February 2021

coordinates of the point of view

50°N, 1°E

media type

image/png

File history

File usage