import numpy as np

# given xx, yy, zz

grad_x, grad_y = np.gradient(zz)
dxx = np.gradient(xx, axis=1)
dyy = np.gradient(yy, axis=0)
# partial derivative
f_x = grad_x / dxx
f_y = grad_y / dyy
# Calculate the magnitude of the gradient
grad_mag = np.sqrt(f_x**2 + f_y**2)

# Find the maximum gradient and its location
max_index = np.unravel_index(np.argmax(grad_mag, axis=None), grad_mag.shape)
steepest_slope = grad_mag[max_index]
steepest_slope_location = (xx[max_index], yy[max_index])

print(f"The steepest slope is {steepest_slope} at location {steepest_slope_location}")

# Compute the magnitude of the surface normal
normal_mag = np.sqrt(f_x**2 + f_y**2 + 1)
# Compute the cosine of the angle between the surface normal and the upward vector
cos_theta = 1 / normal_mag
# Compute the angle in radians
theta_rad = np.arccos(cos_theta)
# Convert the angle to degrees
theta_deg = np.degrees(theta_rad)

steepest_angle = theta_deg[max_index]
steepest_angle_location = (xx[max_index], yy[max_index])

print(f"The steepest angle from the Z axis is {steepest_angle} degrees at location {steepest_angle_location}")

**Output**