### Sidebar

spatial-ecology.org

Trainings:

Learn:

Data:

Community:
teachers
students
projects
Matera 2015
Vancouver 2015
Santa Barbara 2015
Site design
Hands on:
Installations

Donations
USD

GBP

wiki:python:basicpython2

# Geospatial Python crash course, part 2

## Some examples of code

Some example of reasonable complicated programs:

| 01/SimpleGIS.py
```import turtle as t

# DATA MODEL
# All layers will have a name, 1+ points, and population count
NAME = 0
POINTS = 1
POP = 2

# Create the state layer
state = ["COLORADO", [[-109, 37], [-109, 41], [-102, 41], [-102, 37]], 5187582]

# Cities layer list
# city = [name, [point], population]
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 GRAPHICS RENDERING
map_width = 800
map_height = 500

# State Bounding Box
# Use Python min/max function to get bounding box
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
# Get earth distance on each axis
dist_x = maxx - minx
dist_y = maxy - miny

# Scaling ratio each axis
# to map points from world to screen
x_ratio = map_width / dist_x
y_ratio = map_height / dist_y

# Function to convert lat/lon to screen coordinates
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]

# Draw the state
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()
# Go back to the first point
t.goto(first_pixel)
# Label the state
t.up()
t.goto([0,0])
t.write(state[NAME], align="center", font=("Arial",16,"bold"))

# Draw the 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()

# Perform an attribute query
# Question: Which city has the largest population?
# Write the result but make sure it's under the map
biggest_city = max(cities, key=lambda city:city[POP])
t.goto(0, -1*((map_height/2)+20))
t.write("The biggest city is: " +  biggest_city[NAME])

# Perform a spatial query
# Question: Which is the western most city?
# Write the result but make sure it's under the other question
western_city = min(cities, key=lambda city:city[POINTS])
t.goto(0, -1*((map_height/2)+40))
t.write("The western-most city is: " + western_city[NAME])

# Hide our map pen
t.pen(shown=False)
t.done()```
02/fmtDecode.py
```#!/usr/bin/python3

"""
fmtDecode.py - Simple program to manually decode binary formats.
"""

import struct
import pickle
import os

def export():
print("Saving results")
out = None
if cached:
out = open(oname, "a")
else:
out = open(oname, "w")
for record in fileDesc:
for field in record:
out.write("%s\t" % field)
out.write("\n")
out.close()
pickle.dump(cached, open(pickleJar, "w"))

fileDesc = []
files = os.listdir(".")
count = 1
print("Available Files:")

for f in files:
print(" %s. %s" % (count, f))
count += 1

fnum = input("Enter the number of the file to decode: ")
fname = files[int(fnum)-1]
base = os.path.splitext(fname)[0]

pickleJar = "%s.p" % base

cached = []

if os.path.exists(pickleJar):
print("Cached session available.\n")
useCache = input("Use it? Yes (Y), No (N)?")
if "y" in useCache.lower() or useCache == "":
else: cached = []

oname = "%s_decode.txt" % base

f = open(fname, "rb")
loc = f.tell()
f.seek(0,2)
eof = f.tell()
f.seek(0)
prev = 0

if len(cached)>0:
print("Using cache...")
f.seek(cached[-1])
prev = cached[-2]

print("Starting at byte %s..." % f.tell())

try:
formats = {"char":{"format":"c","len":1},
"signed char":{"format":"b","len":1},
"unsigned char":{"format":"B","len":1},
"_Bool":{"format":"?","len":1},
"short":{"format":"h","len":2},
"unsigned short":{"format":"h","len":2},
"int":{"format":"i","len":4},
"unsigned int":{"format":"I","len":4},
"long":{"format":"l","len":4},
"unsigned long":{"format":"L","len":4},
"long long":{"format":"q","len":8},
"unsigned long long":{"format":"Q","len":8},
"float":{"format":"f","len":4},
"double":{"format":"d","len":8}}

while f.tell() < eof:
record = []
start = f.tell()
record.append("%s\t" % start)
cached.append(start)
fields = []
print()
count = 1
try:
# Little endian formats
for fmt in formats:
form = formats[fmt]["format"]
bytes = formats[fmt]["len"]
field = struct.unpack("<%s" % form, f.read(bytes))
print("%s. Little %s: %s" % (count, fmt, field))
count += 1
f.seek(start)
fields.append([str(field[0]), fmt, "little", str(bytes)])
except: pass

try:
# Big endian formats
for fmt in formats:
form = formats[fmt]["format"]
bytes = formats[fmt]["len"]
field = struct.unpack(">%s" % form, f.read(bytes))
print("%s. Big %s: %s" % (count, fmt, field))
count += 1
f.seek(start)
fields.append([str(field[0]), fmt, "big", str(bytes)])
except: pass

print("%s. Go back to previous" % count)
print()
print("Current location: %s" % f.tell())
choice = input("Enter the number of one of the above options: ")
choice = int(choice.strip())
desc = ""
if choice != count:
desc = input("Enter a field description: ")
record.append("%s\t" % desc)
record.append("%s\t" % fields[choice-1][0])
record.append("%s\t" % fields[choice-1][1])
record.append("%s\t" % fields[choice-1][2])
f.seek(start + int(fields[choice-1][3]))
prev = start
fileDesc.append(record)
elif choice == count:
f.seek(prev)
print("Going back to previous field.")
f.close()
export()
except KeyboardInterrupt:
print()
reverse = input("How many records back? ")
for i in range(reverse):
cached.pop()
pickle.dump(cached, open(pickleJar, "w"))
print("The program will exit.  Restart and use cached version.")

except:
export()

```
02/GeoJSON.py
```#!/usr/bin/python3
gc = { "type": "GeometryCollection",
"geometries": [
{ "type": "Point",
"coordinates": [-89.33, 30.0]
},
{ "type": "LineString",
"coordinates": [ [-89.33, 30.30], [-89.36, 30.28] ]
}
]
}
print(gc)

```
02/structDemo.py
```#!/usr/bin/python3
"""
structDemo.py - demonstrate using the struct module
by reading the bounding box from a shapefile.
"""
import struct
# Open the shapefile
f = open("../files/hancock.shp","rb")
# Go to the start of the
# bounding box coordinates
f.seek(36)
# in little endian format
# Read all 4 values at once
f.seek(36)

## Playing with network

Python is a general purpose language, so it can be used for network applications.

04/01_urlretrieve.py
```#!/usr/bin/python3

