Running in the Cascades

Today was my first full day back in Colorado after a week in Washington with Ruth's clan. My phone stopped charging during the trip and so I didn't take as many pictures of the Cascades as I would have liked. I took a few during a short and steep run up a fire road near our rental house on Sunday and this is the best.

Looking down Bear Creek to Cle Elum Lake

I didn't carry my phone on my longest run into the Alpine Lakes Wilderness upstream from Cle Elum Lake and have to content myself with memories. According to Forest Service web pages, much of the area is covered by old growth, never logged, forest. Indeed, I saw many titanic Douglas firs around and above Pete Lake. Such giant trees are very rare in Colorado.

Offseason running

After slacking off for a couple weeks after Quad Rock, I've succeeded at getting back to 30+ miles per week. I'm running 3 times Monday through Friday and doing one long easy pace run on the weekend. Starting July 1, I'll start to build towards races in September and October.

I strained some muscle in my back 3 weeks ago and this has forced me to pay more attention to my running form. Engaging my glutes and hips and trying to get my upper body to float was a good way to minimize the pain in my back and is something I'm going to continue to practice.


This week will be busier than the last! I'm writing this from the neighborhood pool where Bea is training with her swim team. It's sunny and mild and it feels great to be writing outside.

Fix for Shapely's GEOS library loading bug in 1.7a2

Including GDAL and GEOS in binary wheels that go to the Python Package Index (PyPI) has been good for users of fiona, rasterio, and shapely, but has also exposed Python GIS programmers on OS X to a perplexing bug. Shapely's issue #553 has been weighing on me for a long time. It's been hard to understand, hard to explain, and bites very hard. The short program below will trigger the bug on any version of OS X with any version of fiona and any version of shapely before 1.7a2 which is installed using a wheel from PyPI.


import fiona
from shapely.geometry import Point
from shapely.ops import unary_union

print("Buffering points to create polygons")
SHAPES = [Point(i, i).buffer(1.5) for i in range(20)]

print("Computing union of polygons")
union = unary_union(SHAPES)

print("Union: {}".format(union))

For example:

$ python
Buffering points to create polygons
Computing union of polygons
Assertion failed: (!static_cast<bool>("should never be reached")), function itemsTree, file AbstractSTRtree.cpp, line 373.
Abort trap: 6

The unary_union function uses a GEOS STRtree and an assertion that should never be reached is in fact reached and the program aborts.

There has been a work-around – import shapely before importing other modules that are dynamically linked with GEOS, such as fiona and rasterio – but that's smelly and hard to ensure in a complex project.

The problem exists in a small niche: one operating system and one kind of program, one that that loads GEOS twice in different ways. Once when fiona's extension modules are loaded and cause the linked GEOS library to be loaded as with any compiled program, and once again when shapely is imported and it calls dlopen from python to load GEOS. I've found no diagnosis or solution of a similar problem on the web. Was the problem in GEOS? Was it in the library files that I built for the fiona, rasterio, and shapely wheels? Was the bug in delocate, the program I use to find dependencies and bundle them into the wheels? I asked GEOS developers and folks at work, to no avail. I was stumped for a long time.

I asked Even Rouault about the issue and he suggested a solution might be found in being more careful about using the correct mode with dlopen and he offered some example Python code that included this: handle = CDLL(None). My first thought was "can we use the ctypes CDLL class like this?" The docs say that the first argument is a path to a library file and that's how I've been using it in Shapely.

While looking for guidance on passing None to the CDLL constructor, I found some commented code in

# CDLL(None) invokes dlopen(NULL), which loads the currently running
# process - in our case Python itself. Since Python is linked with
# libc, readdir_r will be found there.
# Alternatively, we can just explicitly load ''.
lib = CDLL(None)

The OS X dlopen man page doesn't exactly say that,

dlopen() examines the mach-o file specified by path. If the file is compatible with the current process and has not already been loaded into the current process, it is loaded and linked. After being linked, if it contains any initializer functions, they are called, before dlopen() returns. dlopen() can load dynamic libraries and bundles. It returns a handle that can be used with dlsym() and dlclose(). A second call to dlopen() with the same path will return the same handle, but the internal reference count for the handle will be incremented. Therefore all dlopen() calls should be balanced with a dlclose() call.

If a null pointer is passed in path, dlopen() returns a handle equivalent to RTLD_DEFAULT.

but the POSIX spec is clear that dlopen(NULL) returns a handle to a global symbols object and that's what I observe on my MacBook. Here's an interpreter session where I import fiona and look in the global symbols object for a GDAL C API function.

>>> import fiona
>>> from ctypes import CDLL
>>> handle = CDLL(None)
>>> handle.GDALGetDriverCount
<_FuncPtr object at 0x111a3eb38>
>>> handle.GDALGetDriverCount()

When fiona is imported, Python calls dlopen with the path to fiona's extension (.so) modules, and as these are liked with libgdal, the GDAL library is loaded. Now, libgdal is linked with libgeos_c, so will we find GEOS C API functions in the global symbols object?

>>> handle.initGEOS
<_FuncPtr object at 0x111a3fb38>

