Object-oriented programming is a software development paradigm in which concepts are modeled as objects that have properties and behaviors represented as attributes and methods, respectively. The goal of this paradigm is more modular software in which one object can inherit from one or more other objects to encourage software reuse.

The Python programming language is known for its ability to serve multiple roles as a well-designed, object-oriented language, a procedural scripting language, or even a functional programming language. However, you never completely abandon object-oriented programming in Python because even its native data types are objects and all Python libraries, known as modules, adhere to a basic object structure and behavior.

Geospatial analysis is the perfect activity for object-oriented programming. In most object-oriented programming projects, the objects are abstract concepts such as database connections that have no real-world analogy. However, in geospatial analysis, the concepts modeled are, well, real-world objects! The domain of geospatial analysis is the Earth and everything on it. Trees, buildings, rivers, and people are all examples of objects within a geospatial system.

A common example in literature for newcomers to object-oriented programming is the concrete analogy of a cat. Books on object-oriented programming frequently use some form of the following example.

Imagine that you are looking at a cat. We know some information about the cat, such as its name, age, color, and size. These features are the properties of the cat. The cat also exhibits behaviors such as eating, sleeping, jumping, and purring. In object-oriented programming, objects have properties and behaviors too. You can model a real-world object such as the cat in our example or something more abstract such as a bank account.

Most concepts in object-oriented programming are far more abstract than the simple cat paradigm or even the bank account. However, in geospatial analysis, the objects that are modeled remain concrete, such as the simple cat analogy, and in many cases are living, breathing cats. Geospatial analysis allows you to continue with the simple cat analogy and even visualize it. The following map represents the feral cat population of Australia using data provided by the Atlas of Living Australia (ALA):

As its name suggests, a thematic map portrays a specific theme. A general reference map visually represents features as they relate geographically for navigation or planning. A thematic map goes beyond location to provide the geographic context for information around a central idea. Usually, a thematic map is designed for a targeted audience to answer specific questions. The value of thematic maps lies in what they do not show. A thematic map will use minimal geographic features to avoid distracting the reader from the theme. Most thematic maps include political boundaries such as country or state borders but omit navigational features, such as street names or points of interest beyond major landmarks that orient the reader. The cholera map by Dr. John Snow earlier in this chapter is a perfect example of a thematic map. Common uses for thematic maps are visualizing health issues, such as disease, election results, and environmental phenomena such as rainfall. These maps are also the most common output of geospatial analysis. The following map from the United States Census Bureau shows cancer mortality rates by state:

Thematic maps tell a story and are very useful. However, it is important to remember that while thematic maps are models of reality as any other map, they are also generalizations of information. Two different analysts using the same source of information will often come up with very different thematic maps depending on how they analyze and summarize the data. They may also choose to focus on different aspects of the dataset. The technical nature of thematic maps often leads people to treat them as if they are scientific evidence. However, geospatial analysis is often inconclusive. While the analysis may be based on scientific data, the analyst does not always follow the rigor of the scientific method. In his classic book, How to Lie with Maps, Mark Monmonier, University of Chicago Press, demonstrates in great detail how maps are easily manipulated models of reality, which are commonly abused. This fact doesn't degrade the value of these tools. The legendary statistician, George Box, wrote in his 1987 book, Empirical Model-Building and Response Surfaces that, "Essentially, all models are wrong, but some are useful." Thematic maps have been used as guides to start (and end) wars, stop deadly diseases in their tracks, win elections, feed nations, fight poverty, protect endangered species, and rescue those impacted by disaster. Thematic maps may be the most useful models ever created.

In its purest form, a database is simply an organized collection of information. A database management system (DBMS) is an interactive suite of software that can interact with a database. People often use the word database as a catch-all term referring to both the DBMS and underlying data structure. Databases typically contain alphanumeric data and, in some cases, binary large objects or blobs, which can store binary data such as images. Most databases also allow a relational database structure in which entries in normalized tables can be referenced to each other in order to create many-to-one and one-to-many relationships among data.

Spatial databases, also known as geodatabases, use specialized software to extend a traditional relational database management system (RDBMS) to store and query data defined in two-dimensional or three-dimensional space. Some systems also account for a series of data over time. In a spatial database, attributes about geographic features are stored and queried as traditional relational database structures. The spatial extensions allow you to query geometries using Structured Query Language (SQL) in a similar way as traditional database queries. Spatial queries and attribute queries can also be combined to select results based on both location and attributes.