# Retrieve a file using urllib
import urllib.request
url = "http://spatial-ecology.net/dokuwiki/lib/exe/fetch.php?media=wiki:python:hancock.zip"
fileName = "hancock.zip"
urllib.request.urlretrieve(url, fileName)```
04/02_earthquake-data.py
```#!/usr/bin/python3

# Read data from a url - in this case recent USGS earthquakes
import urllib.request
url = "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_hour.csv"
earthquakes = urllib.request.urlopen(url)
# Read the first two earthquakes
# Iterate through the rest
for record in earthquakes:
print(record)```
04/03_ftp.py
```#!/usr/bin/python3

# Read tsunami monitoring data via ftp
import ftplib
server = "ftp.ngdc.noaa.gov"
dir = "hazards/DART/20070815_peru"
fileName = "21415_from_20070727_08_55_15_tides.txt"
ftp = ftplib.FTP(server)
ftp.cwd(dir)
out = open(fileName, "wb")
ftp.retrbinary("RETR " + fileName, out.write)
out.close()
dart = open(fileName)
for line in dart:
if "LAT," in line:
print(line)
break```
04/04_ftp-url.py
```#!/usr/bin/python3

# Read FTP data via the web
import urllib.request
ftpURL = "ftp://anonymous:anonymous@"
server = "ftp.ngdc.noaa.gov"
dir = "hazards/DART/20070815_peru"
fileName = "21415_from_20070727_08_55_15_tides.txt"
dart = urllib.request.urlopen(ftpURL + server + "/" + dir + "/" + fileName)
for line in dart:
if b"LAT," in line:
print(line)
break```

## Zip files and others

It can be used to manipulate archive files:

04/05_zipfile.py
```#!/usr/bin/python3

# Unzip a shapefile
import zipfile
zip = open("../files/hancock.zip", "rb")
zipShape = zipfile.ZipFile(zip)
shpName, shxName, dbfName = zipShape.namelist()
shpFile = open(shpName, "wb")
shxFile = open(shxName, "wb")
dbfFile = open(dbfName, "wb")
shpFile.close()
shxFile.close()
dbfFile.close()```
04/06_zipfile-loop.py
```#!/usr/bin/python3

# Unzip a shapefile with a for loop
import zipfile
zip = open("../files/hancock.zip", "rb")
zipShape = zipfile.ZipFile(zip)
for fileName in zipShape.namelist():
out = open(fileName, "wb")
out.close()```
04/07_tarfile-create.py
```#!/usr/bin/python3

# Add a shapefile to a tar archive
import tarfile
tar = tarfile.open("hancock.tar.gz", "w:gz")
tar.close()```
04/08_tarfile-extract.py
```#!/usr/bin/python3

# Extract a shapefile from a gzipped tar archive
import tarfile
tar = tarfile.open("hancock.tar.gz", "r:gz")
tar.extractall()
tar.close()```
04/09_cloud-zipfile.py
```#!/usr/bin/python3

# Extract a zipped shapefile via a url
import urllib.request
import zipfile
import io
import struct
url = "http://spatial-ecology.net/dokuwiki/lib/exe/fetch.php?media=wiki:python:hancock.zip"
cloudshape = urllib.request.urlopen(url)
zipshape = zipfile.ZipFile(memoryshape)
# Access Python string as an array
struct.unpack("<dddd", cloudshp[36:68])```

## Playing with XML

Many useful file formats are of XML type (e.g. KML) and Python has a good deals of modules to manage them:

04/10_kml-minidom.py
```#!/usr/bin/python3

# Parse KML and count placemarks
from xml.dom import minidom
kml = minidom.parse("../files/time-stamp-point.kml")
Placemarks = kml.getElementsByTagName("Placemark")
# Count placemarks
print(len(Placemarks))
# Check the first one
print(Placemarks[0])
print(Placemarks[0].toxml())
# Extract the coordinates
coordinates = Placemarks[0].getElementsByTagName("coordinates")
point = coordinates[0].firstChild.data
print(point)
# Extract x,y,z values as floats
x,y,z = point.split(",")
print(x)
print(y)
print(z)
x = float(x)
y = float(y)
z = float(z)
print(x,y,z)
# Use list comprehensions for efficiency
x,y,z = [float(c) for c in point.split(",")]
print(x,y,z)```
04/12_make-kml-strings.py
```#!/usr/bin/python3

# Build KML manually using strings
xml = """<?xml version="1.0" encoding="utf-8"?>"""
xml += """<kml xmlns="http://www.opengis.net/kml/2.2">"""
xml += """  <Placemark>"""
xml += """    <name>Office</name>"""
xml += """    <description>Office Building</description>"""
xml += """    <Point>"""
xml += """      <coordinates>"""
xml += """        -122.087461,37.422069"""
xml += """      </coordinates>"""
xml += """    </Point>"""
xml += """  </Placemark>"""
xml += """</kml>"""
print(xml)```
04/13_make-kml-etree.py
```#!/usr/bin/python3

# Build kml with eTree
try:
import xml.etree.cElementTree as ET
except ImportError:
import xml.etree.ElementTree as ET

