New Mexicans are used to people thinking we're not part of the US.
Every New Mexican has stories, like trying to mail-order
something and being told "We don't ship outside the US".
I had a little spare time and decided I'd follow a tutorial that's been
on my to-do list for a while:
Creating Beautiful River Maps with Python.
It combines river watercourse data from
gaia.geosci.unc.edu
with watershed boundaries from the
HydroSheds
project using Python and GeoPandas, making a map that is,
as promised in the title, beautiful.
I've been using the Wildland Fires map from MappingSupport.com
to keep an eye on the Cerro Pelado fire and the larger (though more
distant from me) Hermit's Peak/Calf Canyon fires raging in the Pecos.
It's an excellent map, but it's a little sporadic in whether it shows
the fire perimeter. In any case, as a data junkie, I wanted to know how
to get the data and make my own display, maybe for a quick viewer that
I can pop up when I sign on in the morning.
Also, Los Alamos County, on its
Cerro
Pelado Information page, has a map showing the "Go" lines (if the
fire crosses these lines, we have to evacuate) for Los Alamos and
White Rock and I'd like to be able to view those lines
on the same map with the fire perimeter and hot spots.
I record track files in OsmAnd
for most of my hikes. But about a quarter of the
time, when I get back to the car, I forget to turn tracking off.
So I end up with a track file that shows the hike plus at least
part of the car trip afterward. Not very useful for purposes like
calculating mileage.
My PyTopo can do simple
track edits, like splitting the track into two or deleting from some
point to the end. But most of the time it's easier just to edit the
GPX file with a text editor.
Eesh, edit GPX directly?? Sounds like something you oughtn't
do, doesn't it? But actually, GPX (a form of
XML) is human readable
and editable. And this specific case, separating the walking from the
car trip in an OsmAnd track file, is particularly easy, because OsmAnd
helpfully adds a speed to every point it saves.
These instructions seem long, but really, once you've done it once,
you'll realize that it's all very straightforward; explaining the steps
is harder than doing them.
San Juan County Council districtsSan Juan County Council voting precincts
For this year's LWV NM Voter Guides at
VOTE411.org,
I've been doing a lot of GIS fiddling, since the system needs to know the
voting districts for each race.
You would think it would be easy to find GIS for voting districts —
surely that's public information? — but counties and the state
are remarkably resistant to giving out any sort of data
(they're happy to give you a PDF or a JPG),
so finding the district data takes a lot of searching.
Often, when we finally manage to get GIS info, it isn't for what we want.
For instance, for San Juan County, there's a file that claims to be
County Commission districts (which would look like the image above left),
but the shapes in the file are actually voting precincts (above right).
A district is made up of multiple precincts; in San Juan, there are 77
precincts making up five districts.
In a case like that, you need some way of combining several shapes (a
bunch of precincts) into one (a district).
GIS "Dissolving"
It turns out that the process of coalescing lots of small
shapes into a smaller number of larger shapes is unintuitively
called "dissolving".
A friend recently introduced me to Folium, a quick and easy way of
making web maps with Python.
The Folium
Quickstart gets you started in a hurry. In just two lines of Python
(plus the import line), you can write an HTML file that you can load
in any browser to display a slippy map, or you can display it inline in a
Jupyter notebook.
Folium uses the very mature Leaflet
JavaScript library under the hood. But it lets you do all the
development in a few lines of Python rather than a lot of lines
of Javascript.
Having run through most of the quickstart, I was excited to try
Folium for showing GeoJSON polygons. I'm helping with a redistricting
advocacy project; I've gotten shapefiles for the voting districts in
New Mexico, and have been wanting to build a map that shows them
which I can then extend for other purposes.
Step 1: Get Some GeoJSON
The easiest place to get voting district data is from TIGER, the
geographic arm of the US Census.
For the districts resulting from the 2010 Decadal Census,
start here:
Cartographic
Boundary Files - Shapefile (you can also get them as KML,
but not as GeoJSON). There's a category called
"Congressional Districts: 116th Congress", and farther down the page,
under "State-based Files", you can get shapefiles for the upper and
lower houses of your state.
You can also likely download them from
at www2.census.gov/geo/tiger/TIGER2010/,
as long as you can figure out how to decode the obscure directory names.
ELSD and POINTLM, so
the first step is to figure out what those mean; I never found anything
that could decode them.
(Before I found the TIGER district data, I took a more roundabout path that
involved learning how to merge shapes; more on that in a separate post.)
Okay, now you have a shapefile (unzip the TIGER file to get a bunch of
files with names like cb_2018_35_sldl_500k.* -- shape "files"
are an absurd ESRI concept that actually use seven separate files for
each dataset, so they're always packaged as a zip archive and programs
that read shapefiles expect that when you pass them a .shp,
there will be a bunch of other files with the same basename but
different extensions in the same directory).
But Folium can't handle shapefiles, only GeoJSON. You can do that
translation with a GDAL command:
Or you can do it programmatically with the GDAL Python bindings:
def shapefile2geojson(infile, outfile, fieldname):
'''Translate a shapefile to GEOJSON.'''
options = gdal.VectorTranslateOptions(format="GeoJSON",
dstSRS="EPSG:4326")
gdal.VectorTranslate(outfile, infile, options=options)
The EPSG:4326 specifier, if you read man ogr2ogr, is supposedly
for reprojecting the data into WGS84 coordinates, which is what most
web maps want (EPSG:4326 is an alias for WGS84). But it has an equally
important function: even if your input shapefile is already in WGS84,
adding that option somehow ensures that GDAL will use degrees as the
output unit. The TIGER data already uses degrees so you don't strictly
need that, but some data, like the precinct data I got from UNM RGIS,
uses other units, like meters, which will confuse Folium and Leaflet.
And the TIGER data isn't in WGS84 anyway; it's in GRS1980 (you can tell
by reading the .prj file in the same directory as the .shp).
Don't ask me about details of all these different geodetic reference systems;
I'm still trying to figure it all out. Anyway, I recommend adding the
EPSG:4326 as the safest option.
Step 2: Show the GeoJSON in a Folium Map
In theory, looking at the Folium Quickstart, all you need is
folium.GeoJson(filename, name='geojson').add_to(m).
In practice, you'll probably want to more, like
color different regions differently
show some sort of highlight when the user chooses a region
Here's a simple style function that chooses random colors:
import random
def random_html_color():
r = random.randint(0,256)
g = random.randint(0,256)
b = random.randint(0,256)
return '#%02x%02x%02x' % (r, g, b)
def style_fcn(x):
return { 'fillColor': random_html_color() }
I wanted to let the user choose regions by clicking, but it turns out
Folium doesn't have much support for that (it may be coming in a
future release). You can do it by reading the GeoJSON yourself,
splitting it into separate polygons and making them all separate Folium
Polygons or GeoJSON objects, each with its own click behavior; but
if you don't mind highlights and popups on mouseover instead of
requiring a click, that's pretty easy. For highlighting in red whenever
the user mouses over a polygon, set this highlight_function:
In this case, 'NAME' is the field in the shapefile that I want to display
when the user mouses over the region.
If you're not sure of the field name, the nice thing about GeoJSON
is that it's human readable. Generally you'll want to look inside
"features", for "properties" to find the fields defined for each polygon.
For instance, if I use jq to prettyprint the JSON generated for the NM
state house districts:
If you still aren't sure which property name means what (for example,
"NAME" could be anything), just keep browsing through the JSON file to
see which fields change from feature to feature and give the values
you're looking for, and it should become obvious pretty quickly.
Here's a working code example:
polidistmap.py,
and here's an example of a working map:
This is Part IV of a four-part article on ray tracing digital elevation
model (DEM) data.
The goal: render a ray-traced image of mountains from a digital
elevation model (DEM).
Part IV, Make a 3D Panorama:
Making a Panorama from Raytraced DEM Images (this article)
Except there are actually several more parts on the way, related to
using GRASS to make viewsheds. So maybe this is actually a five- or
six-parter. We'll see.
The Easy Solution
Skipping to the chase here ... I had a whole long article
written about how to make a sequence of images with povray, each
pointing in a different direction, and then stitch them together
with ImageMagick.
But a few days after I'd gotten it all working, I realized none of it
was needed for this project, because ... ta-DA —
povray accepts this argument inside its camera section:
angle 360
Duh! That makes it so easy.
You do need to change povray's projection to cylindrical;
the default is "perspective" which warps the images.
If you set your look_at to point due south --
the first and second coordinates are the same as your observer coordinate,
the third being zero so it's looking southward -- then povray will
create a lovely strip starting at 0 degrees bearing (due north), and
with south right in the middle.
The camera section I ended up with was:
with the same light_source and height_field as in Part III.
There are still some more steps I'd like to do.
For instance, fitting names of peaks to that 360-degree pan.
The rest of this article discusses some of the techniques I would
have used, which might be useful in other circumstances.
A Script to Spin the Observer Around
Angles on a globe aren't as easy as just adding 45 degrees to the
bearing angle each time. You need some spherical trigonometry to
make the angles even, and it depends on the observer's coordinates.
Obviously, this wasn't something I wanted to calculate by hand, so
I wrote a script for it:
demproj.py.
Run it with the name of a DEM file and the observer's coordinates:
demproj.py demfile.png 35.827 -106.1803
It takes care of calculating the observer's elevation, normalizing
to the image size and all that. It generates eight files, named
outfileN.png, outfileNE.png etc.
Stitching Panoramas with ImageMagick
To stitch those demproj images manually in ImageMagick, this should
work in theory:
But in practice, some of the images have a few-pixel offset,
and I never did figure out why; maybe it's a rounding error
in my angle calculations.
I opened the images as layers in GIMP, and used my GIMP script
Pandora/
to lay them out as a panorama. The cylindrical projection should
make the edges match perfectly, so you can turn off the layer masking.
Then use the Move tool to adjust for the slight errors
(tip: when the Move tool is active, the arrow keys will move
the current layer by a single pixel).
If you get the offsets perfect and want to know what they are
so you can use them in ImageMagick or another program, use
GIMP's Filters->Python-Fu->Console.
This assumes the panorama image is the only one loaded in GIMP,
otherwise you'll have to inspect gimp.image_list() to see where
in the list your image is.
>>> img = gimp.image_list()[0]
>>> for layer in img.layers:
... print layer.name, layer.offsets
This is Part III of a four-part article on ray tracing digital elevation
model (DEM) data.
The goal: render a ray-traced image of mountains from a digital
elevation model (DEM).
In Part II, I showed how the povray camera position and angle
need to be adjusted based on the data, and the position of the light
source depends on the camera position.
In particular, if the camera is too high, you won't see anything
because all the relief will be tiny invisible bumps down below.
If it's too low, it might be below the surface and then you
can't see anything.
If the light source is too high, you'll have no shadows, just a
uniform grey surface.
That's easy enough to calculate for a simple test image like the one I
used in Part II, where you know exactly what's in the file.
But what about real DEM data where the elevations can vary?
Explore Your Test Data
For a test, I downloaded some data that includes the peaks I can see
from White Rock in the local Jemez and Sangre de Cristo mountains.
(or whatever your favorite image view is, if not
pho).
The image at right shows the hillshade for the data I'm using, with a
yellow cross added at the location I'm going to use for the observer.
Sanity check: do the lowest and highest elevations look right?
Let's look in both meters and feet, using the tricks from Part I.
>>> import gdal
>>> import numpy as np
>>> demdata = gdal.Open('mountains.tif')
>>> demarray = np.array(demdata.GetRasterBand(1).ReadAsArray())
>>> demarray.min(), demarray.max()
(1501, 3974)
>>> print([ x * 3.2808399 for x in (demarray.min(), demarray.max())])
[4924.5406899, 13038.057762600001]
That looks reasonable. Where are those highest and lowest points,
in pixel coordinates?
Those coordinates are reversed because of the way numpy arrays
are organized: (1386, 645) in the image looks like
Truchas Peak (the highest peak in this part of the Sangres), while
(175, 1667) is where the Rio Grande disappears downstream off the
bottom left edge of the map -- not an unreasonable place to expect to
find a low point. If you're having trouble eyeballing the coordinates,
load the hillshade into GIMP and watch the coordinates reported at the
bottom of the screen as you move the mouse.
While you're here, check the image width and height. You'll need it later.
>>> demarray.shape
(1680, 2160)
Again, those are backward: they're the image height, width.
Choose an Observing Spot
Let's pick a viewing spot: Overlook Point in White Rock
(marked with the yellow cross on the image above).
Its coordinates are -106.1803, 35.827. What are the pixel coordinates?
Using the formula from the end of Part I:
>>> import affine
>>> affine_transform = affine.Affine.from_gdal(*demdata.GetGeoTransform())
>>> inverse_transform = ~affine_transform
>>> [ round(f) for f in inverse_transform * (-106.1803, 35.827) ]
[744, 808]
Just to double-check, what's the elevation at that point in the image?
Note again that the numpy array needs the coordinates in reverse order:
Y first, then X.
The camera should be at the observer's position, and povray needs
that as a line like
location <rightward, upward, forward>
where those numbers are fractions of 1.
The image size in pixels is
2160x1680, and the observer is at pixel location (744, 808).
So the first and third coordinates of location should
be 744/2160 and 808/1680, right?
Well, almost. That Y coordinate of 808 is measured from the top,
while povray measures from the bottom. So the third coordinate is
actually 1. - 808/1680.
Now we need height, but how do you normalize that? That's another thing
nobody seems to document anywhere I can find; but since we're
using a 16-bit PNG, I'll guess the maximum is 216 or 65536.
That's meters, so DEM files can specify some darned high mountains!
So that's why that location <0, .25, 0> line I got
from the Mapping Hacks book didn't work: it put the camera at
.25 * 65536 or 16,384 meters elevation, waaaaay up high in the sky.
My observer at Overlook Point is at 1,878 meters elevation, which
corresponds to a povray height of 1878/65536. I'll use the same value
for the look_at height to look horizontally. So now we can
calculate all three location coordinates: 744/2160 = .3444,
1878/65536 = 0.0287, 1. - 808/1680 = 0.5190:
location <.3444, 0.0287, .481>
Povray Glitches
Except, not so fast: that doesn't work. Remember how I mentioned in
Part II that povray doesn't work if the camera location is at ground
level? You have to put the camera some unspecified minimum distance
above ground level before you see anything. I fiddled around a bit and
found that if I multiplied the ground level height by 1.15 it worked,
but 1.1 wasn't enough. I have no idea whether that will work in
general. All I can tell you is, if you're setting location to
be near ground level and the generated image looks super dark
regardless of where your light source is, try raising your location a
bit higher. I'll use 1878/65536 * 1.15 = 0.033.
For a first test, try setting look_at to some fixed place in
the image, like the center of the top (north) edge (right .5, forward 1):
That means you won't be looking exactly north, but that's okay, we're
just testing and will worry about that later.
The middle value, the elevation, is the same as the camera elevation
so the camera will be pointed horizontally. (look_at can be at
ground level or even lower, if you want to look down.)
Where should the light source be? I tried to be clever and put the
light source at some predictable place over the observer's right
shoulder, and most of the time it didn't work. I ended up just
fiddling with the numbers until povray produced visible terrain.
That's another one of those mysterious povray quirks. This light
source worked fairly well for my DEM data, but feel free to experiment:
And once I finally got to this point I could immediately see it was correct.
That's Black Mesa (Tunyo) out in the valley a little right of center,
and I can see White Rock
canyon in the foreground with Otowi Peak on the other side of the canyon.
(I strongly recommend, when you experiment with this, that you choose a
scene that's very distinctive and very familiar to you, otherwise you'll
never be sure if you got it right.)
Next Steps
Now I've accomplished my goal: taking a DEM map and ray-tracing it.
But I wanted even more. I wanted a 360-degree panorama of
all the mountains around my observing point.
Povray can't do that by itself, but
in Part IV, I'll show how to make a series of povray renderings
and stitch them together into a panorama.
Part IV,
Making a Panorama
from Raytraced DEM Images
This is Part II of a four-part article on ray tracing digital elevation
model (DEM) data. (Actually, it's looking like there may be five or
more parts in the end.)
The goal: render a ray-traced image of mountains from a digital
elevation model (DEM).
My goal for that DEM data was to use ray tracing to show the elevations
of mountain peaks as if you're inside the image looking out at those peaks.
I'd seen the open source ray tracer povray used for that purpose in
the book Mapping Hacks: Tips & Tools for Electronic Cartography:
Hack 20, "Make 3-D Raytraced Terrain Models", discusses how to use
it for DEM data.
Unfortunately, the book is a decade out of date now, and lots of
things have changed. When I tried following the instructions in Hack
20, no matter what DEM file I used as input I got the same distorted
grey rectangle. Figuring out what was wrong meant understanding how
povray works, which involved a lot of testing and poking since the
documentation isn't clear.
Convert to PNG
Before you can do anything, convert the DEM file to a 16-bit greyscale PNG,
the only format povray accepts for what it calls height fields:
If your data is in some format like ArcGIS that has multiple files,
rather than a single GeoTIFF file, try using the name of the directory
containing the files in place of a filename.
Set up the .pov file
Now create a .pov file, which will look something like this:
The trick is setting up the right values for the camera and light source.
Coordinates like the camera location and look_at,
are specified by three numbers that represent
<rightward, upward, forward>
as a fraction of the image size.
Imagine your DEM tilting forward to lie flat in
front of you: the bottom (southern) edge of your DEM image
corresponds to 0 forward, whereas the top (northern) edge is 1 forward.
0 in the first coordinate is the western edge, 1 is the eastern.
So, for instance, if you want to put the virtual camera at the
middle of the bottom (south) edge of your DEM and look straight
north and horizontally, neither up nor down, you'd want:
location <.5, HEIGHT, 0>
look_at <.5, HEIGHT, 1>
(I'll talk about HEIGHT in a minute.)
It's okay to go negative, or to use numbers bigger than zero; that
just means a coordinate that's outside the height map.
For instance, a camera location of
location <-1, HEIGHT, 2>
would be off the west and north edges of the area you're mapping.
look_at, as you might guess, is the point the camera is looking at.
Rather than specify an angle, you specify a point in three dimensions
which defines the camera's angle.
What about HEIGHT?
If you make it too high, you won't see anything
because the relief in your DEM will be too far below you and will disappear.
That's what happened with the code from the book: it specified
location <0, .25, 0>, which, in current DEM files,
means the camera is about 16,000 feet up in the sky, so high that the
mountains shrink to invisibility.
If you make the height too low, then everything disappears
because ... well, actually I don't know why. If it's 0, then you're
most likely underground and I understand why you can't see anything,
but you have to make it significantly higher than ground level, and
I'm not sure why. Seems to be a povray quirk.
Once you have a .pov file with the right camera and light source,
you can run povray like this:
then take a look at rendered.png in your favorite image viewer.
Simple Sample Data
There's not much documentation for any of this. There's
povray:
Placing the Camera, but it doesn't explain details like which
number controls which dimension or why it doesn't work if you're too
high or too low. To figure out how it worked, I made a silly little
test image in GIMP consisting of some circles with fuzzy edges.
Those correspond to very tall pillars with steep sides: in these
height maps, white means the highest point possible, black means
the lowest.
Then I tried lots of different values for location and look_at
until I understood what was going on.
For my bowling-pin image, it turned out looking northward (upward)
from the south (the bottom of the image) didn't work, because the
pillar at the point of the triangle blocked everything else.
It turned out to be more useful to put the camera beyond the top
(north side) of the image and look southward, back toward the image.
location <.5, HEIGHT, 2>
look_at <.5, HEIGHT, 0>
The position of the light_source is also important.
For instance, for my circles, the light source given in the original hack,
<0, 3000, 0>,
is so high that the pillars aren't visible at all,
because the light is shining only on their tops and not on their sides.
(That was also true for most DEM data I tried to view.)
I had to move the light source much lower, so it illuminated the sides
of the pillars and cast some shadows, and that was true for DEM data
as well.
The .pov file above, with the camera halfway up the field (.5)
and situated in the center of the north end of the field,
looking southward and just slightly up from horizontal (.6),
rendered like this. I can't explain the two artifacts in the middle.
The artifacts at the tops and bottoms of the pillars are presumably
rounding errors and don't worry me.
Finally, I felt like I was getting a handle on povray camera positioning.
The next step was to apply it to real Digital Elevation Maps files.
I'll cover that in Part III, Povray on real DEM data:
Ray-Tracing Digital Elevation Data in 3D with Povray
![[North and Central American rivers in New Mexico]](https://shallowsky.com/blog/images/na-rivers/na-ca-rivers.jpg)
New Mexicans are used to people thinking we're not part of the US.
![[Fire map from mapping support.com]](https://shallowsky.com/blog/images/cerropelado/mappingsupport.jpg)
![[San Juan County Council Districts]](http://shallowsky.com/blog/images/screenshots/qgis/SJCo_commdist.jpg)
![[San Juan County voting precincts]](http://shallowsky.com/blog/images/screenshots/qgis/SJCo_precincts.jpg)
![[360 povray panorama from Overlook Point]](http://shallowsky.com/blog/images/dem-povray/povray360-south.jpg)
![[Hillshade of northern New Mexico mountains]](http://shallowsky.com/blog/images/dem-povray/hillshade-overlook.jpg)
For a test, I downloaded some data that includes the peaks I can see
from White Rock in the local Jemez and Sangre de Cristo mountains.
![[Povray-rendering of Black and Otowi Mesas from Overlook Point]](http://shallowsky.com/blog/images/dem-povray/overlooktest.jpg)
![['bowling pin' sample DEM for testing povray]](http://shallowsky.com/blog/images/dem-povray/bowlingDEM.png)
![[povray ray-traced bowling pin result]](http://shallowsky.com/blog/images/dem-povray/bowlingpov.jpg)
![[RSS feed]](https://shallowsky.com/images/logos/rssicon-25x25.png){kind=small}
![[Twitter]](https://shallowsky.com/images/logos/twitter-bird-white-on-blue-25x25.png){kind=small}
![[envelope]](https://shallowsky.com/images/logos/envelope-25x25.png){kind=small}