Module 03

Geoprocessing II: spatial modeling

Yesterday in the Great Smokies
Webcam archive

GEO 409: Advanced GIS

Announcements

Please see me in office hours
or schedule a time to Zoom
if Notebooks still seem completely alien.

Virtual field trip to

Louisville Waterfront

38.263568,-85.741381

Historic images from ExploreUKy.edu

Lab 2 review

Your boss says, "Shearer, get in here!"
"This looks like a long 20 miles."
You think, "Uh-oh, I know I used 20 miles."

This looks fishy!

0.66%, 10.5 sq mi of Open Water
5.96%, 95.3 sq mi of Developed, Open Space
3.76%, 60.2 sq mi of Developed, Low Intensity
2.09%, 33.5 sq mi of Developed, Medium Intensity
0.85%, 13.7 sq mi of Developed, High Intensity
0.26%, 4.1 sq mi of Barren Land
14.73%, 235.6 sq mi of Deciduous Forest
0.73%, 11.6 sq mi of Evergreen Forest
4.65%, 74.4 sq mi of Mixed Forest
0.13%, 2.1 sq mi of Shrub/Scrub
0.46%, 7.3 sq mi of Herbaceuous
61.94%, 990.7 sq mi of Hay/Pasture
3.7%, 59.1 sq mi of Cultivated Crops
0.04%, 0.7 sq mi of Woody Wetlands
0.03%, 0.5 sq mi of Emergent Herbaceuous Wetlands

Know your tools

Read the docs for functions you use.
At the very least, know the required parameters.

Solution?

If you Notebook already works:
- Delete output geodatabase
- Apply the fix in the Notebook
- Rerun entire notebook
Submit revised Notebook by Thursday, end-of-day.

Goals

Fine tune our AOIs.
Use spatial modelling to predict unobserved conditions.
Use lidar derivatives 🏋️‍♀️.

Review Geoprocessing

Geoprocessing

"A framework and set of tools for processing geographic and related data."
Framework: workspace, environment, sequence of functions.
Tools: functions like buffer, clip, etc.
Input => Output

Automated

Steps are scripted in the Notebook.
Workspace parameters are exposed at the head of the script to control downstream cells.

Automate at scale

Programmatically run an entire workflow
on unlimited amounts of data.
Once coded, work on same data but from different area.

Previous labs

were totally input => output.
This lab will be more exploratory.
Get comfortable changing function parameters, rerunning cell, examining output.

Custom AOIs

Essential first step

to define the area of interest.
Can we use counties, river corridors, whatever we can draw on a map?
Yes!

Refactor point AOI

into a single function.
Use it anywhere to get a buffer around a point.

### Imagine   
					if you could just do this

```py
					# Create AOI
					makeAoi(lat, lon, dist, output, square=True)
					```

---
					### First, set up environment 
					---

### Load ArcPy site package
					```py
					import arcpy
					```
					
					---
					### Local Environment
					Set up our local paths
					
					---
					
					Tune input and output paths to your computer
					```py
					# Input database
					inputDB = 'c:\\BoydsGIS\\data\\ky_small_scale.gdb'
					
					# Create string representing output locations
					myData = 'c:\\BoydsGIS\\data'
					myDB = 'practice-03.gdb'
					outputDB = f'{myData}\\{myDB}'
					
					# Verify
					print(f'Input: {inputDB} Output: {outputDB}')
					
					```
					
					---
					
					### Local Environment
					Set workspace, projections, and overwrite output
					
					
					
					---
					```py
					# Input database
					arcpy.env.workspace = inputDB
					# Define output CRS
					ky = arcpy.SpatialReference(3089)
					# Set the output CRS
					arcpy.env.outputCoordinateSystem = ky
					# Enable overwriting exist layers
					arcpy.env.overwriteOutput = True
					```
					
					---
					
					### Local Environment
					Create the output database
					
					---
					```py
					# Creating output database
					print("Creating workspace geodatabase...")
					# Check if it already exists
					if arcpy.Exists(outputDB):
						print(f"Output database exists: {outputDB}")
					else:
						arcpy.management.CreateFileGDB(myData, myDB)
						print(f"Created output database: {outputDB}")
					
					print(f"Using {arcpy.env.workspace} for input data.")
					```
					
					---
					### Define new function
					to create circular or square AOI

