The form of the “import” command that you definitely should already know is

import <module name>

e.g.,

import arcpy 

What happens here is that the module (either a module from the standard library, a module that is part of another package you installed, or simply another .py file in your project directory) is loaded, unless it has already been loaded before, and the name of the module becomes part of the namespace of the script that contains the import command. As a result, you can now access all variables, functions, or classes defined in the imported module, by writing

<module name>.<variable or function name>

e.g.,

arcpy.Describe(…) 

You can also use the import command like this instead:

import arcpy as ap 

This form introduces a new alias for the module name, typically to save some typing when the module name is rather long, and instead of writing

arcpy.Describe(…)  

, you would now use

ap.Describe(…) 

in your code.

Another approach of using “import” is to directly add content of a module (again either variables, functions, or classes) to the namespace of the importing Python script. This is done by using the form "from … import …" as in the following example:

from arcpy import Describe, Point , …  

... 

Describe(…) 

The difference is that now you can use the imported names directly in our code without having to use the module name (or an alias) as a prefix as it is done in line 5 of the example code. However, be aware that if you are importing multiple modules, this can easily lead to name conflicts if, for instance, two modules contain functions with the same name. It can also make your code a little more difficult to read since

  arcpy.Describe(...)
  

helps you or another programmer recognize that you’re using something defined in arcpy and not in another library or the main code of your script.

You can also use

from arcpy import *

to import all variable, function and class names from a module into the namespace of your script if you don’t want to list all those you actually need. However, this can increase the likelihood of a name conflict.

Next, let’s quickly revisit loops in Python. There are two kinds of loops in Python, the for-loop and the while-loop. You should know that the for-loop is typically used when the goal is to go through a given set or list of items or do something a certain number of times. In the first case, the for-loop typically looks like this

for item in list: 

    # do something with item 

while in the second case, the for-loop is often used together with the range(…) function to determine how often the loop body should be executed:

for i in range(50):  

	# do something 50 times

In contrast, the while-loop has a condition that is checked before each iteration and if the condition becomes False, the loop is terminated and the code execution continues after the loop body. With this knowledge, it should be pretty clear what the following code example does:

import random 

r = random.randrange(100) # produce random number between 0 and 99 

attempts = 1 

while r != 11: 

    attempts += 1 

    r = random.randrange(100) 

print('This took ' + str(attempts) + ' attempts') 

What you may not yet know is that there are two additional commands, break and continue, that can be used in combination with either a for or a while-loop. The break command will automatically terminate the execution of the current loop and continue with the code after it. If the loop is part of a nested loop only the inner loop will be terminated. This means we can rewrite the program from above using a for-loop rather than a while-loop like this:

import random 

attempts = 0 

for i in range(1000):  

    r = random.randrange(100) 

    attempts += 1 

    if r == 11: 

        break  # terminate loop and continue after it 

print('This took ' + str(attempts) + ' attempts') 

When the random number produced in the loop body is 11, the body of the if-statement, so the break command, will be executed and the program execution immediately leaves the loop and continues with the print statement after it. Obviously, this version is not completely identical to the while based version from above because the loop will be executed at most 1000 times here.

If you have experience with programming languages other than Python, you may know that some languages have a "do … while" loop construct where the condition is only tested after each time the loop body has been executed so that the loop body is always executed at least once. Since we first need to create a random number before the condition can be tested, this example would actually be a little bit shorter and clearer using a do-while loop. Python does not have a do-while loop but it can be simulated using a combination of while and break:

import random

attempts = 0  

while True: 

    r = random.randrange(100) 

    attempts += 1 

    if r == 11: 

        break 

print('This took ' + str(attempts) + ' attempts') 

A while loop with the condition True will in principle run forever. However, since we have the if-statement with the break, the execution will be terminated as soon as the random number generator rolls an 11. While this code is not shorter than the previous while-based version, we are only creating random numbers in one place, so it can be considered a little bit more clear.

When a continue command is encountered within the body of a loop, the current execution of the loop body is also immediately stopped, but in contrast to the break command, the execution then continues with the next iteration of the loop body. Of course, the next iteration is only started if, in the case of a while-loop, the condition is still true, or in the case of a for-loop, there are still remaining items in the list that we are looping through. The following code goes through a list of numbers and prints out only those numbers that are divisible by 3 (without remainder).

l = [3,7,99,54,3,11,123,444] 

for n in l: 

    if n % 3 != 0:   # test whether n is not divisible by 3 without remainder 

        continue 

    print(n)

This code uses the modulo operator % to get the remainder of the division of n and 3 in line 5. If this remainder is not 0, the continue command is executed and, as a result, the program execution directly jumps back to the beginning of the loop and continues with the next number. If the condition is False (meaning the number is divisible by 3), the execution continues as normal after the if-statement and prints out the number. Hopefully, it is immediately clear that the same could have been achieved by changing the condition from != to == and having an if-block with just the print statement, so this is really just a toy example illustrating how continue works.

As you saw in these few examples, there are often multiple ways in which for, while, break, continue, and if-else can be combined to achieve the same thing. While break and continue can be useful commands, they can also make code more difficult to read and understand. Therefore, they should only be used sparingly and when their usage leads to a simpler and more comprehensible code structure than a combination of for /while and if-else would do.

You are already familiar with Python binary operators that can be used to define arbitrarily complex expressions. For instance, you can use arithmetic expressions that evaluate to a number, or boolean expressions that evaluate to either True or False. Here is an example of an arithmetic expression using the arithmetic operators – and *:

x = 25 – 2 * 3 

Each binary operator takes two operand values of a particular type (all numbers in this example) and replaces them by a new value calculated from the operands. All Python operators are organized into different precedence classes, determining in which order the operators are applied when the expression is evaluated unless parentheses are used to explicitly change the order of evaluation. This operator precedence table shows the classes from lowest to highest precedence. The operator * for multiplication has a higher precedence than the – operator for subtraction, so the multiplication will be performed first and the result of the overall expression assigned to variable x is 19.

Here is an example for a boolean expression:

x = y > 12 and z == 3

The boolean expression on the right side of the assignment operator contains three binary operators: two comparison operators, > and ==, that take two numbers and return a boolean value, and the logical ‘and’ operator that takes two boolean values and returns a new boolean (True only if both input values are True, False otherwise). The precedence of ‘and’ is lower than that of the two comparison operators, so the ‘and’ will be evaluated last. So if y has the value 6 and z the value 3, the value assigned to variable x by this expression will be False because the comparison on the left side of the ‘and’ evaluates to False.

In addition to all these binary operators, Python has a ternary operator, so an operator that takes three operands as input. This operator has the format

 x if c else y

x, y, and c here are the three operands while ‘if’ and ‘else’ are the keywords making up the operator and demarcating the operands. While x and y can be values or expressions of arbitrary type, the condition c needs to be a boolean value or expression. What the operator does is it looks at the condition c and if c is True it evaluates to x, else it evaluates to y. So for example in the following line of code

 p = 1 if x > 12 else 0

variable p will be assigned the value 1 if x is larger than 12, else p will be assigned the value 0. Obviously what the ternary if-else operator does is very similar to what we can do with an if or if-else statement. For instance, we could have written the previous code as

p = 1 

if x > 12: 

    p = 0 

The “x if c else y” operator is an example of a language construct that does not add anything principally new to the language but enables writing things more compactly or more elegantly. That’s why such constructs are often called syntactic sugar. The nice thing about “x if c else y” is that in contrast to the if-else statement, it is an operator that evaluates to a value and, hence, can be embedded directly within more complex expressions as in the following example that uses the operator twice:

newValue = 25 + (10 if oldValue < 20 else 44) / 100 + (5 if useOffset else 0) 

Using an if-else statement for this expression would have required at least five lines of code.

In GEOG 485, we used the + operator for string concatenation to produce strings from multiple components to then print them out or use them in some other way, as in the following two examples:

print('The feature class contains ' + str(n) + ' point features.') 

queryString = '"'+ fieldName+ '" = ' + "'" + countryName + "'" 

An alternative to this approach using string concatenation is to use the string method format(…). When this method is invoked for a particular string, the string content is interpreted as a template in which parts surrounded by curly brackets {…} should be replaced by the variables given as parameters to the method. Here is how the two examples from above would look in this approach:

print('The feature class contains {0} point features.'.format(n) ) 

queryString = '"{0}" = \'{1}\''.format(fieldName, countryName) 

In both examples, we have a string literal '….' and then directly call the format(…) method for this string literal to give us a new string in which the occurrences of {…} have been replaced. In the simple form {i} used here, each occurrence of this pattern will be replaced by the i-th parameter given to format(…). In the second example, {0} will be replaced by the value of variable fieldName and {1} will be replaced by variable countryName. Please note that the second example will also use \' to produce the single quotes so that the entire template could be written as a single string. The numbers within the curly brackets can also be omitted if the parameters should be inserted into the string in the order in which they appear.

The main advantages of using format(…) are that the string can be a bit easier to produce and read as in particular in the second example, and that we don’t have to explicitly convert all non-string variables to strings with str(…). In addition, format allows us to include information about how the values of the variables should be formatted. By using {i:n}, we say that the value of the i-th variable should be expanded to n characters if it’s less than that. For strings, this will by default be done by adding spaces after the actual string content, while for numbers, spaces will be added before the actual string representation of the number. In addition, for numbers, we can also specify the number d of decimal digits that should be displayed by using the pattern {i:n.df}. The following example shows how this can be used to produce some well-formatted list output:

items = [('Maple trees', 45.232 ),  ('Pine trees', 30.213 ), ('Oak trees', 24.331)]

for i in items: 

    '{0:20} {1:3.2f}%'.format(i[0], i[1]) 

Output:

Maple trees                          45.23% 
Pine trees                           30.21% 
Oak trees                            24.33% 

The pattern {0:20} is used here to always fill up the names of the tree species in the list with spaces to get 20 characters. Then the pattern {1:3.2f} is used to have the percentage numbers displayed as three characters before the decimal point and two digits after. As a result, the numbers line up perfectly.

The format method can do a few more things, but we are not going to go into further details here. Check out this page about formatted output if you would like to learn more about this.

The parameters we have been using so far, for which we only specify a name in the function definition, are called positional parameters or positional arguments because the value that will be assigned to them when the function is called depends on their position in the parameter list: The first positional parameter will be assigned the first value given within the parentheses (…) when the function is called, and so on. Here is a simple function with two positional parameters, one for providing the last name of a person and one for providing a form of address. The function returns a string to greet the person with.

def greet(lastName, formOfAddress): 

      return 'Hello {0} {1}!'.format(formOfAddress, lastName) 


print(greet('Smith', 'Mrs.')) 

Output: 

Hello Mrs. Smith! 

Note how the first value used in the function call (“Smith”) in line 6 is assigned to the first positional parameter (lastName) and the second value (“Mrs.”) to the second positional parameter (formOfAddress). Nothing new here so far.

The parameter list of a function definition can also contain one or more so-called keyword arguments. A keyword argument appears in the parameter list as

<argument name> = <default value>   

A keyword argument can be provided in the function by again using the notation

<name> = <value/expression>

It can also be left out, in which case the default value specified in the function definition is used. This means keyword arguments are optional. Here is a new version of our greet function that now supports English and Spanish, but with English being the default language:

def greet(lastName, formOfAddress, language = 'English'): 

      greetings = { 'English': 'Hello', 'Spanish': 'Hola' }

      return '{0} {1} {2}!'.format(greetings[language], formOfAddress, lastName) 

 
print(greet('Smith', 'Mrs.')) 

print(greet('Rodriguez', 'Sr.', language = 'Spanish')) 

Output: 

Hello Mrs. Smith! 
Hola Sr. Rodriguez! 

Compare the two different ways in which the function is called in lines 8 and 10. In line 8, we do not provide a value for the ‘language’ parameter so the default value ‘English’ is used when looking up the proper greeting in the dictionary stored in variable greetings. In the second version in line 10, the value ‘Spanish’ is provided for the keyword argument ‘language,’ so this is used instead of the default value and the person is greeted with “Hola” instead of "Hello." Keyword arguments can be used like positional arguments meaning the second call could also have been

print(greet('Rodriguez', 'Sr.', 'Spanish'))

without the “language =” before the value.

Things get more interesting when there are several keyword arguments, so let’s add another one for the time of day:

def greet(lastName, formOfAddress, language = 'English', timeOfDay = 'morning'): 

      greetings = { 'English': { 'morning': 'Good morning', 'afternoon': 'Good afternoon' }, 

                    'Spanish': { 'morning': 'Buenos dias', 'afternoon': 'Buenas tardes' } } 

      return '{0}, {1} {2}!'.format(greetings[language][timeOfDay], formOfAddress, lastName) 

 
print(greet('Smith', 'Mrs.')) 

print(greet('Rodriguez', 'Sr.', language = 'Spanish', timeOfDay = 'afternoon')) 

Output:

Good morning, Mrs. Smith! 
Buenas tardes, Sr. Rodriguez!

Since we now have four different forms of greetings depending on two parameters (language and time of day), we now store these in a dictionary in variable greetings that for each key (= language) contains another dictionary for the different times of day. For simplicity reasons, we left it at two times of day, namely “morning” and “afternoon.” In line 7, we then first use the variable language as the key to get the inner dictionary based on the given language and then directly follow up with using variable timeOfDay as the key for the inner dictionary.

The two ways we are calling the function in this example are the two extreme cases of (a) providing none of the keyword arguments, in which case default values will be used for both of them (line 10), and (b) providing values for both of them (line 12). However, we could now also just provide a value for the time of day if we want to greet an English person in the afternoon:

print(greet('Rogers', 'Mrs.', timeOfDay = 'afternoon')) 

Output: 

Good afternoon, Mrs. Rogers! 

This is an example in which we have to use the prefix “timeOfDay =” because if we leave it out, it will be treated like a positional parameter and used for the parameter ‘language’ instead which will result in an error when looking up the value in the dictionary of languages. For similar reasons, keyword arguments must always come after the positional arguments in the definition of a function and in the call. However, when calling the function, the order of the keyword arguments doesn’t matter, so we can switch the order of ‘language’ and ‘timeOfDay’ in this example:

print(greet('Rodriguez', 'Sr.', timeOfDay = 'afternoon', language = 'Spanish'))

Of course, it is also possible to have function definitions that only use optional keyword arguments in Python.

Let us continue with the “greet” example, but let’s modify it to be a bit simpler again with a single parameter for picking the language, and instead of using last name and form of address we just go with first names. However, we now want to be able to not only greet a single person but arbitrarily many persons, like this:

greet('English', 'Jim', 'Michelle')

Output: 

Hello Jim! 
Hello Michelle! 

greet('Spanish', 'Jim', 'Michelle', 'Sam') 

Output: 

Hola Jim! 
Hola Michelle! 
Hola Sam! 

To achieve this, the parameter list of the function needs to end with a special parameter that has a * symbol in front of its name. If you look at the code below, you will see that this parameter is treated like a list in the body of the function:

def greet(language, *names): 

     greetings = { 'English': 'Hello', 'Spanish': 'Hola' } 

     for n in names: 

          print('{0} {1}!'.format(greetings[language], n))    

What happens is that all values given to that function from the one corresponding to the parameter with the * on will be placed in a list and assigned to that parameter. This way you can provide as many parameters as you want with the call and the function code can iterate through them in a loop. Please note that for this example we changed things so that the function directly prints out the greetings rather than returning a string.

We also changed language to a positional parameter because if you want to use keyword arguments in combination with an arbitrary number of parameters, you need to write the function in a different way. You then need to provide another special parameter starting with two stars ** and that parameter will be assigned a dictionary with all the keyword arguments provided when the function is called. Here is how this would look if we make language a keyword parameter again:

def greet(*names, **kwargs): 

     greetings = { 'English': 'Hello', 'Spanish': 'Hola' } 

     language = kwargs['language'] if 'language' in kwargs else 'English' 

     for n in names: 

          print('{0} {1}!'.format(greetings[language], n))      

If we call this function as

greet('Jim', 'Michelle')

the output will be:

Hello Jim! 
Hello Michelle!

And if we use

greet('Jim', 'Michelle', 'Sam', language = 'Spanish')

we get:

Hola Jim! 
Hola Michelle! 
Hola Sam! 

Yes, this is getting quite complicated, and it’s possible that you will never have to write functions with both * and ** parameters, still here is a little explanation: All non-keyword parameters are again collected in a list and assigned to variable names. All keyword parameters are placed in a dictionary using the name appearing before the equal sign as the key, and the dictionary is assigned to variable kwargs. To really make the ‘language’ keyword argument optional, we have added line 5 in which we check if something is stored under the key ‘language’ in the dictionary (this is an example of using the ternary "... if ... else ..." operator). If yes, we use the stored value and assign it to variable language, else we instead use ‘English’ as the default value. In line 9, language is then used to get the correct greeting from the dictionary in variable greetings while looping through the name list in variable names.

When making the transition from a beginner to an intermediate or advanced Python programmer, it also gets important to understand the intricacies of variables used within functions and of passing parameters to functions in detail. First of all, we can distinguish between global and local variables within a Python script. Global variables are defined outside of any function. They can be accessed from anywhere in the script and they exist and keep their values as long as the script is loaded which typically means as long as the Python interpreter into which they are loaded is running.

In contrast, local variables are defined inside a function and can only be accessed in the body of that function. Furthermore, when the body of the function has been executed, its local variables will be discarded and cannot be used anymore to access their current values. A local variable is either a parameter of that function, in which case it is assigned a value immediately when the function is called, or it is introduced in the function body by making an assignment to the name for the first time.

Here are a few examples to illustrate the concepts of global and local variables and how to use them in Python.

def doSomething(x):      # parameter x is a local variable of the function

     count = 1000 * x    # local variable count is introduced

     return count 

 

y = 10            # global variable y is introduced 

print(doSomething(y)) 

print(count)      # this will result in an error 

print(x)          # this will also result in an error

This example introduces one global variable, y, and two local variables, x and count, both part of the function doSomething(…). x is a parameter of the function, while count is introduced in the body of the function in line 3. When this function is called in line 11, the local variable x is created and assigned the value that is currently stored in global variable y, so the integer number 10. Then the body of the function is executed. In line 3, an assignment is made to variable count. Since this variable hasn’t been introduced in the function body before, a new local variable will now be created and assigned the value 10000. After executing the return statement in line 5, both x and count will be discarded. Hence, the two print statements at the end of the code would lead to errors because they try to access variables that do not exist anymore.

Now let’s change the example to the following:

def doSomething(): 

     count = 1000 * y    # global variable y is accessed here

     return count 

 

y = 10          

print(doSomething()) 

This example shows that global variable y can also be directly accessed from within the function doSomething(): When Python encounters a variable name that is neither the name of a parameter of that function nor has been introduced via an assignment previously in the body of that function, it will look for that variable among the global variables. However, the first version using a parameter instead is usually preferable because then the code in the function doesn’t depend on how you name and use variables outside of it. That makes it much easier to, for instance, re-use the same function in different projects.

So maybe you are wondering whether it is also possible to change the value of a global variable from within a function, not just read its value? One attempt to achieve this could be the following:

 def doSomething(): 

     count = 1000  

     y = 5 

     return count * y 

 

y = 10 

print(doSomething()) 

print(y)      # output will still be 10 here

However, if you run the code, you will see that last line still produces the output 10, so the global variable y hasn't been changed by the assignment in line 5. That is because the rule is that if a name is encountered on the left side of an assignment in a function, it will be considered a local variable. Since this is the first time an assignment to y is made in the body of the function, a new local variable with that name is created at that point that will overshadow the global variable with the same name until the end of the function has been reached. Instead, you explicitly have to tell Python that a variable name should be interpreted as the name of a global variable by using the keyword ‘global’, like this:

def doSomething(): 

     count = 1000 

     global y      # tells Python to treat y as the name of global variable

     y = 5         # as a result, global variable y is assigned a new value here

     return count * y 

 

y = 10 

print(doSomething()) 

print(y)       # output will now be 5 here

In line 5, we are telling Python that y in this function should refer to the global variable y. As a result, the assignment in line 7 changes the value of the global variable called y and the output of the last line will be 5. While it's good to know how these things work in Python, we again want to emphasize that accessing global variables from within functions should be avoided as much as possible. Passing values via parameters and returning values is usually preferable because it keeps different parts of the code as independent of each other as possible.

So after talking about global vs. local variables, what is the issue with mutable vs. immutable mentioned in the heading? There is an important difference in passing values to a function depending on whether the value is from a mutable or immutable data type. All values of primitive data types like numbers and boolean values in Python are immutable, meaning you cannot change any part of them. On the other hand, we have mutable data types like lists and dictionaries for which it is possible to change their parts: You can, for instance, change one of the elements in a list or what is stored under a particular key in a given dictionary without creating a completely new object.

What about strings and tuples? You may think these are mutable objects, but they are actually immutable. While you can access a single character from a string or element from a tuple, you will get an error message if you try to change it by using it on the left side of the equal sign in an assignment. Moreover, when you use a string method like replace(…) to replace all occurrences of a character by another one, the method cannot change the string object in memory for which it was called but has to construct a new string object and return that to the caller.

Why is that important to know in the context of writing functions? Because mutable and immutable data types are treated differently when provided as a parameter to functions as shown in the following two examples:

def changeIt(x): 

     x = 5   # this does not change the value assigned to y

 
y = 3 

changeIt(y) 

print(y)     # will print out 3

As we already discussed above, the parameter x is treated as a local variable in the function body. We can think of it as being assigned a copy of the value that variable y contains when the function is called. As a result, the value of the global variable y doesn’t change and the output produced by the last line is 3. But it only works like this for immutable objects, like numbers in this case! Let’s do the same thing for a list:

def changeIt(x): 

     x[0] = 5   # this will change the list y refers to

 
y = [3,5,7] 

changeIt(y)     

print(y)        # output will be [5, 5, 7]

The output [5,5,7] produced by the print statement in the last line shows that the assignment in line 3 changed the list object that is stored in global variable y. How is that possible? Well, for values of mutable data types like lists, assigning the value to function parameter x cannot be conceived as creating a copy of that value and, as a result, having the value appear twice in memory. Instead, x is set up to refer to the same list object in memory as y. Therefore, any change made with the help of either variable x or y will change the same list object in memory. When variable x is discarded when the function body has been executed, variable y will still refer to that modified list object. Maybe you have already heard the terms “call-by-value” and “call-by-reference” in the context of assigning values to function parameters in other programming languages. What happens for immutable data types in Python works like “call-by-value,” while what happens to mutable data types works like “call-by-reference.” If you feel like learning more about the details of these concepts, check out this article on Parameter Passing.

While the reasons behind these different mechanisms are very technical and related to efficiency, this means it is actually possible to write functions that take parameters of mutable type as input and modify their content. This is common practice (in particular for class objects which are also mutable) and not generally considered bad style because it is based on function parameters and the code in the function body does not have to know anything about what happens outside of the function. Nevertheless, often returning a new object as the return value of the function rather than changing a mutable parameter is preferable. This brings us to the last part of this section.

It happens quite often that you want to hand back different things as the result of a function, for instance four coordinates describing the bounding box of a polygon. But a function can only have one return value. It is common practice in such situations to simply return a tuple with the different components you want to return, so in this case a tuple with the four coordinates. Python has a useful mechanism to help with this by allowing us to assign the elements of a tuple (or other sequences like lists) to several variables in a single assignment. Given a tuple t = (12,3,2,2), instead of writing

top = t[0] 

left = t[1] 

bottom = t[2] 

right = t[3] 

you can write

top, left, bottom, right = t 

and it will have the exact same effect. The following example illustrates how this can be used with a function that returns a tuple of multiple return values. For simplicity, the function computeBoundingBox() in this example only returns a fixed tuple rather than computing the actual tuple values from a polygon given as input parameter.

def computeBoundingBox(): 

     return (12,3,41,32) 

 

top, left, bottom, right = computeBoundingBox() # assigns the four elements of the returned tuple to individual variables

print(top)    # output: 12

This section has been quite theoretical, but you will often encounter the constructs presented here when reading other people’s Python code and also in the rest of this course.

Here’s another simple script that finds all cells over 3500 meters in an elevation raster and makes a new raster that codes all those cells as 1. Remaining values in the new raster are coded as 0. By now, you’re probably familiar with this type of “map algebra” operation which is common in site selection and other GIS scenarios.

Just in case you’ve forgotten, the expression Raster(inRaster) tells arcpy that it needs to treat your inRaster variable as a raster dataset so that you can perform map algebra on it. If you didn't do this, the script would treat inRaster as just a literal string of characters (the path) instead of a raster dataset.

# This script uses map algebra to find values in an 
#  elevation raster greater than 3500 (meters). 

import arcpy 
from arcpy.sa import * 

# Specify the input raster 
inRaster = "C:/Data/Elevation/foxlake" 
cutoffElevation = 3500 

# Check out the Spatial Analyst extension 
arcpy.CheckOutExtension("Spatial") 

# Make a map algebra expression and save the resulting raster 
outRaster = Raster(inRaster) > cutoffElevation 
outRaster.save("C:/Data/Elevation/foxlake_hi_10") 

# Check in the Spatial Analyst extension now that you're done 
arcpy.CheckInExtension("Spatial") 

You can probably easily work out what this script is doing but, just in case, the main points to remember on this script are:

• Notice the lines of code that check out the Spatial Analyst extension before doing any map algebra and check it back in after finishing. Because each line of code takes some time to run, avoid putting unnecessary code between checkout and checkin. This allows others in your organization to use the extension if licenses are limited. The extension automatically gets checked back in when your script ends, thus some of the Esri code examples you will see do not check it in. However, it is a good practice to explicitly check it in, just in case you have some long code that needs to execute afterward, or in case your script crashes and against your intentions "hangs onto" the license.
• inRaster  begins as a string, but is then used to create a Raster object once you run Raster(inRaster). A Raster object is a special object used for working with raster datasets in ArcGIS. It's not available in just any Python script: you can use it only if you import the arcpy module at the top of your script.
• cutoffElevation  is a number variable that you declare early in your script and then use later on when you build the map algebra expression for your outRaster.
• The expression  outRaster = Raster(inRaster) > cutoffElevation  is saying, in plain terms, "Make a new raster and call it  outRaster. Do this by taking all the cells of the raster dataset at the path of  inRaster  that are greater than the number assigned to the variable  cutoffElevation."
• outRaster  is also a Raster object, but you have to call the method  outRaster.save()  in order to make it permanent on disk. The  save() method takes one argument, which is the path to which you want to save.

Copy the code above into a file called Lesson1A.py (or similar as long as it has a .py extension) in spyder or your favorite IDE or text editor and then save it.

We don’t need to do anything to this code to get it to work in Python 3, it will be fine just as it is. Feel free to check it against Analyze Tools for Pro if you like. Your results should say “Analyze Tools for Pro Completed Successfully” with the lack of warnings signifying that the code you supplied is compatible with Python 3.

Next, we'll convert the script to a Tool.

32-bit software or hardware can only directly represent and operate with numbers up to 2^32 and, hence, can only address up to a maximum of 4GB of memory (that is 2^32 = 4294967296 bytes). If the file system of your operating system is limited to 32-bit integers as well, this also means you cannot have any single file larger than 4GB either in memory or on disk (you can still page or chain larger files together though).

64-bit architectures don’t have this limit. Instead you can access up to 16 terabytes of memory and this is actually only the limit of current chip architectures which "only" use 44 bits which will change over time as software and hardware architectures evolve. Technically with a 64-bit architecture you could access 16 Exabytes of memory (2^64) and while not wanting to paraphrase Bill Gates, that is probably more than we’ll need for the foreseeable future.

There most likely won't be any innate performance benefits to be gained by moving from 32-bit and 64-bit unless you need that extra memory. While in principle, you can move larger amounts of data per time between memory and CPU with 64-bit, this typically doesn't result in significantly improved execution times because of caching and other optimization techniques used by modern CPUs. However if we start using programming models where we run many tasks at once, you might want more than 4GB allocated to those processes. For example if you had 8 tasks that all needed 500MB of RAM each – that’s very close to the 4GB limit in total (500MB * 8 = 4000MB). If you had a machine with more processors (e.g. 64) you would very easily hit the 32-bit 4GB limit as you would only be able to allocate 62.5MB of RAM per processor from your code.

Even with hardware architectures and operating systems mainly being 64-bit these days, a lot of software still is only available as 32-bit versions. 64-bit operating systems are designed to be backwards compatible with 32-bit applications, and if there is no real expected benefit for a particular software, the developer of the software may just as well decide to stick with 32-bit and avoid the efforts and costs that it would take to make the change to 64-bit or even support multiple versions of the software. ArcGIS Desktop is an example of a software that is only available as 32-bit and this is (most likely) not going to change anymore since ArcGIS Pro, which is 64-bit, fills that role. However, Esri also provides a 64-bit geoprocessing extension for ArcGIS Desktop which will be further described in the next section. However, this section is considered optional. You may read through it and learn how it can be set up and about what performance gain can be achieved using it or you may skip most of it and just have a look at the table with the computation time comparison at the very end of the section. But we strongly recommend that you do not try to install 64-bit geoprocessing and perform the steps yourself before you have worked through the rest of lesson and the homework assignment.

This section is provided for interest only - as it only applies to ArcGIS Desktop - not Pro (which is natively 64 bit). It is recommended that you only read/skim through the section and check out the computation time comparison at the end without performing the described steps yourself and then loop back to it at the end of the lesson if you have free time and are interested in exploring 64-bit Background Geoprocessing in ArcGIS Desktop.

A number of versions ago (since Server 10.1), Esri added support for 64-bit arcpy. Esri also introduced 64-bit geoprocessing using the 64-bit Background Geoprocessing patch which was part of 10.1 Service Pack 1 as an option in ArcGIS Desktop (Pro is entirely 64-bit) to work around these memory issues for large geoprocessing tasks. Not all tools support 64-bit geoprocessing within Desktop and there are some tricks to getting it installed so you can access it in Desktop. There is also a 64-bit arcpy geoprocessing library so you can run your code (any code) from the command line. Background Geoprocessing (64-bit) is still available as a separate installation on top of ArcGIS (see this ArcMap/Background Geoprocessing (64-bit) page) and we’ve provided a link for students to obtain it within Canvas. You'll find this link on the "64-bit Geoprocessing downloads for ArcGIS (optional)" page under Lesson 1 in Canvas.

As Esri hint in their documentation, 64-bit processing within ArcGIS Desktop requires that the tool run in the background. This is because it is running using a separate set of tools which are detached from the Desktop application (which is 32-bit). Personally, I rarely use Background Geoprocessing but I do make use of the 64-bit version of Python that it installs to run a lot of my scripts in 64-bit mode from the command line.

If you’ve typically run your code in the past from within an IDE (such as PythonWin, IDLE or spyder) or from within ArcGIS you might not be aware that you can also run that code from the command line by calling Python directly.

For ArcGIS Desktop you can start a regular command prompt and, using the standard Windows commands, change to that path where your Python script is located. Usually when you open a command window, it will start in your home folder (e.g. c:\users\yourname). We could dedicate an entire class to operating system commands but Microsoft has a good resource at this Windows Commands page for those who are interested.

We just need a couple of the commands listed there :

• cd : change directory. We use this to move around our folders. Full help at this Commands/cd page.
• dir : list the files and folders in my directory. Full help at this Commands/dir page. 

We’ll change the directory to where our code from section 1.6 is (e.g. mine is c:\wcgis\geog489\lesson1) and see how to run the script using the command line versions of 32-bit and 64-bit Python.

cd c:\wcgis\geog489\lesson1

If you downloaded and installed the 64-bit Background Geoprocessing from above you will have both 32-bit and 64-bit Python installed. We’ll use the 32-bit Python first which should be located at c:\python27\arcgis10.6\python.exe (substitute 10.6 by whichever version of ArcGIS you have installed).

There’s a neat little feature built into the command line where you can use the TAB key to autocomplete paths so you could start typing c:\python27\a and then hit TAB and you should see the path cycling through the various ArcGIS folders.

We’ll run our code using:

C:\python27\ArcGIS10.6\python.exe myScriptName.py

Where myScriptName.py is whatever you saved the code from section 1.6 as. You will now see the code run in the command window and pop up all of the same messages you would have seen if you had run it from an IDE (but not the AddMessages messages as they are only interpreted by ArcGIS).

To run the code against the 64-bit version of Python the command is almost identical except that you’ll use the x64 version of Python that has been installed by Background Geoprocessing. In my case that means the command is:

C:\python27\ArcGISx6410.5\python.exe myScriptName.py

Once your script finishes you’ll have a few time stamp values. Running that code from Section 1.6 through the 32-bit and 64-bit versions a few times we have some sample results below. The first runs of each are understandably slower as arcpy is imported for the first time. You have probably witnessed this behavior yourself as your code takes longer the very first time it runs.

We can see a couple of things with these results – they are a little inconsistent depending on what else my PC was doing at the time the code was running and, if you are looking at individual executions of the code, it is difficult to see which pieces of the code are slower or faster from time to time. This is the problem that we will solve later in the lesson when we look at profiling code where we examine how long each line of code takes to run.

You have probably noticed if you have a relatively modern PC (anything from the last several years) that when you open Windows Task Manager (from the bottom of the list when you press CTRL-ALT-DEL) and you click the Performance tab and right click on the CPU graph and choose Change Graph to -> Logical Processors you have a number of processors (or cores) within your PC. These are actually “logical processors“ within your main processor but they function as though they were individual processors – and we’ll just refer to them as processors here for simplicity.

Figure 1.10 Logical Processors 

Credit: ArcGIS Pro

Now because we have multiple processors, we can run multiple tasks in parallel at the same time instead of one at a time. There are two ways that we can run tasks at the same time – multithreaded and multiprocessing. We’ll look at the differences in each in the following but it’s important to know that arcpy doesn’t support multithreading but it does support multiprocessing. In addition, there is a third form of parallelisation called distributed computing, which involves distributing the task over multiple computers, that we will also briefly talk about. 

With all of these approaches to speeding up our code, what are the elements which will cause bottlenecks and slow us down?

Well, there are a few – these include the time to set up each of the processes for multiprocessing. Remember earlier we mentioned that because each process doesn’t share memory it needs a copy of the data to use. This will need to be copied to a memory location. Also as each process runs its own Python.exe instance, it needs to be launched and arcpy needs to be imported for each instance (although fortunately, multiprocessing takes care of this for us). Still, all of that takes time to start so our code won’t appear to do much at first while it is doing this housekeeping - and if we're not starting a lot of processes then we won't see enough of a speed up in processing to make up for those start-up time costs.

Other things that can slow us down are the speed of our RAM. Access times for RAM used to be measured in nanoseconds but now are measured in megahertz (MHz). The method of calculating the speed isn’t especially important but if you’re moving large files around in RAM or performing calculations that require getting a number out of RAM, adding, subtracting, multiplying, etc. and then putting the result into another location in RAM and you’re doing that millions of times very, very small delays will quickly add up to seconds or minutes. Another speedbump is running out of RAM. While we can allocate more than 4GB per process using 64-bit programming, if we don’t have enough RAM to complete all of the tasks that we might launch then our operating system will start swapping between RAM (which is fast) and our hard disk (which isn’t – even if it’s one of the solid state types – SSDs).

Speaking of hard disks, it’s very likely that we’re loading and saving data to them and as our disks are slower than our RAM and our processors, that is going to cause a delay. The less we need to load and save data the better, so good multiprocessing practice is to keep as much data as possible in RAM (see the caveat above about running out of RAM). The speed of disks is governed by a couple of factors; the speed that the motor spins (unless it is an SSD), seek time and the amount of cache that the disk has. Here is how these elements all work together to speed up (or slow down) your code. The hard disk receives a request for data from the operating system, which it then goes looking for. This is the seek time referring to how long it takes the disk to position the read head over the segment of disk the data is located on, which is a function of motor speed as well. Then once the file is found, it needs to be loaded into memory – cache - and then this is sent through to the process that needed the data. When data is written back to your disk, the reverse process takes place. The cache is filled (as memory is faster than disks ) and then the cache is written to the disk. If the file is larger than the cache, the cache gets topped up as it starts to empty until the file is written. A slow spinning hard disk motor or a small amount of cache can both slow down this process.

It’s also possible that we’re loading data from across a network connection (e.g. from a database or remotely stored files) and that will also be slow due to network latency – the time it takes to get to and from the other device on the network with the request and the result.

We can also be slowed down by inefficient code, for example, using too many loops or an inefficient if / else / elif statement that we evaluate too many times or using a mathematical function that is slower than its alternatives. We'll examine these sorts of coding bottlenecks - or at least how to identify them when we look at code profiling later in the lesson.

From the brief description in the previous section, you might have realized that there are generally two broad types of tasks – those that are input/output (I/O) heavy which require a lot of data to be read, written or otherwise moved around; and those that are CPU (or processor) heavy that require a lot of calculations to be done. Because getting data is the slowest part of our operation, I/O heavy tasks do not demonstrate the same improvement in performance from multiprocessing as CPU heavy tasks. The more work there is to do for the CPU the greater the benefit in splitting that workload among a range of processors so that they can share the load.

The other thing that can slow us down is outputting to the screen – although this isn’t really an issue in multiprocessing because printing to our output window can get messy. Think about two print statements executing at exactly the same time – you’re likely to get the content of both intermingled, leading to a very difficult to understand message. Even so, updating the screen with print statements is a slow task.

Don’t believe me? Try this sample piece of code that sums the numbers from 0-100.

# Setup _very_ simple timing. 
import time 
start_time = time.time() 

sum = 0 
for i in range(0,100): 
    sum += i 
    print(sum) 
 