Spatial indexing is a process that organizes the geospatial vector data for faster retrieval. It is a way of prefiltering the data for common queries or rendering. Indexing is commonly used in large databases to speed up returns to queries. Spatial data is no different. Even a moderately-sized geodatabase can contain millions of points or objects. If you perform a spatial query, every point in the database must be considered by the system in order to include it or eliminate it in the results. Spatial indexing groups data in ways that allow large portions of the dataset to be eliminated from consideration by doing computationally simpler checks before going into a detailed and slower analysis of the remaining items.

Metadata is defined as data about data. Accordingly, geospatial metadata is data about geospatial datasets that provide traceability for the source and history of a dataset as well as summary of the technical details. Metadata also provides long-term preservation of information holdings. Geospatial metadata can be represented by several possible standards. One of the most prominent standards is the international standard, ISO 19115-1, which includes hundreds of potential fields to describe a single geospatial dataset. Additionally, the ISO 19115-2 includes extensions for geospatial imagery and gridded data. Some example fields include spatial representation, temporal extent, and lineage. The primary use of metadata is cataloging datasets. Modern metadata can be ingested by geographic search engines making it potentially discoverable by other systems automatically. It also lists points of contact for a dataset if you have questions. Metadata is an important support tool for geospatial analysts and adds credibility and accessibility to your work. The Open Geospatial Consortium (OGC) created the Catalog Service for the Web (CSW) to manage metadata. The pycsw Python library implements the CSW standard. You can learn more about it at http://pycsw.org.

Map projections have entire books devoted to them and can be a challenge for new analysts. If you take any three-dimensional object and flatten it on a plane, such as your screen or a sheet of paper, the object is distorted. Many grade school geography classes demonstrate this concept by having students peel an orange and then attempt to lay the peel flat on their desk in order to understand the resulting distortion. The same effect occurs when you take the round shape of the Earth and project it on a computer screen.

In geospatial analysis, you can manipulate this distortion to preserve common properties, such as area, scale, bearing, distance, or shape. There is no one-size-fits-all solution to map projections. The choice of projection is always a compromise of gaining accuracy in one dimension in exchange for error in another. Projections are typically represented as a set of over 40 parameters as either XML or a text format called Well-Known Text (WKT), which is used to define the transformation algorithm.

The International Association of Oil & Gas Producers (IOGP) maintains a registry of the most known projections. The organization was formerly known as the European Petroleum Survey Group (EPSG). The entries in the registry are still known as EPSG codes. The EPSG maintained the registry as a common benefit for the oil and gas industry, which is a prolific user of geospatial analysis for energy exploration. At the last count, this registry contained over 5,000 entries.

As recently as 10 years ago, map projections were a primary concern for a geospatial analyst. Data storage was expensive, high-speed Internet was rare, and cloud computing didn't really exist. Geospatial data was typically exchanged among small groups working in separate areas of interest. The technology constraints at the time meant that geospatial analysis was highly localized. Analysts would use the best projection for their area of interest. Data in different projections cannot be displayed on the same map because they represent two different models of the Earth. Any time an analyst received data from a third party, it had to be reprojected before using it with the existing data. This process was tedious and time-consuming. Most geospatial data formats do not provide a way to store the projection information. This information is stored in an ancillary file as text or XML usually. As analysts didn't exchange data often, many people wouldn't bother defining projection information. Every analyst's nightmare was to come across an extremely valuable dataset missing the projection information. It rendered the dataset useless. The coordinates in the file are just numbers and offer no clue to the projection. With over 5,000 choices, it was nearly impossible to guess.

Now, thanks to modern software and the Internet making data exchange easier and more common, nearly every data format has added a metadata format that defines the projection or places it in the file header if supported. Advances in technology have also allowed for global basemaps, which allow for more common uses of projections such as the common Google Mercator projection used for Google Maps. This projection is also known as Web Mercator and uses code EPSG:3857 (or the deprecated EPSG:900913). Geospatial portal projects such as OpenStreetMap.org and NationalAtlas.gov have consolidated datasets for much of the world in common projections. Modern geospatial software can also reproject data on the fly saving the analyst the trouble of preprocessing the data before using it. Closely related to map projections are geodetic datums. A datum is a model of the Earth's surface used to match the location of features on the Earth to a coordinate system. One common datum is called WGS 84 that is used by GPS devices.