root = ET.Element("kml")
root.attrib["xmlns"] = "http://www.opengis.net/kml/2.2"
placemark = ET.SubElement(root, "Placemark")
office = ET.SubElement(placemark, "name")
office.text = "Office"
point = ET.SubElement(placemark, "Point")
coordinates = ET.SubElement(point, "coordinates")
coordinates.text = "-122.087461,37.422069"
tree = ET.ElementTree(root)
tree.write("placemark.kml", xml_declaration=True,encoding='utf-8',method="xml")```
04/14_gpx-soup.py
```#!/usr/bin/python3

# Parse broken GPX data with BeautifulSoup
from BeautifulSoup import BeautifulStoneSoup
gpx = open("../files/broken_data.gpx")
# Check the first track point
print(soup.trkpt)
# Find the rest of the track points and count them
tracks = soup.findAll("trkpt")
print(len(tracks))
# Save the fixed xml file
fixed = open("fixed_data.gpx", "w")
fixed.write(soup.prettify())
fixed.close()```

## Playing with some geospatial modules

Some very basic geospatial snippets of code;

04/15_shapely-wkt.py
```#!/usr/bin/python3

import shapely.wkt
wktPoly = "POLYGON((0 0,4 0,4 4,0 4,0 0),(1 1, 2 1, 2 2, 1 2,1 1))"
print(poly.area)
print(poly.wkt)```
04/16_ogr-wkt.py
```#!/usr/bin/python3

# Convert a shapefile to WKT using ogr
from osgeo import ogr
shape = ogr.Open("../files/polygon.shp")
layer = shape.GetLayer()
feature = layer.GetNextFeature()
geom = feature.GetGeometryRef()
wkt = geom.ExportToWkt()
# View the WKT
print(wkt)
# Ingest the WKT we made and check the envelope
poly = ogr.CreateGeometryFromWkt(wkt)
print(poly.GetEnvelope())```
04/17_json.py
```#!/usr/bin/python3

# Parse GeoJson data
jsdata = """{
"type": "Feature",
"id": "OpenLayers.Feature.Vector_314",
"properties": {},
"geometry": {
"type": "Point",
"coordinates": [ 97.03125, 39.7265625 ]
},
"crs": {
"type": "name",
"properties": {
"name": "urn:ogc:def:crs:OGC:1.3:CRS84"
}
}
}"""
# Try to eval() the data
point = eval(jsdata)
print(point["geometry"])
# Use the json module
import json
# Parse and then dump GeoJSON
print(json.dumps(pydata))```
04/18_geojson.py
```#!/usr/bin/python3

# Read and write GeoJson using the geojson module
import geojson
p = geojson.Point([-92, 37])
geojs = geojson.dumps(p)
print(geojs)
# Use __geo_interface__ between geojson and shapely
from shapely.geometry import asShape
point = asShape(p)
print(point.wkt)```
04/19_ogr.py
```#!/usr/bin/python3

#Examine a shapefile with ogr

from osgeo import ogr
# open the shapefile
shp = ogr.Open("../files/point.shp")
# Get the layer
layer = shp.GetLayer()
# Loop through the features
# and print information about them
for feature in layer:
geometry = feature.GetGeometryRef()
print(geometry.GetX(), geometry.GetY(), feature.GetField("FIRST_FLD"))```
04/20_pyshp.py
```#!/usr/bin/python3

# Examine a shapefile with pyshp                                                                                     import shapefile
import shapefile

for feature in shp.shapeRecords():
point = feature.shape.points[0]
rec = feature.record[0]
print(point[0], point[1], rec)```
04/21_dbfpy.py
```#!/usr/bin/python3

# Examine and update a dbf file with dbfpy
from dbfpy import dbf
db = dbf.Dbf("../files/GIS_CensusTract_poly.dbf")
print db[0]
rec = db[0]
field = rec["POPULAT10"]
rec["POPULAT10"] = field
rec["POPULAT10"] = field+1
rec.store()
del rec
print(db[0]["POPULAT10"])```
04/22_shapely.py
```#!/usr/bin/python3

# Create a polygon buffer with shapely

