File:Earth like planet habitability as relation of distance to sun like central star 1.png

3. earth-like planet habitability vs distance to solar twin
   1. plot exoplasim data from different simulations
4. 1. 28.06.2024 0000.0000.0001

import xarray as xr import numpy as np import math import matplotlib.pyplot as plt

iname1="./indata2/land30_94au.nc" iname2="./indata2/land30_97au.nc" iname3="./indata2/land30_1au.nc" iname4="./indata2/land30_103au.nc" iname5="./indata2/land30_106au.nc"

def plot_temperature_au(aus2, temp1):

   #hab1=np.nan_to_num(hab1,0)
   plt.plot(aus2, temp2, lw=10, color="#1f1fff")
   plt.title("Planet distance vs habitability", size=18)
   plt.xlabel("Distance from Sun-like G2V star AU", size=16)
   plt.ylabel("Habitability estimation", size=16)
   plt.xticks(fontsize=16)
   plt.yticks(fontsize=16)
   plt.grid()
   #plt.ylim([0,100])
   plt.legend(fontsize=16)
   #plt.plot(fox2)
   plt.show()

def plot_habitability_au(aus2, hab2):

   hab1=np.nan_to_num(hab2,0)
   plt.plot(aus2, hab2, lw=6, color="#1f1fff")
   plt.title("Earth-like planet distance - habitability", size=18)
   plt.xlabel("Distance from Sun-like G2V star au", size=16)
   plt.ylabel("Habitability estimation", size=16)
   plt.xticks(fontsize=16)
   plt.yticks(fontsize=16)
   plt.grid()
   #plt.ylim([0,100])
   plt.legend(fontsize=16)
   #plt.plot(fox2)
   plt.show()

def habitability_f_T_RH(T, RH):

   # https://phl.upr.edu/projects/standard-primary-habitability-sph
   #https://phl.upr.edu/projects/standard-primary-habitability-sph
   # Surface Temperature (°C) 	-8.3 	24.1 	39.0  	3.0 
   # Relative Humidity (%) 	9.2 	77.5 	 100 	2.0 
   # http://cohemisferico.uprm.edu/prysig/pdfs/pres_amendez09.pdf
   TOPT=273.15+24.1
   TMIN=273.15-8.3
   TMAX=273.15+39
   RHOPT=77.5
   RHMIN=9.2
   RHMAX=100
   ## vegetation
   TOPT=273.15+22.0
   TMIN=273.15-0.6
   TMAX=273.15+43.6
   RHOPT=93.1
   RHMIN=0
   RHMAX=100

   amin=T-TMIN
   amax=T-TMAX
   amean=np.power((T-TOPT),2)
   
   a1=amin*amax
   a2=a1-amean
   HTT=np.power((a1/a2),3)

   #if(T>TMAX): HTT=0
   #if(T<TMIN): HTT=0

   b1=RH-RHOPT
   b2=RH-RHMIN
   b3=np.power((b1/b2),2)
   HRH=np.power((1+b3), -1)

   #if(RH>RHMAX): HRH=0
   #if(RH<RHMIN): HRH=0
   
   fhab=HTT
   #fhab=HRH 
   #fhab=HRH*HTT
   #fhab=HTT

   fhab=np.where(fhab<0,0, fhab)   
   return(fhab)

def plot_rh_habitability_npp(oce2, npp2):

   npp2=np.nan_to_num(npp2,0)
   relative_npp=npp2/np.max(npp2)
   plt.title("Ocean percent vs NPP of earth-like planet", size=18)
   plt.xlabel("Ocean percent", size=16)
   plt.ylabel("Relative NPP", size=16)
   plt.xticks(fontsize=16)
   plt.yticks(fontsize=16)   
   plt.plot(oce2, relative_npp, lw=3, color="green", label="relative npp")
   plt.grid()
   plt.legend(fontsize=16)
   #plt.plot(fox2)
   plt.show()

def calc_spatial_mean(xr_da, lon_name="lon", lat_name="lat"):

   coslat = np.cos(np.deg2rad(xr_da[lat_name]))
   return (xr_da * coslat).sum(dim=[lon_name, lat_name]) / (
       coslat.sum(lat_name) * len(xr_da[lon_name])
   )

