...
 
Commits (3)
......@@ -6,34 +6,8 @@ import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
filename = "lidar_souf_5m_ll.asc"
cellsize = 5
with open(filename, "r") as txt:
data = [line.strip().split() for line in txt.readlines()]
dim = data[:6]
# Number of points
ncols, nrows = [int(d[1]) for d in data[:2]]
# Origin coordinates
Ox, Oy = [float(d[1]) for d in data[2:4]]
if cellsize is None:
cellsize = float(data[4][1])
# Topography
data = np.array(data[6:], dtype="float64")
xsize = (ncols-1)*cellsize # m
X = [Ox + i*cellsize for i in range(ncols)]
Y = [Oy + j*cellsize for j in range(nrows)]
elevation = data - data.min()
widthIRL = 400 # mm
# https://matplotlib.org/examples/pylab_examples/contour_image.html
extent = (min(X), max(X), min(Y), max(Y))
slice_height = 50
nb_slices = int(elevation.max() // slice_height)
def make_image(data, outputname, size=(2400, 1220)):
fig = plt.figure()
......@@ -47,10 +21,42 @@ def make_image(data, outputname, size=(2400, 1220)):
ax.imshow(data, cmap="gray", aspect="equal")
plt.savefig(outputname, dpi=dpi)
# print(fig.get_size_inches()*fig.dpi, fig.dpi, fig.get_dpi())
# print(size, Image.open(outputname).size)
# plt.savefig(outputname)
def laser_cut(X, Y, Z, thickness=3, prefix="laser"):
# https://matplotlib.org/examples/pylab_examples/contour_image.html
extent = (min(X), max(X), min(Y), max(Y))
widthIRL = 400 # mm
scale = (extent[1]-extent[0])/widthIRL
slice_height = thickness * scale
nb_slices = int(Z.max() // slice_height)
levels = np.arange(nb_slices)*slice_height
lim = abs(Z).max()
norm = cm.colors.Normalize(vmax=lim, vmin=-lim)
cmap = cm.PRGn
size=Z.shape[::-1]
fig = plt.figure()
dpi = fig.get_dpi()
fig.set_size_inches(size[0]/float(dpi), size[1]/float(dpi))
#ax.imshow(elevation, cmap="gray", extent=extent)
#ax.imshow(elevation, cmap=cmap, norm=norm, extent=extent)
print("nb couches", nb_slices)
for level in range(1,nb_slices+1):
plt.clf()
# Remove axis and fill the space
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
levels = [level*slice_height]
plt.contour(Z, levels, colors='k', origin='upper', extent=extent)
plt.plot(X[430], Y[707], 'ro')
plt.plot(X[419], Y[466], 'ro')
plt.savefig(f"{prefix}_{level:03d}.png", dpi=dpi)
def asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
"""Convert asc (lidar ASCII format)
......@@ -68,7 +74,7 @@ def asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
width in mm of the final stl. If you choose to output a png, this script will
print some values to set in Cura to have this width.
"""
_formats = ["png", "vts", "vtk"]
_formats = ["png", "vts", "vtk", "laser"]
if not fmt in _formats:
raise ValueError("fmt ({}) must be in {}".format(fmt, _formats))
......@@ -78,10 +84,12 @@ def asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
dim = data[:6]
# Number of points
ncols, nrows = [int(d[1]) for d in data[:2]]
# Origin coordinates (Actually these can be ignored in our case)
Ox, Oy = [float(d[1]) for d in data[2:4]]
if cellsize is None:
cellsize = float(data[4][1])
# Topography
z = np.array(data[6:], dtype="float")
elevation = z - z.min()
......@@ -113,26 +121,16 @@ def asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
)
elif fmt == "png":
xsize = (ncols-1) * cellsize
ysize = (nrows-1) * cellsize
hsize = data.max() - data.min()
print("dim:", xsize, ysize, hsize)
# Image dimension (pixel)
# pixel_height = 1200.
# pixel_scale = pixel_height/ncols
# print("scale 1 pix = {} m".format(pixel_scale))
# print(pixel_scale)
# pixel_depth = pixel_scale * ysize
# make_image(data, outputname+"_001.png", size=(pixel_width, pixel_depth))
# print(data.shape, (ncols,nrows) )
# WARNING data shape and pixel size are reversed
make_image(data, outputname+".png", size=data.shape[::-1])
# make_image(data, outputname+".png", size=(ncols*cellsize, nrows*cellsize))
make_image(elevation, outputname+".png", size=elevation.shape[::-1])
print(outputname+".png")
# dimension cura en mm
xsize = (ncols-1) * cellsize
hsize = elevation.max()
scale = widthIRL/xsize
width = scale*xsize
depth = scale*ysize
# depth = scale*ysize
height = scale*hsize
print("---------------------------")
print("Importer l'image dans Cura.")
......@@ -165,6 +163,9 @@ def asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
# Attention pensez à adapter l'echelle en divisant l'élevation par la longueur de la cellule
np.savetxt(f"{outputname}.dat", elevation.transpose()/cellsize, delimiter=" ")
elif fmt == "laser":
laser_cut(X, Y, elevation, thickness=3, prefix="laser")
def slices(height=3):
for sl in range(nb_slices):
data_tmp = np.ones(elevation.shape)
......@@ -210,6 +211,7 @@ def laser(height=3):
plt.savefig(f"contour_{level:03d}.png", dpi=dpi)
if __name__ == "__main__":
# for w in [400, 100]:
# asc2fmt("lidar_souf_5m_ll.asc", f"Soufriere_{w}", fmt="png", cellsize=5, widthIRL=w)
for w in [400]:
asc2fmt("lidar_souf_5m_ll.asc", f"Soufriere_{w}", fmt="png", cellsize=5, widthIRL=w)
asc2fmt("lidar_souf_5m_ll.asc", "Soufriere", fmt="vts", cellsize=5)
asc2fmt("lidar_souf_5m_ll.asc", "Soufriere", fmt="laser", cellsize=5)