Making shaded relief directly from DEMs projected in degrees?

In a presentation I gave at GeoIgnite 2021, I was explaining the process of using QGIS to make shaded relief from DEMs, more or less as detailed in this post. There wasn’t time to explain the process of re-projection, so I rather boldly asserted that you don’t need to re-project your DEMs before making shaded relief. Instead, I said, you can just use the magic scaling number 111120 to convert from degrees to metres.

Really? Is this true?

Well, yes, in many cases it works pretty well. In QGIS you pull in an unreprojected DEM (its native CRS using degrees as units) but set the Project CRS to some reasonable metric projection, such as UTM. Then, opening the GDAL hillshade tool from the Toolbox, you use the DEM as your Input layer, set scale to 111120, and enter the parameters of your choice for Z factor, Azimuth of the light and Altitude of the light.

And, look: shaded relief:

This seems so much simpler than having to first re-project your DEM to square cells measured in metres. But there are some caveats. This technique fools with the terrain a bit, and could potentially distort how it looks. If you understand how it distorts it, you’ll be able to evaluate if you want to use it in your map.

The magic number

The magic number 111120 comes from the GDAL documentation, specifically the documentation for the gdaldem tool:

The gdaldem generally assumes that x, y and z units are identical. If x (east-west) and y (north-south) units are identical, but z (elevation) units are different, the scale (-s) option can be used to set the ratio of vertical units to horizontal. For LatLong projections near the equator, where units of latitude and units of longitude are similar, elevation (z) units can be converted to be compatible by using scale=370400 (if elevation is in feet) or scale=111120 (if elevation is in meters). For locations not near the equator, it would be best to reproject your grid using gdalwarp before using gdaldem.

https://gdal.org/programs/gdaldem.html

Gdaldem is the tool that I recommend using for creating hillshades in QGIS, and there in the docs they say it: unless your locale is near the equator, it’s best to re-project your grid before using gdaldem.

The obvious problem with entering a scale factor of 111,120 metres per degree is that as you move away from the equator, lines of longitude become closer to one another. For example, in Canada at 54° north, one degree of latitude is close to 111,0000 metres, but a degree of longitude is only about 65,000 metres. In other words, a DEM cell is not square on the ground. So when you enter 111,120 as the scaling factor to convert degrees to metres, you will get a hillshade calculated on the basis of a square cell that has been created by stretching a rectangular cell in an east–west direction.

For example, with a 1″ DEM, each cell is 0.0002777777° (that’s 1 ÷ 3600) on a side. If we scale both axes at 111,120, the hillshading algorithm thinks it has a square cell 30 metres on a side. However the reality at 54° north is that on the ground the cell is about 30 m north–south and 18 metres east–west.

How bad is it?

You’d expect that because the landscape is stretched east–west, the slope would be under-calculated for cells that slope east–west. Basically, these cells would appear to be on a slope less steep than it really is.

And slope is important in the hillshading algorithm: it plays a role in determining how bright a hillshade pixel is. Pixels whose aspect and slope indicate that they are orthogonal to the sun direction are the brightest white, while those facing away from the sun direction are darkest black.

So if, for example, you had a 45° slope facing west, and the sun in the west at an elevation of 45 °, you’d expect pixels on the slope to have the highest possible hillshade value, 255. But if you calculated the hillshade from a degree-DEM using a scale of 111,120, that slope would be under-estimated; maybe it would appear to be 35°. As a result it would not get assigned a 255, but something less, maybe 240. There might be a nearby area, also facing west, whose true slope was 55°, but it now gets assigned the brightest value because it appears to be 45°.

In general then, for east or west-facing surfaces, we’d expect the hillshade built from the un-re-projected degree-DEM to distribute its brightest and darkest patches differently, particularly when the sun is due east or west.

That’s just what you see.

However, when the sun is not due east or west it’s hard to see this difference in slope calculation.