def habitability_f_T_RH(T, RH):

   # https://phl.upr.edu/projects/standard-primary-habitability-sph
   #https://phl.upr.edu/projects/standard-primary-habitability-sph
   # Surface Temperature (°C) 	-8.3 	24.1 	39.0  	3.0 
   # Relative Humidity (%) 	9.2 	77.5 	 100 	2.0 
   # http://cohemisferico.uprm.edu/prysig/pdfs/pres_amendez09.pdf
   TOPT=273.15+24.1
   TMIN=273.15-8.3
   TMAX=273.15+39
   RHOPT=77.5
   RHMIN=9.2
   RHMAX=100
   ## vegetation
   TOPT=273.15+22.0
   TMIN=273.15-0.6
   TMAX=273.15+43.6
   RHOPT=93.1
   RHMIN=0
   RHMAX=100

   amin=T-TMIN
   amax=T-TMAX
   amean=np.power((T-TOPT),2)
   
   a1=amin*amax
   a2=a1-amean
   HTT=np.power((a1/a2),3)

   #if(T>TMAX): HTT=0
   #if(T<TMIN): HTT=0

   b1=RH-RHOPT
   b2=RH-RHMIN
   b3=np.power((b1/b2),2)
   HRH=np.power((1+b3), -1)

   #if(RH>RHMAX): HRH=0
   #if(RH<RHMIN): HRH=0
   
   #fhab=HTT
   #fhab=HRH 
   fhab=HRH*HTT
   #fhab=HTT

   fhab=np.where(fhab<0,0, fhab)   
   return(fhab)
   

def fit_log_poly(inx, iny, outx, order1):

   iny=np.nan_to_num(iny,0)
   logx=np.log10(inx)
   logy=np.log10(iny)
   logx2=np.log10(outx)
   #m,b = np.polyfit(logx, logy, 1)
   coef = np.polyfit(logx,logy,order1)
   poly1d_fn = np.poly1d(coef)     
   logy2=poly1d_fn(logx2)
   logy2=np.where(logy2<1e-3,1e-3,logy2)
   y2=np.power(10, logy2)
   return(y2)

def plot_rain(oce2, rain, landrain, landrain2):

   rain=np.nan_to_num(rain,0)
   plt.plot(oce2*100, rain, lw=12, color="#afafff")
   plt.plot(oce2*100, rain, lw=6, color="#7f7fff")
   plt.plot(oce2*100, rain, lw=3, color="blue", label="annual precipitation mm")
   plt.plot(oce2*100, landrain, lw=12, color="#afffaf")
   plt.plot(oce2*100, landrain, lw=6, color="#7fff7f")
   plt.plot(oce2*100, landrain, lw=3, color="green", label="annual precip. in land mm")
   plt.plot(oce2*100, landrain2, lw=6, color="red", label="precip per land area mm")
   plt.title("Ocean percent vs mean precipitation of Earth-like planet", size=18)
   plt.xlabel("Ocean per cent", size=16)
   plt.ylabel("Annual precipitation mm", size=16)
   plt.xticks(fontsize=16)
   plt.yticks(fontsize=16)
   plt.grid()
   plt.ylim([0,2000])
   plt.legend(fontsize=16)
   #plt.plot(fox2)
   plt.show()

1. 1. main program ....

ds01 = xr.open_dataset(iname1) ds02 = xr.open_dataset(iname2) ds03 = xr.open_dataset(iname3) ds04 = xr.open_dataset(iname4) ds05 = xr.open_dataset(iname5)

dss1=[ds01, ds02, ds03, ds04, ds05]

print(dss1)

lsms1=[] temps1=[] npps1=[] prs1=[] prlands1=[] forests1=[] hurs1=[]

1. 1. clt alb glac prw ps psl sic sit snd ts tas veggpp veglai vagplantc vegsoilc

for n in range(0,5):

   des1=dss1[n]
   lsm0=des1.lsm
   temp0=des1.tas
   npp0=des1.vegnpp
   pr0=des1.pr
   vegf0=des1.vegf

   lsm1=np.array(lsm0[0])
   lsm2=np.mean(np.array(calc_spatial_mean(lsm0[0:11])))
   temp2=np.mean(np.array(calc_spatial_mean(temp0[0:11])))
   mean_hur=np.array(calc_spatial_mean(des1.hur))
   hur2=np.mean(mean_hur[:,9])
   npp2=np.mean(np.array(calc_spatial_mean(npp0[0:11])))*3600*24*365
   pr2=np.mean(np.array(calc_spatial_mean(pr0[0:11])))*3600*24*365*1000
   prland0=lsm0*pr0
   prland2=np.mean(np.array(calc_spatial_mean(prland0[0:11])))*3600*24*365*1000
   forestland0=lsm0*vegf0
   forestland2=np.mean(np.array(calc_spatial_mean(forestland0[0:11])))
   #print((1-lsm2), npp2)
   lsms1.append(lsm2)
   temps1.append(temp2)
   npps1.append(npp2)
   prs1.append(pr2)
   hurs1.append(hur2)
   prlands1.append(prland2)
   forests1.append(forestland2)

