Module 02

Essential geoprocessing and analysis

Yesterday in the Great Smokies
Webcam archive

GEO 409: Advanced GIS

Announcements

Lab is now open after-hours.
Use student ID for swipe access.

Virtual field trip to

Natural Bridge State Resort Park

37.773072,-83.685751

Historic images from ExploreUKy.edu

Goals

Review first lab
Build a Notebook that will clip any number of layers to any AOI.
Perform summary statistics on the clipped layers
Start making a base map for Kentucky (to use in module 4).

Let's make a little ❤️

Vertices

make points, lines, and polygons.
The basis of the vector data model.
Programatically create geometries with this notebook.

Lab review

Lab 1

Find points within search area.
Now we have a density measure.
- 3.14 * 5² = 78.5 sq mi
- 2 places per 10 sq mi.
How does this compare to other locations?
Trivial question because we have a script!

Imagine...

if you had a hundred locations to search and compare.
Wouldn't it be awesome to automate this redundant task?
Will discover the for loop in this lesson.

Hold up!

We found 15 and 16 places.
How do you represent a polygon as a point?

Center of polygon?

Irregular and multipart polygons?
Centroid?
Why not buffer the polygon?

Lab 2

You have a geodatabase with many layers
that need extraction to an AOI.
Perform summary statistics on the layers.
Automate workflow in Notebook.

Go from massive amounts of data

to analysis and visualization.

Python Essentials

In-class task

Create a Notebook called, 'practice-02.ipynb'
Practice along in the following examples.

What you need to know for ArcPy?

Functions
Lists
for loops

Functions

Block of statements that execute when called.
The function might accept input data as arguments.
The function might return data.
Many ArcPy functions output data to the disk.
📚 ➡ 🔨 ➡ 🌟

Most important function?

The print() function
reveals the contents of variables
and tells you if a block of code runs.
Tip: whenever code is indented, print something!

## Syntax
					---
					Void function example
					```python
					# Assign data to a variable
					x = f"Know what's cool? F-strings can do math: {2 + 2} is 4."
					# Call the print() function and pass x as an argument
					y = print(x)
					# outputs: Know what's cool? F-strings can do math: 4 is 4.
					# when print(x) is invoked (called) it prints the output.
					print(y) # outputs: None
					```

---
					Fruitful function example
					```py
					# Import the random module
					import random

# Create a random integer between 1 and 10
					x = random.randint(1, 10)

# print the random integer
					print(x)
					```
					---
					defining functions
					```python 
					def function_name(parameters):
						# Indented block of statements
						"""docstring: this prints when you call help(function_name)"""
						statement(s) using the parameters
					```

---
					```python [1-5|7-8|10-12]
					# define function
					def squareMe(x):
						"""This function squares any number x"""
						print(x**2)
						return x**2

# use function
					squareMe(5) # outputs 25

# store output to variable
					y = squareMe(5)
					print(y) # outputs 25
					```

Parameters

Variables inside function
which allow us to pass information to the function.
Their sequence is critical.

Clip function

Perhaps most used GIS function
Clips spatial data to an Area of Interest (AOI)
Requires three parameters
Most everyone, initially, gets the parameter order mixed up.

Clip function

The spatial data is the cookie dough
The AOI is the cookie cutter
The output is the cookie

ArcGIS Pro example

Documentation

ALWAYS read the documentation.
Do a web search with 'arcpy {function name}'

F-Strings

String formatting that allows concatenation of:

variables
expressions
and numbers.

Build strings from variables.

Lists

Sequence of values, aka items or elements
Also known as arrays.

Enclosed in square brackets [ ]
Elements separated by commas ,
Elements are ordered.

Common way to store data.

for loops

iterate over a sequence of values.
Use with a list
to access each element in the list.

##### Challenge for 2.5 points
					Print the list of elevations and convert them from meters to feet. 
					List is in the next slide. 
					
					---

```py
					elevations = [132, 152, 194, 59, 72, 35, 93, 206]
					```
					---

##### Challenge
					Find the average elevation in a list of elevations.

---
					
					```py
					# Using a for loop
					# and the plus-equal operator
					x = 0
					x += 1 # Alternative to x = x + 1
					print(x) # prints 1
					```

---
					```py
					x = 0
					n = 0
					for e in elevations:
						x += e
						n += 1
						print(x, n, e)
					avg = x / n
					print(f"Of {n} elevations, the average elevation is {avg}.")
					```