from shapely import wkt, geometry
# Create a simple wkt polygon string
wktPoly = "POLYGON((0 0,4 0,4 4,0 4,0 0))"
# Load the polygon into shapely
# Check the area
print(poly.area)
# Create a buffer
buf =(poly.buffer(5.0)
# Get the buffer area
print(buf.area)
# Compute the difference between the two
print(buf.difference(poly).area)```
04/23_gdal.py
```#!/usr/bin/python3

# Open a raster with gal
from osgeo import gdal
raster = gdal.Open("../files/SatImage.tif")
print(raster.RasterCount)
print(raster.RasterXSize)
print(raster.RasterYSize)```
04/24_numpy.py
```#!/usr/bin/python3

# Numpy/gdalnumeric - Read an image, extract a band, save a new image
from osgeo import gdalnumeric
band1 = srcArray[0]
gdalnumeric.SaveArray(band1, "band1.jpg", format="JPEG")```
04/25_PIL.py
```#!/usr/bin/python3

# Rasterize a shapefile with PIL
import shapefile
from PIL import Image, ImageDraw
xdist = r.bbox[2] - r.bbox[0]
ydist = r.bbox[3] - r.bbox[1]
iwidth = 400
iheight = 600
xratio = iwidth/xdist
yratio = iheight/ydist
pixels = []
for x,y in r.shapes()[0].points:
px = int(iwidth - ((r.bbox[2] - x) * xratio))
py = int((r.bbox[3] - y) * yratio)
pixels.append((px,py))
img = Image.new("RGB", (iwidth, iheight), "white")
draw = ImageDraw.Draw(img)
draw.polygon(pixels, outline="rgb(203, 196, 190)", fill="rgb(198, 204, 189)")
img.save("hancock.png")```
04/26_PNGC.py
```#!/usr/bin/python3

# Rasterize a shapefile with PNGCanvas
import shapefile
import pngcanvas
xdist = r.bbox[2] - r.bbox[0]
ydist = r.bbox[3] - r.bbox[1]
iwidth = 400
iheight = 600
xratio = iwidth/xdist
yratio = iheight/ydist
pixels = []
for x,y in r.shapes()[0].points:
px = int(iwidth - ((r.bbox[2] - x) * xratio))
py = int((r.bbox[3] - y) * yratio)
pixels.append([px,py])
c = pngcanvas.PNGCanvas(iwidth,iheight)
c.polyline(pixels)
f = file("hancock_pngcvs.png", "wb")
f.write(c.dump())
f.close()```
04/27_PDF.py
```#!/usr/bin/python3

# Create a pdf with a map image

import fpdf

# PDF constructor:
# Portrait, millimeter units, A4 page size
pdf=fpdf.FPDF("P", "mm", "A4")
# create a new page
# Set font: arial, bold, size 20
pdf.set_font('Arial','B',20)
# Create a new cell: 160 x 25mm, text contents, no border, centered
pdf.cell(160,25,'Hancock County Boundary', border=0, align="C")
pdf.image("hancock.png",25,50,150,160)
pdf.output('map.pdf','F')```

## Using OGR

Let us see some very basic useful snippets of OGR code which could be combined together to build more complicated things:

03/01_ogr.py
```#!/usr/bin/python

import osgeo.ogr
import os.path
import sys

mypath = os.path.dirname(os.path.realpath(sys.argv[0]))
shapefile_name = os.path.join(mypath,"../files/tl_2012_us_cbsa.shp")
shapefile = osgeo.ogr.Open(shapefile_name)
numLayers = shapefile.GetLayerCount()

print ("Shapefile contains %d layers" % numLayers)
print

for layerNum in range(numLayers):
layer = shapefile.GetLayer(layerNum)
spatialRef = layer.GetSpatialRef().ExportToProj4()
numFeatures = layer.GetFeatureCount()
print "Layer %d has spatial reference %s" % (layerNum, spatialRef)
print "Layer %d has %d features" % (layerNum, numFeatures)
print

for featureNum in range(numFeatures):
feature = layer.GetFeature(featureNum)
featureName = feature.GetField("NAME")

print "Feature %d has name %s" % (featureNum, featureName)```
03/02_ogr.py
```#!/usr/bin/python
import osgeo.ogr
import os.path
import sys

mypath = os.path.dirname(os.path.realpath(sys.argv[0]))
shapefile_name = os.path.join(mypath,"../files/tl_2012_us_cbsa.shp")
shapefile = osgeo.ogr.Open(shapefile_name)
layer = shapefile.GetLayer(0)
feature = layer.GetFeature(2)

print "Feature 2 has the following attributes:"
print

attributes = feature.items()

for key,value in attributes.items():
print " %s = %s" % (key, value)

geometry = feature.GetGeometryRef()
geometryName = geometry.GetGeometryName()

print
print "Feature's geometry data consists of a %s" % geometryName```
03/03_ogr.py
```#!/usr/bin/python

import osgeo.ogr
import os.path
import sys

def analyzeGeometry(geometry, indent=0):
s = []
s.append("  " * indent)
s.append(geometry.GetGeometryName())
if geometry.GetPointCount() > 0:
s.append(" with %d data points" % geometry.GetPointCount())
if geometry.GetGeometryCount() > 0:
s.append(" containing:")

print "".join(s)

for i in range(geometry.GetGeometryCount()):
analyzeGeometry(geometry.GetGeometryRef(i), indent+1)

mypath = os.path.dirname(os.path.realpath(sys.argv[0]))
shapefile_name = os.path.join(mypath,"../files/tl_2012_us_cbsa.shp")
shapefile = osgeo.ogr.Open(shapefile_name)
layer = shapefile.GetLayer(0)
feature = layer.GetFeature(55)
geometry = feature.GetGeometryRef()

analyzeGeometry(geometry)```
03/04_ogr.py
```#!/usr/bin/python

import osgeo.ogr
import os.path

import sys

def findPoints(geometry, results):
for i in range(geometry.GetPointCount()):
x,y,z = geometry.GetPoint(i)
if results['north'] == None or results['north'][1] < y:
results['north'] = (x,y)
if results['south'] == None or results['south'][1] > y:
results['south'] = (x,y)

for i in range(geometry.GetGeometryCount()):
findPoints(geometry.GetGeometryRef(i), results)

mypath = os.path.dirname(os.path.realpath(sys.argv[0]))
shapefile_name = os.path.join(mypath,"../files/tl_2012_us_cbsa.shp")
shapefile = osgeo.ogr.Open(shapefile_name)
layer = shapefile.GetLayer(0)
feature = layer.GetFeature(55)
geometry = feature.GetGeometryRef()

results = {'north' : None,
'south' : None}

findPoints(geometry, results)

print "Northernmost point is (%0.4f, %0.4f)" % results['north']
print "Southernmost point is (%0.4f, %0.4f)" % results['south']```

Playing with projection in OGR and UTM modules:

05/06-utm2latlon.py
```"""Convert utm to lat/lon"""
import utm
y = 479747.0453210057
x = 5377685.825323031
zone = 32
band = 'U'
print utm.to_latlon(y,x,zone,band)```
05/07-latlon2utm.py
```"""Convert lat/lon to utm"""
import utm
print utm.from_latlon(48.55199390882121, 8.725555729071763)```

Hands on work: write a couple of programs to convert coordinates taken from input (hint: sys.stdin), command line, etc.

05/08-reprojection.py
```"""Reproject a shapefile"""
import ogr
import osr
import os
import shutil

# Source and target file names
srcName = "NYC_MUSEUMS_LAMBERT.shp"
tgtName = "NYC_MUSEUMS_GEO.shp"

# Target spatial reference
tgt_spatRef = osr.SpatialReference()
tgt_spatRef.ImportFromEPSG(4326)

# Source shapefile
driver = ogr.GetDriverByName("ESRI Shapefile")
src = driver.Open(srcName, 0)
srcLyr = src.GetLayer()

# Source spatial reference
src_spatRef = srcLyr.GetSpatialRef()

# Target shapefile -
# there.
if os.path.exists(tgtName):
driver.DeleteDataSource(tgtName)
tgt = driver.CreateDataSource(tgtName)
lyrName = os.path.splitext(tgtName)[0]
tgtLyr = tgt.CreateLayer(lyrName, geom_type=ogr.wkbPoint)

# Layer definition
featDef = srcLyr.GetLayerDefn()

# Spatial Transform
trans = osr.CoordinateTransformation(src_spatRef, tgt_spatRef)

# Reproject and copy features
srcFeat = srcLyr.GetNextFeature()
while srcFeat:
geom = srcFeat.GetGeometryRef()
geom.Transform(trans)
feature = ogr.Feature(featDef)
feature.SetGeometry(geom)
tgtLyr.CreateFeature(feature)
feature.Destroy()
srcFeat.Destroy()
srcFeat = srcLyr.GetNextFeature()
src.Destroy()
tgt.Destroy()

# Create the prj file
tgt_spatRef.MorphToESRI()
prj = open(lyrName + ".prj", "w")
prj.write(tgt_spatRef.ExportToWkt())
prj.close()

# Just copy dbf contents over rather
# than rebuild the dbf using the
# ogr API since we're not changing
# anything.
srcDbf = os.path.splitext(srcName)[0] + ".dbf"
tgtDbf = lyrName + ".dbf"
shutil.copyfile(srcDbf, tgtDbf)```
```"""Open a shapefile for reading"""
import shapefile
print r
print r.bbox
print r.shapeType
print r.numRecords
print r.fields
print [item[0] for item in r.fields[1:]]
print r.record(2)
print r.record(2)[4]
fieldNames = [item[0] for item in r.fields[1:]]
name10 = fieldNames.index("NAME10")
print name10
print r.record(2)[name10]
zipRec = zip(fieldNames, rec)
print zipRec
for z in zipRec:
if z[0] == "NAME10":
print z[1]
for rec in enumerate(r.records()[:3]):
print rec[0]+1, ": ", rec[1]
counter = 0
for rec in r.iterRecords():
counter += 1
print counter
geom = r.shape(0)
print geom.points```
05/10-change-shp.py
```"""Convert a shapefile from lat/lon to UTM"""
import shapefile
import utm
w = shapefile.Writer(r.shapeType)
w.fields = list(r.fields)
w.records.extend(r.records())
for s in r.iterShapes():
lon,lat = s.points[0]
y,x,zone,band = utm.from_latlon(lat,lon)
w.point(x,y)
w.save("NYC_MUSEUMS_UTM")```
```"""Add a dbf field and column"""
import shapefile
w = shapefile.Writer(r.shapeType)
w.fields = list(r.fields)
w.records.extend(r.records())
w.field("LAT","F",8,5)
w.field("LON","F",8,5)
for i in range(geo.numRecords):
lon, lat = geo.shape(i).points[0]
w.records[i].extend([lat,lon])
w._shapes.extend(r.shapes())
w.save("NYC_MUSEUMS_UTM")```
05/12-merge.py
```"""Merge multiple shapefiles"""
import glob
import shapefile
files = glob.glob("footprints_*shp")
w = shapefile.Writer()
r = None
for f in files:
w._shapes.extend(r.shapes())
w.records.extend(r.records())
w.fields = list(r.fields)
w.save("Merged")```
05/13-split.py
```"""Split a shapefile"""
import shapefile
import utm
w = shapefile.Writer(r.shapeType)
w.fields = list(r.fields)
for sr in r.shapeRecords():
utmPoints = []
for p in sr.shape.points:
x,y,band,zone = utm.from_latlon(p[1],p[0])
utmPoints.append([x,y])
area = abs(shapefile.signed_area(utmPoints))
if  area <= 100:
w._shapes.append(sr.shape)
w.records.append(sr.record)
w.save("footprints_185")
# Verify changes
print r.numRecords
print subset.numRecords```
05/14-pip.py
```def point_in_poly(x,y,poly):
# check if point is a vertex
if (x,y) in poly: return True

# check if point is on a boundary
for i in range(len(poly)):
p1 = None
p2 = None
if i==0:
p1 = poly[0]
p2 = poly[1]
else:
p1 = poly[i-1]
p2 = poly[i]
if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
return True

n = len(poly)
inside = False

p1x,p1y = poly[0]
for i in range(n+1):
p2x,p2y = poly[i % n]
if y > min(p1y,p2y):
if y <= max(p1y,p2y):
if x <= max(p1x,p2x):
if p1y != p2y:
xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
if p1x == p2x or x <= xints:
inside = not inside
p1x,p1y = p2x,p2y

if inside: return True
return False

# Test a point inside the polygon for inclusion
myPolygon = [(-70.593016,-33.416032), (-70.589604,-33.415370),
(-70.589046,-33.417340), (-70.592351,-33.417949),
(-70.593016,-33.416032)]
lon = -70.592000
lat = -33.416000

print point_in_poly(lon, lat, myPolygon)

# test an vertex for inclusion
lon = -70.593016
lat = -33.416032
print point_in_poly(lon, lat, myPolygon)```
05/15-attr-sel.py
```'''
Attribute selection for subset
'''
import shapefile
# Create a writer instance
w = shapefile.Writer(r.shapeType)
# Copy the fields to the writer
w.fields = list(r.fields)
# Grab the geometry and records from all features
# with the correct population
selection = []
for rec in enumerate(r.records()):
if rec[1][15] < 5000:
selection.append(rec)
# Add the geometry and records to the writer
for rec in selection:
w._shapes.append(r.shape(rec[0]))
w.records.append(rec[1])
# Save the new shapefile
w.save("MS_Urban_Subset")```
05/16-dot-density.py
```"""Create a dot-density thematic map"""                                                import shapefile
import random
import pngcanvas

def point_in_poly(x,y,poly):
"""Boolean: is a point inside a polygon?"""
# check if point is a vertex
if (x,y) in poly: return True
# check if point is on a boundary
for i in range(len(poly)):
p1 = None
p2 = None
if i==0:
p1 = poly[0]
p2 = poly[1]
else:
p1 = poly[i-1]
p2 = poly[i]
if p1[1] == p2[1] and p1[1] == y and \
x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
return True
n = len(poly)
inside = False
p1x,p1y = poly[0]
for i in range(n+1):
p2x,p2y = poly[i % n]
if y > min(p1y,p2y):
if y <= max(p1y,p2y):
if x <= max(p1x,p2x):
if p1y != p2y:
xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
if p1x == p2x or x <= xints:
inside = not inside
p1x,p1y = p2x,p2y
if inside: return True
else: return False

def world2screen(bbox, w, h, x, y):
"""convert geospatial coordinates to pixels"""
minx,miny,maxx,maxy = bbox
xdist = maxx - minx
ydist = maxy - miny
xratio = w/xdist
yratio = h/ydist
px = int(w - ((maxx - x) * xratio))
py = int((maxy - y) * yratio)
return (px,py)

# Open the census shapefile
# Set the output image size
iwidth = 600
iheight = 400
# Get the index of the population field
pop_index = None
dots = []
for i,f in enumerate(inShp.fields):
if f[0] == "POPULAT11":
# Account for deletion flag
pop_index = i-1
# Calculate the density and plot points
for sr in inShp.shapeRecords():
population = sr.record[pop_index]
# Density ratio - 1 dot per 100 people
density = population / 100
found = 0
# Randomly distribute points until we
# have the correct density
while found < density:
minx, miny, maxx, maxy = sr.shape.bbox
x = random.uniform(minx,maxx)
y = random.uniform(miny,maxy)
if point_in_poly(x,y,sr.shape.points):
dots.append((x,y))
found += 1

# Set up the PNG output image
c = pngcanvas.PNGCanvas(iwidth,iheight)
# Draw the red dots
c.color = (255,0,0,0xff)
for d in dots:
x,y = world2screen(inShp.bbox, iwidth, iheight, *d)
c.filledRectangle(x-1,y-1,x+1,y+1)
# Draw the census tracts
c.color = (0,0,0,0xff)
for s in inShp.iterShapes():
pixels = []
for p in s.points:
pixel = world2screen(inShp.bbox, iwidth, iheight, *p)
pixels.append(pixel)
c.polyline(pixels)
# Save the image
img = open("DotDensity.png","wb")
img.write(c.dump())
img.close()```
05/17-choropleth.py
```"""Choropleth thematic map"""
import math
import shapefile
import Image
import ImageDraw

def world2screen(bbox, w, h, x, y):
"""convert geospatial coordinates to pixels"""
minx,miny,maxx,maxy = bbox
xdist = maxx - minx
ydist = maxy - miny
xratio = w/xdist
yratio = h/ydist
px = int(w - ((maxx - x) * xratio))
py = int((maxy - y) * yratio)
return (px,py)

# Open our shapefile
iwidth = 600
iheight = 400
# PIL Image
img = Image.new("RGB", (iwidth,iheight), (255,255,255))
# PIL Draw module for polygon fills
draw = ImageDraw.Draw(img)
# Get the population AND area index
pop_index = None
area_index = None
for i,f in enumerate(inShp.fields):
if f[0] == "POPULAT11":
# Account for deletion flag
pop_index = i-1
elif f[0] == "AREASQKM":
area_index = i-1
# Draw the polygons
for sr in inShp.shapeRecords():
density = sr.record[pop_index]/sr.record[area_index]
weight = min(math.sqrt(density/80.0), 1.0) * 50
R = int(205 - weight)
G = int(215 - weight)
B = int(245 - weight)
pixels = []
for x,y in sr.shape.points:
(px,py) = world2screen(inShp.bbox, iwidth, iheight, x, y)
pixels.append((px,py))
draw.polygon(pixels, outline=(255,255,255), fill=(R,G,B))
img.save("choropleth.png")```
05/18-excel.py
```"""Convert a spreadsheet to a shapefile"""
import xlrd
import shapefile

xls = xlrd.open_workbook("NYC_MUSEUMS_GEO.xls")
sheet = xls.sheet_by_index(0)

# Open the shapefile writer
w = shapefile.Writer(shapefile.POINT)

# Move data from spreadsheet to shapefile
for i in range(sheet.ncols):
w.field(str(sheet.cell(0,i).value), "C", 40)
for i in range(1, sheet.nrows):
values = []
for j in range(sheet.ncols):
values.append(sheet.cell(i,j).value)
w.record(*values)
w.point(float(values[-2]),float(values[-1]))
w.save("NYC_MUSEUMS_XLS2SHP")```
05/19-nmea.py
```"""Parse NMEA GPS strings"""
from pynmea.streamer import NMEAStream
nmeaFile = open("nmea.txt")
nmea_stream = NMEAStream(stream_obj=nmeaFile)
next_data = nmea_stream.get_objects()
nmea_objects = []
while next_data:
nmea_objects += next_data
next_data = nmea_stream.get_objects()
# The NMEA stream is parsed!
# Let's loop through the
# Python object types:
for nmea_ob in nmea_objects:
if hasattr(nmea_ob, "lat"):
print "Lat/Lon: (%s, %s)" % (nmea_ob.lat, nmea_ob.lon)```

## Playing with rasters

06/01-swap-bands.py
```"""Swap bands in a raster satellite image"""
# Module within the GDAL python package
import gdalnumeric
# name of our source image
src = "FalseColor.tif"
# load the source image into an array
# swap bands 1 and 2 for a natural color image.
# We will use numpy "advanced slicing" to reorder the bands.
# Using the source image
gdalnumeric.SaveArray(arr[[1,0,2],:], "swap.tif", format="GTiff", prototype=src)```
06/02-histogram.py
```"""Graph a histogram of a remotely sensed image"""
import gdalnumeric
import turtle as t

def histogram(a, bins=range(0,256)):
"""
Histogram function for multi-dimensional array.
a = array
bins = range of numbers to match
"""
fa = a.flat
n = gdalnumeric.numpy.searchsorted(gdalnumeric.numpy.sort(fa), bins)
n = gdalnumeric.numpy.concatenate([n, [len(fa)]])
hist = n[1:]-n[:-1]
return hist

def draw_histogram(hist, scale=True):
t.color("black")
# Draw the axes
axes = ((-355, -200),(355, -200),(-355, -200),(-355, 250))
t.up()
for p in axes:
t.goto(p)
t.down()
# Labels
t.up()
t.goto(0, -250)
t.write("VALUE",font=("Arial,",12,"bold"))
t.up()
t.goto(-400, 280)
t.write("FREQUENCY",font=("Arial,",12,"bold"))
# Tick marks
# x axis
x = -355
y = -200
t.up()
for i in range(1,11):
x = x+65
t.goto(x,y)
t.down()
t.goto(x,y-10)
t.up()
t.goto(x,y-25)
t.write("%s" % (i*25), align="center")
# y axis
x = -355
y = -200
t.up()
pixels = sum(hist[0])
if scale:
max = 0
for h in hist:
hmax = h.max()
if hmax > max:
max = hmax
pixels = max
label = pixels/10
for i in range(1,11):
y = y+45
t.goto(x,y)
t.down()
t.goto(x-10,y)
t.up()
t.goto(x-15,y-6)
t.write("%s" % (i*label), align="right")
# Plot each histogram as a colored line
x_ratio = 709.0 / 256
y_ratio = 450.0 / pixels
# Add more colors to this list if comparing
# more than 3 bands or 1 image
colors = ["red", "green", "blue"]
for j in range(len(hist)):
h = hist[j]
x = -354
y = -199
t.up()
t.goto(x,y)
t.down()
t.color(colors[j])
for i in range(256):
x = i * x_ratio
y = h[i] * y_ratio
x = x - (709/2)
y = y + -199
t.goto((x,y))

im = "swap.tif"
histograms = []
for b in arr:
histograms.append(histogram(b))
draw_histogram(histograms)

# Hide our pen
t.pen(shown=False)
t.done()```
06/03-stretch.py
```"""Perform a histogram stretch on an image"""
import gdalnumeric
import operator

def histogram(a, bins=range(0,256)):
"""
Histogram function for multi-dimensional array.
a = array
bins = range of numbers to match
"""
fa = a.flat
n = gdalnumeric.numpy.searchsorted(gdalnumeric.numpy.sort(fa), bins)
n = gdalnumeric.numpy.concatenate([n, [len(fa)]])
hist = n[1:]-n[:-1]
return hist

def stretch(a):
"""
Performs a histogram stretch on a gdalnumeric array image.
"""
hist = histogram(a)
lut = []
for b in range(0, len(hist), 256):
# step size
step = reduce(operator.add, hist[b:b+256]) / 255
# create equalization lookup table
n = 0
for i in range(256):
lut.append(n / step)
n = n + hist[i+b]
gdalnumeric.numpy.take(lut, a, out=a)
return a

src = "swap.tif"
stretched = stretch(arr)
gdalnumeric.SaveArray(arr, "stretched.tif", format="GTiff", prototype=src)```
06/04-clip.py
```"""Clip a raster image using a shapefile"""
import operator
import gdal, gdalnumeric, osr
import shapefile
import Image, ImageDraw

# Raster image to clip
raster = "stretched.tif"

# Polygon shapefile used to clip
shp = "hancock"

# Name of clipped raster file(s)
output = "clip"

def imageToArray(i):
"""
Converts a Python Imaging Library array to a gdalnumeric image.
"""
a=gdalnumeric.numpy.fromstring(i.tostring(),'b')
a.shape=i.im.size[1], i.im.size[0]
return a

def world2Pixel(geoMatrix, x, y):
"""
Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
the pixel location of a geospatial coordinate
"""
ulX = geoMatrix[0]
ulY = geoMatrix[3]
xDist = geoMatrix[1]
yDist = geoMatrix[5]
rtnX = geoMatrix[2]
rtnY = geoMatrix[4]
pixel = int((x - ulX) / xDist)
line = int((ulY - y) / xDist)
return (pixel, line)

# Load the source data as a gdalnumeric array

# Also load as a gdal image to get geotransform (world file) info
srcImage = gdal.Open(raster)
geoTrans = srcImage.GetGeoTransform()

# Use pyshp to open the shapefile

# Convert the layer extent to image pixel coordinates
minX, minY, maxX, maxY = r.bbox
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)