lsms1=np.array(lsms1) temps1=np.array(temps1) npps1=np.array(npps1) prs1=np.array(prs1) prlands1=np.array(prlands1) forests1=np.array(forests1)/lsms1 hurs1=np.array(hurs1)

oceans1=(1-lsms1)

print ("oceanfrac:",oceans1) print("tsurfair:", temps1) print ("precip:",prs1) print("precipinland:", prlands1) print("rh:", hurs1) print("npp:",npps1) print("forestinland:", forests1)

print("oceanfrac tsair precip precipinland rh npp forestinland") for n in range(0,5):

   print(round(oceans1[n],3),round(temps1[n],3),  round(prs1[n],0), round(prlands1[n]),round(hurs1[n],9),round(npps1[n],9), forests1[n])

1. oceanfrac: [0.70743782 0.70743782 0.70743782 0.70743782 0.70743782]
2. tsurfair: [309.39436378 299.50792493 291.2773634 282.42979126 273.48767201]
3. precip: [1791.17692459 1406.23741651 1148.68153485 920.87646469 719.58338112]
4. precipinland: [433.22410895 324.14415487 270.4668116 216.25549866 161.74408844]
5. rh: [81.2404649 79.13623723 78.28275684 75.45653865 71.86652914]
6. npp: [0.08187775 0.07653408 0.07373737 0.06272353 0.04615242]
7. forestinland: [0.18520579 0.15474801 0.15202314 0.13788886 0.08973894]
8. oceanfrac tsair precip precipinland rh npp forestinland
9. 0.707 309.394 1791.0 433 81.240464904 0.08187775 0.18520579438451157
10. 0.707 299.508 1406.0 324 79.136237233 0.076534077 0.1547480120569912
11. 0.707 291.277 1149.0 270 78.282756842 0.073737367 0.1520231381282169
12. 0.707 282.43 921.0 216 75.456538652 0.062723535 0.1378888614032402
13. 0.707 273.488 720.0 162 71.866529138 0.046152417 0.08973893801103539

1. quit(-1)

oceans2=np.linspace(0,1,100) lands2=1-oceans2

aus1=np.array([0.94, 0.97,1,1.03, 1.06]) aus2=np.linspace(0.9,1.1,10000)

1. print(oceans2)

1. npps2=fit_log_poly(oceans1, npps1, oceans2,2)
2. prs2=fit_log_poly(oceans1, prs1, oceans2,2)
3. prlands2=fit_log_poly(oceans1, prlands1, oceans2,2)

npps2=fit_log_poly(aus1, npps1, aus2,3)

temps2=fit_log_poly(aus1, temps1, aus2,3)

hurs2=fit_log_poly(aus1, hurs1, aus2,3)

hab1=habitability_f_T_RH(temps2, hurs2)

hab1=np.round(((100*hab1)/np.max(hab1)), 3)

ESI_t = np.power(abs(temps2 / 288), 0.7)

PHI_t = np.power(abs(288 / temps2), 0.2) PHI_a = np.exp(-0.2 * np.power((aus2 - 1), 2)) PHI=PHI_t*PHI_a

hab2=ESI_t

1. prlands3=prlands2/lands2

1. plt.plot(aus2,hab1)

1. plt.plot(aus2,hab2)
2. plt.plot(aus2,hurs2)

1. plt.plot(aus2,temps2-273.15)

1. plt.plot(aus2,npps2)

1. plt.show()

1. print(hab1)

plot_habitability_au(aus2, hab1)

1. plot_rain(oceans2,prs2, prlands2, prlands3)

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=5780 # assumption Earth-like

1. 1. Gliese 12 b

planetname="TEST1"

1. S1=1.0

au=0.94 S1=1/(au*au)

flux=S1*Sol year=360 rotationperiod=1 radius=1 ## rel to earth mass=math.pow(radius, 3.7)

eccentricity=0.013 obliquity=23.5

fixedorbit=True synchronous=False

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

1. pressure=(patm1+patm2)/2
2. pressure=10
3. pressure=0.1

pressure=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. map generator #1
3. 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-1W.bmp" ## diku seed 11111 sea level -0.05 wrinky

1. top1=7000
2. floor1=-8000

top1=8400 floor1=-11200

land_height_scaling_exponent_1=5 ## land height scaling exponent 15 16 tex

sealevel1=0.065 ## 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(".")