The exciting part of geospatial analysis is visualization. As geospatial analysis is a computer-based process, it is good to be aware of how geographic data appears on a computer screen.

Geographic data including points, lines, and polygons are stored numerically as one or more points, which come in (x,y) pairs or (x,y,z) tuples. The x represents the horizontal axis on a graph. The y represents the vertical axis. The z represents terrain elevation. In computer graphics, a computer screen is represented by an x and y axis. A z axis is not used because the computer screen is treated as a two-dimensional plane by most graphics software APIs. However, as desktop computing power continues to improve, three-dimensional maps are starting to become more common.

Another important factor is screen coordinates versus world coordinates. Geographic data is stored in a coordinate system representing a grid overlaid on the Earth, which is three-dimensional and round. Screen coordinates, also known as pixel coordinates, represent a grid of pixels on a flat, two-dimensional computer screen. Mapping x and y world coordinates to pixel coordinates is fairly straightforward and involves a simple scaling algorithm. However, if a z coordinate exists, then a more complicated transform must be performed to map coordinates from three-dimensional space to a two-dimensional plane. These transformations can be computationally costly and therefore slow if not handled correctly.

In the case of remote sensing data, the challenge is typically the file size. Even a moderately-sized satellite image that is compressed can be tens, if not hundreds, of megabytes. Images can be compressed using lossless or lossy methods. Lossless methods use tricks to reduce the file size without discarding any data. Lossy compression algorithms reduce the file size by reducing the amount of data in the image while avoiding a significant change in the appearance of the image. Rendering an image on the screen can be computationally-intensive. Most remote sensing file formats allow for the storing of multiple lower-resolution versions of the image—called overviews or pyramids—for the sole purpose of faster rendering at different scales. When zoomed out from the image to a scale where you couldn't see the detail of the full resolution image, a preprocessed, lower-resolution version of the image is displayed quickly and seamlessly.

Raster data is captured digitally as square tiles. This means that the data is stored on a computer as a numerical array of rows and columns. If the data is multispectral, the dataset will usually contain multiple arrays of the same size, which are geospatially referenced together to represent a single area on the Earth. These different arrays are called bands. Any numerical array can be represented on a computer as an image. In fact, all computer data is ultimately numbers. It is important in geospatial analysis to think of images as a numeric array because mathematical formulas are used to process them.

In remotely sensed images, each pixel represents both space (location on the Earth of a certain size) and the reflectance captured as light reflected from the Earth at this location into space. So, each pixel has a ground size and contains a number representing the intensity. As each pixel is a number, we can perform mathematical equations on this data to combine data from different bands and highlight specific classes of objects in the image. If the wavelength value is beyond the visible spectrum, we can highlight features not visible to the human eye. Substances such as chlorophyll in plants can be greatly contrasted using a specific formula called Normalized Difference Vegetation Index (NDVI).

By processing remotely sensed images, we can turn this data into visual information. Using the NDVI formula, we can answer the question, what is the relative health of the plants in this image? You can also create new types of digital information, which can be used as input for computer programs to output other types of information.

Computer screens display images as combinations of Red, Green, and Blue (RGB) to match the capability of the human eye. Satellites and other remote sensing imaging devices can capture light beyond this visible spectrum. On a computer, wavelengths beyond the visible spectrum are represented in the visible spectrum so that we can see them. These images are known as false color images. In remote sensing, for instance, infrared light makes moisture highly visible. This phenomenon has a variety of uses such as monitoring ground saturation during a flood or finding hidden leaks in a roof or levee.

GIS vector data uses coordinates consisting of, at a minimum, an x horizontal value and a y vertical value to represent a location on the Earth. In many cases, a point may also contain a z value. Other ancillary values are possible including measurements or timestamps.

These coordinates are used to form points, lines, and polygons to model real-world objects. Points can be a geometric feature in and of themselves or they can connect line segments. Closed areas created by line segments are considered polygons. Polygons model objects such as buildings, terrain, or political boundaries.

A GIS feature can consist of a single point, line, or polygon or it can consist of more than one shape. For example, in a GIS polygon dataset containing world country boundaries, the Philippines, which is made up of 7,107 islands, would be represented as a single country made up of thousands of polygons.

Vector data typically represents topographic features better than raster data. Vector data has better accuracy potential and is more precise. However, to collect vector data on a large scale is also traditionally more costly than raster data.

