checkpoint ca topo post

content/sundries/a-very-digital-artifact/index.md

@@ -133,26 +133,118 @@ second](https://gisgeography.com/srtm-shuttle-radar-topography-mission/) (which
 "GL**1**"), or roughly 30x30 meters, and the height data is accurate to within 16 meters. Not too
 shabby!
 
-The only problem was that you could only download data covering up to 450,000 square kilometers at a
+They provided the data in the form of [GeoTIFF](https://en.wikipedia.org/wiki/GeoTIFF)s, which are
-time, so I had had to download three or four separate
+basically an image where each pixel represents one data point (so, a 30x30 square meter plot)
-[GeoTIFF](https://en.wikipedia.org/wiki/GeoTIFF) files and then mosaic them together. A GeoTIFF file
-is 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 height is
 mapped to brightness, so the lowest spot's value is `0` (black), and the highest spot is
-`65535`[^16-bit-ints] (brightest white). These files are not small
+`65535`[^16-bit-ints] (brightest white).
 
-## Thanks, California state!
+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:
 
-https://data.ca.gov/dataset/ca-geographic-boundaries
+> 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
 
-## Give it a good smear
+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, that's an image that's 40,757 pixels wide and
+35,418 pixels tall. The "pixel size" is 0.000277777777778 by 0.000277777777778; since each pixel is
+1 arc second, and 1 arc second is 1/3600th of a degree, and 1/3600 is 0.000277777777..., we
+can infer that the unit of size there is 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, and relative
+to the [World Geodetic System 1984](https://en.wikipedia.org/wiki/World_Geodetic_System) datum; this
+refers to height; the very last line is the lowest and highest points in file, in meters from that WGS84
+baseline. 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, 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.
+
+## Thank you, State of California!
+
+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!
 
 # 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 domain 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"
@@ -161,14 +253,18 @@ https://data.ca.gov/dataset/ca-geographic-boundaries
 
 [meshy-cube]: meshy-cube.png "an overly-complicated mesh of a cube"
 
-[^math-computers]: I'm pretty sure this is more "represent shapes with math" than with a computer, but
+[geotiff-files]: geotiff-files.png "the input geotiff files and the resulting 'ca_topo.tif' output file, which is 2.7 gigabytes"
-the computer is just 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
+[small_ca_topo]: small_ca_topo.png "a 'raw' heightmap of california and parts of nevada, arizona, and mexico"
-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
+[^math-computers]: I'm pretty sure this is more "represent shapes with math" than with a computer, but
-cloth sheet with a hole cut out and the edge of the hole hemmed or otherwise "sealed", which is then
+the computer is helping us do the math and it's more relatable.
-a *manifold boundary*. See [this powerpoint
+
+[^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!
@@ -177,3 +273,6 @@ you). If I'm wrong about a 2D sheet with a hole being possibly manifold, I invit
 
 [^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 perfect sphere.