Commit bd1dd181 authored by Cedrick Copol's avatar Cedrick Copol

Sauvegarde rapide. Implémentation a revoir

parent 13ca0ceb
......@@ -5,6 +5,41 @@ __version__ = "1.0"
import sys
import numpy as np
import matplotlib.pyplot as plt
# from PIL import Image
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
print(elevation.max())
nb_slices = int(elevation.max() // slice_height)
def make_image(data, outputname, size=(2400, 1220)):
fig = plt.figure()
......@@ -16,8 +51,11 @@ def make_image(data, outputname, size=(2400, 1220)):
ax.set_axis_off()
fig.add_axes(ax)
ax.imshow(data, cmap="gray", origin="lower")
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 asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
......@@ -51,7 +89,9 @@ def asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
if cellsize is None:
cellsize = float(data[4][1])
# Topography
#data = np.array(data[6:], dtype="float64", order="F")
data = np.array(data[6:], dtype="float64")
print(data.shape, (nrows, ncols))
if fmt == "vts":
#To export data stored in a rectilinear grid (AKA irregular Cartesian grid):
......@@ -68,6 +108,7 @@ def asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
for i, Xi in enumerate(X):
x[i,:,0] = Xi
#x = np.expand_dims(Xi, axis=0)
for j, Yj in enumerate(Y):
y[:,j,0] = Yj
......@@ -106,12 +147,71 @@ def asc2fmt(filename, outputname, fmt="png", cellsize=None, widthIRL=400):
elif fmt == "stl":
try:
import stl_tools
import stl
except:
raise ImportError("stl_tools not found")
#stl_tools.numpy2stl
raise ImportError("stl not found. Please install numpy-stl")
mesh = stl.mesh.Mesh(data)
mesh.save(outputname+".stl")
elif fmt == "scad":
# Surface
with open(f"{outputname}.scad", "w") as txt:
print("//surface.scad", file=txt)
print(f"//width: {ncols}", file=txt)
print(f"//depth: {nrows}", file=txt)
print(f"//height: {elevation.max()}", file=txt)
print(f"widthIRL = {widthIRL};", file=txt)
print(f"ech = {widthIRL}/{ncols};", file=txt)
print("scale(ech*[1,1,1])", file=txt)
print(f'surface(file = "{outputname}.dat", center = true, convexity = 5);', file=txt)
# Attention pensez à adapter l'echelle en divisant l'élevation par la longueur de la cellule
np.savetxt(f"{outputname}.dat", elevation.transpose()/cellsize, delimiter=" ")
def slices(height=3):
for sl in range(nb_slices):
data_tmp = np.ones(elevation.shape)
data_tmp[ elevation<=sl*slice_height ] = 0
#data_tmp[ (elevation<sl*slice_height) | (elevation>(sl+1)*slice_height) ] = 0
make_image(data_tmp, f"Soufriere_{sl:03d}.png", size=data_tmp.shape[::-1])
def laser(height=3):
# https://matplotlib.org/examples/pylab_examples/contour_image.html
extent = (min(X), max(X), min(Y), max(Y))
thickness = 3 # mm
widthIRL = 400 # mm
scale = xsize/widthIRL
slice_height = thickness * scale
nb_slices = int(elevation.max() // slice_height)
levels = np.arange(nb_slices)*slice_height
norm = cm.colors.Normalize(vmax=abs(elevation).max(), vmin=-abs(elevation).max())
cmap = cm.PRGn
size=elevation.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(nb_slices):
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(elevation, levels, colors='k', origin='upper', extent=extent)
plt.plot(X[430], Y[707], 'ro')
plt.plot(X[419], Y[466], 'ro')
plt.savefig(f"contour_{level:03d}.png", dpi=dpi)
if __name__ == "__main__":
asc2fmt("lidar_souf_5m_ll.asc", "Soufriere", fmt="png", cellsize=5, widthIRL=400)
#asc2fmt("lidar_souf_5m_ll.asc", "Soufriere_limit", fmt="png", cellsize=5, widthIRL=400)
for w in [400, 100]:
asc2fmt("lidar_souf_5m_ll.asc", f"Soufriere_{w}", fmt="png", cellsize=5, widthIRL=w)
#asc2fmt("lidar_souf_5m_ll.asc", "Soufriere", xmax=2500, fmt="vts", cellsize=5)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment