Shallow Thoughts : : mapping

Akkana's Musings on Open Source Computing and Technology, Science, and Nature.

Tue, 01 Oct 2019

Making Web Maps using Python, Folium and Shapefiles

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:

ogr2ogr -t_srs EPSG:4326 -f GeoJSON file.json file.shp

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

Each of these requires some extra work.

You can color the regions with a style function:

folium.GeoJson(jsonfile, style_function=style_fcn).add_to(m)

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:

def highlight_fcn(x):
    return { 'fillColor': '#ff0000' }

For tooltips:

tooltip = folium.GeoJsonTooltip(fields=['NAME'])
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:
$ jq . House.json | less
{
  "type": "FeatureCollection",
  "name": "cb_2018_35_sldl_500k",
  "crs": {
    "type": "name",
    "properties": {
      "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
    }
  },
  "features": [
    {
      "type": "Feature",
      "properties": {
        "STATEFP": "35",
        "SLDLST": "009",
        "AFFGEOID": "620L600US35009",
        "GEOID": "35009",
        "NAME": "9",
        "LSAD": "LL",
        "LSY": "2018",
        "ALAND": 3405159792,
        "AWATER": 5020507
      },
      "geometry": {
        "type": "Polygon",
        "coordinates": [
...

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:

Tags: , ,
[ 12:29 Oct 01, 2019    More mapping | permalink to this entry | comments ]

Tue, 23 Jul 2019

360 Panoramas with Povray and/or ImageMagick

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).

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:

camera {
    cylinder 1

    location <0.344444, 0.029620, 0.519048>
    look_at  <0.344444, 0.029620, 0>

    angle 360
}
with the same light_source and height_field as in Part III.

[360 povray panorama from Overlook Point]

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:

convert -size 3600x600 xc:black \
    outfile000.png -geometry +0+0 -composite \
    outfile045.png -geometry +400+0 -composite \
    outfile090.png -geometry +800+0 -composite \
    outfile135.png -geometry +1200+0 -composite \
    outfile180.png -geometry +1600+0 -composite \
    outfile225.png -geometry +2000+0 -composite \
    outfile270.png -geometry +2400+0 -composite \
    outfile315.png -geometry +2800+0 -composite \
    out-composite.png
or simply
convert outfile*.png +smush -400 out-smush.png

Adjusting Panoramas in GIMP

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

Tags: , , ,
[ 15:28 Jul 23, 2019    More mapping | permalink to this entry | comments ]

Wed, 17 Jul 2019

Ray-Tracing Digital Elevation Data in 3D with Povray (Part III)

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

[Hillshade of northern New Mexico mountains] 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.

wget -O mountains.tif 'http://opentopo.sdsc.edu/otr/getdem?demtype=SRTMGL3&west=-106.8&south=35.1&east=-105.0&north=36.5&outputFormat=GTiff'

Create a hillshade to make sure it looks like the right region:

gdaldem hillshade mountains.tif hillshade.png
pho hillshade.png
(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?

>>> np.where(demarray == demarray.max())
(array([645]), array([1386]))
>>> np.where(demarray == demarray.min())
(array([1667]), array([175]))

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.

>>> demarray[808, 744], demarray[808, 744] * 3.28
(1878, 6159.839999999999)

1878 meters, 6160 feet. That's fine for Overlook Point. We have everything we need to set up a povray file.

Convert to PNG

As mentioned in Part II, povray will only accept height maps as a PNG file, so use gdal_translate to convert:

gdal_translate -ot UInt16 -of PNG mountains.tif mountains.png

Use the Data to Set Camera and Light Angles

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):

    location <.3444, 0.033, .481>
    look_at <.5, 0.033, 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:

light_source { <2, 1, -1> color <1,1,1> }

All Together Now

Put it all together in a mountains.pov file:

camera {
    location <.3444, 0.0330, .481>
    look_at <.5, 0.0287, 1>
}

light_source { <2, 1, -1> color <1,1,1> }

height_field {
    png "mountains.png"
    smooth
    pigment {
        gradient y
        color_map {
            [ 0 color <.7 .7 .7> ]
            [ 1 color <1 1 1> ]
        }
    }
    scale <1, 1, 1>
}
[Povray-rendering of Black and Otowi Mesas from Overlook Point] Finally, you can run povray and generate an image!
povray +A +W800 +H600 +INAME_OF_POV_FILE +OOUTPUT_PNG_FILE

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

Tags: , , , ,
[ 16:43 Jul 17, 2019    More mapping | permalink to this entry | comments ]

Fri, 12 Jul 2019

Height Fields in Povray (Ray Tracing Elevation Data, Part II)

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:

gdal_translate -ot UInt16 -of PNG demfile.tif demfile.png

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:

camera {
    location <.5, .5, 2>
    look_at  <.5, .6, 0>
}

light_source { <0, 2, 1> color <1,1,1> }

height_field {
    png "YOUR_DEM_FILE.png"

    smooth
    pigment {
        gradient y
        color_map {
            [ 0 color <.5 .5 .5> ]
            [ 1 color <1 1 1> ]
        }
    }

    scale <1, 1, 1>
}

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:

povray +A +W800 +H600 +Idemfile.pov +Orendered.png
then take a look at rendered.png in your favorite image viewer.

Simple Sample Data

['bowling pin' sample DEM for testing povray] 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>

[povray ray-traced bowling pin result]

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

Tags: , ,
[ 18:02 Jul 12, 2019    More mapping | permalink to this entry | comments ]

Sun, 07 Jul 2019

Working with Digital Elevation Models with GDAL and Python (Ray Tracing Elevation Data, Part I)

Part III of a four-part article:

One of my hiking buddies uses a phone app called Peak Finder. It's a neat program that lets you spin around and identify the mountain peaks you see.

Alas, I can't use it, because it won't work without a compass, and [expletive deleted] Samsung disables the compass in their phones, even though the hardware is there. I've often wondered if I could write a program that would do something similar. I could use the images in planetarium shows, and could even include additions like predicting exactly when and where the moon would rise on a given date.

Before plotting any mountains, first you need some elevation data, called a Digital Elevation Model or DEM.

Get the DEM data

Digital Elevation Models are available from a variety of sources in a variety of formats. But the downloaders and formats aren't as well documented as they could be, so it can be a confusing mess.

USGS

[Typical experience with USGS map tiles not loading] USGS steers you to the somewhat flaky and confusing National Map Download Client. Under Data in the left sidebar, click on Elevation Products (3DEP), select the accuracy you need, then zoom and pan the map until it shows what you need.

Current Extent doesn't seem to work consistently, so use Box/Point and sweep out a rectangle. Then click on Find products. Each "product" should have a download link next to it, or if not, you can put it in your cart and View Cart.

Except that National Map tiles often don't load, so you can end up with a mostly-empty map (as shown here) where you have no idea what area you're choosing. Once this starts happening, switching to a different set of tiles probably won't help; all you can do is wait a few hours and hope it gets better..

Or get your DEM data somewhere else. Even if you stick with the USGS, they have a different set of DEM data, called SRTM (it comes from the Shuttle Radar Topography Mission) which is downloaded from a completely different place, SRTM DEM data, Earth Explorer. It's marginally easier to use than the National Map and less flaky about tile loading, and it gives you GeoTIFF files instead of zip files containing various ArcGIS formats. Sounds good so far; but once you've wasted time defining the area you want, suddenly it reveals that you can't download anything unless you first make an account, and you have to go through a long registration process that demands name, address and phone number (!) before you can actually download anything.

Of course neither of these sources lets you just download data for a given set of coordinates; you have to go through the interactive website any time you want anything. So even if you don't mind giving the USGS your address and phone number, if you want something you can run from a program, you need to go elsewhere.

Unencumbered DEM Sources

Fortunately there are several other sources for elevation data. Be sure to read through the comments, which list better sources than in the main article.

The best I found is OpenTypography's SRTM API, which lets you download arbitrary areas specified by latitude/longitude bounding boxes.

Verify the Data: gdaldem

[Making a DEM visible with GIMP Levels] Okay, you've got some DEM data. Did you get the area you meant to get? Is there any data there? DEM data often comes packaged as an image, primarily GeoTIFF. You might think you could simply view that in an image viewer -- after all, those nice preview images they show you on those interactive downloaders show the terrain nicely. But the actual DEM data is scaled so that even high mountains don't show up; you probably won't be able to see anything but blackness.

One way of viewing a DEM file as an image is to load it into GIMP. Bring up Colors->Levels, go to the input slider (the upper of the two sliders) and slide the rightmost triangle leftward until it's near the right edge of the histogram. Don't save it that way (that will mess up the absolute elevations in the file); it's just a quick way of viewing the data.

[hillshade generated by gdaldem] A better way to check DEM data files is a beautiful little program called gdaldem. It has several options, like generating a hillshade image:

gdaldem hillshade n35_w107_1arc_v3.tif hillshade.png

Then view hillshade.png in your favorite image viewer and see if it looks like you expect. Having read quite a few elaborate tutorials on hillshade generation over the years, I was blown away at how easy it is with gdaldem.

Here are some other operations you can do on DEM data.

Translate the Data to Another Format

gdal has lots more useful stuff beyond gdaldem. For instance, my ultimate goal, ray tracing, will need a PNG:

gdal_translate -ot UInt16 -of PNG srtm_54_07.tif srtm_54_07.png

gdal_translate can recognize most DEM formats. If you have a complicated multi-file format like ARCGIS, try using the name of the directory where the files live.

Get Vertical Limits, for Scaling

What's the highest point in your data, and at what coordinates does that peak occur? You can find the highest and lowest points easily with Python's gdal package if you convert the gdal.Dataset into a numpy array:

import gdal
import numpy as np

demdata = gdal.Open(filename)
demarray = np.array(demdata.GetRasterBand(1).ReadAsArray())
print(demarray.min(), demarray.max())

That gives you the highest and lowest elevations. But where are they in the data? That's not super intuitive in numpy; the best way I've found is:

indices = np.where(demarray == demarray.max())
ymax, xmax = indices[0][0], indices[1][0]
print("The highest point is", demarray[ymax][xmax])
print("  at pixel location", xmax, ymax)

Translate Between Lat/Lon and Pixel Coordinates

But now that you have the pixel coordinates of the high point, how do you map that back to latitude and longitude? That's trickier, but here's one way, using the affine package:

import affine

affine_transform = affine.Affine.from_gdal(*demdata.GetGeoTransform())
lon, lat = affine_transform * (xmax, ymax)

What about the other way? You have latitude and longitude and you want to know what pixel location that corresponds to? Define an inverse transformation:

inverse_transform = ~affine_transform
px, py = [ round(f) for f in inverse_transform * (lon, lat) ]

Those transforms will become important once we get to Part III. But first, Part II, Understand Povray: Height Fields in Povray

Tags: , ,
[ 18:15 Jul 07, 2019    More mapping | permalink to this entry | comments ]

Mon, 15 Apr 2019

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.

[QGIS categorized layers] 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

[Land ownership overlay in 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: , ,
[ 18:13 Apr 15, 2019    More mapping | permalink to this entry | comments ]

Wed, 10 Apr 2019

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:


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

[Georeferenced map overlay 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: , , ,
[ 19:08 Apr 10, 2019    More mapping | permalink to this entry | comments ]

Fri, 26 Aug 2016

More map file conversions: ESRI Shapefiles and GeoJSON

I recently wrote about Translating track files between mapping formats like GPX, KML, KMZ and UTM But there's one common mapping format that keeps coming up that's hard to handle using free software, and tricky to translate to other formats: ESRI shapefiles.

ArcGIS shapefiles are crazy. Typically they come as an archive that includes many different files, with the same base name but different extensions: filename.sbn, filename.shx, filename.cpg, filename.sbx, filename.dbf, filename.shp, filename.prj, and so forth. Which of these are important and which aren't?

To be honest, I don't know. I found this description in my searches: "A shape file map consists of the geometry (.shp), the spatial index (.shx), the attribute table (.dbf) and the projection metadata file (.prj)." Poking around, I found that most of the interesting metadata (trail name, description, type, access restrictions and so on) was in the .dbf file.

You can convert the whole mess into other formats using the ogr2ogr program. On Debian it's part of the gdal-bin package. Pass it the .shp filename, and it will look in the same directory for files with the same basename and other shapefile-related extensions. For instance, to convert to KML:

 ogr2ogr -f KML output.kml input.shp

Unfortunately, most of the metadata -- comments on trail conditions and access restrictions that were in the .dbf file -- didn't make it into the KML.

GPX was even worse. ogr2ogr knows how to convert directly to GPX, but that printed a lot of errors like "Field of name 'foo' is not supported in GPX schema. Use GPX_USE_EXTENSIONS creation option to allow use of the <extensions> element." So I tried ogr2ogr -f "GPX" -dsco GPX_USE_EXTENSIONS=YES output.gpx input.shp but that just led to more errors. It did produce a GPX file, but it had almost no useful data in it, far less than the KML did. I got a better GPX file by using ogr2ogr to convert to KML, then using gpsbabel to convert that KML to GPX.

Use GeoJSON instead to preserve the metadata

But there is a better way: GeoJSON.

ogr2ogr -f "GeoJSON" -t_srs crs:84 output.geojson input.shp

That preserved most, maybe all, of the metadata the .dbf file and gave me a nicely formatted file. The only problem was that I didn't have any programs that could read GeoJSON ...

[PyTopo showing metadata from GeoJSON converted from a shapefile]

But JSON is a nice straightforward format, easy to read and easy to parse, and it took surprisingly little work to add GeoJSON parsing to PyTopo. Now, at least, I have a way to view the maps converted from shapefiles, click on a trail and see the metadata from the original shapefile.

See also:

Tags: , , ,
[ 12:11 Aug 26, 2016    More mapping | permalink to this entry | comments ]