In [ ]:
import numpy as np
import pandas as pd
import cmocean as cmo
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
from scipy.interpolate import RegularGridInterpolator
from mpl_toolkits.axes_grid1 import make_axes_locatable
import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter(action = "ignore", category = FutureWarning)
%config InlineBackend.figure_format = 'svg'
%matplotlib inline
Here we extract several parameters relative to morphometrics analysis. The analysis relies on the combined file (i.e. surface coordinates and LEC values) produced in the notebook 1-computeLEC.
The following suite of geomorphic attributes could be extracted:
In [ ]:
def assignBCs(z,nx,ny):
Zbc = np.zeros((nx + 2, ny + 2))
Zbc[1:-1,1:-1] = z
Zbc[0, 1:-1] = z[0, :]
Zbc[-1, 1:-1] = z[-1, :]
Zbc[1:-1, 0] = z[:, 0]
Zbc[1:-1, -1] = z[:,-1]
Zbc[0, 0] = z[0, 0]
Zbc[0, -1] = z[0, -1]
Zbc[-1, 0] = z[-1, 0]
Zbc[-1, -1] = z[-1, 0]
return Zbc
def cmptParams(x,y,Z):
Zbc = assignBCs(Z,x.shape[0],x.shape[1])
z1 = Zbc[2:, :-2]
z2 = Zbc[2:,1:-1]
z3 = Zbc[2:,2:]
z4 = Zbc[1:-1, :-2]
z5 = Zbc[1:-1,1:-1]
z6 = Zbc[1:-1, 2:]
z7 = Zbc[:-2, :-2]
z8 = Zbc[:-2, 1:-1]
z9 = Zbc[:-2, 2:]
dx = x[0,1]-x[0,0]
zz = z2+z5
r = ((z1+z3+z4+z6+z7+z9)-2.*(z2+z5+z8))/(3. * dx**2)
t = ((z1+z2+z3+z7+z8+z9)-2.*(z4+z5+z6))/(3. * dx**2)
s = (z3+z7-z1-z9)/(4. * dx**2)
p = (z3+z6+z9-z1-z4-z7)/(6.*dx)
q = (z1+z2+z3-z7-z8-z9)/(6.*dx)
u = (5.*z1+2.*(z2+z4+z6+z8)-z1-z3-z7-z9)/9.
with np.errstate(invalid='ignore',divide='ignore'):
grad = np.arctan(np.sqrt(p**2+q**2))
aspect = np.arctan(q/p)
hcurv = -(r*q**2-2.*p*q*s+t*p**2) / \
((p**2+q**2)*np.sqrt(1+p**2+q**2))
vcurv = -(r*p**2+2.*p*q*s+t*q**2) / \
((p**2+q**2)*np.sqrt(1+p**2+q**2))
return grad,aspect,hcurv,vcurv
In [ ]:
def regridDataSet(filename,outputfile):
azimuth=315.0
altitude=45.0
df = pd.read_csv(filename)
x = df['X']
y = df['Y']
z = df['Z']
lec = df['LEC']
dx = (x[1]-x[0])
nx = int((x.max() - x.min())/dx+1)
ny = int((y.max() - y.min())/dx+1)
xi = np.linspace(x.min(), x.max(), nx)
yi = np.linspace(y.min(), y.max(), ny)
xi, yi = np.meshgrid(xi, yi)
xyi = np.dstack([xi.flatten(), yi.flatten()])[0]
XY = np.column_stack((x,y))
zi = np.reshape(z,(ny,nx))
leci = np.reshape(lec,(ny,nx))
# Calculate gradient
Sx, Sy = np.gradient(zi)
rad2deg = 180.0 / np.pi
slope = 90. - np.arctan(np.sqrt(Sx**2 + Sy**2))*rad2deg
slp = np.sqrt(Sx**2 + Sy**2)
aspect = np.arctan2(-Sx, Sy)
deg2rad = np.pi / 180.0
shaded = np.sin(altitude*deg2rad) * np.sin(slope*deg2rad) \
+ np.cos(altitude*deg2rad) * np.cos(slope*deg2rad) \
* np.cos((azimuth - 90.0)*deg2rad - aspect)
shaded = shaded * 255
slp,aspect,hcurv,vcurv = cmptParams(xi,yi,zi)
df2 = pd.DataFrame({'X':xi.flatten(),'Y':yi.flatten(),'Z':zi.flatten(),'LEC':leci.flatten(),
'LEC':leci.flatten(),'SLP':slp.flatten(),'ASPECT':aspect.flatten(),
'HCURV':hcurv.flatten(),'VCURV':vcurv.flatten()})
df2.to_csv(outputfile,columns=['X', 'Y', 'Z', 'LEC','SLP', 'ASPECT', 'HCURV', 'VCURV'], sep=',', index=False ,header=1)
return zi,leci,slp,aspect,hcurv,vcurv
Previous python functions are used to compute the metrics defined above. To make this work for your specific case, you will need to change the values of filename and outputfile in the cell below.
In [ ]:
filename = 'dataset/combined.csv'
outputfile = 'dataset/topoData.csv'
zi,leci,slp,aspect,hcurv,vcurv = regridDataSet(filename,outputfile)
In [ ]:
def visMetrics(data=None,color=cmo.cm.delta,title='metric',figsize=(6,6),extend=None,save=None):
save = None
fig = plt.figure(figsize=figsize)
ax1 = plt.gca()
ax1.set_title(title, fontsize=10)
if extend is None:
im1 = plt.imshow(np.flipud(data), vmin=data.min(), vmax=data.max(),
cmap=color, aspect='auto')
else:
im1 = plt.imshow(np.flipud(data), vmin=extend[0], vmax=extend[1],
cmap=color, aspect='auto')
ax1.tick_params(axis='x', labelsize=8)
ax1.tick_params(axis='y', labelsize=8)
ax1.grid(False)
divider1 = make_axes_locatable(ax1)
cax1 = divider1.append_axes("right", size="5%", pad=0.05)
cbar1 = plt.colorbar(im1,cax=cax1)
cbar1.ax.tick_params(labelsize=9)
plt.tight_layout()
plt.show()
if save is not None:
fig.savefig(save,dpi=200, bbox_inches='tight')
return
In [ ]:
visMetrics(data=zi,color=cmo.cm.delta,title='Landscape elevation',figsize=(6,6),save=None)
In [ ]:
visMetrics(data=leci,color=cmo.cm.balance,title='Landscape elevational connectivity',figsize=(6,6),save=None)
In [ ]:
visMetrics(data=slp,color=cmo.cm.tempo,title='Slope',figsize=(6,6),extend=[1.5,1.6],save=None)
In [ ]:
visMetrics(data=vcurv,color=cmo.cm.balance,title='Vertical curvature',figsize=(6,6),extend=[-2,2],save=None)
In [ ]:
visMetrics(data=hcurv,color=cmo.cm.balance,title='Horizontal curvature',figsize=(6,6),extend=[-2,2],save=None)
In [ ]:
visMetrics(data=aspect,color=cmo.cm.balance,title='Aspect',figsize=(6,6),extend=[-2,2],save=None)
In [ ]: