...
 
Commits (2)
......@@ -5,3 +5,6 @@
*.ipynb
*.vts
Untitled*
*.scad
*.gif
*.dat
__author__ = "Cedrick COPOL"
__version__ = "1.0"
__version__ = "1.1"
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
......@@ -38,7 +33,6 @@ 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)):
......@@ -74,7 +68,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"]
_formats = ["png", "vts", "vtk"]
if not fmt in _formats:
raise ValueError("fmt ({}) must be in {}".format(fmt, _formats))
......@@ -89,30 +83,34 @@ 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))
z = np.array(data[6:], dtype="float")
elevation = z - z.min()
if fmt == "vts":
#To export data stored in a rectilinear grid (AKA irregular Cartesian grid):
from evtk.hl import gridToVTK
X = [Ox + i*cellsize for i in range(ncols)]
Y = [Oy + j*cellsize for j in range(nrows)]
X = [Ox + i*cellsize for i in range(ncols)]
Y = [Oy + j*cellsize for j in range(nrows)]
# X = np.linspace(Ox, Ox + (ncols-1)*cellsize, ncols, dtype="float64")
# Y = np.linspace(Oy, Oy + (nrows-1)*cellsize, nrows, dtype="float64")
x = np.zeros((ncols,nrows,1))
y = np.zeros((ncols,nrows,1))
z = data.flatten().reshape(ncols, nrows, 1)
if fmt in ["vts", "vtk"]:
#To export data stored in a rectilinear grid (AKA irregular Cartesian grid):
try:
from evtk.hl import gridToVTK
except:
raise ImportError("evtk not found. See https://bitbucket.org/pauloh/pyevtk/src/")
x = np.zeros((ncols, nrows, 1))
y = np.zeros((ncols, nrows, 1))
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
z = z.flatten(order="F").reshape(ncols, nrows, 1)
print(gridToVTK(outputname, x, y, z, pointData = {"elevation" : z.flatten(order="F")}))
print(gridToVTK(outputname, y, x, z,
pointData = {"height" : z.flatten(order="F"),
"elevation": elevation.flatten()
}
)
)
elif fmt == "png":
xsize = (ncols-1) * cellsize
......@@ -212,6 +210,6 @@ 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)
#asc2fmt("lidar_souf_5m_ll.asc", "Soufriere", xmax=2500, fmt="vts", cellsize=5)
# 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", fmt="vts", cellsize=5)