File:Ordovician ice age 445ma tas annual S1312 co2 120 tilt22 ecc 0p06 mvelp 90 1.png

1. 3. Exoplasim planet running code
   4. exoplasim example
   5. stepper code
2. 2. 01.06.2023 0000.0004b
3. 1. in ubuntu you must install
4. 1. pip3 install exoplasim[netCDF4]
   2. not
   3. "sudo pip3 install exoplasim[netCDF4]"

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

from matplotlib.pyplot import figure, draw, pause

from scipy.interpolate import interp2d import netCDF4

import exoplasim as exo

NLAT=0 NLON=0

def plot_temperature(model, fig):

   lon = model.inspect("lon")
   lat = model.inspect("lat")
   ts = model.inspect("tas",tavg=True)
   daata1=np.array(ts)-273.15
   print("dim ", np.shape(ts))
   dimx1=64
   dimy1=32
   meantas1=np.mean(daata1)
   #print(meantas1)
   ax = fig.gca()
   ax.clear()
   ax.set_title("Temperature mean "+ str(round(meantas1, 0)), fontsize=22)  
   alons1=np.linspace(-180,180,num=dimx1)
   alats1=np.linspace(-90,90,num=dimy1)
   levels1=[-250,-200,-150,-120,-100,-80,-50,-40,-30,-20,-10,0,10,20,30,40,50,60,80,100,120, 150,200,300]
   ext1=[-180,180,-90, 90]
   vmin1=-60
   vmax1=60
   ax.imshow(daata1, extent=ext1, cmap="coolwarm", origin="lower", vmin=vmin1, vmax=vmax1)
   cs1 = ax.contour( alons1, alats1, daata1,origin="lower", extent=ext1, colors=['k'], levels=levels1)
   ax.clabel(cs1, cs1.levels, inline=True, fmt="%3.1f", fontsize=15)
   ax.xaxis.set_tick_params(labelsize=18)
   ax.yaxis.set_tick_params(labelsize=18)
   #ax.set_xticklabels(fontsize=16)
   #ax.set_xlabel(fontsize=16)
   fig.canvas.draw() ##important for animation
   fig.canvas.flush_events()	## important for animation
   return(0)

def setup_monska(dimx1, dimy1):

   # imgs is a N x X x Y image stack
   plt.ion() ## important for animation
   fig = figure()
   ax = fig.gca()
   #img=np.random.rand(dimx1*dimy1).reshape((dimy1,dimx1))*100
   #hoo = ax.imshow(img, cmap='rainbow')
   plt.show()
   return(fig)

def monitor_temperaturepzka(model): lon = model.inspect("lon") lat = model.inspect("lat") ts = model.inspect("tas",tavg=True) im=plt.pcolormesh(lon,lat,ts,cmap="RdBu_r",vmin=273.15-60.0,vmax=273.15+60.0,shading="Gouraud") plt.contour(lon,lat,ts-273.15,[-30,-20,-10,0,10,20,30],colors=['gray',]) plt.colorbar(im,label="Surface Temperature [K]") plt.xlabel("Longitude [deg]") plt.ylabel("Latitude [deg]") plt.title("Surface air Temperature") plt.show()

def writeSRA(name,kcode,field,NLAT,NLON):

   label=name+'_surf_%04d.sra'%kcode
   header=[kcode,0,20170927,0,NLON,NLAT,0,0]
   fmap = field.reshape((int(NLAT*NLON/8),8))
   sheader = 
   for h in header:
       sheader+=" %11d"%h
   
   lines=[]
   i=0
   while i<NLAT*NLON/8:
       l=
       for n in fmap[i,:]:
           l+=' %9.3f'%n
       lines.append(l)
       i+=1

   text=sheader+'\n'+'\n'.join(lines)+'\n' 

   f=open(label,'w')
   f.write(text)
   f.close()
   print (label)

