User Tools

Site Tools


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 = []
 
# Add Denver
cities.append(["DENVER",[-104.98, 39.74], 634265])
# Add Boulder
cities.append(["BOULDER",[-105.27, 40.02], 98889])
# Add Durango
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")
      out.write(header)
    for record in fileDesc:
        for field in record:
            out.write("%s\t" % field)
        out.write("\n")
    out.close()
    pickle.dump(cached, open(pickleJar, "w"))
 
header = "POSITION\tFIELD\tSAMPLE\tTYPE\tBYTE_ORDER\n"
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 == "":
        cached = pickle.load(open(pickleJar, "r"))
    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)
# Read min-x,min-y,max-x,max-y
# in little endian format
print(struct.unpack("<d", f.read(8)))
print(struct.unpack("<d", f.read(8)))
print(struct.unpack("<d", f.read(8)))
print(struct.unpack("<d", f.read(8)))
# Read all 4 values at once
f.seek(36)
print(struct.unpack("<dddd", f.read(32)))

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
print(earthquakes.readline())
print(earthquakes.readline())
# 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.login()
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.write(zipShape.read(shpName))
shxFile.write(zipShape.read(shxName))
dbfFile.write(zipShape.read(dbfName))
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.write(zipShape.read(fileName))
  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.add("hancock.shp")
tar.add("hancock.shx")
tar.add("hancock.dbf")
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)
memoryshape = io.BytesIO(cloudshape.read())
zipshape = zipfile.ZipFile(memoryshape)
cloudshp = zipshape.read("hancock.shp")
# Access Python string as an array
# read shapefile boundingbox
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")
soup = BeautifulStoneSoup(gpx.read())
# 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))"
poly = shapely.wkt.loads(wktPoly)
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
print(json.loads(jsdata))
# Parse and then dump GeoJSON
pydata = json.loads(jsdata)
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
 
shp = shapefile.Reader("../files/point")
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
poly = wkt.loads(wktPoly)
# 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
srcArray = gdalnumeric.LoadFile("../files/SatImage.tif")
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
r = shapefile.Reader("../files/hancock.shp")
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
r = shapefile.Reader("../files/hancock.shp")
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
pdf.add_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 -
# delete if it's already
# 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)
05/09-read-shp.py
"""Open a shapefile for reading"""
import shapefile
r = shapefile.Reader("MSCities_Geo_Pts")
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
r = shapefile.Reader("NYC_MUSEUMS_GEO")
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")
05/11-add-field.py
"""Add a dbf field and column"""
import shapefile
r = shapefile.Reader("NYC_MUSEUMS_UTM")
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)
geo = shapefile.Reader("NYC_MUSEUMS_GEO")
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:
  r = shapefile.Reader(f)
  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
r = shapefile.Reader("footprints_se")
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
r = shapefile.Reader("footprints_se")
subset = shapefile.Reader("footprints_185")
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 reader instance
r = shapefile.Reader("MS_UrbanAnC10")
# 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
inShp = shapefile.Reader("GIS_CensusTract_poly")
# 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
inShp = shapefile.Reader("GIS_CensusTract_poly")
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
# Shade the census tracts
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
 
# Open the spreadsheet reader
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
arr = gdalnumeric.LoadFile(src)
# 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 = []
arr = gdalnumeric.LoadFile(im)
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"
arr = gdalnumeric.LoadFile(src)
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
srcArray = gdalnumeric.LoadFile(raster)
 
# 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
r = shapefile.Reader("%s.shp" % shp)
 
# 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
mask = imageToArray(rasterPoly)   
 
# Clip the image using the mask
clip = gdalnumeric.numpy.choose(mask, (clip, 0)).astype(gdalnumeric.numpy.uint8)
 
# 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
srcArr = gdalnumeric.LoadFile(src)
 
# 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])):
        rgb[j] = gdalnumeric.numpy.choose(mask, (rgb[j], lut[i][j]))
    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
srcArr = gdalnumeric.LoadFile(src)
 
# 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])):
        rgb[j] = gdalnumeric.numpy.choose(mask, (rgb[j], lut[i][j]))
    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
mask = band
# 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
r = shapefile.Reader("extract.shp")
# 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
ar1 = gdalnumeric.LoadFile(im1).astype(np.int8)                                 
ar2 = gdalnumeric.LoadFile(im2)[1].astype(np.int8)
# 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])):
        rgb[j] = np.choose(mask, (rgb[j], lut[i][j]))
    start = classes[i]+1
# Save the output image
gdalnumeric.SaveArray(rgb, "change.tif", format="GTiff", prototype=im2)

Some introductory readings

wiki/python/basicpython2.txt · Last modified: 2020/07/17 06:39 (external edit)