Two other important terms related to vector data structures are bounding box and convex hull. The bounding box or minimum bounding box is the smallest possible square that contains all of the points in a dataset. The following image demonstrates a bounding box for a collection of points:

The convex hull of a dataset is similar to the bounding box, but instead of a square, it is the smallest possible polygon that can contain a dataset. The bounding box of a dataset always contains its convex hull. The following image shows the same point data as the previous example with the convex hull polygon shown in red:

A buffer operation can be applied to spatial objects including points, lines, or polygons. This operation creates a polygon around the object at a specified distance. Buffer operations are used for proximity analysis, for example, establishing a safety zone around a dangerous area. In the following image, the black shapes represent the original geometry while the red outlines represent the larger buffer polygons generated from the original shape:

A dissolve operation creates a single polygon out of adjacent polygons. A common use for a dissolve operation is to merge two adjacent properties in a tax database that have been purchased by a single owner. Dissolves are also used to simplify data extracted from remote sensing:

Objects that have more points than necessary for the geospatial model can be generalized to reduce the number of points used to represent the shape. This operation usually requires a few attempts to get the optimal number of points without compromising the overall shape. It is a data optimization technique to simplify data for the efficiency of computing or better visualization. This technique is useful in web-mapping applications. Computer screens have a resolution of 72 dots per inch (dpi). Highly-detailed point data, which would not be visible, can be reduced so that less bandwidth is used to send a visually equivalent map to the user:

An intersection operation is used to see if one part of a feature intersects with one or more features. This operation is for spatial queries in proximity analysis and is often a follow-on operation to a buffer analysis:

A merge operation combines two or more non-overlapping shapes in a single multishape object. Multishape objects mean that the shapes maintain separate geometries but are treated as a single feature with a single set of attributes by the GIS:

A fundamental geospatial operation is checking to see whether a point is inside a polygon. This one operation is the atomic building block of many different types of spatial queries. If the point is on the boundary of the polygon, it is considered inside. Very few spatial queries exist that do not rely on this calculation in some way. However, it can be very slow on a large number of points.

The most common and efficient algorithm to detect if a point is inside a polygon is called the ray casting algorithm. First, a test is performed to see if the point is on the polygon boundary. Next, the algorithm draws a line from the point in question in a single direction. The program counts the number of times the line crosses the polygon boundary until it reaches the bounding box of the polygon. The bounding box is the smallest box that can be drawn around the entire polygon. If the number is odd, the point is inside. If the number of boundary intersections is even, the point is outside:

The union operation is less common but very useful to combine two or more overlapping polygons in a single shape. It is similar to dissolve, but in this case, the polygons are overlapping as opposed to being adjacent. Usually, this operation is used to clean up automatically-generated feature datasets from remote sensing operations:

A join or SQL join is a database operation used to combine two or more tables of information. Relational databases are designed to avoid storing redundant information for one-to-many relationships. For example, a U.S. state may contain many cities. Rather than creating a table for each state containing all of its cities, a table of states with numeric IDs is created, while a table for all the cities in every state is created with a state numeric ID. In a GIS, you can also have spatial joins that are part of the spatial extension software for a database. In spatial joins, combine the attributes to two features in the same way that you do in a SQL join, but the relation is based on the spatial proximity of the two features. To follow the previous cities example, we could add the county name that each city resides in using a spatial join. The cities layer could be loaded over a county polygon layer whose attributes contain the county name. The spatial join would determine which city is in which county and perform a SQL join to add the county name to each city's attribute row.

In geospatial analysis, there are several general rules of thumb regarding polygons that are different from mathematical descriptions of polygons:

Different geospatial software packages and libraries handle exceptions to these rules differently and can lead to confusing errors or software behaviors. The safest route is to make sure that your polygons obey these rules. There is one more important piece of information about polygons. A polygon is by definition a closed shape, which means that the first and last vertices of the polygon are identical. Some geospatial software will throw an error if you haven't explicitly duplicated the first point as the last point in the polygon dataset. Other software will automatically close the polygon without complaining. The data format that you use to store your geospatial data may also dictate how polygons are defined. This issue is a gray area and so it didn't make the polygon rules, but knowing this quirk will come in handy someday when you run into an error that you can't explain easily.