def writeSRA2(label,kcode,field,NLAT,NLON):

   #label=name+'_surf_%04d.sra'%kcode
   header=[kcode,0,20170927,0,NLON,NLAT,0,0]
   fmap = field.reshape((int(NLAT*NLON/8),8))
   sheader = 
   for h in header:
       sheader+=" %11d"%h
   
   lines=[]
   i=0
   while i<NLAT*NLON/8:
       l=
       for n in fmap[i,:]:
           l+=' %9.3f'%n
       lines.append(l)
       i+=1

   text=sheader+'\n'+'\n'.join(lines)+'\n' 

   f=open(label,'w')
   f.write(text)
   f.close()
   print (label)

def savenetcdf_single_frommem(outfilename1, outvarname1, xoutvalue1,xoutlats1,xoutlons1): nlat1=len(xoutlats1) nlon1=len(xoutlons1) #indata_set1=indata1 print(outfilename1) ncout1 = netCDF4.Dataset(outfilename1, 'w', format='NETCDF4') outlat1 = ncout1.createDimension('lat', nlat1) outlon1 = ncout1.createDimension('lon', nlon1) outlats1 = ncout1.createVariable('lat', 'f4', ('lat',)) outlons1 = ncout1.createVariable('lon', 'f4', ('lon',)) outvalue1 = ncout1.createVariable(outvarname1, 'f4', ('lat', 'lon',)) outvalue1.units = 'Unknown' outlats1[:] = xoutlats1 outlons1[:] = xoutlons1 outvalue1[:, :] =xoutvalue1[:] ncout1.close() return 0

def loadnetcdf_single_tomem(infilename1, invarname1): global cache_lons1 global cache_lats1 print(infilename1) inc1 = netCDF4.Dataset(infilename1) inlatname1="latitude" inlonname1="longitude" inlats1=inc1[inlatname1][:] inlons1=inc1[inlonname1][:] cache_lons1=inlons1 cache_lats1=inlats1 indata1_set1 = inc1[invarname1][:] dim1=indata1_set1.shape nlat1=dim1[0] nlon1=dim1[1] inc1.close() return (indata1_set1)

def create_sras(topo, seamasklevel1):

global NLAT global NLON

topo2=np.copy(topo)

seamasklevel2=seamasklevel1+1.0 topo2[topo2 < seamasklevel1] = seamasklevel1 masko=np.copy(topo2) masko[masko == seamasklevel1] = -9999999 masko[masko > seamasklevel1] = 1 masko[masko == -9999999 ] = 0

grid=np.flipud(masko) name="Example" writeSRA(name,129,topo,NLAT,NLON) writeSRA(name,172,grid,NLAT,NLON) writeSRA2("topo.sra",129,topo2,NLAT,NLON) writeSRA2("landmask.sra",172,grid,NLAT,NLON) return(0)

def convert_to_t21(infilename1, outfilename1, seamasklevel1):

global NLAT global NLON

indimx=361 indimy=181 #indimx=360 #indimy=360

## t21 64x32 shapex=64 shapey=32 NLAT=shapex NLON=shapey nc = netCDF4.Dataset(infilename1)

inlats=nc['latitude'][:] inlons=nc['longitude'][:] #print(inlats) #print(inlons) latlen=len(inlats) lonlen=len(inlons)

#print(lonlen, latlen)

indimx=lonlen indimy=latlen

dem000=nc['z'] dem=np.flipud(dem000) #dem=dem000 dem2=np.copy(dem) #dem2[dem2 < 0] = 0 #plt.imshow(dem,cmap='gist_earth') #plt.imshow(dem2,cmap='gist_earth') #plt.show() #quit(0) lts0=[85.7606, 80.2688, 74.7445, 69.2130, 63.6786, 58.1430, 52.6065, 47.0696, 41.5325,35.9951, 30.4576, 24.9199, 19.3822, 13.8445, 8.3067, 2.7689, -2.7689, -8.3067, -13.8445, -19.3822, -24.9199, -30.4576, -35.9951, -41.5325, -47.0696, -52.6065, -58.1430, -63.6786, -69.2130, -74.7445, -80.2688, -85.7606]