# Calculate the pixel size of the new image
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)

clip = srcArray[:, ulY:lrY, ulX:lrX]

# Create a new geomatrix for the image
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY

# Map points to pixels for drawing the county boundary
# on a blank 8-bit, black and white, mask image.
pixels = []
for p in r.shape(0).points:
pixels.append(world2Pixel(geoTrans, p[0], p[1]))
rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
# Create a blank image in PIL to draw the polygon.
rasterize = ImageDraw.Draw(rasterPoly)
rasterize.polygon(pixels, 0)
# Convert the PIL image to a NumPy array

# Clip the image using the mask

# Save ndvi as tiff
gdalnumeric.SaveArray(clip, "%s.tif" % output, format="GTiff", prototype=raster)```
06/05-classify.py
```"""Classify a remotely sensed image"""
import gdalnumeric

# Input file name (thermal image)
src = "thermal.tif"

# Output file name
tgt = "classified.jpg"

# Load the image into numpy using gdal

# Split the histogram into 20 bins as our classes
classes = gdalnumeric.numpy.histogram(srcArr, bins=20)[1]

# Color look-up table (LUT) - must be len(classes)+1.
# Specified as R,G,B tuples
lut = [[0,0,0],[191,48,48],[166,0,0],[255,64,64],
[255,115,115],[255,116,0],[191,113,48],[255,178,115],
[0,153,153],[29,115,115],[0,99,99],[166,75,0],
[0,204,0],[51,204,204],[255,150,64],[92,204,204],[38,153,38],[0,133,0],
[57,230,57],[103,230,103],[184,138,0]]
print len(lut)

