File:Fictional exoplanet surface temperature c mean annual 1.png

1. 1. "R" script bmp exoplanet dem to netcdf georaster
   2. uses pixmap, bmp
2. 1. 28.07.2022 0000.0005

1. install.packages("bmp")

library(raster) library(ncdf4) library(rgdal) library(pixmap) library(png) library(bmp) library(stringr)

save_as_nc<-function(outname1,r, outvar1, longvar1, unit1) { ext2<-c(-180, 180, -90,90) extent(r)<-ext2 crs(r)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(r, outname1, overwrite=TRUE, format="CDF", varname=outvar1, varunit=unit1,

       longname=longvar1, xname="lon",   yname="lat")

}

bmp_dem_to_nc<-function(inbmpname1, outncname1, deltaz1, sealevel1) { #bm1=read.bmp(inbmpname1) print(inbmpname1) bm1=readPNG(inbmpname1, native = FALSE, info = FALSE) str(bm1) plot(bm1) #stop(-1)

pr1=pixmapGrey(bm1) data00<-pr1@grey data01<-data00 dim1=dim(data01) dimx=dim1[1] dimy=dim1[2] print (dimx) print (dimy) rin1 <- raster(data01) ext1<-c(-180, 180, -90,90) extent(rin1)<-ext1 min1=minValue(rin1) max1=maxValue(rin1) delta1=max1-min1 rin2=(rin1-min1)/delta1 rin3=(rin1*deltaz1)-sealevel1 #image(rin3) rin3a<-flip(rin3, direction="y")

save_as_nc(outncname1,rin3, "z", "z", "m")

1. rin4=flip(rin3)
2. outncname2=paste0(outncname1,"sand_flipped.nc")

1. save_as_nc(outncname2,rin4, "z", "z", "m")

return(rin3) }

process_base_rasters<-function(infile1, file2, file3, masksealevel1) { ## create absolute sealevel dem and lon, lat ras

ur1<-raster(infile1)

ur1<-ur1-masksealevel1

ur1[ur1[]<1] <- 0

#ur1<-flip(ur1)

lonr1 <- init(ur1, 'x')

latr1 <- init(ur1, 'y')

crs(ur1)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" writeRaster(ur1, file2, overwrite=TRUE, format="CDF", varname="Band1", varunit="m",

       longname="Band1", xname="lon",   yname="lat")

writeRaster(ur1, file3, overwrite=TRUE, format="GTiff", varname="Band1", varunit="m",

       longname="Band1", xname="lon",   yname="lat")
       
  crs(lonr1)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
  crs(latr1)<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
 

1. writeRaster(lonr1, "./origo/lons.nc", overwrite=TRUE, format="CDF", varname="Band1", varunit="deg",
2. longname="Band1", xname="lon", yname="lat")
3. writeRaster(latr1, "./origo/lats.nc", overwrite=TRUE, format="CDF", varname="Band1", varunit="deg",
4. longname="Band1", xname="lon", yname="lat")

}