---
					Function definition
					```py
					def makeAoi(lat, lon, dist, output, square=True):
						"""Function to create round or square area of interest for lat, lon coordinate
						Parameters:
							lat = latitude (Number)
							lon = longitude (Number)
							dist = buffer distance (Number or valid string for Buffer function)
							output = workspace polygon layer name (String)
							square {optional} = square (default) or circular aoi polygon
						"""
						out1 = f'{output}_pt_temp'
						out2 = f'{output}_ro_temp'
						arcpy.env.addOutputsToMap = False
						wgs84 = arcpy.SpatialReference(4326)
						point = arcpy.Point(lon, lat) 
						ptGeometry = arcpy.PointGeometry(point, wgs84)
						arcpy.management.CopyFeatures(ptGeometry, out1)
						if square:
							arcpy.analysis.Buffer(out1, out2, dist)
							arcpy.env.addOutputsToMap = True
							arcpy.management.FeatureEnvelopeToPolygon(out2, output)
							arcpy.management.Delete(out2)
						else:
							arcpy.env.addOutputsToMap = True
							arcpy.analysis.Buffer(out1, output, dist)
						arcpy.management.Delete(out1)
						print(f"Created layer: {output}") 
					```

---
					Then, run the function
					```py
					# Define the variables first
					lat = 38.038125
					lon = -84.503860
					dist = '1 Mile' # 1 mile buffer
					output = 'whitehall_2mi' # What's the workspace?
					# Call the function
					makeAoi(lat, lon, dist, output)
					```

---
					Then, run the function using output database
					```py
					# Define the variables first
					lat = 38.038125
					lon = -84.503860
					dist = '1 Mile' # 1 mile buffer
					output = f'{outputDB}\\whitehall_2mi' # What's the workspace?
					# Call the function
					makeAoi(lat, lon, dist, output)
					```

In-class task

Create new 'module-03' project and 'lesson-03.ipynb' Notebook.
Open lesson docs to 'custom areas of interest'
and create the function and make an AOI.

AOIs from existing features

Use a county or river corridor as an AOI.
Select feature
and use CopyFeatures() and/or Buffer() functions.

Arbitrary AOIs

Any shape you can draw on a map.
Use use CreateFeatureclass() function to create empty layer.
Use Edit ribbon to create shape and save it.

In-class task

Create three new AOIs.
- A county
- A river
- An arbitrary shape

Spatial Modeling

Spatial data encodes observations and models where conditions exist.
Can accurate/precise observed data predict unobserved conditions?
Roads are now mostly GPS'd. Can their density predict population?

Spatial Modeling

Where x_n in y, z exists

Spatial Analyst tools

arcpy.sa module
Leverage Python and map algebra to facilitate modeling and analysis.
Note: if you have a license from another department, you might need an update.

Map algebra

Uses arithmetic operators (+, - <,>, etc.)
on one or many raster layers
to locate where a condition exists.
Produces temporary raster.

Map algebra

Enclose all queries with ()
- z = (elev > 1000)
Logical operators and and or
- & (and)
- | (or)
- z = ((elev < 1000) & (elev > 900))

Raster Review

Resolutions

Spatial: how large is each cell?
Temporal: how often is the same scene captured?
Bit depth: 2ⁿ number of values the cell can store.
Spectral: how much of the spectrum, e.g. light, is captured?

Matrix of cells

of numeric values.
Is it a quantity
or an category?

Continuous

Measures a varying phenomena
like elevation, slope, etc.
in a single band.
Cell values are a quantity.

Discrete

Measures distinct types
like land cover type
in a single band.
Cell values are a category.

Image

Measures reflected or emitted light
like a photograph
in multiple bands.
Cell values are a relative intensity.

Identity

identifies whether a condition exists
like elevations above 1000 feet
in single band.
Cell values are either 0 or 1.

In-class task

Find the area of elevations over 2,000 feet in Ky.
Find the highest and lowest elevations in Ky.

Where are locations over 2,000 feet in elevation?
				```py
				# Create raster object
				elev = arcpy.sa.Raster('DEM_120_M')

# Find elevations over 2,000 feet
				# Group all queries in ()
				z = (elev > 2000)
				```
				---
				```py
				z.save('elevOver2000')
				# Where is this saved?
				```
				---
				Highest and lowest locations
				```py
				# Get statistics for elevation layer
				elev.getStatistics()
				```
				---
				```py
				# Get selected statistics for elevation layer
				print(f'Minimum elevation: {elev.minimum}')
				print(f'Maximum elevation: {elev.maximum}')
				```
				---
				```py
				# Extract locations with the minimum elevation
				min_elev = arcpy.sa.ExtractByAttributes(elev, f'VALUE = {elev.minimum}')
				arcpy.conversion.RasterToPoint(min_elev, f'{outputDB}/min_elev')
				
				# Extract locations with the maximum elevation
				max_elev = arcpy.sa.ExtractByAttributes(elev, f'VALUE = {elev.maximum}')
				arcpy.conversion.RasterToPoint(max_elev, f'{outputDB}/max_elev')
				```

