Conversion utilities for input and output - test

[oculus99]

Python3 script for creating own "earth.txt" input file from png image

-- coding: utf-8 --

v 0001

Python code

to convert ascii map for pycli

test only

input file is large paletted greyscale planet map dem

basic usage python3 aski.py --in map.png --out earth.txt

python3 aski.py --in rapa.png --out earth.txt --cols 360

import sys, random, argparse
import numpy as np
import math
from PIL import Image

gray scale level values from:

http://paulbourke.net/dataformats/asciiart/

70 levels of gray

#gscale1 = "$@b%8&WM#*oahkbdpqwmZO0QLCJUYXzcvunxrjft/|()1{}[]?-_+~<>i!lI;:,"^`'. "

10 levels of gray

#gscale2 = '@%#*+=-:. '
gscale1= '01'
gscale2= '01'

def getAverageL(image):
im = np.array(image)
w,h = im.shape
return np.average(im.reshape(w*h))

def makemap(fileName, cols, scale, moreLevels):
global gscale1, gscale2
image = Image.open(fileName).convert('L')
W, H = image.size[0], image.size[1]
print("input image dims: %d x %d" % (W, H))
w = W/cols
h = w/scale
rows = int(H/h)
print("cols: %d, rows: %d" % (cols, rows))
print("tile dims: %d x %d" % (w, h))
if cols > W or rows > H:
print("Image too small for specified cols!")
exit(0)
aimg = []
for j in range(rows):
y1 = int(j*h)
y2 = int((j+1)*h)
if j == rows-1:
y2 = H
aimg.append("")

    for i in range(cols):
        x1 = int(i*w)
        x2 = int((i+1)*w)
        if i == cols-1:
            x2 = W
        img = image.crop((x1, y1, x2, y2))
        avg = int(getAverageL(img))
        if moreLevels:
            gsval = gscale1[int((avg*69)/255)]
        else:
            #gsval = gscale2[int((avg*9)/255)]
            gsval = gscale2[int((avg*1)/255)]
        aimg[j] += gsval
        aimg[j] += ' '
return aimg

def main():
descStr = "PNG into ASCII map."
parser = argparse.ArgumentParser(description=descStr)
parser.add_argument('--in', dest='imgFile', required=True)
parser.add_argument('--scale', dest='scale', required=False)
parser.add_argument('--out', dest='outFile', required=False)
parser.add_argument('--cols', dest='cols', required=False)
parser.add_argument('--morelevels',dest='moreLevels',action='store_true')
args = parser.parse_args()

imgFile = args.imgFile
outFile = 'out.txt'
if args.outFile:
   outFile = args.outFile
scale = 1
if args.scale:
   scale=1
   scale = float(args.scale)
cols = 180
if args.cols:
   cols = int(args.cols)
print('Generating ASCII map...')
aimg = makemap(imgFile, cols, scale, args.moreLevels)
f = open(outFile, 'w')
for row in aimg:
   f.write(row + '\n')
f.close()
print("ASCII map written to %s" % outFile)

if name == 'main':
main()

Create sample elevation netcdf and mask

create land mask png

end planet map dem

with python3 "noise" perlin noise

v0001

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

from PIL import Image

import netCDF4 as nc

from PIL import Image, ImageOps

import noise

def save_netcdf_elev(fn, dataa1):
shp1=np.shape(dataa1)
dimx1=shp1[1]
dimy1=shp1[0]
ds = nc.Dataset(fn, 'w', format='NETCDF4')
lat = ds.createDimension('lat', dimy1)
lon = ds.createDimension('lon', dimx1)
lats = ds.createVariable('lat', 'f4', ('lat',))
lons = ds.createVariable('lon', 'f4', ('lon',))
data = ds.createVariable('data', 'f4', ('lat', 'lon',))
data.units = 'm'
lats[:] = np.arange(-90, 90, 180/dimy1)
lons[:] = np.arange(-180.0, 180, 360/dimx1)
data[:, :] = np.flipud(dataa1)
ds.close()

def profile_twohalf(ines1):
ones1=np.where(ines1 < 0, np.power(ines1,2), np.power(ines1,0.5) )
return(ones1)

def ringed(ines1):
ones1=2 * (0.5 - np.abs(0.5 - ines1) )
return(ones1)

def savepng_gray_16(outname, I):
maxiima1=(256.0*256.0)-1.0
I16 = (((I - I.min()) / (I.max() - I.min())) * maxiima1).astype(np.uint16)
img = Image.fromarray(I16)
#img2 = img.filter(ImageFilter.GaussianBlur(3))
img.save(outname)
return(0)

def normal1(x, c,d):
q1=-math.pow((x-c),2.0)
q2=math.pow(2*d,2.0)
fu1= 1+ math.exp(q1/q2 )
fu2=1-1/fu1
return ( fu2 )

