File:Late ordovician tas 445ma co2 150 s1318 obl245 ecc001 mvelp45.png

1. 3. Exoplasim planet running code
   4. exoplasim example
   5. stepper code
2. 2. 03.06.2023 0000.0004c
3. 1. in ubuntu you must install
4. 1. pip3 install exoplasim[netCDF4]
   2. not
   3. "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,-35,-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,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)

1. 1. 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 ---")

1. !!! Sun

1. input_dem="dema1.nc" ##dem of exoplanet
2. 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"

1. 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

1. 1. a_baseflux=1365.2
2. a_baseflux=1354
3. a_baseflux=1318

a_baseflux=1312

a_yearlength=360 a_radius=1.0 a_gravity=9.81

1. 1. 1. !!! ROTATION 7 !!!!

a_rotationperiod=0.883 ## ordovician

1. 1. !!!!! earth nowadays
2. a_eccentricity1=0.01671022
3. a_obliquity1=23.4392811
4. a_lonvernaleq1=102.94719

1. a_eccentricity1=0.06
2. a_obliquity1=22
3. a_lonvernaleq1=270

a_eccentricity1=0.01 a_obliquity1=24.5 a_lonvernaleq1=45

a_pCO21=150e-6

print("Exoplasim ...")

1. 1. 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(".")

1. "R" temperature downscaler

1. 1. exoplasim output && paleodem temperature 6.5 C/1km downscaler

1. 1. 4.6.2023 v 0000.0001a

library(stringr)

library(raster) library(ncdf4) library(rgdal) library(viridis)

library(pixmap) library(png) library(bmp) library(OpenImageR) library(graphx)

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")

}

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<-t(mask2) idem3<-t(idem2) mask4 <- apply(mask3, 2, rev) idem4 <- apply(idem3, 2, rev) r <- writePNG(mask4, maskname1) r <- writePNG(idem4, 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
 

1. save_as_nc(demname2,dem2, "Band1", "Band1", "")
2. save_as_nc(maskname1,mask1, "Band1", "Band1", "")
3. save_as_nc(maskname2,mask2, "Band1", "Band1", "")
4. save_as_nc(maskname3,mask3, "Band1", "Band1", "")
5. save_as_nc(maskname4,mask4, "Band1", "Band1", "")

}

downscale_temperature_65 <- function (coarse_rastera, fine_rastera) { ## methods: 0 delta, 1 spatialeco, 2 dissever, 3 temperature lapse 6.5 C/1 km lm

   print ("Downscaling temperature, 6.5 C/1 km")			

coarse_raster<-coarse_rastera fine_raster<-fine_rastera p1<-fine_raster p2<-fine_raster

plot(fine_raster) plot(coarse_raster, col=viridis(200))

coarseoro<- resample(p1, coarse_raster) coarseoro_big<-resample(coarseoro, p1) orodelta<-(p1-coarseoro_big)

baset1 <- resample(coarse_raster, p1)

raster_stack <- stack(p1,p2)

orotemp<-orodelta*0.0065*-1

coarseorotemp<- resample(orotemp, coarse_raster) coarseorotemp_big<-resample(coarseorotemp, p1)

orotempdelta<-orotemp-coarseorotemp_big

outtemp<-baset1+orotempdelta

#outtempr<-rotate(outtemp)

#plot(outtempr)

     return(outtemp)

}

1. 1. warning: check, if you need
   2. flip input MOST raster, You can do this witk "lsm" raster

1. 1. only input files

dataname1<-"./exo_445_omega45_co2_150/MOST.00030.nc" finename1<-"./maps/Map77_PALEOMAP_6min_Late_Ordovician_445Ma.nc"

1. 1. output files ...

coarsename1<-"tassi.nc"

maskname1= "mask_land.png" maskname2= "mask_sea.png"

edgename1="edge.png"

1. maskname3= "mask_land_rotated.nc"
2. maskname4= "mask_sea_rotated.nc"

demname1= "planet_dem.nc"

uxsize=720 uysize=360

ext1<-extent(-180,180,-90,90) crs1<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

1. stop(-1)

nc_data <- nc_open(dataname1)

lon <- ncvar_get(nc_data, "lon") lat <- ncvar_get(nc_data, "lat", verbose = F) t <- ncvar_get(nc_data, "time")

temp.array <- ncvar_get(nc_data, "tas") # store the data in a 3-dimensional array mask.array <- ncvar_get(nc_data, "lsm") ## for ratster flip check only

1. plot(mask.array)

1. stop(-1)

temp.array <-temp.array-273.15

dim(temp.array) nc_close(nc_data)

print(t)

1. temp.slice <- temp.array[, , 6]

print (dim(temp.array))

temp.slice<-apply(temp.array, c(1,2), mean) mask.slice<-apply(mask.array, c(1,2), mean)

1. plot(temp.slice)

1. stop(-1)

r <- raster(t(temp.slice), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))

1. plot(r)

maskr <- raster(t(mask.slice), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))

plot(maskr)

1. stop(-1)

1. 1. check whethet must fli ras

1. 1. r <- flip(r, direction='y')

extent(r)<-ext1 crs(r)<-crs1 extent(maskr)<-ext1 crs(maskr)<-crs1

1. plot(r)

1. plot(maskr)

writeRaster(r, coarsename1, "CDF", overwrite=TRUE)

1. stop(-1)

coarser1<-raster(coarsename1) finer1<-raster(finename1)

plot(coarser1)

1. plot(finer1)

1. stop(-1)

finer2<-finer1

finer2[finer2<0.0]<-0.0

plot(coarser1) plot(finer2)

1. stop(-1)

matrix1 <- matrix(1:(uxsize*uysize), nrow=uysize, ncol=uxsize)

sabluna1 <- raster(matrix1) extent(sabluna1) <- ext1 crs(sabluna1)<-crs1

finer4<-resample(finer2, sabluna1) finer0<-finer4

rout1<-downscale_temperature_65(coarser1, finer4)

plot(rout1)

1. stop(-1)

crs(rout1)<-crs1 extent(rout1 )<-ext1

1. image(rout1)

writeRaster(rout1, "temp_dskaled1.nc", overwrite=TRUE, format="CDF", varname="atm_temp", varunit="degC",

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

finer3<-finer0

finer3[finer3>0]<-255 finer3[finer3<=0]<-0

writeRaster(finer3, "mask1.tiff", overwrite=TRUE, format="GTiff", varname="mask", varunit="unit",

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

dem1<-raster(finename1)

writeRaster(dem1, "demo1.nc", overwrite=TRUE, format="CDF", varname="mask", varunit="unit",

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

masksealevel1=0

create_mask_rasters(finename1,maskname1,demname1)

image1 = readImage(maskname1)

1. image1

1. print (dim(image1))

image2=image1

1. str(image2)

1. stop(-1)

res1 = edge_detection(image2, method = 'Frei_chen', conv_mode = 'same')

1. plot(res1)

res2=res1

res2[res2[]==0]=2 res2=res2-1

writeImage(res2, edgename1)

image2[image2[]==0]=2 image2=image2-1

writeImage(image2, maskname2)

1. createmask(dem1,demname1,maskname1, maskname2, maskname3, maskname4, masksealevel1)