Topographic Exposure Index

Here is my attempt at creating a Topographic Exposure Index. Topographic exposure or TOPEX essentially looks out from a given pixel and scans the surrounding landscape for maximum angles from the given pixel on a digital elevation model. Maximum angles are measured for each of 8 horizontal direction classes and those maximums are summed for the final score. A high TOPEX value would indicate a very sheltered area whereas a value of 0 would indicate you were on a flat plain or mountain top.

I ran it on the Blue River Watershed in Oregon and it took it about 5 days to complete. The input DEM was 7084x5882 pixels and the search kernel looked out 2,500m. The main innovation relative to other implementations I have seen is I created a kernel with concentric rings to scan the landscape for maximum ridge angles within 8 horizontal direction classes. The example below has 50 concentric rings which are separated by 50m. The DEM was at a 5m resolution so I went with a relatively fine spacing for the concentric rings.
 

	#-------------------------------------------------------------------------------
	# Name: Calc_TOPEX.py
	# Purpose: This program will calculate the vertical angle for a given direction and distance
	# on a given DEM. The maximum angle is calculated for each of 8 directions and summed. A greater angle
	# equates to less exposure.
	#
	# Author: olsenk
	#
	# Created: 17/02/2014
	# Copyright: (c) olsenk 2014
	# Licence: <your licence>
	#-------------------------------------------------------------------------------
	#!/usr/bin/env python
	import arcpy
	import numpy as np
	#from scipy.ndimage import filters
	from math import atan2, atan, degrees
	nodataVal = -999
	# funciton to construct a circular kernal from a given radius
	def buildAzimuthKernel(radiusMeters, referenceDesc):
	kernel = []
	# build kernel from radius and input raster cell size
	window = int(radiusMeters / referenceDesc.meanCellHeight)
	for row in range(-window, window + 1):
	for col in range(-window, window + 1):
	Xdist = abs(col * referenceDesc.meanCellHeight)
	Ydist = abs(row * referenceDesc.meanCellHeight)
	totalDist = pow(Xdist**2 + Ydist**2, 0.5)
	if (totalDist > 0) and (totalDist <= radiusMeters) and (totalDist % 50 == 0):
	angleDeg = degrees(atan2(col,row))
	if (angleDeg > -22.5 and angleDeg < 22.5):
	# north
	azClass = 0
	elif (angleDeg >= 22.5 and angleDeg < 67.5):
	# northeast
	azClass = 1
	elif (angleDeg >= 67.5 and angleDeg < 112.5):
	# east
	azClass = 2
	elif (angleDeg >= 112.5 and angleDeg < 157.5):
	# southeast
	azClass = 3
	elif (angleDeg >= 157.5 or angleDeg < -157.5):
	# south
	azClass = 4
	elif (angleDeg >= -157.5 and angleDeg < -112.5):
	# southwest
	azClass = 5
	elif (angleDeg >= -112.5 and angleDeg < -67.5):
	azClass = 6
	else:
	azClass = 7
	kernel.append([row,col,int(totalDist + 0.5),azClass])
	print str(len(kernel)) + " pixels in kernel"
	return(kernel)
	# -------------------------------------------------------------------------------
	def main():
	demImage = "T:\\groups\\clib\\spatial\\HJA\\phy\\dem5m_merge_b"
	outRasterName = "N:\\projects\\HJA\\topex\\hja_topex"
	outTestRasterName = "N:\\projects\\HJA\\topex\\hja_topex_s"
	logFile = open("N:\\projects\\HJA\\topex\\topexLog.txt", 'w')
	# read DEM into memory
	referenceDesc = arcpy.Describe(demImage)
	lowerLeftCorner = arcpy.Point(referenceDesc.Extent.XMin + (referenceDesc.meanCellWidth / 2), referenceDesc.Extent.YMin + (referenceDesc.meanCellHeight / 2))
	myArray = arcpy.RasterToNumPyArray(demImage, lowerLeftCorner, referenceDesc.width, referenceDesc.height, nodataVal)
	demImage = np.ma.masked_values(np.rint(myArray), nodataVal)
	# set workspaces
	arcpy.env.workspace = "c:\\temp"
	arcpy.env.scratchWorkspace = "c:\\temp"
	# pad DEM with nodataVal to avoid boundary issues
	## pad = np.ma.zeros((500,referenceDesc.width), dtype=np.int32)
	## demImage = np.ma.vstack((pad, demImage, pad))
	## referenceDesc.height += 1000
	## pad = np.ma.zeros((referenceDesc.height, 500), dtype=np.int32)
	## demImage = np.ma.hstack((pad, demImage, pad))
	## referenceDesc.width += 1000
	kernel = buildAzimuthKernel(2500, referenceDesc)
	# create empty image to hold topex values
	topexImage = np.ma.zeros((referenceDesc.height, referenceDesc.width), dtype=np.float32)
	topexImage.mask = demImage.mask
	# run topex function over image
	# NOTE: input DEM must have nodata buffer same width as kernel radius to avoid using boundary if statements (below)
	## i = 0
	for y in xrange(3001,referenceDesc.height):
	for x in xrange(referenceDesc.width):
	if (demImage.data[y,x] != nodataVal):
	exposure = [0,0,0,0,0,0,0,0]
	for pixelOffset in kernel:
	# if ((y + pixelOffset[0]) >= 0 and (y + pixelOffset[0]) < referenceDesc.height and (x + pixelOffset[1]) >= 0 and (x + pixelOffset[1]) < referenceDesc.width):
	if (demImage.data[y + pixelOffset[0], x + pixelOffset[1]] != nodataVal):
	exposeRatio = (demImage.data[y + pixelOffset[0],x + pixelOffset[1]] - demImage.data[y,x]) / (pixelOffset[2])
	if exposeRatio > exposure[pixelOffset[3]]:
	exposure[pixelOffset[3]] = exposeRatio
	# set topex values to sum of 8 direction exposure degree maximums
	exposure = np.degrees(np.arctan(exposure))
	topexImage.data[y,x] = sum(exposure)
	# testImage.data[y,x] = exposureAz[0]
	print "Processing row: " + str(y)
	logFile.write("Processing row: " + str(y) + "\n")
	## if (y % 1000 == 0):
	## i += 1
	## outRasterNameNew = outRasterName + str(i)
	## try:
	## lowerLeftCorner = arcpy.Point(referenceDesc.Extent.XMin, referenceDesc.Extent.YMin)
	## logFile.write("Saving numpy array to raster " + str(i) + "\n")
	## rasterToSave = arcpy.NumPyArrayToRaster(topexImage.filled(nodataVal), lowerLeftCorner, referenceDesc.meanCellWidth, referenceDesc.meanCellHeight, nodataVal)
	## rasterToSave.save(outRasterNameNew)
	##
	## # define the projection based on the reference raster
	## arcpy.DefineProjection_management(outRasterNameNew, referenceDesc.spatialReference)
	## except Exception, e:
	## print "\n\n" + e.args[0]
	##
	## except:
	## print "unhandled Error!!"
	# convert the masked array to normal numpy array and fill the masked area with nodataVal
	try:
	logFile.write("Saving numpy array to raster\n")
	rasterToSave = arcpy.NumPyArrayToRaster(topexImage.filled(nodataVal), lowerLeftCorner, referenceDesc.meanCellWidth, referenceDesc.meanCellHeight, nodataVal)
	rasterToSave.save(outRasterName)
	# define the projection based on the reference raster
	arcpy.DefineProjection_management(outRasterName, referenceDesc.spatialReference)
	# arcpy.BuildRasterAttributeTable_management(rasterOutputName, "Overwrite")
	except Exception, e:
	print "\n\n" + e.args[0]
	except:
	print "unhandled Error!!"
	logFile.close()
	print "done."
	if __name__ == '__main__':
	main()