---
					Bonus: find highest and lowest elevations.
					```py
					# Find the number of elements
					len = len(elevations)
					
					# Sort the list
					elevations.sort()
					# sorts highest to lowest on original list

# Print the first and last elements
					print (elevations[0]) # lowest elevation
					print (elevations[len-1]) # highest elevation
					```

Toward completing the lab

In-class task

Create a Notebook called, 'state-base-map.ipynb'
Open lesson documentation and build Notebook cells.
Finish notebook using the section called "Refactored Notebook".

### Setting up the Notebook
						Clip all layers in one database to an area of interest (AOI) and save to a new database.
						
						---
	
						```py
						import arcpy
						```
						---
	
						
						### Local Environment
						Set up our local paths and desired AOI
	
						---
						
						Tune input and output paths to your computer
						```py
						# Input database
						inputDB = 'c:\\BoydsGIS\\data\\input.gdb'
						
						# Create string representing output locations
						myData = 'c:\\BoydsGIS\\data'
						myDB = 'output.gdb'
						outputDB = f'{myData}\\{myDB}'
						
						# Verify
						print(f'Input: {inputDB} Output: {outputDB}')
						
						# AOI in input database
						aoi = 'state_polygon'
						```

---

### 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.")
						```

---
						
						### Analysis
						Clip and save new layers
						
						---
						```py
						# Clip vector layers
						# Create a list of feature class names
						features = arcpy.ListFeatureClasses()
						# Loop over list of features
						for feature in features:
							# Make output name
							output = f"{outputDB}\\{feature}"
							# If aoi, then don't clip; copy
							if feature == aoi:
								arcpy.management.Copy(feature, output)
								print(f"Copying {feature} to {output}")
							else:
								print(f"Clipping {feature} with {aoi} and outputting layer to {output} ...")
								arcpy.analysis.Clip(feature, aoi, output)
						
						# Clip raster layers
						# Create a list of raster names
						rasters = arcpy.ListRasters()
						# Loop over list of rasters
						for raster in rasters:
							output = f"{outputDB}\\{raster}"
							print(f"Clipping {raster} with {aoi} and outputting layer to {output} ...")
							arcpy.management.Clip(raster, "#", output, aoi, "#", "ClippingGeometry", "#")
						```

Notebook output

Layers clipped to AOI (Kentucky)
ready for cartographic excellence
and statistical analysis.

Analyzing attributes

Summary Statistics

Aggregate features based on shared attributes.
How long is the Kentucky River?
How much forested land in Kentucky?
How much of x is in y?

Vector layers

Look for common attributes.
Combine features
and look at Shape_Length and Shape_Area attributes.

Functions

Summary Statistics (Analysis)
- Works on attribute table
Dissolve (Data Management)
- Works on attribute table and geometry
Both functions output a new layer.

In-class task

Explore streams layer
both geometry and attribute table
and use these two functions.

### Summary Statistics
					Function outputs new attribute table with desired summary statistics.
					
					---
					First, switch to output database
					```py
					arcpy.env.workspace = outputDB
					```
					---

### Run function
					What is length of the Kentucky River?

---
					
					```py
					# Parameters: input, output, statistics fields, case field
					arcpy.analysis.Statistics("streams", "streams_stats", [['SHAPE_LENGTH', "Sum"]], "GNIS_NAME")
					# Pseudocode: Sum the length of streams by GNIS_NAME
					```
					```

---

### Dissolve function (default)
					Create new geometry based on shared attributes.
					
				
					---
					```py
					# Parameters: input, output, case field
					arcpy.management.Dissolve('streams', 'streams_joined', 'GNIS_NAME')
					# Pseudocode: Join all lines that share the same GNIS_NAME
					```

---
					
					### Dissolve function (modified)
					Create new geometry based on shared attributes and connectedness.

---
					```py
					# Parameters: input, output, case field, skip option, skip option, unsplit lines
					arcpy.management.Dissolve('streams', 'streams_unique', 'GNIS_NAME', '#', "#", 'UNSPLIT_LINES')
					# Pseudocode: Join all lines that share the same GNIS_NAME and are connected
					```

Statistics on-the-fly

SearchCursor function
iterates over records in a table
accessing any field value.

Raster layers

Matrix of cell values
always a number.
Look for cell size.

Continuous raster

Cell number represent varying quantity.
Elevation, temperature, precipitation, etc.
Digital elevation model (DEM).

Imaging raster

