I released Rasterio 1.0.25 yesterday. It has a few important bug fixes, but the core of the work was writing and testing shims to make the package compatible with GDAL version 3.0 and PROJ version 6. Norman Barker did much of that work and I only had to make sure that we were using the right coordinate axis order strategy everywhere and figure out which of the output changes were new Rasterio bugs and which were actually improvements delivered by PROJ 6.

Please note that the binary wheels for 1.0.25 on PyPI contain GDAL 2.4.2, not 3.0, and that no new features of GDAL 3 and PROJ 6 are intentionally exposed in Rasterio's API. My wheel builds are already running up against time limits on Travis CI, and GDAL 3 and PROJ 6 take even longer to compile. My system is going to need some more hacks before Rasterio wheels with the latest GDAL and PROJ are possible. You might be able to get a combination of Rasterio 1.0.25, GDAL 3, and PROJ 6 from conda-forge. I look forward to hearing how that works for users.

I hope you'll appreciate that I've managed to shrink the size of the GDAL shared libraries in the manylinux wheels by 50%. I wish I knew why they are so much bigger than the OS X libraries. I suspect it's due to the ancient toolchain and glibc used by manylinux1.

In July the the Gnar Runners put on a race in State Forest State Park, the Never Summer 100K. It's a 64 mile loop through the Never Summer and Rawah ranges and the Colorado State Forest. The race requires park rangers, emergency medical technicians, ham radio operators, and lots of volunteer bodies.

Lulu and Thunder mountains, and moose, from Cameron Pass at 7 a.m.

In 2018 I flipped burgers at the finish line so that runners could get a hot meal in the middle of the night after being on the trail for twenty hours or more. In 2019, I spent an afternoon and evening supporting runners at the "Canadian" aid station at the race's 50 mile mark. The aid station gets its name from being near Never Summer Nordic's North Fork Canadian Yurt near the North Fork of the Canadian River, a tributary of the North Platte River. We had Canadian flags and a life-sized cutout of Justin Trudeau. At some point there was consensus that next year Ryan Gosling should join us.

Canadian aid station ready for runners to arrive, 2 p.m.

I had a cold, so I didn't cook or handle food, but I did a little of everything else. I hauled gear, deployed a portable toilet, set up canopies, organized and fetched drop bags, helped runners arrive, recover, and head out again. After the 1 a.m. cutoff, I helped break down the aid station and pack it all up. I was out on the course for 15 hours, 3 hours more than the first placed runner, but 8 hours less than the last finisher.

Last of the day's rain showers, 8 p.m.

We could see rain on the course as early as 11 a.m. None fell at the Canadian aid station until about 4 p.m., but then it rained until just before sunset. We had a pretty solid shower of small hail as well. Runners arriving in this rain were pretty worn out from a 6 mile slog through heavy mud. Some contemplated dropping out and a few did. Most found the energy to continue after some food, a cup of ramen or tea, and a couple minutes out of the rain. It was, after all, only 14 miles to the finish and there was plenty of time left.

Sunset and mud puddles

I felt good helping runners accomplish their goal and had a great time hanging out with other volunteers, many of them experienced ultra-runners, and listening to their stories. I wish I'd spent more time with the ham radio team and learned more about packet radio and running a network in the backcountry.

Race director Nick Clark's official recap of the race is here: http://gnarrunners.com/2019/08/a-recap-of-the-2019-never-summer-100km/.

I'm blogging this a day early because I expect to be too busy tomorrow. My Mom recently sent me her copy of the Detroit Free Press, her hometown paper, from Monday, July 21, 1969, the day of my birth.

I've enjoyed being connected to this milestone in space exploration and have been nuts about space for as long as I can remember. I wouldn't trade my birthday for any other.

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.

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.

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.

# crash.py

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 crash.py
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 https://eli.thegreenplace.net/2013/03/09/python-ffi-with-ctypes-and-cffi.

# 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 'libc.so.6'.
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()
187

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.

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.