# Python & GeoComputation

Longzhu Shen Spatial Ecology Jun 2019

Code avalability at

wget https://github.com/selvaje/spatial-ecology-codes/blob/master/docs/source/PYTHON/02_Geo_Python.ipynb


## Configuration

Adjust the ‘ast_note_interactivity’ setting to print all versus last expression in Jupyter cell.

:

#from IPython.core.interactiveshell import InteractiveShell
#InteractiveShell.ast_node_interactivity = "last_expr" #"all"


## Numpy Intro

Numpy is a foundational package for data analysis in Python. It enables extremely efficient mathematical operations on large data sets. It also can be used to store and manipulate tabular data.

### Creating N-dimensional arrays using NumPy

• There are many ways to create N-dimensional arrays.

Create 2X3 array initialized to all zeroes & check the dimensions of the array

:

import numpy as np

:

a = np.zeros((2,3))
a

:

array([[0., 0., 0.],
[0., 0., 0.]])

:

a.shape

:

(2, 3)


Create array using arange function

:

a = np.arange(6)
a

:

array([0, 1, 2, 3, 4, 5])


Create array initialized by list of lists

:

a = np.array([[0,1,2],[3,4,5]])
a

:

array([[0, 1, 2],
[3, 4, 5]])


### Get values from N-dimensional array

:

a[0,0]

:

0


slicing

:

a[0,:]

:

array([0, 1, 2])

:

a[:,-1]

:

array([2, 5])

:

a[0:2,1:3]

:

array([[1, 2],
[4, 5]])


### Modifying N-dimensional arrays

:

a

:

array([[0, 1, 2],
[3, 4, 5]])

:

a[0,0] = 25.0
a

:

array([[25,  1,  2],
[ 3,  4,  5]])

:

a[0,:] = np.array([10,11,12])
a

:

array([[10, 11, 12],
[ 3,  4,  5]])

:

a[:,-1] = np.array([20,21])
a

:

array([[10, 11, 20],
[ 3,  4, 21]])

:

a[0:2,1:3] = np.array([[10,11],[20,21]])
a

:

array([[10, 10, 11],
[ 3, 20, 21]])


• “Describes how numpy treats arrays with different shapes during arithmetic operations”

• Subject to constraints, smaller array is ‘broadcast’ across larger array so they have compatible shapes

• Looping occurs in C instead of Python - fast!

:

a[0,:] = 10.0
a

:

array([[10, 10, 10],
[ 3,  4, 21]])

:

a[:,:] = np.array([30,31,32])
a

:

array([[30, 31, 32],
[30, 31, 32]])

:

tmp = np.array([40,41]).reshape(2,1)
tmp

:

array([,
])

:

a[:,:] = tmp
a

:

array([[40, 40, 40],
[41, 41, 41]])

:

a[0:2,1:3] = 100.0
a

:

array([[ 40, 100, 100],
[ 41, 100, 100]])


### Math Operations on Arrays

• Arithmetic of arrays vs. scalars applies scalar to all elements of array

• Arithmetic of array vs. array operates elementwise

:

a = np.arange(4)
a

:

array([0, 1, 2, 3])

:

a + 0.5

:

array([0.5, 1.5, 2.5, 3.5])

:

a + a

:

array([0, 2, 4, 6])

:

a * a

:

array([0, 1, 4, 9])

:

np.dot(a, a)

:

14

:

np.dot(a.reshape(4,1), a.reshape(1,4))

:

array([[0, 0, 0, 0],
[0, 1, 2, 3],
[0, 2, 4, 6],
[0, 3, 6, 9]])


### Working with Tabular Data

:

solar_system_file = 'solar_system_abbr.csv'
dtype=['S10', 'S20', 'i4', 'f4' , 'f4'])

:

solar_system_data

:

array([(b'Sun', b'Star',  0, 1.3920e+06, 3.33e+05),
(b'Mercury', b'Terrestrial planet',  1, 4.8780e+03, 5.50e-02),
(b'Venus', b'Terrestrial planet',  2, 1.2104e+04, 8.15e-01),
(b'Earth', b'Terrestrial planet',  3, 1.2756e+04, 1.00e+00),
(b'Mars', b'Terrestrial planet',  4, 6.7870e+03, 1.07e-01),
(b'Jupiter', b'Gas giant',  6, 1.4280e+05, 3.18e+02),
(b'Saturn', b'Gas giant',  7, 1.2000e+05, 9.50e+01),
(b'Uranus', b'Gas giant',  8, 5.1118e+04, 1.50e+01),
(b'Neptune', b'Gas giant',  9, 4.9528e+04, 1.70e+01),
(b'Ceres', b'Dwarf planet',  5, 9.7460e+02, 1.60e-04),
(b'Pluto', b'Dwarf planet', 10, 2.3000e+03, 2.00e-03),
(b'Haumea', b'Dwarf planet', 11, 1.3000e+03, 7.00e-04),
(b'Makemake', b'Dwarf planet', 12, 1.4200e+03, 6.70e-04),
(b'Eris', b'Dwarf planet', 13, 2.3260e+03, 2.80e-03)],
dtype=[('f0', 'S10'), ('f1', 'S20'), ('f2', '<i4'), ('f3', '<f4'), ('f4', '<f4')])

