利用ICESat角点数据通过最小二乘精化DEM

# -*- coding:utf-8 -*-import numpy as npfrom scipy.optimize import leastsqimport xlrddata = xlrd.open_workbook('icesat.xls')table = data.sheets()[0]nrows = table.nrowsncols = table.ncolsusedata = []for i in range(nrows):    #删除0值点    if table.cell(i, 4).value != 0:        usedata.append(table.row_values(i))#删除名字usedata.pop(0)#每n个点采样（n=3）usedata = usedata[ ::3 ]usedata = np.array(usedata)usedata_T = np.transpose(usedata)h_icesat = usedata_T[3]def_topo = usedata_T[4]h_2000 = 1000g_X = usedata_T[5]g_Y = usedata_T[6]funY = h_icesat-h_2000-(h_2000*def_topo)def func(p,x,y):    a,b,c,d,e,f = p    return a+b*x+c*y+d*x*y+e*x**2+f*y**2def error(p,x,y,Z,s):    print( s )    return func(p, x, y)-ZZ0 = np.zeros(6)s = "Test the number of iteration"Para = leastsq(error, Z0, args=(g_X, g_Y, funY, s))a, b, c, d, e, f = Para[0]print(a, b, c, d, e, f)err = []for j in range(len(usedata)):    err.append((func(Para[0], g_X[j], g_Y[j])-funY[j]))print (err)f = open('result.txt', 'w+')for m in Para[0]:    f.write(str(m))    f.write('\n')f1 = open('error.txt', 'w+')for n in err:    f1.write(str(n))    f1.write('\n')