Band math is multidimensional array mathematics. In array math, arrays are treated as single units, which are added, subtracted, multiplied, and divided. However, in an array, the corresponding numbers in each row and column across multiple arrays are computed simultaneously. These arrays are termed matrices and computations involving matrices are the focus of linear algebra.

Change detection is the process of taking two images of the same location at different times and highlighting the changes. A change can be due to the addition of something on the ground, such as a new building or the loss of a feature such as coastal erosion. There are many algorithms to detect changes among images and also determine qualitative factors such as how long ago the change took place. The following image from a research project by the U.S. Oak Ridge National Laboratory (ORNL) shows rainforest deforestation between 1984 and 2000 in the state of Rondonia, Brazil. Colors are used to show how recently the forest was cut. Green represents virgin rainforests, white is a forest cut within two years of the end of the date range, red is within 22 years, and the other colors fall in between as described in the legend:

A histogram is the statistical distribution of values in a dataset. The horizontal axis represents a unique value in a dataset while the vertical axis represents the frequency of this unique value in the raster. A histogram is a key operation in most raster processing. It can be used for everything from enhancing contrast in an image to serving as a basis for object classification and image comparison. The following example from NASA shows a histogram for a satellite image that has been classified into different categories representing the underlying surface features:

Feature extraction is the manual or automatic digitization of features in an image to points, lines, or polygons. This process serves as the basis for the vectorization of images in which a raster is converted to a vector dataset. An example of feature extraction is extracting a coastline from a satellite image and saving it as a vector dataset. If this extraction is performed over several years, you could monitor the erosion or other changes along this coastline.

Objects on the Earth reflect different wavelengths of light depending on the materials that they are made of. In remote sensing, analysts collect wavelength signatures for specific types of land cover (for example, concrete) and build a library for a specific area. A computer can then use this library to automatically locate classes in the library in a new image of the same area.

In an unsupervised classification, a computer groups pixels with similar reflectance values in an image without any other reference information other than the histogram of the image.

As stated earlier, this book assumes that you have some basic knowledge of Python. The examples in this book are based on Python 3.4.3, which you can download here:

https://www.python.org/downloads/release/python-343/

The only module used in the following example is the turtle module that provides a very simple graphics engine based on the Tkinter library included with Python. If you used the installers for Windows or Mac OS X, the Tkinter library should be included already. If you compiled Python yourself or are using a distribution from somewhere besides Python.org, then make sure that you can import the turtle module by typing the following at a command prompt to run the turtle demo script:

python –m turtle

If your Python distribution does not have Tkinter, you can find information on installing it from the following page. The information is for Python 2.3 but the process is still the same:

http://tkinter.unpythonic.net/wiki/How_to_install_Tkinter

The official Python wiki page for Tkinter can be found here:

https://wiki.python.org/moin/TkInter

The documentation for Tkinter is in the Python Standard Library documentation that can be found at https://docs.python.org/2/library/tkinter.html.

If you are new to Python, Dive into Python is a free online book, which covers all the basics of Python and will bring you up to speed. For more information, refer to http://www.diveintopython.net.

The code is divided into two different sections. The first is the data model section and the second is the map renderer that draws this data. For the data model, we will use simple Python lists. A Python list is a native data type, which serves as a container for other Python objects in a specified order. Python lists can contain other lists and are great for simple data structures. They also map well to more complex structures or even databases if you decide you want to develop your script further.

The second portion of the code will render the map using the Python turtle graphics engine. We will have only one function in the GIS that converts the world coordinates, in this case, longitude and latitude, into pixel coordinates. All graphics engines have an origin point of (0,0) that is usually in the top-left or lower-left corner of the canvas. Turtle graphics are designed to teach programming visually. The turtle graphics canvas uses an origin of (0,0) in the center, similar to a graphing calculator. The following image illustrates the type of Cartesian graph that the turtle module uses. In the following graph, some points are plotted in both positive and negative space:

This also means that the turtle graphics engine can have negative pixel coordinates, which is uncommon for graphics canvases. However, for this example, the turtle module is the quickest and simplest way to render our map.

Step by step

You can run this program interactively in the Python interpreter or you can save the complete program as a script and run it. The Python interpreter is an incredibly powerful way to learn new concepts because it gives you real-time feedback on errors or unexpected program behavior. You can easily recover from these issues and try something else until you get the results that you want:

1. In Python, you usually import modules at the beginning of the script so we'll import the turtle module first. We'll use Python's import as feature to assign the module the name t to save space and time when typing turtle commands:

   import turtle as t
   