# Output how long the process took. 
print ("--- %s seconds ---" % (time.time() - start_time)) 

If I run it with the print function in the loop the code takes 0.049 seconds to run on my PC. If I comment that print function out, the code runs in 0.0009 seconds.

4278
4371
4465
4560
4656
4753
4851
4950
--- 0.04900026321411133 seconds ---

runfile('C:/Users/jao160/Documents/Teaching_PSU/Geog489_SU_21/Lesson 1/untitled1.py', wdir='C:/Users/jao160/Documents/Teaching_PSU/Geog489_SU_21/Lesson 1')
--- 0.0009996891021728516 seconds ---

In Penn State's GEOG 485 course, we simulated 10,000 runs of the children's game Cherry-O to determine the average number of turns it takes. If we printed out the results, the code took a minute or more to run. If we skipped all but the final print statement the code ran in less than a second. We’ll revisit that Cherry-O example as we experiment with moving code from the single processor paradigm to multiprocessor. We’ll start with it as a simple, non arcpy example and then move on to two arcpy examples – one raster (our raster calculation example from before) and one vector.

Since you most likely did not take GEOG 485, you may want to have a quick look at the description.

Following is the original Cherry-O code.

# Simulates 10K game of Hi Ho! Cherry-O  
# Setup _very_ simple timing.  
import time 
 
start_time = time.time() 
import random 
 
spinnerChoices = [-1, -2, -3, -4, 2, 2, 10] 
turns = 0 
totalTurns = 0 
cherriesOnTree = 10 
games = 0 
 
while games < 10001: 
    # Take a turn as long as you have more than 0 cherries  
    cherriesOnTree = 10 
    turns = 0 
    while cherriesOnTree > 0: 
        # Spin the spinner  
        spinIndex = random.randrange(0, 7) 
        spinResult = spinnerChoices[spinIndex] 
        # Print the spin result      
        # print ("You spun " + str(spinResult) + ".")     
        # Add or remove cherries based on the result  
        cherriesOnTree += spinResult 
 
        # Make sure the number of cherries is between 0 and 10     
        if cherriesOnTree > 10: 
            cherriesOnTree = 10 
        elif cherriesOnTree < 0: 
            cherriesOnTree = 0 
            # Print the number of cherries on the tree         
        # print ("You have " + str(cherriesOnTree) + " cherries on your tree.")      
        turns += 1 
    # Print the number of turns it took to win the game  
    # print ("It took you " + str(turns) + " turns to win the game.")  
    games += 1 
    totalTurns += turns 
print("totalTurns " + str(float(totalTurns) / games)) 
# lastline = raw_input(">")  
# Output how long the process took.  
print("--- %s seconds ---" % (time.time() - start_time))  

We've added in our very simple timing from earlier and this example runs for me in about 1/3 of a second (without the intermediate print functions). That is reasonably fast and you might think we won't see a significant improvement from modifying the code to use multiprocessor mode but let's experiment. 

The Cherry-O task is a good example of a CPU bound task; we’re limited only by the calculation speed of our random numbers, as there is no I/O being performed. It is also an embarrassingly parallel task as none of the 10,000 runs of the game are dependent on each other. All we need to know is the average number of turns; there is no need to share any other information. Our logic here could be to have a function (Cherry-O) which plays the game and returns to our calling function the number of turns. We can add that value returned to a variable in the calling function and when we’re done divide by the number of games (e.g. 10,000) and we’ll have our average. 

Now that we have completed a non-ArcGIS parallel processing exercise, let's look at a couple of examples using ArcGIS functions. There are a number of caveats or gotchas to using multiprocessing with ArcGIS and it is important to cover them up-front because they affect the ways in which we can write our code.

Esri describe a number of best practices for multiprocessing with arcpy. These include:

• Use “in_memory“ workspaces to store temporary results because as noted earlier memory is faster than disk.
• Avoid writing to file geodatabase (FGDB) data types and GRID raster data types. These data formats can often cause schema locking or synchronization issues. That is because file geodatabases and GRID raster types do not support concurrent writing – that is, only one process can write to them at a time. You might have seen a version of this problem in arcpy previously if you tried to modify a feature class in Python that was open in ArcGIS. That problem is magnified if you have an FGDB and you’re trying to write many feature classes to it at once. Even if all of the feature classes are independent you can only write them to the FGDB one at a time.
• Use 64-bit. This isn’t an issue if we are writing code in ArcGIS Pro (although Esri does recommend using a version of Pro greater than 1.4) because we are already using 64-bit, but if you were planning on using Desktop as well, then you would need to use ArcGIS Server 10.5 (or greater) or ArcGIS Desktop with Background Geoprocessing (64-bit). The reason for this is that as we previously noted 64-bit processes can access significantly more memory and using 64-bit might help resolve any large data issues that don’t fit within the 32-bit memory limits of 4GB.

So bearing the top two points in mind we should make use of in_memory workspaces wherever possible and we should avoid writing to FGDBs (in our worker functions at least – but we could use them in our master function to merge a number of shapefiles or even individual FGDBs back into a single source).

As you may remember from GEOG 485, the simplest method of debugging is to embed print statements in your code to either determine how far your code is running through a loop or to print out the contents of a variable.

A more complicated and detailed method involves using the tools or features of your IDE to create watches for checking the contents of variables and breakpoints for stepping through your code. While we will introduce a range of IDEs at the end of the lesson and you will be investigating IDEs further and examining how their debugging tools work, we will provide a generic overview of the techniques here of setting watches, breakpoints and stepping through code. Don’t focus on the specifics of the interface as we do this. Instead, it is more important to understand the purpose of each of the different methods of debugging. While we could examine debugging using IDLE (or PythonWin as we did in GEOG 485) we’ve chosen to use a more representative IDE – Spyder as it more accurately demonstrates the features available in the majority of IDEs and you should also already have it installed from earlier.

We will start off by looking at Spyder’s debugging functions, which are similar to those of PythonWin from GEOG 485 and all other IDEs. You might remember back in GEOG 485 when we looked at debugging we examined setting up breakpoints, watches, and stepping through code. We’re going to revisit those concepts here briefly using Spyder; the same functionality will be available in all of the IDEs you’ll be investigating later but you might have to dig a little in the menus for it. There's more details in the Spyder help here.

For our debugging and profiling with Spyder, we’re going to use our raster multiprocessing example (the one that involved delineating streams from lidar data) from earlier. After starting Spyder and opening that file, click on the Debug menu and you will see a number of options.

Set the cursor to the line of code which filters our list of rasters :

new_rasters = filter_list(rasters,wildCardList) 

and set a breakpoint by choosing Set/Clear breakpoint from the Debug menu or pressing F12 or double-clicking alongside the line number of the code where you want a breakpoint (removing breakpoints uses the same methods as adding them).

Now run your code using the Debug item from the Debug menu, by pressing CTRL+F5, or using the button. Your code will now start running and then pause prior to running the line of code at your breakpoint. 

One of the nice touches to Spyder is that it has a variable explorer pane, accessible from the View > Panes menu. This will show you a list of your variables, their type, size, and contents once they’re declared, and you can watch them change.

If you bring this pane up you will see two variables already in it, start_time and rasters, our list which has a size of 43 – that is 43 elements - so the code found 43 rasters in my case.

We can step to the next line of code using either Debug>Step, CTRL+F10, or the Step button . Notice that when you do this, an extra variable will appear in your list (new_rasters) and your next executable line of code will be highlighted in a pale grey band.

We can step through our code line by line to monitor its state as well as the contents of our variables to ensure that our code is doing what we expect it to. While we are doing this, we are looking for unexpected results in variables, loops that are not running correctly (either doing too many or too few iterations or not at all) and if/elif/else statements which are not being correctly evaluated.

While we can step through lines individually, there are also two more useful options in most IDEs which are "run current function until it returns" and "run to the next breakpoint." Spyder has these implemented using the buttons and , respectively, and I suggest that you experiment with both to see how your functions are executed – particularly with the multiprocessing code. If you want to run between multiple breakpoints, you will need to add some more breakpoints to your code, otherwise, the code will run to the end.

Lastly, if you don’t want to step through all of the lines of code in a function, you can use the Step Return button  , which will return you to the calling function.

We have experimented with some very simple code profiling by introducing the time() function into our code in previous examples using it to record the start and end times of our code and check the overall performance.

While that gives us a high level view of the performance of our code, we do not know where specific bottlenecks might exist within that code. For example what is the slowest part of our algorithm, is it the file reading or writing or the calculations in between? When we know where these bottlenecks are we can investigate ways of removing them or at least using faster techniques to speed up our code.

In this section, we will focus on basic profiling that looks at how long each function call takes – not each individual line of code. This basic type of code profiling is built into Spyder (and most IDEs). However, we also provide some complementary materials that explain how to visualize profiling results and on how to profile each line of code individually but these parts will be entirely optional because they are quite a bit more complex and require the installation of additional software and packages. It is possible that you will run into some technical/installation issues, so our recommendation is that you only come back and try to run the described steps in these optional subsections yourself if you are done with lesson and homework assignment and still have time left.

Introduction to Online VCS

Some popular online hosting solutions for VCS and DVCS code repositories include: GitHub, Bitbucket, Google Code and Microsoft CodePlex. These online repositories are often used as the main trunk repositories for open-source projects with many developers who may be geographically dispersed. For the purposes of this class, we will focus on GitHub.

Introduction to GitHub

GitHub takes all of Git’s version control components, adds a graphical user interface to repositories, change history, and branch documentation, and adds several social components. Users can add comments, submit issues, and get involved in the bug tracking process. Users can follow other GitHub users or specific projects to be notified of project updates. GitHub can either be used entirely online or with an application download for easily managing and syncing local and online repositories. Optional (not required for class): Click here for the desktop application download.

The following exercise will cover the basics of Git and how they’re used in the GitHub website.

GitHub's change log

GitHub has the ability to display everything that changed with every commit. Take a look at GitHub's Kansasgis/NG911 page. If you click on one of the titles of one of the commits, it displays whatever basic description the developer included of the changes and then as you scroll down, you can see every code change that occurred - red highlighting what was removed and green highlighting what was added. If you mouse over the code, a plus sign graphic shows up, and users can leave comments and such.

Resolving conflicts on GitHub

Conflicts occur if two branches being merged have had different changes in the same places. Git automatically flags conflicts and will not complete the merge like normal; instead, the user will be notified that the conflicts must be resolved. Some conflicts can be resolved inside GitHub, and other types of conflicts have to be resolved in the Git command line [4]. Due to the complexity of resolving conflicts in the command line, it’s best to plan ahead and silo projects as much as possible to avoid conflicts.

Git adds three different markers to the code to flag conflicts:

<<<<<<<HEAD – This marker indicates the beginning of the conflict in the base branch. The code from the base branch is located directly under this marker.

======= – This marker divides the base branch code from the other branch.

>>>>>>> BRANCH-NAME – This marker will have the name of the other branch next to it and indicates the end of the conflict.

Here’s a full example of how Git flags a conflict between branches:

<<<<<<<HEAD
myString = “Monty Python and the Holy Grail is the best. ”
======= 
myString = “John Cleese is hilarious.”
>>>>>>> cleese-branch

To resolve the conflict, the user needs to pick what myString will equal. Possible resolution options-

Keeping the base branch -
myString = “Monty Python and the Holy Grail is the best.”

Using the other branch -
myString = “John Cleese is hilarious.”

Combining branches, in this case combining the options -
myString = “Monty Python and the Holy Grail is the best. John Cleese is hilarious.”

GitHub has an interface that can be activated for resolving basic conflicts by clicking on the “Resolve Conflicts” button under the “Pull Requests” tab. This interface steps through each conflict and the user must decide which version to take, keep their changes, use the other changes, or work out a way to integrate both sets of changes. Inside the GitHub interface, the user must also remove the Git symbols for the conflict. The user steps through every conflict in that particular file to decide how to resolve the conflict and then will eventually click on the “Mark as resolved” button. The next file in the project with conflicts will show up and the user will repeat all of the steps until the conflicts are resolved. At this point, the user will click “Commit merge” and then “Merge pull request.”

For more complex types of conflicts like one branch deleting a file that the other keeps, the resolution has to take place in the Git command line. This process can hopefully be avoided, but basic instructions are available at GitHub Help: Resolving a merge conflict using the command line.

Resources:
[4] https://help.github.com/articles/resolving-a-merge-conflict-on-github/

GitHub is a great fit for managing open source code projects since, with a free account, all repositories are available on the internet at large. For example, the open source GIS software QGIS (see Lesson 4) is housed on GitHub at GitHub's qgis/QGIS page. Take a look at the repository.

On the front page, you can see in the dashboard statistics that (at the time of this writing) there have been over 40,000 commits, 50 branches, 100 releases, and 250 contributors to the QGIS project. Users worldwide can now contribute their ideas, bugs, and code improvements to a central location that can be managed with standard version control workflows.

Some software companies that have traditionally been protective about their code have adopted GitHub to open certain projects. Esri is rather active on GitHub at GitHub's Esri page including the documentation and samples for the ArcGIS API for Python (see Lesson 3). Microsoft also is present at the GitHub Microsoft page with the tagline “Open source, from Microsoft with love.”

While GitHub is open to all digital files and any programming languages, Python is a great fit for use in GitHub for multiple reasons. Unlike other, heavier programming languages, Python doesn’t require extensive libraries with complex dlls and installation structures to get the job done.

Creating Python repositories is as simple as adding the .py files, and then the project can be shared, documented, and updated as needed. GitHub is also a great place to find both Python snippets and entire modules to use. For basic purposes, users can copy/paste just the portions of code off another project they want to try. Otherwise, users can fork an entire repository and tweak it as necessary to fit their purposes.

GitHub strongly recommends that every repository contain a README.txt or README.md file. This file will act as the “home page” for the project and is displayed on the repository page after files and folders are listed. This document should contain specific information about the project, how to use it, licensing, and support.

Text files will show up without formatting, so many users choose to use an .md (markdown) file instead. Markdown notation will be interpreted to show various formatting components like font size, bold, italics, imbedded links, numbered lists, and bullet points.

For more information on markdown formatting, visit GitHub Guide's Mastering Markdown page. We will also use Markdown in Lesson 3, in the context of Jupyter notebooks, and provide a brief introduction there.

Gists

While all free GitHub accounts are required to publish public repositories, all accounts have the ability to create Gists. Gists are single page repositories in GitHub, so they don't support projects with folder structures or multiple files. Since Gists are a single page repository, they are good for storing code snippets or one page projects. Gists can be public or private, even with a free account.

To create a Gist in GitHub, log into GitHub and then click on the plus sign in the upper right hand corner. In the options presented, choose "New gist." Enter a description of the Gist (in figure 1.30 "Delete if Exists Shortcut" is the description) as well as the filename with extension (in figure 1.30 this is DeleteIfExists.py). Enter code or notes in the large portion of the screen or import the code by using the "Add File" button. You have two options for saving your Gist- either "Create secret gist" or "Create public gist." 

"Secret" Gists are only mostly secret since they use the internet philosophy of difficult-to-guess urls. If you create a secret Gist, you can still share the Gist with anyone by sending them the url, but there are no logins required to view the Gist. Along this same philosophy, if someone stumbles across the url, they will be able to see the Gist.

For more information about Gists, see the official GitHub documentation at "About Gists" page on Github's website.

Figure 1.30 Example of creating a Gist, note the secret and public gist options

GeoJson

For GIS professionals, Gists are additionally useful since a Gist can be a single GeoJson file. GeoJson files are essentially a text version of geographic data in json formatting. Other developers can instantly access your GeoJson data and incorporate it from GitHub into their online mapping applications without needing to get a hard copy of the shapefile or geodatabase feature class or rely on some kind of map server. GitHub will automatically display GeoJson files as a map whether the file is a Gist or a part of a larger repository. For example, take a look at GitHub's lyzidiamond/learn-geojson page. At first, you’ll see the GeoJson file interpreted as a map. If you click the “Raw” button located on the upper right-ish side of the map, you will see what the GeoJson file looks like in text form. GeoJson can be easily used in Python since after reading in the file, Python can work with the text as if it is one giant dictionary.

Using GitHub in This Course

In GEOG 489, using GitHub to store the sample code and exercise code from the lessons can be a great way to practice and gain experience with a new software tool. Using GitHub is not required and we don't recommend that you store your completed projects on there. GitHub is an encouraged platform for students to learn since many organizations use GitHub or other VCS. 

Conclusion

Git and GitHub provide fast and convenient ways to track projects, whether the project is by one individual or a team of software developers. Although GitHub has many complex features available, it’s easily accessible for individual and small projects that need some kind of tracking mechanism. In addition to version control, GitHub provides users with a social platform for project management as well as the ability for users to create Gists and store GeoJson.  

We already mentioned some of the main visible components that serve as the construction blocks for building GUIs in typical GUI libraries, like windows and buttons. The image below shows a few more that we commonly encounter in today’s software applications including group boxes, labels, check boxes, radio buttons, combo boxes, line input fields, text input areas, tab views, and list views. Others that you are probably familiar with are tool bars, tool buttons, menu bars, context menus, status bars, and there are many, many more!

Figure 2.3 Different widgets we commonly find in GUIs

Typically, the GUI library contains classes for each of these visible elements, and they are often referred to as the widgets. Certain widgets can serve as containers for other widgets and, as a result, widgets tend to be organized hierarchically within a concrete graphical interface. For instance, a dialog box widget can contain many other widgets including a tab widget that in turn contains labels and buttons on each of its tab areas. If a widget A directly contains another widget B, we say that B is a child of A and A is B’s parent. A widget without a parent is a window that will be displayed independently on the screen. Widgets can have many different attributes for controlling their visual appearance, their layout behavior, and how they operate. Methods defined in the respective widget class allow for accessing and modifying these attributes. The most common operations performed with widgets in program code are:

• Creating the widget
• Adding the widget to another widget (widget becomes the child of that other widget)
• Adding another widget to the widget (the widget becomes the parent containing the other widget)
• Changing an attribute of the widget (for instance, you may change the text displayed by a label widget)
• Reading an attribute of the widget (for instance, you may need to get the text that a user entered into a line input widget)
• Setting the layout management method for the widget; this determines how the child widgets of that widget will be arranged to fill the widget’s content area
• Linking an event that the widget can trigger to event handler code that should be executed in that case (for instance, you may want that a particular function in your code be called when a particular button is clicked)

We will explain the ideas of layout management and event handling hinted at in the last two bullet points above in more detail in the next sections. From the user's perspective, widgets can be interacted with in many different ways depending on the type of the widget, including the following very common forms of interactions:

• The user can click on the widget with the mouse to start some action, change the state of the widget, or open a context menu.
• The user can give focus to a widget either by clicking it or using the TAB key for moving the focus to the next widget in a specified order. At any time, only a single widget can have focus, meaning it will receive keyboard input which allows for typing into a text input widget or "clicking" a button by pressing ENTER. Similarly, the user makes a widget lose focus when giving focus to another widget.
• The user can enter some text into the widget.
• The user can drag the widget and drop it onto another widget, or drop something on the widget.

In addition, there are complex widgets that allow the user to interact with them and change their state by clicking on particular parts of a widget. Examples are the user unfolding a combo box to select a different value, the user clicking on a menu in the menu bar to open that menu and select an item in it, the user moving the slider component of a widget to adapt some value, or the user selecting a color by clicking on a particular location in a widget with a chromatic circle. The events caused by such user interactions are what drives the order of code execution in the underlying program code as will be further explained in Section 2.4.3.

Widgets can be freely arranged within the content area of their parent widget (for instance a window widget). This can be done with the help of pixel coordinates that are typically measured from the top left corner of the content area. However, it would be very cumbersome to create GUIs by specifying the x and y coordinates and width and height values for each widget. More importantly, such a static layout will only look good for the particular size of the parent it has been designed for. When the containing window is resized by the user, you would expect that the dimensions and positions of the contained widgets adapt accordingly which will not be the case for such a static coordinate-based layout.

Therefore, in modern GUI libraries the task of arranging the child widgets within the parent widget is taken care of by so-called layout managers. This happens dynamically, so if the window is resized, all content will be rearranged again down the hierarchy of widgets by the different layout managers involved. The GUI library defines different layout classes to create the layout manager objects from. The three most common layout types are:

• Horizontal layout: All child widgets are arranged horizontally in the order in which they have been added to the parent
• Vertical layout: All child widgets are arranged vertically in the order in which they have been added to the parent
• Grid layout: Child widgets are arranged in a table or grid layout consisting of columns and rows. The child widgets are added to a particular cell in the layout by providing row and column indices and can potentially span multiple rows and/or columns.

The images below illustrate these three basic layout types for a set of three label and three push button widgets. Because of their layout preferences, both labels and buttons remain at their preferred height in the horizontal layout, but are expanded to fill the available space horizontally in the vertical layout. In the grid layout, the labels have been set up to form the first column of a grid with 3 rows and 2 columns, while the buttons have been set up to occupy the cells in the second column.

Figure 2.4 Horizontal layout of three labels and three buttons

Figure 2.5 Vertical layout of the same label and button widgets

Figure 2.6 Grid layout with two columns and three rows

In many cases, these three layout types are already sufficient to arrange the widgets as desired because of the way layouts can be nested. For instance, you can have the components of a window organized in a grid layout and then use a vertical layout to arrange several widgets within one of the cells of the grid. As indicated above, widgets have attributes that affect their layout behavior and, as a result, how much space the layout manager will assign to them in the vertical and horizontal dimensions taking into account the other child widgets of the same parent and their attributes. These attributes can, for instance, define the minimum, preferred, and/or maximum dimensions or general size policies for width and height. For instance, as we saw in the examples from the previous figures, a standard push button widget often has a size policy for its height that says that the height should not be increased beyond its default height even if there is space available, while in the horizontal dimension the button may be expanded to fill available space that the layout manager would like to fill. We will talk more about layout management when we will start to work with the QT library later on in this lesson.

So far, we are used to the code in a script file being executed line-by-line from the top to the bottom with the order of execution only being affected by loops, if-else and function calls. GUI based applications operate a bit differently. They use what is called an event-driven programming approach. In event-driven programming, the code is organized as follows:

• Initalization phase:
  • The GUI is created by instantiating the widgets (= creating objects of the widget classes) and organizing them in parent-child hierarchies using suitable layout manager objects to achieve the desired arrangement.
  • Event handling code is defined for dealing with events from user interactions (like clicking a button) or other types of events.
  • Different events are associated with the corresponding event handling code.
• Execution phase:
  • An infinite loop is started that waits for GUI events and only terminates if the application is closed. In the loop, whenever an event occurs, the respective event handling code is executed, then the waiting continues until the next event happens.

The order of the first two points of the initialization phase can sometimes be swapped. The code for running the event processing loop is something you do not have to worry about when programming GUI based applications because that part is being taken care of by the GUI library code. You just have to add a command to start the loop and be aware that this is happening in the background. Your main job is to produce the code for creating the GUI and defining the event handlers in the initialization part of the program.

As indicated above, widgets can be interacted with in different ways and such interactions cause certain types of events that can be reacted on in the code. For instance, a button widget may cause a “button pressed” event when the user presses the left mouse button while having the mouse cursor over that button and a “button released” event when the mouse button is released again. In addition, it will cause a “button triggered” event after the release but this event can also be triggered by pressing the RETURN key, while the button has focus (e.g. when the button is hovered over with the mouse). The functionality of the GUI is realized by setting up the event handler code. That code typically consists of the definitions of event handler functions that are invoked when a certain event is caused and that contain the code that determines what should happen in this case. For instance, we may set up an event handler function for the “button triggered” event of the mentioned button. The code of that function may, for example, open a dialog box to get further information from the user or start some computations followed by displaying the results.

Precisely how events are linked to event handling functions depends on the GUI library used. We will see quite a few examples of this later in this lesson. However, we already want to mention that in the QT library we are going to use, the event-based approach is covered under what QT calls the signals & slots approach. When an event occurs for a particular QT widget (e.g., user clicks button), that widget emits a signal specific for that event. A slot is a function that can be called in response to a signal (so essentially an event handler function). QT’s widgets have predefined slots so that it is possible to directly connect a signal of one widget to a slot of another widget. For instance, the “clicked” signal of a button can be connected to the “clear” slot of a text widget, such that the text of that widget is cleared whenever the button is clicked. In addition, you can still write your own slot functions and connect them to signals to realize the main functionality of your application. No worries if this all sounds very abstract at the moment; it will get clear very quickly as soon as we look at some concrete examples.

As we already mentioned, Tkinter is the standard GUI for Python, but only in the sense that it is a package in the Python standard library, so it is available for all platforms without requiring any additional installation. Its name stands for “Tk interface”. It is certainly possible that you have not heard about Tk and Tcl before but Tk is one of the oldest free and open-source, cross-platform GUI toolkits (written in the Tcl scripting language and initially released in 1991) and has been adopted for building GUIs in many programming languages. Tkinter has been written by Fredrik Lundh and is essentially a set of wrapper classes and functions that use a Tcl interpreter embedded into the Python interpreter to create and manage the Tk GUI widgets.

To get an impression of how tkinter is used to build a GUI in Python, let us look at the example of creating a simple miles-to-kilometers conversion tool. The tool is shown in the figure below. It has a single window with five different widgets: two label widgets, two widgets for entering or displaying single lines of text, and a button in the middle. The user can enter a number of miles into the line input field at the top, then press the button, and then the entered number of miles will be converted and displayed as kilometers in the line input field at the bottom.

Figure 2.7 Tkinter version of the miles-to-kilometer converter GUI

We are using a line input field to display the resulting distance in kilometers just to make things more symmetrical. Since we do not want the user to be able to enter anything into this text field, it has been disabled for input and we could just as easily have used another label widget. The Python code to create this tool with the help of tkinter is shown below and the explanation of the code follows.

from tkinter import Tk, Label, Entry, Button, DISABLED, StringVar 

def convert(): 
    """Takes miles entered, converts them to km, and displays the result""" 
    miles = float(entryMiles.get()) 
    kilometers.set(str(miles * 1.60934)) 

# create the GUI 

rootWindow = Tk() # create main window 
rootWindow.title("Miles to kilometers") 
rootWindow.geometry('500x200+0+0') 
rootWindow.grid_columnconfigure(1, weight = 1) 

labelMiles = Label(rootWindow, text='Distance in miles:') # create label for miles field 
labelMiles.grid(row=0, column=0) 

labelKm = Label(rootWindow, text='Distance in kilometers:') # create label for km field 
labelKm.grid(row=2, column=0) 

entryMiles = Entry(rootWindow) # create entry field for miles 
entryMiles.grid(row=0, column=1, sticky='w,e') 

kilometers = StringVar() # create entry field for displaying km 
entryKm = Entry(rootWindow, textvariable = kilometers, state=DISABLED) 
entryKm.grid(row=2, column=1, sticky='w,e') 

convertButton = Button(rootWindow, text='Convert', command = convert) # create button for running conversion 
convertButton.grid(row=1, column=1) 

# run the event processing loop  

rootWindow.mainloop() 

Let us ignore the first few lines of Python code for a moment and first look at lines 10 to 29. This is where the GUI of our little program is produced starting with the root window widget in lines 10 to 13. The widget is created by calling the function Tk() defined in tkinter and the created object is stored in variable rootWindow. We then use different methods of the widget to set its title, initial size, and some properties for its grid layout that we are going to use to arrange the child widgets within the content area of the root window.

Next, the label saying “Distance in miles:” is created. The tkinter widget class for labels is called Label and we provide rootWindow as a parameter to Label(…), so that the widget knows what its parent widget is. As mentioned, we will be using a grid layout, namely one with three rows and two columns. We place the created label in the cell in the first row and first column of its parent by calling the grid(…) method with row = 0 and column = 0. We then take the exact same steps to create the other label and place it in the third row of the first column.

In the next steps, the two text input fields are created as widget objects of the tkinter Entry class. An additional parameter sticky=’w,e’ is used for placing these widgets in the grid. This parameter says that the widgets should expand horizontally (west and east) to fill the entire cell. This is required to make the layout fill out the window horizontally and have the text field grow and shrink when the window is resized. Moreover, the Entry widget for displaying the distance in kilometers is set to DISABLED so that the user cannot enter text into it, and it is associated with a variable kilometers of tkinter class StringVar which is needed for us to be able to change the text displayed in the widget from code.

Finally, the button is created as a widget of tkinter class Button. What is new here is what happens with the ‘command’ parameter given to Button(…) in line 28. Here we are saying that if this button is clicked, the function convert() that we are defining at the top of our code should be executed to deal with this event. So this is an example of connecting an event to an event handler function. What happens in convert() is very simple: With the help of the get() method, we get the current text from the Entry widget for the distance in miles, multiply it with a constant to convert it to kilometers, and then use the set() method of the StringVar object in variable kilometers to change the text displayed in the Entry widget to the distance in kilometers associated with that variable.

In the last line of the code, we call the mainloop() method of our root window to start the infinite event processing loop. The program execution will only return from this call when the user closes the root window, in which case the program execution will be terminated.

The only part of the code we haven’t talked about is the first line where we simply import the widget classes and other auxiliary classes from tkinter that we need in our code.

Hopefully, it is clear that this is just a very prototypical implementation of a miles-to-kilometers conversion tool focusing on the GUI. We have neither implemented any sort of checking whether input values are valid nor any sort of error handling. It is therefore very easy to make the tool crash, e.g. by entering something that is not a number into the field for distance in miles. If you haven’t already done so, we suggest you create a Python script with the code from above and try out the tool yourself and see how the layout adapts if you resize the window. Feel free to experiment with making small changes to the code, like adapting the text shown by the labels or adding another button widget to the currently still empty second row of the first column; then make the button call another event handler function you write to, for instance, just print some message to the console.

Don’t worry if some of the details happening here don’t seem entirely clear at the moment. A real introduction to creating GUIs with Python will follow later in this lesson. Here we just wanted to give you a general idea of how the different concepts we discussed in Section 2.4 are realized in tkinter: You saw how different widgets are created, how they were arranged in a grid layout by placing them in different cells of the layout, how to connect an event (button clicked) with a self-defined event handler function (convert()), and how to execute the application by starting the event processing loop (rootWindow.mainloop()). Now let’s move on and talk about QT as an alternative to tkinter and see how this same example would look like when produced with the PyQt instead of tkinter.

We already mentioned a few things about QT in this lesson. It is a widely used cross-platform library written in C++, modern and under very active development. In addition to the GUI functionality, the library provides support for internationalization, Unicode, database and network access, XML and JSON code processing, thread management, and more. That’s why it is also called an application framework, not just a GUI library. QT was originally developed by the company Trolltech and its initial release was in 1995. KDE, one of the early GUIs for the Linux operating system, was based on QT and that triggered a lot of discussion and changes to the license and organization QT was published under. These days, the company developing QT is called The QT Company, a successor of Trolltech, and QT is published in four different editions, including the Community edition that is available under different open source licenses GPL 3.0, LGPL 3.0, and LPGL 2.1 with a special QT exception. QT is very commonly used for both open source and commercial software, and if you have worked with QT in one programming language, it is typically relatively easy to learn to use it in a different language. QT5 was released in 2012 and the current version of QT at the time of this writing is 5.10.

To familiarize ourselves with PyQT5 and get to know the most common QT5 widgets, let’s go through three smaller examples.

While it’s good and useful to understand how to write Python code to create a GUI directly, it’s obviously a very laborious and time consuming approach that requires writing a lot of code. So together with the advent of early GUI frameworks, people also started to look into more visual approaches in which the GUI of a new software application is clicked and dragged together from predefined building blocks. The typical approach is that first the GUI is created within the graphical GUI building tool and then the tool translates the graphical design into code of the respective programming language that is then included into the main code of the new software application.

The GUI building tool for QT that we are going to use in this lesson is called QT Designer. QT Designer is included in the PyQT5 Python package. The tool itself is platform and programming language independent. Instead of producing code from a particular programming language, it creates a .ui file with an xml description of the QT-based GUI. This .ui file can then be translated into Python code with the pyuic5 GUI compiler tool that also comes with PyQT5. There is also support for directly reading the .ui file from Python code in PyQT5 and then generating the GUI from its content, but we here will use the approach of creating the Python code with pyui5 because it is faster and allows us to see and inspect the Python code for our application. However, you will see an example of reading in the content of the .ui file directly in Lesson 4. In the following, we will take a quick look at how the QT Designer works.

Since you already installed PyQt5 from the ArcGIS Pro package manager, the QT Designer executable will be in your ArcGIS Pro default Python environment folder, either under "C:\Users\<username>\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\Library\bin\designer.exe" or under “C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Library\bin\designer.exe”. It might be good idea to create a shortcut to this .exe file on your desktop for the duration of the course, allowing you to directly start the application. After starting, QT Designer will greet you as shown in the figure below:

Figure 2.21 QT Designer

QT Designer allows for creating so-called “forms” which can be the GUI for a dialog box, a main window, or just a simple QWidget. Each form is saved as a single .ui file. To create a new form, you pick one of the templates listed in the “New Form” dialog. Go ahead and double-click the “Main Window” template. As a result, you will now see an empty main window form in the central area of the QT Designer window.

Figure 2.22 Empty form for a new QMainWindow

Let’s quickly go through the main windows of QT Designer. On the left side, you see the “Widget Box” pane that lists all the widgets available including layout widgets and spacers. Adding a widget to the current form can be done by simply dragging the widget from the pane and dropping it somewhere on the form. Go ahead and place a few different widgets like push buttons, labels, and line edits somewhere on the main window form. When you do a right-click in an empty part of the central area of the main window form, you can pick “Lay out” in the context menu that pops up to set the layout that should be used to arrange the child widgets. Do this and pick “Lay out horizontally” which should result in all the widgets you added being arranged in a single row. See what happens if you instead change to a grid or vertical layout. You can change the layout of any widget that contains other widgets in this way.

On the right side of QT Designer, there are three different panes. The one at the top called “Object Inspector” shows you the hierarchy of all the widgets of the current form. This currently should show you that you have a QMainWindow widget with a QWidget for its central area, which in turn has several child widgets, namely the widgets you added to it. You can pretty much perform the same set of operations that are available when interacting with a widget in the form (like changing its layout) with the corresponding entry in the “Object Inspector” hierarchy. You can also drag and drop widgets onto entries in the hierarchy to add new child widgets to these entries, which can sometimes be easier than dropping them on widgets in the form, e.g., when the parent widget is rather small.

The “Object” column lists the object name for each widget in the hierarchy. This name is important because when turning a GUI form into Python code, the object name will become the name of the variable containing that widget. So if you need to access the widget from your main code, you need to know that name and it’s a good idea to give these widgets intuitive names. To change the object name to something that is easier to recognize and remember, you can double-click the name to edit it, or you can use “Change objectName” from the context menu when you right-click on the entry in the hierarchy or the widget itself.

Below the “Object Inspector” window is the “Property Editor”. This shows you all the properties of the currently selected widget and allows you to change them. The yellow area lists properties that all widgets have, while the green and blue areas below it (you may have to scroll down to see these) list special properties of that widget class. For instance, if you select a push button you added to your main window form, you will find a property called “text” in the green area. This property specifies the text that will be displayed on the button. Click on the “Value” column for that property, enter “Push me”, and see how the text displayed on the button in the main window form changes accordingly. Some properties can also be changed by double-clicking the widget in the form. For instance, you can also change the text property of a push button or label by double-clicking it.

If you double-click where it says “Type Here” at the top, you can add a menu to the menu bar of the main window. Give this a try and call the menu “File”.

Figure 2.23 Adding a File menu to the menu bar by double-clicking on "Type Here"

The last pane on the right side has three tabs below it. “Resource browser” allows for managing additional resources, like files containing icons to be used as part of the GUI. “Action editor” allows for creating actions for your GUI. Remember that actions are for things that can be initiated via different GUI elements. If you click the “New” button at the top, a dialog for creating a new action opens up. You can just type in “test” for the “Text” property and then press “OK”. The new action will now appear as actionTest in the list of actions. You can drag it and drop it on the File menu you created to make this action an entry in that menu.

Figure 2.24 A new action created in the Action Editor window

Finally, the “Signal/Slot Editor” tab allows for connecting signals of the widgets and actions created to slots of other widgets. We will mainly connect signals with event handler functions in our own Python code rather than in QT Designer, but if some widgets’ signals should be directly connected to slots of other widgets this can already be done here.

You now know the main components of QT Designer and a bit about how to place and arrange widgets in a form. We cannot teach QT Designer in detail here but we will show you an example of creating a larger GUI in QT Designer as part of the following first walkthrough of this lesson. If you want to learn more, there exist quite a few videos explaining different aspects of how to use QT Designer on the web. For now, here is a short video (12:21min) that we recommend checking out to see a few more basic examples of the interactions we briefly described above before moving on to the walkthrough.

Once you have created the form(s) for the GUI of your program and saved them as .ui file(s), you can translate them into Python code with the help of the pyuic5 tool, e.g. by running a command like the following from the command line using the tool directly with Python :

"C:\Users\<username>\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\python.exe" –m PyQt5.uic.pyuic mainwindow.ui –o mainwindow.py

or

"C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe" –m PyQt5.uic.pyuic mainwindow.ui –o mainwindow.py