You can also experiment with other, lower scale factors, such as 90,000, 80,000 or 70,000. These, in effect, cause slopes to be less under-estimated on east–west aspects, but over-estimated on north–south aspects. In other words you are distributing the error around.

This suggests that in many cases, making a hillshade directly from a DEM projected in degrees is fine.

Other questions

A second distortion arising from creating a hillshade using a degree-DEM and a scale factor of 111,120 concerns the sun’s azimuth. If you are mapping in a local UTM projection, you know that UTM grid north is usually different from true north. This was visible in the figure above with the contour lines: the cell that was “square” in degrees turns out to be rectangular in a UTM projection, and also a bit tilted to the side.

With UTM projections, this difference is not more than 3°. At worst, then, you might ask for a sun azimuth of 315° and instead get one of 312°. I doubt whether you can notice this difference.

The last problem we might wonder about is whether shadows calculated by the hillshading algorithm are distorted. Are they foreshortened in an east–west direction when the hillshade is re-projected into a metric projection? The answer to this may be surprising: the GDAL hillshading algorithm does not calculate shadows. It calculates shadowed surfaces, but this is purely on the basis of the slope and aspect of the cell it is looking at (and its eight neighbours).

In other words, a piece of terrain at the bottom of a deep gorge will be rendered purely on the basis of its slope and aspect, as if the sun could shine on it right through the surrounding mountains.

Hence, no, we do not need to worry that shadows will be shortened.

So…

All of these things considered, I come to the surprising (for me) conclusion that yes, you actually can often just hillshade a DEM that’s still projected in degrees, using the scale factor 111,120. It seems heretical, but when I can’t tell the difference in the results, I have to admit that for practical purposes it’s a good solution.

QGIS 3 and Vector map tiles

If you’ve ever tried printing a map from the Open Street Map tile service (available in QGIS through QuickMapServices plugin) you’ve experienced the annoyance of having all your labels come out tiny. In other words, you see this on the screen…

But then you get this out of your 300 dpi print layout….

The reason those labels and icons come out so small with raster tiles is that your screen display is something like 96 pixels/inch, but the print layout has to pack them down to 300 pixels/inch. The print layout also winds up asking for tiles at a higher zoom level. Text labels and icons, like everything else, are “baked” into the raster tiles at a reasonable size for display at 96 dpi. At 300 dpi, everything is much smaller.

In 2018 I described, at the end of a post, a way to work around this, but now there’s a better option. If you work instead with vector tiles, the text labels are not baked into the tile: they are rendered by QGIS. This means that if the label is supposed to be 2.5 mm high, it’ll be 2.5 mm high no matter what dpi you export at out of the print layout.

QGIS 3.14 and above can handle vector tile servers as data sources. This means the server sends tiles of vector data, which are then styled at the user’s end by QGIS.

How they work

Vector map tiles (VMTs) were not developed for GIS use. Their main use is for people viewing maps through browsers. You commonly have a map that’s been embedded in a web page that includes the javascript libraries for mapbox rendering, such as here.

But vector tiles offer the cartographer a lot. As raster tiles, OpenStreetMap data was ideal for areas where you just could not get decent data—if, say, you were mapping the streets of a city in Turkmenistan. But now, as vector tiles, this data allows you to control colours, line weight, font, and, of course the visibility of elements. Imagine your own version of OpenStreetMap, with pink highways, no railroads and emphasized parks, and you begin to get the idea.

There are multiple formats for tiles, but here we’ll just look at the Mapbox vector tiles format, since that is what QGIS can support.

Because vector map tiles are designed to be viewed over the internet, the most basic principle of VMTs is that you typically to have a tile server, which is firing off vector tiles of data, and a GL style file that understands what data layers are delivered in this server’s tiles, and tells the browser how to render that data. Often these are at different URLs.

The URLs of vector tile servers look like this:

The z, y and x in curly brackets get replaced by the zoom level, row and column, and the data comes as a PBF (Protobuf) file or a MVT (Mapbox vector tile) file.

