import gdal
import gdalconst
import numpy
from stl import mesh
np = numpy
!pip install pymesh
DEPRECATION: Python 2.7 will reach the end of its life on January 1st, 2020. Please upgrade your Python as Python 2.7 won't be maintained after that date. A future version of pip will drop support for Python 2.7.
Collecting pymesh
Downloading https://files.pythonhosted.org/packages/94/ca/5a49647246c189856a5fb7d36382893330788035d7734b1447e05641baa6/pymesh-1.0.2.tar.gz
Requirement already satisfied: numpy in ./.venv/lib/python2.7/site-packages (from pymesh) (1.16.3)
Building wheels for collected packages: pymesh
Building wheel for pymesh (setup.py) ... done
Stored in directory: /Users/jared/Library/Caches/pip/wheels/a6/d7/76/ed7121207ca5d4310797e72ae0921f54935d5afaf08bd67855
Successfully built pymesh
Installing collected packages: pymesh
Successfully installed pymesh-1.0.2
def convert_file(fname, outfile, sample=5, chunk_config=None):
file = gdal.Open(fname, gdalconst.GA_ReadOnly)
data = file.ReadAsArray(0, 0, file.RasterXSize, file.RasterYSize)
if chunk:
data = chunk(data, *chunk_config)
scaled = data[::sample, ::sample]
m = geo_mesh2(scaled, 0, 0, len(scaled) - 1, len(scaled[0]) - 1)
rescale(m)
print 'ok'
mesh.Mesh(m).save(outfile)
# Timp at 1/3 arc second
convert_file('./raw_data/USGS_NED_13_n41w112_ArcGrid/grdn41w112_13/', 'timp_3.stl', sample=1, chunk_config=(1-.38, 1-.66, .25, .25))
2702 2702 ok
# Overview files, sampled at 5 arc seconds
# Wasatch
convert_file('./raw_data/N039W111/N039W111_AVE_DSM.tif', 'n39w111.stl')
convert_file('./raw_data/N039W112/N039W112_AVE_DSM.tif', 'n39w112.stl')
convert_file('./raw_data/N040W112/N040W112_AVE_DSM.tif', 'n40w112.stl')
# Grand Canyon
convert_file('./raw_data/N036W113/N036W113_AVE_DSM.tif', 'n36w113.stl')
719 719 ok 719 719 ok 719 719 ok 719 719 ok
def rescale(m):
m['vectors'][:,:,0] /= m['vectors'][:,:,0].max() / 10.0
m['vectors'][:,:,1] /= m['vectors'][:,:,1].max() / 10.0
m['vectors'][:,:,2] /= 4000
def load(fname):
file = gdal.Open(fname, gdalconst.GA_ReadOnly)
return file.ReadAsArray(0, 0, file.RasterXSize, file.RasterYSize)
data = load('./raw_data/N040W112/N040W112_AVE_DSM.tif')
from matplotlib import pyplot
import matplotlib
matplotlib.image.
scaled = chunk(data, 1-0.37, 1-0.68, 0.26, 0.26)
m = geo_mesh2(scaled, 0, 0, len(scaled) - 1, len(scaled[0]) - 1)
rescale(m)
print 'ok'
mesh.Mesh(m).save('./timp.stl')
935 935
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-139-68683d94b92d> in <module>() 1 scaled = chunk(data, 1-0.37, 1-0.68, 0.26, 0.26) ----> 2 m = geo_mesh2(scaled, 0, 0, len(scaled) - 1, len(scaled[0]) - 1) 3 rescale(m) 4 print 'ok' 5 mesh.Mesh(m).save('./timp.stl') <ipython-input-128-ebce252437c5> in geo_mesh2(geo, x, y, w, h) 23 total = p0 # numpy.concatenate([p0, p1]) 24 m = numpy.zeros((w - 1) * (h - 1), dtype=mesh.Mesh.dtype) ---> 25 m['vectors'] = total 26 return m ValueError: could not broadcast input array from shape (872356,4,3) into shape (872356,3,3)
# Grand Canyon?
data = load('./raw_data/N036W113/N036W113_AVE_DSM.tif')
scaled = chunk(data, 1-0.2, 1-0.2, 0.2, 0.2)
m = geo_mesh2(scaled, 0, 0, len(scaled) - 1, len(scaled[0]) - 1)
rescale(m)
print 'ok'
mesh.Mesh(m).save('./grand_canyon.stl')
719 719
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-129-2c80c9b49419> in <module>() 2 data = load('./N036W113/N036W113_AVE_DSM.tif') 3 scaled = chunk(data, 1-0.2, 1-0.2, 0.2, 0.2) ----> 4 m = geo_mesh2(scaled, 0, 0, len(scaled) - 1, len(scaled[0]) - 1) 5 rescale(m) 6 print 'ok' <ipython-input-128-ebce252437c5> in geo_mesh2(geo, x, y, w, h) 23 total = p0 # numpy.concatenate([p0, p1]) 24 m = numpy.zeros((w - 1) * (h - 1), dtype=mesh.Mesh.dtype) ---> 25 m['vectors'] = total 26 return m ValueError: could not broadcast input array from shape (515524,4,3) into shape (515524,3,3)
scaled = data[::5, ::5]
def chunk(data, x, y, w, h):
fw = len(data)
fh = len(data[0])
x0 = int(fw * (x - w/2))
x1 = int(fw * (x + w/2))
y0 = int(fh * (y - h/2))
y1 = int(fh * (y + h/2))
return data[x0:x1,y0:y1]
# (-3.7, -6.8), w 2.4 to a side
def geo_mesh2(geo, x, y, w, h):
data = numpy.zeros(w * h * 3).reshape((w, h, 3))
print len(data)
i0 = numpy.indices((w, h))
print len(data)
indices = i0.swapaxes(0,2).swapaxes(0,1)
data[:,:,:2] = indices
min = geo.min()
data[:, :, 2] = (geo[x:x+w,y:y+h] - min)
# now we have x, y, z just fine
# return data
p0 = numpy.array([
data[:-1,:-1].reshape((-1, 3)),
data[1:, :-1].reshape((-1, 3)),
data[1:, 1:].reshape((-1, 3))
]).swapaxes(0, 1)
p1 = numpy.array([
data[:-1, :-1].reshape((-1, 3)),
data[1:, 1:].reshape((-1, 3)),
data[:-1, 1:].reshape((-1, 3))
]).swapaxes(0, 1)
total = numpy.concatenate([p0, p1])
m = numpy.zeros((w - 1) * (h - 1) * 2, dtype=mesh.Mesh.dtype)
m['vectors'] = total
return m
def geo_pymesh(geo, x, y, w, h):
data = numpy.zeros(w * h * 3).reshape((w, h, 3))
print len(data)
i0 = numpy.indices((w, h))
print len(data)
indices = i0.swapaxes(0,2).swapaxes(0,1)
data[:,:,:2] = indices
min = geo.min()
data[:, :, 2] = (geo[x:x+w,y:y+h] - min)
vertices = data.reshape((-1, 3))
faces = numpy.array([
])
faces = numpy.array([
data[:-1,:-1].reshape((-1, 3)),
data[1:, :-1].reshape((-1, 3)),
data[1:, 1:].reshape((-1, 3)),
data[:-1, 1:].reshape((-1, 3)),
]).swapaxes(0, 1)
m = numpy.zeros((w - 1) * (h - 1) * 2, dtype=mesh.Mesh.dtype)
m['vectors'] = total
return m
def show_numpy(data):
# Optionally render the rotated cube faces
from matplotlib import pyplot
from mpl_toolkits import mplot3d
# Create a new plot
figure = pyplot.figure()
axes = mplot3d.Axes3D(figure)
# Render the cube faces
m = mesh.Mesh(data.copy())
axes.add_collection3d(mplot3d.art3d.Poly3DCollection(m.vectors))
# Auto scale to the mesh size
# scale = m.points.flatten()
# axes.auto_scale_xyz(scale, scale, scale)
# axes.auto_scale_xyz(m.points[::3].flatten(), m.points[1::3].flatten(), m.points[2::3].flatten())
axes.auto_scale_xyz(m.vectors[:,:,0].flatten(), m.vectors[:,:,1].flatten(), m.vectors[:,:,2].flatten())
# Show the plot to the screen
pyplot.show()
m = geo_mesh2(data, 0, 0, len(data) - 1, len(data[0]) - 1)
print 'ok'
mesh.Mesh(m).save('./mountains.stl')
3599 3599 ok
# Create 3 faces of a cube
mdata = numpy.zeros(6, dtype=mesh.Mesh.dtype)
# Top of the cube
mdata['vectors'][0] = numpy.array([[0, 1, 1],
[1, 0, 1],
[0, 0, 1]])
mdata['vectors'][1] = numpy.array([[1, 0, 1],
[0, 1, 1],
[1, 1, 1]])
# Front face
mdata['vectors'][2] = numpy.array([[1, 0, 0],
[1, 0, 1],
[1, 1, 0]])
mdata['vectors'][3] = numpy.array([[1, 1, 1],
[1, 0, 1],
[1, 1, 0]])
# Left face
mdata['vectors'][4] = numpy.array([[0, 0, 0],
[1, 0, 0],
[1, 0, 1]])
mdata['vectors'][5] = numpy.array([[0, 0, 0],
[0, 0, 1],
[1, 0, 1]])
# Since the cube faces are from 0 to 1 we can move it to the middle by
# substracting .5
mdata['vectors'] -= .5
# Generate 4 different meshes so we can rotate them later
meshes = [mesh.Mesh(mdata.copy()) for _ in range(4)]
# Rotate 90 degrees over the Y axis
meshes[0].rotate([0.0, 0.5, 0.0], math.radians(90))
# Translate 2 points over the X axis
meshes[1].x += 2
# Rotate 90 degrees over the X axis
meshes[2].rotate([0.5, 0.0, 0.0], math.radians(90))
# Translate 2 points over the X and Y points
meshes[2].x += 2
meshes[2].y += 2
# Rotate 90 degrees over the X and Y axis
meshes[3].rotate([0.5, 0.0, 0.0], math.radians(90))
meshes[3].rotate([0.0, 0.5, 0.0], math.radians(90))
# Translate 2 points over the Y axis
meshes[3].y += 2
show_numpy(mdata)