## lns0=[0, 5.6250, 11.2500, 16.8750, 22.5000, 28.1250, 33.7500 ,39.3750, 45.0000, 50.6250, 56.2500, 61.8750, 67.5000, 73.1250, 78.7500, 84.3750, 90.0000, 95.6250, 101.2500, 106.8750, 112.5000, 118.1250, 123.7500, 129.3750, 135.0000, 140.6250, 146.2500, 151.8750, 157.5000, 163.1250, 168.7500, 174.3750, 180.0000, 185.6250, 191.2500, 196.8750, 202.5000, 208.1250, 213.7500, 219.3750, 225.0000, 230.6250, 236.2500, 241.8750, 247.5000, 253.1250, 258.7500, 264.3750, 270.0000, 275.6250, 281.2500, 286.8750, 292.5000, 298.1250, 303.7500, 309.3750, 315.0000, 320.6250, 326.2500, 331.8750, 337.5000, 343.1250, 348.7500, 354.3750]

lts1=np.array(lts0) lns1=np.array(lns0)

lns=lns1 lts=np.flip(lts1)

ly2=len(lts) lx2=len(lns) shapex=lx2 shapey=ly2

#print("sheip") #print(shapex, shapey)

lons, lats = np.meshgrid(lns,lts) #print (lts) #print (lns) new_W, new_H = (shapey,shapex) xrange = lambda x: np.linspace(0, 360, x) f2 = interp2d(xrange(indimx), xrange(indimy), dem2, kind="linear") #f2 = interp2d(range(indimx), range(indimy), dem2, kind="cubic") demo = f2(xrange(shapex), xrange(shapey)) #plt.imshow(demo) #plt.show() #quit(0) f3 = interp2d(xrange(indimx), xrange(indimy), dem2, kind="linear") #masko = f3(xrange(shapex), xrange(shapey)) #topo=np.flipud(demo) topo=np.copy(demo) topo2=np.copy(topo) masko=np.copy(topo)

seamasklevel2=seamasklevel1+1.0

topo2[topo2 < seamasklevel1] = seamasklevel1

masko=np.copy(topo2) masko[masko == seamasklevel1] = -9999999 masko[masko > seamasklevel1] = 1 masko[masko == -9999999 ] = 0 #plt.imshow(demo) #plt.imshow(masko) #plt.imshow(topo2) #plt.show()

#grid=np.fliplr(masko) #def savenetcdf_single_frommem(outfilename1, outvarname1, xoutvalue1,xoutlats1,xoutlons1): savenetcdf_single_frommem(outfilename1, "z", topo,lts,lns) savenetcdf_single_frommem("mapt21.nc", "z", topo2,lts,lns) savenetcdf_single_frommem("maskt21.nc", "z", masko,lts,lns) return(topo,lons,lats)

1. 1. stepper exoplasim ...

