Starting with what we had done by the end of the last post, I would like to continue on to doing some actual analyses using QGIS, GRASS, and my other GIS-type tools. This post has been a long time coming, partly as always due to my laziness, but also because I was having a difficult time getting a lot of my software up and running on my “new” computer. Rather than trying to do everything from inside of QGIS, I decided to work separately in QGIS and GRASS(mostly…).
Although I’m not confident in the realism of the elevations file I composed in the last post. The elevations are quite low. Univariate analysis shows that the range of values is -3m(I’m not going to concern myself with what’s under the oceans except to avoid problems with basin fill, flow algorithms and the like) to 971m(which is genuinely very low). I’m not sure whether it needs to have less emphasis on the 1.0 exponent values or not(as it is, the low areas appear very relatively rugged, but until I have higher altitudes(5-9km), its difficult to say). At the very least, I probably need to add a large, high-exponent pass(perhaps 5-8,000m of exponent 5.0). Anyway, I’m using it as is, for now. We’ll see where this gets us.
Okay, SAGA’s hydro-modeling tools look very appealing, but I can’t set up Boot Camp with Windows XP, I’m in no hurry to buy a modern version of Windows, SAGA works poorly on WINE(on my Mac – I’ve heard tell of it working well, but not around here…), and the QGIS-internal versions of SAGA, GRASS, etc. have not been working at all well for me. Grrr… Looks like I’m back to the drawing table.
The most important thing I want is flow mapping that takes differential rainfall into account. Next, I want d-inf flow mapping. After messing around with a lot of GRASS’s raster hydrology modules, and producing a lot of unsatisfactory results, I settled on the r.watershed module. This requires an input elevation raster. I have that. It also uses a “flow” raster, which is described in the manual thus,”flow Input map: amount of overland flow per cell. This map indicates the amount of overland flow units that each cell will contribute to the watershed basin model. Overland flow units represent the amount of overland flow each cell contributes to surface flow. If omitted, a value of one (1) is assumed. “. Yeah. This sounds good. Unfortunately, it doesn’t do Tarboton’s d-inf flow mapping. That’s unfortunate, but it does provide both multiple flow direction and deterministic 8-directional flow models, so I’m hopeful I can at least create better incise flow erosion than Wilbur.
The available water map(rainfall-evaporation, my model doesn’t even try to deal with variable runoff infiltration) was generated using a method derived from this. Roper’s sech(lat) model looked good, but first I had to convert from watts per square meter of energy absorption/release due to evaporation/precipitation to meters of rainfall/evaporation per year. The input flow map only really needs relative contributions of each cell to overland flow, but having proper units seems like it might be useful down the line.
Once I figured this out(correctly, I hope), I entered the following formulae into r.mapcalc:
evaporation = 2.307 / (1 + ((y() – 15)/17)^2 + 2.307 / (1 + ((y() + 15)/17)^2)
This would generate the “evaporation” raster map.
precipitation = 1.957 / (1 + ((y() – 33)/15)^2) + 6.292 / (1 + (y()/4)^2) + 1.957 / (1 + ((y()+33)/15)^2)
This would generate the “precipitation” map. As my maps were already in latlong format and correctly located(more or less: more about that if I ever get SAGA working…), the y() internal variable represents latitude.
Looking at the results, rainfall seemed awfully heavy under the subtropical high pressure zone, but I decided to go with it for now. At least the net_moisture(generated in r.mapcalc with the formula: net_moisture = precipitation – evaporation) seemed properly dry under the STHZ, although it actually seemed a little too dry generally and the width of the desert band(r.mapcalc desert = net_moisture < 0.0) seemed… excessive. I’m following Carl Davidson’s climate-modeling efforts with great interest, but I’ll play with this for now. Later biome development might demand knowledge of precipitation and temperature(and by extension, evaporation) for summer and winter, rather than just annual averages and totals, but this will have to do for now. I’m largely just trying to see how this will work.
In an attempt to further refine the effect a bit I created a further version of the map with slightly higher precipitation at higher altitudes ( modified_precipitation = precipitation * (1 + 0.001*elevation*sin(2*y())) ). Then I added a bit of precipitation to a blurred version of the sea areas( coastal_precipitation = modified_precipitation + 0.05*sea_raster_blurred ). The exact details of developing that sea raster are embarrassing and given the limited quality of the results, I’m not going to share them. Suffice to say, you can generate an initial sea raster with r.mapcalc( sea_raster = elevation<= 0.0 ), you can dilate the result with r.neighbors and a method of ‘maximum’, and blur the result of that with r.resamp.filter with filter “gauss,hermite”. In retrospect, I probably should have blurred the elevation raster as well. Ehhh…
Thinking about that coastal_net_moisture map, I don’t think I want it propagating
negative flow contributions through and out of desert areas. Logically, it would be nice to have the negative values taken out of flow passing through a cell(evaporation), but having negative flow accumulations propagate into neighboring cells would not be reasonable. So I think I need to create a version of coastal_net_moisture floored at zero. In r.mapcalc, use the formula r.mapcalc expression=”coastal_net_moisture_floor = ( coastal_net_moisture@PERMANENT >= 0.0 ) ? coastal_net_moisture@PERMANENT : 0.0″. That won’t take into account extreme evaporation of the stream itself, but it should still be more valid than dealing with negative flows.
In addition, I really don’t want to calculate stream flows under the ocean. The submarine areas are not accurately represented anyway, as I did not have bathymetric data available, and for convenience simply set the ocean areas to -3 meters on the elevation map. So, here, I will create a version of the basin filled elevation map with all submarine areas set to null. There is definitely a module to do this, but I forget what it is. Anyway, the Map Calculator is sooo versatile, and it’s good to know how to get around in there. I should research more standard workflows, but for now I will use the formula r.mapcalc expression=”elevation_set1_basinfill_land = elevation_set1@PERMANENT >= 0.0 ? elevation_set1_basinfill@PERMANENT : null()”.
I already did all of this before I started writing this. Call it a flash-back. On to r.watershed!
I set the input elevation raster map toleration_set1_basinfill_land, and the input flow to coastal_net_moisture_floor. Minimum size of exterior watershed basin I will leave at default(though I may need to rejigger it if the results aren’t satisfactory, the manual says it’s a sensitive parameter😟), the rest of the inputs I’ll leave at default or off.
Next we set the outputs. The accumulation output I will name accumulation_MFD. This is to differentiate it from the possible future accumulation_D8 I might create in the future. I will also create topographical_index_MFD, stream_power_index_MFD, and all the rest I will name the same as the parameter name with an _MFD subscript appended. The stream output parameter I will name stream_segments_MFD for clarity. I will use the ‘b’ option to beautify flat ares, and I will leave the convergence at its default 5 for now.
I’m pretty satisfied with the results. I have a few good streams heading through the desert. I have quite a few decently long, but not too crowded streams running through the moister northern regions. And the wet, but rugged deep southern peninsula has several short streams. The cumulative flow in the desert regions goes negative even where streamm segments are found, which is not a little bit odd. My best guess is because the MFD is routing some flow away from the streams and the stream segment generator is trying to force downhill flow to go all the way to the ‘sea’. I don’t see any nasty straight segments, which is awesome. The desert regions, though they seem overly large to my eye(probably due to the weaknesses inherent in my very lame attempt at a precipitation and evaporation model, are pleasantly devoid of streams. Because this was generated from a real world DEM of a fairly moist area, the deserts are pretty extensively incised with stream-eroded features. This is a downside of using re-purposed real-world elevation data. It might work well to simplify the elevations with gaussian blur and apply erosion to the result. One good tool, which I really don’t have, would be a good aeolian erosion filter. Convert some of the flatter desert regions into saharan fields of marching barchans.
Some of the good behaviors, as well, probably, as some of the bad behaviors could easily be lain at the feet of the MFD model rather than a single deterministic flow direction. For better or worse, some can be blamed on the use of real-world data, even heavily massaged. A better test might be to use this on generated noise-based elevations.
I neglected a lot of the layers that were created in this process to create the images shown here. Also, beyond the narrow band of yellow to represent the desert regions, there is no climatic data here. Only a shaded relief based on the un-basin filled elevation map, a nicely massaged desert map in yellows at 25% opacity, a basin-filled land elevation-colored map at 43% opacity, a landmass vector map at 50% opacity, and the generated stream segments map at full opacity, with a monochrome blue color table are shown overlaid.
Future directions in research would be a much(much, MUCH, MUCH)better climate model, use of less natural initial elevations, use of accumulated flow or stream power index raised to a fractional power and possibly blurred to erode those elevations into a better approximation of naturalistic shapes, and you know, some actual biome data generated from climate and elevations. With this, I might be able to create something akin to the ‘satellite view’ of tectonics.js, only hopefully a bit more aesthetic.
I hope you found this to be a good introduction to using some of the available free GIS tools, and that you will download GRASS and QGIS and other open source GIS tools and use them to help your worldbuilding efforts. Thank you for reading,