# Starting value for classification
start = 1

# Set up the RGB color JPEG output image
rgb = gdalnumeric.numpy.zeros((3, srcArr.shape[0], srcArr.shape[1],), gdalnumeric.numpy.float32)

# Process all classes and assign colors
for i in range(len(classes)):
mask = gdalnumeric.numpy.logical_and(start <= srcArr, srcArr <= classes[i])
for j in range(len(lut[i])):
start = classes[i]+1

# Save the image
gdalnumeric.SaveArray(rgb.astype(gdalnumeric.numpy.uint8), tgt, format="JPEG")```
06/06-threshold.py
```"""Threshold an image to black and white"""
import gdalnumeric

# Input file name (thermal image)
src = "islands.tif"

# Output file name
tgt = "islands_classified.tiff"

# Load the image into numpy using gdal

# Split the histogram into 20 bins as our classes
classes = gdalnumeric.numpy.histogram(srcArr, bins=2)[1]

lut = [[255,0,0],[0,0,0],[255,255,255]]

# Starting value for classification
start = 1

# Set up the output image
rgb = gdalnumeric.numpy.zeros((3, srcArr.shape[0], srcArr.shape[1],), gdalnumeric.numpy.float32)

# Process all classes and assign colors
for i in range(len(classes)):
mask = gdalnumeric.numpy.logical_and(start <= srcArr, srcArr <= classes[i])
for j in range(len(lut[i])):
start = classes[i]+1

