Shallow Thoughts

Finding Astronomical Alignments in Ancient Monuments (or anywhere else)

Dave and I will be presenting a free program on Stonehenge at the Los Alamos Nature Center tomorrow, June 14.

The nature center has a list of programs people have asked for, and Stonehenge came up as a topic in our quarterly meeting half a year ago. Remembering my seventh grade fascination with Stonehenge and its astronomical alignments -- I discovered Stonehenge Decoded at the local library, and built a desktop model showing the stones and their alignments -- I volunteered. But after some further reading, I realized that not all of those alignments are all they're cracked up to be and that there might not be much of astronomical interest to talk about, and I un-volunteered.

But after thinking about it for a bit, I realized that "not all they're cracked up to be" makes an interesting topic in itself. So in the next round of planning, I re-volunteered; the result is tomorrow night's presentation.

The talk will include a lot of history of Stonehenge and its construction, and a review of some other important or amusing henges around the world. But this article is on the astronomy, or lack thereof.

The Background: Stonehenge Decoded

Stonehenge famously aligns with the summer solstice sunrise, and that's when tens of thousands of people flock to Salisbury, UK to see the event. (I'm told that the rest of the time, the monument is fenced off so you can't get very close to it, though I've never had the opportunity to visit.)

Curiously, archaeological evidence suggests that the summer solstice wasn't the big time for prehistorical gatherings at Stonehenge; the time when it was most heavily used was the winter solstice, when there's a less obvious alignment in the other direction. But never mind that.

