File:World 200ma mean temperature of year 2.png

author       = {Scotese, Christopher R and
                Wright, Nicky M},
title        = {{PALEOMAP Paleodigital Elevation Models (PaleoDEMS) 
                 for the Phanerozoic}},
month        = aug,
year         = 2018,
publisher    = {Zenodo},
doi          = {10.5281/zenodo.5460860},
url          = {https://doi.org/10.5281/zenodo.5460860}

}

https://zenodo.org/record/5460860/files/Scotese_Wright_2018_Maps_1-88_6minX6min_PaleoDEMS_nc.zip?download=1}}

1. 1. take 16-bit dem for koppenpasta
   2. uses gdal, imagemagick

1. inputdem1=$1

inputdem1="./indata/Map43_PALEOMAP_6min_Late_Triassic_200Ma.nc"

demdimx=360 demdimy=180

1. demdimx=1600
2. demdimy=800

rmdir origo mkdir origo gdal_translate -of GTiff $inputdem1 ./origo/dem0.tif gdal_calc.py -A ./origo/dem0.tif --outfile=./origo/dem6.tif --calc="A*(A>0)"

gdal_translate -of NetCDF ./origo/dem6.tif ./origo/dem6.nc

gdalwarp -of GTiff -overwrite -ts $demdimx $demdimy ./origo/dem6.tif ./origo/dem1600.tif

gdal_translate -of NetCDF ./origo/dem1600.tif ./origo/dem1600.nc gdal_translate -of png -ot Uint16 ./origo/dem1600.tif ./origo/dem1600.png

gdal_translate -ot Int16 -of GTiff ./origo/dem1600.tif ./origo/dem16.tif

convert ./origo/dem16.tif -compose Copy_Opacity -auto-level ./origo/koppendem0.tif

convert ./origo/koppendem0.tif ./origo/koppendem0.png

gdal_translate -of png -ot Uint16 ./origo/koppendem0.tif ./origo/koppendem.tif

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3
2. 1. python3 calculate mannual mean temperature from exoplasim output
4. NOTE:
5. uses average of layers, assumes uuniform p sttructure in atmopshere
7. 30.6.2022 v 0000.0000
9. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3

import netCDF4 as nc import numpy as np from scipy import interpolate from scipy.interpolate import griddata

import matplotlib.pyplot as plt from matplotlib.pylab import * import matplotlib.mlab as mlab

1. 1. std atmosphere util libs

from ambiance import Atmosphere import metpy.calc as mpcalc from metpy.units import units

1. import skimage

from skimage.transform import resize from scipy import misc

def ncread_float_array(fn1, var1): ds2 = nc.Dataset(fn1) varr1=ds2[var1][:] return(varr1.astype(float))

def ncread(fn1, var1): ds2 = nc.Dataset(fn1) varr1=ds2[var1] return(varr1)

def ncread_12months_avg(fn1, var1): ds2 = nc.Dataset(fn1) varr0=ds2[var1][:] varr1=np.array(np.average(varr0,axis=0)).astype(float) return(varr1)

def savenetcdf_single_frommem(outfilename1, outvarname1, xoutvalue1,xoutlats1,xoutlons1): nlat1=len(xoutlats1) nlon1=len(xoutlons1) #indata_set1=indata1 print(outfilename1) ncout1 = nc.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 ncsave_T21(oname1, ovar1, ovals1): newx=64 newy=32 X = np.arange(-180, 180, 360/newx) Y = np.arange(-90, 90, 180/newy) savenetcdf_single_frommem(oname1, ovar1,ovals1,Y, X)

def ncsave_360(oname1, ovar1, ovals1): newx=360 newy=180 X = np.arange(-180, 180, 360/newx) Y = np.arange(-90, 90, 180/newy) savenetcdf_single_frommem(oname1, ovar1,ovals1,Y, X)

