diff --git a/content/sundries/a-very-digital-artifact/geotiff-files.png b/content/sundries/a-very-digital-artifact/geotiff-files.png index a7ef7d4..ae73afc 100644 Binary files a/content/sundries/a-very-digital-artifact/geotiff-files.png and b/content/sundries/a-very-digital-artifact/geotiff-files.png differ diff --git a/content/sundries/a-very-digital-artifact/index.md b/content/sundries/a-very-digital-artifact/index.md index 38ba049..2125708 100644 --- a/content/sundries/a-very-digital-artifact/index.md +++ b/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 -time, so I had had to download three or four separate -[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) +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 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] +*
last_little_piece_i_swear_final_final2.tif
* + +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] +*
if you squint, you can kinda see the mountains
* + +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 -the computer is just helping us do the math and it's more relatable. +[geotiff-files]: geotiff-files.png "the input geotiff files and the resulting 'ca_topo.tif' output file, which is 2.7 gigabytes" -[^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 and the edge of the hole hemmed or otherwise "sealed", which is then -a *manifold boundary*. See [this powerpoint +[small_ca_topo]: small_ca_topo.png "a 'raw' heightmap of california and parts of nevada, arizona, and mexico" + +[^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! @@ -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. diff --git a/content/sundries/a-very-digital-artifact/small_ca_topo.png b/content/sundries/a-very-digital-artifact/small_ca_topo.png new file mode 100644 index 0000000..bf06855 Binary files /dev/null and b/content/sundries/a-very-digital-artifact/small_ca_topo.png differ