Multiple layer queries

Surface Parameters

Spatial Analyst function to model elevation surface.
Slope: rate of elevation change in steepest direction. Measured in degrees.
Aspect: direction of steepest slope. Measured in degrees 0-360 with 0 being north.

In-class task

Find the north-facing steep slopes in Ky.
Slope greater than 16 degrees.

KyFromAbove

Authoritative source of high-resolution aerial imagery and elevation data for Kentucky.
Download individual tiles or load as an ArcGIS service.

Use as an ArcGIS service

Add a new ArcGIS Server and use this URL:
https://kyraster.ky.gov/arcgis/services

Data

DSM (Digital Surface Model) shows the first return of the laser pulse: elevation of above ground features.
DEM (Digital Elevation Model) shows hydrologically enforced (e.g., bridges removed) bare-earth elevations.
NAIP (National Agriculture Imagery Program) leaf-on 4-band imagery.

Height raster

Subtracting the DEM from the DSM gives us a height raster for above ground features
such as buildings, bridges, and trees.
How high is the bridge connecting Whitehall Classroom Building with POT?

Tree canopy

Where are the trees?
Model tree location and height.

Where are green things?

A derivative of NAIP imagery assesses the health of vegetation.
Uses 4-band imagery to calculate a vegetation index.
A higher value indicates more vegetation.

NDVI

Normalized difference vegetation index (NDVI) compares reflected red to near-infrared light (which we cannot see with our eyes but can with optical sensors)
estimate the health and density of green vegetation.

Inspection

After creating the NDVI raster, which value is the minimum threshold for vegetation?

In-class task

Find cells that have vegetation and are over 5 feet tall.

Location?

We might have location but not height.
Use Con function to create a new raster.
Where this query is true, use the cell value from this raster.

Analysis

Group height values into discrete classes.
5-15 feet, 15-35 feet, etc.
Count the number of pixels in each class.

Reclassify
					```py
					# Reclassify old values to new range
					# format [old_lo, old_high, new_value], next new value, etc. 
					reramp = arcpy.sa.RemapRange([[5,15,1],[15,35,2],[35,75,3],[75,140,4],[140,300,5]])
					trees_classed = arcpy.sa.Reclassify(trees, "Value", reramp, "NODATA")
					trees_classed.save('trees_classed')
					```
					---
					Create fancy output
					```py
					# Lookup table for reclassified values
					lookup = {1: '5-15 ft', 2: '15-35 ft', 3: '35-75 ft', 4: '75-140 ft', 5: '140-300 ft (potential error)'}
					```
					---
					Summarize
					```py
					# Get stats for each new class
					resolution = 5
					sumCover = 0
					n = 1
					with arcpy.da.SearchCursor('trees_classed' , ['Count']) as cursor:
						for row in cursor:
							sumCover += row[0]
					with arcpy.da.SearchCursor('trees_classed', ['Count', 'Value']) as cursor:
						for row in cursor:
							a = (row[0]/sumCover) * 100
							b = (row[0] * resolution**2) / 43560
							a = round(a, 2)
							b = round(b, 1)
							c = f'{a}%, {b} acres of {lookup[n]}'
							n += 1
							print(c)
					```

Lab 3

Lab

Find a desired condition in your AOI.
Notebook should produce an analysis layer.
This 'recipe' should run on any ArcGIS Pro instance with minimal modification.

Exploration

The conditions you set require some exploration of your data.
What are the ranges of values in your raster?
Do they correspond to your on-the-ground observations?
No right answer – just how well does your analysis match the desired condition.

In-class task

Visualize data in new Local Scene
add DSM as ground elevation layer
symbolize NDVI to visualize vegetation
and explore 3D map.

Addendum

Mini apps

Workflow

Self-contained Notebooks
- using data from KyFromAbove.
Run top to bottom.
Rerun cells testing different arguments.
Save final layers to geodatabase.

Viz

Visualize your output
- over different base maps.
Explore!
We'll do cartography in the next module.

Contouring