view raw Calc_TOPEX.py hosted with ❤ by GitHub

The ESRI International Developers Summit

Well, I'm dipping my toes into the world of web map development. Here is a rundown of some of the things I was exposed to at the ESRI developers summit in Palm Springs, CA. I skipped the Plenary session so here is a link to the videos.

Joshua Tree National Park

The keynote by Chris Wanstrath was excellent. As the co-owner of GitHub Chris talked about his experience developing new businesses and also the frontier of collaborative code development. He did not really evangelize about open source (GitHub has a pay per private option), but rather described how open sharing of code is useful to private and public operations. GitHub shares code by community 'forking' the master branch code. The new branch is cloned to your desktop and then improvements are synced back up to your account. Once you are happy with the changes you make a pull request to the original developer and the code is probably merged back into the main branch. Forking used to be an insult, now it's a declaration of love.

Some sessions I attended:

1. Intro to Javascript preconference seminar. This two day workshop was a little slow at first, but was ultimately very helpful. I learned all about manipulating the DOM and building pseudo class objects in JS. Namespace is tricky as anything outside of a function is a global variable. We learned how to nest functions to shield their namespace. Finally, we did a little intro into DOJO which extends the JS framework.
2. ArcGIS API for Javascript. There is a wealth of info at developers.arcgis.com. Many samples and there is a new sandbox where you can tryout playing with your data. Many styling and widgets (DIJITS) were also in the samples. I was able to hook into the HJA rest service easily. Some suggested IDEs were Aptana, Sublime, Webstorm, Eclipse, or Visual Studio. Code cleaners for JS include Lint and Hint. Also checkout Grunt and Gulp.
3. ESRI on GitHub. ESRI has many open source projects on GitHub. David Spriggs did a great session on DOJO the good parts. He has a configurable JS API viewer on GitHub.
4. iOS SDK. You need a mac for development. It appears to be written in Xcode. Over my head.
5. Android SDK. There is an Eclipse plugin for Android. Checkout developer.android.com/tools. The best emulator is called Genymotion. Android development is in JAVA. It allows for async tasks that can run in the background. 
6. Web vs. Native app development debate. Pros for web apps: One app to rule them all, smooth delivery of tiles, canvas blending, faster debug time and updates, cheaper to develop. Pros for native development: access to device features such as touch gestures, disk, accelerometer,  text. Also, native can work offline, process faster and access multiple monitors.
7. ESRI and Leaflet. This session was done by Patrick Arlt (@patrickarlt). His notes are here. Leaflet is simple and clear. It has 100's of plugins for labeling, charting, heatmaps, ect. It can load tiles from ArcGIS server. Works great to integrate into existing apps. Works with non-web Mercator projections. ESRI leaflet is on GitHub. There are multiple projects including a geocoder and demographic popups. Cool project by Code For America in Chattanooga.
8. ArcGIS Pro, the road ahead. This seems to be a branch off of arcMap. It will have integrated 2D and 3D capabilities, multi-threaded draw engine, 3D editing, contextual menu ribbon, 64 bit, and spreadsheet like editing of attributes and their properties. It will have more CAD like options, 3D hotspot analysis (density and frequency) and is available at 10.3. Python scripts will run but ArcObjects and not exposed.
9. Big data visualization with Hadoop. Hadoop is an open source project that can reduce 100s of billions of points with MapReduce. Also introduced was YARN (yet another resource negotiator) and Hdfs which is a distributed file system. GIS tools for Hadoop on GitHub include geoprocessing tool and geometry API. Hadoop features distributed processing that is resilient to failure.
10. Ruby on rails is being used in the PDX ESRI shop. It seems to be an alternative to JS. Many GitHub projects for Ruby including JSONatra and Angelo. A cool geo-trigger app was demonstrated.

Web mapping nerds can be penned with glass walls!

Some general impressions and things I learned:

1. About 80% of the people there were developing for web mapping applications. There was very little about desktop applications.
2. The flex and silverlight booths were small and lonely.
3. Most developers prefer jQuery to DOJO, but learn DOJO because they want their apps to work with ESRI services.
4. This dev summit marked a decided opening of ESRI products to the open source world. It seems they could not hold back the tide any longer. Most everything was focused on JS and JS frameworks. Data is being passed with geoJSON. 
5. Most folks seemed to be using the chrome debugging environment to work with Javascript. Firebug works too, but I did not see anyone use it. 
6. There is an ArcGIS editor for Open Street Map on GitHub. 
7. ESRI talks to developers on twitter. Everyone has a twitter handle. 
8. A good book to have for Javascript is called "Javascript, the good parts". 
9. Web developers love to play dodge ball. I signed up too late. There were 32 teams with names like "SQL Servers" and "Fork You!"