:

solar_system_data.dtype.names = ('planet', 'type', 'order', 'diameter', 'mass')
solar_system_data

:

array([(b'Sun', b'Star',  0, 1.3920e+06, 3.33e+05),
(b'Mercury', b'Terrestrial planet',  1, 4.8780e+03, 5.50e-02),
(b'Venus', b'Terrestrial planet',  2, 1.2104e+04, 8.15e-01),
(b'Earth', b'Terrestrial planet',  3, 1.2756e+04, 1.00e+00),
(b'Mars', b'Terrestrial planet',  4, 6.7870e+03, 1.07e-01),
(b'Jupiter', b'Gas giant',  6, 1.4280e+05, 3.18e+02),
(b'Saturn', b'Gas giant',  7, 1.2000e+05, 9.50e+01),
(b'Uranus', b'Gas giant',  8, 5.1118e+04, 1.50e+01),
(b'Neptune', b'Gas giant',  9, 4.9528e+04, 1.70e+01),
(b'Ceres', b'Dwarf planet',  5, 9.7460e+02, 1.60e-04),
(b'Pluto', b'Dwarf planet', 10, 2.3000e+03, 2.00e-03),
(b'Haumea', b'Dwarf planet', 11, 1.3000e+03, 7.00e-04),
(b'Makemake', b'Dwarf planet', 12, 1.4200e+03, 6.70e-04),
(b'Eris', b'Dwarf planet', 13, 2.3260e+03, 2.80e-03)],
dtype=[('planet', 'S10'), ('type', 'S20'), ('order', '<i4'), ('diameter', '<f4'), ('mass', '<f4')])

:

solar_system_data['mass'].mean()

:

23817.64

:

import math
volume_of_earth = 4/3*math.pi*(solar_system_data['diameter']/2.)**3.
volume_of_earth

:

1086781292542.8892


## Interfacing with Network

:

# Retrieve a file
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)

:

('hancock.zip', <http.client.HTTPMessage at 0x152e96691390>)

:

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

b'time,latitude,longitude,depth,mag,magType,nst,gap,dmin,rms,net,id,updated,place,type,horizontalError,depthError,magError,magNst,status,locationSource,magSource\n'
b'2019-05-31T17:17:57.190Z,34.0203323,-117.6094971,19.12,1.96,ml,13,134,0.1116,0.28,ci,ci38613536,2019-05-31T17:20:00.408Z,"4km SSE of Ontario, CA",earthquake,2.14,1.35,0.412,2,automatic,ci,ci\n'
b'2019-05-31T17:05:15.100Z,33.4946667,-116.795,1.59,0.49,ml,19,57,0.01821,0.23,ci,ci38613504,2019-05-31T17:08:56.963Z,"9km NE of Aguanga, CA",earthquake,0.45,0.67,0.127,10,automatic,ci,ci\n'
b'2019-05-31T16:58:50.720Z,33.4896667,-116.7948333,2.37,1.59,ml,56,36,0.02311,0.2,ci,ci38613480,2019-05-31T17:10:18.350Z,"8km NE of Aguanga, CA",earthquake,0.19,0.34,0.151,27,automatic,ci,ci\n'
b'2019-05-31T16:39:45.470Z,45.867,-111.361,4.94,1.57,ml,13,88,0.302,0.13,mb,mb80346349,2019-05-31T17:05:26.100Z,"2km WNW of Manhattan, Montana",earthquake,0.39,8.54,0.255,6,reviewed,mb,mb\n'
b'2019-05-31T16:37:53.690Z,33.5103333,-116.5108333,11.17,0.7,ml,25,82,0.04997,0.18,ci,ci38613456,2019-05-31T16:41:30.754Z,"16km ESE of Anza, CA",earthquake,0.3,0.53,0.161,16,automatic,ci,ci\n'
b'2019-05-31T16:36:34.830Z,34.0466667,-117.5038333,2.53,1.07,ml,31,64,0.06375,0.2,ci,ci38613448,2019-05-31T16:40:14.360Z,"4km NNW of Glen Avon, CA",earthquake,0.32,0.47,0.168,30,automatic,ci,ci\n'
b'2019-05-31T16:22:07.250Z,38.8214989,-122.8526688,2.31,1.13,md,20,64,0.006873,0.03,nc,nc73190601,2019-05-31T17:06:03.673Z,"10km WNW of The Geysers, CA",earthquake,0.23,0.56,0.03,4,automatic,nc,nc\n'


## Zipping Files

:

# Unzip a shapefile
import zipfile
zipShape = zipfile.ZipFile(zip)
for fileName in zipShape.namelist():
out = open(fileName, "wb")
out.close()

:

# Add a shapefile to a tar archive
import tarfile
tar = tarfile.open("hancock.tar.gz", "w:gz")
tar.close()

tar

:

<tarfile.TarFile at 0x152e9608d3c8>

:

# Extract a shapefile from a gzipped tar archive
import tarfile
tar = tarfile.open("hancock.tar.gz", "r:gz")
tar.extractall()
tar.close()

