509 lines
26 KiB
Markdown
509 lines
26 KiB
Markdown
+++
|
||
title = "A Thoroughly Digital Artifact"
|
||
slug = "a-thoroughly-digital-artifact"
|
||
date = "2023-01-17"
|
||
[taxonomies]
|
||
tags = ["3dprinting", "CAD", "GIS", "CNC", "art", "sundries", "proclamation"]
|
||
+++
|
||
|
||
![A plywood slab carved with CNC into a topographic representation of California][main_image]
|
||
|
||
# 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 examined[^introspection], 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.
|
||
|
||
![the programmer's creed: we do these things not because they are easy, but because we thought they
|
||
were going to be easy -- from twitter user @unoservix, 2016-08-05][programmers_creed]
|
||
*<center><sup><sub>me, every. single. time.</sub></sup></center>*
|
||
|
||
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.[^math-computers] I first said I could provide an *STL
|
||
file*. [STL](https://en.wikipedia.org/wiki/STL_(file_format)) 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
|
||
donut[^manifold_holes], 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
|
||
[ST**P**](https://www.fastradius.com/resources/everything-you-need-to-know-about-step-files/). `.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: foreshadowing*[^chekhovs-ram]).
|
||
|
||
![an overly-complicated mesh of a cube][meshy-cube]
|
||
*<center><sup><sub>this cube's mesh has too many vertices and edges, I hope my computer has enough
|
||
RAM to work with it</sub></sup></center>*
|
||
|
||
But at that moment, 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](https://usgs.gov)'s website for
|
||
[heightmap](https://en.wikipedia.org/wiki/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/](https://apps.nationalmap.gov/downloader/)
|
||
|
||
Did I just accidentally miss it then? Did I find it and not recognize its utility 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](https://portal.opentopography.org/raster?opentopoID=OTSRTM.082015.4326.1), from
|
||
OpenTopography, an organization run by the UCSD Supercomputer Center, under a grant from the
|
||
National Science Foundation. So, still hooray for public data!
|
||
|
||
That particular page is for a particular dataset; in this case, "[SRTM
|
||
GL1](http://www2.jpl.nasa.gov/srtm/) Global 30m". "SRTM" stands for "[Shuttle Radar Topography
|
||
Mission](https://en.wikipedia.org/wiki/Shuttle_Radar_Topography_Mission)", which was a Space Shuttle
|
||
mission in February, 2000, where it did a [fancy radar
|
||
scan](https://en.wikipedia.org/wiki/Interferometric_synthetic-aperture_radar) 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](https://gisgeography.com/srtm-shuttle-radar-topography-mission/) (which is why it's
|
||
"GL**1**"), 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 [GeoTIFF](https://en.wikipedia.org/wiki/GeoTIFF)s, 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](https://gdal.org/faq.html#what-does-gdal-stand-for). 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:
|
||
|
||
![listing of tif files with sizes][geotiff-files]
|
||
*<center><sup><sub>last_little_piece_i_swear_final_final2.tif</sub></sup></center>*
|
||
|
||
Using [another tool](https://gdal.org/programs/gdalinfo.html) called `gdalinfo`, we can examine the
|
||
metadata of the mosaic we just created:
|
||
|
||
``` text
|
||
$ 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 Earth[^wgs-ellipsoid], 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](https://en.wikipedia.org/wiki/World_Geodetic_System) vertical datum
|
||
(distances from this reference surface in the dataset are within 16 meters of the true distance in
|
||
reality); the very last line is the lowest and highest points in file, which are <a
|
||
name="minmax-height"></a>-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:
|
||
|
||
![the ca_topo image; it's hard to make out details and very dark][small_ca_topo]
|
||
*<center><sup><sub>if you squint, you can kinda see the mountains</sub></sup></center>*
|
||
|
||
This is because the highest possible value an image like that could have for a pixel is
|
||
65,535[^16-bit-ints], 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 not be represented in our dataset; 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](https://en.wikipedia.org/wiki/Shapefile) that described the
|
||
California state border. Luckily, [that exact
|
||
thing](https://data.ca.gov/dataset/ca-geographic-boundaries) 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](https://en.wikipedia.org/wiki/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](https://en.wikipedia.org/wiki/EPSG_Geodetic_Parameter_Dataset) 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](https://en.wikipedia.org/wiki/Web_Mercator_projection), 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](https://gdal.org/programs/gdalwarp.html), `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:
|
||
|
||
``` text
|
||
$ 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:
|
||
|
||
``` python
|
||
#!/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](https://geojson.org/) or anything, the
|
||
[Fiona](https://pypi.org/project/Fiona/1.4.2/) 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 done[^the-real-treasure-is-the-gd-treasure].
|
||
|
||
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 geotiff digital elevation maps, which use 16-bit numbers as previously mentioned, that
|
||
maximum possible value is 65,535. But unlike a geotiff, a generic heightmap 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](https://gdal.org/programs/gdal_translate.html) I used,
|
||
`gdal_translate`. This is something that can read in a geotiff, and write out a different image
|
||
format. When in doubt, [PNG](https://en.wikipedia.org/wiki/Portable_Network_Graphics) is fine, I
|
||
always say. It's a simple format that nearly everything can read, and is compressed so it should be
|
||
a much smaller file on disk, even if it's the same number of pixels. Smaller file size is always
|
||
easier.
|
||
|
||
> `gdal_translate -of PNG -ot UInt16 -scale -130 4412 0 65535 masked_ca_topo.tif heightmap.png`
|
||
|
||
Like we saw <a href="#minmax-height">earlier</a>, the lowest point we had in our data was -130
|
||
meters, and the highest was 4,412. The `-scale -130 4412 0 65535` arguments are saying, "anything
|
||
with a height of -130 should be totally dark, and anything with a height of 4,412 should be as
|
||
bright as possible, and anything in-between should be set proportionately." This is a linear
|
||
mapping, and preserves the relationships between vertical features (that is, if something is twice
|
||
as tall as another thing, that will still be true after being scaled), so in a sense, it's
|
||
"accurate" (*note: more foreshadowing*).
|
||
|
||
Once I had the PNG file, I used the [ImageMagick](https://imagemagick.org/script/convert.php) `convert`
|
||
command to scale the file down to a reasonable size. Finally, I had something I could use to make a
|
||
mesh:
|
||
|
||
![the heightmap made by doing a linear scale of height to brightness][scaled_heightmap]
|
||
|
||
Pretty cool, right? I thought so! The detail is pretty great; that bright spot near the top is
|
||
[Mt. Shasta](https://en.wikipedia.org/wiki/Mount_Shasta), for example;
|
||
[Mt. Whitney](https://en.wikipedia.org/wiki/Mount_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
|
||
|
||
## Back to the realm of the image
|
||
|
||
# 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)?
|
||
|
||
---
|
||
|
||
[main_image]: PXL_20220723_214758454.jpg "A plywood slab carved with CNC into a topographic representation of California"
|
||
|
||
[programmers_creed]: /images/programmers_creed.jpg "jfk overlaid with the programmer's creed: we do these things not because they are easy, but because we thought they were going to be easy"
|
||
|
||
[meshy-cube]: meshy-cube.png "an overly-complicated mesh of a cube"
|
||
|
||
[geotiff-files]: geotiff-files.png "the input geotiff files and the resulting 'ca_topo.tif' output file, which is 2.7 gigabytes"
|
||
|
||
[small_ca_topo]: small_ca_topo.png "a 'raw' heightmap of california and parts of nevada, arizona, and mexico"
|
||
|
||
[scaled_heightmap]: scaled_heightmap.png "the heightmap made by doing a linear mapping of height to brightness"
|
||
|
||
[^introspection]: The conclusion upon examination was, "I just wasn't thinking".
|
||
|
||
[^math-computers]: 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.
|
||
|
||
[^manifold_holes]: 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](https://pages.mtu.edu/~shene/COURSES/cs3621/SLIDES/Mesh.pdf) 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!
|
||
|
||
[^chekhovs-ram]: A textbook example of *Chekhov's Scarce Computational Resource*.
|
||
|
||
[^16-bit-ints]: Each pixel is 16 bits, so the possible values are from 0 to 2^16 - 1. 2^16 is 65536,
|
||
so there you go.
|
||
|
||
[^wgs-ellipsoid]: 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.
|
||
|
||
[^the-real-treasure-is-the-gd-treasure]: 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.
|