Notebook: create-smooth-contours.ipynb
Lines of equal elevation interpolated from an elevation surface.
Find cliff line and estimate height (if cliffs are present).

Labels and symbology

Symbolize lines with index field
Give labels a halo effect that matches the color of the base map.
ArcGIS Pro labeling basics

Solar mapping

Notebook: solar-mapping.ipynb
Estimate relative rate of direct insolation

Can I see sunshine?

Use the Hillshade function with shadows enabled
to estimate the amount of direct sunlight
based on slope of surface.

Where is the sun?

Use NOAA's solar calculator
to find the sun's position
on any given day and time.

Symbology

Symbolize output layer over NDVI
with a high values bright yellow
and use Layer blending mode: Screen.

Bonus

Create animation of sun's illumination over time
with open-source command line tool ffmpeg.
Example in addendum.

Review

Python III

Videos
- lists watch 🎥
- dictionaries watch 🎥
- tuples watch 🎥

Python collections

Collections

of data.
Ordered vs. unordered
mutable vs. immutable
Nested collections (collection of collections)
Looping through collections

Lists

An mutable, sequence of items or elements
separated by a comma
and enclosed in square brackets.
Items can be any value.

```python
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch'] 
					```

---
					### Index bracket notation
					Address items in sequence with integer, just like strings.

---
					```py
					# where item is an integer
					item = 1
					archesList[item] # returns 'Castle Arch'
					```

---
					## Challenge: 
					### Loop over archesList
					```py
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					```

---
					```python
					# archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch'] 
					n = 0
					for i in archesList:
						print(i)
						print(archesList[n])
						n += 1
					```

---
					### Add items
					```python
					archesList[3]= "Rainbow Arch"
					# ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch', "Rainbow Arch"] 
					```

---
					### Change items
					```python
					archesList[3]= "Sky Bridge"
					# ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch', "Sky Bridge"] 
					```

---
					### Nested collections
					Because lists can contain any type, you can have lists inside lists. Gang together brackets to access nested collection.

---
					```py
					myList[outerItem][innerItem]
					# where outerItem and innerItem are integers
					```

---
					```py
					archCoords = [[37.7873, -83.6243], [37.8123, -83.5796], [37.7686, -83.5582]]
					```

---
					```python
					# archCoords = [[37.7873, -83.6243], [37.8123, -83.5796], [37.7686, -83.5582]]
					n = 0
					for i in archCoords:
						print(i)
						print(archesList[n])
						print(i[0])
						print(archesList[n][1])
						n += 1
					```

---
					### List methods
					Changes original list!
					```py
					.sort() # alphabetically
					.append() # add one element
					.extend() # add elements from one iterable
					.remove() # remove first occurrence
					```

---
					### `.sort()`
					Alphabetize elements in list
					```py
					# archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archesList.sort() 
					# ['Castle Arch', 'Rock Bridge Arch', 'Silvermine Arch']
					archesList.sort(reverse=True)  # keyword argument
					```

---
					### `.append()`
					One argument. Adds element to end of list
					```py
					# archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archesList.append('Grays Arch') 
					# ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch', 'Grays Arch']
					```

---
					### `.extend()`
					One argument: an iterable. Adds element(s) to end of list
					```py
					# archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archesList.extend(['Grays Arch', 'Natural Bridge']) 
					# ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch', 'Grays Arch', 'Natural Bridge']
					```

---
					### Useful list functions
					Returns a value
					```py
					len() # number of elements in the list
					range() # range of index values
					```

---
					### `len()`
					Returns number of characters in string or elements in list
					```py
					# archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					x = len(archesList)
					print(x) 
					# 3
					```

---
					### `range()`
					Returns sequence of numbers (over which to iterate)
					```py
					for i in range(6):
						print(i)
					for i in range(len(archesList)):
						print(i, archesList[i])
					```

---
					### multiple iterators
					```
					archCoords = [[37.7873, -83.6243], [37.8123, -83.5796], [37.7686, -83.5582]]
					for x, y in archCoords:
						print(y, x)
					```

---
					## Challenge: 
					### Print a sorted list of arches
					```py
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					```

---
					## Challenge: 
					### combine and sort arch names and locations
					```py
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archCoords = [[37.7873, -83.6243], [37.8123, -83.5796], [37.7686, -83.5582]]
					```

---
					## Challenge: 
					### Sort by east to west
					```py
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archCoords = [[37.7873, -83.6243], [37.8123, -83.5796], [37.7686, -83.5582]]
					```