def run_exoplasim_wyb(a_input_dem1,s_seamasklevel1, a_gridtype, a_layers, a_years,a_timestep,a_snapshots,a_ncpus,a_eccentricity,a_obliquity,a_lonvernaleq,a_pCO2,a_startemp,a_baseflux,a_yearlength,a_radius,a_gravity,a_rotationperiod):

   print("Exoplasim runner .")

   pressure1=1.0
   a_pO2=(1-a_pCO2-0.7808)*pressure1
   a_pN2=(1-0.2095-a_pCO2)*pressure1
   a_pCO2=a_pCO2*pressure1
   
   output_format=".nc"

   print("Process input grid, to type ",a_gridtype)

   if(a_gridtype=="T21"):
       print("T21")
       topo, lons, lats=convert_to_t21(a_input_dem1,"demT21.nc", a_seamasklevel1)    

   print("SRA")
   create_sras(topo, a_seamasklevel1)
   print("Creating exoplasim object ")
   fig=setup_monska(64, 32)

   testplanet= exo.Model(workdir="planet_run",modelname="PLANET",ncpus=a_ncpus,resolution=a_gridtype,layers=a_layers, outputtype=output_format, crashtolerant=True)

   glaciers1= {

"toggle": True, "mindepth":1, "initialh":-1 }

   # !!!! fixed rotation!
   testplanet.configure(

startemp=a_startemp, flux=a_baseflux,# Stellar parameters eccentricity=a_eccentricity, obliquity=a_obliquity, lonvernaleq=a_lonvernaleq, year=a_yearlength, fixedorbit=False, # Orbital parameters rotationperiod=a_rotationperiod, # Check if otation IS FIXED!!! synchronous=False, # CHECK if rotation IS FIXED!!! topomap="topo.sra", landmap="landmask.sra", radius=a_radius, gravity=a_gravity, # Bulk properties #seaice=False, #maxsnow=False, #glaciers=False, stormclim=True, #vegetation=0, wetsoil=True, #alters albedo of soil based on how wet it is

               vegetation=1,                               #toggles vegetation module; 1 for static vegetation, 2 to allow growth
               vegaccel=1, 

seaice=True, maxsnow=-1, glaciers=glaciers1, #stormclim=True, #vegetation=0, pN2=a_pN2, pCO2=a_pCO2, pO2=a_pO2, ozone=True, # Atmosphere timestep=a_timestep, snapshots=0, ## jos a_snapshots, vie muistia! #wetsoil=True, physicsfilter="gp|exp|sp") # Model dynamics

   # ru
   testplanet.exportcfg()
   looplen1=a_runsteps1
   peen=0
   runc1=1

   print("Phase 2 !!! Stepper runner.")	
   for n in range(0,looplen1):
       print("Exoplasim runner year ",n)
       a_years2=1
       runc1=1
       testplanet.run(years=1,crashifbroken=True)
       savename = 'planet_run_'+str(runc1)
       testplanet.finalize(savename,allyears=False,clean=False,keeprestarts=True)
       testplanet.save(savename)
       tas=testplanet.inspect("tas")
       mint=testplanet.inspect("mint")
       maxt=testplanet.inspect("maxt")
       tas1=np.ravel(tas)
       maxt1=np.ravel(maxt)
       mint1=np.mean(mint)
       meantas1=np.mean(tas1)-273.15
       mintas1=np.min(tas1)-273.15
       maxtas1=np.max(tas1)-273.15
       meanmint1=np.mean(mint1)-273.15
       minmint1=np.min(mint1)-273.15
       maxmint1=np.max(mint1)-273.15
       meanmaxt1=np.mean(maxt1)-273.15
       minmaxt1=np.min(maxt1)-273.15
       maxmaxt1=np.max(maxt1)-273.15
       dailydeltat1=meanmaxt1-meanmint1
       globedeltat1=maxtas1-mintas1
       print("Tas mean min max             : ", round(meantas1,2),round(mintas1,2),round(maxtas1,2) )
       print("Pole-eqv deltaT Daily deltaT : ", round(globedeltat1,2),round(dailydeltat1,2)  )
       plot_temperature(testplanet,fig)
   
   print("Return.")

   return(0)

print(" Exoplasim simulation code ---")

1. !!! Sun

1. input_dem="dema1.nc" ##dem of exoplanet
2. input_dem="elev1.nc" ##dem of exoplanet

input_dem="map445Ma.nc" ##dem of exoplanet

a_modelname1="planet" a_workdir1="planet_run"

a_runsteps1=300 a_years1=a_runsteps1 a_timestep1=15 a_snapshots1=0 a_ncpus1=4 a_layers1=4 a_outputtype1=".nc"

1. a_resolution1="T42"

a_resolution1="T21" a_precision1=4 a_crashtolerant1=True a_landmap1="landmask.sra" a_topomap1="topo.sra" a_seamasklevel1=0.0

a_startemp=5772

1. 1. a_baseflux=1365.2
2. a_baseflux=1354
3. a_baseflux=1318

a_baseflux=1312

a_yearlength=360 a_radius=1.0 a_gravity=9.81

1. 1. 1. !!! ROTATION 7 !!!!

a_rotationperiod=0.883

1. 1. !!!!! earth nowadays
2. a_eccentricity1=0.01671022
3. a_obliquity1=23.4392811
4. a_lonvernaleq1=102.94719

a_eccentricity1=0.06 a_obliquity1=22 a_lonvernaleq1=90

a_pCO21=120e-6

print("Exoplasim ...")

1. 1. attempt to run exoplasim stepper code

run_exoplasim_wyb(input_dem, a_seamasklevel1, a_resolution1, a_layers1, a_years1,a_timestep1,a_snapshots1,a_ncpus1,a_eccentricity1,a_obliquity1,a_lonvernaleq1,a_pCO21,a_startemp,a_baseflux,a_yearlength,a_radius,a_gravity,a_rotationperiod)

print(".")

2. 1. temperature plotter
3. python 3
   1. 06.06.2023 0000.0004

import matplotlib.mlab as mlab import matplotlib.pyplot as plt import netCDF4 from netCDF4 import Dataset from matplotlib.pylab import * import numpy as np from scipy import interpolate from scipy.interpolate import griddata

from sklearn.linear_model import LinearRegression from pygam import LogisticGAM from pygam import LinearGAM from sklearn.ensemble import RandomForestRegressor from sklearn.preprocessing import StandardScaler from sklearn.preprocessing import PolynomialFeatures from sklearn.pipeline import make_pipeline

1. from skimage.transform import resize

import skimage

1. 1. random forest setting

RF_estimators1=1 RF_features1=1

def gam_multiple_vars( x1, y1, x2): print ("GAM") xa1=np.array(x1) ya1=np.array(y1) xa2=np.array(x2)

#print (xa1)

1. 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 random_forest_multiple_vars( x1, y1, x2): print ("RF")

global RF_estimators1 global RF_features1

print(RF_estimators1, RF_features1)

#quit(-1)

xa1=np.array(x1) ya1=np.array(y1) xa2=np.array(x2)

#print (xa1)

1. quit(-1)

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)