def logistic2(x, a, b, c):
return ( 1 / ( c + math.exp((-x*a)+b) ) )

def gaus(X,C,X_mean,sigma):
return (Cnp.exp(-(X-X_mean)**2/(2sigma**2)))

def siirros(X,C,X_mean,sigma):
return( 1 - gaus(X,C,X_mean,sigma) )

def normalize(arra):
mini=np.min(arra)
maxi=np.max(arra)
deltai=maxi-mini
return((arra-mini)/deltai)

def normalizeminus1plus1(arra):
mini=np.min(arra)
maxi=np.max(arra)
deltai=maxi-mini
basis1=(arra-mini)/deltai
basis2=basis1*2-1
return(basis2)

def create_perlin_noise(dimx1, dimy1, seed1, scale, ar, octaves,persistence,lacunarity):

shape = (dimy1,dimx1)
pi2=math.pi*2
pi=math.pi
fiis=np.linspace(0,pi,dimy1)
phiis=np.linspace(0,pi2,dimx1)
leni=len(fiis)
lenj=len(phiis)
noise0 = np.zeros(shape)

for i in range(shape[0]):
	for j in range(shape[1]):
		theta=fiis[i]
		phi=phiis[j]
		x=ar * math.sin(theta) * math.cos(phi)
		y=ar * math.sin(theta) * math.sin(phi)
		z=ar * math.cos(theta)
		noise0[i][j] = noise.pnoise3(x,y,z,octaves=octaves,persistence=persistence,lacunarity=lacunarity, repeatx=1024, repeaty=1024,repeatz=1024,base=seed1)
		
noise0=normalizeminus1plus1(noise0)   	
return(noise0)

def create_noise(dimx1, dimy1, seed1, skala1):
scale = dimx1/2
ar=0.66
octaves = 13
persistence = 0.5
lacunarity = 2.0
shape=(dimy1, dimx1)
noise0000=create_perlin_noise(dimx1, dimy1, seed1, scale, ar, octaves,persistence,lacunarity)
return(noise0000)

############################

betw 0 ...1

continentlevel1=0.8

maskname1="mask.png"
demname1="dembase.png"
demname2="demland.png"
elevname1="demland.nc"
elevname2="demlandsea.nc"

seed1=2
seed2=1
skala1=1/4
dimx1=7201
dimy1=3601

sealevel1=-7000
depthm1=-1500
heightm1=3000

print ("Creating noises for map .. ")

perlin noise

noise0=normalize(create_noise(dimx1, dimy1, seed1, skala1))

ringed noise

noise1=normalize(ringed(create_noise(dimx1, dimy1, seed2, skala1)))

two exponents a - attempt to make more realistic profile

##qnoise1=normalize(np.exp(noise1))
qnoise1=noise1-continentlevel1

#plt.imshow(qnoise1)
#plt.show()

#quit(-1)

qnoise2=np.where(qnoise1 <0, qnoise1, np.exp(np.exp(qnoise1)))

qnoise3=normalize(qnoise2)

map1=np.copy(qnoise2)
map1=map1

dem5=np.copy(map1)
dem5=np.where(dem5 <0, (1-dem5)depthm1, dem5heightm1)
#dem5=np.where(dem5 >0, dem5/10, dem5)

#plt.imshow(dem5)
#plt.show()

#quit(-1)

dem5=dem5+sealevel1

#plt.imshow(dem5)
#plt.show()

dem1=np.copy(dem5)
dem2=np.copy(dem5)
dem3=np.copy(dem5)
dem4=np.copy(dem5)

dem1[dem1<0]=0
dem4[dem4<0]=0

dem1png=np.uint8(dem1*250)

dem3[dem3<0]=0
dem3[dem3>0]=1

landpixels1=np.count_nonzero(dem3)
allpixels1=dimx1dimy1
landpros1=int((landpixels1100)/allpixels1)
minelev1=int(np.min(dem5))
maxelev1=int(np.max(dem5))
landmean1=int(np.mean(dem4))

print("Land % ", landpros1)
print("Z min Max ", minelev1, maxelev1)
print("Landmean ", landmean1)

dem3png=np.uint8(dem3*255)

#plt.imshow(dem3)

#plt.imshow(dem1)
#plt.imshow(dem3png)
plt.imshow(dem3png)
##plt.imshow(dem4)
#plt.imshow(dem5, cmap="jet")
#plt.imshow(dem5, cmap="viridis")

#plt.imshow(dem5,cmap='gist_earth')

#plt.contour(dem5)

shaded=0

