I’ve started experimenting with rendering representations of LIDAR-measured Digital Surface Model map datasets, which in contrast to Digital Elevation Model or Digital Terrain Model datasets - which are more common, and only consist of the raw terrain elevation data - have human-built structures in the height data (i.e. buildings). Previous Map Renderings (A full list of ones I’ve done so far can be found here) have involved DEM model data which just consists of the natural terrain, and in terms of scale, I’ve focused there on rendering entire countries or islands.
So I’ve been curious to try generating more detailed imagery of more localised areas, in particular of cities where human-made buildings and architecture are clearly visible, and here are some very early initial attempts with the raw DSM data for London.
The above render (Full 2.5k Render) is just a plain render of the surface being displaced by the LIDAR dataset values, using data from the UK’s Environment Agency from the LIDAR Composite DSM 2022 - 1m dataset showing central London.
(Note that when downloading the tilesets from the website, the category defaults to the ‘DTM’ Digital Terrain Model version which doesn’t have buildings, so ensure you switch to the ‘DSM’ version if that’s the one you want.)
When zoomed out and viewed from above, the fidelity seems pretty good: showing buildings, bridges and trees in nice detail.
Using a slightly less boring shader look - an occlusion shader, driving a diffuse surface colour gradient - gives what I think is a quite pleasing effect, accentuating the streets:
(Full 2.5k Render)
The height scale I’ve used is not physically-based to real scale (to the horizontal extents) currently: I’ve just eyeballed something which looks fairly reasonable, but I think it’s probably a bit too high still.
Once you start to look at the resulting generated surface in a bit more detail from closer up or at a more horizontal viewing angle, then the understandable limitations of the fidelity and format (2D single height values) of the data becomes a bit more apparent, especially in situations like overhangs, where single DEM/DSM values per point on a 2D surface are obviously limited. This can be seen in representations of bridges (there’s no gap underneath them), Tower Bridge and the London Eye in the below alternative view render:
London did seem to previously have some 25 cm resolution DSM datasets available at some point, but data.gov.uk doesn’t seem to have those any more, although they may still be part of other overall datasets, so I’ll try and find them. That won’t solve the ‘overhang’ problem (you’d need to use a full 3D point-cloud representation for that), but it should provide extra detail which might be interesting.
I’d also like to try and colour in water areas (and maybe foliage-heavy areas like parks) specific colours to add a bit more contrast and produce some more “artistic” versions, which I’ll look into doing over the next few months. Looking at rendering other cities that have more hilly terrain than London might be interesting as well, in order to have a combination of terrain and human-built buildings. I did look to see if I could find LIDAR DSM models of Wellington, NZ (where I’m currently living), which has an interesting combination of the two (although most buildings are quite small on the hills here), but I could only find DEM terrain models (without structures) for NZ. Cities like San Francisco are likely to have good data in this category, and I’ll have a look at other cities as well to see what’s available.
In a third instalment of attempting to copy images I’ve seen online with my own code, I recently saw some images generated by Johannes Kröger, whereby he ‘integrated’ or averaged a satellite image taken every day from the Suomi VIIRS Satellite into a final image which approximated the median average of cloud cover over the year.
He had an original blog post in 2019 here, and a follow-up in 2021 with more technical details here.
I liked the look of the imagery, and was curious how easy it would be to generate myself, and on top of that, was also interested in generating per-quarter/season images rather than ones only for the entire year, in order to try to see obvious variations between seasons.
It should be noted that these will be approximations: the source imagery is taken once a day - generally around noon (although it varies per day per location due to the satellite orbits, as can be seen when comparing adjacent per-day images) - and these processed imagery will include snow/ice cover as well, as shown in this preview of the North Pole area for the approximate average pixel colour of all 366 days of 2020:
Johannes Kröger’s 2021 blog post contained a bash script example which used the gdal_translate command of the Geospatial Data Abstraction Library (GDAL) suite of tools to download the source imagery from NASA’s Global Imagery Browse Services (GIBS) using a web API which provides tiled images, allowing the download of entire images from source tilesets.
I needed to modify his script to get it working (the ‘TileLevel’ needed to be changed), but I didn’t really want to use bash shell anyway, so I wrote a Python script to do the same thing, but added the functionality to also download imagery for multiple days at a time as a date range, and to also use multiple threads to download multiple images in parallel (downloading a single set of tiles for a single date is quite slow), and also added a ‘cubic’ resize filter to the gdal_translate command line args.
The Python script in its final form (albeit slightly sanitised - the save path will need to be changed in order to use it) can be downloaded here.
Note that the size of the images are quite large on disk given their fairly high resolution.
Johannes Kröger did give instructions on how to use available software (GDAL in his case) to perform the ‘averaging’ operation, but this was the bit I wanted to fully implement in code myself: I already had fairly comprehensive Image reading and processing infrastructure code of my own, so modifying it to perform ‘mean’ averages was pretty trivial: to just loop through the entire planar image of each final .tif file for each day’s imagery, and add them all together, and then divide each pixel value by the total number of images. This worked, however the result of using ‘mean’ average pixel values produces an image which does not really represent (at least directly) pixel values that actually occurred in terms of cloud cover: it’s an interpolation, and doesn’t show the pixel values that were most common (i.e. the colour values which occurred the most over the duration of the year for each day).
To find the most common pixel values over the course of the year for each pixel position in the imagery, the ‘median’ average needs to be used, and to calculate this was more development work, as a ‘median’ average requires having all the values for a pixel sorted in order (of luminance/brightness generally), and doing that for 365 16k images at full float32 precision in linear space (despite the source data being in sRGB 8-bit space, it’s generally a good idea to pull pixel values into ‘linear’ colourspaces in order to do computation on them) would take around 1.61 GB of memory (16,384 x 8,192 pixels x 3 channels x 4 bytes) per image, and so it was not going to be feasible to store all 365 entire planar images in memory at once (that would take at least 587 GB of RAM). I could have quantised the pixel values a bit whilst still keeping them in linear space (say to half float16 precision) or something lower with fixed-point, but that still wouldn’t have been anywhere near enough of a reduction, so it was clearly going to require breaking the image up into chunks, so I decided to process the ‘median’ average values in tiled regions, given the source TIFF images were in tiled form anyway, and so reading the individual tile regions for each source image would be easy and pretty fast to iterate through them.
I ended up with an algorithm that would for each tile region (256 x 256 size for the images I had downloaded) of all images (they all obviously have to be the same resolution), iterate through all images for the year, but just for that single tile region at a time, and accumulate all pixel values for all images into arrays per pixel position within the tile region. This way, the total memory usage was “just” the tile size dimensions (256 x 256) x 3 x 4 bytes = 786 KB (plus a bit extra for data structure overhead). Total memory cost for 365 tiles would then be around 287 MB, which is much more reasonable. Then for each pixel position within the tile region, all the pixel values for that pixel position from all the source images needed to be sorted (by luminance), and then the middle ‘median’ value picked. This single RGB value per pixel position within the tile region could then be baked down to a single final image buffer for the tile region, and the memory allocation of all previous pixel values for all 365 tile images could be freed, and then the next tile region could be processed in the same way, for all tile regions in the source images.
Then, finally, these per-tile-region final images would be re-assembled based off their tile position into a final image of the resolution of the original full source images, and this result saved to a final full output image. Given enough memory (my main Linux desktop has 32 GB), it was also trivial enough to process the per-tile-region reading of all source tile images for that tile region and ‘median’ sorting and evaluation in multiple threads, as each tile region could be completely independent from one another, speeding up processing considerably.
Processing 365 16k images into a final output image took around 14 minutes using 12 threads on a Ryzen 5900X, which wasn’t too bad, and there’s still a bit of room for further optimisation I think.
After experimentation with the output, I also added thresholding so as to not accumulate pixel values that were black: the poles of the earth in the satellite imagery were occasionally black, depending on the orbit and that affected the output values a bit.
I had tried to produce average images for the 2022 year, but it turns out the Suomi VIIRS Satellite was missing imagery for late July 2022 and the first half of August 2022, so I used 2021 instead which seemed to have full imagery, and also did 2020 for comparison purposes.
The output result of this for all days in 2021 is this image (full 4k version link):
which is using the same WGS 84 projection that the source images used. Reprojected to a more “true-size preserving” projection - Robinson Projection - provides this image (full 4k version link):
Below is a table containing links to 4k versions of per-quarter images for 2021:
The per-quarter versions do clearly show (as expected) obvious differences in seasons, although there may well be yearly variation as well, and snow/ice cover changes will also be included in the changes.
I’m keen to produce more of these in the future - ideally at higher resolution for more localised regions with better projections - in addition to attempting to generate (mostly) cloudless imagery similar to the famous “Blue Marble” images, to see how easy it is to detect clouds vs snow/ice on the ground: either with large vs. small changes day-to-day between images, or I wonder if it’s possible to use Infrared imagery to detect if colours are likely clouds or not, or by using some of the other output info from the VIIRS sensors.
Last week I spent a few days down in the South Island around Queenstown and Wanaka, mainly to try and photograph the Autumn Colours. Apparently the term for this is a “Leaf Peeper”!
I’ve been several times before, but always in Spring or Summer, so this was the first time I’d seen any of the South Island in Autumn.
Whilst there is a bit of colour in Autumn in the North Island in places (in parks north of Wellington for example), there’s not much that I’ve seen in general (at least in comparison to the Northern Hemisphere), so it was a nostalgic memory of Autumn in the Northern Hemisphere, seeing the very widespread colours that occur there and I grew up with.
I think I’m correct in saying that almost all of the trees in New Zealand that have leaves which change colour in the Autumn are “exotic”, non-native species that have been imported from Europe, North America and Japan, with there being very few native deciduous trees (most are evergreen), and of those very few of them actually change colour.
Maples, English Beech (the native New Zealand Beech species are evergreen from what I can tell), Gum, Cherry and Horse chestnut trees seem almost certain to have been imported by early settlers, and the fact they seem to mostly be found along rivers and near settlements rather than being found out further away in more remote parts of the landscape seems to re-enforce that theory, but it’s difficult to know for certain.
Regardless, they do provide a very noticeable splash of colour on top of New Zealand’s already beautiful scenery in the area, and I did get some nice photos I’m happy with.
I’ve been progressively trying to learn the Rust programming language for around a year and a half now and as well as porting some of my existing C/C++ apps I have to Rust as learning exercises (but with the C/C++ versions still mostly being the main ones for the moment) I have also started to write some new from-scratch apps in Rust when I think that makes sense, rather than defaulting to C/C++ as I previously would have.
I’d recently had the need to record some basic app process stats (CPU usage and RSS memory usage) over the duration of its running, and while there are existing applications out there - i.e. psrecord - that would largely do what I wanted, I was tempted to try and write my own version in a “native” language - at least for the recording part: the plotting / visualisation part is a bit more tricky. This was partly so as to have another small project with which to gain more experience with Rust, but also because I wanted additional features like the ability to control whether to have “normalised” or “absolute” CPU usage, and the ability to separate CPU time into user and system time - so I’ve written an initial version of psrec which is my own equivalent to psrecord written in Rust.
This initial version is pretty basic so far, and doesn’t yet support all the additional features I wanted: I’m making use of the psutil crate to extract the process information for the moment, but its functionality is incomplete as it doesn’t support extracting info on a process’ children, or info like the number of active threads a process has, which is functionality I will likely want at some point, so I’m no doubt going to have to get down in the weeds with the /proc/<pid>/ file system interface, which can sometimes seem a bit primitive and messy in my experience (it would be nice to have a first-class API to access the data efficiently, rather than doing string parsing, although it does make things very visual and easy to debug).
For the moment, I’m using Python and matplotlib to plot the resulting data, which to some extent is a bit at odds with writing the main app in a compiled language like Rust, but I think it’s okay for the moment, and it allows the recording app to theoretically be more efficient and low-overhead (although polling the /proc/ file system and recording samples every second isn’t really that much overhead), whilst using Python for things it’s very good at.
Below is a basic example of the chart plot of a quick process recording:
Last week I received a work-provided (on loan) new Apple MacBook Pro (14-inch, 2021) with the Apple M1 Pro 10-core processor, so I was very curious to see how it performed with some of my code, given the generally excellent reviews and reports I’d seen by people with it online.
I have my own MacBook Pro 15-inch 2015 (Intel i7) model which I’ll be comparing it against as a basic benchmark performance test, although it won’t be a completely fair test as the two machines are running different versions of MacOS and the clang compiler, and I didn’t want to upgrade my own machine’s MacOS version to a newer version (for various reasons). So the code being executed on the different machines will be generated by different versions of clang, but given the micro-architecture is totally different anyway, I think it’s still a meaningful test and at least will show the performance difference between the two laptops.
My old MacBook Pro 15-inch 2015 is a 2.5 GHz Intel Core i7, with four cores, eight threads, and 1600 Mhz DDR3 memory, running macOS 10.14.5 (Mojave), with Xcode 10.2.1 installed which has clang version 10.0.1.
The new MacBook Pro 14-inch 2021 is the M1 Pro 10-core version, with 8 performance cores and 2 efficiency cores, running macOS 12.3 (Monterey), with Xcode 13.4.1 installed which has clang version 13.1.6.
Tests:
I’ll be using two of my apps as benchmarks: my Mint interpreter language VM (originally based off Robert Nystrom’s excellent Crafting Interpreters Lox language tutorial) but with additional functionality and performance improvements, which I’ll use to benchmark two Mint scripts as single-threaded tests, and also my Imagine pathtracing renderer, which has native SSE intrinsics support for Intel and native Neon intrinsics support for ARM, which I’ll run in both single- and multi-threaded scenarios.
Both apps will be compiled with -march=native on the Intel side and -mcpu=native on the Apple M1 Pro ARM side, using the clang version on the machine in question, as well as optimisation level: -O3.
The two Mint script tests will be loop value calculation as Test 1:
var a = 1;
for (var i = 0; i < 100000000; i += 1)
{
a = (i + i + 3 * 2 + i + 1 - 0.42) / a;
}
print a;
and a variation of Project Euler 21 to calculate the sum of all Amicable numbers under 15,000 as Test 2.
The Imagine rendering tests will render the below example scene of a basic maze with 12 spherical area lights and one hemisphere light, with four Next Event Estimation samples being taken each bounce to evaluate the lighting, and using Cone sampling to perfectly sample each spherical area light from each hitpoint, with 256 samples-per-pixel being done, using stratified sampling and a ray depth of 6 diffuse bounces. Depending on the setup of the tests (single-threaded, multi-threaded, etc), I’ll render the scene at resolutions of: 150 x 150, 300 x 300 and 600 x 600.
For all tests, I’ll run them on battery power and then mains power, in case that affects the clocking of the processors, and will also wait between test runs for the processor temperatures to be below 50 degC before running the next test, to try and reduce the impact of thermal throttling. However, for the longer-running multi-threaded Imagine rendering tests I’m fairly sure thermal throttling will start to take place on the Intel 2015 machine given how fast the fans start to spin in longer-running scenarios.
All tests will be run three times, and results below will show the mean average of those numbers.
Results
Single-threaded Mint VM interpreter:
Single-threaded Mint interpreter VM benchmarks, smaller values are better:
For these single-threaded tests, the M1 Pro has fairly modest (~11% faster in Test 1, ~16% faster in Test 2) improvements over the seven year old Intel i7, however as the Mint VM execution is often branch-prediction constrained within the main VM bytecode interpreter loop, I expect this doesn’t really allow the M1 Pro’s eight-wide design to really stretch its legs running these benchmarks, as there’s a limit to the amount of Instruction-Level Parallelism that’s achievable.
It’s interesting to see that there’s no apparent difference in performance between running on battery vs running connected to the mains: I was expecting to see a difference - even if only on the Intel i7 MBP - as I know my Thinkpad T480s running Linux shows significant performance differences between the two states, but Apple clearly don’t seem to be implementing that strategy.
Single-threaded Imagine rendering:
Single-threaded Imagine rendering benchmarks, smaller values are better:
Things start to look more impressive for the M1 Pro in the single-threaded rendering tests, as here there are a lot more floating point calculations being done with a lot less branching, and there’s a fair amount of SIMD use (i.e. ray intersection / traversal code, which uses native SSE or Neon instructions for the respective platform architectures) allowing the M1 Pro to show what it can do.
For both the 150 x 150 resolution and 300 x 300 resolution renders, the M1 Pro machine completed the tasks in ~46% less time, so it’s close to twice as fast as the older Intel machine. Again, there was no apparent performance difference between running powered off the battery compared to running with the machines plugged into the mains power.
Multi-threaded Imagine rendering:
Multi-threaded Imagine rendering benchmarks, smaller values are better:
The multi-threaded rendering tests show that for the shorter-duration 300 x 300 resolution render, the M1 Pro was more than three times faster (render time was less than a 1/3 of the Intel MBP’s i7 render time), and for the longer-duration 600 x 600 resolution render this ratio was slightly better still - I suspect because of the fact the Apple M1 Pro has significantly better thermal performance than the Intel i7, so is not throttling as much - if at all - although watching the temps on the M1 Pro machine towards the end of the test shows it did reach 84 degC, so it’s possible throttling is happening a bit.
Conclusion
So for most code that’s likely not branch/dependency-bound, it would seem the almost one year-old M1 Pro being tested here is around three times faster than the seven year-old Intel i7 model (although the base M1 Pro model with two less performance cores will not be as fast), and for multi-threaded code that parallelises well it can be more than three times faster, with much better thermal performance (which means less fan noise) and battery life, which is exactly what you’d want from a laptop, and I think shows that the Apple M1 processors are something to be excited about.
Over the past year I’ve become somewhat interested in visualisations of maps and terrain representations, and in particular “artistic” renders of relief maps. I’ve tried making a couple of different types of renderings of maps in some form or another over the past few years with varying degrees of success, and over the past few months I’ve been progressively creating what I’m terming “Minimalist” relief maps, where a surface is displaced by a Digital Elevation Model terrain height map (often to an exaggerated extent compared to real-word scale) with fairly simple shading, but using light and shadows as a key element to provide a sense of the terrain topology and height in a stylistic way.
These maps are a lot less manual-labour intensive to set up than the previous historical topographic ones I tried which required a lot of manual warping of the DEM images to match up with the historical map, so there’s a lot less effort (on my part anyway) to generate these more minimalistic style ones, and I still find them very nice to look at.
I’ve experimented with various different types of shading, with the main two types I settled on being the diffuse colour mapped to a colour gradient, driven by:
the terrain height value (the Y coordinate height in worldspace)
the occlusion ratio in a hemisphere around the shading point
The latter of which I think I like best (and is what both renders above show) as it allows nice highlighting of the edges and gradients of steeper terrain, although it means the stronger colours are generally in darker more occluded areas which hides it a bit, and the renders are slower as well, as additional raytracing to calculate the surface point occlusion needs to be done during shading.
I’ve just got back from a fantastic week and a bit in the South Island, and was very lucky to have perfect weather for most of the trip, as well as getting a lucky view of the Aurora Australis (Southern Lights) one night.
I was in Queenstown at the time, and a phone app I have which gives notifications when Aurora might be visible at your location alerted me that it was likely visible, so I grabbed my camera and tripod to try and photograph it, starting during the later stages of dusk. There were two different “types” that were fairly faint but just visible to the naked eye: a green glow that really looked like light pollution - although there were no large towns in exactly that direction, and it morphed over time changing shape - and blue streaks, gently pulsing over time.
With a long exposure, they are much more clearly visible.
A very magical experience.
I also managed to get some excellent views of Aoraki / Mt Cook.
The above render (Full 2.5k Render) is just a plain render of the surface being displaced by the LIDAR dataset values, using data from the UK’s Environment Agency from the LIDAR Composite DSM 2022 - 1m dataset showing central London.
(Full 2.5k Render)
