Although this code is far from complete, it can still serve as an example of starting an exploratory scientific project. It is well-known in some parts of the geophysics community that power laws are ubiquitous in the relations. For example, treating the altitude in a landscape as a function, the spectrum reveals a power law (the correlation function as well) (See Chapter 7 of Fractals and Chaos in Geophysics and Geology by Turcotte). One thought to generate an artificial landscape is to start with white noise and then filter its spectrum with a power law filter, thereby smoothing it (though there are some limitations which might be covered in a future blog post, which I covered during my M.Sc. ca. 10 years ago …). You can read up more on this here and here.

Here is sample code to show that power scaling of altitude in a simulated landscape is obeyed:

# pieces of code adapted from:
# https://stackoverflow.com/questions/50314243/fourier-transform-in-python-2d
# import libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft2, ifft2
from mpl_toolkits.mplot3d import Axes3D
import random

N_x, N_y = 128, 128
range_x, range_y = np.arange(N_x), np.arange(N_y)
# real space grid vectors
xv, yv = (range_x - 0.5 * N_x), (range_y - 0.5 * N_y)
x, y = np.meshgrid(range_x, range_y, sparse=False, indexing='ij')
f = np.random.normal(0,1, [N_x, N_y])
F= fft2(f)
"""PLOTTING"""
fig = plt.figure()
ax = Axes3D(fig)
surf = ax.plot_surface(x,y, np.abs(f), cmap='viridis')
# for other plots I changed to
fig2 = plt.figure()
ax2 =Axes3D(fig2)
surf = ax2.plot_surface(x, y, np.abs(F), cmap='viridis')
plt.show()

Example result of white noise and corresponding spectrum

Example corresponding spectrum of instance of white noise

k_x, k_y = np.roll(np.arange(-N_x/2+1, N_x/2+1), int(N_x/2+1)), np.roll(np.arange(-N_y/2+1, N_y/2+1), int(N_x/2+1))

# reformat this; it did not copy and paste as expected
def scale_Matrix(matrix, N1, N2, beta):
new_matrix = np.array([[matrix[i,j].real if k_x[i] == 0 and k_y[j] == 0 else
matrix[i,j].real*(math.pow(1./math.sqrt((k_x[i]**2+k_y[j]**2)),beta))
for j in range(0,int(N_y))] for i in range(0,int(N_x))], np.float32)
new_matrix = new_matrix + complex(0,1) * np.array([[matrix[i,j].imag if k_x[i] == 0 and k_y[j] == 0 else
matrix[i,j].imag*(math.pow(1./math.sqrt((k_x[i]**2+k_y[j]**2)),beta))
for j in range(0,int(N_y))] for i in range(0,int(N_x))], np.float32)
return new_matrix

H = 0.25
F_scaled = scale_Matrix(F,N_x, N_y, 2 * H + 1)
f_trans = ifft2(F_scaled)

"""PLOTTING"""
fig = plt.figure()
ax = Axes3D(fig)
surf = ax.plot_surface(x,y, np.abs(f_trans), cmap='viridis')
# for other plots I changed to
fig2 = plt.figure()
ax2 =Axes3D(fig2)
surf = ax2.plot_surface(x,y, np.abs(F_scaled), cmap='viridis')
fig3 = plt.figure()
ax3 =Axes3D(fig3)
surf = ax3.plot_surface(np.log(x + 1e-2),np.log(y + 1e-2), np.log(np.abs(F_scaled)), cmap='viridis')
plt.show()

Simulated topography (coloured Gaussian noise (beta = 1.5))

Spectrum of coloured white noise (little visible)

Coloured noise spectrum on log-log-log scale

(Note that the plots are much smoother than the white-noise ones.)

#1d sections fourier transform
F_scaled_y = [np.sum(np.abs(F_scaled[i,1:int(N_y)])**2) for i in range(1,N_x//2)]
F_scaled_x = [np.sum(np.abs(F_scaled[1:int(N_x),i])**2) for i in range(1,N_y//2)]

# match exponents
plt.loglog(range_x[1:int(N_x/2)],F_scaled_x)
plt.loglog(range_x[1:int(N_x/2)],np.reciprocal([0.00004*math.pow(1.0*range_x[i],2.0) for i in range(1,int(N_x/2))]))
plt.show()
#and match exponents ... do a line of best fit
plt.loglog(range_y[1:int(N_y/2)],F_scaled_y)
plt.loglog(range_y[1:int(N_y/2)],np.reciprocal([0.00004*math.pow(1.0*range_y[i],2.0) for i in range(1,int(N_y/2))]))
plt.show()

Averaged spectrum along x (note the k^-1.5 power law)

Averaged spectrum along y (note the k^-1.5 power law)

The limitations of what has been outlined here (ignoring the limitations of the physical model) are as follows:

- There is no scale indicated in the code. This was done for simplicity. https://stackoverflow.com/questions/50314243/fourier-transform-in-python-2d has some details about how scale can be implemented.
- The code needs to be adapted for easy automation, so that it can be used (the scaleMatrix function is a step on the way there, but there needs to be more)
- The coding practices leave room for improvement. This code was written in a Jupyter notebook. It can be adapted to *.py files/modules. Making them work well-together, which includes concepts like “encapsulation.” For example, k_x is defined as a global variable and only used in the scaleMatrix function.

- An explanation of the physical meaning is lacking. I was planning on using the concepts outlined in this post to make clear what the filtering is doing, which could follow subsequently. As well, there is a chapter in a Delphi book whose name I’ve forgotten on generating fractals in a much more “physical” way (iteratively raise and lower elements in the topography) which I can compare with this (hopefully an upcoming post) to illustrate what’s happening on a physical level better.