File:Ordovician ice age 445ma tas annual S1312 co2 120 tilt22 ecc 0p06 mvelp 90 1.png
![File:Ordovician ice age 445ma tas annual S1312 co2 120 tilt22 ecc 0p06 mvelp 90 1.png](https://upload.wikimedia.org/wikipedia/commons/thumb/0/0b/Ordovician_ice_age_445ma_tas_annual_S1312_co2_120_tilt22_ecc_0p06_mvelp_90_1.png/800px-Ordovician_ice_age_445ma_tas_annual_S1312_co2_120_tilt22_ecc_0p06_mvelp_90_1.png?20230606120642)
Original file (1,600 × 800 pixels, file size: 448 KB, MIME type: image/png)
Captions
Captions
Summary
[edit]DescriptionOrdovician ice age 445ma tas annual S1312 co2 120 tilt22 ecc 0p06 mvelp 90 1.png |
English: Mean temperature, Ordovician, 445 Ma if CO2 250 ppm, S=1318 obl 22 ecc 0.06 mvelp 90 |
Date | |
Source | Own work |
Author | Merikanto |
This image is based on Exoplasim and Paleodem
Simulation with exoplasim
https://github.com/alphaparrot/ExoPlaSim
Scotese paloedem downscaling
Scotese paloedem downscaling
PaleoDEM Resource – Scotese and Wright (2018)
11 August, 2018 by Sabin Zahirovic
Download PaleoDEM rasters as a zip file from here:
August 1, 2018 Dataset Open Access
PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic
Scotese, Christopher R; Wright, Nicky M
https://zenodo.org/record/5460860
https://zenodo.org/record/5460860/files/Scotese_Wright_2018_Maps_1-88_6minX6min_PaleoDEMS_nc.zip?download=1
For more information about this resource, contact Christopher R. Scotese at cscotese@gmail.com.
Exoplasim simulation code
- Exoplasim planet running code
- exoplasim example
- stepper code
-
- 01.06.2023 0000.0004b
-
- in ubuntu you must install
-
- pip3 install exoplasim[netCDF4]
- not
- "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)
- 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 ---")
- !!! Sun
- input_dem="dema1.nc" ##dem of exoplanet
- 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"
- 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
- a_baseflux=1365.2
- a_baseflux=1354
- a_baseflux=1318
a_baseflux=1312
a_yearlength=360
a_radius=1.0
a_gravity=9.81
- !!! ROTATION 7 !!!!
a_rotationperiod=0.883
- !!!!! earth nowadays
- a_eccentricity1=0.01671022
- a_obliquity1=23.4392811
- a_lonvernaleq1=102.94719
a_eccentricity1=0.06
a_obliquity1=22
a_lonvernaleq1=90
a_pCO21=120e-6
print("Exoplasim ...")
- 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(".")
PaleoDEM Resource – Scotese and Wright (2018)
11 August, 2018 by Sabin Zahirovic
Download PaleoDEM rasters as a zip file from here:
August 1, 2018 Dataset Open Access
PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic
Scotese, Christopher R; Wright, Nicky M
https://zenodo.org/record/5460860
https://zenodo.org/record/5460860/files/Scotese_Wright_2018_Maps_1-88_6minX6min_PaleoDEMS_nc.zip?download=1
For more information about this resource, contact Christopher R. Scotese at cscotese@gmail.com.
Exoplasim code to produce data
File:Late_ordovician_tas_445ma_co2_250_s_1318_o_22_ecc006_mvelp_90_1.png#%7B%7Bint%3Afiledesc%7D%7D
Python3 code to plot image
-
- temperature plotter
- python 3
- 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
- from skimage.transform import resize
import skimage
- 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)
- 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)
- 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)
- main
- kaption1="Temperature of July 67 million years ago, degrees Celsius"
kaption1="Mean temperature of year degrees Celsius, final Ordovician 445 Ma "
- kaption2="S=1312 CO2=180 ppm obliq=22 ecc=0.06 mvelp=90"
- 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"
- outname1="cretaceous_67ma_co2_4x_tjanuary_1.png"
year1=0
month1=0 ## if 0, mean surf temp of year
- month1=1
vari1="tas"
tsa_offset=0
- fintemp1="./exoplasim_445ma_co2_125/MOST.00040.nc"
- fintemp1="./exoplasim_445ma_co2_180/MOST.00050.nc"
- fintemp1="./exo_445_omega45_co2_150/MOST.00030.nc"
- fintemp1="./exo_445_omega270_co2_120/MOST.00024.nc" ## nok
fintemp1="./exo_445_omega90_co2_120/MOST.00020.nc"
- findem1="Map22_PALEOMAP_6min_Mid-Cretaceous_95Ma.nc"
findem1="./maps/Map77_PALEOMAP_6min_Late_Ordovician_445Ma.nc"
- downscaleing size
dimy3=360
dimx3=720
extent1=(-180,180,-90,90)
- extent1=(-180,180,90,-90)
- print(num1, kk1)
- print(fintemp1)
kindex1=1
ncinu1 = Dataset(fintemp1, "r+", format = "NETCDF4")
print(ncinu1)
- quit(-1)
ncdem1 = Dataset(findem1, "r+", format = "NETCDF4")
- print(ncinu1)
- quit(-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])
- snd = np.array(ncinu1.variables['snd'])
tsa=np.mean(tsa00, axis=0)
- if(month1==0): tsa=np.mean(tsa, axis=0)
- 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)
- plt.imshow(maxsnowdepth)
- plt.contour(minsnowdepth, levels=[0])
- plt.imshow(mask00, cmap="terrain")
- plt.contour(minsnowdepth, levels=[0])
- plt.imshow(tsa_january, cmap="jet")
- plt.imshow(tsa_july, cmap="jet")
- plt.imshow(tsa, cmap="jet")
- plt.show()
- quit(-1)
dem00=np.array(ncdem1.variables['z'])
dem01=np.copy(dem00)
dem1=skimage.transform.resize(dem01,(dimy3, dimx3))
- warning
- dem1=np.flipud(dem1)
- plt.imshow(np.array(ncinu1.variables['tso'][1]), extent=extent1)
- plt.imshow(tsa_january, cmap="terrain")
- plt.imshow(dem1, cmap="terrain")
- plt.show()
- quit(-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
- plt.imshow(smalons)
- plt.imshow(smalats)
- plt.imshow(biglons)
- plt.imshow(biglats)
- plt.show()
- resize input dem
- print (len(dem00))
kmap2 = plt.cm.get_cmap('RdBu_r')
kmap1=plt.get_cmap('YlGnBu')
dmap1=plt.get_cmap('binary')
- kmap2 = plt.cm.get_cmap('bwr')
kmap2 = plt.cm.get_cmap('coolwarm') ## ok
- kmap2 = plt.cm.get_cmap('jet')
- kmap2 = plt.cm.get_cmap('turbo') ##ok
- kmap2 = plt.cm.get_cmap('gist_rainbow_r')
lon1=np.array(lon0)
lat1=np.array(lat0)
- print(lon1)
- tsa=np.array(tsa00)
- mask1=np.array(mask00) ## ok
- print(np.shape(tsa[month1]))
- tsa=tsa+tsa_offset
- plt.imshow(tsa, cmap=kmap2, interpolation="bicubic")
- plt.imshow(dem00, cmap=kmap2, interpolation="bicubic")
- plt.show()
- quit(-1)
dem2=np.copy(dem1)
dem2[dem2>-1]=1
dem2[dem2<0]=0
- sla=[]
- sla.append(tsa)
- apoints1=
- bpoints1=list(zip(*sla))
dem4=np.copy(dem1)
dem4[dem4<0]=0
masco1=np.copy(dem1)
masco1[masco1<0]=0
masco1[masco1>0]=255
- plt.imshow(masco1)
- plt.imshow(mask1)
- plt.imshow(tsa)
- plt.imshow(dem4)
- plt.show()
- 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
- 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 )
- 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)
- plt.imshow(mask1)
- plt.imshow(bigtsa)
- plt.imshow(bigsmalldem)
- plt.imshow(bigtsa, cmap=kmap2)
- plt.imshow(smalldem)
- plt.imshow(bigsmalldem)
- plt.imshow(dstsa, cmap=kmap2)
- plt.imshow(tempdifference, cmap=kmap2)
- plt.show()
- print(mask1)
- quit(-1)
avgtemp1=np.mean(tsa)
print("Tavg: ",avgtemp1)
- 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)
- 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)
- 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)
- plt.margins(0.0015, tight=True)
- var1=ncinu1.variables['tso'][6]
- plt.imshow(dstsa, extent=extent1,interpolation="bicubic", cmap="coolwarm")
- plt.show()
- quit()
- selectvar=tsa
selectvar=dstsa
- 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(".")
Licensing
[edit]![w:en:Creative Commons](https://upload.wikimedia.org/wikipedia/commons/thumb/7/79/CC_some_rights_reserved.svg/90px-CC_some_rights_reserved.svg.png)
![attribution](https://upload.wikimedia.org/wikipedia/commons/thumb/1/11/Cc-by_new_white.svg/24px-Cc-by_new_white.svg.png)
![share alike](https://upload.wikimedia.org/wikipedia/commons/thumb/d/df/Cc-sa_white.svg/24px-Cc-sa_white.svg.png)
- 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 | 12:06, 6 June 2023 | ![]() | 1,600 × 800 (448 KB) | Merikanto (talk | contribs) | update |
06:50, 3 June 2023 | ![]() | 1,600 × 800 (330 KB) | Merikanto (talk | contribs) | Update of data and layout | |
14:45, 1 June 2023 | ![]() | 1,600 × 800 (861 KB) | Merikanto (talk | contribs) | Change of layout | |
16:45, 31 May 2023 | ![]() | 1,600 × 800 (958 KB) | Merikanto (talk | contribs) | update | |
16:42, 31 May 2023 | ![]() | 1,600 × 800 (1.16 MB) | Merikanto (talk | contribs) | Update | |
18:41, 29 May 2023 | ![]() | 1,600 × 800 (1.21 MB) | Merikanto (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
The following page uses this file:
- File:Ordovician 445ma tas 1.png (file redirect)
Metadata
This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Software used | |
---|---|
Horizontal resolution | 39.37 dpc |
Vertical resolution | 39.37 dpc |