tar

:

<tarfile.TarFile at 0x152e9608dac8>


## Raster Images

:

# Open a raster with gal
import gdal

raster = gdal.Open("SatImage.tif")

print (raster.RasterCount)
print (raster.RasterXSize)
print (raster.RasterYSize)

3
2592
2693

:

# Numpy/gdalnumeric - Read an image, extract a band, save a new image
from osgeo import gdalnumeric

band1 = srcArray
gdalnumeric.SaveArray(band1, "band1.jpg", format="JPEG")

:

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x7ff5036390f0> >


### Rasterization

:

# Rasterize a shapefile with PIL
import shapefile
from PIL import Image, ImageDraw

xdist = r.bbox - r.bbox
ydist = r.bbox - r.bbox
iwidth = 400
iheight = 600
xratio = iwidth/xdist
yratio = iheight/ydist
pixels = []
for x,y in r.shapes().points:
px = int(iwidth - ((r.bbox - x) * xratio))
py = int((r.bbox - 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")

:

# 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')


:

''


### Raster Processing

:

"""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)


:

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x7ff501aacb10> >

:

"""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)

# 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, diff.shape,), 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])):
start = classes[i]+1

# Save the output image
gdalnumeric.SaveArray(rgb, "change.tif", format="GTiff", prototype=im2)

:

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x7ff503644a50> >

[ ]:




## Shapefiles

:

#Examine a shapefile with ogr

from osgeo import ogr

# open the shapefile
shp = ogr.Open("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"))

1.0 1.0 First
3.0 1.0 Second
4.0 3.0 Third
2.0 2.0 Fourth
0.0 0.0 Appended

:

# Examine a shapefile with pyshp
import shapefile

for feature in shp.shapeRecords():
point = feature.shape.points
rec = feature.record
print (point, point, rec)

1.0 1.0 First
3.0 1.0 Second
4.0 3.0 Third
2.0 2.0 Fourth
0.0 0.0 Appended

:

"""Choropleth thematic map"""
import math
import shapefile
from PIL import Image, ImageDraw

'''convert geospatial coordinates to pixels'''
def world2screen(bbox, w, h, x, y):
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 == "POPULAT11":
# Account for deletion flag
pop_index = i-1
elif f == "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))

# Save the image to file
img.save("choropleth.png")

:

from osgeo import ogr

shapefile_name = "tl_2012_us_cbsa.shp"
shapefile = ogr.Open(shapefile_name)
numLayers = shapefile.GetLayerCount()

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

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\n" % (layerNum, numFeatures))

for featureNum in range(10):
feature = layer.GetFeature(featureNum)
featureName = feature.GetField("NAME")
print ("Feature %d has name %s" % (featureNum, featureName))


Shapefile contains 1 layers

Layer 0 has spatial reference +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
Layer 0 has 955 features

Feature 0 has name Hilo, HI
Feature 1 has name Bozeman, MT
Feature 2 has name Magnolia, AR
Feature 3 has name Kingston, NY
Feature 4 has name Astoria, OR
Feature 5 has name Kennewick-Pasco-Richland, WA
Feature 6 has name Fayetteville-Springdale-Rogers, AR-MO
Feature 7 has name Seattle-Tacoma-Bellevue, WA
Feature 8 has name Ellensburg, WA
Feature 9 has name Pierre, SD

:

from osgeo import ogr

shapefile_name = "tl_2012_us_cbsa.shp"
shapefile = ogr.Open(shapefile_name)
layer = shapefile.GetLayer(0)
feature = layer.GetFeature(2)

print ("Feature 2 has the following attributes:")

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)

Feature 2 has the following attributes:
CSAFP = None
CBSAFP = 31620
GEOID = 31620
NAME = Magnolia, AR
NAMELSAD = Magnolia, AR Micro Area
MEMI = 2
MTFCC = G3110
ALAND = 1984070931
AWATER = 1809509
INTPTLAT = +33.2230377
INTPTLON = -093.2328433

Feature's geometry data consists of a POLYGON

:

s = []
s.append("  " * 0)

:

s = []
s.append("  " * 0)
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))

POLYGON containing:

:

from osgeo import ogr

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)

shapefile_name = "tl_2012_us_cbsa.shp"
shapefile = osgeo.ogr.Open(shapefile_name)
layer = shapefile.GetLayer(0)
feature = layer.GetFeature(55)
geometry = feature.GetGeometryRef()

analyzeGeometry(geometry)

POLYGON containing:
LINEARRING with 9193 data points

:

# Examine a shapefile with pyshp
import shapefile

for feature in shp.shapeRecords():
point = feature.shape.points
rec = feature.record
print (point, point, rec)

1.0 1.0 First
3.0 1.0 Second
4.0 3.0 Third
2.0 2.0 Fourth
0.0 0.0 Appended

:

from osgeo import ogr

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

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

shapefile_name = "tl_2012_us_cbsa.shp"
shapefile = 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'])


Northernmost point is (-93.9197, 33.0195)
Southernmost point is (-93.8752, 31.8442)


## Sources : this material was adopted from the following sources

[ ]: