26 KiB
+++ title = "A Thoroughly Digital Artifact" slug = "a-thoroughly-digital-artifact" date = "2023-01-11" [taxonomies] tags = ["3dprinting", "CAD", "GIS", "CNC", "art", "sundries", "proclamation"] +++
A birthday wish
Last summer, I wanted to get my wife something nice for her birthday. For many years, she had expressed an occasional and casual desire for a topographic carving of the state of California, where we live, and I thought it might be something I could figure out how to get her. In the end, after many dozens of hours of work, five weeks, and several hundred dollars paid to a professional CNC machine shop, I had the artifact shown in the picture above. This is the story of its creation, starting from knowing almost nothing about GIS, cartography, or CNC machining.
First steps
Before you ask, I did not do a ton of research before embarking on this. As I write this, about six months later, it only now occurred to me to do a basic search for an actual physical thing I could buy, and luckily it seems that CNC-carved wooden relief maps of the whole state are not trivially easy to come by, so, phew!
No, my first step was to see if there were any shops in the area that could carve something out of nice plywood, about a week before the intended recipient's birthday. I found one that was less than ten minutes away, and filled out their web contact form. They had a field for material, and I said, "some nice plywood between 0.75 and 1.0 inches thick or similar" (I didn't know exactly what was available and wanted to give broad acceptable parameters), and under "project description", I wrote,
A relief map of California, carved from wood. Height exaggerated enough to visibly discern the Santa Monica mountains. I can provide an STL file if needed.
For some [incorrect] reason that I only later examined1, I just sort of assumed that the shop would have a library of shapes available for instantiating into whatever material medium you might need. But just in case, I included that hedge about being able to provide an STL file. Needless to say, that was a bluff.
Also needless to say, my bluff was immediately called, and I had the following exchange with the shop:
CNC Shop: STL can work but I can’t manipulate it, which could save some money. If possible can it be exported to an .igs or .iges or .stp format?
Me: Yeah, STP should be no problem. Can you give a rough estimate of the cost for 1x2-foot relief carving?
Shop: Without seeing the drawings, I can’t give even a close price but in the past they range from a few hundred dollars to several thousand dollars.
Me: That's totally fair! I'll get you some files in a few days.
"STP should be no problem ... I'll get you some files in a few days," was an even harder lean into the bluff; my next communication with the shop was nearly four weeks later. But that's getting ahead of things.
Meshes and solid bodies
First off, let's talk about file formats and how to represent shapes with a computer.2 I first said I could provide an STL file. STL is a pretty bare-bones format that describes the outside surface of a shape as a mesh of many, many triangles, each of which is described by three 3D points, where each point (but not necessarily each edge) of the triangle lies on the surface of the shape of the thing you're modeling. This format is popular with 3D printers, which is how I became familiar with it.
STL is simple to implement and easy for a computer to read, but if you have a model in that format that you need to manipulate, like you want to merge it with another shape, you won't have a good time. In order to actually do things like that, it needs to be converted into a CAD program's native representation of a "solid body", which is pretty much what it sounds like: a shape made of a finite volume of "stuff", and NOT just an infinitesimally thin shell enclosing an empty volume, which is what a mesh is.
In order for the CAD program to convert a mesh into a solid body, the mesh must be manifold, meaning, no missing faces (triangles), and with a clearly-defined interior and exterior (all triangles are facing in one direction relative to their interior). When there are no missing faces, it's called "water tight". You can still have "holes" in a mesh, like if you have a model of a donut3, but the surface of the donut can't have any missing faces. A valid STL file's meshes are manifold.
The CNC shop had requested a model in a format called
STP. .stp
is the extension for a "STEP" file; STEP is supposed to be short for "standard for the exchange of
product data", so someone was playing pretty fast and loose with their initialisms, but I
digress. The main thing about STEP files is that CAD programs can really easily convert them
into their native internal solid body representation, which allows easy manipulation. Another thing
about them is that a CAD program can usually turn a manifold mesh into an STP file, unless the mesh is
too complicated and your computer doesn't have enough RAM (note: foreshadowing4).
But so far, I had nothing at all. Time to get some data and see if I can turn it into a model.
Public data
My first impulse was to search USGS's website for heightmap data, but I wound up not finding anything appropriate. Searching now with the wisdom of experience and hindsight, I found this, which would have been perfect:
https://apps.nationalmap.gov/downloader/
Did I just accidentally miss it then? Did I not know how to recognize it because I didn't know what I was doing at all? The world may never know, but at least now you can benefit from my many, many missteps.
From space?
Anyway, having not found anything I could really use from the USGS, I found this site, from OpenTopography, an organization run by the UCSD Supercomputer Center, under a grant from the National Science Foundation. So, hooray for public data!
That particular page is for a particular dataset; in this case, "SRTM GL1 Global 30m". "SRTM" stands for "Shuttle Radar Topography Mission", which was a Space Shuttle mission in February, 2000, where it did a fancy radar scan of most of the land on Earth. Though, it's hard to verify that the data was not synthesized with other datasets of more recent, non-space origin, especially in places like California. But probably space was involved in some way.
In Australia, it's pronounced "g'dal"
Anyway, I'd found an open source of public data. This dataset's horizontal resolution is 1 arc second (which is why it's "GL1"), or roughly 30x30 meters, and the height data is accurate to within 16 meters. Not too shabby!
They provided the data in the form of GeoTIFFs, which are basically an image where each pixel represents one data point (so, a 30x30 square meter plot) centered at a particular location on the Earth's surface. It's a monochrome image, where absolute height is mapped to absolute brightness of each pixel, and each pixel represents an exact location in the world.
The only problem was that you could only download data covering up to 450,000 square kilometers at a
time, so I had had to download a bunch of separate files and then mosaic them together. Luckily,
there's a whole suite of open source tools called
GDAL. Among that suite is a tool called
gdal_merge.py (yes, the .py is part of the name of the tool that gets installed to your system
when you install the GDAL tools), which does exactly what I wanted:
gdal_merge.py -o ca_topo.tif norcal_topo.tif centcal_topo.tif socal_topo.tif so_cent_cal_topo.tif norcal_topo_redux.tif last_bit.tif east_ca.tif
This produced a file called ca_topo.tif. It was very large, in every sense:
Using another tool called gdalinfo, we can examine the
metadata of the mosaic we just created:
$ gdalinfo -mm ca_topo.tif
Driver: GTiff/GeoTIFF
Files: ca_topo.tif
Size is 40757, 35418
Coordinate System is:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-125.109583333326071,42.114305555553187)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-125.1095833, 42.1143056) (125d 6'34.50"W, 42d 6'51.50"N)
Lower Left (-125.1095833, 32.2759722) (125d 6'34.50"W, 32d16'33.50"N)
Upper Right (-113.7881944, 42.1143056) (113d47'17.50"W, 42d 6'51.50"N)
Lower Right (-113.7881944, 32.2759722) (113d47'17.50"W, 32d16'33.50"N)
Center (-119.4488889, 37.1951389) (119d26'56.00"W, 37d11'42.50"N)
Band 1 Block=40757x1 Type=Int16, ColorInterp=Gray
Computed Min/Max=-130.000,4412.000
If I may draw your attention to a couple things there, the image is 40,757 pixels wide and 35,418 pixels tall. The "pixel size" is 0.000277777777778 by 0.000277777777778; the units, given by the "angleunit", is degrees; 1 arc second is 1/3600th of a degree, which is 0.01754... They're degrees of arc along the surface of the Earth5, at a distance measured from the center of the planet. As previously mentioned, that translates into a size of roughly 30 meters. So if you were ever curious about how many 100-ish-foot squares you'd need to fill a rectangle that fully enclosed the entire border of California, then one billion, four-hundred-forty-three million, five-hundred-thirty-one thousand, and four-hundred-twenty-six (40,757 times 35,418) is pretty close.
The other units in there are under the "Coordinate System is" section, and are meters relative to the World Geodetic System 1984 vertical datum; the very last line is the lowest and highest points in file, which are -130 meters and 4,412 meters respectively, relative to the baseline height defined by the WGS84 ellipsoid. If you were to view the file as though it were an image, it would look like this:
This is because the highest possible value an image like that could have for a pixel is 65,5356, and the highest point in our dataset is only 4,412, which is not that much in comparison. Plus, it includes portions of not-California in the height data, and ideally, we want those places to be 0; we have a little more processing to do before we can use this.
Cartography is complicated
The first order of business is to mask out everything that's not California, and the first thing I needed for that was a shapefile that described the California state border. Luckily, that exact thing is publicly available from the state's website; thank you, State of California!
There was only one issue: the shapefile was in a different map
projection than the data in our geotiff file. A "map
projection" is just the term for how you display a curved, 3D shape (like the border of a state on the
curved surface of the Earth) on a flat, 2D surface, like a map. If you look at the line in the
output of gdalinfo above that says, ID["EPSG",4326], that is telling us the particular
projection used. EPSG 4326 uses
latitude and longitude, expressed in degrees, covers the entire Earth including the poles, and
references the WGS84 ellipsoid as the ground truth.
The shapefile was in a projection called EPSG 3857, or "Web Mercator". This is similar to EPSG 4326, except instead of using the WGS84 ellipsoid, it pretends the Earth is a perfect sphere. It only covers +/- 85-ish degrees of latitude (so not the poles), and it uses meters instead of degrees of lat/long. It's popular with online map services (like Google Maps and Open Street Maps) for displaying maps, hence the name, "Web Mercator", so you'd probably recognize the shapes of things in it.
Once again, there's a handy GDAL tool, gdalwarp, which
is for reprojecting geotiffs. So all we have to do is take our 4326-projected geotiff, use
gdalwarp to project it to 3857/Web Mercator, and then we can use the shapefile to mask off all
other height data outside the border of California.
It's almost too easy.
gdalwarp -t_srs EPSG:3857 ca_topo.tif ca_topo_mercator.tif
This gives us a 3857-projected file called ca_topo_mercator.tif. It still has over a billion
pixels in it (it's a little bigger overall, but the aspect is
much wider, with the different projection); scaling it down will be a very last step, since at that
point, it will no longer be a digital elevation map, it will just be an image. We'll get there,
just not yet.
Cracking open gdalinfo, we get:
$ gdalinfo ca_topo_mercator.tif
Driver: GTiff/GeoTIFF
Files: ca_topo_mercator.tif
Size is 36434, 39852
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-13927135.110024485737085,5178117.270359318703413)
Pixel Size = (34.591411839078859,-34.591411839078859)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-13927135.110, 5178117.270) (125d 6'34.50"W, 42d 6'51.50"N)
Lower Left (-13927135.110, 3799580.326) (125d 6'34.50"W, 32d16'33.21"N)
Upper Right (-12666831.611, 5178117.270) (113d47'17.10"W, 42d 6'51.50"N)
Lower Right (-12666831.611, 3799580.326) (113d47'17.10"W, 32d16'33.21"N)
Center (-13296983.361, 4488848.798) (119d26'55.80"W, 37d21'21.69"N)
Band 1 Block=36434x1 Type=Int16, ColorInterp=Gray
You can see that the PROJCRS[ID] value is "EPSG,3857", as expected. The "pixel size" is
"34.591411...." and the "lengthunit" is "metre". But the number of pixels is different, and the
shape is different, yet the coordinates of the bounding corners are the same as the original file's
(the latitude and longitude given as the second tuple). This is all from the Web Mercator's different
projection causing the aspect ratio to stretch horizontally, but it still represents the same area
of the planet.
The one custom script
So, the next step was use the shapefile to mask out the California border in the geotiff. Here is where GDAL failed me, and looking around now as I write this, I still can't find a specific GDAL tool for doing this. Given how useful I found all the other tools, I can't really complain, so I won't! It wasn't that hard to write something that would do it with other open source tools; I didn't even bother checking this into a git repo or anything:
#!/usr/bin/env python3
import fiona # for reading the shapefile
import rasterio # for working with the geotiff
import rasterio.mask as rmask
import sys
def main():
tif = sys.argv[1]
msk = sys.argv[2]
out = sys.argv[3]
print("input: {tif}\nmask: {msk}\noutput: {out}".format(tif=tif, msk=msk, out=out))
if input("Enter 'y' to continue: ").lower() != 'y': # double-check I don't stomp something I wanted to keep
print("See ya.")
return
with fiona.open(msk, "r") as shapefile:
shapes = [feature["geometry"] for feature in shapefile]
with rasterio.open(tif) as in_tif:
out_image, out_xform = rmask.mask(in_tif, shapes, filled=True, crop=True)
out_meta = in_tif.meta
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_xform})
for k, v in out_meta.items():
print("{}: {}".format(k, v)) # just outta curiosity
with rasterio.open(out, "w", **out_meta) as dest:
dest.write(out_image)
print("Wrote masked tif to {}".format(out))
return
if __name__ == "__main__":
main()
I include that just in case anyone else ever needs to do this, and doesn't find one of the hundreds of other examples out there already. This one is nice because you don't need to pre-process the shapefile into GeoJSON or anything, the Fiona package handles things like that transparently for you, but don't think this is great Python or anything; it's the dumbest, quickest thing I could crap out to do the task I needed to be done7.
After running that script, I had a Web Mercator-projected geotiff that included data only for places inside the state border of California. It was still enormous; the mask didn't change the shape and you can't have non-rectangular images anyway, but at this point, I had the final desired dataset. It was time to turn it into a heightmap that we could use to make a mesh.
A usable heightmap
I've been trying to be careful about referring to the image file as a "dataset" or "geotiff", vs. a "heightmap". A geotiff file is not a regular image file, it includes particular metadata and data that is meant to be interpreted as a real map of the land; each pixel in it says something about an exact, actual location in the real world.
A "heightmap" is an image file, like a geotiff, where each pixel's monochromatic intensity is meant to represent height above some lowest plane. The difference is that the height values are normalized so that the lowest value is 0, and the highest is the maximum possible value in the number range. For our geotiff, which uses 16-bit numbers as previously mentioned, that maximum possible value is 65,535. It also has no exact correspondence with anything else; it's not necessarily an accurate dataset, and won't include the GIS stuff like what projection it is, what the coordinate bounding boxes are, etc. But it is useful for turning into a mesh.
And here I get to the final GDAL tool I used,
gdal_translate. This is something that can read in a geotiff, and write out a different image
format. When in doubt, PNG is fine, I
always say. It's a simple format that nearly everything can read, and is compressed so it should be
much smaller file.
gdal_translate -of PNG -ot UInt16 -scale -130 4412 0 65535 masked_ca_topo.tif heightmap.png
Like we saw earlier, the lowest point we had was -130, and the highest
was 4,412. The -scale -130 4412 0 65535 arguments are saying, "anything with a height of -130
should be set to 0, and anything with a height of 4,412 should be set to 65,535, and anything
in-between should be set proportionately." This is a linear mapping, and preserves the relationships
between vertical features (so, if something is twice as tall as another thing, that will
still be true after being scaled), so it's "accurate" in a sense (note: more foreshadowing).
Once I had the PNG file, I used the ImageMagick convert
command to scale the file down to a reasonable size. Finally, I had something I could use to make a
mesh:
Pretty cool, right? I thought so! The detail is pretty great; that bright spot near the top is Mt. Shasta, for example; Mt. Whitney is slightly taller, but not by much, and is part of a range so it doesn't stand out the way Shasta does. It was time to start making some 3D geometry with the heightmap!
A mesh is born
Test prints
Give it a good smear
Final cut
Thank yous, lessons learned, and open questions
thank you: wife, steve at the shop, friends indulging my oversharing, NASA, Blender, FreeCAD
lesson: I started basically knowing nothing, but now retracing my steps I feel like I understand everything much better than when I was first racing through the material trying to get to a point where I could just produce the artifact.
lesson: pipeline of geotiff -> mask -> scaled heightmap -> mesh -> solid body
lesson: see whole article about GIS stuff
question: could I do it better in blender? freecad? could I have used gdal_polygonize (doubtful, given how much I had to blur and decimate)?
-
The conclusion upon examination was, "I just wasn't thinking". ↩︎
-
I'm pretty sure this is more "represent shapes with math" than with a computer, but the computer is helping us do the math and it's more relatable. ↩︎
-
I think you could also have a 2D sheet with a hole cut out of it represented by a mesh that is manifold, as long as the connectivity was correct in terms of how many shared edges and vertices there were (though this would not be a valid STL file). Imagine a cloth sheet with a hole cut out in the middle, and the edge of the hole hemmed or otherwise "sealed", which is then a manifold boundary. See this powerpoint deck for a pretty math-y overview of "mesh basics" (but not really that basic, that's just academics trolling us, don't let it bother you). If I'm wrong about a 2D sheet with a hole being possibly manifold, I invite correction! ↩︎
-
A textbook example of Chekhov's Scarce Computational Resource. ↩︎
-
Technically, it's an arc along the WGS84 ellipsoid, which is a perfectly smooth smushed sphere, which more closely matches the real shape of the Earth vs. a perfectly round sphere. ↩︎
-
Each pixel is 16 bits, so the possible values are from 0 to 2^16 - 1. 2^16 is 65536, so there you go. ↩︎
-
A friend posited at one point that my circuitous journey to the end product was the point, but I assured him that every step I took was trying to get to the end product as quickly and straightforwardly as possible. Still, I did in fact wind up learning a whole shitload of stuff, which is nice, I GUESS. ↩︎