Yes. This changes everything. The problem that has bedeviled shapely users is caused by using dlopen to load a library that's already been loaded. We can avoid the problem by detecting whether libgeos_c has already been loaded and skipping the troublesome dlopen call. I've done this for shapely and made a new 1.7 pre-release. Please considering running pip install -U --pre shapely to try it out.

The OS X man page for dlsym says that looking up functions in the global handle is the slowest possible approach, but ctypes does cache the addresses of functions it finds and my test programs don't run measurably slower, so this just might work out.

I'd like to understand why this problem doesn't also occur on Linux. Is it because of differences in the library loaders or the executable file formats themselves? Is there subtle platform-specific bug in GEOS? I'll write more if I get to the bottom of this.


Ruth is out of town this week and the next and on top of solo parenting I am dog sitting. I'll be extra busy through the end of June. Please forgive me if I don't respond immediately to emails or tickets on GitHub.

Arthur's Rock

Arthur's Rock is one of our local landmarks and hiking destinations. It's a big lump of pink Precambrian granite named after Arthur Howard, one of the previous owners of the land that would become Lory State Park in 1967. The Quad Rock and Black Squirrel courses go by it. I've been back twice since Quad Rock, once on a run and once with my family.

View south from Arthur's Rock Trail, May 26, 2019

May was cool and wet and this weather has continued to some degree in June. Lory State Park is as green as I've ever seen it.

Arthur's Rock from Howard Trail, June 1, 2019

The ponderosa pines in the second photo were charred during a small lightning-sparked fire in June 2016. I'm glad it wasn't worse. The small forest on the east slope of Lory State Park is priceless.

JavaScript mapping in 2018

It's been five-and-a-half years since I last blogged about making maps in a browser using JavaScript. It's cool that my old Syriac Places demo still works. Yay for web standards. A lot of things have changed in the meanwhile and this is my first chance to acknowledge some of those changes.

I've been writing HTML and JavaScript every other day at work for the past 3 weeks, building a web app to give Sales and sales-supporting teams at Mapbox insight into our imagery basemap. It's HTML and mostly Vanilla JS and has been a fun exercise. Since it's for internal customers only I can use new browser and JavaScript features.

Mapbox GL JS, the JavaScript mapping library that's maintained by my coworkers, is the closest thing to a framework in my app. My current project is the first in which I've used Mapbox GL JS. I'd like to point out some features and aspects of the framework that I find especially notable.

Mapbox GL JS has excellent examples. This is a result of I cribbed freely from and everything worked. Not all projects have good examples. John Firebaugh (who has since moved on from Mapbox), Colleen McGinnis, and the GL JS and Documentation teams put a ton of work into these and it shows.

Mapbox GL JS makes a big and important case of dynamic maps simple. I have a GeoJSON Source defined in my app's JavaScript. It is registered with a layer and a map. When I remove features from the source, they disappear from the map. When I add features to the source, they appear on the map. The map is updated fast enough that I can do this on every mousemove event.

Mapbox GL JS Expressions give a developer ample control over styling and filtering while keeping the styling API abstract and simple. In several parts of its API, Mapbox GL JS can evaluate Lisp-like expressions composed of JavaScript arrays. ["+", 4, ["*", 3, 2]] evaluates to 4 + (3 * 2), or 10. In addition to the standard JavaScript operators, the namespace of GL JS expressions operators contains feature property and map attribute getters. An advantages of using JavaScript objects instead of conventional Lisp (as with snuggs and rio-calc) is that developers get syntax highlighting help when writing expressions, and can use JavaScript to build expressions.

Another first for me is that I'm using Turf, the well-known simple features geometry library created by Morgan Herlocker. In my app I'm creating features on the fly derived from the intersections of features in a Mapbox tileset. Turf is fast and its functions take and return GeoJSON objects, which makes it easy for me to use and easy to integrate with Mapbox GL JS and other libraries or plugins. Turf is indispensable for a web app like the one I'm writing.

Lastly, JavaScript itself has changed a lot in the last 5 years. I'm using features of ES6 for the first time ever. You know this, I'll bet, but ES6 arrow functions are like Scheme's lambda expressions and Python's lambda forms. I'm using them exclusively, instead of the old anonymous functions, when filtering and sorting arrays of features. For example, the following sorts features by a date, in descending date order, without needing to type "function" or "return". Handy, and to my eyes, more readable.

features.sort((a, b) => -

These are obviously not the only recent advances in the JavaScript mapping field. Or even the most recent, which is why I chose the title "JavaScript mapping in 2018". The app that I'm building could have been written a year ago using the same software. For a glimpse at JavaScript mapping in 2019-2020, I suggest

Still snowing

The season's 19th snowfall is heavier than the 18th. We're now 24 days after the average last snowfall date 🤷.

May 21, 2019


I haven't run or hiked at all this week. Not a single step.

My strained calf is feeling much better after four days of rest and I'm looking forward to getting back out on the trail for an hour this weekend. My plan is to resume running around 25-30 miles a week, starting next week, and then start building up for races in September and October.