Cell number represent reflected/emitted light.
Multiple bands.
Photography, satellite imagery, renderings, etc.
Shaded relief (artificial illumination).

We'll look at these in the next module

Integer raster

Cell number represent type.
Discrete categories.
Has raster attribute table (RAT)
National land cover dataset (NLCD).

RAT

Count and value fields
with total count * area of cell = total area

In-class task

Explore NLCD layer
and find total area of forest.

Lab 2: Land use in the Bluegrass

Lab 2 Goals

Find percent of total and area of NLCD classes
in a custom area of interest.

Lab 2 Requirements

Check Deliverables.
Basic flow:
- Set up local environment
- Define area of interest
- Clip
- Summary Statistics

Let's jump in

Tip: work with your neighbor.
Code is very sharable.

Custom area of interest?

Use lab 1 to create circular AOI
then use FeatureEnvelopeToPolygon function, convert this circular polygon to square polygon.
Check the Parameters > Python tab in the documentation.

Let's look at a sample Notebook

showing a completed lab.

### Set up local environment
					Set the local file paths and environment variables.
					
					---
					```py
					import arcpy
					```

---
					```py
					# Input database
					inputDB = 'c:\\BoydsGIS\\data\\ky_small_scale.gdb'
					
					# Create string representing output locations
					myData = 'c:\\BoydsGIS\\data'
					myDB = 'bluegrass-crops.gdb'
					outputDB = f'{myData}\\{myDB}'
					
					# Verify
					print(f'Input: {inputDB} Output: {outputDB}')
					```

---
					```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
					```

---
					```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.")
					```
					---
					### Custom AOI
					Find center of AOI and create a point geometry.
				
					---
					```py
					# CRS for WGS84, the standard for GPS
					wgs84 = arcpy.SpatialReference(4326)
					# set the location of AOI's center from GPS coordinates
					lon = -83.682834
					lat = 37.783013
					# create a single point using GPS coordinates
					point = arcpy.Point(lon, lat)
					# create the point geometry using the WGS84 CRS
					ptGeometry = arcpy.PointGeometry(point, wgs84)
					# show point in the Notebook
					ptGeometry
					# save the point to a feature class in our workspace
					arcpy.management.CopyFeatures(ptGeometry, "aoi_center")
					# input is ptGeometry, output is a point layer in database
					```
					---
					### Custom AOI
					Buffer and convert to square.
					
					---
					```py
					# create a one-mile buffer around the point
					arcpy.analysis.Buffer() # Parameters?
					# convert to square
					arcpy.management.FeatureEnvelopeToPolygon() # Parameters?
					```
					---
					### Analysis
					Set new AOI and clip layers to square AOI.
					
					---
					```py
					# Declare aoi
					aoi = 'aoi'
					```

---
					```py
					# Clip vector layers
					# Create a list of feature class names
					features = arcpy.ListFeatureClasses()
					# Loop over list of features
					for feature in features:
						# Make output name
						output = f"{outputDB}\\{feature}"
						# If aoi, then don't clip; copy
						if feature == aoi:
							arcpy.management.Copy(feature, output)
							print(f"Copying {feature} to {output}")
						else:
							print(f"Clipping {feature} with {aoi} and outputting layer to {output} ...")
							arcpy.analysis.Clip(feature, aoi, output)
					
					# Clip raster layers
					# Create a list of raster names
					rasters = arcpy.ListRasters()
					# Loop over list of rasters
					for raster in rasters:
						output = f"{outputDB}\\{raster}"
						print(f"Clipping {raster} with {aoi} and outputting layer to {output} ...")
						arcpy.management.Clip(raster, "#", output, aoi, "#", "ClippingGeometry", "#")
					```

---
					### Summary Statistics
					Find total count of cells and area of each class.

---
					```py
					# Switch to output database
					arcpy.env.workspace = outputDB
					```
					
					---
					```py
					# Find total count of all cells
					total = 0
					with arcpy.da.SearchCursor('NLCD_2016_land_cover', ['Count']) as cursor:
						for row in cursor:
							print(row)
							total += row[0]
					# Create area measurement
					sqMi = (total * 98.42**2) / 2.788e+7
					sqMiRounded = round(sqMi, 1)
					print(f'{sqMiRounded} square miles')
					```
					---
					```py
					# Find percent and area for each land cover class
					with arcpy.da.SearchCursor('NLCD_2016_land_cover', ['Count', 'NLCD_Land_Cover_Class']) as cursor:
						for row in cursor:
							# find percent
							a = (row[0]/total) * 100
							# find area in sq miles
							b = (row[0] * 98.4**2) / 2.788e+7
							# round and redeclare a and b
							a = round(a, 2)
							b = round(b, 1)
							c = f'{a}%, {b} sq mi of {row[1]}'
							print(c)
					```