2. Next, we'll set up the data model starting with some simple variables that allow us to access list indexes by name instead of numbers to make the code easier to follow. Python lists index the contained objects starting with the number 0. So, if we want to access the first item in a list called myList, we would reference it as follows:

   myList[0]
   

3. To make our code easier to read, we can also use a variable name assigned to commonly used indexes:

   firstItem = 0
   myList[firstItem]
   

   Note

   Downloading the example code

   You can download the example code files for all the Packt books that you have purchased from your account at http://www.packtpub.com. If you purchased this book elsewhere, you can visit http://www.packtpub.com/support and register in order to have the files e-mailed to you directly.

   In computer science, assigning commonly used numbers to an easy-to-remember variable is a common practice. These variables are called constants.

   So, for our example, we'll assign constants for some common elements used for all the cities. All cities will have a name, one or more points, and a population count:

   NAME = 0
   POINTS = 1
   POP = 2
   

4. Now, we'll set up the data for Colorado as a list with name, polygon points, and population. Note that the coordinates are a list within a list:

   state = ["COLORADO", [[-109, 37],[-109, 41],[-102, 41],[-102, 37]], 5187582]
   

5. The cities will be stored as nested lists. Each city's location consists of a single point as a longitude and latitude pair. These entries will complete our GIS data model. We'll start with an empty list called cities and then append the data to this list for each city:

   cities = []
   cities.append(["DENVER",[-104.98, 39.74], 634265])
   cities.append(["BOULDER",[-105.27, 40.02], 98889])
   cities.append(["DURANGO",[-107.88,37.28], 17069])
   

6. We will now render our GIS data as a map by first defining a map size. The width and height can be anything that you want depending on your screen resolution:

   map_width = 400
   map_height = 300
   

7. In order to scale the map to the graphics canvas, we must first determine the bounding box of the largest layer, which is the state. We'll set the map's bounding box to a global scale and reduce it to the size of the state. To do so, we'll loop through the longitude and latitude of each point and compare it with the current minimum and maximum x and y values. If it is larger than the current maximum or smaller than the current minimum, we'll make this value the new maximum or minimum, respectively:

   minx = 180
   maxx = -180
   miny = 90
   maxy = -90 
   forx,y in state[POINTS]:
   if x < minx: minx = x
     elif x > maxx: maxx = x
     if y < miny: miny = y
     elif y > maxy: maxy = y
   

8. The second step to scaling is to calculate a ratio between the actual state and the tiny canvas that we will render it on. This ratio is used for coordinate to pixel conversion. We get the size along the x and y axes of the state and then we divide the map width and height by these numbers to get our scaling ratio:

   dist_x = maxx - minx
   dist_y = maxy - miny
   x_ratio = map_width / dist_x
   y_ratio = map_height / dist_y
   

9. The following function, called convert(), is our only function in SimpleGIS. It transforms a point in the map coordinates from one of our data layers to pixel coordinates using the previous calculations. You'll notice that, at the end, we divide the map width and height in half and subtract it from the final conversion to account for the unusual center origin of the turtle graphics canvas. Every geospatial program has some form of this function:

   def convert(point):
     lon = point[0]
     lat = point[1]
     x = map_width - ((maxx - lon) * x_ratio)
     y = map_height - ((maxy - lat) * y_ratio)
     # Python turtle graphics start in the 
     # middle of the screen
     # so we must offset the points so they are centered
     x = x - (map_width/2)
     y = y - (map_height/2)
     return [x,y]
   

10. Now for the exciting part! We're ready to render our GIS as a thematic map. The turtle module uses the concept of a cursor called a pen. Moving the cursor around the canvas is exactly the same as moving a pen around a piece of paper. The cursor will draw a line when you move it. So, you'll notice that throughout the code, we use the t.up() and t.down() commands to pick the pen up when we want to move to a new location and put it down when we're ready to draw. We have some important steps in this section. As the border of Colorado is a polygon, we must draw a line between the last point and first point to close the polygon. We can also leave out the closing step and just add a duplicate point to the Colorado dataset. Once we draw the state, we'll use the write() method to label the polygon:

    t.up()
    first_pixel = None
    for point in state[POINTS]:
      pixel = convert(point)
      if not first_pixel:
        first_pixel = pixel
      t.goto(pixel)
      t.down()
    t.goto(first_pixel)
    t.up()
    t.goto([0,0])
    t.write(state[NAME], align="center", font=("Arial",16,"bold"))
    