In 1963, Gerald Hawkins wrote an article in Nature, which he followed up two years later with a book entitled Stonehenge Decoded. Hawkins had access to an IBM 7090, capable of a then-impressive 100 Kflops (thousand floating point operations per second; compare a Raspberry Pi 3 at about 190 Mflops, or about a hundred Gflops for something like an Intel i5). It cost $2.9 million (nearly $20 million in today's dollars).

Using the 7090, Hawkins mapped the positions of all of Stonehenge's major stones, then looked for interesting alignments with the sun and moon. He found quite a few of them. (Hawkins and Fred Hoyle also had a theory about the fifty-six Aubrey holes being a lunar eclipse predictor, which captured my seventh-grade imagination but which most researchers today think was more likely just a coincidence.)

But I got to thinking ... Hawkins mapped at least 38 stones if you don't count the Aubrey holes. If you take 38 randomly distributed points, what are the chances that you'll find interesting astronomical alignments?

A Modern Re-Creation of Hawkins' Work

Programmers today have it a lot easier than Hawkins did. We have languages like Python, with libraries like PyEphem to handle the astronomical calculations. And it doesn't hurt that our computers are about a million times faster.

Anyway, my script, skyalignments.py takes a GPX file containing a list of geographic coordinates and compares those points to sunrise and sunset at the equinoxes and solstices, as well as the full moonrise and moonset nearest the solstice or equinox. It can find alignments among all the points in the GPX file, or from a specified "observer" point to each point in the file. It allows a slop of a few degrees, 2 degrees by default; this is about four times the diameter of the sun or moon, but a half-step from your observing position can make a bigger difference than that. I don't know how much slop Hawkins used; I'd love to see his code.

My first thought was, what if you stand on a mountain peak and look around you at other mountain peaks? (It's easy to get GPS coordinates for peaks; if you can't find them online you can click on them on a map.) So I plotted the major peaks in the Jemez and Sangre de Cristo mountains that I figured were all mutually visible. It came to 22 points; about half what Hawkins was working with.

My program found (114 alignments.

Yikes! Way too many. What if I cut it down? So I tried eliminating all but the really obvious ones, the ones you really notice from across the valley. The most prominent 11 peaks: 5 in the Jemez, 6 in the Sangres.

That was a little more manageable. Now I was down to only 22 alignments.

Now, I'm pretty sure that the Ancient Ones -- or aliens -- didn't lay out the Jemez and Sangre de Cristo mountains to align with the rising and setting sun and moon. No, what this tells us is that pretty much any distribution of points will give you a bunch of astronomical alignments.

And that's just the sun and moon, all Hawkins was considering. If you look for writing on astronomical alignments in ancient monuments, you'll find all people claiming to have found alignments with all sorts of other rising and setting bodies, like Sirius and Orion's belt. Imagine how many alignments I could have found if I'd included the hundred brightest stars.

So I'm not convinced. Certainly Stonehenge's solstice alignment looks real; I'm not disputing that. And there are lots of other archaeoastronomy sites that are even more convincing, like the Chaco sun dagger. But I've also seen plenty of web pages, and plenty of talks, where someone maps out a collection of points at an ancient site and uses alignments among them as proof that it was an ancient observatory. I suspect most of those alignments are more evidence of random chance and wishful thinking than archeoastronomy.

Tags: astro, programming, python, mapping, speaking
[ 14:54 Jun 13, 2019    More science/astro | permalink to this entry | comments ]

Styling GTK3 in Python with CSS

Lately I've been running with my default python set to Python 3. Debian still uses Python 2 as the default, which is reasonable, but adding a ~/bin/python symlink to /usr/bin/python3 helps me preview scripts that might become a problem once Debian does switch. I thought I had converted most of my Python scripts to Python 3 already, but this link is catching some I didn't convert.

Python has a nice script called 2to3 that can convert the bulk of most scripts with little fanfare. The biggest hassles that 2to3 can't handle are network related (urllib and urllib2) and, the big one, user interfaces. PyGTK, based on GTK2 has no Python 3 equivalent; in Python 3, the only option is to use GObject Introspection (gi) and GTK3. Since there's almost no documentation on python-gi and gtk3, converting a GTK script always involves a lot of fumbling and guesswork.

A few days ago I tried to play an MP3 in my little musicplayer.py script and discovered I'd never updated it. I have enough gi/GTK3 scripts by now that I thought something with such a simple user interface would be easy. Shows how much I know about GTK3!

I got the basic window ported pretty easily, but it looked terrible: huge margins everywhere, and no styling on the text, like the bold, large-sized text I had previously use to highlight the name of the currently playing song. I tried various approaches, but a lot of the old methods of styling have been deprecated in GTK3; you're supposed to use CSS. Except, of course, there's no documentation on it, and it turns out the CSS accepted by GTK3 is a tiny subset of the CSS you can use in HTML pages, but what the subset is doesn't seem to be documented anywhere.

How to Apply a Stylesheet

The first task was to get any CSS at all working. The GNOME Journal: Styling GTK with CSS was helpful in getting started, but had a lot of information that doesn't work (perhaps it did once). At least it gave me this basic snippet:

    css = '* { background-color: #f00; }'
    css_provider = gtk.CssProvider()
    css_provider.load_from_data(css)
    context = gtk.StyleContext()
    screen = Gdk.Screen.get_default()
    context.add_provider_for_screen(screen, css_provider,
                                    gtk.STYLE_PROVIDER_PRIORITY_APPLICATION)

Built-in Class Names

Great! if all you want to do is turn the whole app red. But in reality, you'll want to style different widgets differently. At least some classes have class names:

    css = 'button { background-color: #f00; }'

I found other pages suggesting using 'GtkButton in CSS, but that didn't work for me. How do you find the right class names? No idea, I never found a reference for that. Just guess, I guess.

User-set Class Names

What about classes -- for instance, make all the buttons in a ButtonBox white? You can add classes this way:

    button_context = button.get_style_context()
    button_context.add_class("whitebutton")

If you need to change a class (for instance, turn a red button green), first remove the old class:

    button_context = button.get_style_context()
    entry_style_context.remove_class("red")

Widget Names, like CSS ID

For single widgets, you can give the widget a name and use it like an ID in CSS. Like this:

    label = gtk.Label()
    label.set_use_markup(True)
    label.set_line_wrap(True)
    label.set_name("red_label")
    mainbox.pack_start(label, False, False, 0)
    css = '#red_label { background-color: #f00; }'
[ ... ]

Properties You Can Set

There is, amazingly, a page on which CSS properties GTK3 supports. That page doesn't mention it, but some properties like :hover are also supported. So you can write CSS tweaks like

.button { border-radius: 15; border-width: 2; border-style: outset; }
.button:hover { background: #dff; border-color: #8bb; }

And descendants work, so you can say somthing like

    buttonbox = gtk.ButtonBox(spacing=4)
    buttonbox.set_name("buttonbox")
    mainbox.pack_end(buttonbox, False, False, 0)

    btn = gtk.Button(label="A")
    buttonbox.add(btn)

    btn = gtk.Button(label="B")
    buttonbox.add(btn)

and then use CSS that affects all the buttons inside the buttonbox:

#buttonbox button { color: red; }

No mixed CSS Inside Labels

My biggest disappointment was that I couldn't mix styles inside a label. You can't do something like

label.set_label('Headline'
                'Normal text')

and expect to style the different parts separately. You can use very simple markup like <b>bold</b> normal, but anything further gives errors like "error parsing markup: Attribute 'class' is not allowed on the <span> tag" (you'll get the same error if you try "id"). I had to make separate GtkLabels for each text size and style I wanted, which is a lot more work. If you wanted to mix styles and have them reflow as the content length changed, I don't know how (or if) you could do it.

Fortunately, I don't strictly need that for this little app. So for now, I'm happy to have gotten this much working.

Tags: programming, python, gtk, gtk3
[ 14:49 Jun 05, 2019    More programming | permalink to this entry | comments ]

Plotting a Sequence of Graphs in Matplotlib 3D

A friend and I were talking about temperature curves: specifically, the way the temperature sinks in the evening but then frequently rises again before it really starts cooling off.

I thought it would be fun to plot the curve of temperature as a function of time over successive days, as a 3-D plot. I knew matplotlib had a way to do 3D plots, but I've never actually generated one.

Well, it turns out there are lots of examples, but they all start by generating mysterious data blobs, and none of them explain the structure of the data they're using, and the documentation has mysterious parameters like "zs" that aren't explained anywhere. So getting something that worked was a fiddly process. Creating a color version, to distinguish the graphs better, was even more fiddly.

So I wrote an example that I hope will make it a little clearer for anyone trying to use this library. It can plot using just lines:

... or it can plot in color, cycling colors manually because by default matplotlib makes adjacent colors similar, exactly the opposite of what you'd want:

Here's the demo: multiplot3d.py on GitHub.

... Except there's a Bug

All is not perfect. Axes3D gets a bit confused sometimes about which layer is supposed to be in front of which other layer. You can see that on the two plots: in both cases, the fourth and fifth layers from the front are reversed, so the fifth layer is drawn in front of the fourth layer. I haven't yet found anyone in the matplotlib organization who seems to know much about Axes3D; eventually I'll file a bug but I want to write a shorter, clearer test case to illustrate the problem. Still, even with the bugs it's a useful technique to know.

Tags: python, programming, data, matplotlib
[ 09:57 May 30, 2019    More programming | permalink to this entry | comments ]

Making a Land Ownership overlay: Categorized Styles in QGIS

Now that I know how to make a map overlay for OsmAnd, I wanted a land ownership overlay. When we're hiking, we often wonder whether we're on Forest Service, BLM, or NPS land, or private land, or Indian land. It's not easy to tell.

Finding Land Ownership Data

The first trick was finding the data. The New Mexico State Land Office has an interactive New Mexico Land Status map, but that's no help when walking around, and their downloadable GIS files only cover the lands administered by the state land office, which mostly doesn't include any areas where we hike. They do have some detailed PDF maps of New Mexico Lands if you have a printer capable of printing enormous pages, which most of us don't.

In theory I could download their 11" x 17" Land Status PDF, convert it to a raster file, and georeference it as I described in the earlier article; but since they obviously have the GIS data (used for the interactive map) I'd much rather download the data and save myself all that extra work.

Eventually I found New Mexico ownership data at UNM's RGIS page, which has an excellent collection of GIS data available for download. Click on Boundaries, then download Surface Land Ownership. It's available in a variety of formats; I chose the geojson format because I find it the most readable and the easiest to parse with Python, though ESRI shapefiles arguably might have been easier in QGIS.

Colorizing Polygons in QGIS

You can run qgis on a geojson file directly. When it loads it shows the boundaries, and you can use the Info tool to click on a polygon and see its metadata -- ownership might be BLM, DOE, FS, I, or whatever. But they're all the same color, so it's hard to get a sense of land ownership just clicking around.

To colorize the polygons differently, right-click on the layer name and choose Properties. For Style, choose Categorized. For Column, pick the attribute you want to use to choose colors: for this dataset, it's "own", for ownership.

Color ramp is initially set to random. Click Classify to generate an initial color ramp, then click Apply to see what it looks like on the map.

Then you can customize the colors by doubleclicking on specific color swatches. For instance, by unstated convention most maps show Forest Service land as green, BLM and Indian land as various shades of brown. Click Apply as you change colors, until you're happy with the result.

Exporting to GeoTIFF

You can export the colored layer to GeoTIFF using QGIS' confusing and poorly documented Print Composer. Create one with: Project > New Print Composer, which will open with a blank white canvas.

Zoom and pan in the QGIS window so the full extent of the image you want to export is visible. Then, in the Print Composer, Layout > Add Map. Click and drag in the blank canvas, going from one corner to the opposite corner, and some portion of the map should appear.

There doesn't seem to be any way to Print Composer to import your whole map automatically, or for you to control what portion of the map from the QGIS window will show up in the Print Composer when you drag. If you guess wrong and don't get all of your map, hit Delete, switch to the QGIS window and drag and/or zoom your map a little, then switch back to Print Composer and try adding it again.

You can also make adjustments by changing the Extents in the Item Properties tab, and clicking the Set to map canvas extent button in that tab will enlarge your extents to cover approximately what's currently showing in the QGIS window.

It's a fiddly process and there's not much control, but when you decide it's close enough, Composer > Export as Image... and choose TIFF format. (Print Composer offers both TIFF and TIF; I don't know if there's a difference. I only tried TIFF with two effs.) That should write a GeoTIFF format; to verify that, go to a terminal and run gdalinfo on the saved TIFF file and make sure it says it's GeoTIFF.

Load into OsmAnd

Finally, load the image into OsmAnd's tiles folder as discussed in the previous article, then bring up the Configure map menu and enable the overlay.

I found that the black lines dividing the various pieces of land are a bit thicker than I'd like. You can't get that super accurate "I'm standing with one foot in USFS land and the other foot in BLM land" feeling because of the thick black DMZ dividing them. But that's probably just as well: I suspect the data doesn't have pinpoint accuracy either. I'm sure there's a way to reduce the thickness of the black line or eliminate it entirely, but for now, I'm happy with what I have.

Update: Here's another, easier, way to show land use on OsmAnd using overlay tiles from the BLM (in the US): Adding BLM Land Use Maps to Osmand on Android. It isn't as general (you can only show something you can get from an online tiled source) and it updates in real-time, meaning it might use cellphone data rather than working entirely offline, but it's still a great option to know about.

Tags: mapping, qgis, osmand
[ 18:13 Apr 15, 2019    More mapping | permalink to this entry | comments ]

Making Overlay Maps for OsmAnd on Linux

For many years I've wished I could take a raster map image, like a geology map, an old historical map, or a trail map, and overlay it onto the map shown in OsmAnd so I can use it on my phone while walking around. I've tried many times, but there are so many steps and I never found a method that worked.

Last week, the ever helpful Bart Eisenberg posted to the OsmAnd list a video he'd made: Displaying web-based maps with MAPC2MAPC: OsmAnd Maps & Navigation. Bart makes great videos ... but in this case, MAPC2MAPC turned out to be a Windows program so it's no help to a Linux user. Darn!

But seeing his steps laid out inspired me to try again, and gave me some useful terms for web searching. And this time I finally succeeded. I was also helped by a post to the OsmAnd list by A Thompson, How to get aerial image into offline use?, though I needed to change a few of the steps. (Note: click on any of the screenshots here to see a larger version.)

Georeference the Image Using QGIS

The first step is to georeference the image: turn the plain raster image into a GeoTiff that has references showing where on Earth its corners are. It turns out there's an open source program that can do that, QGIS. Although it's poorly documented, it's fairly easy once you figure out the trick.

I started with the tutorial Georeferencing Basics, but it omits one important point, which I finally found in BBRHUFT's How to Georeference a map in QGIS. Step 11 is the key: the Coordinate Reference System (CRS) must be the same in the georeferencing window as it is in the main QGIS window. That sounds like a no-brainer, but in practice, the lists of possible CRSes shown in the two windows don't overlap, so unless you follow BBRHUFT's advice and type 3857 into the filter box in both windows, you'll likely end up with CRSes that don't match. It'll look like it's working, but the resulting GeoTiff will have coordinates nowhere near where they should be

Instead, follow BBRHUFT's advice and type 3857 into the filter box in both windows. The "WGS 84 / Pseudo Mercator" CRS will show up and you can use it in both places. Then the GeoTiff will come out in the right place.

If you're starting from a PDF, you may need to convert it to a raster format like PNG or JPG first. GIMP can do that.

So, the full QGIS steps are:

• Start qgis, and bring up Project > Project Properties.
• Click on CRS, type 3857 in the filter box, and make sure the CRS is set to WGS 84 / Pseudo Mercator.
• (optional) While you're in Project Properties, you may want to click on General and set Display Coordinates Using: Decimal degrees. By default qgis always uses "meters" (measured from where?) which seems like a bizarre and un-useful default.
• In Plugins > Manage and Install Plugins..., install two plugins: Quick Map Services and Georeferencer GDAL.
• Back in the main QGIS window, go to Web > QuickMapServices (that's one of the plugins you just installed) and pick a background map, probably something from OpenStreetMap. Once you find one that loads, navigate to the area you're interested in.
  
  (All the OpenStreetMap tile sets print "API Key Required" all over everything. You can get an API key from Thunderforest -- I use one for PyTopo -- but it won't help you here, because nobody seems to know how to get QGIS to use an API key, though I found dozens of threads with people asking. If you find a way, let me know.)
  
• Start the georeferencer: Raster > Georeferencer > Georeferencer.... It may prompt you for a CRS; if so, choose WGS 84 / Pseudo Mercator if it's available, and if it's not, type 3857 in the filter box to find it.
• In the Georeferencer window, File > Open Raster...
  
• Select matching points. This part is time consuming but fun.
  
  Zoom in to your raster map in the Georeferencer window and pick a point that's easy to identify: ideally the intersection of two roads, or a place where a road crosses a river, or similar identifying features. Make sure you can identify the same point in both the Georeferencer window and the main QGIS window (see the yellow arrows in the diagram).
  
• Click on the point in the Georeferencer window, then click "from map canvas".
• In the main QGIS window, click on the corresponding point. You can zoom with the mousewheel and pan by middle mouse dragging. Repeat for a total of four to eight points.
  
  You'll see yellow boxes pop up in both windows showing the matching points. If you don't, there's something wrong; probably your CRSes don't match.
  
• In the Georeferencer window, call up Settings > Transformation Settings. Choose a Transformation type (most people seem to like Thin Plate Spline), a Resampling method (Lanczos is good), a target SRS (make sure it's 3857 WGS 84 / Pseudo Mercator to match the main window). Set Output raster to your desired filename, ending in .tiff. Be sure to check the "Load in QGIS when done" box at the bottom. Click OK.
  
• File > Start Georeferencing. It's quite fast, and will create your .tiff file when it's finished, and should load the TIFF as a new layer in the main QGIS window. Check it and make sure it seems to be lined up properly with geographic features. If you need to you can right-click on the layer and adjust its transparency.
• In the Georeferencer window, File > Save GCP Points as... in case you want to go back and do this again. Now you can close the Georeferencer window.
• Project > Save to save the QGIS project, so if you want to try again you don't have to repeat all the steps.
  
  Now you can close qgis. Whew!

Convert the GeoTiff to Map Tiles

The ultimate goal is to convert to OsmAnd's sqlite format, but there's no way to get there directly. First you have to convert it to map tiles in a format called mbtiles.

QGIS has a plug-in called QTiles but it didn't work for me: it briefly displayed a progress bar which then disappeared without creating any files. Fortunately, you can do the conversion much more easily with gdal_translate, which at least on Debian is part of the gdal-bin package.

gdal_translate filename.tiff filename.mbtiles

That will create tiles for a limited range of zoom levels (maybe only one zoom level). gdalinfo will tell you the zoom levels in the file. If you want to be able to zoom out and still see your overlay, you might want to add wider zoom levels, which you can do like this:

gdaladdo -r nearest filename.mbtiles 2 4 8 16

Incidentally, gdal can also create a directory of tiles suitable for a web slippy map, though you don't need that for OsmAnd. For that, use gdal2tiles, which on Debian is part of the python-gdal package:

mkdir tiles
gdal2tiles filename.tiff tiles

Not only does it create tiles, it also includes multiple HTML files you can use to display those tiles using the Leaflet, OpenLayers or Google Maps JavaScript libraries. Very cool!

Create the OsmAnd sqlite file

Tarwirdur has written a nice simple Python script to translate from mbtiles to OsmAnd sqlite: mbtiles2osmand.py. Download it then run

mbtiles2osmand.py filename.mbtiles filename.sqlitedb

So easy to use! Most of the other references I saw said to use Mobile Atlas Creator (MOBAC) and that looked a lot more complicated.

Incidentally, Bart's video says MAPC2MAPC calls the format "Locus/Rmaps/Galileo/OSMAND (sqlite)", which might be useful to know for web search purposes.

Install in OsmAnd

Once you have the .sqlitedb file, copy it to OsmAnd's tiles folder in whatever way you prefer. For me, that's adb push file.sqlitedb $androidSD/Android/data/net.osmand.plus/files/tiles where $androidSD is the /storage/whatever location of my device's SD card.

Then start OsmAnd and tap on the icon in the upper left for your current mode (car, bike, walking etc.) to bring up the Configure map menu. Scroll down to Overlay or Underlay map, enable one of those two and you should be able to choose your newly installed map.

You can adjust the overlay's transparency with a slider that's visible at the bottom of the map (the blue slider just above the distance scale), so you can see your overlay and the main map at the same time.

The overlay disappears if you zoom out too far, and I haven't yet figured out what controls that; I'm still working on those details.

Sure, this process is a lot of work. But the result is worth it. Check out the geologic layers we walked through on a portion of a recent hike in Rendija Canyon (our hike is the purple path).

Tags: mapping, linux, osmand, qgis
[ 19:08 Apr 10, 2019    More mapping | permalink to this entry | comments ]

Antler in the Yard

We see mule deer often enough that I've wondered if we would ever find a shed antler in the yard.

A couple of days ago Dave found one. It's just a 4-pointer, but it's a big 4-pointer. We're still hoping its mate might be somewhere nearby, but no luck so far.

It feels enormously heavy, though the scale says it's just shy of two pounds. Still, carrying around four pounds on your head all the time ... sounds like a recipe for a headache.

The day after the antler turned up, five bucks visited our garden in the evening. Four of them still had substantial buttons where antlers had recently been. The fifth still had antlers (a ten-pointer).

I have to wonder, what do they think about that? Is the guy who still has antlers the macho king of the yard? Or are the other four saying "Gosh, I'm so glad mine dropped last week, so sorry you still have to carry that rack around"?

Tags: nature
[ 15:51 Mar 28, 2019    More nature | permalink to this entry | comments ]

Natural Topiary

Some months after we moved into our new house, we came home to notice a really ugly job of tree trimming on some junipers in the driveway.

We hadn't done it, nor had we requested tree trimming service. Yet we'd been in the house for months, and we were both pretty sure that the junipers hadn't looked like that before. They were nearly denuded at the bottom, belling out higher up.

Was there some rogue topiarist wandering the neighborhood, defoliating people's trees without asking them? When I was growing up, occasionally my grandmother would show up and slash branches off trees in our backyard without asking. But it made no sense that anyone in our neighborhood would do that.

We hoped the trees would grow back to their former bushiness, but several years later, they don't look much better. And all that time it has remained a mystery how the trees came to look like that.

Until a party a few months ago, when a visiting friend cast a knowing look at the trees and said, "I see the deer have been visiting you."

Of course! Somehow that had never occurred to us. We went out to the driveway and checked -- and sure enough, the trimmed parts of the trees go up to roughly the height a deer could easily reach.

So sure enough -- there is a rogue topiarist wandering the neighborhood after all. Lots of them, in fact. Most of the time they eat more tempting fare; but when the weather gets inclement and there isn't much to eat, they'll go after the junipers. And I guess these non-native junipers in the driveway are a little tastier than the native ones that are all around the property.

We don't feel quite so resentful about our unwanted tree trimming now. Sure, we're still not crazy about the look of our oddly shaped junipers, and they've gotten even worse during the current exceptionally snowy winter; but now that we know it's a natural process, not some crazed shear-wielding neighbor, it's hard to be too upset by it.

[ 11:35 Feb 23, 2019    More nature | permalink to this entry | comments ]

Reading the Output of a Weather Station Using Software Defined Radio

A while back, Dave ordered a weather station. His research pointed to the Ambient Weather WS-2000 as the best bang for the buck as far as accuracy (after it's calibrated, which is a time consuming and exacting process where you compare readings to a known-good mercury thermometer, a process that I suspect most weather station owners don't bother with).

It comes with a little 7" display console that sits indoors and reads the radio signal from the outside station as well as a second thermometer inside, then displays all the current weather data. It also uses wi-fi to report the data upstream to Ambient and, optionally, to a weather site such as Wunderground. (Which we did for a while, but now Wunderground is closing off their public API, so why give them data if they're not going to make it easy to share it?)

Having the console readout and the Ambient "dashboard" is all very nice, but of course, being a data geek, I wanted a way to get the data myself, so I could plot it, save it or otherwise process it. And that's where Ambient falls short. The console, though it's already talking on wi-fi, gives you no way to get the data. They sell a separate unit called an "Observer" that provides a web page you can scrape, and we actually ordered one, but it turned out to be buggy and difficult to use, giving numbers that were substantially different from what the console showed, and randomly failing to answer, and we ended up returning the observer for a refund.

The other way of getting the data is online. Ambient provides an API you can use for that purpose, if you email them for a key. It mostly works, but it sometimes lags considerably behind real time, and it seems crazy to have to beg for a key and then get data from a company website that originated in our own backyard.

What I really wanted to do was read the signal from the weather station directly. I'd planned for ages to look into how to do that, but I'm a complete newbie to software defined radio and wasn't sure where to start. Then one day I noticed an SDR discussion on the #raspberrypi IRC channel on Freenode where I often hang out. I jumped in, asked some questions, and got pointed in the right direction and referred to the friendly and helpful #rtlsdr Freenode channel.

An Inexpensive SDR Dongle

On the experts' advice, I ordered a RTL-SDR Blog R820T2 RTL2832U 1PPM TCXO SMA Software Defined Radio with 2x Telescopic Antennas on Amazon. This dongle apparently has better temperature compensation than cheaper alternatives, it came with a couple of different antenna options, and I was told it should work well with Linux using a program called rtl_433.

Indeed it did. The command to monitor the weather station is

rtl_433 -f 915M

rtl_433 already knows the protocol for the WS-2000, so I didn't even need to do any decoding or reverse engineering; it produces a running account of the periodic signals being broadcast from the station. rtl_433 also helpfully offers -F json and -F csv options, along with a few other formats. What a great program!

JSON turned out to be the easiest for me to use; initially I thought CSV would be more compact, but rtl_433's CSV format includes fields for every possible quantity a weather station could ever broadcast. When you think about it, that makes sense: once you're outputting CSV you can't add a new field in mid-stream, so you'd better be ready for anything. JSON, on the other hand, lets you report just whatever the weather station reports, and it's easy to parse from nearly any programming language.

Testing the SDR Dongle

Full disclosure: at first, rtl_433 -f 915M wasn't showing me anything and I still don't know why. Maybe I had a loose connection on the antenna, or maybe I got unlucky and the weather station picked the exact wrong time to take a vacation. But while I was testing, I found another program that was very helpful in testing whether my SDR dongle was working: rtl_fm, which plays radio stations. The only trick is finding the right arguments, since the example from the man page just played static. Here's what worked for me:

rtl_fm -f 101.1M -M fm -g 20 -s 200k -A fast -r 32k -l 0 -E deemp | play -r 32k -t raw -e s -b 16 -c 1 -V1 -

That command plays the 101.1 FM radio station. (I had to do a web search to give me some frequencies of local radio stations; it's been a long time since I listened to normal radio.)

Once I knew the dongle was working, I needed to verify what frequency the weather station was using for its broadcasts. What I really wanted was something that would scan frequencies around 915M and tell me what it found. Everyone kept pointing me to a program called Gqrx. But it turns out Gqrx on Linux requires PulseAudio and absolutely refuses to work or install without it, even if you have no interest in playing audio. I didn't want to break my system's sound (I've never managed to get sound working reliably under PulseAudio), and although it's supposedly possible to build Gqrx without Pulse, it's a difficult build: I saw plenty of horror stories, and it requires Boost, always a sign that a build will be problematic. I fiddled with it a little but decided it wasn't a good time investment.

I eventually found a scanner that worked: RTLSDR-Scanner. It let me set limiting frequencies and scan between them, and by setting it to accumulate, I was able to verify that indeed, the weather station (or something else!) was sending a signal on 915 MHz. I guess by then, the original problem had fixed itself, and after that, rtl_433 started showing me signals from the weather station. It's not super polished, but it's the only scanner I've found that works without requiring PulseAudio.

That Puzzling Rainfall Number

One mystery remained to be solved. The JSON I was getting from the weather station looked like this (I've reformatted it for readablility):

{
    "time" : "2019-01-11 11:50:12",
    "model" : "Fine Offset WH65B",
    "id" : 60,
    "temperature_C" : 2.200,
    "humidity" : 94,
    "wind_dir_deg" : 316,
    "wind_speed_ms" : 0.064,
    "gust_speed_ms" : 0.510,
    "rainfall_mm" : 90.678,
    "uv" : 324,
    "uvi" : 0,
    "light_lux" : 19344.000,
     "battery" : "OK",
    "mic" : "CRC"
}

This on a day when it hadn't rained in ages. What was up with that "rainfall_mm" : 90.678 ?

I asked on the rtl_433 list and got a prompt and helpful answer: it's a cumulative number, since some unspecified time in the past (possibly the last time the battery was changed?) So as long as I make a note of the rainfall_mm number, any change in that number means new rainfall.

This being a snowy winter, I haven't been able to test that yet: the WS-2000 doesn't measure snowfall unless some snow happens to melt in the rain cup.

Some of the other numbers, like uv and uvi, are in mysterious unknown units and sometimes give results that make no sense (why doesn't uv go to zero at night? You're telling me that there's that much UV in starlight?), but I already knew that was an issue with the Ambient. It's not rtl_433's fault.

I notice that the numbers are often a bit different from what the Ambient API reports; apparently they do some massaging of the numbers, and the console has its own adjustment factors too. We'll have to do some more calibration with a mercury thermometer to see which set of numbers is right.

Anyway, cool stuff! It took no time at all to write a simple client for my WatchWeather web app that runs rtl_433 and monitors the JSON output. I already had WatchWeather clients collecting reports from Raspberry Pi Zero Ws sitting at various places in the house with temperature/humidity sensors attached; and now my WatchWeather page can include the weather station itself.

Meanwhile, we donated another weather station to the Los Alamos Nature Center, though it doesn't have the SDR dongle, just the normal Ambient console reporting to Wunderground.

Tags: SDR, weather, programming, raspberry pi
[ 13:20 Feb 12, 2019    More tech | permalink to this entry | comments ]