import laspy
import numpy


source = "Grid/lidar/lidar.las"
output = "Grid/lidar/lidar_asc.asc"

cell = 1.0  # 网格尺寸（数据单位）
NODATA = 0
# 打开LIDAR的LAS文件
las = laspy.file.File(source, mode="r")
# xyz的极值
min = las.header.min
max = las.header.max

# 以m为单位获取x轴的距离
xdist = max[0] - min[0]
# 以y为单位获取y轴的距离
ydist = max[1] - min[1]

cols = int(xdist)/cell
rows = int(ydist)/cell

cols += 1
rows += 1

# 统计平均高程值数量
count = numpy.zeros((int(rows), int(cols))).astype(numpy.float32)
# 平均高程值
zsum = numpy.zeros((int(rows), int(cols))).astype(numpy.float32)
# y分辨率为负值
ycell = -1 * cell

# 将x, y的值投影到网络
projx = (las.x - min[0]) / cell
projy = (las.y - min[1]) / ycell
# 将数据转换为整数并提取部分用作索引
ix = projx.astype(numpy.int32)
iy = projy.astype(numpy.int32)
# 遍历x,y,z数组，将其添加到网格形状中并添加平均值
for x, y, z in numpy.nditer([ix, iy, las.z]):
    count[y, x] += 1
    zsum[y, x] += z

# 修改0值为1避免numpy报错以及数组出现NaN值
nonzero = numpy.where(count > 0, count, 1)
# 平均化z值
zavg = zsum / nonzero
# 将0值插入数组中避免网格出现缺口
mean = numpy.ones((int(rows), int(cols))) * numpy.mean(zavg)
left = numpy.roll(zavg, -1, 1)
lavg = numpy.where(left > 0, left, mean)
right = numpy.roll(zavg, 1, 1)
ravg = numpy.where(right > 0, right, mean)
interpolate = (lavg + ravg) / 2
fill = numpy.where(zavg > 0, zavg, interpolate)

header = "ncols        {}\n".format(fill.shape[1])
header += "nrows        {}\n".format(fill.shape[0])
header += "xllcorner        {}\n".format(min[0])
header += "yllcorner        {}\n".format(min[1])
header += "cellsize        {}\n".format(cell)
header += "NODATA        {}\n".format(NODATA)  # 0 为占位之，无意义

with open(output,"wb") as f:
    f.write(bytes(header,"UTF-8"))
    numpy.savetxt(f, fill, fmt="%1.2f")