The GL style file, on the other hand, is a text file in JavaScript Object Notation (JSON). It might begin like this, with information on the sources for tiles, sprites and glyphs (VMT lingo for markers and fonts):

{
  "version": 8,
  "name": "Positron",
  "metadata": {...  },
  "sources": {
    "openmaptiles": {
      "type": "vector",
      "url": "https://api.maptiler.com/tiles/v3/tiles.json?key={key}"
    }
  },
  "sprite": "https://openmaptiles.github.io/positron-gl-style/sprite",
  "glyphs": "https://api.maptiler.com/fonts/{fontstack}/{range}.pbf?key={key}",

and then it will continue with a list of the layers in the tiles, and how they should be styled.

  "layers": [
    {
      "id": "background",
      "type": "background",
      "paint": {"background-color": "rgb(242,243,240)"}
    },
    {
      "id": "park",
      "type": "fill",
      "source": "openmaptiles",
      "source-layer": "park",
      "filter": ["==", "$type", "Polygon"],
      "layout": {"visibility": "visible"},
      "paint": {"fill-color": "rgb(230, 233, 229)"}
    },
    {
      "id": "water",
      "type": "fill",
      "source": "openmaptiles",
      "source-layer": "water",
      "filter": [
        "all",
        ["==", "$type", "Polygon"],
        ["!=", "brunnel", "tunnel"]
      ],
      "layout": {"visibility": "visible"},
      "paint": {"fill-antialias": true, "fill-color": "rgb(194, 200, 202)"}
    },     ...

The styling information here is written in something called Mapbox GL style, and here are some GL Style file URLs:

Vector tiles are relatively new, and almost no one gives them away for free, but we do have a few heroic companies offering this. Most require API-keys or access-tokens, but these are frequently free.

The first question you might have is what you get if you ask for vector tiles from the tile server, but don’t combine them with a GL Style file. VMTs without a style file aren’t useless—but what you get is data where all points share a single style, as do all lines and polygons.

Data from the Nextgen vector tile server without styling

A simple VMT service brought in using QGIS3’s native abilities

Our best candidate at the moment for free OpenStreetMap (OSM) data delivered through vector tiles is the ESRI OSM vector tile layer. (Thank you ESRI!). The tile server is here:

https://basemaps.arcgis.com/arcgis/rest/services/OpenStreetMap_v2/VectorTileServer/tile/z/y/x.pbf

and the style file is here:

https://www.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/styles/root.json?f=pjson

We can add this tile server to QGIS by

  1. opening the Browser panel
  2. Right-clicking the new Vector Tiles category and choosing New ARCGIS vector tile service connection
  3. Filling it out like this:

Notice that the Service URL ends with “VectorTileServer.” There’s no “tile/{z}/{y}/{x}.pbf.” QGIS will append that automatically before making a tile request.

When we add this layer to our map canvas, and zoom to Baku, Azerbaijan at 1:8,028 scale, we get this pleasing map that looks remarkably like OpenStreetMap:

However, we also get a pile of errors, under the banner “Vector tiles: Style could not be completely converted.” If you click on “Details,” some will be like this:

Referenced font Arial Unicode MS Regular is not available on system

You can repair these font errors by installing the missing font – or simply accepting the local font that has been substituted.

But there are also a pile of errors along the lines of

amenity area/bowling alley: Could not retrieve sprite 'amenity area/bowling alley'

I will explore how to solve these sprite problems below. First let’s talk about what you’ve got so far.

Styling the vector tile layers

If I open the styling panel on the right of the map I get something truly new.

Wow, that’s the longest list of rule-based styling I’ve ever seen!

Can I literally uncheck or change the style for every type of layer? Are these layers within a layer?

Yes and yes. Vector map tiles contain many layers. (The set of layers is called the schema of a tile service.) Out of each layer (e.g., “landcover”), features are selected using a filter, which is quite similar to a rule in rule-based styling. For example, in the image above there is a specific fill style for features in the landcover layer for which _symbol = 48,and the zoom is 9 or higher. It’s called landcover/park.

I can uncheck landcover/park and get a map without parks…

Or I can double-click on the sample fill swatch and change the parks to a pink/white gradient fill.

In fact I can apply any QGIS style to these polygons: shapeburst fill, line fill, marker fill… I don’t have full control—I can’t change the min and max zooms that define this style, or the filter rule—but within the predefined filters, the possibilities are mind-boggling.

Once I like the way in which I have tweaked the styling, I can save the style for this layer to a .qml file. (Open the layer’s Properties, find the Style button, and choose Save Style…) This will allow me to apply this custom styling to the same vector tile layer, but in another QGIS project.

My powers over the styling of the VMTs are clumsy—a GL Style file typically defines scores of styles, and if I want to change, say, all my water styles from blue to grey, I will have to do that to each one individually. But it’s still very powerful.

There are a few other things to point out. Note at the bottom of the styling panel you have a checkbox for Visible rules only. Visible rules only cuts out styles that are not displayed at the current zoom. (It does not remove styling for features not visible in the current map canvas.)

The current zoom is also displayed there. This is always shown as an integer, but with vector map tiles you are not restricted to integer zoom levels. There’s no “Zoom to Native Resolution (100%)” like there is when you right-click a raster tile layer. The zoom level shows as 14, but the map canvas is at 1:10,000, which is somewhere between 1:17,061 (the break point for vector tiles changing from zoom 13 to 14) and 1:8,530 (the break point for vector tiles changing from zoom 14 to 15). If you do the math, you can conclude we’re at something like zoom level 14.17. This might make you queasy! And, welcome to the world of vector tiles.

[Note that zoom levels for vector tiles are one more than the zoom levels for typical raster map tiles (e.g., Google Maps, OSM ). This is a consequence of MapBox vector tiles being 512 pixels on a side. As well, I have a suspicion that, in QGIS, the detail shown on ESRI OSM vector tiles at zoom level 14 is the same as shown on OSM raster tiles at zoom level 13. But I’m not quite sure why this is.]

Styles are implemented in order, from the top of the list down, and you can drag them to reorder.

Finally, it’s important to understand that the styling rules you see here in QGIS are not the GL Style document itself. They are a translation of that document. You cannot (at, least, as far as I know) translate back from QGIS styling to a GL Style file.

Labels

If you move to the labelling tab in the style panel, you see this:

So, you also have extensive control over the labelling. Double click on a label’s style swatch, and change the font, size, colour… all of the usual label parameters are available.

For example, if I would like street names to be blue (bold-italic) and subway stations to be labelled in red, I can just do it by altering three styles.

Finally, note that while you can’t open an attribute table on your vector map tiles layer, you can use the Information tool. You will get results for every single layer under the mouse click.

I can now print at 300 dpi from this, and I get just what I see on the screen!

Sprites

All that might have been enough for you, and now you’re off to try customizing OpenStreetMap yourself. However, if you’re interested in fixing those sprite errors we encountered above, continue reading.

When you add a vector tile service to QGIS you frequently get this:

When you check out the Details, many of the errors refer to “could not retrieve sprite.”

At the same time, in the QGIS’s Network Log window, you’ll see something like this:

2020-12-31T20:27:24     WARNING    Error transferring https://cdn.arcgis.com/sharing/rest/content/items/684018c92e9b424ca8e900da8dad8b23/resources/styles/../sprites/sprite@2x.json - server replied: Not Found
2020-12-31T20:27:24     WARNING    Error transferring https://cdn.arcgis.com/sharing/rest/content/items/684018c92e9b424ca8e900da8dad8b23/resources/styles/../sprites/sprite.json - server replied: Not Found

So what’s going on?

If I paste the URL for the GL Style file into my browser (I’m using Firefox), the result is that whole code for the styling gets delivered to me in an unformatted blob.

It’s a bit dense reading. But near the top I can pick out this piece:

"sprite": "https://cdn.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/styles/../sprites/sprite"

Sprites, in vector map tile parlance, are markers. QGIS needs to use this URL look up the sprite image (a PNG file with all the sprites packed into it) and its partner sprite JSON file (which explains which piece of the image a given sprite is). But, for QGIS, the “styles/..” part of this particular URL is confusing. We might understand it as meaning descend into the styles sub-directory and then pop back up out of it, but QGIS doesn’t know this.

I can make a better address for the sprites by saving the GL Style file locally, editing it so the URL of the sprites is simpler, and then telling QGIS to now use my local version of the GL Style file.

I save the GL Style file locally by selecting all the text out of my browser, and copying and pasting into my text editor. Then I change the sprite URL by removing the “styles/..” part, so it now reads

"sprite": "https://cdn.arcgis.com/sharing/rest/content/items/3e1a00aeae81496587988075fe529f71/resources/sprites/sprite"

Once I’ve saved the file as ESRI_OSM_Style.json, I return to QGIS, where I edit the Vector Tiles connection I made for this server, and specify this new file as the “Style URL.” I could post it on a server somewhere and give a http:// URL to it, but, with Ubuntu, I can also just use the file:// protocol:

Sprites make quite a difference at high zoom levels.

No sprites, zoom level 17
With sprites, zoom level 17

A few final points

There are three different types of support for vector tiles in QGIS. The original one was the Vector Tiles Reader plugin, which received the vector tiles, converted them to geoJSON data, and made a different layer for each kind of features in the tiles. There is the new native support, which I have used here: it receives the tiles and applies filters, or rules, to style the many different kinds of features in them. Finally, there is the MapTiler plugin, which has a great deal of potential but requires a free account with MapTiler/OpenMapTiles, and this has an unfortunate data cap of 100,000 tile requests/month. (It seems like a lot, but while learning to use and style vector map tiles, I went through my 100,000 tiles in eleven days.)

The QGIS implementation of vector tiles permits you to add “New generic connection” if you have URLs for the server and the GL Style file. Yet outside of MapTiler and ESRI I have found precious few servers that also provide GL Style files. Here’s a quick survey:

  • The city of Wien (Vienna)—the data is limited to Austria. [Server URL] [GL Style URL]
  • Ordnance Survey UK—the data is limited to the U. K. and requires a free API key. [Server URL to which you have to append your free API key] [GL Style URL to which you have to append your free API key] Because of the API-key, QGIS does not yet seem to be able to make the request for the sprites correctly.
  • MapTiler/OpenMapTiles (OSM data)—beware the 100,000 tile limit. (Literally every time you zoom, pan or resize the screen you can trigger a tile request.) [Server URL to which you have to append your free access key] [GL Style URL for the “Basic” style]
  • Nextzen (OSM data)—lovely data, but I haven’t found a GL Style file. [Server URL to which you append your API key.]
  • ESRI generously offers another vector tile server based on more data sources than just OpenStreetMap. [Server URL] It comes with many different styles which you can find here. To obtain the GL Style file for one of these, you typically Open the style in Map Viewer and then in the upper right is a button called View style. Copy the link out if it and use it, as we did above, in the Style URL field. (Generally speaking, you’ll have to also do the same hack with the sprites as I showed above.)
  • Mapbox, if you can decode their mysterious URLs (mapbox://mapbox.mapbox-streets-v11), serves up vector tiles, and the public access key is free (without data limits), but I have not been able to find a GL Style file for its schema that QGIS can fully translate. [Server URL for Mapbox-streets-v8] [GL Style URL for MapBox-streets-v11]
  • Thunderforest has nice vector tile servers, such as their Transport layer, and their Outdoors layer, that you can have access to with a free API key—but I have not been able to find GL Style files for the schemas of these servers. [Server URL to which you have to append your free API key]

I understand that converting data to vector tiles, and keeping it up to date, costs money, so I don’t begrudge these companies charging fees or imposing limits. On the other hand, I really appreciate those which don’t!

As time passes there will no doubt be more vector map tile servers, and the way schemas and styling work will change. Nonetheless working with vector map tiles is addictive, and it will be exciting to watch this evolve.

Shaded Relief using Skymodels, courtesy of Raster Chunk Processing

A couple of weeks ago I watched some excellent presentations from the How To Do Map Stuff day organized by Daniel Huffman. One that I particularly enjoyed was Jake Adams’s talk on building shaded relief out of hillshades. Toward the end of his talk he brought in something called a skymodel.

In this post, I’ll explain what Skymodels are, and how to get Jake’s Raster Chunk Processing software running in a linux environment so you can use the Skymodel feature to make your own versions of this unique kind of shaded relief.

Skymodels were introduced in a paper by Pat Kennelly and James Stewart in 2014. The basic idea behind the skymodel is that, in the real world, we see terrain illuminated by light coming from all parts of the sky in various amounts, not just from the sun. The usual hillshading algorithm, on the other hand, calculates what’s illuminated and what’s in shadow based entirely on light coming from a single point.

An overcast sky is one example of a skymodel: in the overcast sky, light comes mostly from above, but not from any particular direction.  In James Stewart’s SkyLum software, this is represented by the redder, “hotter” dots in the dome of the model, while we have yellow and then cooler blue dots down toward the horizon representing less light coming in from those directions.

skymodel_overcast

Contrast that with this model for what SkyLum calls a Type 6 sky, “partly cloudy, no gradation towards zenith, slight brightening towards the sun.”

Type 6 partly cloudy, no gradation towards zenith, slight brightening towards the sun

Or this, the Type 13, “CIE standard clear sky, polluted atmosphere.”

type 13 CIE cstandard lear sky, polluted atmosphere

Each of these will produce a different kind of shaded relief.

Baltoro_assortment
The Baltoro glacier, Pakistan. Left: overcast. Middle: type 6, “partly cloudy, no gradation towards zenith, slight brightening towards the sun.” Right: type 13, “CIE standard clear sky, polluted atmosphere.”

In Blender you can position multiple light sources, but SkyLum writes a CSV file (which I will call the illumination file) that defines 250 different light sources and assigns a relative weight to each. This is new level of complexity in lighting.

180,45,0.000962426
175,37,0.000929834
193,47,0.00118523
187,38,0.000949383
167,47,0.00143157
181,54,0.00168069
169,41,0.00141631
...

Each line in the illumination file above gives the azimuth, elevation, and weight for one illumination source. The weights for all 250 points will add up to 1.

So how does one use this to generate shaded relief?

Well, Jake Adams (thank you, Jake!) has written a clever piece of code called Raster Chunk Processing (RCP, hereafter), which he presents in a three-part blog post. RCP divides up large DEMs into smaller tiles (“chunks”), each of which is processed separately. All of the results are then merged back together for a final result. The point of RCP is that it allows you to work with DEMs that ordinarily would max out your RAM and cause processing to grind to a halt.

This is similar to the way I have divided large DEMS into tiles for processing by Blender, but Jake’s RCP code allows you to use this “divide-and-conquer” strategy to do a whole host of things. One of them is to build a skymodel hillshade, given a DEM and an illumination file from SkyLum.

The RCP code is written in python, which is platform independent, so although Jake gives us instructions for installing it under Windows, we can get it running under linux with just a few changes.  In this case I did this on a system running Ubuntu 18 Bionic, and since I was using this machine for mapping I already had QGIS and GDAL installed.

To begin, head over to Jake’s Git repository, and download the RCP software using the Clone or Download button. Once unzipped, this will produce a directory called rcp-master. Within it you will see, among other things, the main program (raster_chunk_processing.py), a subdirectory called SkyLum (which contains the SkyLum program) and a couple of other python files that will be important to us, like settings.py and methods.py.

RCP is written to be run under python 3, so from a Ubuntu point of view these are the packages you will need to install:

  • python3-numba
  • python3-astropy
  • python3-gdal
  • python3-skimage
  • python3-numpy

(Note that python3-scimage replaces the module Jake calls scikit-image.)

Another thing that RCP needs is a working copy of mdenoise, the software that implements Sun’s algorithm for denoising topographic data. Or rather, RCP needs to at least think it has a copy. So you have a choice: if you want to be able use RCP to denoise DEMs, you should compile yourself a copy of the mdenoise binary; there are instructions here, and it’s not too painful. If not, just use sudo to place an empty file called mdenoise in /bin.

Then use a text editor on the file settings.py to alter its single line about where the mdenoise executable is.

MDENOISE_PATH = '/bin/mdenoise'

The last piece you have to put in place is a way to run SkyLum. SkyLum only comes as a Windows binary, SkyLum.exe, so to run this you need to have wine installed. The good news is that SkyLum runs quite well under wine.

Right-click SkyLum.exe, and choose Open With… >Other Application. In the list of application choose, “winebrowser.”

winebrowser

SkyLum should open right up and show you a piece of terrain with a sky dome over it.

fresh skylum

The complete instructions for SkyLum are in the README file included with it, but I will give a summary here of what I find useful:

  • By default when you open SkyLum the sun (sun position is shown at the bottom) is at 45° elevation and azimuth 180° (due South). 0° is north, and azimuths increase clockwise, as is standard
  • Hit ? for a help screen
  • Move the sun with the arrow keys.
  • Rt-click and choose an illumination model. See Kennely and Stewart’s paper for more on these.
  • Hit p to have SkyLum.exe calculate points after you have positioned the sun and chosen a skymodel
  • Hit o to have SkyLum write a sky model file. It’s conventional to give it a .csv extension.
  • Use a text editor to delete the header lines. All you want left are the comma-separated lines with the azimuth, elevation and weight, as shown above.

Now, how to run RCP?

Let’s assume you have a DEM  called myDEM.tif, and an illumination file you made with SkyLum called illum.csv. You’ve already deleted the header lines from illum.csv. You want to divide your DEM into 1000x1000px chunks, with an overlap of 200 px. We’ll also assume you have 4 processors and you want to use 3 of them for this operation.

Drop both myDEM.tif  and illum.csv in the rcp_master directory. Open a terminal there.

The general form of the RCP command line is:

python3 raster_chunk_processing.py -m [method] {general options} {method-specific options} input_file output_file

Notice that the first element in the command line is python3, not just python.

For my skymodel in this case the command line will be…

 python3 raster_chunk_processing.py -m skymodel -s 1000 -o 200 -p 3 -l illum.csv --verbose myDEM.tif myDEM_skymodel.tif

What you’ll see in the terminal is a bunch of information about the job to be done, and then you’ll see RCP submitting the sub-jobs (the chunks) to be skymodelled. You can walk away, or watch, fascinated, as the chucks get worked on…

Preparing output file myDEM_skymodel.tif...
Output dimensions: 4046 rows by 4515 columns.
Output data type: <class 'numpy.float32'>
Output size: 69.7 MiB
Output NoData Value: -32767.0

Processing chunks...
Tile 0-0: 1 of 25 (4.000%) started at 0:00:00.701178 Indices: [0:1200, 0:1200] PID: 8777
Tile 0-1: 2 of 25 (8.000%) started at 0:00:00.701253 Indices: [0:1200, 800:2200] PID: 8778
Tile 0-2: 3 of 25 (12.000%) started at 0:00:00.701284 Indices: [0:1200, 1800:3200] PID: 8779

All of the switches and parameters are explained quite well in Jake’s post. My additional notes are:

  • the input file is generally a geotiff, but it can also be a TIF/TFW pair. It has to have a NoData value defined or you will get an error No NoData value set in input DEM.
  • RCP will exit with an error if the output file already exists
  • By default RCP applies a vertical exaggeration of 5X to the DEM, because this was what Kennely and Stewart did in their paper. However, you can change this if you would prefer a different vertical exaggeration. Open methods.py in your text editor and go to line 272. This line originally says  in_array *= 5. You can change that 5 to whatever number you would prefer. (No recompile necessary.)
  • the overlap (-o) parameter is based on how far you think shadows may stretch. The shadowing algorithm checks outward for 600 pixels to see if a given point is shadowed by anything. For this reason, there is no point in making overlap larger then 600.
  • the chunk size parameter (-s) is based on how much RAM each process requires while running. You can experiment with this and watch on the system monitor to see how close you are to spilling over into swap.
  • When shadows are being calculated in the skymodel,

The output file can be dragged into QGIS, where I often find I want to increase brightness.

adjusdting_brightness
Left: hillshade image fresh out of RCP. Right: with brightness +100

Another idea that seems to produce some nice results is to combine two or more skymodels with different transparencies and styling. You might also want to check out a gallery of the results of applying all of the different skymodels in SkyLum to the same piece of terrain.

Now that you have RCP running, of course you’ll need to try all of the different skymodels on your favourite DEM to see which you like best. But it’s also worth checking out the other kinds of processing RCP can do on DEMs. It is a very fast conventional hillshader (method hillshade). It runs a nice, quick gaussian blur (method blur_gauss), much faster than the SAGA Gaussian filter module in QGIS. And I haven’t even tried the Clahe and TPI processing yet!

 

 

 

How do we understand the size of Syria?

Syrian_comparison_inset_2004_CIA

This inset appears on a CIA map of Syria from 2004. We can assume it’s meant to give the map reader a sense of the size of Syria by comparing it to a region that he or she is familiar with.

I’m fascinated, first, by the assumption that the decision-maker reading the map is located in the mid-Atlantic states, (Washington, presumably). This kind of regionality runs throughout American politics, a geography of where important people live, and where they don’t, that one carries in the memory and consults without even knowing it. 

This is an attempt to make the map *personal* but I almost think it should be captioned, “This might be helpful to you if you happen to live in the northeast.”

The second thing I’m fascinated by is, assuming we have to stick to the northeast USA, would it have been a smarter decision  to have aligned Damascus with Washington, DC? These are the two national capitals. Such a juxtaposition would, in the context of the Syrian Civil War, allow the President to imagine New York no longer doing his bidding, and Providence functioning like an independent state.

CIA inset Washington and Damascus aligned

Perhaps this juxtaposition puts too much of Syria out to sea, but it does put New York City more or less where Raqqa, the ISIS capital, was.

Of course one must always be careful when constructing these comparison maps. In GIS software a polygon (Syria) whose coordinates are in degrees, dragged to another latitude, will be the wrong size. The safer method, which I have used here, is to make separate maps at the same scale of both Syria and the northeast US, and then juxtapose them in a photo editing program (the GIMP, in my case).  Syria is about 780 km from end to end on its long axis, about the distance from Boston to Richmond, Virginia.

This brings up the question of what it means to understand distance. In the original map, what does the distance from, say, Philadelphia to the eastern tip of Tennessee, mean to people living in Washington? I suspect they rarely go very far to the southwest from the city, and when they do they experience slow, twisting roads that have difficulty passing through the Appalachian mountains. So if the Washington-dweller says “It would take me six hours to drive to the tip of Tennessee!” is there anything meaningful at all there for his understanding of Syria?

I kind of prefer this juxtaposition, both for the similarity of desert landscape, and the sense of distance.

CIA inset LA and Damascus aligned

Los Angeles stands in for Damascus here, and Las Vegas finds itself somewhere up on the Euphrates. (“Las Vegas on the Euphrates” is not yet, but may someday be, the tourism slogan of Deir Ez-Zur or Raqqa.) If you are familiar with the distances and deserts of the American southwest, it is remarkable how much smaller Syria looks when shown like this.

Of course, as Appalachian Trail hikers know, there happens to be a town in Virginia called Damascus. It’s actually right there, just across the border from that eastern tip of Tennessee. So perhaps the best juxtaposition of all aligns Damascus, Virginia, with Damascus, Syria.

CIA inset Damascus and Damascus aligned