if(shaded==1):
x, y = np.gradient(dem4)
slope = np.pi/2. - np.arctan(np.sqrt(xx + yy))
aspect = np.arctan2(-x, y)
altitude = np.pi / 4.
azimuth = np.pi / 2.
shaded = np.sin(altitude) * np.sin(slope)+ np.cos(altitude) * np.cos(slope)* np.cos((azimuth - np.pi/2.) - aspect)
plt.imshow(shaded, cmap='Greys', origin='lower', alpha=0.5)

plt.show()

#quit(-1)

print ("Save noise mask.")

im = Image.fromarray(dem3png, "L")
im.save(maskname1)
#print ("Save demland")
im = Image.fromarray(dem1png, "L")
im.save(demname2)

save_netcdf_elev(elevname1, dem4)
save_netcdf_elev(elevname2, dem5)

Create temperature neycdf from putput file temp_f

pycli output to netcdf file

test only v 0001c

import matplotlib.pyplot as plt
from netCDF4 import Dataset
import numpy as np
import sys

#input_file0 = np.loadtxt('temp_f.txt', delimiter=' ')

indata0=np.genfromtxt("temp_f.txt", delimiter=" ", usemask=False)

#print(input_file0)

#input_file=np.array(input_file0).astype(float)

print(indata0)

#print (np.shape(input_file))

sheip1=np.shape(indata0)
dimx=sheip1[1]
dimy=sheip1[0]

print(dimx, dimy)

indata1=np.array(indata0).astype(float)

indata=np.rot90(indata1, axes=(0,1))

note flipud!

#indata=np.flipud(indata)

#plt.imshow(indata)

#plt.show()

#quit(-1)

fn = 'temp.nc'
ds = Dataset(fn, 'w', format='NETCDF4')

lat = ds.createDimension('lat', dimy)
lon = ds.createDimension('lon', dimx)

lats = ds.createVariable('lat', 'f4', ('lat',))
lons = ds.createVariable('lon', 'f4', ('lon',))
Temperature = ds.createVariable('Temperature', 'f4', ('lat', 'lon',))
Temperature.units = 'deg C'

#lats[:] = np.arange(-90, 90, 180/dimy)
lats[:] = np.arange(90, -90, -180/dimy)

#quit(-1)

lons[:] = np.arange(-180.0, 180, 360/dimx)

Temperature[:, :] = indata1-273.15

ds.close()

Simple -6.5 c/100m downscaling with dem netcdf

v0003

downscale flat land temperature

with dem

in temperatures are from level 0

6.5 C/1000 m

import numpy as np
##import scipy
import skimage
import matplotlib.pyplot as plt
import netCDF4 as nc

def save_netcdf_temp(fn, dataa1):
shp1=np.shape(dataa1)
dimx1=shp1[1]
dimy1=shp1[0]
ds = nc.Dataset(fn, 'w', format='NETCDF4')
lat = ds.createDimension('lat', dimy1)
lon = ds.createDimension('lon', dimx1)
lats = ds.createVariable('lat', 'f4', ('lat',))
lons = ds.createVariable('lon', 'f4', ('lon',))
data = ds.createVariable('Temperature', 'f4', ('lat', 'lon',))
data.units = 'deg C'
lats[:] = np.arange(-90, 90, 180/dimy1)
lons[:] = np.arange(-180.0, 180, 360/dimx1)
data[:, :] = np.flipud(dataa1)
ds.close()

fn1 = 'temp.nc'
fn2 = 'demland.nc'

ds1 = nc.Dataset(fn1)
ds2 = nc.Dataset(fn2)

land1=ds2['data'][:]
temp1=ds1['Temperature'][:]

land2=np.flipud(land1)
temp2=temp1
shapeland2=np.shape(land2)
dimx1=shapeland2[1]
dimy1=shapeland2[0]
lats1= np.arange(-90, 90, 180/dimy1)
lons1= np.arange(-180.0, 180, 360/dimx1)

#plt.imshow(land1)
#plt.contour(land1)

#plt.show()

deltemp2=land2*-6.5/1000
##t2=np.resize(temp2, shapeland2)
t2=new_image = skimage.transform.resize(temp1, shapeland2, order=3)

t3=t2+deltemp2
##t3=np.flipud(t3)

#plt.imshow(temp2)

#plt.imshow(land2)

#plt.contour(temp2, alpha=0.75, cmap="jet")
#plt.contour(t2, alpha=0.75, cmap="jet")

#plt.contour(t3, alpha=0.75, cmap="OrRd")

plt.imshow(t3, cmap="jet", extent=[-180,180,-90,90])
#plt.imshow(t3)

plt.contour(t3, alpha=0.75, cmap="jet", extent=[-180,180,90,-90])

plt.show()

t4=t3

save_netcdf_temp("temp_ds.nc", t4)