def read_ta_levels_m(inf1): #ts=ncread(inf1, "ta") ds = nc.Dataset(inf1) #ta=ds["ta"] #print(ta) #tim1=ds.dimensions['lev'] #print(tim1) #lev=ds["lev"] lev=ds["lev"][:] #print(lev[7]) #print (lev) sealevel = Atmosphere(0) p0=sealevel.pressure #print (p0) pressures=p0*lev #print(pressures) #levms=metpy.calc.pressure_to_height_std(1000) #press = np.array([1000., 500.]) * units.hPa press=np.array(pressures)*units.Pa levms=mpcalc.pressure_to_height_std(press) levms=np.array(levms[:]).astype(float)*1000.0 #print(levms) return(levms)

def ncread_avg_annual_levels(fn1, varname1): ds2 = nc.Dataset(fn1) var0=np.array(ds2[varname1][:])-273.15 #var0s=ds2[varname1] #print(var0s) #quit(-1) var1=var0.mean(axis=0) #print(var1, np.shape(var1)) return(var1)

def calculate_temperature_of_point(inx,iny,levstak, tempstak, dema1):

height1=dema1[iny,inx]

dt65=height1*-6.5 temp65=tempstak[9, iny,inx]+dt65

#print(levstak) #print(np.shape(tempstak))

#quit(-1)

#plt.imshow(tempstak[6]) ##plt.imshow(dema1) #plt.show() #print(levstak) #print(inx,iny,height1) idx2=np.argwhere(np.array(levstak)<=height1)[0] idx1=idx2-1 #print(idx1,idx2)

a=levstak[idx1] ## upper b=levstak[idx2] ## lower

te1=tempstak[idx1, iny,inx] te2=tempstak[idx2, iny,inx]

difab=a-b dift=te1-te2

ko=dift/difab

indh=a-height1

dt1=ko*indh

tempp1=te1-dt1

#print(te1)

return(tempp1)

infilename1="./indata/input.nc"

kaption1="Mean temperature °C, 200 Ma" demname1="./origo/dem1600.nc" savename1="./output/mean_temperature_250ma_1.png"

dimx, dimy = 360, 180

1. 1. read dem

dem1=np.flipud(ncread_float_array(demname1, "Band1"))

1. 1. read average levels, assume std atmosphere

levm=read_ta_levels_m(infilename1)

temp2m=ncread_12months_avg(infilename1, "tas")-273.15 temp0m=ncread_12months_avg(infilename1, "tsa")-273.15

templevs1=ncread_avg_annual_levels(infilename1, "ta")

levm2=np.copy(levm) templevs2=np.copy(templevs1)

levm2=np.append(levm2,2.0) levm2=np.append(levm2,0.0)

at2=np.atleast_3d(temp2m).reshape(1,32,64) at3=np.atleast_3d(temp2m).reshape(1,32,64) print (np.shape(at2))

templevs2=np.append(templevs2, at2, axis=0) templevs2=np.append(templevs2, at3, axis=0)

redemps=[]

for n in range(0,10): print(n) temp0=templevs2[n] temp1=resize(temp0,(dimy, dimx),order=3) redemps.append(temp1) #print ("N ",n, " ... ", levm2[n])

levm3=np.array(levm2).astype(float) redemps3=np.array(redemps).astype(float)

nn=9

temp1=redemps[nn] h1=levm2[nn]

print(h1)

resultemp1=np.empty((dimy, dimx))

print (shape(resultemp1))

1. plt.imshow(resultemp1)

1. plt.show()

1. inx=180
2. iny=90

for iy in range(0,(dimy-1)): for ix in range(0,(dimx-1)): tem1=calculate_temperature_of_point(ix,iy,levm3, redemps3, dem1) resultemp1[iy,ix]=tem1

ncsave_360("temp_annual_dskaledx1.nc", "tas", np.flipud(resultemp1))

plt.imshow(resultemp1)

plt.show()

quit(-1)