Gotchas

Input/output

Where's my data? What database are you using?
Functions have input AND output.

Errors?

If you are not getting them, then you are
- lucky
- or not practicing. 🤦
Errors are your friend.

Syntax errors

Typeos
in file path names, variables, etc.
Using functions properly? Parameters correct?
Throws an error, look at last line.

Solutions

If the error has a code, ERROR 000670, search it.
Remove layers from Contents pane.
Rebuild your Notebook from scratch.

Big Picture

Once you have a working Notebook
you can answer how much of x is in y *
for any area in Kentucky.

Review

Python II

Video lessons
- functions watch 🎥
- strings watch 🎥
- for and while loops watch 🎥

Challenge

Build a function

Store and reuse statements
def keyword defines a function
followed by unique name and():
Parameters? i.e., does the function use data?
Indent function body four spaces

```python
					# Make a function that takes one string argument
					greeting("Yoda")
					# and outputs: "Hello, Yoda!"
					```
					---

```python
					def greeting(name):
						print(f"Hello, {name}!")
					n ="Yoda"
					greeting(n)
					```

---
					```python
					def greeting(name):
						print(f"Hello, {name}!")
					n ="Yoda"
					greeting("Luke")
					# prints "Hello, Luke!" because
					# parameter is used in function definition.
					```
					---

### Convert function to a value

`return` statement returns a value from a function

---

```python
					def greeting(name):
						print(f"Hello, {name}!")
					x = greeting("Yoda")
					print(x) # returns None
					```
					---

```python
					def greeting(name):
						return f"Hello, {name}!"
					x = greeting("Yoda")
					print(x) # returns "Hello, Yoda!"
					```
					---

```python
					# What is this function doing?
					function(2, 3)
					# prints 8
					```
					---
					```python
					# define function
					def powersXtoY(x, y): # function has arbitrary name
						z = x**y
						print(z) # prints to screen

powersXtoY(2,3) # Need documentation!
					```
					---
					```python
					# define function
					def powersXtoY(x, y): # function has arbitrary name
						"""x (1 arg) to y (2 arg) power""" # Docstring provides info
						z = x**y
						print(z) # prints to screen

powersXtoY(2,3) # outputs 8

# Get docstring
					help(powersXtoY) # outputs x (1 arg) to y (2 arg) power
					```

Documentation

Whenever you use a function
look up its documentation.
ArcGIS Pro has excellent docs 🚀

Methods

Functions that only work on certain types. 🔮🐉
String type, integer type, csv type, etc.
Uses dot notation, i.e., myString.upper()

Strings

String Methods

Perform common tasks on strings
Find substrings, change case, etc.
Most return new values, others return `True` or `False`
They do not change the original string.
Built-in methods

Example

justWords = "If you live to be a hundred, I want to live to be a hundred minus one day so I never have to live without you."

String slice

Access substring using [start : end] indexes
Two index values separated by a colon : inside bracket [] notation
End index is up-to-but-not-including
"Let's Hike"[0:2] returns "Le"

String stride

Optional third number in [start:end:stride]
Third number returns every nth character
By default, 1
Negative numbers reverse the string

Challenge

Slice a URL

A uniform resource locator is the unique address of anything on the Web.

Loops

Loop

Interrupts sequential execution of code
runs same block of code over and over
modifying it with each new run
until a desired condition is found.

Syntax

Loop statement ends with semicolon :
loop block is indented four spaces
runs block until a condition is met.

`for` loop

Iterate over sequence
value by value
with an iterating variable
until all values have been accessed
or until a break occurs in the block.

`while` loop

Runs while statement remains true
or until a break occurs in the block.
Can run indefinite number of times.

Python lists

A collection of values in a sequence.
One of the most common data collections.
Values enclosed in [ ] square brackets
and separated by commas.