# Save the image
gdalnumeric.SaveArray(rgb.astype(gdalnumeric.numpy.uint8), tgt, format="GTIFF", prototype=src)```
06/07-extract.py
```"""Automatically extract features of a thresholded image to a shapefile"""
import gdal
import ogr, osr

# Thresholded input raster name
src = "islands_classified.tiff"
# Output shapefile name
tgt = "extract.shp"
# OGR layer name
tgtLayer = "extract"
# Open the input raster
srcDS = gdal.Open(src)
# Grab the first band
band = srcDS.GetRasterBand(1)
# Force gdal to use the band as a mask
# Set up the output shapefile
driver = ogr.GetDriverByName("ESRI Shapefile")
shp = driver.CreateDataSource(tgt)
# Copy the spatial reference
srs = osr.SpatialReference()
srs.ImportFromWkt(srcDS.GetProjectionRef())
layer = shp.CreateLayer(tgtLayer, srs=srs)
# Set up the dbf file
fd = ogr.FieldDefn( "DN", ogr.OFTInteger )
layer.CreateField(fd)
dst_field = 0
# Automatically extract features from an image!
extract = gdal.Polygonize(band, mask, layer, dst_field, [], None)```
06/08-draw.py
```"""Rasterize a shapefile and account for polygon holes"""
import shapefile
import pngcanvas
# Open the extracted islands
# Setup the world to pixels conversion
xdist = r.bbox[2] - r.bbox[0]
ydist = r.bbox[3] - r.bbox[1]
iwidth = 800
iheight = 600
xratio = iwidth/xdist
yratio = iheight/ydist
polygons = []
# Loop through all shapes
for shape in r.shapes():
# Loop through all parts to catch
# polygon holes!
for i in range(len(shape.parts)):
pixels=[]
pt = None
if i<len(shape.parts)-1:
pt = shape.points[shape.parts[i]:shape.parts[i+1]]
else:
pt = shape.points[shape.parts[i]:]
for x,y in pt:
px = int(iwidth - ((r.bbox[2] - x) * xratio))
py = int((r.bbox[3] - y) * yratio)
pixels.append([px,py])
polygons.append(pixels)
# Set up the output canvas
c = pngcanvas.PNGCanvas(iwidth,iheight)
# Loop through the polygons and draw them
for p in polygons:
c.polyline(p)
# Save the image
f = file("extract.png", "wb")
f.write(c.dump())
f.close()```
06/09-change.py
```"""Perform a simple difference image change detection on a
'before' and 'after' image."""
import gdal, gdalnumeric
import numpy as np

# "Before" image
im1 = "before.tif"
# "After" image
im2 = "after.tif"
# Load before and after into arrays
# Perform a simple array difference on the images
diff = ar2 - ar1
# Set up our classification scheme to try
# and isolate significant changes
classes = np.histogram(diff, bins=5)[1]
# The color black is repeated to mask insignificant changes
lut = [[0,0,0],[0,0,0],[0,0,0],[0,0,0],[0,255,0],[255,0,0]]
# Starting value for classification
start = 1
# Set up the output image
rgb = np.zeros((3, diff.shape[0], diff.shape[1],), np.int8)
# Process all classes and assign colors
for i in range(len(classes)):
mask = np.logical_and(start <= diff, diff <= classes[i])
for j in range(len(lut[i])):