sc = StandardScaler() xa = sc.fit_transform(x) xba = sc.transform(xb)

## orig ##model = RandomForestRegressor(n_estimators=10, max_features=2).fit(xa,y) model = RandomForestRegressor(n_estimators=RF_estimators1, max_features=RF_features1).fit(xa,y)

y2= model.predict(xba) 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 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 = 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 savenetcdf_single_frommem(outfilename1, outvarname1, xoutvalue1,xoutlats1,xoutlons1): nlat1=len(xoutlats1) nlon1=len(xoutlons1) #indata_set1=indata1 print(outfilename1) ncout1 = netCDF4.Dataset(outfilename1, 'w', format='NETCDF4') outlat1 = ncout1.createDimension('lat', nlat1) outlon1 = ncout1.createDimension('lon', nlon1) outlats1 = ncout1.createVariable('lat', 'f4', ('lat',)) outlons1 = ncout1.createVariable('lon', 'f4', ('lon',)) outvalue1 = ncout1.createVariable(outvarname1, 'f4', ('lat', 'lon',)) outvalue1.units = 'Unknown' outlats1[:] = xoutlats1 outlons1[:] = xoutlons1 outvalue1[:, :] =xoutvalue1[:] ncout1.close() return 0

def loadnetcdf_single_tomem(infilename1, invarname1): global cache_lons1 global cache_lats1 print(infilename1) inc1 = netCDF4.Dataset(infilename1) inlatname1="latitude" inlonname1="longitude" inlats1=inc1[inlatname1][:] inlons1=inc1[inlonname1][:] cache_lons1=inlons1 cache_lats1=inlats1 indata1_set1 = inc1[invarname1][:] dim1=indata1_set1.shape nlat1=dim1[0] nlon1=dim1[1] inc1.close() return (indata1_set1)