, depending on where your default Python environment is located (don't forget to replace <username> with your actual user name in the first version and if you get an error with this command try typing it in, not copying/pasting). mainwindow.ui, here, is the name of the file produced with QT Designer, and what follows the -o is the name of the output file that should be produced with the Python version of the GUI. If you want, you can do this for your own .ui file now and have a quick look at the produced .py file and see whether or not you understand some of the things happening in it. We will demonstrate how to use the produced .py file to create the GUI from your Python code as part of the walkthrough in the next section.

The idea of the little software tool we are going to build is that sometimes you would like to be able to quickly create a feature class with points for certain places based on their names. Let’s say you want to produce a shapefile showing a few places in the city you are living in. What you could do is open a basemap in ArcGIS or some other GIS software and then digitize point features for these particular locations. Or you could look up the coordinates of the places in some external application such as Google Maps and then create new point features based on these coordinates manually or by exporting the coordinates in a format that can be imported by your GIS. Both these options are relatively time consuming and require several steps. We therefore want to exploit the fact that there exist many web services that allow for querying for geographic places based on name, type, and other properties, and getting back the coordinates. The two web services we are going to work with are the Nominatim query interface to OpenStreetMaps (OSM) and the place query interface of the online geographic gazetteer GeoNames. We will discuss these services in more detail below. Our tool should allow us to enter a query term such as a place name, pick a query service and set some additional query options, and then run the query to get a list of candidate places with coordinates. We can then look at the result list and pick one or more candidates from the list and either

• add them to a layer currently open in ArcGIS Pro (assuming that the software is run as a script tool inside ArcGIS),
• add them to a shapefile on disk, or
• add the result to a CSV text file (basically an alternative that still works even when arcpy is not available).

Figure 2.20 and the video below show the GUI of the software and further illustrate how the tool works. The use case in the shown example is that we use the tool to populate a layer open in ArcGIS Pro with some of the main sightseeing locations in Paris including the Eiffel Tower, the Louvre, etc.

Figure 2.25 The interface of the Locations from Web Services tool

Please watch the following video, Location from Web Service Tool Overview (5:25min):

The GUI is organized into four main parts:

• The text field for entering the query term. We are currently querying for "Eiffel tower".
• The section where you pick a web query service and set some additional options. The figure shows the interface for the Nominatim OSM service. In addition to Nominatim and GeoNames, we also have a “Direct Input” tab for adding places directly by providing a name and lat/lon coordinates.
• Since the query has already been run, the "Results" section shows different candidate toponyms returned by the Nominatim service for our query. The third one looks like the result we want, so we selected just that one. The selection can be changed with the help of the buttons below or by directly checking or unchecking the checkboxes. Next to the result list is a browser widget that is used to display a Leaflet based web map of the results.
• The final section is for adding the selected features to some dataset. The figure shows the tab for adding them to a currently open layer in ArcGIS Pro.  

We can run as many queries as we wish with the tool and collect the results we want to keep in one of the possible output options. While we will only present a basic version of the tool in this walkthrough supporting two different query interfaces and direct input, the tool could easily be extended to provide access to other web portals so that it can be used for running queries to, for instance, get locations of all Starbucks located in a particular city or create locations from a list of addresses.

Before we continue with exploring the GUI and code of the application, let us talk about the two web services we are going to use:

Nominatim – Surely, you know about OpenStreetMaps (OSM) and how it collects geographic data from volunteers all over the world to create a detailed map of the world. OSM data is freely available and can be directly exported from the OSM web site. In addition, there exist quite a few web services built around the OSM data, some of them created with the purpose of allowing for querying the data to only obtain information about particular features in the data. One such example is the Nominatim web service provided by OSM themselves. The website for Nominatim Open Street Maps provides an easy to use interface to the Nominatim query engine. You can type in a search term at the top and then will get to see a list of results displayed on the left side of the page and an OSM based map on the right that shows the currently selected entity from the result list. Give it a try, for instance by running a query for “Eiffel tower, France”.

Note: With some recent changes, it seems Nominatim has become much more restrictive and will often only return a single result rather than multiple options. If you leave out the 'France' in the query, the only result returned will actually not be the Eiffel Tower in Paris. However, you will still get multiple results if you, for instance, enter 'Washington' as the query string. Due to these changes, the results you will get when using Nominatim in the Locations from Web Services Tool will partially deviate from what is shown in the figures and videos in this section (for instance, the list of options shown in Figure 2.25 when only querying for 'Eiffel tower' without country name.

Nominatim provides an HTTP based web API that can be used for running queries from your own code and getting the results back, for instance as JSON or XML data. The web API is explained on this wiki page here. Query terms and additional parameters are encoded in a URL that has the general format:

https://nominatim.openstreetmap.org/search?

where parameters are specified as <parameter name>=<parameter value> pairs and multiple parameters have to be separated by an & symbol. The parameter name for the query string is simply 'q' (so q=...). To run the query, the client sends an HTTP GET request with this URL to the Nominatim server who processes the query and parameters, derives the result, and sends back the results to the client. The format parameter controls how the result is presented and encoded. Without specifying that parameter, you will get the kind of HTML page that you already saw above with the result list and map. When using format=json as we will do in the following, we get the result as a list of entities encoded as JSON. Here is an example query URL, querying for “Eiffel tower, France” again, that you can test out simply by clicking on the link to open the URL in your browser:

https://nominatim.openstreetmap.org/search?q=Eiffel%20tower%2c%20France&format=json

Have a look at the result shown in your browser. As we explained at the beginning of the lesson, JSON uses [...] to denote lists of entities

[ <entity1, entity2>, ... ]

where each entity is described by its properties like in a Python dictionary:

{<property1>: <value>, <property2>: <value>, ...}

Due to the changes, we mentioned above, the result will be a list with just a single entity, looking like this  ...

[{
    "place_id":115316817,
    "licence":"Data © OpenStreetMap contributors, ODbL 1.0. http://osm.org/copyright",
    "osm_type":"way",
    "osm_id":5013364,
    "lat":"48.8582599",
    "lon":"2.2945006358633115",
    "class":"man_made",
    "type":"tower",
    "place_rank":30,
    "importance":0.5868325701744196,
    "addresstype":"man_made",
    "name":"Eiffel Tower",
    "display_name":"Eiffel Tower, 5, Avenue Anatole France, Quartier du Gros-Caillou, 7th Arrondissement, Paris, Ile-de-France, Metropolitan France, 75007, France",
    "boundingbox":["48.8574753","48.8590453","2.2933119","2.2956897"]
}]

We can see that the most important properties for us will be the ‘display_name’ property and the ‘lat’ and ‘lon’ properties (appearing in bold above) in order to create point features and add them to an existing data set with our tool. Feel free to compare this result to what you get when querying for 'Washington' where you will get a list of multiple results.

The following query uses a few more parameters to query for places called London in Canada (countrycodes=CA) and asking for only a single entity be returned (limit=1).

https://nominatim.openstreetmap.org/search?q=London&format=json&countrycodes=CA&limit=1

If you look at the result, you will see that it lists London in Ontario as the only result. Without using the ‘countrycodes’ parameter the result would have been London in the UK because Nominatim uses a ranking scheme to order entities by likelihood/prominence. Without the 'limit' parameter, we would get a list of multiple options in the JSON result. ‘format’, ‘countrycodes’ and ‘limit’ will be the only parameters we will be using in our tool but please have a look at the other parameters and examples given on the Nominatim wiki page to get an idea of what other kinds of queries could be implemented.

GeoNames – You have already seen the Nominatim examples, so we can keep the section about GeoNames a bit shorter, since the URLs for running queries are somewhat similar. GeoNames is essentially an online geographic gazetteer, so a directory of geographic place names with associated information including coordinates. The main page can be used to type in queries directly but we will be using their REST web API that is documented here. Instead of a parameter for specifying the output format, the API uses a special URL for running queries that are replied to using JSON, as in the following example:

http://api.geonames.org/searchJSON?name=Springfield&maxRows=10&username=demo

Please note that a parameter (name=) is used to specify the query term. In addition, GeoNames requires a user name to be provided with the ‘username’ parameter. In case you tried out the link above, you probably got the reply that the daily request limit for user ‘demo’ has been reached. So you will have to create your own account at http://www.geonames.org/login and then use that user name instead of ‘demo’ in the query. The JSON sent back by GeoNames as the result will start like this:

{
    "totalResultsCount": 3308,
    "geonames": [
        {
            "adminCode1": "IL",
            "lng": "-89.64371",
            "geonameId": 4250542,
            "toponymName": "Springfield",
            "countryId": "6252001",
            "fcl": "P",
            "population": 116565,
            "countryCode": "US",
            "name": "Springfield",
            "fclName": "city, village,...",
            "adminCodes1": {
                "ISO3166_2": "IL"
            },
            "countryName": "United States",
            "fcodeName": "seat of a first-order administrative division",
            "adminName1": "Illinois",
            "lat": "39.80172",
            "fcode": "PPLA"
        },
        {
            "adminCode1": "MO",
            "lng": "-93.29824",
            "geonameId": 4409896,
            "toponymName": "Springfield",
            "countryId": "6252001",
            "fcl": "P",
            "population": 166810,
            "countryCode": "US",
            "name": "Springfield",
            "fclName": "city, village,...",
            "adminCodes1": {
                "ISO3166_2": "MO"
            },
            "countryName": "United States",
            "fcodeName": "seat of a second-order administrative division",
            "adminName1": "Missouri",
            "lat": "37.21533",
            "fcode": "PPLA2"
        }
    ]
}

Here the list of entities is stored under the property called ‘geonames’. Each entity in the list has the properties ‘toponymName’ with the entity name, ‘lng’ with the longitude coordinate, and ‘lat’ with the latitude coordinate. Query parameters we will be using in addition to ‘name’ and ‘username’ are ‘maxRows’ to determine the number of results sent back, ‘country’ to restrict the search to a single country, and ‘featureClass’ to look only for features of a particular type (codes A,H,L,P,R,S,T,U,V of the GeoNames feature class codes defined here).

The GUI of our software application is not too complex but still uses many of the most common GUI elements such as a toolbar with a tool button, normal push buttons, group boxes, tabs, labels, text input fields, checkboxes, combo boxes, list views, and a status bar. You already got to know most of these in Section 2.6 but there are also some new ones. The GUI consists of two parts, the GUI for the main window and the GUI for the dialog box that is shown when the user clicks on the “Create new shapefile…” button. Therefore, we will be producing two different .ui files with QT Designer called gui_main.ui and gui_newshapefile.ui. The GUI of the main window uses a vertical layout to organize the different elements into four different rows as shown in the figure below.

Each of the rows is formed by a QGroupBox widget and then other widgets are arranged hierarchically within these group boxes using a combination of grid and horizontal layouts. The tool doesn’t have a menu bar but a toolbar at the top with a button to exit the program and a status bar at the bottom. When creating the GUI in QT Designer, it will be important to name the widgets we need to refer to from our main code as indicated by the orange labels in the two figures below. As we already pointed out, the object names given to the widgets in QT Designer will be the names of the variables used for storing the widgets when the .ui file is compiled into a .py file.

Figure 2.26 Location from Web Services main GUI with labels showing the object names of the different widgets

Figure 2.27 GUI of the "Create shapefile" dialog

The series of videos linked below shows the process of creating the GUI in QT Designer. A zip file with the resulting .ui files is available for download here. We recommend that you work along with the video, pausing it as needed, to create the GUI yourself to get some more practice with QT Designer. The downloadable .ui files are mainly intended as a fallback position in case you experience any difficulties while creating the GUI or later on when compiling the .ui files and using the produced Python code. If you cannot replicate what is shown in the video, please ask for help on the discussion forums.

Create the GUI along with this series of videos [~45 min of video materials]

At this point, you either have saved your own two .ui files or, if you ran into any issues, will continue with these files downloaded from the link posted above. We now need to compile the files into Python code with pyuic5. We do this by running the ArcGIS Pro python.exe from the Windows command line with the pyuic module. The python.exe file is located in the folder “C:\Users \<username>\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\" (unless the ArcGIS Pro Python environment you are using is located in a different folder). So open a command shell, navigate to the folder containing the .ui files and then run the following two commands (again picking the correct version depending on where your default Python environment is installed and replacing <username> as needed):

"C:\Users\<username>\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\python.exe" –m PyQt5.uic.pyuic gui_main.ui –o gui_main.py

"C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe" –m PyQt5.uic.pyuic gui_main.ui –o gui_main.py

and

"C:\Users\<username>\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\python.exe" –m PyQt5.uic.pyuic gui_newshapefile.ui –o gui_newshapefile.py

"C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe" –m PyQt5.uic.pyuic gui_newshapefile.ui –o gui_newshapefile.py

Figure 2.28 Pyuic5 commands to compile the two .ui files into Python files

The parameters given to pyuic5 are the name of the input .ui file and then –o followed by the name of the output file. You should now have the two files gui_main.py and gui_newshapefile.py in the project folder. Let us have a quick look at the produced code. Open the produced file gui_main.py in your preferred Python IDE and see whether you recognize and understand how the different elements are created and how their properties are set. Without going into the details, the code defines a class Ui_MainWindow with a method setupUi(…). The parameter MainWindow is for passing a QMainWindow widget to the method. The rest of the code of the method then either...

• changes properties of MainWindow,
• creates new widgets and layouts storing them as attributes of the Ui_MainWindow object and sets their properties, or
• adds widgets to MainWindow or to other widgets to create the hierarchical organization of the widgets.

class Ui_MainWindow(object): 

    def setupUi(self, MainWindow): 
        … 
        MainWindow.resize(605, 685)   # changes property of main window widget 
        self.centralwidget = QtWidgets.QWidget(MainWindow) # creates child widgets and stores them in attributes of the Ui_MainWindow object 
        … 
        MainWindow.setCentralWidget(self.centralwidget) # adds widgets to main window 
        … 

This all means that we can create a new QMainWindow widget in our code with ...

mainWindow = QMainWindow()

... create an object of the UI_MainWindow class with

ui = gui_main.Ui_MainWindow()

... and then create the GUI for the main window by calling ui.setupUi(…) with mainWindow as the parameter:

ui.setupUi(mainWindow) 

Whenever we need to access a widget created by setupUi(…), we can do so by using the expression

ui.<object name of the widget>

where the object name is the name we gave the widget in QT Designer, e.g.

ui.queryTermLE

for the QLineEdit widget we created for entering the query term.

Now we are going to develop the main code for our tool that imports the gui_main.py and gui_newshapefile.py files, sets up the application, and connects the different GUI elements with event handler functions and functions that realize the actual functionality of the tool. We will organize the code into several modules. In particular, we will keep the functions that realize the main functionality, such as querying the different web portals, creating a new shapefile, etc. in a separate Python script file called core_functions.py . These functions are completely independent of the GUI of our tool and the few global variables we will need, so we might want to use them in other projects. Separating the project cleanly into GUI dependent and GUI independent code fosters reusability of the GUI independent code. Overall, the project will consist of the following Python files:

• gui_main.py – compiled version of gui_main.ui for setting up the GUI of the main window,
• gui_newshapefile.py – compiled version of gui_newshapefile.ui for setting up the GUI of the dialog box for creating a new shapefile,
• core_functions.py – contains definitions of functions for main functionality that are completely independent of GUI and global variables
• main.py – contains the event handler functions and the code for wiring up the GUI as well as setting up and running the application

In the following, we will focus on the code from main.py but we will start with a brief look at core_functions.py, so please download the file core_functions.py and open it so that you can read the code for the different functions. Most of the functions defined in the script should be rather easy to understand from the comments in the code and from your experience with using arcpy to work with shapefiles. Here is an overview of the functions with a few additional explanations:

• queryNominatim(query, limit=10, countryCodes=''): The purpose of this function is to query the Nominatim query interface of OSM for a query string and parameters for the limit of results returned and country code for the country to search in as discussed in Section 2.7.1. The code should be easy to understand: urllib.parse(…) is used to encode the query string to be used as part of a URL, before the final query URL is put together. The get(…) method of the requests package is used to run the actual query and the returned result is translated into JSON before being returned. Since no error handling is done in this function, we will have to deal with potential exceptions raised within the function in the code that calls this function.
• queryGeonames(query, maxRows = 10, username='', country='', featureClass=''): This function does the same as queryNominatim(...) but for GeoNames and for a different set of query parameters. Since GeoNames returns JSON code with the list of candidate features stored under the attribute 'geonames', we return json['geonames'] at the end.
• getStringFieldsForDescribeObject(desc): This is the first of a few auxiliary functions for getting information from shapefiles and layers. An arcpy.Describe object of a data set needs to be passed to it and it then returns a list of all editable string fields of that layer.
• getValidFieldsForShapefile(fileName): This function is used to get a list of editable string fields for a shapefile whose name is provided as a parameter. It relies on getStringFieldsForDescribeObject(…) to do most of the work but before that checks that the shapefile exists and is a Point vector data set. If not, it will return an empty list.
• createPointWGS1984Shapefile(fileName, fieldName): This function creates a new Point shapefile with the name provided as parameter, using WGS 1984 coordinates, and with a single additional string field whose name is also provided as a parameter.
• getPointLayersFromArcGIS(): This function is for getting a list of Layer objects for the layers currently open in ArcGIS Pro but only including Point layers.
• importArcpyIfAvailable(): This function returns True if arcpy is available for import, else False. It attempts to import arcpy within a try-except construct so that, if the operation fails, the resulting exception is caught and False is returned.
• runningAsScriptTool(): This function returns True if the program is run as a script tool inside ArcGIS, else False. This is accomplished by calling arcpy.mp.ArcGISProject("CURRENT") and catching the resulting exception if this operation fails, meaning the program is run as a standalone program outside of ArcGIS.
• webMapFromDictionaryList(features): This function produces and returns the HTML code for displaying the web map as part of the GUI. It gets the features that should be displayed on the map as parameter in the form of a list of dictionaries with name, lat, and lon attributes.

Now that you know the functions we have available for realizing the different operations that we will need, let’s develop the code for main.py together. Open a new file main.py in your IDE, then follow the steps listed on the next few pages.

As was mentioned earlier, when you use the import statement you can import a single .py file. In addition, the import statement can point to a folder containing a set of .py files, or a library written in a different programming language, such as C++.

You may be wondering how Python finds the module or package you specified since you only specify the name. Your own modules may have been in the current directory with the program using it, but arcpy, for example, isn’t. What happens is that Python has a list of locations that it uses in order to find the necessary packages. It traverses the list in the specific order, and, once it finds all the packages it needs, it stops looking. Here is the search order¹ that Python uses:

1. The home directory where the currently executing .py file is located 
2. PYTHONPATH directories: PYTHONPATH is a variable that is optionally set in the operating system. For example, on a Windows machine you would set it in Environmental Variables in System Settings. 
3. Standard library directories: The location where all standard libraries are installed on the local machine – if you have ArcGIS 10.6 Desktop installed to a default location the standard libraries can be found at C:\Python27\ArcGIS10.6\Lib. Browse to the folder, and take a quick look of all the packages that are installed. 
4. The contents of any .pth files: These are text files that can be created to add additional library paths; this option is used by ArcGIS Desktop 10.6. In a standard installation you can find the .pth file at C:\Python27\ArcGIS10.6\Lib\site-packages\Desktop10.6.pth. ArcGIS pro has its own .pth file: C:\Users\<username>\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone\Lib\site-packages\ArcGISPro.pth or C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\Lib\site-packages\ArcGISPro.pth depending on your version of Pro.
5. The site-package home or third-party extensions: Packages placed in Libs\site-packages directory. In a standard install of ArcGIS 10.6 Desktop that would be in the C:\Python27\ArcGIS10.6\Lib\site-packages folder. 

Because of the way Python finds the code it needs to import, you need to be careful how you name your modules and packages, and where you place them. For example, if you were to create an arcpy.py module and put it in the home directory, the ArcGIS arcpy package would not be loaded.

This list above may look intimidating, but the good news is that packages you are likely to need will be packaged with special Python utilities (either pip or conda) and thus setup to place themselves in the appropriate paths without any manual intervention on your part, beyond the installation step. The other good news is that both pip and conda are fairly straightforward to use when it comes to installing packages and managing Python environments. Creating your own pip or conda packages can be a bit more involved though as you will also see in this section but still provides a convenient way for deploying and sharing your own Python applications.

¹Mark Lutz: Learning Python, 5th Edition

There are many Python packages available for use, and there are a couple of different ways to effectively manage (install, uninstall, update) packages. The two package managers that are commonly used are pip and conda. In the following sections, we will discuss each of them in more detail. At the end of the section, we will discuss the merits of the two tools and make recommendations for their use.

We will be doing some more complicated technical "stuff" here so the steps might not work as planned because everyone’s PC is configured a little differently. If you get stuck please check in with the instructor sooner rather than later. A quick troubleshooting / debugging process can involve testing to see if running the command or Command Prompt as Administrator resolves the issue, trying the Windows Command prompt instead of the Python Command prompt (or vice versa), and, if none of that has helped, trying the tech support staple of restarting your PC.

Another option for packaging and distributing your Python programs is to use conda (we will discuss Anaconda a bit later in the lesson). Just like pip, it is a package manager. In addition, it is also an environment manager. What that means is that you can use conda to create virtual environments for Python, while specifying the packages you want to have available in that environment. A little more about that in a moment. Conda comes installed with ArcGIS Pro. While conda should be installed if you were able to install spyder in Lesson 1, we can doublecheck that it is by opening the Python Command Prompt and then typing in:

cd Scripts

followed by:

conda –-version

The output should show the conda version.

In order to find out what packages are installed type in:

conda list

Your output should look something like Figure 2.34:

     

Figure 2.34 Conda Package List

The first column shows the package name, the second the version of the package. The third column provides clues on how the package was installed. You will see that for some of the packages installed, Esri is listed, showing they are related to the Esri installation. The list option of conda is useful, not only to find out if the package you need is already installed but also to confirm that you have the appropriate version.

Conda has the functionality to create different environments. Think of an environment as a sandbox – you can set up the environment with a specific Python version and different packages. That allows you to work in environments with different packages and Python versions without affecting other applications. The default environment used by conda is called base environment. We do not need to create a new environment, but, should you need to, the process is simple – here is an example:

conda create -n gisenv python=3.6 arcpy numpy

the –n flag is followed by the name of the environment (in this case gisenv), then you would choose the Python version which matches the one you already have installed (3.5, 3.6 etc.) and follow that up with a list of packages you want to add to it. If you later find out you need other packages to be added, you could use the install option of conda, for example:

conda install –n gisenv matplotlib 

To activate an environment, you would run:

activate gisenv

And to deactivate an environment, simply:

deactivate 

There are other options you can use with environments – you can clone them and delete them, for example. A great resource for the different options is Conda's Managing Environments page.

If you worked or read through the exercises we worked through for packaging our application in the optional parts of this section, you might have gotten the impression that pip is easier to use than conda, and, since they both facilitate software packaging and distribution, why bother packaging your applications with conda? Here are a number of reasons:

1. Conda has a built-in virtual environment functionality
2. Conda further aids in package installs as it not only checks for dependencies but installs any dependent packages
3. Conda is closely bound with Anaconda and Anaconda Cloud, which is set up to use different channels, providing a finer-grained package organization. For example, if you are interested in packages ESRI has published, you can go to the company’s channel - https://anaconda.org/esri.
4. While pip can only be used with Python, conda can be used with other languages as well – so if you decided to develop something with R, conda can work with it too.

If you are interested in a more detailed comparison of the two tools, a great article to reference is Conda Myths and Misconceptions.  

features = [ { "id": "A", "lat": "23.32", "lon": "-54.22" }, 
             { "id": "B", "lat": "24.39", "lon": "53.11" }, 
             { "id": "C", "lat": "27.98", "lon": "-54.01" } ]
             
featuresAsTuples = [ ("Feature " + feat['id'], float(feat['lat']), float(feat['lon']) ) for feat in features if float(feat['lon']) < 0 ] 

print(featuresAsTuples) 

Let's look at the components of the list comprehension starting with the middle part:

for feat in features

This means we will be going through the list features using a variable feat that will be assigned one of the dictionaries from the features list. This also means that both the if-condition on the right and the expression for the 3-tuples on the left need to be based on this variable feat.

if float(feat['lon']) < 0 

Here we implement the condition that we only want 3-tuples in the new list for dictionaries that contain a lon value that is < 0. 

("Feature " + feat['id'], float(feat['lat']), float(feat['lon']) )

Finally, this is the part where we construct the 3-tuples to be placed in the new list based on the dictionaries contained in variable feat. It should be clear that this is an expression for a 3-tuple with different expressions using the values stored in the dictionary in variable feat to derive the three elements of the tuple. The output produced by this code will be:

Output:

[('Feature A', 23.32, -54.22), ('Feature C', 27.98, -54.01)] 

import requests
from bs4 import BeautifulSoup

url = 'https://www.e-education.psu.edu/geog489/node/2269'

response = requests.get(url)
soup = BeautifulSoup(response.text, 'html.parser')

divElement = soup.select('article#node-book-2269 > div > div')[0] 

wordLengths = [ len(word) for word in divElement.text.split() ] 
print(wordLengths)

After loading the html page and creating the BeautifulSoup structure for it as in the examples you already saw in this lesson, the select(…) method is used in line 9 to get the <div> elements within the <div> element within the <article> element with the special id we are looking for. Since we know there will only be one such element, we can use the index [0] to get that element from the list and store it in variable divElement.

With divElement.text.split() we create a list of all the words in the text and then use this inside the list comprehension in line 11 where we convert the word list into a list of word lengths by applying the len(…) function to each word.

The image below shows the hierarchy of widgets created in QT Designer and the names chosen. You can also download the .ui file and the compiled .py version here and open the .ui file in QT Designer to compare it to your own version. Note that we are using a vertical layout as the main layout, a horizontal layout within that layout for the +, -, =, and Clear button, and a grid layout within the vertical layout for the digit buttons.

Calculator hierarchy

The following code can be used to set up the GUI and run the application but without yet implementing the actual calculator functionality. You can run the code and main window with the GUI you created will show up.

import sys 

from PyQt5.QtWidgets import QApplication, QWidget 

import calculator_gui 

# create application and gui 

app = QApplication(sys.argv) 

mainWindow = QWidget()        # create an instance of QWidget for the main window 
ui = calculator_gui.Ui_Form() # create an instance of the GUI class from calculator_gui.py 
ui.setupUi(mainWindow)        # create the GUI for the main window 

# run app 

mainWindow.show() 
sys.exit(app.exec_()) 

The main .py file for the full calculator can be downloaded here. Function digitClicked(…) is the auxiliary function from step 5, called by the ten event handler functions digitPBClicked(…) defined later in the code. Function evaluateResult() is the auxiliary function from step 6. Function operatorClicked(…) is called from the two event handler functions plusPBClicked() and minusPBClicked(), just with different strings representing the mathematical operation. There are plenty of comments in the code, so it should be possible to follow along and understand what is happening rather easily.

It would be impossible to introduce or even just list all the packages available for conducting spatial data analysis projects in Python here, so the following is just a small selection of those that we consider most important.

numpy

numpy (Python numpy page, Wikipedia numpy page) stands for “Numerical Python” and is a library that adds support for efficiently dealing with large and multi-dimensional arrays and matrices to Python together with a large number of mathematical operations to apply to these arrays, including many matrix and linear algebra operations. Many other Python packages are built on top of the functionality provided by numpy.

matplotlib

matplotlib (Python matplotlib page, Wikipedia matplot page) is an example of a Python library that builds on numpy. Its main focus is on producing plots and embedding them into Python applications. Take a quick look at its Wikipedia page to see a few examples of plots that can be generated with matplotlib. We will be using matplotlib a few times in this lesson’s walkthrough to quickly create simple map plots of spatial data.

SciPy

SciPy (Python SciPy page, Wikipedia SciPy page) is a large Python library for application in mathematics, science, and engineering. It is built on top of both numpy and matplotlib, providing methods for optimization, integration, interpolation, signal processing and image processing. Together numpy, matplotlib, and SciPy roughly provide a similar functionality as the well known software Matlab. While we won’t be using SciPy in this lesson, it is definitely worth checking out if you're interested in advanced mathematical methods.

pandas

pandas (Python pandas page, Wikipedia pandas software page) provides functionality to efficiently work with tabular data, so-called data frames, in a similar way as this is possible in R. Reading and writing tabular data, e.g. to and from .csv files, manipulating and subsetting data frames, merging and joining multiple data frames, and time series support are key functionalities provided by the library. A more detailed overview on pandas will be given in Section 3.8. 

Shapely

Shapely (Python Shapely page, Shapely User Manual) adds the functionality to work with planar geometric features in Python, including the creation and manipulation of geometries such as points, polylines, and polygons, as well as set-theoretic analysis capabilities (intersection, union, …). It is based on the widely used GEOS library, the geometry engine that is used in PostGIS, which in turn is based on the Java Topology Suite  (JTS) and largely follows the OGC’s Simple Features Access Specification.

geopandas

geopandas (Python geopandas page, GeoPandas page) combines pandas and Shapely to facilitate working with geospatial vector data sets in Python. While we will mainly use it to create a shapefile from Python, the provided functionality goes significantly beyond that and includes geoprocessing operations, spatial join, projections, and map visualizations.

GDAL/OGR

GDAL/OGR (Python GDAL page, GDAL/OGR Python) is a powerful library for working with GIS data in many different formats widely used from different programming languages. Originally, it consisted of two separated libraries, GDAL (‘Geospatial Data Abstraction Library‘) for working with raster data and OGR (used to stand for ‘OpenGIS Simple Features Reference Implementation’) for working with vector data, but these have now been merged. The gdal Python package provides an interface to the GDAL/OGR library written in C++. In Section 3.9 and the lesson’s walkthrough, you will see some examples of applying GDAL/OGR.

ArcGIS API for Python

As we already mentioned at the beginning, Esri provides its own Python API (ArcGIS for Python page) for working with maps and GIS data via their ArcGIS Online and Portal for ArcGIS web platforms. The API allows for conducting administrative tasks, performing vector and raster analyses, running geocoding tasks, creating map visualizations, and more. While some services can be used autonomously, many are tightly coupled to Esri’s web platforms and you will at least need a free ArcGIS Online account. The Esri API for Python will be further discussed in Section 3.10.

In this lesson, we will start to work with some software that you probably are not familiar with and we will be using Python packages extensively that we have not used before to demonstrate how a complex GIS project can be solved in Python by combining different languages and packages within a Jupyter Notebook. Therefore, it is probably a good idea to prepare you a bit with an overview of what will happen in the remainder of the lesson.

• We already discussed the idea of using Jupyter Notebooks for data analysis projects. We will start this part of the lesson by introducing you to Jupyter Notebook and explaining to you the basic functionality (Section 3.6) so that you will be able to use it for the remainder of the lesson and future Python projects.
• The R programming language has its roots in statistical computing but also comes with a large library of packages providing data analysis methods for many specialized areas. One such package is the ‘dismo’ package for species distribution modeling. We will use the task of generating a species distribution model for the Solanum Acaule plant species as the data analysis task for this lesson’s walkthrough with the goal of showing you how Python and R functions can be combined within a Jupyter Notebook to solve some pretty complex analysis problem. The species distribution modeling application will be discussed further together with a brief overview on R and the ‘dismo’ package in Section 3.7.
• Using pandas for the manipulation of tabular data will be a significant part of this lesson’s walkthrough. We will use it to clean up the somewhat messy observation data available for Solanum Acaule. As a preparation, we will teach you the basics of manipulating table data with pandas in Section 3.8.
• GDAL/OGR will be the main geospatial extension of Python that we will use in this lesson (a) to perform additional data cleaning based on spatial querying and (b) to prepare additional input data (raster data sets for different climatic variables). We, therefore, provide an overview on its functionality and typical patterns of using GDAL/OGR in Section 3.9.
• We will mainly use the Esri ArcGIS API for Python to create an interactive map visualization within a Jupyter Notebook. However, the API has much more to offer and provides an interesting bridge between the FOSS Python Data Science ecosystem and the proprietary Esri world. We, therefore, provide an overview of the API in Section 3.10.
• The lesson’s walkthrough in Section 3.11 will show you a solution to the task of creating a species distribution model for Solanum Acaule combining both Python and R and making use of the different Python packages introduced in the lesson. The walkthrough will be provided as a Jupyter notebook that you can download and run on your own computer.

Juypter Notebook is already installed in the Anaconda environment AC36 or AC37 we created in Section 3.2. If you have the Anaconda Navigator running, make sure it shows the “Home” page and that our AC36 or AC37 environment is selected. Then you should see an entry for Juypter Notebook with a Launch button (Figure 3.3). Click the ‘Launch’ button to start up the application. This will ensure that your notebook starts with the correct environment. Starting Jupyter Notebook this way will create the link to shortcuts for a Jupyter Notebook to that conda environment so you can use it in the way we describe below.

Figure 3.3: Opening Jupyter Notebook from the Anaconda Navigator using the correct conda environment

Alternatively, you can start up Jupyter directly without having to open Anaconda first: You will find the Juypter Notebook application in your Windows application list as a subentry of Anaconda. Be sure that you start the Jupyter Notebook for the recently created conda environment (which will only be created if you change the dropdown in Anaconda Navigator above). Alternatively, simply press the Windows key and then type in the first few characters of Jupyter until Jupyter (with the correct conda environment) shows up in the search results.

When you start up Jupyter, two things will happen: The server component of the Jupyter application will start up in a Windows command line window showing log messages, e.g. that the server is running locally under the address http://localhost:8888/ (see Figure 3.4 (a)). When you start Jupyter from the Anaconda Navigator, this will actually happen in the background and you won't get to see the command line window with the server messages.  In addition, the web-based client application part of Jupyter will open up in your standard web browser showing you the so-called Dashboard, the interface for managing your notebooks, creating new ones, and also managing the kernels. Right now it will show you the content of the default Jupyter home folder, which is your user’s home folder, in a file browser-like interface (Figure 3.4 (b)).

Figure 3.4 (a): Shell window for the Jupyter server

The file tree view allows you to navigate to existing notebooks on your disk, to open them, and to create new ones. Notebook files will have the file extension .ipynb . Let’s start by creating a new notebook file to try out the things shown in the next sections. Click the ‘New…’ button at the top right, then choose the ‘Python 3’ option. A new tab will open up in your browser showing an empty notebook page as in Figure 3.5.

Figure 3.5: New notebook file in the browser

Before we are going to explain how to edit and use the notebook page, please note that the page shows the title of the notebook above a menu bar and a toolbar that provide access to the main operations and settings. Right now, the notebook is still called ‘Untitled...’, so, as a last preparation step, let’s rename the notebook by clicking on the title at the top and typing in ‘MyFirstJupyterNotebook’ as the new title and then clicking the ‘Rename’ button (Figure 3.6).

Figure 3.6: Giving your Jupyter Notebook a name

If you go back to the still open ‘Home’ tab with the file tree view in your browser, you can see your new notebook listed as MyFirstJupyterNotebook.ipynb and with a green ‘Running’ tag indicating that this notebook is currently open. You can also click on the ‘Running’ tab at the top to only see the currently opened notebooks (the ‘Clusters’ tab is not of interest for us at the moment). Since we created this notebook in the Juypter root folder, it will be located directly in your user’s home directory. However, you can move notebook files around in the Windows File Explorer if, for instance, you want the notebook to be in your Documents folder instead. To create a new notebook directly in a subfolder, you would first move to that folder in the file tree view before you click the ‘New…’ button.

Figure 3.7: The new notebook in the file tree view

We will now explain the basics of editing a Jupyter Notebook. We cannot cover all the details here, so if you enjoy working with Jupyter and want to learn all it has to offer as well as all the little tricks that make life easier, the following resources may serve as good starting points:

• Jupyter Documentation
• The Jupyter Notebook
• 28 Jupyter Notebook tips, tricks, and shortcuts
• Jupyter Notebook Tutorial: The Definitive Guide

A Jupyter notebook is always organized as a sequence of so called ‘cells’ with each cell either containing some code or rich text created using the Markdown notation approach (further explained in a moment). The notebook you created in the previous section currently consists of a single empty cell marked by a blue bar on the left that indicates that this is the currently active cell and that you are in ‘Command mode’. When you click into the corresponding text field to add or modify the content of the cell, the bar color will change to green indicating that you are now in ‘Edit mode’. Clicking anywhere outside of the text area of a cell will change back to ‘Command mode’.

Let’s start with a simple example for which we need two cells, the first one with some heading and explaining text and the second one with some simple Python code. To add a second cell, you can simply click on the  symbol. The new cell will be added below the first one and become the new active cell shown by the blue bar (and frame around the cell’s content). In the ‘Insert’ menu at the top, you will also find the option to add a new cell above the currently active one. Both adding a cell above and below the current one can also be done by using the keyboard shortcuts ‘A’ and ‘B’ while in ‘Command mode’. To get an overview on the different keyboard shortcuts, you can use Help -> Keyboard Shortcuts in the menu at the top.

Both cells that we have in our notebook now start with “In [ ]:” in front of the text field for the actual cell content. This indicates that these are ‘Code’ cells, so the content will be interpreted by Jupyter as executable code. To change the type of the first cell to Markdown, select that cell by clicking on it, then change the type from ‘Code’ to ‘Markdown’ in the dropdown menu  in the toolbar at the top. When you do this, the “In [ ]:” will disappear and your notebook should look similar to Figure 3.8 below. The type of a cell can also be changed by using the keyboard shortcuts ‘Y’ for ‘Code’ and ‘M’ for ‘Markdown’ when in ‘Command mode’.

Figure 3.8 Notebook with two cells with the second cell being a 'Code' cell

Let’s start by putting some Python code into the second(!) cell of our notebook. Click on the text field of the second cell so that the bar on the left turns green and you have a blinking cursor at the beginning of the text field. Then enter the following Python code:

from bs4 import BeautifulSoup 
import requests 

documentURL = 'https://www.e-education.psu.edu/geog489/l1.html' 

html = requests.get(documentURL).text 
soup = BeautifulSoup(html, 'html.parser') 

print(soup.get_text())

This brief code example is similar to what you already saw in Lesson 2. It uses the requests Python package to read in the content of an html page from the URL that is provided in the documentURL variable. Then the package BeautifulSoup4 (bs4) is used for parsing the content of the file and we simply use it to print out the plain text content with all tags and other elements removed by invoking its get_text() method in the last line.

While Jupyter by default is configured to periodically autosave the notebook, this would be a good point to explicitly save the notebook with the newly added content. You can do this by clicking the disk symbol or simply pressing ‘S’ while in ‘Command mode’. The time of the last save will be shown at the top of the document, right next to the notebook name. You can always revert back to the last previously saved version (also referred to as a ‘Checkpoint’ in Jupyter) using File -> Revert to Checkpoint. Undo with CTRL-Z works as expected for the content of a cell while in ‘Edit mode’; however, you cannot use it to undo changes made to the structure of the notebook such as moving cells around. A deleted cell can be recovered by pressing ‘Z’ while in ‘Command mode’ though.

Now that we have a cell with some Python code in our notebook, it is time to execute the code and show the output it produces in the notebook. For this you simply have to click the run symbol button or press ‘SHIFT+Enter’ while in ‘Command mode’. This will execute the currently active cell, place the produced output below the cell, and activate the next cell in the notebook. If there is no next cell (like in our example so far), a new cell will be created. While the code of the cell is being executed, a * will appear within the squared brackets of the “In [ ]:”. Once the execution has terminated, the * will be replaced by a number that always increases by one with each cell execution. This allows for keeping track of the order in which the cells in the notebook have been executed.

Figure 3.9 below shows how things should look after you executed the code cell. The output produced by the print statement is shown below the code in a text field with a vertical scrollbar. We will later see that Jupyter provides the means to display other output than just text, such as images or even interactive maps.

Figure 3.9 Notebook with output produced by running the cell with the code example

In addition to running just a single cell, there are also options for running all cells in the notebook from beginning to end (Cell -> Run All) or for running all cells from the currently activated one until the end of the notebook (Cell -> Run All Below). The produced output is saved as part of the notebook file, so it will be immediately available when you open the notebook again. You can remove the output for the currently active cell by using Cell -> Current Outputs -> Clear, or of all cells via Cell -> All Output -> Clear.

Let’s now put in some heading and information text into the first cell using the Markdown notation. Markdown is a notation and corresponding conversion tool that allows you to create formatted HTML without having to fiddle with tags and with far less typing required. You see examples of how it works by going Help -> Markdown in the menu bar and then clicking the “Basic writing and formatting syntax” link on the web page that opens up. This page here also provides a very brief overview on the markdown notation. If you browse through the examples, you will see that a first level heading can be produced by starting the line with a hashmark symbol (#). To make some text appear in italics, you can delimit it by * symbols (e.g., *text*), and to make it appear in bold, you would use **text** . A simple bullet point list can be produced by a sequence of lines that start with a – or a *.

Let’s say we just want to provide a title and some bullet point list of what is happening in this code example. Click on the text field of the first cell and then type in:

# Simple example of reading a web page and converting it to plain text 
How the code works: 
* package **requests** is used to load web page from URL given in variable *documentURL* 
* package **BeautifulSoup4 (bs4)** is used to parse content of loaded web page 
* the call of *soup.get_text()* in the last line provides the content of page as plain text 

While typing this in, you will notice that Jupyter already interprets the styling information we are providing with the different notations, e.g. by using a larger blue font for the heading, by using bold font for the text appearing within the **…**, etc. However, to really turn the content into styled text, you will need to ‘run the cell’ (SHIFT+Enter) like you did with the code cell. As a result, you should get the nicely formatted text shown in Figure 3.10 below that depicts our entire first Jupyter notebook with text cell, code cell, and output. If you want to see the Markdown code and edit it again, you will have to double-click the text field or press ‘Enter’ to switch to ‘Edit mode’.

Figure 3.10 Notebook with styled text explanation produced with Markdown

If you have not worked with Markdown styling before, we highly recommend that you take a moment to further explore the different styling options from the “Basic writing and formatting syntax” web page. Either use the first cell of our notebook to try out the different notations or create a new Markdown cell at the bottom of the notebook for experimenting.

This little example only covered the main Jupyter operations needed to create a first Jupyter notebook and run the code in it. The ‘Edit’ menu contains many operations that will be useful when creating more complex notebooks, such as deleting, copying, and moving of cells, splitting and merging functionality, etc. For most of these operations, there also exist keyboard shortcuts. If you find yourself in a situation in which you can’t figure out how to use any of these operations, please feel free to ask on the forums.

Jupyter provides a number of so-called magic commands that can be used in code cells to simplify common tasks. Magic commands are interpreted by Jupyter and, for instance, transformed into Python code before the content is passed on to the kernel for execution. This happens behind the scenes, so you will always only see the magic command in your notebook. Magic commands start with a single % symbol if they are line-oriented meaning they should be applied to the remaining content of the line, and with %% if they are cell-oriented meaning they should be applied to the rest of the cell. As a first example, you can use the magic command %lsmagic to list the available magic commands (Figure 3.11). To get the output you have to execute the cell as with any other code cell.

Figure 3.11 Output produced by the magic command %lsmagic that lists all availble magic commands

The %load_ext magic command can be used for loading IPython extension which can add new magic commands. The following command loads the IPython rpy2 extension. If that code gives you a long list of errors then the rpy2 package isn't installed and you will need to go back to Section 3.2 and follow the instructions there.

We recently had cases where loading rpy2 failed on some systems due to the R_HOME environment variable not being set correctly. We therefore added the first line below which you will have to adapt to point to the lib\R folder in your AC Python environment.

import os, rpy2
os.environ['R_HOME'] = r'C:\Users\username\anaconda3\envs\AC37\lib\R' # workaround for R.dll issue occurring on some systems
%load_ext rpy2.ipython

Using a ? symbol in front of a magic command will open a subwindow with the documentation of that command at the bottom of the browser window. Give it a try by executing the command

     ?%R 

?BeautifulSoup

or

?soup.get_text()

Give it a try and see if you understand what the documentation is telling you.

Jupyter notebooks can also include interactive elements, referred to as widgets as in Lesson 2, like buttons, text input fields, sliders, and other GUI elements, as well as visualizations, plots, and animations. Figure 3.12 shows an example that places three button widgets and then simply prints out which button has been pressed when you click on them. The ipywidgets and IPython.display packages imported at the beginning are the main packages required to place the widgets in the notebook. We then define a function that will be invoked whenever one of the buttons is clicked. It simply prints out the description attribute of the button (b.description). In the for-loop we create the three buttons and register the onButtonClick function as the on_click event handler function for all of them. 

from ipywidgets import widgets 
from IPython.display import display 

def onButtonClick(b):
    print("Button " + b.description + " has been clicked")

for i in range(1,4):
    button = widgets.Button(description=str(i))
    display(button)
    button.on_click(onButtonClick)

Figure 3.12 Notebook example using three button widgets and an event handler function that prints out which button has been clicked

If you get an error with this code "Failed to display Jupyter Widget of type Button" that means the widgets are probably not installed which we can potentially fix in our Anaconda prompt:

conda install -n base -c conda-forge widgetsnbextension
conda install -n AC37 -c conda-forge ipywidgets

After installing the packages, exit your Jupyter notebook and restart it and try to re-run your code. It's possible you will receive the error again as the widget tries to run before the Javascript library that runs the widgets has opened. In that case try to select your code, wait a few more seconds and then click Run.

If you're still getting an error, it's likely that your packages didn't install properly (or in a way that Jupyter/Anaconda could find them). The fix for this is to close Jupyter Notebook, return to Anaconda Navigator, click Environments (on the left), choose your environment and then search for "ipy", you may need to either change the "Installed" dropdown to "Not Installed" if they are missing or perhaps they should be updated (by clicking on the upward point arrow or the blue text).

Figure 3.13 Anaconda Navigator showing how to install / update packages

It is easy to imagine how this example could be extended to provide some choices on how the next analysis step in a longer Data Science project should be performed. Similarly, a slider or text input field could be used to allow the notebook user to change the values of important input variables.

Let’s close this brief introduction to Jupyter with a few more things that are good to know when starting to write longer and more complex notebooks. Like normal development environments, Juypter has an autocomplete feature that helps with the code writing and can save a lot of typing: while editing code in a code cell, you can press the TAB key and Jupyter will either automatically complement the name or keyword you are writing or provide you with a dropdown list of choices that you can pick from. For instance, type in soup.ge and then press TAB and you get the list of options, as shown in Figure 3.13 including the get_text() function that we used in our code.

Figure 3.13 Autocompletion of code by pressing TAB

Another useful keyboard command to remember is SHIFT+TAB. When you place the cursor on a variable name or function call and press this key combination, a window will open up showing helpful information like the type of the variable and its current content or the parameters of the function as in Figure 3.14. This is of great help if you are unsure about the different parameters of a function, their order or names. Try out what you get when you use this key combination for different parts of the code in this notebook.

Figure 3.14 Help text provided for a function when pressing SHIFT+TAB

As in all programming, it may occasionally happen that something goes completely wrong and the execution of the code in a cell won’t terminate in the expected time or not at all. If that happens, the first thing to try is to use the “Interrupt Kernel” button located to the right of the “Execute cell” button. This should stop the execution of the current cell and you can then modify the code and try to run the cell again. However, sometimes even that won’t help because the kernel has become unresponsive. In that case, the only thing you can do is restart the kernel using the “Restart Kernel” button to the right of the “Interrupt Kernel” button. Unfortunately, that means that you will have to start the execution of the code in your notebook from the beginning because all imports and variable assignments will be lost after the restart.

Once you have finished your notebook, you may want to publish or share it. There are many options by wich to do so. In the File menu, there exists the “Download as…” option for obtaining versions of your notebook in different formats. The .ipynb format, as we mentioned, is the native format in which Jupyter saves the notebooks. If you make the .ipynb file available to someone who has access to Juypter, that person can open the notebook and run it or modify the content. The .py option allows for exporting content as a Python script, so that the code can be run outside of Jupyter. If you want a version of your notebook that others can read even without access to Jupyter, there are several options like exporting the notebook as HTML, as Latex, or as PDF. Some of these options require additional software tools to be installed on your computer and there are some limitations. For instance, if you export your notebook as HTML to add it to your personal web page, interactive widgets such as the interactive web map you will see later in Section 3.10 will not be included.

To close this section, we want to again refer to the links we provided at the beginning of Section 3.6.2 if you want to keep reading about Jupyter and learn tricks that we weren't able to cover in this brief introduction. In the remainder of this lesson, please use Jupyter to try out the code examples by entering them into a Jupyter notebook and running the code there to get some more practice with Jupyter.

In our examples, we will be using pandas in combination with the numpy package, the package that provides many fundamental scientific computing functionalities for Python and that many other packages are built on. So we start by importing both packages:

import pandas as pd
import numpy as np

A data frame consists of cells that are organized into rows and columns. The rows and columns have names that serve as indices to access the data in the cells. Let us start by creating a data frame with some random numbers that simulate a time series of different measurements (columns) taken on consecutive days (rows) from January 1, 2017 to January 6, 2017. The first step is to create a pandas series of dates that will serve as the row names for our data frame. For this, we use the pandas function date_range(…):

dates = pd.date_range('20170101' , periods=6, freq='D')
dates

The first parameter given to date_range is the starting date. The ‘periods’ parameter tells the function how many dates we want to generate, while we use ‘freq’ to tell it that we want a date for every consecutive day. If you look at the output from the second line we included, you will see that the object returned by the function is a DatetimeIndex object which is a special class defined in pandas.

Next, we generate random numbers that should make up the content of our date frame with the help of the numpy function randn(…) for creating a set of random numbers that follow a standard normal distribution:

numbers = np.random.randn(6,4)
numbers

The output is a two-dimensional numpy array of random numbers normally distributed around 0 with 4 columns and 6 rows. We create a pandas data frame object from it with the following code:

df = pd.DataFrame(numbers, index=dates, columns=['m1', 'm2', 'm3', 'm4'])
df

Note that we provide our array of random numbers as the first parameter, followed by the DatetimeIndex object we created earlier for the row index. For the columns, we simply provide a list of the names with ‘m1’ standing for measurement 1, ‘m2’ standing for measurement 2, and so on. Please also note how the resulting data frame is displayed as a nicely formatted table in your Jupyter Notebook because it makes use of IPython widgets. Please keep in mind that because we are using random numbers for the content of the cells, the output produced by commands used in the following examples will look different in your notebook because the numbers are different.

Now that we have a data frame object available, let’s quickly go through some of the basic operations that one can perform on a data frame to access or modify the data.

The info() method can be used to display some basic information about a data frame such as the number of rows and columns and data types of the columns:

df.info()

The output tells us that we have four columns, all for storing floating point numbers, and each column has 6 rows with values that are not null. If you ever need the number of rows and columns, you can get them by applying the len(…) operation to the data frame and to the columns property of the data frame:

print(len(df))          # gives you the number of rows
print(len(df.columns))  # gives you the number of columns

We can use the head(…) and tail(…) methods to get only the first or last n rows of a data frame:

firstRows = df.head(2)
print(firstRows)
lastRows = df.tail(2)
print(lastRows)

We can also just get a subset of consecutive rows by applying slicing to the data frame similar to how this can be done with lists or strings:

someRows = df[3:5]    # gives you the 4th and 5th row
print(someRows)

This operation gives us rows 4 and 5 (those with index 3 and 4) from the original data frame because the second number used is the index of the first row that will not be included anymore.

If we are just interested in a single column, we can get it based on its name:

thirdColumn = df.m3
print(thirdColumn)

The same can be achieved by using the notation df['m3'] instead of df.m3 in the first line of code. Moreover, instead of using a single column name, you can use a list of column names to get a data frame with just these columns and in the specified order:

columnsM3andM2 = df[ ['m3', 'm2'] ]
columnsM3andM2

The column subsetting and row slicing operations shown above can be concatenated into a single expression. For instance, the following command gives us columns ‘m3’ and ‘m2’ and only the rows with index 3 and 4:

someSubset = df[['m3', 'm2']][3:5]
someSubset

The order here doesn’t matter. We could also have written df[3:5][['m3', 'm2']] .

The most flexible methods for creating subsets of data frame are via the loc and .iloc index properties of a data frame. .iloc[…] is based on the numerical indices of the columns and rows. Here is an example:

someSubset = df.iloc[2:4,1:3] 
print(someSubset)

The part before the comma determines the rows (rows with indices 2 and 3 in this case) and the part after the comma, the columns (columns with indices 1 and 2 in this case). So we get a data frame with the 3rd and 4th rows and 2nd and 3rd columns of the original data frame. Instead of slices we can also use lists of row and column indices to create completely arbitrary subsets. For instance, using iloc in the following example

someSubset = df.iloc[ [0,2,4], [1,3] ]
print(someSubset)

...gives us a data frame with the 1st, 3rd, and 5th row and 2nd and 4th column of the original dataframe. Both the part before the comma and after the comma can just be a colon symbol (:) in which case all rows/columns will be included. For instance,

allRowsSomeColumns = df.iloc[ : , [1,3] ]
print(allRowsSomeColumns)

...will give you all rows but only the 2nd of 4th column.

In contrast to iloc, loc doesn’t use the row and column numbers but instead is based on their labels, while otherwise working in the same way as iloc. The following command produces the same subset of the 1st, 3rd, and 5th rows and 2nd and 4th columns as the iloc code from two examples above:

someSubset = df.loc[ [pd.Timestamp('2017-01-01'), pd.Timestamp('2017-01-03'), pd.Timestamp('2017-01-05')] , ['m2', 'm4'] ]
print(someSubset)

Please note that, in this example, the list for the column names at the very end is simply a list of strings but the list of dates for the row names has to consist of pandas Timestamp objects. That is because we used a DatetimeIndex for the rows when we created the original data frame. When a data frame is displayed, the row names show up as simple strings but they are actually Timestamp objects. However, a DatetimeIndex for the rows has many advantages; for instance, we can use it to get all rows for dates that are from a particular year, e.g. with

df.loc['2017' , ['m2', 'm4'] ]

...to get all dates from 2017 which, of course, in this case, are all rows. Without going into further detail here, we can also get all dates from a specified time period, etc.

The methods explained above for accessing subsets of a data frame can also be used as part of an assignment to change the values in one or several cells. In the simplest case, we can change the value in a single cell, for instance with

df.iloc[0,0] = 0.17 

...or

df.loc['2017-01-01', 'm1'] = 0.17 

...to change the value of the top left cell to 0.17. Please note that this operation will change the original data frame, not create a modified copy. So if you now print out the data frame with

df 

you will see the modified value for 'm1' of January 1, 2017. Even more importantly, if you have used the operations explained above for creating a subset, the data frame with the subset will still refer to the original data, so changing a cell value will change your original data. If you ever want to make changes but keep the original data frame unchanged, you need to explicitly make a copy of the data frame by calling the copy() method as in the following example:

dfCopy = df.copy()
dfCopy.iloc[0,0] = 0
print(df)
print(dfCopy)

Check out the output and compare the top left value for both data frames. The data frame in df still has the old value of 0.17, while the value will be changed to 0 in dfCopy. Without using the copy() method in the first line, both variables would still refer to the same underlying data and both would show the 0 value. Here is another slightly more complicated example where we change the values for the first column of the 1st and 5th rows to 1.2 (intentionally modifying the original data):

df.iloc[ [0,4] , [0] ] = 1.2
print(df)

If you ever need to explicitly go through the rows in a data frame, you can do this with a for-loop that uses the itertuples(…) method of the data frame. itertuples(…) gives you an object that can be used to iterate through the rows as named tuples, meaning each element in the tuple is labeled with the respective column name. By providing the parameter index=False to the method, we are saying that we don’t want the row name to be part of the tuple, just the cell values for the different columns. You can access the elements of the tuple either via their index or via the column name:

for row in df.itertuples(index=False):
    print(row)     # print entire row tuple
    print(row[0])  # print value from column with index 0
    print(row.m2)  # print value from column with name m2
    print('----------')

Try out this example and have a look at the named tuple and the output produced by the other two print statements.

Pandas also provides operations for sorting the rows in a data frame. The following command can be used to sort our data frame by the values in the ‘m2’ column in decreasing order:

dfSorted = df.sort_values(by='m2', ascending=False)
dfSorted

The ‘by’ argument specifies the column that the sorting should be based on and, by setting the ‘ascending’ argument to False, we are saying that we want the rows to be sorted in descending rather than ascending order. It is also possible to provide a list of column names for the ‘by’ argument, to sort by multiple columns. The sort_values(...) method will create a new copy of the data frame, so modifying any cells of dfSorted in this example will not have any impact on the data frame in variable df.

Adding a new column to a data frame is very simple when you have the values for that column ready in a list. For instance, in the following example, we want to add a new column ‘m5’ with additional measurements and we already have the numbers stored in a list m5values that is defined in the first line of the example code. To add the column, we then simply make an assignment to df['m5'] in the second line. If a column ‘m5’ would already exist, its values would now be overwritten by the values from m5values. But since this is not the case, a new column gets added under the name ‘m5’ with the values from m5values.

m5values = [0.432523, -0.123223, -0.231232, 0.001231, -0.23698, -0.41231]
df['m5'] = m5values
df

For adding new rows, we can simply make assignments to the rows selected via the loc operation, e.g. we could add a new row for January 7, 2017 by writing

df.loc[pd.Timestamp('2017-01-07'),:] = [ ... ]

where the part after the equal sign is a list of five numbers, one for each of the columns. Again, this would replace the values in the case that there already is a row for January 7. The following example uses this idea to create new rows for January 7 to 9 using a for loop:

for i in range(7,10):
    df.loc[ pd.Timestamp('2017-01-0'+str(i)),:] = [ np.random.rand() for j in range(5) ]
df

In the body of the for loop, the part on the left of the equal sign uses loc(...) to refer to a row for the new date based on loop variable i, while the part on the right side simply uses the numpy rand() method inside a list comprehension to create a list of five random numbers that will be assigned to the cells of the new row.

If you ever want to remove columns or rows from a data frame, you can do so by using df.drop(...). The first parameter given to drop(...) is a single column or row name or, alternatively, a list of names that should be dropped. By default, drop(...) will consider these as row names. To indicate these are column names that should be removed, you have to specify the additional keyword argument axis=1 . We will see an example of this in a moment.

The following short example demonstrates how pandas can be used to merge two data frames based on a common key, so to perform a join operation in database terms. Let’s say we have two tables, one listing capitals of states in the U.S. and one listing populations for each state. For simplicity we just define data frames for these with just entries for two states, Washington and Oregon:

df1 = pd.DataFrame( {'state': ['Washington', 'Oregon'], 'capital': ['Olympia', 'Salem']} )
print(df1)
df2 = pd.DataFrame( {'name': ['Washington', 'Oregon'], 'population':[7288000, 4093000]} )
print(df2)

The two data frames produced by this code look like this:

We here use a new way of creating a data frame, namely from a dictionary that has entries for each of the columns where the keys are the column names (‘state’ and ‘capital’ in the case of df1, and ‘name’ and ‘population’ in case of df2) and the values stored are lists of the values for that column. We now invoke the merge(...) method on df1 and provide df2 as the first parameter meaning that a new data frame will be produced by merging df1 and df2. We further have to specify which columns should be used as keys for the join operation. Since the two columns containing the state names are called differently, we have to provide the name for df1 through the ‘left_on’ argument and the name for df2 through the ‘right_on’ argument.

merged = df1.merge(df2, left_on='state', right_on='name')
merged

The joined data frame produced by the code will look like this:

As you see, the data frames have been merged correctly. However, we do not want two columns with the state names, so, as a last step, we remove the column called ‘name’ with the previously mentioned drop(...) method. As explained, we have to use the keyword argument axis=1 to indicate that we want to drop a column, not a row.

newMerged = merged.drop('name', axis=1)
newMerged

Result:

If you print out variable merged, you will see that it still contains the 'name' column. That is because drop(...) doesn't change the original data frame but rather produces a copy with the column/row removed.

When working with tabular data, it is very common that one wants to do something with particular data entries that satisfy a certain condition. For instance, we may want to restrict our analysis to rows that have a value larger than a given threshold for one of the columns. Pandas provides some powerful methods for this kind of filtering, and we are going to show one of these to you in this section, namely filtering with Boolean expressions.

The first important thing to know is that we can use data frames in comparison expressions, like df > 0, df.m1 * 2 < 0.2, and so on. The output will be a data frame that only contains Boolean values (True or False) indicating whether the corresponding cell values satisfy the comparison expression or not. Let’s try out these two examples:

df > 0

The result is a data frame with the same rows and columns as the original data frame in df with all cells that had a value larger than 0 set to True, while all other cells are set to False:

df.m1 * 2 < 0.2

Here we are doing pretty much the same thing but only for a single column (‘m1’) and the comparison expression is slightly more complex involving multiplication of the cell values with 2 before the result is compared to 0.2. The result is a one-dimensional vector of True and False values corresponding to the cells of the ‘m1’ column in the original data frame:

Just to introduce another useful pandas method, we can apply the value_counts() method to get a summary of the values in a data frame telling how often each value occurs:

(df.m1 * 2 < 0.2).value_counts()

The expression in the parentheses will give us a boolean column vector as we have seen above, and invoking its value_counts() method tells us how often True and False occur in this vector. (The actual numbers will depend on the random numbers in your original data frame).

The second important component of Boolean indexing is that we can use Boolean operators to combine Boolean data frames as illustrated in the next example:

v1 = df.m1 * 2 < 0.2
print(v1)
v2 = df.m2 > 0
print(v2)
print(~v1)
print(v1 & v2)
print(v1 | v2)

This will produce the following output data frames:

What is happening here? We first create two different Boolean vectors using two different comparison expressions for columns ‘m1’ and ‘m2’, respectively, and store the results in variables v1 and v2. Then we use the Boolean operators ~ (not), & (and), and | (or) to create new Boolean vectors from the original ones, first by negating the Boolean values from v1, then by taking the logical AND of the corresponding values in v1 and v2 (meaning only cells that have True for both v1 and v2 will be set to True in the resulting vector), and finally by doing the same but with the logical OR (meaning only cells that have False for both v1 and v2 will be set to False in the result). We can construct arbitrarily complex Boolean expressions over the values in one or multiple data frames in this way.

The final important component is that we can use Boolean vectors or lists to select rows from a data frame. For instance,

df[ [True, False, True, False, True, False, True, False, True] ]

... will give us a subset of the data frame with only every second row:

Taken these three things together means we can use arbitrarily logical expressions over the values in a data frame to select a subset of rows that we want to work with. To continue the examples from above, let’s say that we want only those rows that satisfy both the criteria df.m1 * 2 < 0.2 and df.m2 > 0, so only those rows for which the value of column ‘m1’ times 2 is smaller than 0.2 and the value of column ‘m2’ is larger than 0. We can use the following expression for this:

df[v1 & v2 ]

Or even without first having to define v1 and v2:

df[ (df.m1 * 2 < 0.2)  & (df.m2 > 0)  ]

Here is the resulting data frame:

Hopefully you are beginning to see how powerful this approach is and how it allows for writing very elegant and compact code for working with tabular data. You will get to see more examples of this and of using pandas in general in the lesson’s walkthrough. There we will also be using GeoPandas, an extension built on top of pandas that allows for working with data frames that contain geometry data, e.g. entire attribute tables of an ESRI shapefile.

OGR and GDAL exist as separate modules in the osgeo package together with some other modules such as OSR for dealing with different projections and some auxiliary modules. As with previous packages you will probably need to install gdal (which encapsulates all of them in its latest release) to access these packages. So typically, a Python project using GDAL would import the needed modules similar to this:

from osgeo import gdal
from osgeo import ogr
from osgeo import osr

Let’s start with taking a closer look at the ogr module for working with vector data. The following code, for instance, illustrates how one would open an existing vector data set, in this case an Esri shapefile. OGR uses so-called drivers to access files in different formats, so the important thing to note is how we first obtain a driver object for ‘ESRI Shapefile’ with GetDriverByName(...) and then use that to open the shapefile with the Open(...) method of the driver object. The shapefile we are using in this example is a file with polygons for all countries in the world (available here) and we will use it again in the lesson’s walkthrough. When you download it, you may still have to adapt the path in the first line of the code below.

shapefile = r'C:\489\L3\TM_WORLD_BORDERS-0.3.shp'
drv =ogr.GetDriverByName('ESRI Shapefile')
dataSet = drv.Open(shapefile)

dataSet now provides access to the data in the shapefile as layers, in this case just a single layer, that can be accessed with the GetLayer(...) method. We then use the resulting layer object to get the definitions of the fields with GetLayerDefn(), loop through the fields with the help of GetFieldCount() and GetFieldDefn(), and then print out the field names with GetName():

layer = dataSet.GetLayer(0)
layerDefinition = layer.GetLayerDefn()
for i in range(layerDefinition.GetFieldCount()):
	print(layerDefinition.GetFieldDefn(i).GetName())

Output:
FIPS
ISO2
ISO3
UN
NAME
...
LAT

If you want to loop through the features in a layer, e.g., to read or modify their attribute values, you can use a simple for-loop and the GetField(...) method of features. It’s important that, if you want to be able to iterate through the features another time, you have to call the ResetReading() method of the layer after the loop. The following loop prints out the values of the ‘NAME’ field for all features:

for feature in layer:
	print(feature.GetField('NAME'))
layer.ResetReading()

Output:
Antigua and Barbuda
Algeria
Azerbaijan
Albania
....

We can extend the previous example to also access the geometry of the feature via the GetGeometryRef() method. Here we use this approach to get the centroid of each country polygon (method Centroid()) and then print it out in readable form with the help of the ExportToWkt() method. The output will be the same list of country names as before but this time, each followed by the coordinates of the respective centroid.

for feature in layer:
    print(feature.GetField('NAME') + '–' + feature.GetGeometryRef().Centroid().ExportToWkt())
layer.ResetReading()

We can also filter vector layers by attribute and spatially. The following example uses SetAttributeFilter(...) to only get features with a population (field ‘POP2005’) larger than one hundred million.

layer.SetAttributeFilter('POP2005 > 100000000')
for feature in layer:
	print(feature.GetField('NAME'))
layer.ResetReading()

Output:
Brazil
China
India
...

We can remove the selection by calling SetAttributeFilter(...) with the empty string:

layer.SetAttributeFilter('')

The following example uses SetSpatialFilter(...) to only get countries overlapping a given polygonal area, in this case an area that covers the southern part of Africa. We first create the polygon via the CreateGeometryFromWkt(...) function that creates geometry objects from Well-Known Text (WKT) strings. Then we apply the filter and use a for-loop again to print the selected features:

wkt = 'POLYGON ( (6.3 -14, 52 -14, 52 -40, 6.3 -40, 6.3 -14))'  # WKT polygon string for a rectangular area
geom = ogr.CreateGeometryFromWkt(wkt) 
layer.SetSpatialFilter(geom)
for feature in layer:
	print(feature.GetField('NAME'))
layer.ResetReading()

Output:
Angola
Madagascar
Mozambique
...
Zimbabwe

Access to all features and their geometries together with the geometric operations provided by GDAL allows for writing code that realizes geoprocessing operations and other analyses on one or several input files and for creating new output files with the results. To show one example, the following code takes our selection of countries from southern Africa and then creates 2 decimal degree buffers around the centroids of each country. The resulting buffers are written to a new shapefile called centroidbuffers.shp. We add two fields to the newly produced buffered centroids shapefile, one with the name of the country and one with the population, copying over the corresponding field values from the input country file. When you will later have to use GDAL in one of the parts of the lesson's homework assignment, you can use the same order of operations, just with different parameters and input values.

To achieve this task, we first create a spatial reference object for EPSG:4326 that will be needed to create the layer in the new shapefile. Then the shapefile is generated with the shapefile driver object we obtained earlier, and the CreateDataSource(...) method. A new layer inside this shapefile is created via CreateLayer(...) and by providing the spatial reference object as a parameter. We then create the two fields for storing the country names and population numbers with the function FieldDefn(...) and add them to the created layer with the CreateField(...) method using the field objects as parameters. When adding new fields, make sure that the name you provide is not longer than 8 characters or GDAL/OGR will automatically truncate the name. Finally, we need the field definitions of the layer available via GetLayerDefn() to be able to later add new features to the output shapefile. The result is stored in variable featureDefn.

sr = osr.SpatialReference()   # create spatial reference object
sr.ImportFromEPSG(4326)       # set it to EPSG:4326
outfile = drv.CreateDataSource(r'C:\489\L3\centroidbuffers.shp') # create new shapefile
outlayer = outfile.CreateLayer('centroidbuffers', geom_type=ogr.wkbPolygon, srs=sr)  # create new layer in the shapefile 

nameField = ogr.FieldDefn('Name', ogr.OFTString)        # create new field of type string called Name to store the country names
outlayer.CreateField(nameField)                         # add this new field to the output layer
popField = ogr.FieldDefn('Population', ogr.OFTInteger) # create new field of type integer called Population to store the population numbers
outlayer.CreateField(popField)                         # add this new field to the output layer

featureDefn = outlayer.GetLayerDefn()  # get field definitions

Now that we have prepared the new output shapefile and layer, we can loop through our selected features in the input layer in variable layer, get the geometry of each feature, produce a new geometry by taking its centroid and then calling the Buffer(...) method, and finally create a new feature from the resulting geometry within our output layer.

for feature in layer:                                              # loop through selected features
    ingeom = feature.GetGeometryRef()                              # get geometry of feature from the input layer
    outgeom = ingeom.Centroid().Buffer(2.0)                        # buffer centroid of ingeom

    outFeature = ogr.Feature(featureDefn)                          # create a new feature
    outFeature.SetGeometry(outgeom)                                # set its geometry to outgeom
    outFeature.SetField('Name', feature.GetField('NAME'))          # set the feature's Name field to the NAME value of the input feature
    outFeature.SetField('Population', feature.GetField('POP2005')) # set the feature's Population field to the POP2005 value of the input feature 
    outlayer.CreateFeature(outFeature)                             # finally add the new output feature to outlayer
    outFeature = None

layer.ResetReading()
outfile = None         # close output file

The final line “outfile = None” is needed because otherwise the file would remain open for further writing and we couldn’t inspect it in a different program. If you open the centroidbuffers.shp shapefile and the country borders shapefile in a GIS of your choice, the result should look similar to the image below. If you check out the attribute table, you should see the Name and Populations columns we created populated with values copied over from the input features. 

Figure 3.16 Shapefile with buffered centroids produced by our example code overlaid on the world borders shapefile in a GIS

Centroid() and Buffer() are just two examples of the many availabe methods for producing a new geometry in GDAL. In this lesson's homework assignment, you will instead have to use the ogr.CreateGeometryFromWkt(...) function that we used earlier in this section to create a clip polygon from a WKT string representation but, apart from that, the order of operations for creating the output feature will be the same. The GDAL/OGR cookbook contains many Python examples for working with vector data with GDAL, including how to create different kinds of geometries from different input formats, calculating envelopes, lengths and areas of a geometry, and intersecting and combining geometries. We recommend that you take a bit of time to browse through these online examples to get a better idea of what is possible. Then we move on to have a look at raster manipulation with GDAL.

To open an existing raster file in GDAL, you would use the Open(...) function defined in the gdal module. The raster file we will use in the following examples contains world-wide bioclimatic data and will be used again in the lesson’s walkthrough. Download the raster file here.

raster = gdal.Open(r'c:\489\L3\wc2.0_bio_10m_01.tif')

We now have a GDAL raster dataset in variable raster. Raster datasets are organized into bands. The following command shows that our raster only has a single band:

raster.RasterCount

Output: 
1

To access one of the bands, we can use the GetRasterBand(...) method of a raster dataset and provide the number of the band as a parameter (counting from 1, not from 0!):

band = raster.GetRasterBand(1)  # get first band of the raster

If your raster has multiple bands and you want to perform the same operations for each band, you would typically use a for-loop to go through the bands:

for i in range(1, raster.RasterCount + 1):
    b = raster.GetRasterBand(i)
    print(b.GetMetadata())      # or do something else with b

There are a number of methods to read different properties of a band in addition to GetMetadata() used in the previous example, such as GetNoDataValue(), GetMinimum(), GetMaximum(), andGetScale().

print(band.GetNoDataValue())
print(band.GetMinimum())
print(band.GetMaximum())
print(band.GetScale())

Output:
-1.7e+308
-53.702073097229
33.260635217031
1.0

GDAL provides a number of operations that can be employed to create new files from bands. For instance, the gdal.Polygonize(...) function can be used to create a vector dataset from a raster band by forming polygons from adjacent cells that have the same value. To apply the function, we first create a new vector dataset and layer in it. Then we add a new field ‘DN’ to the layer for storing the raster values for each of the polygons created:

drv = ogr.GetDriverByName('ESRI Shapefile')
outfile = drv.CreateDataSource(r'c:\489\L3\polygonizedRaster.shp') 
outlayer = outfile.CreateLayer('polygonized raster', srs = None )
newField = ogr.FieldDefn('DN', ogr.OFTReal)
outlayer.CreateField(newField)

Once the shapefile is prepared, we call Polygonize(...) and provide the band and the output layer as parameters plus a few additional parameters needed:

gdal.Polygonize(band, None, outlayer, 0, [])
outfile = None

With the None for the second parameter we say that we don’t want to provide a mask for the operation. The 0 for the fourth parameter is the index of the field to which the raster values shall be written, so the index of the newly added ‘DN’ field in this case. The last parameter allows for passing additional options to the function but we do not make use of this, so we provide an empty list. The second line "outfile = None" is for closing the new shapefile and making sure that all data has been written to it. The result produced in the new shapefile rasterPolygonized.shp should look similar to the image below when looked at in a GIS and using a classified symbology based on the values in the ‘DN’ field.

Figure 3.17 Polygonized raster file produced by the previous code example shown in a GIS

Polygonize(...) is an example of a GDAL function that operates on an individual band. GDAL also provides functions for manipulating raster files directly, such as gdal.Translate(...) for converting a raster file into a new raster file. Translate(...) is very powerful with many parameters and can be used to clip, resample, and rescale the raster as well as convert the raster into a different file format. You will see an example of Translate(...) being applied in the lesson’s walkthrough. gdal.Warp(...) is another powerful function that can be used for reprojecting and mosaicking raster files.

While the functions mentioned above and similar functions available in GDAL cover many of the standard manipulation and conversion operations commonly used with raster data, there are cases where one directly wants to work with the values in the raster, e.g. by applying raster algebra operations. The approach to do this with GDAL is to first read the data of a band into a GDAL multi-dimensional array object with the ReadAsArray() method, then manipulate the values in the array, and finally write the new values back to the band with the WriteArray() method.

data = band.ReadAsArray()
data

If you look at the output of this code, you will see that the array in variable data essentially contains the values of the raster cells organized into rows. We can now apply a simple mathematical expression to each of the cells, like this:

data = data * 0.1
data

The meaning of this expression is to create a new array by multiplying each cell value with 0.1. You should notice the change in the output from +308 to +307. The following expression can be used to change all values that are smaller than 0 to 0:

print(data.min())
data [ data < 0 ] = 0
print(data.min())

data.min() in the previous example is used to get the minimum value over all cells and show how this changes to 0 after executing the second line. Similarly to what you saw with pandas data frames in Section 3.8.6, an expression like data < 0 results in a Boolean array with True for only those cells for which the condition <0 is true. Then this Boolean array is used to select only specific cells from the array with data[...] and only these will be changed to 0. Now, to finally write the modified values back to a raster band, we can use the WriteArray(...) method. The following code shows how one can first create a copy of a raster with the same properties as the original raster file and then use the modified data to overwrite the band in this new copy:

drv = gdal.GetDriverByName('GTiff')     # create driver for writing geotiff file
outRaster = drv.CreateCopy(r'c:\489\newRaster.tif', raster , 0 )   # create new copy of inut raster on disk
newBand = outRaster.GetRasterBand(1)                               # get the first (and only) band of the new copy
newBand.WriteArray(data)                                           # write array data to this band 
outRaster = None                                                   # write all changes to disk and close file

This approach will not change the original raster file on disk. Instead of writing the updated array to a band of a new file on disk, we can also work with an in-memory copy instead, e.g. to then use this modified band in other GDAL operations such as Polygonize(...) . An example of this approach will be shown in the walkthrough of this lesson. Here is how you would create the in-memory copy combining the driver creation and raster copying into a single line:

tmpRaster = gdal.GetDriverByName('MEM').CreateCopy('', raster, 0) # create in-memory copy

The approach of using raster algebra operations shown above can be used to perform many operations such as reclassification and normalization of a raster. More complex operations like neighborhood/zonal based operators can be implemented by looping through the array and adapting cell values based on the values of adjacent cells. In the lesson’s walkthrough you will get to see an example of how a simple reclassification can be realized using expressions similar to what you saw in this section.

While the GDAL Python package allows for realizing the most common vector and raster operations, it is probably fair to say that it is not the most easy-to-use software API. While the GDAL Python cookbook contains many application examples, it can sometimes take a lot of search on the web to figure out some of the details of how to apply a method or function correctly. Of course, GDAL has the main advantage of being completely free, available for pretty much all main operating systems and programming languages, and not tied to any other GIS or web platform. In contrast, the Esri ArcGIS API for Python discussed in the next section may be more modern, directly developed specifically for Python, and have more to offer in terms of visualization and high-level geoprocessing and analysis functions, but it is tied to Esri’s web platforms and some of the features require an organizational account to be used. These are aspects that need to be considered when making a choice on which API to use for a particular project. In addition, the functionality provided by both APIs only overlaps partially and, as a result, there are also merits in combining the APIs as we will do later on in the lesson’s walkthrough.

The central class of the API is the GIS class defined in the arcgis.gis module. It represents the GIS you are working with to access and publish content or to perform different kinds of management or analysis tasks. A GIS object can be connected to AGOL or Enterprise/Portal for ArcGIS. Typically, the first step of working with the API in a Python script consists of constructing a GIS object by invoking the GIS(...) function defined in the argis.gis module. There are many different ways of how this function can be called, depending on the target GIS and the authentication mechanism that should be employed. Below we are showing a couple of common ways to create a GIS object before settling on the approach that we will use in this class, which is tailored to work with your pennstate organizational account.

Most commonly, the GIS(...) function is called with three parameters, the URL of the ESRI web GIS to use (e.g. ‘https://www.arcgis.com’ for a GIS object that is connected to AGOL), a username, and the login password of that user. If no URL is provided, the default is ArcGIS Online and it is possible to create a GIS object anonymously without providing username and password. So the simplest case would be to use the following code to create an anonymous GIS object connected to AGOL (we do not recommend running this code - it is shown only as an example):

from arcgis.gis import GIS
gis = GIS()

Due to not being connected to a particular user account, the available functionality will be limited. For instance, you won’t be able to upload and publish new content. If, instead, you'd want to create the GIS object but with a username and password for a normal AGOL account (so not an enterprise account), you would use the following way with <your username> and <your password> replaced by the actual user name and password (we do not recommend running this code either- it is shown only as an example):

from arcgis.gis import GIS
gis = GIS('https://www.arcgis.com', '<your username>', '<your password>')
gis?

This approach unfortunately does not work with the pennstate organizational account you had to create at the beginning of the class. But we will need to use this account in this lesson to make sure you have the required permissions. You already saw how to connect with your pennstate organizational account if you tried out the test code in Section 3.2. The URL needs to be changed to ‘https://pennstate.maps.arcgis.com’ and we have to use a parameter called client_id to provide the ID of an app that we created for this class on AGOL. When using this approach, calling the GIS(...) function will open a browser window where you can authenticate with your PSU credentials and/or you will be asked to grant the permission to use Python with your pennstate AGOL account. After that, you will be given a code that you have to paste into a box that will be showing up in your notebook. The details of what is going on behind the scenes are a bit complicated but whenever you need to create a GIS object in this class to work with the API, you can simply use these exact few lines and then follow the steps we just described. The last line of the code below is for testing that the GIS object has been created correctly. It will print out some help text about the GIS object in a window at the bottom of the notebook. The screenshot below illustrates the steps needed to create the GIS object and the output of the gis? command.  (for safety it would be best to restart your Notebook kernel prior to running the code below if you did run the code above as it can cause issues with the below instructions working properly)

import arcgis
from arcgis.gis import GIS
gis = GIS('https://pennstate.maps.arcgis.com', client_id='lDSJ3yfux2gkFBYc')
gis?

Figure 3.18 Steps for creating GIS object and output produced by the gis? command in the previous code example

The GIS object in variable gis gives you access to a user manager (gis.users), a group manager (gis.groups) and a content manager (gis.content) object. The first two are useful for performing administrative tasks related to groups and user management. If you are interested in these topics, please check out the examples on the Accessing and Managing Users and Batch Creation of Groups pages. Here we simply use the get(...) method of the user manager to access information about a user, namely ourselves. Again, you will have to replace the <your username> tag with your PSU AGOL name to get some information display related to your account:

user = gis.users.get('<your username>')
user

The user object now stored in variable user contains much more information about the given user in addition to what is displayed as the output of the previous command. To see a list of available attributes, type in

user.

and then press the TAB key for the Jupyter autocompletion functionality. The following command prints out the email address associated with your account:

user.email

We will talk about the content manager object in a moment. But before that, let’s talk about maps and how easy it is to create a web map like visualization in a Jupyter Notebook with the ArcGIS API. A map can be created by calling the map(...) method of the GIS object. We can pass a string parameter to the method with the name of a place that should be shown on the map. The API will automatically try to geocode that name and then set up the parameters of the map to show the area around that place. Let’s use State College as an example:

map = gis.map('State College, PA')
map

Figure 3.19 Map widget in Jupyter Notebook produced by the previous code example

The API contains a widget for maps so that these will appear as an interactive map in a Jupyter Notebook. As a result, you will be able to pan around and use the buttons at the top left to zoom the map. The map object has many properties and methods that can be used to access its parameters, affect the map display, provide additional interactivity, etc. For instance, the following command changes the zoom property to include a bit more area around State College. The map widget will automatically update accordingly.

map.zoom = 11

With the following command, we change the basemap tiles of our map to satellite imagery. Again the map widget automatically updates to reflect this change.

map.basemap = 'satellite'

As another example, let’s add two simple circle marker objects to the map, one for Walker Building on the PSU campus, the home of the Geography Department, and one for the Old Main building. The ArcGIS API provides a very handy geocode(...) function in the arcgis.geocoding module so that we won’t have to type in the coordinates of these buildings ourselves. The properties of the marker are defined in the pt_sym dictionary.

from arcgis.geocoding import geocode

pt_sym = {
    "type": "esriSMS",
    "style": "esriSMSCircle",
    "color": [255,0,0,255],        
}

walker = geocode('Walker Bldg, State College, PA')[0]
oldMain = geocode('Old Main Bldg, State College, PA')[0]
map.draw(walker, {'title':'Walker Bldg', 'content': walker['attributes']['LongLabel']}, symbol=pt_sym)
map.draw(oldMain, {'title':'Old Main Bldg', 'content': oldMain['attributes']['LongLabel']}, symbol=pt_sym)

Figure 3.20 Map widget in Jupyter Notebook produced by the previous code example

The map widget should now include the markers for the two buildings as shown in the image above (the markers may look different). A little explanation on what happens in the code: the geocode(...) function returns a list of candidate entities matching the given place name or address, with the candidate considered most likely appearing as the first one in the list. We here simply trust that the ArcGIS API has been able to determine the correct candidate and hence just take the first object from the respective lists via the “[0]” in the code. What we have now stored in variables walker and oldMain are dictionaries with many attributes to describe the entities. When adding the markers to the map, we provide the respective variables as the first parameter and the API will automatically access the location information in the dictionaries to figure out where the marker should be placed. The second parameter is a dictionary that specifies what should be shown in the popup that appears when you click on the marker. We use short names for the title of the popups and then use the ‘LongLabel’ property from the dictionaries, which contains the full address, for the content of the popups.

Let's have a look at a last map example demonstrating how one can use the API to manually draw a route from Walker building to Old Main and have the API calculate the length of that route. The first thing we do is import the lengths(...) function from the arcgis.geometry module and define a function that will be called once the route has been drawn to calculate the length of the resulting polyline object that is passed to the function in parameter g:

from arcgis.geometry import lengths

def calcLength(map,g):
    length = lengths(g['spatialReference'], [g], "", "geodesic")
    print('Length: '+ str(length[0]) + 'm.')

Once the length has been computed, the function will print out the result in meters. We now register this function as the function to be called when a drawing action on the map has terminated. In addition, we define the symbology for the drawing to be a dotted red line:

map.on_draw_end(calcLength)

line_sym = {
"type": "esriSLS",
"style": "esriSLSDot",
"color": [255,0,0,255],
"width": 3
}

With the next command we start the drawing with the ‘polyline’ tool. Before you execute the command, it's a good idea to make sure the map widget is zoomed in to show the area between the two markers but you will still be able to pan the map by dragging with the left mouse button pressed and to zoom with the mouse wheel during the drawing. You need to do short clicks with the left mouse button to set the points for the start point and end points of the line segments. To end the polyline, make a double click for the last point while making sure you don't move the mouse at the same time. After this, the calculated distance will appear under the map widget.

map.draw('polyline', symbol=line_sym)

You can always restart the drawing by executing the previous code line again. Using the command map.clear_graphics() allows for clearing all drawings from the map but you will have to recreate the markers after doing so.

Figure 3.21 Map widget produced by the previous code example showing the drawn route

Now let’s get back to the content manager object of our GIS. The content manager object allows us to search and use all the content we have access to. That includes content we have uploaded and published ourselves but also content published within our organization and content that is completely public on AGOL. Content can be searched with the search(...) method of the content manager object. The method takes a query string as parameter and has additional parameters for setting the item type we are looking for and the maximum number of entries returned as a result for the query. For instance, try out the following command to get a list of available feature services that have PA in the title:

featureServicesPA = gis.content.search(query='title:PA', item_type='Feature Layer Collection', max_items = 50)
featureServicesPA

This will display a list of different AGOL feature service items available to you. The list is probably rather short at the moment because not much has been published in the new pennstate organization yet, but there should at least be one entry with municipality polygons for PA. Each feature service from the returned list can, for instance, be added to your map or used for some spatial analysis. To add a feature service to your map as a new layer, you have to use the add_layer(...) method of the map object. For example, the following command takes the first feature service from our result list in variable featureServicesPA and adds it to our map widget:

map.add_layer(featureServicesPA[0], {'opacity': 0.8})

The first parameter specifies what should be added as a layer (the feature service with index 0 from our list in this case), while the second parameter is an optional dictionary with additional attributes (e.g., opacity, title, visibility, symbol, renderer) specifying how the layer should be symbolized and rendered. Feel free to try out a few other queries (e.g. using "*" for the query parameter to get a list of all available feature services) and adding a few other layers to the map by changing the index used in the previous command. If the map should get too cluttered, you can simply recreate the map widget with the map = gis.map(...) command from above.

The query string given to search(...) can contain other search criteria in addition to ‘title’. For instance, the following command lists all feature services that you are the owner of (replace <your agol username> with your actual Penn State AGOL/Pro user name):

gis.content.search(query='owner:<your agol username>', item_type='Feature Service')

Unless you have published feature services in your AGOL account, the list returned by this search will most likely be empty ([]). So how can we add new content to our GIS and publish it as a feature service? That’s what the add(...) method of the content manager object and the publish(...) method of content items are for. The following code uploads and publishes a shapefile of larger cities in the northeast of the United States.

To run it, you will need to first download the zipped shapefile and then slightly change the filename on your computer. Since AGOL organizations must have unique services names, edit the file name to something like ne_cities_[initials].zip or ne_cities_[your psu id].zip so the service name will be unique to the organization. Adapt the path used in the code below to refer to this .zip file on your disk.

neCities = gis.content.add({'type': 'Shapefile'}, r'C:\489\L3\ne_cities.zip')
neCitiesFS = neCities.publish()

The first parameter given to gis.content.add(...) is a dictionary that can be used to define properties of the data set that is being added, for instance tags, what type the data set is, etc. Variable neCitiesFS now refers to a published content item of type arcgis.gis.Item that we can add to a map with add_layer(...). There is a very slight chance that this layer will not load in the map. If this happens for you - continue and you should be able to do the buffer & intersection operations later in the walkthrough without issue. Let’s do this for a new map widget:

cityMap = gis.map('Pennsylvania')
cityMap.add_layer(neCitiesFS, {})
cityMap

Figure 3.22 Map widget produced by the previous code example

The image shows the new map widget that will now be visible in the notebook. If we rerun the query from above again...

gis.content.search(query='owner:<your agol username>', item_type='Feature Service')

...our newly published feature service should now show up in the result list as:

<Item title:"ne_cities" type:Feature Layer Collection owner:<your AGOL username>>

The ArcGIS Python API also provides the functionality to access individual features in a feature layer, to manipulate layers, and to run analyses from a large collection of GIS tools including typical GIS geoprocessing tools. Most of this functionality is available in the different submodules of arcgis.features. To give you a few examples, let’s continue with the city feature service we still have in variable neCitiesFS from the previous section. A feature service can actually have several layers but, in this case, it only contains one layer with the point features for the cities. The following code example accesses this layer (neCitiesFE.layers[0]) and prints out the names of the attribute fields of the layer by looping through the ...properties.fields list of the layer:

for f in neCitiesFS.layers[0].properties.fields:
	print(f)

When you try this out, you should get a list of dictionaries, each describing a field of the layer in detail including name, type, and other information. For instance, this is the description of the STATEABB field:

{
 "name": "STATEABB",
 "type": "esriFieldTypeString",
 "actualType": "nvarchar",
 "alias": "STATEABB",
 "sqlType": "sqlTypeNVarchar",
 "length": 10,
 "nullable": true,
 "editable": true,
 "domain": null,
 "defaultValue": null
 }

The query(...) method allows us to create a subset of the features in a layer using an SQL query string similar to what we use for selection by attribute in ArcGIS itself or with arcpy. The result is a feature set in ESRI terms. If you look at the output of the following command, you will see that this class stores a JSON representation of the features in the result set. Let’s use query(...) to create a feature set of all the cities that are located in Pennsylvania using the query string "STATEABB"='US-PA' for the where parameter of query(...):

paCities = neCitiesFS.layers[0].query(where='"STATEABB"=\'US-PA\'')
print(paCities.features)
paCities

Output: 
[{"geometry": {"x": -8425181.625237303, "y": 5075313.651659228}, 
"attributes": {"FID": 50, "OBJECTID": "149", "UIDENT": 87507, "POPCLASS": 2, "NAME": 
"Scranton", "CAPITAL": -1, "STATEABB": "US-PA", "COUNTRY": "USA"}}, {"geometry": 
{"x": -8912583.489066456, "y": 5176670.443556941}, "attributes": {"FID": 53, 
"OBJECTID": "156", "UIDENT": 88707, "POPCLASS": 3, "NAME": "Erie", "CAPITAL": -1, 
"STATEABB": "US-PA", "COUNTRY": "USA"}}, ...  ]
<FeatureSet> 10 features

Of course, the queries we can use with the query(...) function can be much more complex and logically connect different conditions. The attributes of the features in our result set are stored in a dictionary called attributes for each of the features. The following loop goes through the features in our result set and prints out their name (f.attributes['NAME']) and state (f.attributes['STATEABB']) to verify that we only have cities for Pennsylvania now:

for f in paCities:
	print(f.attributes['NAME'] + " - "+ f.attributes['STATEABB'])

Output: 
Scranton - US-PA
Erie - US-PA
Wilkes-Barre - US-PA
...
Pittsburgh - US-PA

Now, to briefly illustrate some of the geoprocessing functions in the ArcGIS Python API, let’s look at the example of how one would determine the parts of the Appalachian Trail that are within 30 miles of a larger Pennsylvanian city. For this we first need to upload and publish another data set, namely one that’s representing the Appalachian Trail. We use a data set that we acquired from PASDA for this, and you can download the DCNR_apptrail file here. As with the cities layer earlier, there is a very small chance that this will not display for you - but you should be able to continue to perform the buffer and intersection operations. The following code uploads the file to AGOL (don’t forget to adapt the path and add your initials or ID to the name!), publishes it, and adds the resulting trail feature service to your cityMap from above:

appTrail = gis.content.add({'type': 'Shapefile'}, r'C:\489\L3\dcnr_apptrail_2003.zip')
appTrailFS = appTrail.publish()
cityMap.add_layer(appTrailFS, {})

Figure 3.23 Map widget now showing the cities and trail layers

Next, we create a 30 miles buffer around the cities in variable paCities that we can then intersect with the Appalachian Trail layer. The create_buffers(...) function that we will use for this is defined in the arcgis.features.use_proximity module together with other proximity based analysis functions. We provide the feature set we want to buffer as the first parameter, but, since the function cannot be applied to a feature set directly, we have to invoke the to_dict() method of the feature set first. The second parameter is a list of buffer distances allowing for creating multiple buffers at different distances. We only use one distance value, namely 30, here and also specify that the unit is supposed to be ‘Miles’. Finally, we add the resulting buffer feature set to the cityMap above.

from arcgis.features.use_proximity import create_buffers
bufferedCities = create_buffers(paCities.to_dict(), [30], units='Miles')
cityMap.add_layer(bufferedCities, {})

Figure 3.24 Map widget with buffers created by the previous code example

As the last step, we use the overlay_layers(...) function defined in the arcgis.features.manage_data module to create a new feature set by intersecting the buffers with the Appalachian Trail polylines. For this we have to provide the two input sets as parameters and specify that the overlay operation should be ‘Intersect’.

from arcgis.features.manage_data import overlay_layers
trailPartsCloseToCity = overlay_layers(appTrailFS, bufferedCities, overlay_type='Intersect')

We show the result by creating a new map ...

resultMap = gis.map('Pennsylvania')
resultMap

... and just add the features from the resulting trailPartsCloseToCity feature set to this map:

resultMap.add_layer(trailPartsCloseToCity)

The result is shown in the figure below.

Figure 3.25 New map widget showing the final result of our analysis

This section gave you a brief introduction to ESRI's ArcGIS API for Python and showed you some of the most important methods and patterns that you will typically encounter when using the API. The API is much too rich to cover everything here, so, to get a better idea of the geoprocessing and other capabilities of the API, we recommend that you check and try out some of the sample notebooks that ESRI has published, such as:

• Chennai Floods 2015 - a geographic analysis
• Fighting California forest fires using spatial analysis

A lot of valuable information can also be gained from the ESRI's API guide, and, of course, the API documentation. 

import re 

datePattern = re.compile('\d\d\d\d-(0[1-9]|1[0-2])-(0[1-9]|[12]\d|3[01])$') 

def isValidDate(s): 
    return datePattern.match(s) != None 

Explanation: Since we are using match(…) to compare the compiled pattern in variable datePattern to the string in parameter s given to our function isValidDate(…), we don’t have to worry about additional characters before the start of the date because match(…) will always try to match the pattern to the start of the string. However, we use $ as the last character in our pattern to make sure there are no additional characters following the date. That means the pattern has the form

“…-…-…$”

where the dots have to be replaced with some regular expression notation for the year, month, and day parts. The year part is easy, since we allow for any 4-digit number here. So we can use \d\d\d\d here, or alternatively \d{4,4} (remember that \d stands for the predefined class of all digits).

For the month, we need to distinguish two cases: either it is a 0 followed by one of the digits 1-9 (but not another 0) or a 1 followed by one of the digits 0-2. We therefore write this part as a case distinction (…|…) with the left part 0[1-9] representing the first option and the second part 1[0-2] representing the second option.

For the day, we need to distinguish three cases: (1) a 0 followed by one of the digits 1-9, (2) a 1 or 2 followed by any digit, or (3) a 3 followed by a 0 or a 1. Therefore we use a case-distinction with three options (…|…|…) for this part. The first part 0[1-9] is for option (1), the second part [12]\d for option (2), and the third part 3[01] for the third option.

Part 1:

import operator 

from functools import reduce 

l = [True, False, True] 

r = reduce(operator.and_, l, True) 

print(r)  #  output will be False in this case 

To check whether or not at least one element is True, the call has to be changed to:

r = reduce(operator.or_, l, False) 

Part 2:

import operator 

from functools import reduce 

l = [-4, 2, 1, -6 ] 

r = reduce(operator.and_, map(lambda n: n > 0, l), True) 

print(r)   # will print False in this case 

We use map(…) with a lambda expression for checking whether or not an individual element from the list is >0. Then we apply the reduce(…) version from part 1 to the resulting list of Boolean values we get from map(…) to check whether or not all elements are True.

Part 3:

import operator

l = [True, False, True] 

def myReduce(f, l, i): 
	intermediateResult = i 
	for element in l: 
		intermediateResult = f(intermediateResult, element) 
	return intermediateResult 

r = myReduce(operator.and_, l, True) 
print(r)  #  output will be False in this case

Maybe you were expecting that an implementation of reduce would be more complicated, but it’s actually quite simple. We set up a variable to always contain the intermediate result while working through the elements in the list and initialize it with the initial value provided in the third parameter i. When looping through the elements, we always apply the function given in parameter f to the intermediate result and the element itself and update the intermediate result variable with the result of this operation. At the end, we return the value of this variable as the result of the entire reduce operation.

import pandas as pd 

# create the data frame from a list of tuples 
data = pd.DataFrame( [('Mike',7,10,5.5),
     ('Lisa', 6.5, 9, 8),
     ('George', 4, 3, 7),
     ('Maria', 7, 9.5, 4),
     ('Frank', 5, 5, 5) ] )
     
# set column names
data.columns = ['Name', 'Assignment 1', 'Assignment 2', 'Assignment 3']

# set row names
data.index = range(1,len(data)+1)

# show table 
print(data)
 
# add column with averages
data['Average'] = (data['Assignment 1'] + data['Assignment 2'] + data['Assignment 3']) / 3
 
# part a (all students with a1 score < 7)
print(data[ data['Assignment 1'] < 7])
 
# part b (all students with a1 and a2 score > 6)
print(data[ (data['Assignment 1'] > 6) & (data['Assignment 2'] > 6)])

# part c (at least one assignment < 5) 
print( data[ data[ ['Assignment 1', 'Assignment 2', 'Assignment 3'] ].min(axis = 1) < 5 ] )
 
# part d (name starts with M, only Name and Average columns)
print(data [ data [ 'Name' ].map(lambda x: x.startswith('M')) ] [ ['Name','Average'] ])

# sort by Name 
print(data.sort_values(by = ['Name']))

If any of these steps is unclear to you, please ask for further explanation on the forums.

Sets are another built-in collection in Python in addition to lists, tuples, and dictionaries. The idea is that of a mathematical set, meaning that there is no order between the elements and an element can only be contained in a set once (in contrast to lists). Sets are mutable like lists or dictionaries.

The following code example shows how we can create a set using curly brackets {…} to delimit the elements (similar to a dictionary but without the : part separating keys from values) and that all duplicates are automatically removed. The comparison between elements is done using the == operator. One restriction of sets is that they can only contain immutable values.

s = {3,4,1,3,4,1} # create set 
print(s) 

Output: 
{1, 3, 4} 

Since sets are unordered, it is not possible to access their elements via an index but we can use the “in” operator to test whether or not a set contains an element as well as use a for-loop to iterate through the elements:

x = 3 
if x in s: 
     print("already contained") 

for e in s:
    print(e)

Output: 
already contained 
1 
3 
4 

One of the nice things about sets is that they provide the standard set theoretical operations union, intersection, etc. as shown in the following code example:

group1 = { "Jim", "Maria", "Frank", "Susan"} 
group2 = { "Sam", "Steve", "Jim" }

print( group1 | group2 )  # or group1.union(group2) 

print( group1 & group2 )  # or group1.intersection(group2) 

print( group1 - group2 )  # or group1.difference(group2)  

print( group1 ^ group2 )  # or group1.symmetric_difference(group2)

Output: 
{'Frank', 'Sam', 'Steve', 'Susan', 'Maria', 'Jim'} 
{'Jim'} 
{'Susan', 'Frank', 'Maria'} 
{'Frank', 'Sam', 'Steve', 'Susan', 'Maria'} 

The difference between the last and second-to-last operation here is that group1 - group2 returns the elements of the set in group1 that are not also elements of group2, while the symmetric difference operation group1 ^ group2 returns a set with all elements that are only contained in one of the groups but not in both.

One common operation on collections is sorting the elements in a collection. Python provides a function sorted(…) to sort the elements in a collection that allows for iterating over the elements (e.g. a list, tuple, or dictionary). The result is always a list. Here are two examples:

l1 = [9,3,5,1,-2]
print(sorted(l1))

l2 = ("Maria", "Frank", "Sam", "Mike")
print(sorted(l2))

Output: 
[-2, 1, 3, 6, 9] 
['Frank', 'Maria', 'Mike', 'Sam'] 

If our collection is a list, we can also use the list method sort() instead of the function sorted(…), e.g. l1.sort() instead of sorted(l1) . Both work exactly the same.

sorted(…) and sort() by default sort the elements in ascending order based on the < comparison operator. This means numbers are sorted in increasing order and strings are sorted in lexicographical order. When we define our own classes (Section 4.6) and want to be able to sort objects of a class based on their properties, we have to define the < operator in a suitable way in the class definition.

The keyword argument ‘reverse’ can be used to sort the elements in descending order instead:

print( sorted(l2, reverse = True) )

Output:

['Sam', 'Mike', 'Maria', 'Frank']

In addition, we can use the keyword argument ‘key’ to provide a function that will be applied to the elements and they will then be sorted based on the values returned by this function. For instance, the following example uses a lambda expression for the ‘key’ parameter to sort the names from l2 based on their length (in descending order) rather than based on their lexicographical order:

print( sorted(l2, reverse = True, key = lambda x: len(x)) )

Output:

['Maria', 'Frank', 'Mike', 'Sam']

Sorting can be a somewhat time-consuming operation for larger collections. Therefore, if you mainly need to access the elements in a collection in a specific order based on their properties (e.g. always the element with the lowest value for a certain attribute), it is advantageous to use a specialized data structure that keeps the collection sorted whenever an element is added or removed. This can save a lot of time compared to frequently re-sorting the collection. An example of such a data structure is the so-called priority queue or heap. The heapq module of the Python standard library implements an algorithm for realizing a priority queue in a Python list and we are going to discuss it in the final part of this section.

The idea of a priority queue is that items in the collection are always kept ordered based on their < relation so that when we take the first item from the queue it will always be the one with lowest (or highest) value.

For instance, let’s get back to the example we started this section with of managing a queue of assignments or tasks that need to be performed. Let’s say that instead of performing the assignments in the order in which they arrive (first-in-first-out), the assignments have a priority value between 1 and 9 with 1 meaning highest and 9 meaning lowest priority. That means we need to make sure we keep the assignments in the queue ordered based on their priority so that taking the first assignment from the queue will be that with the highest priority.

The heapq module among other things provides a set of functions for adding elements to a list (function heappush(…)) and for removing the first item with highest priority (function heappop(…)). In the following code, we again use strings for representing our assignments and encode the priority in the strings themselves so that their lexicographical order corresponds to their priority, i.e. “Assignment 1” < “Assignment 2” < … < “Assignment 9”. The reason we defined the highest priority to be given by the number 1 and the lowest priority by the number 9 is that heapq implements a min heap in which heappop(…) always returns the lowest value element according to the < relation in contrast to a max heap in which heappop(…) would always return the highest value element. The code starts with an empty list in variable pQueue and then simulates the arrival of 100 assignments with random priority using heappush(…) to add a new assignment to the queue.

import heapq 
import random 

pQueue = [] 

for count in range(1,101): 
    priority = random.randint(1,9) 
    heapq.heappush(pQueue, 'Assignment ' + str(priority)) 

print(pQueue) 

When you look at the output produced by the print statement in the last line, you may be disappointed because it doesn’t look like the list is really ordered based on the priority numbers of the assignments. However, the list also does not reflect the order in which the assignments have been added to the queue. The list is actually a “flattened” representation of a binary tree, the data structure that heapq is using to make the push and pop operations as efficient as possible, while making sure that heappop(…) always gives you the lowest value element from the queue.

Now add the following code that calls heappop(…) 100 times to remove all assignments from the queue and print out their names including their priority value:

for count in range(1,101): 
    assignment = heapq.heappop(pQueue) 
    print(assignment)

Output: 

Assignment 1 
Assignment 1 
…
Assignment 2 
…
Assignment 9 

As you can see, by using heappop(…) we indeed get the assignments in the right order from the queue. Of course, this is a simplified example in which we first fill the queue completely and then empty it again, but it works in the same way if we add and remove assignments in any arbitrary order. Using heapq for this task is much, much faster than any simple approach such as always searching through the entire list to find the element with the lowest value or, slightly better, always searching for the correct position when inserting a new assignment into the list to keep the list sorted. If you don't believe it, try to implement your own method and do some profiling to see how it compares to the priority queue based approach. 

In the walkthrough of this lesson, we will employ this notion of a priority queue for keeping a number of bus track GPS observations sorted based on their timestamps. For this, we will have to define the < method for our observation points in a suitable way to work with heapq. This will allow us to process the observation points in chronological order.

This section gave you a bit of a taste of the idea of efficient data structures for collections, the algorithms behind them, and the trade-offs involved (a data structure that is very efficient for certain operations will be suboptimal for other operations). Computer science students spend a lot of time studying the implementation and properties of such data structures and the time and space complexities of the operations involved. We were only able to scratch the surface of this topic here, but, as indicated above, there are many books and other resources on this topic, including some specifically written for Python.

In this section, we will provide instructions for installing QGIS via the OSGeo4W distribution manager and for setting up your system to be prepared for the QGIS programming work, we are going to do in this lesson. The OSGeo4W/QGIS installation includes its own Python 3 environment and you will have make sure that you use this Python installation for running the qgis based examples from the later sections. One way to achieve this is by executing the scripts via commands in the OSGeo4W shell, after executing some commands that make sure that all environment variables are set up correctly. This will also be explained below.

OSGeo4W and QGIS installation

To install the OSGeo4W environment with QGIS 3.x, please follow the steps below:

1. Go to the qgis website
2. Click "Download now"
3. Pick the Windows "OSGeo4W Network Installer" close to the top (not the green button for the standalone installer!)
        

   Figure 4.1 Downloading the OSGeo4W installer
4. Click on the link that says "Download OSGeo4W Installer and start it" to downloaded a file called osgeo4w-setup.exe . Then run that file to start the installation.
5. Select the "Express Desktop Install" option (Unless you already have OSGeo4W/QGIS installed, then use the Advanced option to ensure you get the latest versions of all packages/tools or pick a different installation folder. If you get errors attempting to run QGIS, you may have to delete your c:\OSGEO4W folder and re-run the installation.)
   

   Figure 4.2 Selecting the Express Desktop Install option
6. If asked, select a site from which to install (does not matter which) 
7. When asked which packages to install, select the options shown in the screenshot below. Apart from that you can accept the default settings on the next pages.
   

   Figure 4.2b Selecting the packages to install
8. Accept the different licenses and run the installation (ignore warnings about corrupted packages if you get them and just click 'Next' if you get warnings about missing dependencies)

Where to find what?

After the installation has finished, you should have a folder called OSGeo4W in the root folder of your C: drive (unless you picked a different folder for the installation). Here we list the main programs from this installation folder that you will need in this lesson:

• C:\OSGeo4W\OSGeo4W.bat - This opens the OSGeo4W shell that can be used for executing python scripts from the command line.
• C:\OSGeo4W\bin\qgis-ltr-bin.exe - This is the main QGIS executable that you need to run for starting QGIS 3.
• C:\OSGeo4W\apps\Qt5\bin\designer.exe - This is the QT Designer executable that you can use for creating Qt5 GUIs in this lesson. If you simply double-click the .exe file in the Windows File Explorer, you will most likely get an DLL related error message because some environment variables won't be set correctly. But you should be able to run the program by opening the OSGeo4W shell and typing the command  designer there.
• python-qgis-ltr - As explained further below, you can use this command in the OSGeo4W shell for setting the path and environment variables for running qgis and PyQT5 based Python 3 code as well as executing scripts directly.

Running OSGeo4W shell and commands for qgis and PyQt5 development

When you run OSGeo4W.bat, the OSGeo4W shell will show up looking similar to the normal Windows command line but providing some additional commands that can be listed by typing in "o-help".

Figure 4.3 OSGeo4W shell after the o-help command

When using the OSGeo4W shell in this lesson, it is best to always execute the command

python-qgis-ltr

first to make sure all environment variables are set up correctly for running qgis and PyQt5 based Python code. The command will start a Python interpreter (recongnizable by the >>> prompt) that you can immediately leave again by typing the command  quit() . You can also directly run Python scripts with python-qgis-ltr by writing

python-qgis-ltr xyz.py

rather than just

python xyz.py

You can also use the command  pyrcc5  in the OSGeo4W shell for compiling QT5 resource files that we will need later on in this lesson.

Installing geopy package and pandas

Most of the Python packages we will need in this lesson (like PyQt5) are already installed in the Python environment that comes with OSGeo4W/QGIS, but a few additional pieces are necessary. There is one package that we will use for performing distance calculations between WGS84 points in the two walkthroughs of the lesson. The package is called geopy and it needs to be installed first. To do this, please open the OSGeo4W shell and change to Python 3 by running the python-qgis-ltr command followed by quit() as described above, and then run the following pip installation command:

python -m pip install geopy

The package is small, so the installation should only take a couple of seconds. The output you are getting may look slightly different than what is shown in the image below but should indicate that geopy has been installed successfully.

Figure 4.4 Installing geopy in the OSGeo4W shell

In the practice exercise for this lesson, we will also use pandas. In earlier versions of QGIS/OSGeo4W, pandas wasn't installed by default. To make sure, simply run the following command for installing pandas; most likely it's going to tell you that pandas is already installed:

python -m pip install pandas

Installing the required QGIS plugins

We will need a few QGIS plugins in this lesson, so let's install those as well. Some of these are for the optional part at the end but they are small and installation should be quick, so let's install all of them now. Please follow the instructions below for this:

1. Start QGIS
2. Go to Plugins -> Manage and Install Plugins in the main menu bar
3. Under "Settings", make sure the box next to "Show also experimental plugins" is checked.
4. Under "Not installed" look for the following three plugins and install them:
   • QuickMapServices (allows for quickly adding basemaps like OSM to a project)
   • Plugin Builder 3 (for creating templates for new plugins)
   • Plugin Reloader (for reloading a plugin after modifying the code)

Figure 4.5 QGIS Plugin installation

If you now click on "Installed", all three plugins should appear in the list of installed plugins with a checkmark on the left, which indicates that the plugin is activated.

Important note: This lesson has a lot of content and this is one of its less important sections. We included it so that, if you have not worked with QGIS before, you get an idea of where to find what and how things work in QGIS in general. However, since we will mainly be using the QGIS programming API rather than doing things in QGIS itself, we recommend that you go through this section quickly and then maybe come back at the end of the lesson if you have an interest in learning more about QGIS and its interface. 

When you open QGIS 3 for the first time, it will look similar to the image below. The main elements are the main menu bar at the top, a number of horizontal toolbars with buttons for different operations below the menu bar, a smaller vertical toolbar on the left side with buttons for adding or creating layers, and then three main windows: a panel with a file browser, a panel that lists the layers in your project (currently empty), and then the main window for displaying the current project. At the very bottom, you can find a status bar displaying information related to the project window such as the scale and coordinate reference system used. Overall, this all looks somewhat similar to ArcGIS Desktop or Pro. All toolbars and panels can be freely moved around, undocked and docked back again, and there are many additional panels and toolbars that can be enabled/disabled either from the main menu under View -> Panels/Toolbars or by doing a right-click on one of the panel title bars or toolbar areas at the top and left.

Figure 4.6 QGIS after starting it for the first time

There are several ways to add a data set to a project:

1. By navigating to a file in the file browser panel and then double-clicking it.
2. By dragging a file from the Windows File Explorer onto the project window or layers panel.
3. By clicking the “Open Data Source Manager” button, which will open up a dialog with a list of different types of sources on the left including local files, datasets from different databases, and also data sets provided as web services (WMS, WFS, ArcGIS Map Server or Feature Server, etc.).

Figure 4.7 Dialog for adding data layers to the project

Feel free to try out adding different data sets to the project. Similar to ArcGIS, the coordinate reference system used for the project and project window will be that of the first source added, but of course this can be changed, e.g. by going Project -> Project Properties…. in the menu bar or by left-clicking the CRS field in the status bar. Dragging the layers and the buttons at the top of the Layers panel can be used to arrange the layers in a certain order and group or filter them. We here add the world borders layer from Lesson 3 to the project. The layer now shows up in the project window and the Layers panel. Right-clicking the layer in the Layers panel will provide a number of options for that layer. Double-clicking the layer will directly open the “Layer Properties” dialog with a lot of options to change rendering or other properties of the layer.

Figure 4.8 Right-clicking a layer in the Layers panel opens the Layer context menu

Figure 4.9 Layer Properties dialog box

The properties you will most commonly work with are the Symbology and Labels properties. When coming from ArcGIS, working with these dialogs requires a bit of getting used to. Give it a try by attempting to show the world borders layer with a Graduated scheme based on the “AREA” attribute of the layer using a Natural Breaks classification with 8 classes and with labels based on the “NAME” attribute. The result should look somewhat similar to the image below. If you have any problems achieving this, please post on the Lesson 4 discussion forum.

Figure 4.10 World borders layer after changing the symbology

If you want to select features from a layer based on attribute, the Query Builder dialog can be opened by doing a right-click -> Filter … on the layer in the Layers panel. The dialog itself works roughly similar to the corresponding component of ArcGIS. You can check out the attribute table of the layer by doing a right-click -> Open Attribute Table. Working with the attribute table again is roughly similar to ArcGIS. If you want to export a layer as a new data set, you do a right-click -> Export -> Save Features as… . This, for instance, allows for saving only the currently selected features and/or saving the layer in a different format or using a different CRS.

Looking at the main menu bar, we find the main tools for working with Vector and Raster data under the respective submenus. They include typical geoprocessing, data manipulation, and analysis tools. Additional tools can be accessed by opening the Processing Toolbox panel under Processing -> Toolbox. Moreover, QGIS has a plugin interface that allows for writing extensions to QGIS. Plugins can be managed and new plugins can be installed under Plugins -> Manage and Install Plugins, and they can add new entries to menu bar and tool bars. QGIS plugins are written in Python, and you will learn how to do so later on in this lesson. QGIS also has a Python Console (Plugins -> Python Console) that allows for entering and executing Python code that uses the QGIS Python API.

A QGIS project is saved as a .qgz file using Project -> Save or Project -> Save As…. From this menu, you can also open a new project, export the project map in different formats, etc. 

One thing that works a bit differently than in ArcGIS is the layout composer component for creating map views of your project including additional elements such as a legend, scale bar, etc. By going Project -> New Print Layout, you can create a new map layout document. This opens up a new window with its own interface that allows you to arrange maps and other elements like images and text in the same way as in a vector graphics or publishing tool. The created layout can just be a single page or span multiple pages and contain different maps. Elements are added to the page with the buttons from the toolbar on the left. A list of all elements is shown in the panel on the top right. The properties of the currently selected element can be accessed and changed with the panel on the bottom right. The simple layout in the image below was created by adding our current map with the add map  button, adding a text element with the add text button , and then adding a legend for the current map with the add legend  button.

Figure 4.11 Map view created in the QGIS layout composer

Layouts can be exported as images or PDF files and previously created layouts can be accessed via the Layout Manager under Project -> Layout Manager… or directly be accessed from Project -> Layouts -> … .

This short overview should be enough to get you started but, of course, only covers the basics. This lesson will focus on the QGIS Python API and using it to write programs or plugins for QGIS, rather than on working with the QGIS interface directly. Nevertheless, if you want to learn more about QGIS at some point, the following tutorials covering certain tasks in more detail can be used as a starting point. 

• Making a map
• Working with attributes
• Basic vector styling
• Raster styling and analysis

More tutorials are available at this QGIS Tutorials and Tips page.

Let’s start this introduction by writing some code directly in the QGIS Python console and talking about how you can access the layers currently open in QGIS and add new layers to the currently open project. If you don’t have QGIS running at the moment, please start it up and open the Python console from the Plugins menu in the main menu bar.

When you open the Python console in QGIS, the Python qgis package and its submodules are automatically imported as well as other relevant modules including the main PyQt5 modules. In addition, a variable called iface is set up to provide an object of the class QgisInterface¹ to interact with the running QGIS environment. The code below shows how you can use that object to retrieve a list of the layers in the currently open map project and the currently active layer. Before you type in and run the code in the console, please add a few layers to an empty project including the TM_WORLD_BORDERS-0.3.shp shapefile that we already used in Section 3.9.1 on GDAL/ORG. We will recreate some of the steps from that section with QGIS here so that you also get a bit of a comparison between the two APIs. The currently active layer is the one selected in the Layers window; please select the world borders layer by clicking on it before you execute the code.

layers = iface.mapCanvas().layers() 
for layer in layers:
	print(layer)
	print(layer.name())
	print(layer.id())
	print('------')
# If you copy/paste the code - run the part above
# before you run the part below 
# otherwise you'll get a syntax error.
activeLayer = iface.activeLayer()
print('active layer: ' + activeLayer.name())

Output (numbers will vary): 
...
<qgis._core.QgsVectorLayer object at 0x000000666CF22D38>
TM_WORLD_BORDERS-0.3 
TM_WORLD_BORDERS_0_3_2e5a7cd5_591a_4d45_a4aa_cbba2e639e75 
------ 
...
active layer:  TM_WORLD_BORDERS-0.3 

The layers() method of the QgsMapCanvas object we get from calling iface.mapCanvas() returns the currently open layers as a list of objects of the different subclasses of QgsMapLayer. Invoking the name() method of these layer objects gives us the name under which the layer is listed in the Layers window. layer.id() gives us the ID that QGIS has assigned to the layer which in contrast to the name is unique. The iface.activeLayer() method gives us the currently selected layer.

The type() function of a layer can be used to test the type of the layer:

if activeLayer.type() == QgsMapLayer.VectorLayer: 
	print('This is a vector layer!')

Depending on the type of the layer, there are other methods that we can call to get more information about the layer. For instance, for a vector layer we can use wkbType() to get the geometry type of the layer:

if activeLayer.type() == QgsMapLayer.VectorLayer: 
	if activeLayer.wkbType() == QgsWkbTypes.MultiPolygon: 
		print('This layer contains multi-polygons!')

The output you get from the previous command should confirm that the active world borders layer contains multi-polygons, meaning features that can have multiple polygonal parts.

QGIS defines a function dir(…) that can be used to list the methods that can be invoked for a given object. Try out the following two applications of this function:

dir(iface) 
dir(activeLayer)

To add or remove a layer, we need to work with the QgsProject object for the project currently open in QGIS. We retrieve it like this:

currentProject = QgsProject.instance() 
print(currentProject.fileName())

The output from the print statement in the second row will probably be the empty string unless you have saved the project. Feel free to do so and rerun the line and you should get the actual file name.

Here is how we can remove the active layer (or any other layer object) from the layer registry of the project (you may have to resize/refresh the map canvas afterwards for the layer to disappear there):

currentProject.removeMapLayer(activeLayer.id())

The following command shows how we can add the world borders shapefile again (or any other feature class we have on disk). Make sure you adapt the path based on where you have the shapefile stored. We first have to create the vector layer object providing the file name and optionally the name to be used for the layer. Then we add that layer object to the project via the addMapLayer(…) method:

layer = QgsVectorLayer(r'C:\489\TM_WORLD_BORDERS-0.3.shp', 'World borders') 
currentProject.addMapLayer(layer)

Lastly, here is an example that shows you how you can change the symbology of a layer from your code:

renderer = QgsGraduatedSymbolRenderer() 
renderer.setClassAttribute('POP2005') 
layer.setRenderer(renderer) 
layer.renderer().updateClasses(layer, QgsGraduatedSymbolRenderer.Jenks, 5) 
layer.renderer().updateColorRamp(QgsGradientColorRamp(Qt.white, Qt.red)) 
iface.layerTreeView().refreshLayerSymbology(layer.id())
iface.mapCanvas().refreshAllLayers()

Here we create an object of the QgsGraduatedSymbolRenderer class that we want to use to draw the country polygons from our layer using a graduated color approach based on the population attribute ‘POP2005’. The name of the field to use is set via the renderer’s setClassAttribute() method in line 2. Then we make the renderer object the renderer for our world borders layer in line 3. In the next two lines, we tell the renderer (now accessed via the layer method renderer()) to use a Jenks Natural Breaks classification with 5 classes and a gradient color ramp that interpolates between the colors white and red. Please note that the colors used as parameters here are predefined instances of the Qt5 class QColor. Changing the symbology does not automatically refresh the map canvas or layer list. Therefore, in the last two lines, we explicitly tell the running QGIS environment to refresh the symbology of the world borders layer in the Layers tree view (line 6) and to refresh the map canvas (line 7). The result should look similar to the figure below (with all other layers removed).

     

Figure 4.12 World borders layer after changing the symbology from Python

You will get to see another example of interacting with the layers open in QGIS and setting the symbology (for point and line layers in this case) in Section 4.12 where we take the code from this lesson's walkthrough and turn it into a QGIS plugin.

[1] The qgis Python module is a wrapper around the underlying C++ library. The documentation pages linked in this section are those of the C++ version but the names of classes and available functions and methods are the same.

Let’s keep working with the world borders layer open in QGIS for a bit, looking at how we can access the individual features in a layer and select features by attribute. The following piece of code shows you how we can loop through all the features with the help of the layer’s getFeatures() method:

for feature in layer.getFeatures(): 
    print(feature) 
    print(feature.id()) 
    print(feature['NAME']) 
    print('-----') 

Output: 
<qgis._core.QgsFeature object at 0x...> 
0 
Antigua and Barbuda 
----- 
<qgis._core.QgsFeature object at 0x...> 
1 
Algeria 
----- 
<qgis._core.QgsFeature object at 0x...> 
2 
Azerbaijan 
----- 
<qgis._core.QgsFeature object at 0x...> 
3 
Albania 
----- 
... 

Features are represented as objects of the class QgsFeature in QGIS. So, for each iteration of the for-loop in the previous code example, variable feature will contain a QgsFeature object. Features are numbered with a unique ID that you can obtain by calling the method id() as we are doing in this example. Attributes like the NAME attribute of the world borders polygons can be accessed using the attribute name as the key as also demonstrated above.

Like in most GIS software, a layer can have an active selection. When the layer is open in QGIS, the selected features are highlighted. The layer method selectAll() allows for selecting all features in a layer and removeSelection() can be used to clear the selection. Give this a try by running the following two commands in the QGIS Python console and watch how all countries become selected and then deselected again.

layer.selectAll() 
layer.removeSelection()

The method selectByExpression() allows for selecting features based on their properties with a SQL query string that has the same format as in ArcGIS. Use the following command to select all features from the layer that have a value larger than 300,000 in the AREA column of the attribute table. The result should look as in the figure below.

layer.selectByExpression('"AREA" > 300000')

     

Figure 4.13 World borders layer with countries selected by selectByExpression(...) selected

While there can only be one active selection for a layer, you can create as many subgroups of features from a layer as you want by calling getFeatures(…) with a parameter that is an object of the class QgsFeatureRequest and that has been given a filter expression via its setFilterExpression(…) method. The filter expression can be again an SQL query string. The following code creates a subgroup that will only contain the polygon for Canada. When you run it, this will not change the active selection that you see for that layer in QGIS but variable selectionName now provides access to the subgroup with just that one polygon. We get that first (and only) polygon by calling the __next__() method of selectionName and then print out some information about this particular polygon feature.

selectionName = layer.getFeatures(QgsFeatureRequest().setFilterExpression('"NAME" = \'Canada\'')) 
feature = selectionName.__next__() 
print(feature['NAME'] + "-" + str(feature.id()))
print(feature.geometry()) 
print(feature.geometry().asWkt())

Output: 
Canada – 23 
<qgis._core.QgsGeometry object at 0x...> 
MultiPolygon (((-65.61361699999997654 43.42027300000000878,...)))

The first print statement in this example works in the same way as you have seen before to get the name attribute and id of the feature. The method geometry() gives us the geometric object for this feature as an instance of the QgsGeometry class and calling the method asWkt() gives us a WKT string representation of the multi-polygon geometry. You can also use a for-loop to iterate through the features in a subgroup created in this way. The method rewind() can be used to reset the iterator to the beginning so that when you call __next__() again, it will again give you the first feature from the subgroup.

When you have the geometry object and know what type of geometry it is, you can use the methods asPoint(), asPolygon(), asPolyline(), asMultiPolygon(), etc. to get the geometry as a Python data structure, e.g. in the case of multi-polygons as a list of lists of lists with each inner list containing tuples of the point coordinates for one polygonal component.

print(feature.geometry().asMultiPolygon())

[[[(-65.6136, 43.4203), (-65.6197,43.4181), … ]]]

Here is another example to demonstrate that we can work with several different subgroups of features at the same time. This time we request all features from the layer that have a POP2005 value larger than 50,000,000.

selectionPopulation = layer.getFeatures(QgsFeatureRequest().setFilterExpression('"POP2005" > 50000000'))

If we ever want to use a subgroup like this to create the active selection for the layer from it, we can use the layer method selectByIds(…) for this. The method requires a list of feature IDs and will then change the active selection to these features. In the following example, we use a simple list comprehension to create the ID list from the subgroup in our variable selectionPopulation:

layer.selectByIds([f.id() for f in selectionPopulation]) 

When running this command you should notice that the selection of the features in QGIS changes to look like in the figure below.

     

Figure 4.14 Selection created by calling selectByIds(...)

Let’s save the currently selected features as a new file. We use the GeoPackage format (GPKG) for this, which is more modern than the shapefile format, but you can easily change the command below to produce a shapefile instead; simply change the file extension to “.shp” and replace “GPKG” with “ESRI Shapefile”. The function we will use for writing the layer to disk is called writeAsVectorFormat(…) and it is defined in the class QgsVectorFileWriter. Please note that this function has been declared "deprecated", meaning it may be removed in future versions and it is recommended that you do not use it anymore. In versions up to QGIS 3.16 (the current LTR version that most likely you are using right now), you are supposed to use writeAsVectorFormatV2(...) instead; however, there have been issues reported with that function and it is already replaced by writeAsVectorFormatV3(...) in versions >3.16 of QGIS. Therefore, we have decided to stick with writeAsVectorFormat(…) while things are still in flux. The parameters we give to writeAsVectorFormat(…) are the layer we want to save, the name of the output file, the character encoding to use, the spatial reference to use (we simply use the one that our layer is in), the format (“GPKG”), and True for signaling that only the selected features should be saved in the new data set. Adapt the path for the output file as you see fit and then run the command:

QgsVectorFileWriter.writeAsVectorFormat(layer, r'C:\489\highPopulationCountries.gpkg', 'utf-8', layer.crs(),'GPKG', True)

If you add the new file produced by this command to your QGIS project, it should only contain the polygons for the countries we selected based on their population values.

For changing the attribute values of a feature, we need to work with the “data provider” object of the layer. We can access it via the layer’s dataProvider() method:

dataProvider = layer.dataProvider()

Let’s say we want to change the POP2005 value for Canada to 1 (don’t ask what happened!). For this, we also need the index of the POP2005 column which we can get by calling the data provider’s fieldNameIndex() method:

populationColumnIndex = dataProvider.fieldNameIndex('POP2005')

To change the attribute value we call the method changeAttributeValues(…) of the data provider object providing a dictionary as parameter that maps feature IDs to dictionaries which in turn map column indices to new values. The inner dictionary that maps column indices to values is defined in a separate variable newValueDictionary.

newValueDictionary = { populationColumnIndex : 1 } 
dataProvider.changeAttributeValues( { feature.id(): newValueDictionary } )

In this simple example, the outer dictionary contains only a single key-value pair with the ID of the feature for Canada as key and another dictionary as value. The inner dictionary also only contains a single key-value pair consisting of the index of the population column and its new value 1. Both dictionaries can have multiple entries to simultaneously change multiple values of multiple features. After running this command, check out the attributes of Canada, either via the QGIS Identify tool or in the attribute table of the layer. You will see that the population value in the layer now has been changed to 1 (the same holds for the underlying shapefile). Let’s set the value back to what it was with the following command:

dataProvider.changeAttributeValues( { feature.id(): { populationColumnIndex : 32270507 } } )

For the final part of this section, let’s switch from the Python console in QGIS to writing a standalone script that uses qgis. You can use your editor of choice to write the script and then execute the .py file from the OSGeo4W shell (see again Section 4.4.1) with all environment variables set correctly for a qgis and QT5 based program.

We are going to repeat the task from Section 3.9.1 of creating buffers around the centroids of the countries within a rectangular (in terms of WGS 84 coordinates) area around southern Africa. We will produce two new vector GeoPackage files: a point based one with the centroids and a polygon based one for the buffers. Both data sets will only contain the country name as their only attribute.

We start by importing the modules we will need and creating a QApplication() (handled by qgis.core.QgsApplication) for our program that qgis can run in (even though the program does not involve any GUI).

Important note: When you later write you own qgis programs (e.g. in the L4 homework assignment), make sure that you always "import qgis" first before using any other qgis related import statements such as "import qgis.core". We are not sure why this is needed, but the other imports will most likely fail tend to fail without "import qgis" coming first.

import os, sys 
import qgis 
import qgis.core

To use qgis in our software, we have to initialize it and we need to tell it where the actual QGIS installation is located. To do this, we use the function getenv(…) of the os module to get the value of the environmental variable “QGIS_PREFIX_PATH” which will be correctly defined when we run the program from the OSGeo4W shell. Then we create an instance of the QgsApplication class and call its initQgis() method.

qgis_prefix = os.getenv("QGIS_PREFIX_PATH")      
qgis.core.QgsApplication.setPrefixPath(qgis_prefix, True) 
qgs = qgis.core.QgsApplication([], False)
qgs.initQgis()

Now we can implement the main functionality of our program. First, we load the world borders shapefile into a layer (you may have to adapt the path!).

layer = qgis.core.QgsVectorLayer(r'C:\489\TM_WORLD_BORDERS-0.3.shp')

Then we create the two new layers for the centroids and buffers. These layers will be created as new in-memory layers and later written to GeoPackage files. We provide three parameters to QgsVectorLayer(…): (1) a string that specifies the geometry type, coordinate system, and fields for the new layer; (2) a name for the layer; and (3) the string “memory” which tells the function that it should create a new layer in memory from scratch (rather than reading a data set from somewhere else as we did earlier).

centroidLayer = qgis.core.QgsVectorLayer("Point?crs=" + layer.crs().authid() + "&field=NAME:string(255)", "temporary_points", "memory") 
bufferLayer = qgis.core.QgsVectorLayer("Polygon?crs=" + layer.crs().authid() + "&field=NAME:string(255)", "temporary_buffers", "memory")

The strings produced for the first parameters will look like this: “Point?crs=EPSG:4326&field=NAME:string(255)” and “Polygon?crs=EPSG:4326&field=NAME:string(255)”. Note how we are getting the EPSG string from the world border layer so that the new layers use the same coordinate system, and how an attribute field is described using the syntax “field=<name of the field>:<type of the field>". When you want your layer to have more fields, these have to be separated by additional & symbols like in a URL.

Next, we set up variables for the data providers of both layers that we will need to create new features for them. The new features will be collected in two lists, centroidFeatures and bufferFeatures.

centroidProvider = centroidLayer.dataProvider() 
bufferProvider = bufferLayer.dataProvider() 

centroidFeatures = []
bufferFeatures = []

Then, we create the polygon geometry for our selection area from a WKT string as in Section 3.9.1:

areaPolygon = qgis.core.QgsGeometry.fromWkt('POLYGON ( (6.3 -14, 52 -14, 52 -40, 6.3 -40, 6.3 -14) )')

In the main loop of our program, we go through all the features in the world borders layer, use the geometry method intersects(…) to test whether the country polygon intersects with the area polygon, and, if yes, create the centroid and buffer features for the two layers from the input feature.

for feature in layer.getFeatures(): 
    if feature.geometry().intersects(areaPolygon): 
        centroid = qgis.core.QgsFeature() 
        centroid.setAttributes([feature['NAME']]) 
        centroid.setGeometry(feature.geometry().centroid()) 
        centroidFeatures.append(centroid) 

        buffer = qgis.core.QgsFeature() 
        buffer.setAttributes([feature['NAME']]) 
        buffer.setGeometry(feature.geometry().centroid().buffer(2.0,100)) 
        bufferFeatures.append(buffer)

Note how in both cases (centroids and buffers), we first create a new QgsFeature object, then use setAttributes(…) to set the NAME attribute to the name of the country, and then use setGeometry(…) to set the geometry of the new feature either to the centroid derived by calling the centroid() method or to the buffered centroid created by calling the buffer(…) method of the centroid point. As a last step, the new features are added to the respective lists. Finally, all features in the two lists are added to the layers after the for-loop has been completed. This happens with the following two commands:

centroidProvider.addFeatures(centroidFeatures) 
bufferProvider.addFeatures(bufferFeatures)

Lastly, we write the content of the two in-memory layers to GeoPackage files on disk. This works in the same way as in previous examples. Again, you might want to adapt the output paths.

qgis.core.QgsVectorFileWriter.writeAsVectorFormat(centroidLayer, r'C:\489\centroids.gpkg', "utf-8", layer.crs(), "GPKG") 
qgis.core.QgsVectorFileWriter.writeAsVectorFormat(bufferLayer, r'C:\489\buffers.gpkg', "utf-8", layer.crs(), "GPKG")

Since we are now done with using QGIS functionalities (and actually the entire program), we clean up by calling the exitQgis() method of the QgsApplication, freeing up resources that we don’t need anymore.

qgs.exitQgis()

If you run the program from the OSGeo4W shell and then open the two produced output files in QGIS, the result should look as shown in the image below.

     

Figure 4.15 Centroid and buffer layers created with the previous Python code

QGIS has a toolbox system and visual workflow building component somewhat similar to ArcGIS and its Model Builder. It is called the QGIS processing framework and comes in the form of a plugin called Processing that is installed by default. You can access it via the Processing menu in the main menu bar. All algorithms from the processing framework are available in Python via a QGIS module called processing. They can be combined to solve larger analysis tasks in Python and can also be used in combination with the other qgis methods discussed in the previous sections.

We can get a list of all processing algorithms currently registered with QGIS with the command QgsApplication.processingRegistry().algorithms(). Each processing object in the returned list has an identifying name that you can get via its id() method. The following command, which you can try out in the QGIS Python console, uses this approach to print the names of all algorithms that contain the word “clip”:

[x.id() for x in QgsApplication.processingRegistry().algorithms() if "clip" in x.id()] 

Output: 

['gdal:cliprasterbyextent', 'gdal:cliprasterbymasklayer','gdal:clipvectorbyextent', 'gdal:clipvectorbypolygon', 'native:clip', 'saga:clippointswithpolygons', 'saga:cliprasterwithpolygon', 'saga:polygonclipping'] 

As you can see, there are processing versions of algorithms coming from different sources, e.g. natively built into QGIS vs. algorithms based on GDAL. The function algorithmHelp(…) allows you to get some documentation on an algorithm and its parameters. Try it out with the following command:

processing.algorithmHelp("native:clip")

To run a processing algorithm, you have to use the run(…) function and provide two parameters: the id of the algorithm and a dictionary that contains the parameters for the algorithm as key-value pairs. run(…) returns a dictionary with all output parameters of the algorithm. The following example illustrates how processing algorithms can be used to solve the task of clipping a points of interest shapefile to the area of El Salvador, reusing the two data sets from homework assignment 2 (Section 2.10). This example is intended to be run as a standalone program again and most of the code is required to set up the QGIS environment needed, including initializing the Processing environment.

The start of the script looks like in the example from the previous section:

import os,sys 
import qgis 
import qgis.core

qgis_prefix = os.getenv("QGIS_PREFIX_PATH")      
qgis.core.QgsApplication.setPrefixPath(qgis_prefix, True) 
qgs = qgis.core.QgsApplication([], False)
qgs.initQgis()

After, creating the QGIS environment, we can now initialize the processing framework. To be able to import the processing module we have to make sure that the plugins folder is part of the system path; we do this directly from our code. After importing processing, we have to initialize the Processing environment and we also add the native QGIS algorithms to the processing algorithm registry.

# Be sure to change the path to point to where your plugins folder is located
# it may not be the same as this one.
sys.path.append(r"C:\OSGeo4W\apps\qgis-ltr\python\plugins")
import processing 
from processing.core.Processing import Processing 
Processing.initialize() 
qgis.core.QgsApplication.processingRegistry().addProvider(qgis.analysis.QgsNativeAlgorithms())

Next, we create input variables for all files involved, including the output files we will produce, one with the selected country and one with only the POIs in that country. We also set up input variables for the name of the country and the field that contains the country names.

poiFile = r'C:\489\L2\assignment\OSMpoints.shp' 
countryFile = r'C:\489\L2\assignment\countries.shp' 
pointOutputFile = r'C:\489\L2\assignment\pointsInCountry.shp' 
countryOutputFile = r'C:\489\L2\assignment\singleCountry.shp' 

nameField = "NAME" 
countryName = "El Salvador"

Now comes the part in which we actually run algorithms from the processing framework. First, we use the qgis:extractbyattribute algorithm to create a new shapefile with only those features from the country data set that satisfy a particular attribute query. In the dictionary with the input parameters for the algorithm, we specify the name of the input file (“INPUT”), the name of the query field (“FIELD”), the comparison operator for the query (0 here stands for “equal”), and the value to which we are comparing (“VALUE”). Since the output will be written to a new shapefile, we don’t really need the output dictionary that we get back from calling run(…) but the print statement shows how this dictionary in this case contains the name of the output file under the key “OUTPUT”.

output = processing.run("qgis:extractbyattribute", { "INPUT": countryFile, "FIELD": nameField, "OPERATOR": 0, "VALUE": countryName, "OUTPUT": countryOutputFile  }) 
print(output['OUTPUT'])

To perform the clip operation with the new shapefile from the previous step, we use the “native:clip” algorithm. The input paramters are the input file (“INPUT”), the clip file (“OVERLAY”), and the output file (“OUTPUT”). Again, we are just printing out the content stored under the “OUTPUT” key in the returned dictionary. Finally, we exit the QGIS environment.

output = processing.run("native:clip", { "INPUT": poiFile, "OVERLAY": countryOutputFile, "OUTPUT": pointOutputFile }) 
print(output['OUTPUT']) 

qgs.exitQgis()

Below is how the resulting two layers should look when shown in QGIS in combination with the original country layer.

     

Figure 4.16 Output files produced by the previous Python code

In this section, we showed you how to perform common GIS operations with the QGIS Python API. Once again we have to say that we are only scratching the surface here; the API is much more complex and powerful, and there is hardly anything you cannot do with it. What we have shown you will be sufficient to understand the code from the two walkthroughs of this lesson, but, if you want more, below are some links to further examples. Keep in mind though that since QGIS 3 is not that old yet, some of the examples on the web have been written for QGIS 2.x. While many things still work in the same way in QGIS 3, you may run into situations in which an example won’t work and needs to be adapted to be compatible with QGIS 3.

• PyQGIS Developer Cookbook
• QGIS Tutorial and Tips (Section on Python Scripting (PyQGIS))
• "How To" in QGIS

Let’s recapitulate a bit: the underlying perspective of object-oriented programming is that the domain modeled in a program consists of objects belonging to different classes. If your software models some part of the real world, you may have classes for things like buildings, vehicles, trees, etc. and then the objects (also called instances) created from these classes during run-time represent concrete individual buildings, vehicles, or trees with their specific properties. The classes in your software can also describe non real-world and often very abstract things like a feature layer or a random number generator.

Class definitions specify general properties that all objects of that class have in common, together with the things that one can do with these objects. Therefore, they can be considered blueprints for the objects. Each object at any moment during run-time is in a particular state that consists of the concrete values it has for the properties defined in its class. So, for instance, the definition of a very basic class Car may specify that all cars have the properties owner, color, currentSpeed, and lightsOn. During run-time we might then create an object for “Tom’s car” in variable carOfTom with the following values making up its state:

carOfTom.owner = "Tom" 
carOfTom.color = "blue" 
carOfTom.currentSpeed = 48   (mph) 
carOfTom.lightsOn = False

While all objects of the same class have the same properties (also called attributes or fields), their values for these properties may vary and, hence, they can be in different states. The actions that one can perform with a car or things that can happen to a car are described in the form of methods in the class definition. For instance, the class Car may specify that the current speed of cars can be changed to a new value and that lights can be turned on and off. The respective methods may be called changeCurrentSpeed(…), turnLightsOn(), and turnLightsOff(). Methods are like functions but they are explicitly invoked on an object of the class they are defined in. In Python this is done by using the name of the variable that contains the object, followed by a dot, followed by the method name:

carOfTom.changeCurrentSpeed(34) # change state of Tom’s car to current speed being 34mph 

carOfTom.turnLightsOn()         # change state of Tom’s car to lights being turned on

The purpose of methods can be to update the state of the object by changing one or several of its properties as in the previous two examples. It can also be to get information about the state of the car, e.g. are the lights turned on? But it can also be something more complicated, e.g. performing a certain driving maneuver or fuel calculation.

In object-oriented programming, a program is perceived as a collection of objects that interact by calling each other’s methods. Object-oriented programming adheres to three main design principles:

• Encapsulation: Definitions related to the properties and methods of any class appear in a specification that is encapsulated independently from the rest of the software code and properties are only accessible via a well-defined interface, e.g. via the defined methods.
• Inheritance: Classes can be organized hierarchically with new classes being derived from previously defined classes inheriting all the characteristics of the parent class but potentially adding specialized properties or specialized behavior. For instance, our class Car could be derived from a more general class Vehicle adding properties and methods that are specific for cars.
• Polymorphism: Inherited classes can change the behavior of methods by overwriting them and the code executed when such a method is invoked for an object then depends on the class of that object.

We will talk more about inheritance and polymorphism in section 4.8. All three principles aim at improving reusability and maintainability of software code. These days, most software is created by mainly combining parts that already exist because that saves time and costs and increases reliability when the re-used components have already been thoroughly tested. The idea of classes as encapsulated units within a program increases reusability because these units are then not dependent on other code and can be moved over to a different project much more easily.

For now, let’s look at how our simple class Car can be defined in Python.

 class Car(): 

     def __init__(self): 
          self.owner = 'UNKNOWN' 
          self.color = 'UNKNOWN' 
          self.currentSpeed = 0 
          self.lightsOn = False 

     def changeCurrentSpeed(self,newSpeed): 
          self.currentSpeed = newSpeed 

     def turnLightsOn(self): 
          self.lightsOn = True 

     def turnLightsOff(self): 
          self.lightsOn = False 

     def printInfo(self): 
          print('Car with owner = {0}, color = {1}, currentSpeed = {2}, lightsOn = {3}'.format(self.owner, self.color, self.currentSpeed, self.lightsOn)) 

Here is an explanation of the different parts of this class definition: each class definition in Python starts with the keyword ‘class’ followed by the name of the class (‘Car’) followed by parentheses that may contain names of classes that this class inherits from, but that’s something we will only see later on. The rest of the class definition is indented to the right relative to this line.

The rest of the class definition consists of definitions of the methods of the class which all look like function definitions but have the keyword ‘self’ as the first parameter, which is an indication that this is a method. The method __init__(…) is a special method called the constructor of the class. It will be called when we create a new object of that class like this:

carOfTom = Car()    # uses the __init__() method of Car to create a new Car object

In the body of the constructor, we create the properties of the class Car. Each line starting with “self.<name of property> = ...“ creates a so-called instance variable for this car object and assigns it an initial value, e.g. zero for the speed. The instance variables describing the state of an object are another type of variable in addition to global and local variables that you already know. They are part of the object and exist as long as that object exists. They can be accessed from within the class definition as “self.<name of the instance variable>” which happens later in the definitions of the other methods, namely in lines 10, 13, 16 and 19. If you want to access an instance variable from outside the class definition, you have to use <name of variable containing the object>.<name of the instance variable>, so, for instance:

print(carOfTom.lightsOn)    # will produce the output False because right now this instance variable still has its default value

The rest of the class definition consists of the methods for performing certain actions with a Car object. You can see that the already mentioned methods for changing the state of the Car object are very simple. They just assign a new value to the respective instance variable, a new speed value that is provided as a parameter in the case of changeCurrentSpeed(…) and a fixed Boolean value in the cases of turnLightsOn() and turnLightsOff(). In addition, we added a method printInfo() that prints out a string with the values of all instance variables to provide us with all information about a car’s current state. Let us now create a new instance of our Car class and then use some of its methods:

carOfSue = Car() 
carOfSue.owner = 'Sue' 
carOfSue.color = 'white' 
carOfSue.changeCurrentSpeed(41) 
carOfSue.turnLightsOn() 
carOfSue.printInfo()

Output: 

Car with owner = Sue, color = white, currentSpeed = 41, lightsOn = True

Since we did not define any methods to change the owner or color of the car, we are directly accessing these instance variables and assigning new values to them in lines 2 and 3. While this is okay in simple examples like this, it is recommended that you provide so-called getter and setter methods (also called accessor and mutator methods) for all instance variables that you want the user of the class to be able to read (“get”) or change (“set”). The methods allow the class to perform certain checks to make sure that the object always remains in an allowed state. How about you go ahead and for practice create a second car object for your own car (or any car you can think of) in a new variable and then print out its information?

A method can call any other method defined in the same class by using the notation “self.<name of the method>(...)”. For example, we can add the following method randomSpeed() to the definition of class Car:

def setRandomSpeed(self): 
    self.changeCurrentSpeed(random.randint(0,76))

The new method requires the “random” module to be imported at the beginning of the script. The method generates a random number and then uses the previously defined method changeCurrentSpeed(…) to actually change the corresponding instance variable. In this simple example, one could have simply changed the instance variable directly but in more complex cases changes to the state can require more code so that this approach here actually avoids having to repeat that code. Give it a try and add some lines to call this new method for one of the car objects and then print out the info again.

It can be a bit cumbersome to use methods or assignments to set all the instance variables to the desired initial values after a new object has been created. Instead, one would rather like to pass initial values to the constructor and get back an object with these values for the instance variables. It is possible to do so in Python by adding additional parameters to the constructor. Go ahead and change the definition of the constructor in class Car to the following version:

def __init__(self, owner = 'UNKNOWN', color = 'UNKNOWN', currentSpeed = 0, lightsOn = False): 
	self.owner = owner 
	self.color = color 
	self.currentSpeed = currentSpeed 
	self.lightsOn = lightsOn
    

Please note that we here used identical names for the instance variables and corresponding parameters of the constructor used for providing the initial values. However, these are still distinguishable because instance variables always have the prefix “self.”. In this new version of the constructor we are using keyword arguments for each of the properties to provide maximal flexibility to the user of the class. The user can now use any combination of providing their own initial values or using the default values for these properties. Here is how to re-create Sue’s car by providing values for all the properties:

carOfSue = Car(owner='Sue', color='white', currentSpeed = 41, lightsOn = True) 
carOfSue.printInfo()

Output: 

Car with owner = Sue, color = white, currentSpeed = 41, lightsOn = True

Here is a version in which we only specify the owner and the speed. Surely you can guess what the output will look like.

carOfSue = Car(owner='Sue', currentSpeed = 41) 
carOfSue.printInfo() 

In addition to __init__(…) for the constructor, there is another special method called __str__(). This method is called by Python when you either explicitly convert an object from that class to a string using the Python str(…) function or implicitly, e.g. when printing out the object with print(…). Try out the following two commands for Sue’s car and see what output you get:

print(str(carOfSue)) 
print(carOfSue)

Now add the following method to the definition of class Car:

	def __str__(self): 
            return 'Car with owner = {0}, color = {1}, currentSpeed = {2}, lightsOn = {3}'.format(self.owner, self.color, self.currentSpeed, self.lightsOn)

Now repeat the two commands from above and look at the difference. The output should now be the following line repeated twice:

Car with owner = Sue, color = UNKNOWN, currentSpeed = 41, lightsOn = False

For implementing the method, we simply used the same string that we were printing out from the printInfo() method. In principal, this method is not really needed anymore now and could be removed from the class definition.

Objects can be used like any other value in Python code. Actually, everything in Python is an object, even primitive data types like numbers and Boolean values. That means we can …

• use objects as parameters of functions and methods (you will see an example of this with the function stopCar(…) defined below),
• return objects as the return value of a function or method,
• store objects inside sequences or containers, for instance in lists like this: carList = [ carOfTom, carOfSue, Car(owner = 'Mike'],
• store objects in instance variables of other objects.

To illustrate this last point, we can add another class to our car example, one for representing car manufacturers:

class Manufacturer(): 
     def __init__(self, name): 
          self.name = name

Usually such a class would be much more complex, containing additional properties for describing a concrete car manufacturer. But we keep things very simple here and say that the only property is the name of the manufacturer. We now modify the beginning of the definition of class Car so that another instance variable is created called self.manufacturer. This is used for storing an object of class Manufacturer inside each Car object for representing the manufacturer of that particular car. For parameters that are objects of classes, it is common to use the special value None as the default value when the parameter is not provided.

class Car(): 

     def __init__(self, manufacturer = None, owner = 'UNKNOWN', color = 'UNKNOWN', currentSpeed = 0, lightsOn = False): 
        self.manufacturer = manufacturer 
        self.owner = owner 
        self.color =  color 
        self.currentSpeed = currentSpeed 
        self.lightsOn = lightsOn

The rest of the class definition can stay the same although we would typically change the __str__(...) method to include this new instance variable. The following code shows how to create a new Car object by first creating a Manufacturer object with name 'Chrysler'. This object could also come from a predefined list or dictionary of car manufacturer objects if we want to be able to use the same Manufacturer object for several cars. Then we use this object for the manufacturer keyword argument of the Car constructor. As a result, this object gets assigned to the manufacturer instance variable of the car as reflected by the output from the final print statement.

m = Manufacturer('Chrysler') 
carOfFrank = Car(manufacturer = m, owner = 'Frank', currentSpeed = 70) 
print(carOfFrank.manufacturer.name)

Output: 
Chrysler

Note how in the last line of the example above, we chain things together via dots starting from the variable containing the car object (carOfFrank), followed by the name of an instance variable (manufacturer) of class Car, followed by the name of an instance variable of class Manufacturer (name): carOfFrank.manufacturer.name . This is also something you have probably seen before, for instance as “describeObject.SpatialReference.Name” when accessing the name of the spatial reference object that is stored inside an arcpy Describe object.

We briefly discussed in Section 4.2 when talking about collections that when defining our own classes we may have to provide definitions of comparison operators like == and < for them to work as we wish when placed into a collection. So a question for instance would be, when should two car objects be considered to be equal? We could take the standpoint that they are equal if the values of all instance variables are equal. Or it could make sense for a particular application to define that two Car objects are equal if the name of the owner and the manufacturer are equal. If our instance variables would include the license plate number that would obviously make for a much better criterion. Similarly, let us say we want to keep our Car objects in a priority queue sorted by their current speed values. In that case, we need to define the < comparison operator so that car A < car B holds if the value of the currentSpeed variable of A is smaller than that of B.

The meaning of the == operator is defined via a special method called __eq__(…) for “equal”, while that of the < operator is defined in a special method called __lt__(…) for “less than”. The following code example extends the most recent version of our class Car with a definition of the __eq__(…) method based on the idea that cars should be treated as equal if owner and manufacturer are equal. It then uses a Python list with a single car object and another car object with the same owner and manufacturer but different speed to illustrate that the new definition works as intended for the list operations “in” and index(…).

class Car(): 

    … # just add the method below to the previous definition of the class 

    def __eq__(self, otherCar): 
        return self.owner == otherCar.owner and self.manufacturer == otherCar.manufacturer 

m = 'Chrysler' 
carList = [ Car(owner='Sue', currentSpeed = 41, manufacturer = m) ] 
car = Car(owner='Sue', currentSpeed = 0, manufacturer = m) 

if car in carList: 
    print('Already contained in the list') 
print(carList.index(car))

Output: 
Already contained in the list 
0

Note that __eq__(…) takes another Car object as parameter and then simply compares the values of the owner and manufacturer instance variables of the Car object the method was called for with the corresponding values of that other Car object. The output shows that Python considers the car to be already located in the list as the first element, even though these are actually two different car objects with different speed values. This is because these operations use the new definition of the == operator for objects of our class Car that we provided with the method __eq__(...).

You now know the basics of writing own classes in Python and how to instantiate them and use the created objects. To wrap up this section, let’s come back to a topic that we already discussed in Section 1.4 of Lesson 1. Do you remember the difference between mutable and immutable objects when given as a parameter to functions? Mutable objects like lists used as parameters can be changed within the function. All objects that we create from classes are also mutable, so you can in principle write code like this:

def stopCar(car): 
	car.currentSpeed = 0 

stopCar(carOfFrank) 
print(carOfFrank)

When stopCar(...) is called, the parameter car will refer to the same car object that variable carOfFrank is referring to. Therefore, all changes made to that object inside the function referring to variable car will be reflected by the final print statement for carOfFrank showing a speed of 0. What we have not discussed so far is that there is a second situation where this is important, namely when making an assignment. You may think that when you write something like

anotherCar = carOfFrank 

a new variable will be created and a copy of the car object in variable carOfFrank will be assigned to that variable so that you can make changes to the instance variables of that object without changing the object in carOfFrank. However, that is only how it works for immutable values. Instead, after the assignment, both variables will refer to the same Car object in memory. Therefore, when you add the following commands

anotherCar.color = 'green' 
anotherCar.changeCurrentSpeed(12) 
print(carOfFrank)

The output will be:

Car with owner = Frank, color = green, currentSpeed = 12, lightsOn = False

It works in the same way for all mutable objects, so also for lists for example. If you want to create an independent copy of a mutable object, the module copy from the Python standard library contains the functions copy(…) and deepcopy(…) to explicitly create copies. The difference between the two functions is explained in the documentation and only plays a role when the object to be copied contains other objects, e.g. if you want to make a copy of a list of Car objects.

Let’s change our example so that both Circle and Rectangle are derived from such a general class called Geometry. This class will be an abstract class in the sense that it is not intended to be used for creating objects from. Its purpose is to introduce properties and templates for methods that all geometric classes in our project have in common.

class Geometry(): 

	def __init__(self, x = 0.0, y = 0.0): 
		self.x = x 
		self.y = y 

	def computeArea(self): 
    	pass 

    def computePerimeter(self): 
		pass 

    def move(self, deltaX, deltaY): 
		self.x += deltaX 
		self.y += deltaY 

    def __str__(self): 
		return 'Abstract class Geometry should not be instantiated and derived classes should override this method!'

The constructor of class Geometry looks pretty normal, it just initializes the instance variables that all our geometry objects have in common, namely x and y coordinates to describe their location in our 2D coordinate system. This is followed by the definitions of the methods computeArea(), computePerimeter(), move(…), and __str__() that all geometry objects should support. For move(…), we can already provide an implementation because it is entirely based on the x and y instance variables and works in the same way for all geometry objects. That means the derived classes for Circle and Rectangle will not need to provide their own implementation. In contrast, you cannot compute an area or perimeter in a meaningful way just from the position of the object. Therefore, we used the keyword pass to indicate that we are leaving the body of the computeArea() and computePerimeter() methods intentionally empty. These methods will have to be overridden in the definitions of the derived classes with implementations of their specialized behavior. We could have done the same for __str__() but instead we return a warning message that this class should not have been instantiated.

It is worth mentioning that, in many object-oriented programming languages, the concepts of an abstract class (= a class that cannot be instantiated) and an abstract method (= a method that must be overridden in every subclass that can be instantiated) are built into the language. That means there exist special keywords to declare a class or method to be abstract and then it is impossible to create an object of that class or a subclass of it that does not provide an implementation for the abstract methods. In Python, this has been added on top of the language via a module in the standard library called abc (for abstract base classes). Although we won’t be using it in this course, it is a good idea to check it out and use it if you get involved in larger Python projects. This Abstract Classes page is a good source for learning more.

Here is our new definition for class Circle that is now derived from class Geometry. We also use a few commands at the end to create and use a new Circle object of this class to make sure everything is indeed working as before:

import math  

class Circle(Geometry): 

	def __init__(self, x = 0.0, y = 0.0, radius = 1.0): 
		super(Circle,self).__init__(x,y) 
		self.radius = radius 

	def computeArea(self): 
		return math.pi * self.radius ** 2 

	def computePerimeter (self): 
		return 2 * math.pi * self.radius 

	def __str__(self): 
		return 'Circle with coordinates {0}, {1} and radius {2}'.format(self.x, self.y, self.radius) 

circle1 = Circle(10, 10, 10) 
print(circle1.computeArea()) 
print(circle1.computePerimeter()) 
circle1.move(2,2) 
print(circle1)

Here are the things we needed to do in the code:

• In line 3, we had to change the header of the class definition to include the name of the base class we are deriving Circle from (‘Geometry’) within the parentheses.
• The constructor of Circle takes the same three parameters as before. However, it only initializes the new instance variable radius in line 7. For initializing the other two variables it calls the constructor of its base class, so the class Geometry, in line 6 with the command “super(Circle,self).__init__(x,y)”. This is saying “call the constructor of the base class of class Circle and pass the values of x and y as parameters to it”. It is typically a good idea to call the constructor of the base class as the first command in the constructor of the derived class so that all general initializations are taken care off.
• Then we provide definitions of computeArea() and computePerimeter() that are specific for circles. These definitions override the “empty” definitions of the Geometry base class. This means whenever we invoke computeArea() or computePerimeter() for an object of class Circle, the code from these specialized definitions will be executed.
• Note that we do not provide any definition for method move(…) in this class definition. That means when move(…) will be invoked for a Circle object, the code from the corresponding definition in its base class Geometry will be executed.
• We do override the __str__() method to produce the same kind of string with information about all instance variables that we had in the previous definition. Note that this function accesses both the instance variables defined in the parent class Geometry as well as the additional one added in the definition of Circle.

The new definition of class Rectangle, now derived from Geometry, looks very much the same as that of Circle if you replace “Circle” with “Rectangle”. Only the implementations of the overridden methods look different, using the versions specific for rectangles.

class Rectangle(Geometry): 

	def __init__(self, x = 0.0, y = 0.0, width = 1.0, height = 1.0): 
		super(Rectangle, self).__init__(x,y) 
		self.width = width 
        self.height = height 

	def computeArea(self): 
		return self.width * self.height 

	def computePerimeter (self): 
		return 2 * (self.width + self.height) 

	def __str__(self): 
		return 'Rectangle with coordinates {0}, {1}, width {2} and height {3}'.format(self.x, self.y, self.width, self.height ) 

rectangle1 = Rectangle(15,20,4,5) 
print(rectangle1.computeArea()) 
print(rectangle1.computePerimeter()) 
rectangle1.move(2,2) 
print(rectangle1)

Overall, the new definitions of Circle and Rectangle have gotten shorter and redundant code like the implementation of move(…) only appears once, namely in the most general class Geometry. Let’s add another class to the hierarchy, a class for axis-aligned Square objects. Of course, you could argue that our class Rectangle is already sufficient to represent such squares. That is correct but we want to illustrate how it would look if you specialize a class already derived from Geometry further and one could well imagine a more complex version of our toy GIS example in which squares would add some other form of specialization. The resulting class hierarchy will then look like in the image below. The new class Square is a derived class of class Rectangle (so Rectangle is its base class) but it is also indirectly derived from class Geometry. Therefore, we say both Geometry and Rectangle are superclasses of Square and Square is a subclass of both these classes. Please note that the way we have been introducing these terms here, the terms base and derived class desribe the relationship between two nodes directly connected by a single arrow in the hierarchy graph, while superclass and subclass are more general and describe the relationship between two classes that are connected via any number of directed arrows in the graph.

Figure 4.18 Class hierarchy after adding class Square. Square is a subclass of both Geometry and Rectangle. Both these classes are superclasses of Square.

Here is the code for class Square:

class Square(Rectangle): 

	def __init__(self, x = 0.0, y = 0.0, sideLength = 1.0): 
		super(Square,self).__init__(x, y, sideLength, sideLength) 

	def __str__(self): 
		return 'Square with coordinates {0}, {1} and sideLength {2}'.format(self.x, self.y, self.width ) 

square1 = Square(5, 5, 8) 
print(square1.computeArea()) 
print(square1.computePerimeter()) 
square1.move(2,2)
print(square1)

Right, the definition of Square is really short; we only define a new constructor that only takes x and y coordinates and a single sideLength value rather than width and height values. In the constructor we call the constructor of the base class Rectangle and provide sideLength for both the width and height parameters of that constructor. There are no new instance variables to initialize, so this is all that needs to happen in the constructor. Then the only other thing we have to do is override the __str__() method to produce some square-specific output message using self.width for the side length information for the square. (Of course, we could have just as well used self.height here.) The implementations of methods computeArea() and computePerimeter() are inherited from class Rectangle and the implementation of move(…) indirectly from class Geometry.

Now that we have this class hierarchy consisting of one abstract and three instantiable classes, the following code example illustrates the power of polymorphism. Imagine that in our toy GIS we have created a layer consisting of objects of the different geometry types. If we now want to implement a function computeTotalArea(…) that computes the combined area of all the objects in a layer, this can be done like this:

layer = [ circle1, rectangle1, square1, Circle(3,3,9),  Square(30, 20, 5) ] 

def computeTotalArea(geometryLayer): 
	area = 0 
	for geom in geometryLayer: 
		area +=  geom.computeArea() 
	return area 

print(computeTotalArea(layer))

Output: 
677.6282702997526

In line 1, you see how we can create a list of objects of the different classes from our hierarchy to represent the layer. We included objects that we already created previously in variables circle1, rectangle1, and square1 but also added another Circle and another Square object that we are creating directly within the square brackets […]. The function computeTotalArea(…) then simply takes the layer list, loops through its elements, and calls computeArea() for each object in the list. The returned area values are added up and returned as the total area.

The code for this is really compact and elegant without any need for if-else to realize some case-distinction based on the geometry type of the given object in variable geom. Let’s further say we would like to add another class to our hierarchy, a class Polygon that – since polygons are neither specialized versions of circles or rectangles – should be derived from the root class Geometry. Since polygons are much more complex than the basic shapes we have been dealing with so far (e.g. when it comes to computing their area), we will not provide a class definition here. But, once we have written the class, we can include polygons in the layer list from the previous example …

layer = [  Polygon(…), circle1, rectangle1, square1, Circle(3,3,9),  Square(30, 20, 5) ]

… and the code for computing the total area will immediately work without further changes. All changes required for making this addition are nicely contained within the class definition of Polygon because of the way inheritance and polymorphism are supported in Python.

Understanding this idea of building reusable GUI components via inheritance is so important that we should look at another example. While doing so, we will also learn how you can actually programmatically draw on a widget to display your own content. What we are going to do is take the classes from our Geometry hierarchy from the previous section and create a widget that actually draws the instances of the classes we have stored in a list to the screen. To make this a bit more interesting, we also want all the objects of the different geometry types to have a “color” attribute that determines in which color the object should be drawn. Before we look at what changes need to be made to the different geometry classes, here is a quick introduction to drawing with PyQt5.

Every widget in QT5 has a method called paintEvent(…) that is called when the widget needs to be drawn (for instance, when its drawn for the first time or when the size of the widget has changed). The only parameter passed to this method is an event object that can be used to get the current dimensions of the content area that we can draw on by calling its rect() method. That means when we want to use a widget for drawing something on it, we derive a new class from the respective widget class and override the paintEvent(…) method with our own implementation that takes care of the drawing. To do the actual drawing, we need to create an object of the class QPainter and then use the drawing methods it provides. Here is a simple example; the details will be explained below:

import sys 
from PyQt5 import QtGui, QtWidgets 
from PyQt5.QtCore import Qt, QPoint 

class MyWidget(QtWidgets.QWidget): 

    def paintEvent(self, event): 
        qp = QtGui.QPainter() 
        qp.begin(self) 

        qp.setPen(QtGui.QColor(200,0,0)) 
        qp.drawText(20,20, "Text at fixed coordinates") 
        qp.drawText(event.rect(), Qt.AlignCenter, "Text centered in the drawing area") 
        qp.setPen(QtGui.QPen(Qt.darkGreen, 4)) 
        qp.drawEllipse(QPoint(50,60),30,30) 
        qp.setPen(QtGui.QPen(Qt.blue, 2, join = Qt.MiterJoin)) 
        qp.drawRect(20,60,50,80) 

        qp.end() 

app = QtWidgets.QApplication(sys.argv) 

window = MyWidget() 
window.show() 

sys.exit(app.exec_())

When you run this small script, you should see the following window on your screen:

     

Fig 4.19 Window produced by the previous code example illustrating some drawing methods of QPainter

Let’s look at the coarse structure first: we are defining a new class derived from QWidget and only overriding the paintEvent(…) method, meaning in all other aspects this widget will behave like an instance of QWidget. In the main code, we simply create an instance of our new widget class and make it show up on the screen. Now, let’s look at the body of method paintEvent(…): The first thing to note here is that all drawing needs to be preceded by the creation of the QPainter object (line 8) and the call of its begin(…) method using “self” as the parameter standing for the widget object itself because that is what we want to draw on (line 9). To conclude the drawing, we need to call the end() method of the QPainter object (line 19).

Next, let us look at the methods of the QPainter object we are invoking that all start with “draw…”. These are the methods provided by QPainter to draw different kinds of entities like text, circles or ellipses, rectangles, images, etc. We here use the method drawText(…) twice to produce the two different lines of text (lines 12 and 13). The difference between the two calls is that in the first one we use absolute coordinates, so the text will be drawn at pixel coordinates 20, 20 counting from the top left corner of the widget’s content area. The second call takes a rectangle (class QRectF) as the first parameter and then draws the text within this rectangle based on the additional text options given as the second parameter which here says that the text should be centered within the rectangle. This is an example where a class provides several methods with the same name but different parameters, something that is called overloading. If you check out the documentation of QPainter, you will see that most methods come in different versions. Now go ahead and resize the window a bit and see how the text produced by the first call always remains at the same absolute position, while that from the second call always stays centered within the available area.

In line 15, we use the method drawEllipse(…) to produce the circle. There is no special circle drawing method, so we use this one and then provide the same number for the two radii. To draw the rectangle, we use the method drawRect(…) in the version that takes the coordinates of the corner plus width and height values as parameters.

The remaining calls of methods of the QPainter object are there to affect the way the objects are drawn, e.g. their color. Colors inPyQt5 are represented by instances of the class QColor. In line 11, we create a new QColor object by providing values between 0 and 255 for the color’s red, green, and blue values. Since the red value is 200 and both green and blue are zero, the overall color will be the kind of red that the text appears in. QT5 also has a number of predefined colors that we are using in lines 14 (Qt.darkGreen) and 16 (Qt.blue).

QPainter uses objects of class QPen and QBrush to draw the boundary and inside of a shape. In line 11, it is stated that a pen with red color should be used for the following drawing operations. As a result, both text lines appear in red. In line 14, we create a new QPen object to be used by the QPainter and specify that the color should be dark green, and the line width should be 4. This is used for drawing the circle. In line 16, we do the same with color blue and line width 2, and, in addition, we say that sharp corners should be used for the connection between to adjacent line segments of the shape’s border. This is used for drawing the rectangle. We won’t go further into the details of the different pen and brush properties here but the documentation of the QPen and QBrush classes provides some more examples and explanations. In addition, you will see more use cases in the walkthrough in the next section.

We are now going to revist the geometry example we saved in Section 4.8.

To prepare our geometry classes to be drawn on the screen, we first need to modify their definitions by

1. introducing the new “color” instance variable for storing a QColor object,
2. defining a method called “paint” for drawing the respective object on a screen. This method will have one parameter called “painter” for passing an object of the QPainter class that can be used for drawing onto the widget.

As an exercise, think about where in the class hierarchy you would need to make changes to address points (1) and (2). Once you have thought about this for a bit, read on.

The new attribute “color” is something that all our geometry classes have in common. Therefore, the best place to introduce it is in the root class Geometry. In all subclasses (Circle, Rectangle, Square), you then only have to adapt the constructor to include an additional keyword parameter for the color. Regarding point (2): As we saw above, drawing a circle or a rectangle requires different methods of the QPainter object to be called with different kinds of parameters specific to the particular geometry type. Therefore, we define the method paint(…) in class Geometry but then override it in the subclasses Circle and Rectangle. For class Square, the way it is based on class Rectangle allows us to directly use the implementation of paint(…) from the Rectangle class, so we do not have to override the method in the definition of Square. Here are the changes made to the four classes.

Class Geometry: as discussed, in class Geometry we introduce the “color” variable, so we need to change the constructor a bit. In addition, we add the method paint(…) but without an implementation. The rest of the definition remains unchanged:

class Geometry: 

	def __init__(self, x = 0.0, y = 0.0, color = Qt.black): 
       	 self.x = x 
         self.y = y 
         self.color = color 

    ...

    def paint(self, painter): 
         pass 

Classes Circle and Rectangle: for the classes Circle and Rectangle, we adapt the constructors to also include a keyword argument for “color”. The color is directly passed on to the constructor of the base class, while the rest of the constructor remains unchanged. We then provide a definition for method paint(…) that sets up the Pen object to use the right color and then uses the corresponding QPainter method for drawing the object (drawEllipse(…) for class Circle and drawRect(…) for class Rectangle) providing the different instance variables as parameters. The rest of the respective class definitions stay the same:

class Circle (Geometry): 

    def __init__(self, x = 0.0, y = 0.0, radius = 1.0, color = Qt.black): 
         super(Circle,self).__init__(x,y,color) 
         self.radius = radius 

    ...

    def paint(self, painter): 
         painter.setPen(QtGui.QPen(self.color, 2)) 
         painter.drawEllipse(QPoint(self.x, self.y), self.radius, self.radius) 


class Rectangle(Geometry): 

    def __init__(self, x = 0.0, y = 0.0, width = 1.0, height = 1.0, color = Qt.black): 
         super(Rectangle,self).__init__(x,y, color) 
         self.width = width 
         self.height = height 

    ...

    def paint(self, painter): 
         painter.setPen(QtGui.QPen(self.color, 2,  join = Qt.MiterJoin)) 
         painter.drawRect(self.x, self.y, self.width, self.height)

Class Square: sor the Square class, we just adapt the constructor to include the color, the rest remains unchanged:

class Square(Rectangle): 

    def __init__(self, x = 0.0, y = 0.0, sideLength = 1.0, color = Qt.black): 
         super(Square,self).__init__(x, y, sideLength, sideLength, color)

Now that we have the geometry classes prepared, let us again derive a specialized class GeometryDrawingWidget from QWidget that stores a list of geometry objects and in the paintEvent(…) method sets up a QPainter object for drawing on its content area and then invokes the paint(…) methods for all objects from the list. This is bascially the same thing we did towards the end of Section 4.7.2 but now happens inside the new widget class. The list of objects is supposed to be given as a parameter to the constructor of GeometryDrawingWidget:

import math, sys 
from PyQt5 import QtGui, QtWidgets 
from PyQt5.QtCore import Qt, QPoint

class GeometryDrawingWidget(QtWidgets.QWidget): 

     def __init__(self, objects): 
         super(GeometryDrawingWidget,self).__init__() 
         self.objectsToDraw = objects 

     def paintEvent(self, event): 
         qp = QtGui.QPainter() 
         qp.begin(self) 
         for obj in self.objectsToDraw: 
             obj.paint(qp) 
         qp.end()

Finally, we create another new class called MyMainWindow that is derived from QMainWindow for the main window containing the drawing widget. The constructor takes the list of objects to be drawn and then creates a new instance of GeometryDrawingWidget passing the object list as a parameter and makes it its central widget in line 6:

class MyMainWindow(QtWidgets.QMainWindow): 

    def __init__(self, objects): 
	    super(MyMainWindow, self).__init__() 
	    self.resize(300,300) 
	    self.setCentralWidget(GeometryDrawingWidget(objects))

In the main code of the script, we then simply create an instance of MyMainWindow with a predefined list of geometry objects and then call its show() method to make the window appear on the screen:

app = QtWidgets.QApplication(sys.argv) 

objects = [ Circle(93,83,45, Qt.darkGreen), Rectangle(10,10,80,50, QtGui.QColor(200, 0, 250)), Square(30,70,38, Qt.blue)] 

mainWindow = MyMainWindow (objects) 
mainWindow.show() 

sys.exit(app.exec_())

When you run the program, the produced window should look like in the figure below. While this is not really visible, every time you resize the window, the paintEvent(…) method of GeometryDrawingWidget will be called and the content consisting of the three geometric objects will be redrawn. While in this simple widget, we only use fixed absolute coordinates so that the drawn content is independent of the size of the widget, one could easily implement some logic that would scale the drawing based on the available space.

     

Figure 4.20 Window produced by the previous code example with the drawings of the different geometric objects

Since this project involves quite a bit of code, we have tried to cleanly organize it into different class definitions and multiple files. We don’t expect you to type in the code yourself in this walkthrough but rather study the files carefully and use the explanations provided in this text to make sure you understand how everything plays together. Here is an overview on what each of the involved files contains:

core_classes.py – This file contains most of the basic classes for our project that are not derived from other classes.

• class Timepoint: a Timepoint in this project will represent a point in space and time and is therefore defined by latitude and longitude properties which are both instance variables of type float and a time property that is an instance of the class datetime from the datetime module of the Python standard library (that you worked with in homework assignment 3). We will use Timepoints to represent the bus observations from the GPS data as well as events happening at a particular location at a particular time.
• class Bus: the class Bus represents a single vehicle appearing in the GPS data and its GPS points extracted from the input data. It is defined by its vehicle ID, the line it belongs to, and a list of chronologically-ordered Timepoints representing the GPS observation points read from the input data.
• class BusTracker: this class represents the current status of an individual bus vehicle during the analysis and event detection process. A BusTracker object will be maintained for each bus vehicle. It links to the corresponding object of class Bus and, in addition, is used to keep track of whether the vehicle’s status is currently “driving”, “stopped”, or “in depot”, of its current speed estimate, of the last Timepoint that has already been processed for this bus, and of some other information.
• class Observation: in the analysis stage of our program, we will maintain a priority queue (Section 4.2.3) in which we store the next observations for each bus. These observations need to be processed in order of their time of occurrence. The objects we are storing in this queue will be of class Observation that combines the BusTracker for the bus this observation is about, the Timepoint of this observation, and the index of that Timepoint in the list of Timepoints of the bus. Please note that the class has to define its own __lt__(...) method for the < comparison operator to achieve that the Obervation objects in the priority queue will be ordered based on their time attribute of the Timepoint object they contain. 
• class Depot: the class Depot is used to represent a single bus depot, simply defined by its name and bounding box which is represented as a 4-tuple of floats as in the original input data in file dublin_depots.csv.

bus_events.py - This file defines the hierarchy of bus events starting with the abstract root class BusEvent from which we derive three more specialized (but still abstract) classes SingleBusEvent, MultipleBusesEvent, and BusDepotEvent. The classes for the events that we are actually trying to detect in the data are derived from these three intermediate classes. The overall bus event hierarchy is depicted in the figure below. We are keeping things somewhat simple here. One could certainly imagine other kinds of events that could be of interest and easily added to this hierarchy.

     

Figure 4.23 Class hierarchy of bus event classes used in this project

• class BusEvent: each bus event has an associated Timepoint. Therefore, the instance variable for this Timepoint object is already introduced in the abstract root class BusEvent. In addition, each bus event should have a method description() to produce a string with a brief description of the event suitable to be stored in an attribute table of a feature class, and a class function (Section 4.8) called detect(…) that takes an object of class Observation and checks whether the respective event occurs at this observation point and if so, creates and returns one or multiple event objects of that event class. We define templates for both these functions in BusEvent but they need to be overwritten by the derived subclasses.
• class SingleBusEvent: this class is intended as a superclass for all events involving just a single bus, so those that are entirely based on properties of the bus and its GPS data. The class adds an instance variable to refer to the Bus object this event is about.
• class MultipleBusesEvent: this class is the superclass for all bus events that involve two or more buses (like bus encounter events). It adds an instance variable for storing a list of Bus objects, namely those involved in the particular event.
• class BusDepotEvent: all event classes that involve a single bus and a bus depot are supposed to be derived from this intermediate class. It adds instance variables for the involved bus (class Bus) and involved depot (class Depot).
• class BusStoppedEvent: this is the first instantiable class in our hierarchy and it is derived from the SingleBusEvent class because it is only about the properties of an individual bus: when the given Timepoint on the bus doesn’t move more than 3 meters for at least a minute, we consider the bus to be stopped, generate an event of this class, and update the status of the BusTracker for this bus accordingly. This all happens in the function detect(…) that is overwritten in the definition of class BusStoppedEvent. This class adds another instance variable called duration for keeping track of how long the bus remained stopped.
• class BusEncounterEvent: this event derived from class MultipleBusesEvent represents situations in which two buses encounter each other along their routes. The detect(…) function overwritten in the class definition creates an event of this class if the given bus is observed within 20 meters of the last accounted position of another bus while both buses are “driving” (so not stopped and not currently in a depot).
• class LeavingDepotEvent: this class is derived from BusDepotEvent and an event of this class will be generated (again by the overwritten function detect(…)) if the bus still has the status “in depot” but for the given observation now is located outside any depot. The status of the corresponding BusTracker will be changed from “in depot” to “driving”.
• class EnteringDepotEvent: this is the counter part to LeavingDepotEvent for when a bus that still has the status “driving” now for the current observation is located inside one of the depots.

bus_track_analyzer.py – This file contains just a single class definition, the definition of class BusTrackAnalyzer that is our main class for performing the analysis and event detection over the data read in from the two input files. Its constructor takes two input parameters: a dictionary that maps a bus vehicle ID to the corresponding object of class Bus (created from the data from the GPS input file) and a list of Depot objects (created from the data in the depot input file). While its code could also have become the main program for this project, it is advantageous to have this all encapsulated into a class definition with methods for performing a single step of the analysis, resetting the analysis to start from the beginning, and for producing output vector data sets of the bus tracks created and the events detected so far. This way we can use this analyzer differently in different contexts and have full control over when and in which order the individual analysis steps and other actions will be performed. When we turn the project into a QGIS plugin in Section 4.12, we will make use of this by linking the methods of this class to a media player like control GUI with buttons for starting, pausing, and resetting the analysis. We will explain how this main class of the project works in more detail in a moment.

bus_tracker_widget.py – This file also defines just a single class, BusTrackerWidget, which is for visualizing the current status of the buses during the analysis and, therefore, is a bit of an optional component in this project to practice what you learned about drawing on a widget some more in this lesson. A BusTrackerWidget object is directly linked to a BusTrackAnalyzer object that is given to it as a parameter to the constructor. Whenever the content is supposed to be drawn, it accesses the analyzer object and, in particular, the list of BusTracker objects maintained there and depicts the status of the different buses as shown in the image below with each line representing one of the buses:

     

Figure 4.24 GUI for the event detection process implemented in class BusTrackerWidget

The class uses the three images from the files status_driving.png, status_indepot.png, and status_stopped.png to show the current status as an icon in the leftmost column of each row. This is followed by a colored circle depicting the vehicles current estimated speed using red for speeds below 15 mph, orange below 25 mph, yellow below 40 mph, and green for speeds larger than 40 mph. Then it displays the bus ID and line information followed by a short text description providing more details on the status like the exact speed estimate or the name of the depot the bus is currently located in. The widget also shows the time of the last processed observation in green at the top. We will discuss how this class has been implemented in more detail later in this section.

main.py – Lastly, this file contains the main program for this project in which we put everything together. Since most of the analysis functionality is implemented in class BusTrackAnalyzer and the detect(…) functions of the different bus event classes, this main program is comparatively compact. It reads in the data from the two input files, creates a BusTrackAnalyzer for the resulting Bus and Depot objects, and sets up a QWidget for a window that hosts an instance of the BusTrackerWidget class. When the main button in this QWidget is pressed, the run method is executed and processes the data step-by-step by iteratively calling the nextStep() method of the analyzer object until all observations have been processed. The method also makes sure that the BusTrackerWidget is repainted after each step so that we can see what is happening during the analysis. Finally, it saves the detected events and bus tracks as vector data sets on the disk.

Now that you have a broad overview of the classes and files involved, let’s look at the code in more detail. The code is too long to explain every line and there are parts that should be easy to understand with the knowledge you have now and the comments included in the source code, so we will only be picking out the main points. Nevertheless, please make sure to study carefully the class definitions and how they work together, and if something is unclear, please ask on the forums.

The plugin builder 3 plugin should already be installed in your QGIS version. If not, please go back to Section 4.4.1 and follow the installation instructions there. Now run the plugin by going Plugins -> Plugin Builder in the main menu.

To create a template for your plugin, you have to work through the dialog boxes of Plugin Builder and fill out the information there. Clicking the “Help” button will open a local .html page with detailed information on the purpose and meaning of the different fields. We fill out the first page as shown in the figure below. Here is a brief overview:

• Class name: this is the name of the main class for our plugin. It needs to be a valid Python class name.
• Plugin name: this is the name of the plugin in readable form. In contrast to the class name, it can contain spaces and other special characters.
• Description: a short description of what the plugin is for.
• Module name: the name of the .py file that will contain the main class for your plugin. Typically, the class name written with underscores rather than in CamelCase is used here.
• Version number & Minimum QGIS version: these are used to specify a version number for your plugin and the minimum QGIS version required to run it.
• Author/Company & Email address: here you provide information about the author of the plugin that will be used to create the copyright information at the beginning of the different files. We have no intentation of publishing this plugin, so we just use “489” here.

     

Figure 4.29 Filling out the first dialog page of Plugin Builder

On the next page, you can enter a longer description of your plugin. Since this is just a toy example, we don’t bother with this here and leave the text as it is.

On the next page, you can choose between different kinds of templates for your plugin, e.g. a simple dialog box or a dock widget, meaning an actual panel that can be docked and moved inside the QGIS GUI like the other panels, e.g. the Layers panel. We here will go with the dialog option. In the next section, we will then use the dock widget option. With the “Text for the menu item” option, we specify which text should show up in the menu for our plugin in the main menu bar. With the “Menu” option we pick in which menu of the menu bar this entry should be located. We will fill out this page as shown below:

     

Figure 4.30 Picking the template for our plugin

The checkboxes on the next page allow for determining which files Plugin Builder is going to create. It’s ok to leave all options checked.

The next page is specifying information that is relevant if you plan to publish your plugin, e.g. on Github. Since we are not planning this, we just leave the page unchanged.

On the last page, we can determine in which folder the new folder with the files for our plugin will be created. By default this is the default plugin folder of our QGIS installation, meaning the plugin will immediately be listed in QGIS when we start it next. If, instead of the path, you just see a dot (.), please browse to the plugins folder yourself, replacing the part directly after "C:/Users/" with your Windows user name. It is possible that the “AppData” folder in your user’s home directory is not visible in which case you will have to change your settings to show hidden files.

     

Figure 4.31 Final dialog page of Plugin Builder for picking the path where the folder for the plugin will be created

We now click the “Generate” button and Plugin Builder will then create the folder and different files for our plugin. It's possible that you will get a warning message that about Plugin Builder not being able to compile the resources.qrc file; that's ok, we will take care of that in a moment. Plugin Builder will now show us a summary page like in the figure below with some valuable information about where the plugin has been created and what the next steps should be. Even though we won’t be following these exactly, it’s a good idea to take a screenshot of this information or note down where the folder for our plugin is located. You may also want to add the plugins folder under “Quick access” in your Windows File Explorer since we will need to access it quite a lot in the remainder of this lesson.

     

Figure 4.32 Plugin Builder summary page after creating the new plugin

If we now open the new folder “random_number_generator” in the QGIS default plugin folder, you will see the following file structure:

     

Figure 4.33 Content of the folder created by Plugin Builder for our new plugin 

In this introduction to QGIS plugin development, we won’t go into the details of features like internationalization and test code generation, so you can ignore the different subfolders and also some of the files in the main folder. The important files are:

• random_number_generator.py: this file defines the main class RandomNumberGenerator for our plugin. If you look at the code that Plugin Builder 3 generated for this class, you will see that this class has to define a number of methods to interact with QGIS and its plugin interface. One central method is the method run() that in this case here shows our dialog box stored in instance variable self.dlg. If we wanted some code to be executed only if the dialog box is accepted (e.g. closed with the Ok button), then we would put that code here after the “if result:” line.
• random_number_generator_dialog.py: this is the class that defines the main GUI widget for our plugin derived from QDialog. We will talk about this class in more detail in a moment because this is the class that we are going to modify to implement what needs to happen when the user interacts with the GUI of our plugin.
• random_number_generator_dialog_base.ui: this is the class in which the GUI for our dialog box is defined. We will modify it in QT Designer in a moment to adapt the GUI for our purposes.
• __init__.py: A file with this name needs to be contained in the main folder of each plugin and it lets QGIS know about the plugin. QGIS calls the classFactory(…) function defined in this file to create an object of the main class of the plugin (class RandomNumberGenerator in this case).
• metadata.txt: this file contains most of the info about our plugin that we entered into Plugin Builder.

When you now restart QGIS and open the plugin manager under Plugins -> Manage and Install Plugins… , the Random Number Generator plugin should appear in the list of installed plugins but it still needs to be activated. However, if we try to activate it now, we will get an error message that no module called “randomnumbergenerator.resources” can be found. This is a file that we have to generate ourselves by compiling the file called resources.qrc located in our plugin folder with the pyrcc5 tool. The resources file contains information about all additional GUI related resources needed for our plugin, like additional icons for example.

Usually, you would edit the .qrc file first, e.g. in QT Designer, and then compile it. But we don’t need any additional resources for this project, so we just compile it directly with the following command in the OSGeo4W shell after first moving to our plugin folder with the cd command. (The screenshot below contains some commands that were needed in earlier versions of OSGeo4W; you will most likely only need the cd and pyrcc5 commands shown at the bottom.)

pyrcc5 resources.qrc –o resources.py

     

Figure 4.34 Commands for compiling the resource.qrc file in the OSGeo4W shell

After running this command, we can now activate the plugin by clicking on the checkbox:

There will now be a submenu called “Random Number Generator” with an entry with the same name in the Plugins menu:

In addition, a toolbar with a single button has been added for our plugin to the toolbar section at the top. Since we didn’t make any changes, the default button  is used. Either clicking the menu entry or the button will open the dialog box for our plugin:

     

Figure 4.35 Empty dialog box appearing after starting the new plugin in QGIS

Currently, the dialog box only contains the default elements, namely two buttons for accepting and rejecting the dialog. We are now going to change the GUI of the dialog box and add the random number generation functionality.

The file random_number_generator_dialog.py defines the class RandomNumberGeneratorDialog derived from QDialog. This is the widget class for our dialog box. The GUI itself is defined in the .ui file random_number_generator_dialog_base.ui. In contrast to first compiling the .ui file into a Python .py file and then using the Python code resulting from this to set up the GUI, this class uses the previously briefly mentioned approach of directly reading in the .ui file. That means to change the GUI, all we have to do is modify the .ui file in QT Designer, no compilation needed.

We open random_number_generator_dialog_base.ui with QT Designer and make the following small changes:

• We delete the button box with the Ok and cancel buttons.
• The layout of the dialog box is changed to a vertical layout.
• We add a push button and a label to the dialog box with vertical spacers in between.
• The button text is set to “Press here to generate random number between 1 and 100”.
• The label text is set to “The number is:” and its horizontal alignment property is changed to be horizontally centered.

The image below shows the produced layout in QT Designer.

     

Figure 4.36 Creating the GUI for our plugin in QT Designer

After saving the form, we want to see if we now get the new version of the GUI in QGIS as well. However, for this to happen, the plugin needs to be reloaded. For this, we use the Plugin Reloader plugin that we also installed at the beginning of the lesson. We pick the Plugins -> Plugin Reloader -> Configure option. Then we pick our random_number_generator plugin from the list and press OK.

     

Figure 4.37 Configuring Plugin Reloader to reload the code for our plugin

Now we just have to go Plugins -> Plugin Reloader -> Reload: random_number_generator, whenever we want to reload the plugin after making some changes. After doing so and starting our plugin again, it should show up like this in QGIS:

     

Figure 4.38 Running the plugin in QGIS now shows the dialog box with the GUI we just created

Now we just have to implement the functionality for our plugin, which for this toy example is extremely simple. We just need to write some code for rolling a random number between 1 and 100 and updating the text of the label accordingly when the button is pressed. We implement this in the definition of our class RandomNumberGeneratorDialog in file random_number_generator_dialog.py. The changes we make are the following:

• In the line where the os module is imported, we also import the random module from the Python standard library.
• We add a new method to the class definition with the event handler code:

  def generateNewNumber(self): 
  	r = random.randint(1,100) 
  	self.label.setText("The number is: " + str(r))

• Finally, we connect the button to this new method by adding the following line at the end of the __init__() method:

  self.pushButton.clicked.connect(self.generateNewNumber)

Below is the entire code part of random_number_generator_dialog.py after making these changes:

import os, random 

from PyQt5 import uic 
from PyQt5 import QtWidgets 

FORM_CLASS, _ = uic.loadUiType(os.path.join( os.path.dirname(__file__), 'random_number_generator_dialog_base.ui')) 
 
class RandomNumberGeneratorDialog(QtWidgets.QDialog, FORM_CLASS): 
    def __init__(self, parent=None): 
        """Constructor.""" 
        super(RandomNumberGeneratorDialog, self).__init__(parent) 
        # Set up the user interface from Designer. 
        # After setupUI you can access any designer object by doing 
        # self.<objectname>, and you can use autoconnect slots - see 
        # http://qt-project.org/doc/qt-4.8/designer-using-a-ui-file.html 
        # #widgets-and-dialogs-with-auto-connect 
        self.setupUi(self) 

        self.pushButton.clicked.connect(self.generateNewNumber) 

    def generateNewNumber(self): 
        r = random.randint(1,100) 
        self.label.setText("The number is: " + str(r))

Now the only thing we have to do is reload the plugin again in QGIS, and when we then start the plugin one more time, the button should work and allow us to create random numbers. Please note that the dialog box is opened non-modally from our main plugin class so that you can keep the dialog box open in the background while still interacting with the rest of the QGIS GUI.

     

Figure 4.39 The Random Number Generator plugin in QGIS after adding the main functionality

In this small example, all code changes needed were about wiring up the GUI of the dialog box. In such cases, modifying the file with the class definition for the main widget can be sufficient. In the next section, we will implement a plugin that is significantly more complex and interacts with QGIS in more sophisticated ways, e.g., by adding new layers to the currently open project and continuously updating the content of these layers. However, the general approach is the same in that we are mainly going to modify the .ui file and the file with the main widget class from the template files created by Plugin Builder 3. If you want to learn more about QGIS plugin development and the role the other files play, the official documentation is a good starting point.

To create a folder with a template version for this plugin, please follow the steps below.

1. Run the Plugin Builder plugin and fill out the first page as follows before clicking “Next >”:
   

   Figure 4.40 First Plugin Builder dialog page with information about the new plugin
2. Just click “Next >” again as we will skip the description part and then on the next page pick the “Tool button with dock widget” template at the top and fill out the rest of the information as shown below before pressing “Next >” again.
   

   Figure 4.41 Picking the "Tool button with dock widget" template in Plugin Builder
3. We will keep the default settings for the following setup pages, so just keep pressing “Next >” until you arrive at the last page where the button is called “Generate” instead. Take a look at the paths where Plugin Builder is going to create the folder for the new plugin. It might be a good idea to copy the path to somewhere from which you can retrieve it at any time. Now click on “Generate” and the folder for the plugin will be created.
4. If you now navigate to the plugins folder from the previous step, it will contain a new folder called “bus_track_analyzer_for_qgis”. Enter that folder and check out the files that have been created. The content should look like this:
   

   Figure 4.42 Newly created plugin folder with files for our plugin
5. Please close QGIS for a moment. We still need to compile the file resources.qrc into a Python .py file with the help of pyrcc5 as we did in Section 4.11. For this, run OSGeo4W.bat from your OSGeo4W/QGIS installation and, in the shell that opens up, navigate to the folder containing your plugin, e.g. type in the following command but adapt the path to match the path from steps 3 and 4:

   cd C:\Users\xyz\AppData\roaming\QGIS\QGIS3\profiles\default\python\plugins\bus_track_analyzer_for_qgis

6. Finally, compile the file with the command below. There should now be a file called resources.py in your plugin folder.

   pyrcc5 resources.qrc –o resources.py

7. You can now restart QGIS and open the plugin manager where “Bus Track Analyzer for QGIS” should now appear in the list of installed plugins. Enable it and then check the Plugins menu to make sure there is now an entry “Bus Track Analyzer for QGIS” there. Start the plugin and a rather empty dock widget will appear in the right part of the QGIS window. Try moving this dock widget around. You will see that you can move it to the other areas of the main QGIS window or completely undock it and have it as an independent window on the screen.

For the next steps, it’s best if you again close QGIS for a bit. In case you made any changes to the files during the bus tracking project in Section 4.10, it would be best if you re-download them from here. Then copy the following files from the Section 4.10 project folder (if you didn't edit anything) or the fresh download into the folder for the Bus Tracker plugin:

bus_events.py 
bus_track_analyzer.py 
bus_tracker_widget.py 
core_classes.py 
dublin_bus_data.csv
dublin_depots.csv
status_driving.png 
status_indepot.png 
status_stopped.png 

Please note that we are intentionally not including main.py. Also, it wouldn’t really be necessary to include the input data sets (csv files), but there is also no harm in doing so, and it means that we have everything needed to create and run the plugin together in the same folder.

     

Figure 4.43 Plugin folder with files from the Bus Track Analyzer project copied in 

The GUI we will be using for our dock widget is shown in the image below. It has an area at the top where the user can select the GPS and depot input files and a button “Read and init” for reading in the data from the selected files. The central area contains a QScrollArea widget that will host our BusTrackerWidget in the same way as we had it embedded into the main window in the original project. In the area at the bottom, we have the controls for running the analysis consisting of three buttons “Stop and reset”, “Pause”, and “Start” and a QSlider widget for setting the delay between consecutive analysis steps. The image also shows the object names of the important GUI elements that will become instance variables of class BusTrackAnalyzerForQGISDockWidget that we can access and connect to.

     

Figure 4.44 GUI for the dock widget of our plugin with object names of the different GUI elements as labels

If you look at the files in the folder for our plugin, you will see that Plugin Builder has created a file called bus_track_analyzer_for_qgis_dockwidget.py. This file contains the definition of class BusTrackAnalyzerForQGISDockWidget derived from QDockWidget with the GUI for our plugin. The class itself directly reads the GUI specification from the file bus_track_analyzer_for_qgis_dockwidget_base.ui as explained in Section 4.11.

So the next thing we are going to do is open that .ui file in QT Designer and modify it so that we get the GUI shown in the previous image. The image below shows the new GUI and its widget hierarchy in QT Designer. You don’t have to create this yourself. It is okay if you download the resulting .ui file and extract it into the plugin folder overwriting the default file that is already there (you might need to rename the downloaded file to match the default file). Then open the .ui file in QT Designer for a moment and have a closer look at how the different widgets have been arranged.

     

Figure 4.45 GUI for our dock widget and widget hierarchy in QT Designer

The next thing we are going to do is make a few smaller changes to the files we copied over from the original project. First of all, it is unfortunately required that we adapt all import statements in which we are importing .py files located in the project folder. The reason is that when we write something like

from core_classes import BusTracker 

, this will work fine when the file core_classes.py is located in the current working directory when the program is executed. This is usually the same folder in which the main Python script is located. Therefore, we didn’t have any problems when executing main.py since main.py and the other .py files we wrote are all in the same directory. However, when being run as a QGIS plugin, the working directory will not be the folder containing the plugin code. As a result, you will get error messages when trying to import the other .py files like this. What we have to do is adapt the import statements to start with a dot which tells Python to look for the file in the same folder in which the file in which the import statement appears is located. So the previous example needs to become:

from .core_classes import BusTracker 

Here is quick overview of where we have to make these changes:

1. In bus_events.py change the import statement at the beginning to:

   from .core_classes import Depot, BusTracker 

2. In bus_track_analyzer.py change the two import statement at the beginning to:

   from .bus_events import BusEvent 
   from .core_classes import BusTracker, Observation, Depot

3. Lastly, in bus_tracker_widget.py change the import statement at the beginning to:

   from .core_classes import BusTracker 

In addition to adapting the import statement, we are going to slightly adapt the bus_track_analyzer.py code to better work in concert with the GUI related classes of our plugin code: we are going to add the functionality to emit signals that we can connect to using the QT signal-slot approach. The two signals we are going to add are the following:

• observationProcessed: This signal will be emitted at the end of the nextStep() method, so whenever the BusTrackAnalyzer object is done with processing one observation, the observation object that has just been processed will be included with the signal.
• eventDetected: This signal will be emitted for each new event that is added to the allEvents list. The emitted signal includes the new bus event object.

Both signals will be used to connect an object of a new class we are going to write in Section 4.12.6 that has the purpose of showing the developing bus tracks and detected events live in QGIS. For this, it is required that the object be informed about newly processed observations and newly detected events, and this is what we are going to facilitate with these signals. Luckily, adding these signals to bus_track_analyzer.py just requires you to make a few small changes:

4. Add the following PyQt5 import statement somewhere before the line that starts the class definition with “class BusTrackAnalyzer …”:

   from PyQt5.QtCore import QObject, pyqtSignal

5. To be able to emit signals we need the class BusTrackerAnalyzer to be derived from the QT class QObject we are importing with the newly added import statement. Therefore, please change the first line of the class definition to:

   class BusTrackAnalyzer(QObject):

6. Then directly before the start of the constructor with “def __init__...", add the following two lines indented relative to the “class BusTrackAnalyzer…” but not part of any method:

   observationProcessed = pyqtSignal(Observation) 
   eventDetected = pyqtSignal(BusEvent)

   With these two lines we are defining the two signals that can be emitted by this class and the types of the parameters they will include.

7. Since we are now deriving BusTrackAnalyzer from QObject, we should call the constructor of the base class QObject from the first line of our __init__(…) method. So please add the following line there:

   super(BusTrackAnalyzer, self).__init__()

8. The last two changes we need to make are both in the definition of the nextStep() method. We therefore show the entire new version of that method here:

     def nextStep(self): 
            """performs next step by processing Observation at the front of the Observations priority queue""" 
            observation = heapq.heappop(self._observationQueue)           # get Observation that is at front of queue 
   
            # go through list of BusEvent subclasses and invoke their detect() method; then collect the events produced 
            # and add them to the allEvents lists 
            for evClass in self._eventClasses: 
                eventsProduced = evClass.detect(observation, self._depotData, self.allBusTrackers) # invoke event detection method 
                self.allEvents.extend(eventsProduced)  # add resulting events to event list 
                for event in eventsProduced:            
                    self.eventDetected.emit(event) 
   
            # update BusTracker of Observation that was just processed 
            observation.busTracker.lastProcessedIndex += 1 
            observation.busTracker.updateSpeed() 
   
            if observation.busTracker.status == BusTracker.STATUS_STOPPED: # if  duration of a stopped event has just expired, change status to "DRIVING" 
                if observation.timepoint.time > observation.busTracker.statusEvent.timepoint.time + observation.busTracker.statusEvent.duration: 
   		        observation.busTracker.status = BusTracker.STATUS_DRIVING 
                   observation.busTracker.statusEvent = None 
   
            # if this was not the last GPS Timepoint of this bus, create new Observation for the next point and add it to the Observation queue 
            if observation.timepointIndex < len(observation.busTracker.bus.timepoints) - 1:  # not last point 
                heapq.heappush(self._observationQueue, Observation(observation.busTracker, observation.timepointIndex + 1)  )  
   
            # update analyzer status 
            self.lastProcessedTimepoint = observation.timepoint 
   
            self.observationProcessed.emit(observation)

   In lines 10 and 11 of this new version we added a for-loop that goes through the events produced by the previous call of detect(…) and emit an eventDetected signal for each using the bus event object as a parameter. In the last line of the method, we do the same with the observationProcessed signal including the just processed Observation object.

At this point, our plugin is still missing the code that ties everything together, that is, the code that reads in the data from the input files when the “Read and init” button is clicked and reacts to the control buttons at the bottom of the BusTrackAnalyzerForQGISDockWidget widget by starting to continuously call the analyzer’s nextStep() method, pausing that process, or completely resetting the analysis to start from the beginning. We are going to place the code for this directly in the definition of class BusTrackAnalyzerForQGISDockWidget in bus_track_analyzer_for_qgis_dockwidget.py, so you should open that file for editing now. Here is the code that needs to be added, together with some explanations.

1. First, we need to important quite a few classes from our own .py files and also a few additional PyQt5 classes; so please change the import statements at the beginning of the code to the following:

   import os 
   
   from PyQt5 import QtGui, QtWidgets, uic 
   from PyQt5.QtCore import pyqtSignal, QTimer, QCoreApplication 
   
   # added imports 
   from .bus_track_analyzer import BusTrackAnalyzer 
   from .bus_tracker_widget import BusTrackerWidget 
   #from .qgis_event_and_track_layer_creator import QGISEventAndTrackLayerCreator 
   from .core_classes import Bus, Depot 
   from .bus_events import LeavingDepotEvent, EnteringDepotEvent, BusStoppedEvent, BusEncounterEvent 
   

   Note that we again have to use the notation with the dot at the beginning of the module names here. Also, please note that there is one import statement that is still commented out because it is for a class that we have not yet written. That will happen a bit later in Section 4.12.6 and we will then uncomment this line.

2. Next, we are going to add some initialization code to the constructor, directly after the last line of the __init__(…) method saying “self.setupUi(self)”:

            # own code added to template file 
   
            self.running = False # True if currently running analysis 
            self.delay = 0       # delay between steps in milliseconds 
   
            self.eventClasses = [ LeavingDepotEvent, EnteringDepotEvent, BusStoppedEvent, BusEncounterEvent ] # list of event classes to detect 
            self.busFileColumnIndices = { 'lon': 8, 'lat': 9, 'busId': 5, 'time': 0, 'line': 1 } # dictionary of column indices for required info 
   
            # create initial BusTrackAnalyzer and BusTrackerWidget objects, and add the later to the scroll area of this widget 
            self.analyzer = BusTrackAnalyzer({}, [], self.eventClasses) 
            self.trackerWidget = BusTrackerWidget(self.analyzer) 
            self.busTrackerContainerWidget.layout().addWidget(self.trackerWidget,0,0) 
            self.layerCreator = None    # QGISEventAndTrackLayerCreator object, will only be initialized when input files are read 
   
            # create a QTimer user to wait some time between steps  
            self.timer = QTimer() 
            self.timer.setSingleShot(True) 
            self.timer.timeout.connect(self.step) 
   
            # set icons for play control buttons 
            self.startTB.setIcon(QCoreApplication.instance().style().standardIcon(QtWidgets.QStyle.SP_MediaPlay)) 
            self.pauseTB.setIcon(QCoreApplication.instance().style().standardIcon(QtWidgets.QStyle.SP_MediaPause)) 
            self.stopAndResetTB.setIcon(QCoreApplication.instance().style().standardIcon(QtWidgets.QStyle.SP_MediaSkipBackward)) 
   
            # connect play control buttons and slide to respetive methods defined below 
            self.startTB.clicked.connect(self.start) 
            self.pauseTB.clicked.connect(self.stop) 
            self.stopAndResetTB.clicked.connect(self.reset) 
            self.delaySlider.valueChanged.connect(self.setDelay) 
   
            # connect edit fields and buttons for selecting input files to respetive methods defined below 
            self.browseTrackFileTB.clicked.connect(self.selectTrackFile) 
            self.browseDepotFileTB.clicked.connect(self.selectDepotFile) 
            self.readAndInitPB.clicked.connect(self.readData)

   What happens in this piece of code is the following:

   • First, we introduce some new instance variables for this class that are needed for reading data and setting up a BusTrackAnalyzer object (you probably recognize the variables eventClasses and busFileColumnIndices from the original main.py file) and for controlling the analysis steps. Variable running will be set to True when the analyzer’s nextStep() method is supposed to be called continuously with small delays between the steps whose length is given by variable delay. It is certainly not optimal that the column indices for the GPS data file are hard-coded here in variable busFileColumnIndices but we decided against including an option for the user to specify these in this version to keep the GUI and code as simple as possible.
   • In the next block of code, we create an instance variable to store the BusTrackAnalyzer object we will be using, and we create the BusTrackerWidget widget and place it within the scroll area in the center of the GUI. Since we have not read any input data yet, we are using an empty bus dictionary and empty depot list to set these things up.
   • An instance variable for a QTimer object is created that will be used to queue up the execution of the method step() defined later when the widget is in continuous execution mode.
   • The next code block simply sets the icons for the different control buttons at the bottom of the widget using standard QT icons.
   • Next, we connect the “clicked” signals of the different control buttons to different methods of the class that are supposed to react to these signals. We do the same for the “valueChanged” signal of the QSlider widget in our GUI to be able to adapt the value of the variable delay when the slider is moved by the user.
   • Finally, we link the three buttons at the top of the GUI for reading in the data to the respective event handler methods.

3. Now the last thing that needs to happen is adding the different event handler methods we have already been referring to in the previously added code. This is another larger chunk of code since there are quite a few methods to define. Please add the definitions at the end of the file after the definition of the method closeEvent(…) that is already there by default.

        # own methods added to template file 
   
        def selectTrackFile(self): 
            """displays open file dialog to select bus track input file""" 
            fileName, _ = QtWidgets.QFileDialog.getOpenFileName(self,"Select CSV file with bus track data", "","(*.*)") 
            if fileName: 
                self.trackFileNameLE.setText(fileName) 
   
        def selectDepotFile(self): 
            """displays open file dialog to select depot input file""" 
            fileName, _ = QtWidgets.QFileDialog.getOpenFileName(self,"Select CSV file with depot data", "","(*.*)") 
            if fileName: 
                self.depotFileNameLE.setText(fileName) 
   
        def readData(self): 
            """reads bus track and depot data from selected files and creates new analyzer and creates analyzer and layer creator for new input""" 
            if self.running: 
                self.stop() 
            try:   # read data 
                depotData = Depot.readFromCSV(self.depotFileNameLE.text())   
                busData = Bus.readFromCSV(self.trackFileNameLE.text(), self.busFileColumnIndices)  
            except Exception as e: 
                QtWidgets.QMessageBox.information(self, 'Operation failed', 'Could not read data from files provided: '+ str(e.__class__) + ': ' + str(e), QtWidgets.QMessageBox.Ok) 
                busData = {} 
                depotData = [] 
   
            # create new analyzer and layer creator objects and connect them 
            self.analyzer = BusTrackAnalyzer(busData, depotData, self.eventClasses) 
            self.trackerWidget.analyzer = self.analyzer 
   #         self.createLayerCreator() 
            self.trackerWidget.updateContent() 
   
        def stop(self): 
            """halts analysis but analysis can be continued from this point""" 
            self.timer.stop() 
            self.running = False 
   
        def reset(self): 
            """halts analysis and resets analyzer to start from the beginning""" 
            self.stop() 
            self.analyzer.reset() 
   #         self.createLayerCreator() 
            self.trackerWidget.updateContent() 
   
        def start(self): 
            """starts analysis if analysis isn't already running""" 
            if not self.running: 
                self.running = True 
                self.step() 
   
        def step(self): 
            """performs a single analysis step of the BusTrackAnalyzer but starts singleshot timer after each step to call itself again""" 
            if self.running: 
                if self.analyzer.isFinished(): 
                    self.stop() 
                else: 
                    self.analyzer.nextStep()               # perform next analysis step 
                    self.trackerWidget.updateContent()     # redraw tracker widget  
                    self.timer.start(max([5,self.delay])) # start timer to call this method again after delay 
   
        def setDelay(self): 
            """adapt delay when slider has been moved""" 
            self.delay = 10 * self.delaySlider.value() 
            if self.running:                           # if analysis is running, change to the new delay immediately 
                self.timer.stop() 
                self.timer.start(max([5,self.delay]))
   

The first two methods selectTrackFile() and selectDepotFile() are called when the “…” buttons at the top are clicked and will open file dialog boxes for picking the input files. The method readData() is invoked when the “Read and init” button is clicked. It stops all ongoing executions of the analyzer, attempts to read the data from the selected files, and then creates a new BusTrackAnalyzer object for this input data and connects it to the BusTrackerWidget in our GUI. The code of this function contains another two lines that are commented out and that we will uncomment later.

The other methods we define in this piece of code are the event handler functions for the control buttons at the bottom:

• stop() is called when the “Pause” button in the middle is clicked and simply stops the timer if its currently running and sets the variable running to False so that the analyzer’s nextStep() method won’t be called anymore.
• reset() is called when the control button on the left is called and it also stops the execution but, in addition, resets the analyzer so that the analysis will start from the beginning when the “Start” button is clicked again.
• start() is called when the right control button (“Play”) for continuously performing analysis steps is clicked. It sets variable running to True and then invokes method step() to run the first analysis step.
• step() is either called from start() or by the timer it sets up itself after each call of the analyzer’s nextStep() method to perform the next analysis step after a certain delay until all observations have been processed. Please note that we are using minimum delay of 5 milliseconds to make sure that the QGIS GUI remains responsive while we are running the analysis.
• setDelay() is called when the slider is moved and it translates the slider position into a delay value between 0 and 1 second. It also immediately restarts the timer to use this new delay value.

Our plugin is operational now and you can open QGIS and run it. In case you already have QGIS running or you encounter any errors that need to be fixed, don’t forget to reload the plugin code with the help of the Plugin Reloader plugin. Once the dock widget appears at the right side of the QGIS window (as shown in the figure below), do the following:

1. Use the “…” buttons to select the input files for the bus GPS and depot data.
2. Click “Read and init”; after that you should see the initial bus tracker configuration in the central area of the dock widget.
3. Press the “Play” button to start the analysis; the content of the BusTrackerWidget should now continuously update to reflect the current state of the analysis
4. Test out the “Pause” and “Rewind” buttons as well as changing the delay between steps with the slider to control the speed the analysis is run at.

Figure 4.46 The plugin running in QGIS showing the current status of the analysis and individual buses

So far, so good. We have now turned our original standalone project into a QGIS plugin and even added in some extra functionality allowing the user to pause and restart the analysis and control the speed. However, typically a QGIS plugin in some way interacts with the content of the project that is currently open in QGIS, for instance by taking some of its layers as input or adding new layers to the project. We will add this kind of functionality in the next section and this will be the final addition we are making to our plugin.

In addition to showing the current bus tracker states in the dock widget, we want our plugin to add two new layers to the currently open QGIS project that show the progress of the analysis and, once the analysis is finished, contain its results. The two layers will correspond to the two output files we produced in the original project in Section 4.10:

• a polyline layer showing routes taken by the different bus vehicles with the bus IDs as attribute,
• a point layer showing the detected events with the type of the event and a short description as attributes.

Since we don’t want to just produce these layers at the end of the analysis but want these to be there from the start of the analysis and continuously update whenever a new bus GPS observation is processed or an event is detected, we are going to write some code that reacts to the observationProcessed and eventDetected signals emitted by our class BusTrackAnalyzer (see the part of Section 4.12.4 where we added these signals). We will define a new class for all this that will be called QGISEventAndTrackLayerCreator and it will be defined in a new file qgis_event_and_track_layer_creator.py. The class definition consists of the constructor and two methods called addObservationToTrack(…) and addEvent(…) that will be connected to the corresponding signals of the analyzer.

Let’s start with the beginning of the class definition and the constructor. All following code needs to be placed in file qgis_event_and_track_layer_creator.py that you need to create in the plugin folder.

import qgis 

class QGISEventAndTrackLayerCreator(): 

     def __init__(self): 
         self._features = {}    # dictionary mapping bus id string to polyline feature in bus track layer  
         self._pointLists = {}  # dictionary mapping bus id string to list of QgsPointXY objects for creating the poylines from 

         # get project currently open in QGIS 
         currentProject = qgis.core.QgsProject.instance() 

         # create track layer and symbology, then add to current project   
         self.trackLayer = qgis.core.QgsVectorLayer('LineString?crs=EPSG:4326&field=BUS_ID:integer', 'Bus tracks' , 'memory') 
         self.trackProv = self.trackLayer.dataProvider() 

         lineMeta = qgis.core.QgsApplication.symbolLayerRegistry().symbolLayerMetadata("SimpleLine") 
         lineLayer = lineMeta.createSymbolLayer({'color': '0,0,0'}) 
         markerMeta = qgis.core.QgsApplication.symbolLayerRegistry().symbolLayerMetadata("MarkerLine") 

         markerLayer = markerMeta.createSymbolLayer({'width': '0.26', 'color': '0,0,0', 'placement': 'lastvertex'}) 
         symbol = qgis.core.QgsSymbol.defaultSymbol(self.trackLayer.geometryType()) 
         symbol.deleteSymbolLayer(0) 
         symbol.appendSymbolLayer(lineLayer) 
         symbol.appendSymbolLayer(markerLayer) 

         trackRenderer = qgis.core.QgsSingleSymbolRenderer(symbol) 
         self.trackLayer.setRenderer(trackRenderer) 

         currentProject.addMapLayer(self.trackLayer) 

         # create event layer and symbology, then add to current project 
         self.eventLayer = qgis.core.QgsVectorLayer('Point?crs=EPSG:4326&field=TYPE:string(50)&field=INFO:string(255)', 'Bus events' , 'memory') 
         self.eventProv = self.eventLayer.dataProvider() 

         colors = { "BusEncounterEvent": 'yellow', "BusStoppedEvent": 'orange', "EnteringDepotEvent": 'blue', "LeavingDepotEvent": 'green' } 

         categories = [] 
         for ev in colors: 
             categories.append( qgis.core.QgsRendererCategory( ev, qgis.core.QgsMarkerSymbol.createSimple({'name': 'square', 'size': '3.0', 'color': colors[ev]}), ev )) 

         eventRenderer = qgis.core.QgsCategorizedSymbolRenderer("TYPE", categories) 
         self.eventLayer.setRenderer(eventRenderer) 

         currentProject.addMapLayer(self.eventLayer) 

To be able to build polylines for the bus tracks and update these whenever a new observation has been processed by the analyzer, we need to maintain dictionaries with the QGIS features and point lists for each bus vehicle track. These are created at the beginning of the constructor code in lines 6 and 7. In addition, the constructor accesses the currently open QGIS project (line 10) and adds the two new layers called “Bus tracks” and “Bus events” to it (lines 29 and 44). The rest of the code is mainly for setting the symbology of these two layers: For the track layer, we use black lines with a red circle marker at the end to indicate the current location of the vehicle (lines 16 to 27) as shown in the image below. For the events, we use square markers in different colors based on the TYPE of the event (lines 35 to 42).

     

Figure 4.47 QGIS symbology used for the bus track and event layers

Now we are going to add the definition of the method addObservationToTrack(…) that will be connected to the observationProcessed signal emitted when the analyzer object has completed the execution of nextStep().

	def addObservationToTrack(self, observation): 
         """add new vertex to a bus polyline based on the given Observation object""" 
         busId = observation.busTracker.bus.busId; 

         # create new point for this observation 
         p =  qgis.core.QgsPointXY(observation.timepoint.lon, observation.timepoint.lat) 

         # add point to point list and (re)create polyline geometry 
         if busId in self._features:      # we already have a point list and polyline feature for this bus 
             feat = self._features[busId] 
             points = self._pointLists[busId] 
             points.append(p) 

             # recreate polyline geometry and replace in layer 
             polyline = qgis.core.QgsGeometry.fromPolylineXY(points) 
             self.trackProv.changeGeometryValues({feat.id(): polyline}) 

         else:                            # new bus id we haven't seen before 
             # create new  polyline and feature     
             polyline = qgis.core.QgsGeometry.fromPolylineXY([p]) 
             feat = qgis.core.QgsFeature() 
             feat.setGeometry(polyline) 
             feat.setAttributes([int(busId)]) 
             _, f = self.trackProv.addFeatures([feat]) 
             # store point list and polyline feature in respective dictionaries 
             self._features[busId] = f[0] 
             self._pointLists[busId] = [p] 

         # force redraw of layer 
         self.trackLayer.triggerRepaint() 
         qgis.utils.iface.mapCanvas().refresh()

The Observation object given to this method as a parameter provides us with access to all the relevant information we need to update the polyline feature for the bus this observation is about. First, we extract the ID of the bus (line 3) and create a new QgsPointXY object from the Timepoint stored in the Observation object (line 6). If we already have a polyline feature for this vehicle, we get the corresponding feature and point lists from the features and pointLists dictionaries, add the new point to the point list and create a new polyline geometry from it, and finally change the geometry of that feature in the bus track layer to this new geometry (lines 10 to 16). If instead this is the first observation of this vehicle, we create a point list for it to be stored in the pointList dictionary as well as a new polyline geometry with just that single point, and we then set up a new QgsFeature object for this polyline that is added to the bus track layer and also the features dictionary (lines 19 to 27). At the very end of the method, we make sure that the layer is repainted in the QGIS map canvas.

Now we add the code for the addEvent(…) method completing the definition of our class QGISEventAndTrackLayerCreator:

		def addEvent(self, busEvent): 
         	"""add new event point feature to event layer based on the given BusEvent object""" 
            # create point feature with information from busEvent 
            p =  qgis.core.QgsPointXY(busEvent.timepoint.lon, busEvent.timepoint.lat) 
            feat = qgis.core.QgsFeature() 
            feat.setGeometry(qgis.core.QgsGeometry.fromPointXY(p)) 
            feat.setAttributes([type(busEvent).__name__, busEvent.description()]) 

            # add feature to event layer and force redraw 
            self.eventProv.addFeatures([feat]) 

            self.eventLayer.triggerRepaint() 
            qgis.utils.iface.mapCanvas().refresh()

This method is much simpler because we don’t have to modify existing features in the layer but rather always add one new point feature to the event layer. All information required for this is taken from the bus event object given as a parameter: The coordinates for the next point feature are taken from the Timepoint stored in the event object (line 4), for the TYPE field of the event we take the type of the event object (line 7), and for the INFO field we take the string returned by calling the event object’s description(…) method (also line 7).

To incorporate this new class into our current plugin, we need to make a few more modifications to the class BusTrackAnalyzerForQGISDockWidget in file bus_track_analyzer_for_qgis_dock_widget.py. Here are the instructions for this:

1. Remove the # at the beginning of the line “from .qgis_event_and_track_layer_creator import …” to uncomment this line and import our new class.
2. Now we add a new auxiliary method to the class definition that has the purpose of creating an object of our new class BusTrackAnalyzerForQGISDockWidget and connecting its methods to the corresponding signals of the analyzer object. Add the following definition as the first method after the “closeEvent” method, so directly after the comment “# own methods added to the template file”.

   	def createLayerCreator(self): 
            """creates a new QGISEventAndTrackLayerCreator for showing events and tracks on main QGIS window and connects it to analyzer""" 
            self.layerCreator = QGISEventAndTrackLayerCreator() 
            self.analyzer.observationProcessed.connect(self.layerCreator.addObservationToTrack) 
            self.analyzer.eventDetected.connect(self.layerCreator.addEvent) 
   

   We already set up an instance variable layerCreator in the constructor code and we are using it here for storing the newly created layer creator object. Then we connect the signals to the two methods of the layer creator object.

3. Now we just need to uncomment two lines to make sure that the method createLayerCreator() is called at the right places, namely after we read in data from the input files and when the analyzer is reset. In both cases, we want to create a new layer creator object that will set up two new layers in the current QGIS project. The lines that you need to uncomment for this by removing the leading # are:
   • the line “# self.createLayerCreator()” in method readData()
   • the line “# self.createLayerCreator()” in method reset()

   You should make sure that the indentation is correct after removing the hashmarks.

That’s it, we are done with the code for our plugin!

To try out this new version of the plugin, close the dock widget if it’s currently still open in QGIS and then run the Plugin Reloader plugin to load this updated version of our plugin. Add a basemap to your map project and zoom it to the general area of Dublin. When you then load the input files, you will see the two new layers appear in the layer list of the current QGIS project with the symbology we are setting up in the code. When you now start the analysis, what you see should look like the video below with the bus tracker widget continuously updating the bus status information, the bus tracks starting to appear in the QGIS map window, and square symbols starting to pop up for the events detected.

[NOTE:  This video (3:04) does NOT contain sound]

The delay slider can be used to increase the breaks between two analysis steps, which can be helpful if the map window doesn’t seem to update properly because QGIS has problems catching up with the requests to repaint layers. This in particular is a good idea if you want to pan and zoom the map in which case you may notice the basemap tiles not appearing if the delay is too short.

Overall, there is quite a bit that could be optimized to make sure that QGIS remains responsive while the plugin and analysis are running as well as other improvements and extensions that could be made. But all this would increase the amount of code needed quite a bit and this has already been a rather long project to begin with, requiring you to read and understand a lot of code. As we said before, it is not required to understand each line in the code; the crucial points to understand are how we are using classes, objects, and inheritance in this project and make use of the other techniques and concepts taught in this lesson. Reading and understanding other people’s code is one of the main ways to become a better programmer and since we are approaching the end of this course, this was a good place to practice this a bit and maybe provide some inspiration for your term project. However, we certainly don’t expect your term project to be nearly as complex as the plugin created in this section!

Solution

import qgis 
import qgis.core
import sys, os 
import pandas as pd 


# read data into pandas data frame 
data = pd.read_csv(r"C:\489\L4\exercise\L4exercise_data.csv") 

# class definition for PointObject 
class PointObject(): 

    def __init__(self, lat, lon):  
        self.lon = lon 
        self.lat = lat 

    def toQgsFeature(self):  
         feat = qgis.core.QgsFeature() 
         feat.setGeometry(qgis.core.QgsGeometry.fromPointXY(qgis.core.QgsPointXY(self.lon, self.lat))) 
         return feat 

firstObservations = []          # for storing objects of class pointWithID 
firstObservationsFeatures = []  # for storing objects of class QgsFeature 

# code for creating QgsApplication and initializing QGIS environment 
qgis_prefix = os.getenv("QGIS_PREFIX_PATH")                
qgis.core.QgsApplication.setPrefixPath(qgis_prefix, True) 
qgs = qgis.core.QgsApplication([], False)
qgs.initQgis() 
 
# definition of class PointWithID derived from PointObject  
# to represent animal observation from the input data 
class PointWithID(PointObject): 
 
    def __init__(self, pID, lat, lon): 
        super(PointWithID, self).__init__(lat, lon) 
        self.pID = pID  # instance variable for storing animal ID 

    # overwriting the == operator to be based on the animal ID 
    def __eq__(self, other): 
        return self.pID == other.pID 

    # overwriting this method to include animal ID as attribute of QgsFeature created 
    def toQgsFeature(self): 
         feat = super(PointWithID, self).toQgsFeature() 
         feat.setAttributes([self.pID]) 
         return feat 

# create list of PointWithID object with first observations for each animal in the data frame 
for row in data.itertuples(index=False): 
    pointWithID = PointWithID(row[0], row[1], row[2])  
    if not pointWithID in firstObservations: # here __eq__() is used to do the comparison 
        firstObservations.append(pointWithID) 

# list comprehension for creating list of features from firstObservations list 
firstObservationsFeatures = [ o.toQgsFeature() for o in firstObservations ] 

# create new point layer with field for animal ID 
layer = qgis.core.QgsVectorLayer("Point?crs=EPSG:4326&field=AnimalID:string(255)", 'animal first observations' ,'memory') 

# add features to new layer 
prov = layer.dataProvider() 
prov.addFeatures(firstObservationsFeatures) 

# save layer as GeoPackage file 
qgis.core.QgsVectorFileWriter.writeAsVectorFormat( layer, r"C:\489\L4\exercise\firstobservations.gpkg", "utf-8",layer.crs(), "GPKG") 

# clean up 
qgs.exitQgis() 

Explanation

• In line 34 the new class PointWithID is declared to be derived from class PointObject.
• The constructor in line 36 calls the constructor of the base class and then adds the additional instance variable for the animal ID (called pID).
• __eq__() is defined in lines 43 and 44 to compare the pID instance variables of both involved objects and return True if they are equal.
• In lines 48 and 49, we call toQgsFeature() of the base class PointObject and then take the returned feature and add the animal ID as the only attribute.
• In lines 53 to 56, we use itertuples() of the data frame to loop through its rows. We create a new PointWithID object in variable pointWithID from the cell values of the row and then test whether an object with that ID is already contained in list firstObservations. If not, we add this object to the list.
• In line 59, we use toQgsFeature() to create a QgsFeature object from each PointWithID object in the list. These are then added to the new layer in line 66.
• When creating the new layer in line 62, we have to make sure to include the field for storing the animal ID in the string for the first parameter.