from__future__importdivision#Avoids the floor of the mathematical result of division if the args are ints or longsimportnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlib.colorsasmcimportSSFM2Dimportrandom#Parameter assignmentsmax_level=8#This is the exponent controlling the grid size. In this case N=2^8=256. Use only integers.N=2**max_level
sigma=1#Variation for random Gauss generation (standardised normal distribution)H=0.8#Hurst exponent (0.8 is a recommended value for natural phenomena)seed=random.random()#Setting the seed for random Gauss generationprint('The lattice size is '+str(N)+'x'+str(N))#Lattice initializationLattice=np.zeros((256,256))#Calling Spectral fBm functionLattice=SSFM2D.SpectralSynthesisFM2D(max_level,sigma,H,seed,normalise=True,bounds=[0,1])#Plotting the original 256x256 latticeM=np.zeros((257,257))foriinrange(0,256):forjinrange(0,256):M[i][j]=Lattice[i][j]#Normalizing the output matrixprint('Original sum: '+str(round(M[-257:,-257:].sum(),3)))M=M/M[-257:,-257:].max()#Matrix normalization with respect to max#Lattice statisticsprint('Normalized sum: '+str(round(M[-257:,-257:].sum(),3)))print('Normalized max: '+str(round(M[-257:,-257:].max(),3)))print('Normalized min: '+str(round(M[-257:,-257:].min(),3)))print('Normalized avg: '+str(round(M[-257:,-257:].mean(),3)))#Determining the footprintfootprint=0#Initializing the footprint count variableforiinrange(M.shape[0]):forjinrange(M.shape[1]):ifM[i][j]>0.15:# Change here to set the failure thresholdfootprint=footprint+1print('Event footprint: '+str(round(footprint*100/(256*256),2))+'%')#Plotting the 256x256 "mother" matrixplt.imshow(M[-257:,-257:].T,origin='lower',interpolation='nearest',cmap='Reds',norm=mc.Normalize(vmin=0,vmax=M.max()))title_string=('Spatial Loading of a Generic Natural Hazard')subtitle_string=('Inverse FFT on Spectral Synthesis')plt.suptitle(title_string,y=0.99,fontsize=17)plt.title(subtitle_string,fontsize=8)plt.show()#Making a custom list of tick mark intervals for color bar (assumes minimum is always zero)numberOfTicks=5ticksListIncrement=M.max()/(numberOfTicks)ticksList=[]foriinrange((numberOfTicks+1)):ticksList.append(ticksListIncrement*i)plt.tick_params(axis='x',labelsize=8)plt.tick_params(axis='y',labelsize=8)cb=plt.colorbar(orientation='horizontal',format='%0.2f',ticks=ticksList)cb.ax.tick_params(labelsize=8)cb.set_label('Loading',fontsize=12)plt.show()plt.xlim(0,255)plt.xlabel('Easting (Cells)',fontsize=12)plt.ylim(255,0)plt.ylabel('Northing (Cells)',fontsize=12)plt.annotate('fractional Brownian motion on a 256x256 lattice | H=0.8 | Dim(f)= '+str(3-H),xy=(0.5,0),xycoords=('axes fraction','figure fraction'),xytext=(0,0.5),textcoords='offset points',size=10,ha='center',va='bottom')plt.annotate('Max: '+str(round(M[-257:,-257:].max(),3))+' | Min: '+str(round(M[-257:,-257:].min(),3))+' | Avg: '+str(round(M[-257:,-257:].mean(),3))+' | Footprint: '+str(round(footprint*100/(256*256),2))+'%',xy=(0.5,0),xycoords=('axes fraction','figure fraction'),xytext=(0,15),textcoords='offset points',size=10,ha='center',va='bottom')#Producing the 101x101 image(s)numfig=int(raw_input('Insert the number of 101x101 windows to produce: '))N=np.zeros((101,101))x=0#Counts the iterations towards the chosen number of imagesforxinrange(1,numfig+1):print('Image no. '+str(x)+' of '+str(numfig))north=int(raw_input('Northing coordinate (0 thru 155, integer): '))east=int(raw_input('Easting coordinate (0 thru 155, integer): '))foriinrange(101):forjinrange(101):N[i][j]=M[north+i][east+j]#Writing X, Y and values to a .csv file from scratchimportnumpyimportcsvwithopen('C:\\Users\\Francesco\\Desktop\\Python_files\\csv\\fBm_101x101_'+str(x)+'of'+str(numfig)+'.csv','w')asf:#Change directory if necessarywriter=csv.writer(f)writer.writerow(['X','Y','Value'])for(x,y),valinnumpy.ndenumerate(M):writer.writerow([x,y,val])#Plotting the 101x101 "offspring" matricesplt.imshow(N[-101:,-101:].T,origin='lower',interpolation='nearest',cmap='Reds',norm=mc.Normalize(vmin=0,vmax=M.max()))title_string=('Spatial Loading of a Generic Natural Hazard')subtitle_string=('Inverse FFT on Spectral Synthesis | Origin in the 256x256 matrix: '+str(east)+' East; '+str(north)+' North')plt.suptitle(title_string,y=0.99,fontsize=17)plt.title(subtitle_string,fontsize=8)plt.show()#Making a custom list of tick mark intervals for color bar (assumes minimum is always zero)numberOfTicks=5ticksListIncrement=M.max()/(numberOfTicks)ticksList=[]foriinrange((numberOfTicks+1)):ticksList.append(ticksListIncrement*i)plt.tick_params(axis='x',labelsize=8)plt.tick_params(axis='y',labelsize=8)cb=plt.colorbar(orientation='horizontal',format='%0.2f',ticks=ticksList)cb.ax.tick_params(labelsize=8)cb.set_label('Loading',fontsize=12)plt.show()plt.xlim(0,100)plt.xlabel('Easting (Cells)',fontsize=12)plt.ylim(100,0)plt.ylabel('Northing (Cells)',fontsize=12)plt.annotate('fractional Brownian motion on a 101x101 lattice | H=0.8 | Dim(f)= '+str(3-H),xy=(0.5,0),xycoords=('axes fraction','figure fraction'),xytext=(0,0.5),textcoords='offset points',size=10,ha='center',va='bottom')plt.annotate('Max: '+str(round(N[-101:,-101:].max(),3))+' | Min: '+str(round(N[-101:,-101:].min(),3))+' | Avg: '+str(round(N[-101:,-101:].mean(),3))+' | Footprint: '+str(round(footprint*100/(101*101),2))+'%',xy=(0.5,0),xycoords=('axes fraction','figure fraction'),xytext=(0,15),textcoords='offset points',size=10,ha='center',va='bottom')