Fetching OpenStreetMap Details with OSMPythonTools (Shallow Thoughts)

Fetching OpenStreetMap Details with OSMPythonTools

I was talking to a friend about LANL's proposed new powerline. A lot of people are opposing it because the line would run through the Caja del Rio, an open-space piñon-juniper area adjacent to Santa Fe which is owned by the US Forest Service. The proposed powerline would run from the Caja across the Rio Grande to the Lab. It would carry not just power but also a broadband fiber line, something Los Alamos town, if not the Lab, needs badly. On the other hand, those opposed worry about road-building and habitat destruction in the Caja.

I'm always puzzled reading accounts of the debate. There already is a powerline running through the Caja and across the Rio via Powerline Point. The discussions never say (a) whether the proposed line would take a different route, and if so, (b) Why? why can't they just tack on some more lines to the towers along the existing route?

For instance, in the slides from one of the public meetings, the map on slide 9 not only doesn't show the existing powerline, but also uses a basemap that has no borders and NO ROADS. Why would you use a map that doesn't show roads unless you're deliberately trying to confuse people?

I wanted a map of the current powerline, so I could compare it to the proposed new route. OpenStreetMap shows power lines if you zoom in far enough (Google Maps doesn't). But I have to zoom out to see where the powerline is actually going — and at that zoom level the powerline disappears. I needed a way to fetch the data for the powerline so I can examine the route myself.

Get a Feature ID

To make sure I was really looking at the powerline, I zoomed in and clicked the Query features button (an arrow with a question mark next to it) in the right panel at the bottom of the right-side toolbar. It goes green when you're in Query mode, as you can see on the screenshot. Then I clicked on what I thought was the powerline. Tip: the question-mark cursor has its hot spot at the bottom dot, so use the dot as your pointer.

In Nearby features it showed one Way (Way is the term OSM uses for any linear feature, like a road). I clicked on the link for Way 14598545 and verified that it was the powerline. OSM helpfully highlighted the Way object for me — which showed that this feature was only a small piece of the powerline. So I would need to keep clicking, and clicking, to find out where the powerline actually goes.

At this point I thought: there must be a better way, where I could fetch details about Way 14598545 and then follow all other Ways that connect to it and which are also powerlines.

I've been looking for an excuse to learn more about extracting feature details from OpenStreetMap, and this looked like a perfect excuse.

Getting Details From OSM's APIs

I poked around among a few options, and OSMPythonTools looked like the best option. It's not available as a Debian package, so you have to install it with pip.

There's a wiki quickstart page that shows the basics of the three APIs OSMPythonTools covers: the OSM API, which is mostly geared toward making changes to the map data; Nominatim, mostly used for searches by name; and Overpass, which is more general but uses its own arcane database language. Fortunately OSMPythonTools makes it Overpass little easier to deal with.

Get One Object

Getting a single way whose ID you already know is straightforward:

>>> from OSMPythonTools.api import Api
>>> api = Api()
>>> way = api.query('way/14598545')
[api] downloading data: way/14598545
>>> way
<OSMPythonTools.api.ApiResult object at 0x7fe150fd33e0>
>>> way.nodes()
[
    <OSMPythonTools.api.ApiResult object at 0x7fa576732420>,
    <OSMPythonTools.api.ApiResult object at 0x7fa576732600>,
...
    <OSMPythonTools.api.ApiResult object at 0x7fa57676e300>
]
>>> way.geometry()
[api] downloading data: way/14598545/full
{"coordinates": [[-106.188748, 35.708005], [-106.188221, 35.710108], [-106.187753, 35.711926], [-106.187296, 35.713701],
... ], "type": "LineString"}

The nodes are the points that define the powerline way; the geometry is the actual coordinates of the nodes.

I initially hoped that from those nodes, I could find which nodes belonged to more than one way, and thus trace the segments of the powerline. But after several days of searching, I was unable to find a way to go from a Node to a list of the Ways that include it. There's probably a way using the Overpass query language, but I found that impenetrable. So I took another, less elegant, approach.

Find All Powerlines Within a Bounding Box

A bounding box (bbox) contains four coordinates: south latitude, west longitude, north latitude, east longitude.

How do you define a bounding box and get the coordinates? The OpenStreetMap.org website apparently used to have a tool for that, but it's been removed. But there are other websites that offer that, like bbox tool, or bboxfinder (use the rectangle tool). or the XAPI Query Builder which is harder to use, but which also can generate Overpass queries. Or you can use my PyTopo program: right-click on the lower left point then choose the coordinates at the top of the context menu; do the same for the upper right point. PyTopo prints the coordinates to standard output, so run it from a terminal.

So, given a bounding box, let's find some powerlines inside it. Initially, I made a box from the southwest corner of Rio Rancho (but not including Albuquerque) to northeast of Santa Fe.

Reading the documentation on the Overpass API and how to use it in OSMPythonTools, I can build up a query for all the powerlines inside the bounding box:

from OSMPythonTools.overpass import Overpass
overpass = Overpass()
from OSMPythonTools.overpass import overpassQueryBuilder

query = overpassQueryBuilder(bbox=[35.236646, -106.710205, 36.151182, -105.759888],
                             elementType='way',
                             selector='"power"="line"',
                             out='body')
powerlines = overpass.query(query)

How many powerlines did it find?

>>> len(powerlines.ways())
235

It's interesting browsing the tags.

>>> for way in powerlines.ways():
...     print(way.tags())
{'power': 'line'}
{
    'operator': 'Public Service Company of New Mexico',
    'operator:short': 'PNM',
    'operator:wikidata': 'Q112271592',
    'power': 'line',
    'voltage': '115000'
}
{'name': '115kV transmission line', 'power': 'line', 'voltage': '115000'}
{
    'cables': '3',
    'frequency': '60',
    'operator': 'Public Service Company of New Mexico',
    'operator:short': 'PNM',
    'operator:wikidata': 'Q112271592',
    'power': 'line',
    'ref': 'CE',
    'voltage': '115000',
    'wires': 'single'
}
[ ... ]
{'power': 'line', 'tiger:cfcc': 'C20', 'tiger:county': 'Los Alamos, NM'}
{'layer': '1', 'power': 'line'}
{'power': 'line'}
{'power': 'line', 'tiger:cfcc': 'C20', 'tiger:county': 'Los Alamos, NM'}
{'power': 'line'}
{'name': 'Buckman - Pajarito Line', 'power': 'line', 'voltage': '115000'}
{'power': 'line'}
{'power': 'line'}
[ ... ]

Most power lines just have {'power': 'line'}, but some have lots of other useful information, sometimes including a "name" tag.

That 'Buckman - Pajarito Line' is part of the line I'm after, so let's find that to use as a starting point.

for way in powerlines.ways():
    try:
        if 'Buckman - Pajarito' in way.tag('name'):
            buckman = way
            break
    except:
        continue

It found this line:

>>> buckman
<OSMPythonTools.element.Element object at 0x7fa57671faa0>
>>> buckman.tags()
{'name': 'Buckman - Pajarito Line', 'power': 'line', 'voltage': '115000'}
>>> buckman.id()
370732405

To find ways that connect to the current way, search through ways comparing IDs of the nodes they contain. For this powerline case, a simplification is that powerlines only seem to connect at the beginning and end, so you don't need to test all the nodes, only the first and last.

Once I have a list of ways, I can call .geometry() on each way to get a list of [longitude, latitude] pairs, then I can dump that to a GPX file that I can view in a mapping app. But there's one trick: I had to call

from OSMPythonTools.api import Api

even though I don't use the Api interface myself, or calling .geometry() will die with Exception: module 'OSMPythonTools' has no attribute 'api'. I guess the import also sets up some variables that the geometry call depends on.

If you want the details, you'll find them in trace_ways.py.

I was able to trace the powerline down to Bernalillo. But as I was writing this, a new wrinkle appeared in the LANL powerline story, which required more analysis and a bigger bounding box. I'm now more confused than ever about the proposed powerline project. But I'll save that for another article.

Tags: mapping, GIS, programming, python, Los Alamos
[ 12:14 Aug 01, 2024    More mapping | permalink to this entry | ]