2. main

1. kaption1="Temperature of July 67 million years ago, degrees Celsius"

kaption1="Mean temperature of year degrees Celsius, final Ordovician 445 Ma "

1. kaption2="S=1312 CO2=180 ppm obliq=22 ecc=0.06 mvelp=90"

1. kaption2="S=1312 CO2=150 ppm obliq=22 ecc=0.06 mvelp=90"

kaption2="S=1312 CO2=120 ppm obliq=22 ecc=0.06 mvelp=90"

outname1="ordovician_445ma_tas_1.png"

1. outname1="cretaceous_67ma_co2_4x_tjanuary_1.png"

year1=0 month1=0 ## if 0, mean surf temp of year

1. month1=1

vari1="tas"

tsa_offset=0

1. fintemp1="./exoplasim_445ma_co2_125/MOST.00040.nc"
2. fintemp1="./exoplasim_445ma_co2_180/MOST.00050.nc"
3. fintemp1="./exo_445_omega45_co2_150/MOST.00030.nc"
4. fintemp1="./exo_445_omega270_co2_120/MOST.00024.nc" ## nok

fintemp1="./exo_445_omega90_co2_120/MOST.00020.nc"

1. findem1="Map22_PALEOMAP_6min_Mid-Cretaceous_95Ma.nc"

findem1="./maps/Map77_PALEOMAP_6min_Late_Ordovician_445Ma.nc"

1. 1. downscaleing size

dimy3=360 dimx3=720

extent1=(-180,180,-90,90)

1. extent1=(-180,180,90,-90)

1. print(num1, kk1)

1. print(fintemp1)

kindex1=1

ncinu1 = Dataset(fintemp1, "r+", format = "NETCDF4")

print(ncinu1)

1. quit(-1)

ncdem1 = Dataset(findem1, "r+", format = "NETCDF4")

1. print(ncinu1)

1. quit(-1)

1. navalue=9.969209968386869e+36

lon0 = ncinu1.variables['lon'] lat0 = ncinu1.variables['lat'] lon1 = ncdem1.variables['longitude'] lat1 = ncdem1.variables['latitude']

dimx1=len(lon0) dimy1=len(lat0)

dimx2=len(lon1) dimy2=len(lat1)

tsa00 = np.array(ncinu1.variables[vari1])

1. snd = np.array(ncinu1.variables['snd'])

tsa=np.mean(tsa00, axis=0)

1. if(month1==0): tsa=np.mean(tsa, axis=0)
2. else:tsa=tsa[month1-1]

tsa=tsa-273.15

tsa_january=np.array(ncinu1.variables['tsa'][0])-273.15 tsa_july=np.array(ncinu1.variables['tsa'][6])-273.15

mask00=ncinu1.variables['lsm'][1] mask1=np.array(mask00)

snowdepths=np.array(ncinu1.variables['snd'][:]) minsnowdepth=np.min(snowdepths, axis=0) maxsnowdepth=np.max(snowdepths, axis=0)

1. plt.imshow(maxsnowdepth)

1. plt.contour(minsnowdepth, levels=[0])

1. plt.imshow(mask00, cmap="terrain")

1. plt.contour(minsnowdepth, levels=[0])

1. plt.imshow(tsa_january, cmap="jet")

1. plt.imshow(tsa_july, cmap="jet")
2. plt.imshow(tsa, cmap="jet")

1. plt.show()

1. quit(-1)

dem00=np.array(ncdem1.variables['z'])

dem01=np.copy(dem00)

dem1=skimage.transform.resize(dem01,(dimy3, dimx3))