create_mask_rasters<-function(infile1,maskname1,demname1) {

ur1<-raster(infile1) ur1[ur1[]<1] <- 0 ## jn warning flip! #rdem1=flip(ur1) r=flip(ur1) #r=ur1 rdem1=ur1

mini=minValue(rdem1) maxi=maxValue(rdem1) delta=maxi-mini rdem2=(rdem1/delta)*255

dims<-dim(r)

r[r[]<1] <- 0 r[r[]>0] <- 1

zeros1<-sum(r[]==0) ones1<-sum(r[]==1)

print(zeros1) print(ones1) stop(-1)

image(r)

print (dims[1]) print (dims[2])

rows=dims[2] cols=dims[1]

mask0<-r mask1<-mask0[]

idem1<-rdem2[] mask2<-matrix(mask1, ncol=cols, nrow=rows ) idem2<-matrix(idem1, ncol=cols, nrow=rows ) #mask3<-t(mask2) #idem3<-t(idem2) mask3<-mask2 idem3<-idem2 r <- writePNG(mask3, maskname1) r <- writePNG(idem3, demname1) }

createmask <- function(dem1, demname2, maskname1, maskname2,maskname3, maskname4, masksealevel1) {

 dem2<-(dem1+abs(masksealevel1))
 
 mask2<-dem2
 mask1<-dem2
 
 mask2[mask2 > 1 ]=255
 mask2[mask2 < 0 ]=0 
 
 mask1[mask1 > -1 ]=0
 mask1[mask1 < 0 ]=255  
 
 	zeros1<-sum(mask2[]==0)

ones1<-sum(mask2[]==255) all1=zeros1+ones1

landsearatio1=(100*zeros1)/all1

print(zeros1) print(ones1) print (all1) print(800*1600) print("Sea per cent of total surface:") print(landsearatio1)

#print("TEBUK STOP !!!") #stop(-1)

 image(dem1)
 image(mask1)
 image(mask2)
   
 #dem2<-flip(dem2, direction= "y")
 #mask1<-flip(mask1, direction="y") 
 #mask2<-flip(mask2, direction="y")    
 
 ext1<- extent(0, 360, -90, 90)
 ext2<- extent(-180, 180, -90, 90) 
 
 mask30<-mask1
 extent(mask30) <- ext1
 mask3 <- rotate(mask30)
 extent(mask3) <- ext1
 mask40<-mask2
 extent(mask40) <- ext1
 mask4 <- rotate(mask40)
 extent(mask4) <- ext1
 
 save_as_nc(demname2,dem2, "Band1", "Band1", "")
 save_as_nc(maskname1,mask1, "Band1", "Band1", "")
 save_as_nc(maskname2,mask2, "Band1", "Band1", "")
 save_as_nc(maskname3,mask3, "Band1", "Band1", "")
 save_as_nc(maskname4,mask4, "Band1", "Band1", "")  
 

}

1. inbmpname1="./indata/dryplanet1.bmp" ## bmp planet dem map
2. inbmpname1="./indata/Map-111.bmp" ## bmp planet dem map
3. inbmpname1="./indata/Map-111.png" ## bmp planet dem map

1. 1. here
2. https://topps.diku.dk/torbenm/maps.msp
3. seed 1111111, rectangular, greyscale
   1. reso 6400x3200

inbmpname1="./indata/inmap.png" ## bmp planet dem map

outncname1="./indata/dryplanet1.nc"

1. deltaz1=19600

1. 1. near same heights than in earth

deltaz1=11200+8400 ## difference of output min, max height

1. 1. sealevel, that used to create dem

sealevel1=11200 #desert planet, sealevel 0

1. 1. sealevel, that used in landsea mask

1. 1. "real sea level"

masksealevel1=-8000.0

file2="./indata/dryplanet_dem.nc" file3="./indata/dryplanet_dem.tif"

maskname1= "./indata/dryplanet_mask_land.nc" maskname2= "./indata/dryplanet_mask_sea.nc" maskname3= "./indata/dryplanet_mask_land_rotated.nc" maskname4= "./indata/dryplanet_mask_sea_rotated.nc"

demname1= "./indata/dryplanet_dem.nc" demname2= "./indata/dryplanet_dem_sealevel.nc"

dem1<-bmp_dem_to_nc(inbmpname1, outncname1, deltaz1, sealevel1)

createmask(dem1,demname2,maskname1, maskname2, maskname3, maskname4, masksealevel1)

1. process_base_rasters(outncname1, file2, file3, masksealevel1)
2. create_mask_rasters(file2,maskname1,demname1)

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

import numpy as np import math as math import matplotlib.pyplot as plt from scipy.interpolate import interp2d import netCDF4

import exoplasim as exo

NLAT=0 NLON=0

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 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['lat'][:] inlons=nc['lon'][:] #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. exoplasim ,,,

def run_exoplasim_wx(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):

#output_format=".npz" output_format=".nc"

a_pO2=1-a_pCO2-0.79 a_pN2=(1-0.21-a_pCO2)

print("Precess 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)

create_sras(topo, a_seamasklevel1)

print("Creating exoplasim object ")

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

## earth 21000 BP glaciers1= { "toggle": True, "mindepth":2, "initialh":-1 }

testplanet.configure( startemp=5772.0, flux=1367,# Stellar parameters eccentricity=a_eccentricity, obliquity=a_obliquity, lonvernaleq=a_lonvernaleq, year=365, fixedorbit=True, # Orbital parameters rotationperiod=1, # Rotation topomap="topo.sra", landmap="landmask.sra", radius=1.0, gravity=9.80665, # Bulk properties #seaice=False, #maxsnow=False, #glaciers=False, #stormclim=False, #vegetation=0, wetsoil=True, #alters albedo of soil based on how wet it is

               vegetation=2,                               #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

testplanet.exportcfg()

print("Running ExoPlasim ... ")

#testplanet.run(years=a_years,crashifbroken=True)

testplanet.run(years=a_years,crashifbroken=True)

return(0)

1. 1. stepper exoplasim ...

def run_exoplasim_wy(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):

print("Exoplasim runner")

a_startemp=5772.0 a_baseflux=1367 a_yearlength=365 a_radius=1.0 a_gravity=9.80665 a_rotationperiod=1

a_pO2=1-a_pCO2-0.79 a_pN2=(1-0.21-a_pCO2)

output_format=".nc"

print("Precess 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)

create_sras(topo, a_seamasklevel1)

print("Creating exoplasim object ")

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

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

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

               vegetation=2,                               #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

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") tas2=np.ravel(tas) tmean1=np.mean(tas2) tmean2=math.floor( (tmean1-273.15)*100) tmean3=tmean2/100 print("Mean tas "+str(tmean3) )

print("Return.")

return(0)

print(" Exoplasim simulation code ---")

input_dem="./indata/dryplanet1.nc" ##dem of exoplanet

a_modelname1="planet" a_workdir1="planet_run"

a_runsteps1=300 a_years1=a_runsteps1 a_timestep1=30 a_snapshots1=0 a_ncpus1=4 a_layers1=8 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=-8000

1. 1. earth nowadays

a_eccentricity1=0.01671022 a_obliquity1=23.44 a_lonvernaleq1=102.7 a_pCO21=360e-6

print("Exoplasim ...")

1. 1. attempt to run exoplasim stepper code

run_exoplasim_wy(input_dem, a_seamasklevel1, a_resolution1, a_layers1, a_years1,a_timestep1,a_snapshots1,a_ncpus1,a_eccentricity1,a_obliquity1,a_lonvernaleq1,a_pCO21)

1. run_exoplasim_w(input_dem, a_seamasklevel1, a_resolution1, a_layers1, a_years1,a_timestep1,a_snapshots1,a_ncpus1,a_eccentricity1,a_obliquity1,a_lonvernaleq1,a_pCO21)

print(".")