Python: A Simple GIS Model of Colorado

Hello Readers,


Since I have not posted with Python in awhile, this is a great opportunity to explore geographic information systems (GIS) with Python! We shall start mapping with a simple graphic of the Colorado state with a few cities using turtle graphics.

Though we do not use fancy graphics here, the resulting map below introduces the GIS capability of pure Python. For further GIS work with a variety of data formats such as raster and vector data, we will use the geospatial data abstraction library (GDAL). However, with the simple turtle graphics, we can create the map showing the 3 cities in Colorado and their populations:

In addition to the city names and their populations, we also print the largest city and the western most city by longitude from querying the population and coordinates for the maximum and minimum, respectively. The largest city has the largest population, and the western most city has the lowest longitude, being negative in the western hemisphere.

The code is displayed below as a reference:

 # GIS in Python Part 1
 
import turtle as t
 
## Section 1: Data Model
NAME = 0
POINTS = 1
POP = 2
state = ["COLORADO", [[-109, 37],[-109, 41],[-102,41],[-102,37]], 5187582]
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])
map_width = 400
map_height = 300
 
 # bounding box for COLORADO
minx = 180
maxx = -180
miny = 90
maxy = -90
for x,y in state[POINTS]:
    if x < minx: minx = x
    elif x > maxx: maxx = x
    if y < miny: miny = y
    elif y > maxy: maxy = y
dis_x = maxx - minx
dis_y = maxy - miny
x_ratio = map_width / dis_x
y_ratio = map_height / dis_y
 
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 middel of screen
    # must offset points to center
    x = x - (map_width/2)
    y = y - (map_height/2)
    return [x,y]
 
## Section 2: Map Renderer
# COLORADO
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"))

# CITIES
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

# attribute query, use keyword lambda
biggest_city = max(cities, key=lambda city:city[POP])
t.goto(0,-200)
t.write("The biggest city is: " + biggest_city[NAME])
western_city = min(cities, key=lambda city:city[POINTS])
t.goto(0,-220)
t.write("The western-most city is: " + western_city[NAME])
t.goto(0, 200)
t.write("MAP OF COLORADO", align="center", font=("Arial", 20, "bold"))
t.pen(shown=False)
t.done()

There will be more python posts with GDAL and colorful visualizations! Stay tuned, and in the meantime, take a look at 'greatest ever' infographic :



Thanks for reading,

Wayne
@beyondvalence