```python
					x = [3123, 5342, 2432, 5433, 1933, 4322]

# loop statement
					for i in x:
						print(i)
					```
					---
					```python
					x = [3123, 5342, 2432, 5433, 1933, 4322]

# loop statement
					for i in x:
						if i > 5280:
							print("We went a mile high!")
							break
					```
					---
					```python
					x = [3123, 5342, 2432, 5433, 1933, 4322]
					lenX = len(x)
					y = 0
					z = 0

# loop statement
					for i in x:
						if i > y:
							y = i
						z += i
					print(f"The average elevation is {z/lenX} and the highest is {y}.")
					```
					---
					```python
					x = [
					"ftp://ftp.kymartian.ky.gov/kyaped/DEMs_5FT/N091E298_DEM.zip",
					"ftp://ftp.kymartian.ky.gov/kyaped/DEMs_5FT/N091E299_DEM.zip",
					"ftp://ftp.kymartian.ky.gov/kyaped/DEMs_5FT/N091E300_DEM.zip",
					"ftp://ftp.kymartian.ky.gov/kyaped/DEMs_5FT/N091E301_DEM.zip",
					"ftp://ftp.kymartian.ky.gov/kyaped/DEMs_5FT/N091E302_DEM.zip"]

for i in x:
						y = i[43:-4]
						print (y)
					```
					---
					##### Challenge
					Build a function that finds the highest elevation in a list of elevations.

---
					```py
					elevations = [1232, 452, 1224, 599, 745, 355, 899]
					```

---
					```py
					def findStats(dataList):
						x = 0
						for number in dataList:
							if number > x:
								x = number
						print(x)
					```

---
					### `return` value

* To access value we need to return it in the function definition

---
					```py
					def findStats(dataList):
						"""Find the highest number in a list of numbers"""
						x = 0
						for number in dataList:
							if number > x:
								x = number
						print(f"{x} is the highest elevation.")
						return x
					```

---
					```py
					# Use function on this data 
					elevations = [1232, 452, 1224, 599, 745, 355, 899]
					highest = findStats(elevations) # outputs 1232 is the highest elevation.

# Print variable contents
					print(highest) # outputs 1232
					```

### Built-in functions

* Python comes with about [60 functions](https://docs.python.org/3.8/library/functions.html)
					* Don't need to `import` anything
					
					
					---
					### Round function
					
					```python
					round(x, y)
					# x is float type, y is integer type
					print(round(7.555, 2))
					```
					
					---
					# What!?!
					
					---
					### Adding functions
					
					* Python's sprawling [standard library](https://docs.python.org/3/library/) of modules
					* Adds new functions and data types
					* The `decimal` module to correctly round numbers
					* Use dot notation to address functions
					
					
					---
					```python
					import decimal # use the module's namespace to access its functions
					x = decimal.Decimal('7.555')
					
					from decimal import Decimal # access the module functions directly
					x = Decimal('7.555')
					
					print(round(x, 2))
					```

---
					### decimal properties
					
					* Control precision and number of significant digits
					* Rounding rules
					* Useful in mapping and money applications
					
					
					---
					### Lat/lon coords and precision
					
					* 38.038015, -84.5046852!
					* 1° of latitude = 69 miles
					* 1.0°  6.9 mi of precision
					* 1.00°  3600 ft
					* 1.000°  360 ft
					* 1.0000°  36 ft
					* 1.000000°  4 in

Module 02

Essential geoprocessing and analysis

Announcements

Virtual field trip to

Natural Bridge State Resort Park

Goals

Let's make a little ❤️

Vertices

Lab review

Lab 1

Imagine...

Hold up!

Center of polygon?

Lab 2

Go from massive amounts of data

to analysis and visualization.

Python Essentials

In-class task

What you need to know for ArcPy?

Functions

Most important function?

Parameters

Clip function

Clip function

ArcGIS Pro example

Documentation

F-Strings

Lists

for loops

Toward completing the lab

In-class task

Notebook output

Analyzing attributes

Summary Statistics

Vector layers

Functions

In-class task

Statistics on-the-fly

Raster layers

Continuous raster

Imaging raster

We'll look at these in the next module

Integer raster

RAT

In-class task

Lab 2: Land use in the Bluegrass

Lab 2 Goals

Lab 2 Requirements

Let's jump in

Custom area of interest?

Let's look at a sample Notebook

Gotchas

Input/output

Errors?

Syntax errors

Solutions

Big Picture

Review

Python II

Challenge

Build a function

Build a function

Documentation

Methods

Strings

String Methods

Example

String slice

String stride

Challenge

Slice a URL

Loops

Loop

Syntax

for loop

while loop

Python lists

STOP

`for` loop

`while` loop