python - Algorithm for finding approximate area of the ground via elevation grid? -

i have 200x200 grid/array (i.e. 40,000 values) of 10mx10m "pixels" values indicate average elevation of land @ point.

across grid there vast connected regions elevation values 0m because represent actual sea.

question: there quick algorithm obtaining approximate area of land? understand may multiply 200^2*10^2 obtain rough approximation of area, of values vary quite bit.

i think know rather expensive way, i.e. summing of triangles vertices lie @ elevation. there faster/easier way?

numpy , scipy tools kind of problem. here's synthetic 200×200 landscape points on 10 metre grid , height ranging 40 metres above sea level:

>>> import numpy np >>> xaxis = yaxis = np.arange(0, 2000, 10) >>> x, y = np.meshgrid(xaxis, yaxis) >>> z = np.maximum(40 * np.sin(np.hypot(x, y) / 350), 0) 

we can have @ in matplotlib:

>>> import matplotlib.pyplot plt >>> import mpl_toolkits.mplot3d.axes3d axes3d >>> _, axes = plt.subplots(subplot_kw=dict(projection='3d')) >>> axes.plot_surface(x, y, z, cmap=plt.get_cmap('winter')) >>> plt.show() 

now, number of points on land (that is, height greater 0) trivial compute, , can multiply size of grid square (100 m² in question) estimate of land area:

>>> (z > 0).sum() * 100 1396500 

but question, understand want more accurate estimate, 1 takes account slope of land. 1 way make mesh of triangles covering land, , add area of triangles.

first, turn coordinate arrays array of points (a point cloud):

>>> points = np.vstack((x, y, z)).reshape(3, -1).t >>> points array([[  0.000000e+00,   0.000000e+00,   0.000000e+00],        [  1.000000e+01,   0.000000e+00,   1.142702e+00],        [  2.000000e+01,   0.000000e+00,   2.284471e+00],        ...,         [  1.970000e+03,   1.990000e+03,   3.957136e+01],        [  1.980000e+03,   1.990000e+03,   3.944581e+01],        [  1.990000e+03,   1.990000e+03,   3.930390e+01]]) 

second, use scipy.spatial.delaunay triangulate in 2 dimensions, getting surface mesh:

>>> scipy.spatial import delaunay >>> tri = delaunay(points[:,:2]) >>> len(tri.simplices) 79202 >>> tri.simplices array([[39698, 39899, 39898],        [39899, 39698, 39699],        [39899, 39700, 39900],        ...,         [19236, 19235, 19035],        [19437, 19236, 19237],        [19436, 19236, 19437]], dtype=int32) 

the values each triangle in triangulation indexes in points array of 3 points in triangle.

third, select triangles have land in them:

>>> land = (points[tri.simplices][...,2] > 0).any(axis=1) >>> triangles = tri.simplices[land] >>> len(triangles) 27858 

fourth, compute areas of these triangles:

>>> v0, v1, v2 = (points[triangles[:,i]] in range(3)) >>> areas = np.linalg.norm(np.cross(v1 - v0, v2 - v0), axis=1) / 2 >>> areas array([ 50.325028,  50.324343,  50.32315 , ...,  50.308673,  50.313157, 50.313649]) 

finally, add them up:

>>> areas.sum() 1397829.2847141961 

that's not different original estimate, expected because slopes shallow.