File:Gliese 12b example simulation sea land 2 r 1 1 1 1.png

https://github.com/alphaparrot/ExoPlaSim

Planet map generator in "diku" mapm 77777, adjusted curves with image processor: mone lower areas for sea

https://topps.diku.dk/torbenm/maps.msp

Nasa Panoply

https://www.giss.nasa.gov/tools/panoply/

2. 1. exoplasim planet runner
   2. python 3, exoplasim
4. basic run settings T21
5. 1. 29.10.2023 0000.0002a2

import math import numpy as np import exoplasim as exo

timestep=30.0 years=100

Sol=1361.5

1. 1. basic input params

startemp=3296 # assumption Earth-like

1. 1. Gliese 12 b

planetname="Gliese 12 b" S1=1.63 flux=S1*Sol year=12.76140 rotationperiod=year radius=0.958 ## rel to earth mass=math.pow(radius, 3.7)

eccentricity=0 obliquity=0

fixedorbit=True synchronous=True

1. substellarlon=0

substellarlon=180 aquaplanet=False desertplanet=False vegetation=2 ## 0 no 1 static 2 dynamic stormclim=False aerosol=False co2weathering=False outgassing=0 evolveco2=False maxsnow=-1 ## False ? ozone=False vegaccel=1 seaice=False wetsoil=False

glaciers1= { "toggle": True, "mindepth":2, "initialh":-1 }

1. landmap="Alderaan_surf_0129.sra"
2. topomap="Alderaan_surf_0172.sra"

1. landmap=None
2. topomap=None

landmap="mapa_surf_0129.sra" topomap="mapa_surf_0172.sra"

1. relCO2=0.95
2. relN2=1-relCO2

relCO2=280e-6 relO2=0.21 relN2=1-relCO2-relO2

1. relCO2=0.95
2. relO2=0.045
3. relN2=1-relCO2-relO2

density=mass/np.power(radius,3)*5.519 vesc=np.sqrt(mass/radius)*11.186 geese=mass/(radius*radius) geeseg=geese*9.80665

patm1=geese*mass*radius*radius ## estimation patm2=math.pow(radius, 2.4)

gravity=geeseg

pressure=(patm1+patm2)/2

1. pressure=10
2. pressure=0.1

pN2=pressure*relN2 pCO2=pressure*relCO2 pO2=pressure*relO2

print(" Name ", planetname) print(" insol S0 ",round(S1,3) ) print(" insol Wm2 ",round(flux,2) ) print(" mass ",round(mass,2) ) print(" radius ",round(radius,2) ) print(" radius km ",round((radius*6371),2) ) print(" density ",round(density,2) ) print(" vesc ",round(vesc,2) ) print(" gee_e ",round(geese,2) ) print(" gee_ms2 ",round((geese*9.81),2) ) print(" Patm ",round(pressure,6) )

print("Running wait long time ...")

1. quit(-1)

planeta = exo.Model(workdir="planeta_run",modelname="Planeta",ncpus=4,resolution="T21",outputtype=".nc")

planeta.configure(year=year, wetsoil=wetsoil,pO2=pO2, ozone=ozone, vegaccel=vegaccel, seaice=seaice, glaciers=glaciers1, maxsnow=maxsnow, evolveco2=evolveco2, outgassing=outgassing, co2weathering=co2weathering, vegetation=vegetation, stormclim=stormclim, aerosol=aerosol,landmap=landmap, topomap=topomap,startemp=startemp, flux=flux, eccentricity=eccentricity,obliquity=obliquity,fixedorbit=fixedorbit,synchronous=synchronous,substellarlon=substellarlon,rotationperiod=rotationperiod,radius=radius,gravity=gravity,aquaplanet=aquaplanet,desertplanet=desertplanet,pN2=pN2,pCO2=pCO2,timestep=timestep,snapshots=False,physicsfilter="gp|exp|sp")

print(" Run, wait a long time...") planeta.run(years=years,crashifbroken=True)

planeta.exportcfg()

planeta.run(years=years,crashifbroken=True)

planeta.finalize("Planeta",allyears=True,keeprestarts=True)

planeta.save()

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3
2. 1. diku "planet map generator" map converter to exoplasim input test
3. maps from https://topps.diku.dk/torbenm/maps.msp
   1. must use here grayscale palette output option !

1. 1. normalizes map to 0...1, then np.exp(z), then normalize to depest and highest points
   2. python3 source code
2. 1. 24.6.2024 v 0000.0000c2

import os import math as math import numpy as np

from scipy.interpolate import interp2d import matplotlib.pyplot as plt

from PIL import Image, ImageFilter

from imageio import imwrite import imageio

import netCDF4

1. import exoplasim as exo

1. 1. main ctrls

iname1="Map-777.bmp"

top1=7000 floor1=-8000

land_height_scaling_exponent_1=15 ## land height scaling exponent

sealevel1=0.001 ## if land scaling exponent hight mush use here small value! level in original bittmap 0...1

def t21_meanz(inmap):

   latitudes0=np.array([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])

   radlats1=np.radians(latitudes0)
   coslats1=np.cos(radlats1)
   height1=len(radlats1)
   seas1=0
   all1=0
   landmeanz1=0

   for iy in range(0,(height1)):
       line1=inmap[iy,:]
       lek1=abs(coslats1[iy])
       seadd0=np.count_nonzero(line1 == 0)
       #print(seadd0)
       seadd1=seadd0/len(line1)
       seadd2=lek1*seadd1
       seas1=seas1+seadd2
       all1=all1+lek1
       lamez0=line1[line1>0].mean()
       if (math.isnan(lamez0)): lamez0=0
       lamez1=lamez0*lek1
       landmeanz1=landmeanz1+lamez1

   landmeanz2=round((landmeanz1/height1),2)

   seapercent2=round(((seas1*100)/all1),2)
   landpercent2=100-seapercent2
   #print(seas1)
   #print(all1)
   #print(seapercent2)
   #print(landpercent2)
   return(landpercent2, seapercent2, landmeanz2)