1. 1. warning
   2. dem1=np.flipud(dem1)

1. plt.imshow(np.array(ncinu1.variables['tso'][1]), extent=extent1)

1. plt.imshow(tsa_january, cmap="terrain")

1. plt.imshow(dem1, cmap="terrain")

1. plt.show()

1. quit(-1)

1. print(tsa00)

smalons=np.zeros((dimy1, dimx1))*255 smalats=np.zeros((dimy1, dimx1))*255

for iy in range(dimy1): for ix in range(dimx1): sx=ix/dimx1 sy=iy/dimy1 smalons[iy, ix]=sx smalats[iy, ix]=sy

biglons=np.zeros((dimy3, dimx3))*255 biglats=np.zeros((dimy3, dimx3))*255

for iy in range(dimy3): for ix in range(dimx3): sx=ix/dimx3 sy=iy/dimy3 biglons[iy, ix]=sx biglats[iy, ix]=sy

1. plt.imshow(smalons)
2. plt.imshow(smalats)

1. plt.imshow(biglons)
2. plt.imshow(biglats)

1. plt.show()

1. 1. resize input dem

1. print (len(dem00))

kmap2 = plt.cm.get_cmap('RdBu_r') kmap1=plt.get_cmap('YlGnBu')

dmap1=plt.get_cmap('binary')

1. kmap2 = plt.cm.get_cmap('bwr')

kmap2 = plt.cm.get_cmap('coolwarm') ## ok

1. kmap2 = plt.cm.get_cmap('jet')

1. kmap2 = plt.cm.get_cmap('turbo') ##ok
2. kmap2 = plt.cm.get_cmap('gist_rainbow_r')

lon1=np.array(lon0) lat1=np.array(lat0)

1. print(lon1)

1. tsa=np.array(tsa00)
2. mask1=np.array(mask00) ## ok

1. print(np.shape(tsa[month1]))

1. tsa=tsa+tsa_offset

1. plt.imshow(tsa, cmap=kmap2, interpolation="bicubic")
2. plt.imshow(dem00, cmap=kmap2, interpolation="bicubic")

1. plt.show()

1. quit(-1)

dem2=np.copy(dem1)

dem2[dem2>-1]=1 dem2[dem2<0]=0

1. sla=[]

1. sla.append(tsa)

1. apoints1=
2. bpoints1=list(zip(*sla))

dem4=np.copy(dem1) dem4[dem4<0]=0

masco1=np.copy(dem1) masco1[masco1<0]=0 masco1[masco1>0]=255

1. plt.imshow(masco1)
2. plt.imshow(mask1)

1. plt.imshow(tsa)

1. plt.imshow(dem4)
2. plt.show()

1. quit(-1)

smalldem = skimage.transform.resize(dem4,(dimy1, dimx1)) bigsmalldem=skimage.transform.resize(smalldem,(dimy3, dimx3), anti_aliasing=True)

demdifference=bigsmalldem-dem4

tempdifference=demdifference*10/1000 ## simple assumption lapse rate 10 or 6.5 / 1000 m

1. 1. add some random noise

mu, sigma = 0, 0.3 # mean and standard deviation

random0 = np.random.normal(mu, sigma, (int(dimy3/10), int(dimx3/10))) random1=skimage.transform.resize(random0, (dimy3, dimx3), anti_aliasing=True ) mu1, sigma1 = 0, 0.1 random0 = np.random.normal(mu1, sigma1, (int(dimy3/3), int(dimx3/3))) random2=skimage.transform.resize(random0, (dimy3, dimx3), anti_aliasing=True )

1. 1. tempdifference=tempdifference+random1+random2

bigtsa=skimage.transform.resize(tsa, (dimy3, dimx3), anti_aliasing=True )

dstsa=bigtsa+tempdifference

machinelearn=0

