Solar Radiation: free method

This is a follow-up post to Tuesday’s Garden Placement Using Publicly Available LiDAR, which relied on ArcGIS’ Spatial Analyst extension. Huge hat tip to the LAStools Linkedin conversation, and particularly Daniel Grohmann, for mentioning the freely available SAGA GIS tool, Potential Incoming Solar Radiation. I should note that this is truly my first attempt at using SAGA, so apologies in advance if I’m leading you astray in some way.

To recap, we had created a digital elevation model (DEM) from publicly available LiDAR from the City of Prince George, and were ready to calculate solar radiation values for a city block in order to optimally locate a garden plot.

  1. Download, install, and open SAGA GIS (it’s free)
  2. Import your DEM into SAGA. I used the GDAL: Import Raster module (Modules -> File -> GDAL/OGR -> GDAL: Import Raster). Double-click the layer in the Data workspace window to view the raster and ensure it imported correctly.
  3. Open the Potential Incoming Solar Radiation tool (Modules -> Terrain Analysis -> Lighting -> Potential Incoming Solar Radiation). Here is how I set up my run:
  4. After changing the colour ramp to the default (blue -> red), we are left with a very similar, if not identical, solar radiation map comapred to that produced by Spatial Analyst:

Garden Placement Using Publicly Available LiDAR

Here’s a quick method for finding the sunniest spot on your property in Prince George, BC (and elsewhere if you can find the data). Unfortunately, this method is not free from start to finish (it requires ArcGIS Spatial Analyst), but there are a fair number of free tools and data used.

  1. Get the data. Specifically, you need the raw LiDAR (.las file) for your area, which are not available for download from the City of PG’s Open Data Catalogue. While you’re at it, also ask for the orthophotos (aerial photos) and cadastral data (parcel boundaries, road files, etc.). Get in touch with the City of Prince George’s GIS department to make your data request:
  2. Load your orthophoto into your GIS program (I’m using ArcGIS, but you can use QGIS for the time being [it's free]).
  3. Create a new polygon feature class and draw a polygon to narrow down your area of interest.Image
  4. Download LAStools (it’s free for the tools you need). Bring the LAStools toolbox into your GIS (there are ArcGIS and QGIS versions).
  5. Run the lasclip tool to ignore the millions of LiDAR points that you’re not interested in.
  6. Run las2dem to create an elevation raster of your LiDAR data.Image
  7. Consult your favourite source, like the Farmer’s Almanac, to determine the timing for your growing season, .
  8. Here is where I used a Spatial Analyst tool called Area Solar Radiation, which is not free (in fact, it’s darn expensive). Run the tool, using your latitude (PG is about 53.914), and frost free start and end dates (for PG, June 4 to Sept. 3). You should end up with something like this, which my legend tells me ranges from blue (353 watt hours per square metre) to red (883061 WH/m^2): Image
  9. Now, it’s a matter of overlaying parcel boundaries, finding your property, and seeing what kind of sunlight you can expect. I’ve circled a few good candidates on this block with lots of sun in the backyard, and (surprise!) some of them are existing gardens.ImageImage

Bonus: another cool thing you can do (for free) is load your monochrome elevation model in Blender, extrude the terrain heights, drape the aerial imagery overtop, and create a 3D animation like this one.

PaperJS Voronoi Diagrams


I’ve been tinkering around with Voronoi diagrams, inspired by the example on the PaperJS site. My example (here) uses Raymond Hill’s Voronoi Javascript implementation and PaperJS to derive the diagram from a point distribution. Note: the redraw animation is a bit backwards, as the point locations are determined first, driving the placement of the Voronoi diagram.

Another, more festive example, can be found here, which shows the distribution of lights on our Christmas tree. Falalalala-la-la-la-la.

Refractometer: ABV Calculator


I got a mysterious package in the mail yesterday – by the time I figured out that it was a Christmas present from my brother, it was too late.

Turns out the gift was a refractometer, a device used to measure an index of refraction, or how light gets bent going through another substance. Few people outside laboratories have much use for such a thing, but along the way some science-savvy beer brewer realized that you could measure the amount of dissolved sugar, and thus the amount of potential alcohol, in freshly brewed beer. By comparing the amount of sugar at the beginning and end of fermentation, you can infer the amount of alcohol by volume (ABV) in the finished beer.

Most homebrewers use a hydrometer in a similar manner, but this requires several chances of contamination (each time the wine-thief is filled, often several times) and wastes a good deal of beer along the way (about 100ml for each measurement). Refractometer measurements use a mere couple of drops which are collected only once per measurement.

Unfortunately, once alcohol is present in the liquid (i.e. any time after the initial measurement), it messes with the refraction measurement and a mathematical correction must be applied in order to compare the two samples. Using formulas from Northern Brewer, Primetab, Realbeer, I made an interactive calculator for calculating ABV from initial and final refractometer readings (in Brix units), and you can find it here. Drag the bars up and down to match your initial and final measurements.

Notes: there are tons of other online refractometer calculators, and I’ve seen some that use different formulas (or at least different coefficients), so beware. Also, if you’ve got Matlab, my brother made a library for charting and doing the calculations for yourself (blog post here).



I don’t usually write about what goes on at our staff meetings (“what happens at the staff meeting, stays at the staff meeting”), but today the possibility of a staff lottery for the prime parking spot behind the office piqued my interest (and I don’t think this a trade secret). So, I spent a grand total of five minutes finding the original KineticJS code (sorry, HTML5 browsers only), modifying it with staff initials, and uploading for live use the “EDI Wheel of Parking”. Give it a spin here.

Least Cost Path: Dijkstra’s algorithm


I’ve recently started walking to work, and as a result I’ve started thinking about the easiest way to get there. I’m not sure how much other people think about it, or consciously think about it, but every step and other larger movement decision (where to cross the street, what neighbourhood to visit, etc.) has an associated cost (effort, danger, time, etc.). When you find the optimal route according to your cost assignments, you have found what’s called the Least Cost Path (at least, by ArcGIS), or solved the Shortest Path problem in Graph Theory (which I won’t pretend to fully understand).

Your brain tends to do this fairly easily. We know that it’s easier and safer to walk down a sidewalk than the middle of a street or over top of a building. Perhaps you choose the north sidewalk because it gets more sunlight than the south sidewalk. We also avoid certain routes because they are poorly lit at night.

Our computers have more trouble “just knowing things.” Using GIS software, we can assign costs to a continuous cost surface. In my example, I’ve assigned costs to a section of Prince George, BC (spatial data courtesy of the City of Prince George Open Data Catalogue) as follows: building = 100, street = 15, alley = 5, alley/sidewalk intersection = 3, and sidewalk = 1. The higher the cost, the less you’d want to walk there.

Next, we need to evaluate our route options. You can do this in ArcGIS using the Spatial Analyst toolbox (which ain’t free), but it’s cheaper and more fun to make our own using javascript, Google Maps API (the map), PaperJS (access pixel values from cost raster), and Dijkstra’s algorithm (the math). I won’t go into the details, but you can read the code (Ctrl+u) and try to figure it out for yourself.

My example is here. A word of warning: it’s slow. It takes about 80 seconds to evaluate 1000 potential points in a 10 pixel grid (that is, 10 pixels in the original raster, before it gets blown up to fit the map). The example starts in the top left corner and tries to travel to pixel 250 (right), 250 (down). Potential, rejected sites are grey circles; selected sites are purple filled circles connected by the least cost line.