nebkor/blog
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

checkpoint

Browse source
This commit is contained in:
Joe Ardent 2023-01-13 17:59:22 -08:00
parent 845b336dd8
commit 212bc2e226
1 changed files with 181 additions and 12 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

193
content/sundries/a-very-digital-artifact/index.md
View file

@ -14,7 +14,8 @@ Last summer, I wanted to get my wife something nice for her birthday. For many y
expressed an occasional and casual desire for a topographic carving of the state of California, 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, 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 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. 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 # First steps
@ -32,7 +33,7 @@ available and wanted to give broad acceptable parameters), and under "project de
> A relief map of California, carved from wood. Height exaggerated enough > 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. to visibly discern the Santa Monica mountains. I can provide an STL file if needed.
For some [incorrect] reason that I only later examined, I just sort of assumed that the shop would 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 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 need. But just in case, I included that hedge about being able to provide an STL file. Needless to
say, that was a bluff. say, that was a bluff.
@ -101,8 +102,8 @@ But so far, I had nothing at all. Time to get some data and see if I can turn it
My first impulse was to search [USGS](https://usgs.gov)'s website for 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 [heightmap](https://en.wikipedia.org/wiki/Heightmap) data, but I wound up not finding anything
appropriate. Once again, now that I'm looking after I'm done, I found this, which would have been appropriate. Searching now with the wisdom of experience and hindsight, I found this, which would
perfect: have been perfect:
[https://apps.nationalmap.gov/downloader/](https://apps.nationalmap.gov/downloader/) [https://apps.nationalmap.gov/downloader/](https://apps.nationalmap.gov/downloader/)
@ -203,9 +204,9 @@ fill a rectangle that fully enclosed the entire border of California, then one b
four-hundred-forty-three million, five-hundred-thirty-one thousand, and four-hundred-twenty-six four-hundred-forty-three million, five-hundred-thirty-one thousand, and four-hundred-twenty-six
(40,757 times 35,418) is pretty close. (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 The other units in there are under the "Coordinate System is" section, and are meters relative to
to the [World Geodetic System 1984](https://en.wikipedia.org/wiki/World_Geodetic_System) datum; this the [World Geodetic System 1984](https://en.wikipedia.org/wiki/World_Geodetic_System) vertical
refers to height; the very last line is the lowest and highest points in file, in meters from that WGS84 datum; 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: 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] ![the ca_topo image; it's hard to make out details and very dark][small_ca_topo]
@ -216,13 +217,179 @@ the highest point in our dataset is only 4,412, which is not that much in compar
includes portions of not-California in the height data, and ideally, we want those places to be 0; 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. we have a little more processing to do before we can use this.
## Thank you, State of California! ## Cartography is complicated
The first order of business is to mask out everything that's not California, and the first thing I 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 needed for that was a [shapefile](https://en.wikipedia.org/wiki/Shapefile) that described the
California state border. Luckily, [that exact California state border. Luckily, [that exact
thing](https://data.ca.gov/dataset/ca-geographic-boundaries) is publicly available from the state's thing](https://data.ca.gov/dataset/ca-geographic-boundaries) is publicly available from the state's
website! 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 and it's still almost the same size (though it's slightly different now, with the
different projection); scaling it down is a last step, since at that point, it will no longer be a
digital elevation map, it will just be an image. But we'll get there.
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...." since the "lengthunit" is "metre". But the number of pixels is different; it's not
as wide, yet the coordinates of the bounding corners are the same as the original file (the latitude
and longitude given as the second tuple).
## The one custom script
So, the next step was use our shapefile to mask out the California border in our 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 something, it's the dumbest, quickest thing I crapped
out to do the task I needed to be done.
## A usable heightmap
# A mesh is born
# Test prints # Test prints
@ -235,7 +402,7 @@ website!
thank you: wife, steve at the shop, friends indulging my oversharing, NASA, Blender, FreeCAD 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 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 much better than when I was first racing through the material trying to get to a point where I could
just produce the artifact. just produce the artifact.
lesson: pipeline of geotiff -> mask -> scaled heightmap -> mesh -> solid body lesson: pipeline of geotiff -> mask -> scaled heightmap -> mesh -> solid body
@ -257,6 +424,8 @@ given how much I had to blur and decimate)?
[small_ca_topo]: small_ca_topo.png "a 'raw' heightmap of california and parts of nevada, arizona, and mexico" [small_ca_topo]: small_ca_topo.png "a 'raw' heightmap of california and parts of nevada, arizona, and mexico"
[^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 [^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. the computer is helping us do the math and it's more relatable.
@ -269,10 +438,10 @@ deck](https://pages.mtu.edu/~shene/COURSES/cs3621/SLIDES/Mesh.pdf) for a pretty
"mesh basics" (but not really that basic, that's just academics trolling us, don't let it bother "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! you). If I'm wrong about a 2D sheet with a hole being possibly manifold, I invite correction!
[^chekhovs-ram]: A classic example of Chekhov's Scarce Computational Resource. [^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, [^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. so there you go.
[^wgs-ellipsoid]: Technically, it's an arc along the WGS84 ellipsoid, which is a perfectly smooth [^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. *smushed* sphere, which more closely matches the real shape of the Earth vs. a perfectly round sphere.