if(machinelearn==1): sla=[] sla.append(smalldem) sla.append(smalats) sla.append(smalons) slb=[] slb.append(dem4) slb.append(biglats) slb.append(biglons) apoints1=list(zip(*sla)) ## check this! bpoints1=tsa cpoints1=list(zip(*slb)) dows=linear_regression_multiple_vars(apoints1, bpoints1, cpoints1) #dows=gam_multiple_vars(apoints1, bpoints1, cpoints1) #dows=random_forest_multiple_vars(apoints1, bpoints1, cpoints1) #print (shape(dows)) dows2=np.reshape(dows, (dimy3, dimx3)) #dows3=dows2-np.min(dows2) #lapses=(dows3/dem4)*1000 #print (np.mean(lapses)) #plt.imshow(dows2) #plt.show() #quit(-1)

1. plt.imshow(mask1)
2. plt.imshow(bigtsa)
3. plt.imshow(bigsmalldem)
4. plt.imshow(bigtsa, cmap=kmap2)

1. plt.imshow(smalldem)
2. plt.imshow(bigsmalldem)

1. plt.imshow(dstsa, cmap=kmap2)

1. plt.imshow(tempdifference, cmap=kmap2)

1. plt.show()

1. print(mask1)

1. quit(-1)

avgtemp1=np.mean(tsa)

print("Tavg: ",avgtemp1)

1. quit(-1)

newx=64 newy=32

X = np.arange(-180, 180, 360/newx) Y = np.arange(-90, 90, 180/newy)

oname1="atm_temp.nc" ovar1="Band1"

savenetcdf_single_frommem(oname1, ovar1, tsa,Y, X)

X2 = np.arange(-180, 180, 360/dimx3) Y2 = np.arange(-90, 90, 180/dimy3)

oname1="tempds.nc" ovar1="Band1"

savenetcdf_single_frommem(oname1, ovar1, np.flipud(dstsa),Y2, X2)

1. contourrange1=[-50,-30,-20,-10,0,10,20,30,40,50]

contourrange1=[-50,-45,-40,-35,-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,40,45,50]

contourrange2=[-50,0,50]

contourrange3=np.arange(-10,40,1)

1. range3=np.arange(-50,50,0.1)

range3=np.arange(-10,40,0.1)

koloro1=['#0000FF', '#A0A0A0', '#FF0000']

imsiz1=(16,8)

fig,ax = plt.subplots(figsize=imsiz1)

size1=20

plt.suptitle(kaption1, fontsize=size1) ax.set_title(kaption2, fontsize=size1-3) size2=17

ax.set_xlabel("Longitude",fontsize=size2 ) ax.set_ylabel("Latitude", fontsize=size2 ) ax.tick_params(axis='x', labelsize= size2) ax.tick_params(axis='y', labelsize= size2)

tsa=np.flipud(tsa) dem1=np.flipud(dem1) mask1=np.flipud(mask1)

dstsa=np.flipud(dstsa) masco1=np.flipud(masco1)

1. plt.margins(0.0015, tight=True)

1. var1=ncinu1.variables['tso'][6]

1. plt.imshow(dstsa, extent=extent1,interpolation="bicubic", cmap="coolwarm")

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

1. selectvar=tsa

selectvar=dstsa

1. selectvar=tsa

ax.imshow( selectvar, origin="lower", extent=extent1, cmap=kmap2, vmin=-70, vmax=50 ,interpolation="bicubic")

cs = ax.contour(selectvar, origin="lower",extent=extent1, inline=True, colors='#7f0000', alpha=0.5, levels=contourrange1)

ax.clabel(cs, fontsize=20,fmt = '%3.0f', colors='#7f0000')

csdem2 = ax.contour(masco1, origin="lower",extent=extent1, linestyles='-',alpha=0.75, linewidths=[4], inline=True, colors=["#003f00"], levels=[-15000,0,15000])

fig.savefig(outname1)

print("plot")

plt.show()

print(".")