Find and download the bare earth (digital terrain model = DTM) and full feature (digital surface model = DSM) LiDAR zip files for a quarter section of your choice from the LiDAR viewer (I used the files here and here). To familiarize yourself with the concepts of DSMs and DTMs, read this answer.
Download LAStools and add the LAStools toolbox to ArcMap (or ArcCatalog)
Run the las2dem tool within the LAStools toolbox for each of your downloaded .las files. You can accept all default parameters.
The tool should output a .tif image for each surface, projected to UTM Zone 11. If not, run Define Projection to assign the projection.
The names of each layer will be the same, by default (e.g. “710622NW.tif”). Change them to be unique, so geoprocessing tools (and you) can easily tell them apart (e.g. “710622NW_BE.tif” and “710622NW_FF.tif”).
The DTM and DSM each report elevations in meters above sealevel. Since we are interested in the difference between those elevations (i.e. the difference between rooftop and the ground pixels) we need to subtract the two images. You can do this using the Minus tool, or more generally, Raster Calculator. I used the expression: Int(“710622NW_FF.tif”-“710622NW_BE.tif”) to return the difference between the DSM and DTM, as integer numbers (this will simplify things later on).
Eventually, we are going to want polygons to load into Cesium (on the order of hundreds or perhaps a few thousand), not individual pixels (there are 700,448 pixels in my image). At this point, we have what is called a classified image, classified by height in meter increments. We want to group this classified image by similar height values, and we can do that by processing our classified output. Following the instructions in this link, we are guided through filtering, smoothing, and generalizing our output.
Cesium is going to want lat/long coordinates, so you may as well run the Project tool on your polygons to switch the data to WGS84.
As far as I know, ArcGIS does not support GeoJSON (which we are going to use in Cesium). Conveniently, QGIS does. So, load your polygons in QGIS and “Save As…” GeoJSON.
Upload your GeoJSON to your host.
Switching to Cesium, create a map and load your GeoJSON (general example here). Cycling through the GeoJSON polygons, assign the polygon height to 0, and the extruded height to match the elevation attribute in your data (for me, that was stored in the “GRIDCODE” field).
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.
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.
Open the Potential Incoming Solar Radiation tool (Modules -> Terrain Analysis -> Lighting -> Potential Incoming Solar Radiation). Here is how I set up my run:
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:
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.
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: firstname.lastname@example.org
Load your orthophoto into your GIS program (I’m using ArcGIS, but you can use QGIS for the time being [it’s free]).
Create a new polygon feature class and draw a polygon to narrow down your area of interest.
Download LAStools (it’s free for the tools you need). Bring the LAStools toolbox into your GIS (there are ArcGIS and QGIS versions).
Run the lasclip tool to ignore the millions of LiDAR points that you’re not interested in.
Run las2dem to create an elevation raster of your LiDAR data.
Consult your favourite source, like the Farmer’s Almanac, to determine the timing for your growing season, .
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):
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.
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.