I advise using Interpolation and NIntegrate if your data is smooth. Otherwise, if there is noise in your data, the derivatives will be noisy and you will probably want to smooth the data and use BSplineFunction.

Here is an example using Interpolation.

1. You will need to arrange your data so that it corresponds to what Interpolation operates on:

   yourData = {{0.00505051, 0., -18.164}, {0.010101, 
      0., -14.5719}, {0.0151515, 0., -9.9892}, {0.020202, 
      0., -5.792}, {0.0252525, 0., -1.7631}, {0.030303, 0., 0.9439}}
   
   dataForInterpolate = yourData /. {x_, y_, z_} :> {{x, y}, z}
   

Unfortunately, we don't have enough of your data to continue, so here is an example with "example data"

exampleData = 
 Flatten[
  Table[ {{x, y}, Sin[x] Cos[y]}, {x, 0, 1, .01}, {y, 0, 1, .01}], 1]

zInterpolation = Interpolation@exampleData

Get interpolators for derivatives:

dzdxInterpolation = Derivative[1, 0][zInterpolation]
dzdyInterpolation = Derivative[0, 1][zInterpolation]

Numerically integrate:

NIntegrate[
 Sqrt[1 + dzdxInterpolation[x, y]^2 + dzdyInterpolation[x, y]^2 ], {x,
   0, 1}, {y, 0, 1}]

Compare:

With[{dx = D[Sin[x] Cos[y], x ], dy = D[Sin[x] Cos[y], y ] },
 NIntegrate[Sqrt[1 + dx^2 + dy^2], {x, 0, 1}, {y, 0, 1}]
 ]

exampleData = Flatten[Table[{{x, y}, f {x, y}}, {x, 0, 1, .01}, {y, 0, 1, .01}], 1]

Has a syntax error f{x,y} should be f[x,y] if f is a function that takes two arguments. I will check your data sometime later.

Thank you so much for your help. Have applied your method and now everything is working just fine. I have one more question and I would be grateful if you can advise me on this as well. I have to calculate a parameter that in order to do so my data must be transformed into polar coordinates (data is lots of points in x,y,z format) and it is using Angular Spectrum. The result is shown in the polar spectrum graph (I Guess the other name of the graph is rosette graph but I am not sure). Do you have any suggestions on how the code for plotting the graph might be?

Thank you so much in advance.

Thank you so much for your help. I am developing a platform based on the data that I have to calculate some surface parameters. so I don't have a specific number for A and Ar. that is why I used dimx and dimy to limit the surface as desired. the data file was not readable in here so you can download it from the following link. I appreciate your help in advance. enter link description here

You must first reconstruct the surface $z = f(x,y)$. One way to do it is via Radial Basis Function (RBF) interpolation. One you are satisfied with your interpolation, use it in the double integral. You may use NIntegrate to evaluate the integral.