---
					## Solutions

---
					```py
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archesList.sort()
					for arch in archesList:
						print (arch)
					```

---
					```py
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archCoords = [[37.7873, -83.6243], [37.8123, -83.5796], [37.7686, -83.5582]]

aList = []
					i = 0

for arch in archesList:
						aList.append([arch, archCoords[i]])
						i += 1
					aList.sort()
					print(aList)
					```

---
					### a more Pythonic way
					```py
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archCoords = [[37.7873, -83.6243], [37.8123, -83.5796], [37.7686, -83.5582]]

aList = []
					i = 0

# Using zip() function
					for x, y in zip(archesList, archCoords):
						aList.append([x, y])
						i += 1
					aList.sort()
					print(aList)
					```

---
					```py
					aList = []
					archesList = ['Silvermine Arch', 'Castle Arch', 'Rock Bridge Arch']
					archCoords = [[37.7873, -83.6243], [37.8123, -83.5796], [37.7686, -83.5582]]
					n = 0
					for i in archesList:
						aList.append([[archCoords[n][1], archCoords[n][0]], i])
						n += 1
					aList.sort(reverse=True) # reverse needed?
					print(aList)
					```

Tuples

An immutable, sequence of items or elements
separated by a comma
and enclosed in round brackets.
Items can be any value.

Tuples

Think of tuples like records in a database,
a collection or related information
often of a different type.
Tuples are fast.

Dictionaries

An unordered collection of items
stored as key:value pairs
separated by a colon
and enclosed in curly brackets.

Dictionaries

Access key:value pairs is like looking in a dictionary.
Iterating over...
- keys
- values
- key:value pairs

```python
					archesDictionary = {
						'Castle Arch': [37.8123, -83.5796],
						'Red-byrd Arch': [37.8148, -83.5504],
						'Rock Bridge Arch': [37.7686, -83.5582],
						'Silvermine Arch': [37.7873, -83.6243]
						}
					```
					---
					### Iterate through dictionary
					```python
					for i in archesDictionary:
						print(i, archesDictionary[i])
					```

---
					### Iterate just keys
					```py
					for k in archesDictionary.keys():
						print(k)
					```

---
					### Iterate just values
					```py
					for v in archesDictionary.values():
						print(v)
					```
					---
					### Iterate each item with key:value as tuple
					```py
					for i in archesDictionary.items(): # returns tuple of items!
						print(i)
					```
					---
					## Challenge
					### In two lines, print: 
					'Castle Arch is located at 37.8123 lat and -83.5796 lon' for each arch
					```py
					archesDictionary = {
					 'Castle Arch': [37.8123, -83.5796],
					 'Red-byrd Arch': [37.8148, -83.5504],
					 'Rock Bridge Arch': [37.7686, -83.5582],
					 'Silvermine Arch': [37.7873, -83.6243]
					 }
					```
					---
					## Solution
					---
					```py
					for (x, y) in archesDictionary.items():
						print(f'{x} is located a {y[0]} lat and {y[1]} long.')
					```

Advantages of each?

Lists are easy to build/modify
Tuples are fast to access and secure
Dictionaries are fast, tidy, and easy to build/modify

Module 03

Geoprocessing II: spatial modeling

Announcements

Virtual field trip to

Louisville Waterfront

Lab 2 review

This looks fishy!

Know your tools

Solution?

Goals

Review Geoprocessing

Geoprocessing

Automated

Automate at scale

Previous labs

Custom AOIs

Essential first step

Refactor point AOI

In-class task

AOIs from existing features

Arbitrary AOIs

In-class task

Spatial Modeling

Spatial Modeling

Spatial Modeling

Where xn in y, z exists

Spatial Analyst tools

Map algebra

Map algebra

Raster Review

Resolutions

Matrix of cells

Continuous

Discrete

Image

Identity

In-class task

Multiple layer queries

Surface Parameters

In-class task

KyFromAbove

Use as an ArcGIS service

Data

Height raster

Tree canopy

Where are green things?

NDVI

Inspection

In-class task

Location?

Analysis

Lab 3

Lab

Exploration

In-class task

Addendum

Mini apps

Workflow

Viz

Contouring

Labels and symbology

Solar mapping

Can I see sunshine?

Where is the sun?

Symbology

Bonus

Review

Python III

Python collections

Collections

Lists

Tuples

Tuples

Dictionaries

Dictionaries

Advantages of each?

STOP

Where x_n in y, z exists