def count_spherical_land_sea_meanz(inmap):

   radlats1=np.linspace(-np.pi/2,np.pi, height1)
   coslats1=np.cos(radlats1)
   seas1=0
   all1=0
   landmeanz1=0

   for iy in range(0,(height1)):
       line1=inmap[iy,:]
       lek1=abs(coslats1[iy])
       seadd0=np.count_nonzero(line1 == 0)
       #print(seadd0)
       seadd1=seadd0/len(line1)
       seadd2=lek1*seadd1
       seas1=seas1+seadd2
       all1=all1+lek1
       lamez0=line1[line1>0].mean()
       if (math.isnan(lamez0)): lamez0=0
       lamez1=lamez0*lek1
       landmeanz1=landmeanz1+lamez1

   landmeanz2=round((landmeanz1/height1),2)

   seapercent2=round(((seas1*100)/all1),2)
   landpercent2=100-seapercent2
   #print(seas1)
   #print(all1)
   #print(seapercent2)
   #print(landpercent2)
   return(landpercent2, seapercent2, landmeanz2)

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="lat" inlonname1="lon" 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 conv_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['lat'][:] inlons=nc['lon'][:] #print(inlats) #print(inlons) indimx=len(inlons) indimy=len(inlats)

latlen=len(inlats) lonlen=len(inlons)

#print(lonlen, latlen)

indimx=lonlen indimy=latlen

dem000=nc['z'] #dem=np.flipud(dem000) dem=np.copy(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) landpercent2, seapercent2, landmeanz2=t21_meanz(topo2) print("") print("Output T21 dem \n Land % ", landpercent2, "\n Mean land z ", landmeanz2) #savenetcdf_single_frommem("maskt21.nc", "z", masko,lts,lns) writeSRA("mapa",129,np.flipud(topo2),NLAT,NLON) writeSRA("mapa",172,np.flipud(masko),NLAT,NLON) return(topo,lons,lats)

def savenetcdf2(outfilename1, outvarname1, map1):

shape1=np.shape(map1) width1=shape1[1] height1=shape1[0] xoutlons1=np.linspace(0,360,width1) xoutlats1=np.linspace(-90,90, height1)

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[:, :] =np.flipud(map1[:]) ncout1.close() return 0

1. 1. main program

NLAT=0 NLON=0

imin1=Image.open(iname1)

img1 = np.asarray(imin1).astype(float)

max1=np.max(img1) min1=np.min(img1) del1=max1-min1

map0=(img1-min1)/del1

1. map1=np.exp(map0)/np.exp(1)

map1a=np.copy(map0)

1. 1. land_height_scaling_exponent_1=15

map1a=np.power(map1a, land_height_scaling_exponent_1)

max2=np.max(map1a) min2=np.min(map1a) del2=max2-min2

map1=(map1a-min2)/del2

map1=map1[:,:,0]

shape1=np.shape(map1)

width1=shape1[1] height1=shape1[0]

imsize1=width1*height1

1. quit(-1)

map2=np.copy(map1) map3=np.copy(map1)

map2=np.where(map2<sealevel1,0,map2)

map2=np.where(map2<sealevel1,map2*floor1,map2*top1)

1. 1. abyss mat too ...

map3=np.where(map3<sealevel1,map3*floor1,map3*top1)

landpixels1=np.count_nonzero(map2)

landpercent1=round( ((landpixels1*100)/imsize1),1)

1. quit(-1)

landpercent2, seapercent2, landmeanz2=count_spherical_land_sea_meanz(map2)

1. map2a=(np.abs(map2)*256*256).astype(np.uint16)

map2a=map2

1. print (np.max(map2a))

1. quit(-1)

1. print(map2a)

1. quit(-1)

imout1 = Image.fromarray(map2a)

imageio.imwrite('map2.png', map2a.astype(np.uint16))

1. imout1.save("map1.png")

1. quit(-1)

map3=np.roll(map2a, int(width1/2), axis=1)

plt.imshow(map3) plt.show()

1. imout2 = Image.fromarray(map3)

imageio.imwrite('map1.png', map3.astype(np.uint16))

mask2=np.copy(map2a) mask3=np.copy(map3)

mask2=np.where(mask2<1,0, 65535) mask3=np.where(mask3<1,0, 65535)

1. imout3 = Image.fromarray(mask2)
2. imout4 = Image.fromarray(mask3)

imageio.imwrite('mask2.png', mask2.astype(np.uint16))

imageio.imwrite('mask1.png', mask3.astype(np.uint16))

image1 = Image.open("mask1.png")

image1 = image1.convert("L")

image1 = image1.filter(ImageFilter.FIND_EDGES)

image1.save("edge1.png")

image2 = Image.open("mask2.png")

image2 = image2.convert("L")

image2 = image2.filter(ImageFilter.FIND_EDGES)

image2.save("edge2.png")

savenetcdf2("map1.nc", "z", map2a) print() print() print() print("Output coarse map params:") conv_t21("map1.nc", "inmap", 0)

1. print(width1, height1)
2. print(" Land %" , landpercent1)

print() print("Accurate dem") print(" Land  % ", landpercent2) print(" Sea  % ", seapercent2)

print(" Land mean height ", landmeanz2)

print(".")