11. If we were to run the code at this point, we would see a simplified map of the state of Colorado, as shown in the following screenshot:
    

    Note

    If you do try to run the code, you'll need to temporarily add the following line at the end, or the Tkinter window will close as soon as it finishes drawing:

    t.done()
    

12. Now, we'll render the cities as point locations and label them with their names and population. As the cities are a group of features in a list, we'll loop through them to render them. Instead of drawing lines by moving the pen around, we'll use the turtle dot() method to plot a small circle at the pixel coordinate returned by our SimpleGISconvert() function. We'll then label the dot with the city's name and add the population. You'll notice that we must convert the population number to a string in order to use it in the turtle write() method. To do so, we use Python's built-in function, str():

      #for city in cities:
      pixel = convert(city[POINTS])
      t.up()
      t.goto(pixel)
      # Place a point for the city
      t.dot(10)
      # Label the city
      t.write(city[NAME] + ", Pop.: " + str(city[POP]), align="left")
      t.up()
    

13. Now we will perform one last operation to prove that we have created a real GIS. We will perform an attribute query on our data to determine which city has the largest population. Then, we'll perform a spatial query to see which city lies the furthest west. Finally, we'll print the answers to our questions on our thematic map page safely out of the range of the map.
14. As our query engine, we'll use Python's built-in min() and max() functions. These functions take a list as an argument and return the minimum and maximum values of this list. These functions have a special feature called a key argument that allows you to sort complex objects. As we are dealing with nested lists in our data model, we'll take advantage of the key argument in these functions. The key argument accepts a function that temporarily alters the list for evaluation before a final value is returned. In this case, we want to isolate the population values for comparison and then the points. We could write a whole new function to return the specified value, but we can use Python's lambda keyword instead. The lambda keyword defines an anonymous function that is used inline. Other Python functions can be used inline, for example, the string function, str(), but they are not anonymous. This temporary function will isolate our value of interest.
15. So our first question is, which city has the largest population?

    biggest_city = max(cities, key=lambda city:city[POP])
    t.goto(0,-200)
    t.write("The biggest city is: " +  biggest_city[NAME])
    

16. The next question is, which city lies the furthest west?

    western_city = min(cities, key=lambda city:city[POINTS])
    t.goto(0,-220)
    t.write("The western-most city is: " + western_city[NAME])
    

17. In the preceding query, we use Python's built in min() function to select the smallest longitude value and this works because we represented our city locations as longitude and latitude pairs. It is possible to use different representations for points including possible representations where this code would need modification to work correctly. However, for our SimpleGIS, we are using a common point representation to make it as intuitive as possible.

These last two commands are just for clean up purposes. First, we hide the cursor. Then, we call the turtle done() method, which will keep the turtle graphics window with our map open until we choose to close it using the close handle at the top of the window.

t.pen(shown=False)
t.done()

Whether you followed along using the Python interpreter or you ran the complete program as a script, you should see the following map rendered in real time:

Congratulations! You have followed in the footsteps of Paleolithic hunters, the father of GIS Dr. Roger Tomlinson, geospatial pioneer Howard Fisher, and game-changing humanitarian programmers to create a functional, extensible, and technically complete geographic information system. It took less than 60 lines of pure Python code! You will be hard-pressed to find a programming language that can create a complete GIS using only its core libraries in such a finite amount of readable code as Python. Even if you did, it is highly unlikely that the language would survive the geospatial Python journey that you'll take through the rest of this book.

As you can see, there is lots of room for expansion of SimpleGIS. Here are some other ways that you might expand this simple tool using the reference material for Tkinter and Python linked at the beginning of this section:

• Create an overview map in the top-right corner with a U.S. border outline and Colorado's location in the U.S.
• Add color for visual appeal and further clarity
• Create a map key for different features
• Make a list of states and cities and add more states and cities
• Add a title to the map
• Create a bar chart to compare population numbers visually

The possibilities are endless. SimpleGIS can also be used as a way to quickly test and visualize geospatial algorithms that you come across. If you want to add more data layers, you can create more lists, but these lists would become difficult to manage. In this case, you can use another Python module included in the standard distribution. The SQLite module provides a SQL-like database in Python that can be saved to disk or run in memory.