import geometryIO
import shapely
import shapely.geometry
import pysal
import pylab as pl
from shapely.geometry import Point, LineString, Polygon
def plot(shapelyGeometries):
'Plot shapelyGeometries'
figure = pl.figure(num=None, figsize=(4, 4), dpi=180)
axes = pl.axes()
axes.set_aspect('equal', 'datalim')
axes.xaxis.set_visible(False)
axes.yaxis.set_visible(False)
draw(shapelyGeometries)
def draw(gs):
'Draw shapelyGeometries'
# Handle single and multiple geometries
try:
gs = iter(gs)
except TypeError:
gs = [gs]
# For each shapelyGeometry,
for g in gs:
gType = g.geom_type
if gType.startswith('Multi') or gType == 'GeometryCollection':
draw(g.geoms)
else:
draw_(g)
def draw_(g):
'Draw a shapelyGeometry; thanks to Sean Gilles'
gType = g.geom_type
if gType == 'Point':
pl.plot(g.x, g.y, 'k,')
elif gType == 'LineString':
x, y = g.xy
pl.plot(x, y, 'b-')
elif gType == 'Polygon':
x, y = g.exterior.xy
pl.fill(x, y, color='#cccccc', aa=True)
pl.plot(x, y, color='#666666', aa=True, lw=1.0)
for hole in g.interiors:
x, y = hole.xy
pl.fill(x, y, color='#ffffff', aa=True)
pl.plot(x, y, color='#999999', aa=True, lw=1.0)
import geometryIO
import datetime
import itertools
import os
from osgeo import ogr
from shapely import geometry
polygonsPath = os.path.expandvars('$HOME/Downloads/polygons.shp.zip')
pointsPath = os.path.expandvars('$HOME/Downloads/points.shp.tar.gz')
geometryIO.save(
targetPath=polygonsPath,
sourceProj4=geometryIO.proj4LL,
shapelyGeometries=[
geometry.Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)]),
geometry.Polygon([(10, 0), (10, 10), (20, 10), (20, 0), (10, 0)]),
],
fieldPacks=[
('xxx', 11111, 44444.44, datetime.datetime(1980, 1, 1)),
('yyy', 22222, 88888.88, datetime.datetime(1990, 1, 1)),
],
fieldDefinitions=[
('Name', ogr.OFTString),
('Population', ogr.OFTInteger),
('GDP', ogr.OFTReal),
('Date', ogr.OFTDate),
],
driverName='ESRI Shapefile',
targetProj4=geometryIO.proj4SM)
proj4, shapelyGeometries, fieldPacks, fieldDefinitions = geometryIO.load(
sourcePath=polygonsPath,
targetProj4=geometryIO.proj4LL)
for shapelyGeometry, fPack in itertools.izip(shapelyGeometries, fieldPacks):
print
for fValue, (fName, fType) in itertools.izip(fPack, fieldDefinitions):
print '%s = %s' % (fName, fValue)
print shapelyGeometry
geometryIO.save_points(
# Save to a compressed shapefile
targetPath=pointsPath,
# Declare that source coordinates are in longitude and latitude
sourceProj4=geometryIO.proj4LL,
# Specify coordinates
coordinateTuples=[
(0, +1),
(+1, 0),
(0, -1),
(-1, 0),
])
print geometryIO.load_points(pointsPath)[1]
Transform between spatial references.
transformGeometry = geometryIO.get_transformGeometry(geometryIO.proj4LL, geometryIO.proj4SM)
print transformGeometry(LineString(((3, 4), (10, 50), (20, 25))))
transformPoint = geometryIO.get_transformPoint(geometryIO.proj4LL, geometryIO.proj4SM)
print transformPoint(0, 0)
Create a geometry from coordinates.
point1 = Point((6, 10))
lineString1 = LineString(((3, 4), (10, 50), (20, 25)))
polygon1 = Polygon(
((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell
[
((2, 2), (2, 3), (3, 3), (3, 2), (2, 2)), # Hole
((3, 3), (3, 4), (4, 4), (4, 3), (3, 3)), # Hole
]
)
plot(lineString1)
Create a geometry from a geometry.
point2 = Point((10, 6))
lineString2 = LineString(((3, 14), (10, 60), (20, 35)))
polygon2 = Polygon(
((1, 11), (5, 11), (5, 15), (1, 15), (1, 11))
)
from shapely.geometry import MultiPoint, MultiLineString, MultiPolygon
multiPoint = MultiPoint([point1, point2])
multiLineString = MultiLineString([lineString1, lineString2])
multiPolygon = MultiPolygon([polygon1, polygon2])
point1B = point1.buffer(1)
lineString1B = lineString1.buffer(1)
polygon1B = polygon1.buffer(1)
plot(lineString1B)
Create a geometry from WKT.
wkt = 'LINESTRING(3 4, 10 50, 20 25)'
import shapely.wkt
shapelyGeometry = shapely.wkt.loads(wkt)
Create a geometry from WKB.
wkb = '\x01\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x18@\x00\x00\x00\x00\x00\x00$@'
import shapely.wkb
shapelyGeometry = shapely.wkb.loads(wkb)
Create a geometry from a class.
class City(object):
def __init__(self, name, x, y):
self.name = name
self.x = x
self.y = y
def __repr__(self):
return '<City(name={})>'.format(self.name.encode('ascii', 'replace'))
@property
def coordinates(self):
return self.x, self.y
@property
def __geo_interface__(self):
return {'type': 'Point', 'coordinates': self.coordinates}
class Flight(object):
def __init__(self, city1, city2):
self.city1 = city1
self.city2 = city2
@property
def coordinates(self):
return self.city1.coordinates, self.city2.coordinates
@property
def __geo_interface__(self):
return {'type': 'LineString', 'coordinates': self.coordinates}
def __repr__(self):
return '<Flight({}, {})>'.format(self.city1, self.city2)
# Initialize
city1 = City('Kumasi', -1.62196, 6.68711)
city2 = City('Accra', -0.20092, 5.55856)
flight1 = Flight(city1, city2)
# Make shape adapter
shape = shapely.geometry.asShape(flight1)
print type(shape)
print 'length = {}'.format(shape.length)
# Make shape
shape = shapely.geometry.shape(flight1)
print type(shape)
print 'length = {}'.format(shape.length)
Make sure that the geometries you want to simplify are in a planar coordinate system such as spherical mercator. Strange distortions can appear if you try to simplify geometries that are in an angular coordinate system such as latitude and longitude. An added benefit of simplifying in spherical mercator is that the simplification tolerance has units in meters.
geometry = shapely.geometry.Point(0, 0).buffer(0.75)
figure = pylab.figure(num=None, figsize=(4, 4), dpi=180)
for figureIndex, tolerance in enumerate([0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.2, 0.5]):
axes = figure.add_subplot(3, 3, figureIndex + 1, frame_on=False)
axes.xaxis.set_visible(False)
axes.yaxis.set_visible(False)
draw(geometry.simplify(tolerance))
# Import system modules
import shapely.ops
# Merge line segments to make a polygon
lineSegments = [
((0, 0), (1, 1)), ((1, 1), (2, 0)), ((2, 0), (0, 0)),
((3, 0), (4, 1)), ((4, 1), (5, 0)), ((5, 0), (3, 0)),
]
for geometry in shapely.ops.polygonize(lineSegments):
print geometry.wkt
# Import system modules
import shapely.ops
# Merge line segment to make a line segment
lineSegments = [
((0, 0), (1, 1)),
((1, 1), (2, 2)),
((2, 2), (3, 3)),
]
geometry = shapely.ops.linemerge(lineSegments)
print geometry.wkt
print geometry.simplify(0.5).wkt
# Import system modules
import shapely.ops
# Merge two geometries
geometry = point1.union(lineString1)
# Merge multiple geometries efficiently
geometry = shapely.ops.cascaded_union([
point1,
lineString1,
polygon1,
])
# Initialize
templateStatement = 'filter(geometry.intersects, lineStrings)'
templateSetup = """
from shapely.geometry import LineString
from shapely.prepared import prep
import random
lineStrings = [LineString([(
random.choice(xrange(100)),
random.choice(xrange(100)),
), (
random.choice(xrange(100)),
random.choice(xrange(100)),
)]) for x in xrange(100)]
{}
"""
# Compare performance
import timeit
print 'With prepared statement: {}'.format(timeit.timeit(
stmt=templateStatement,
setup=templateSetup.format('geometry = prep(lineStrings[0])'),
number=100))
print 'Without prepared statement: {}'.format(timeit.timeit(
stmt=templateStatement,
setup=templateSetup.format('geometry = lineStrings[0]'),
number=100))
import numpy as np
from itertools import product, izip
from scipy.spatial import KDTree
pointsXY1 = np.array(list(product(xrange(3), xrange(3))))
# print pointsXY1
k1 = KDTree(pointsXY1)
pointsXY2 = np.array(list(product(np.arange(0.5, 3), np.arange(0.5, 3))))
# print pointsXY2
k2 = KDTree(pointsXY2)
# X, Y
from scipy.spatial import KDTree
distances, indices = k1.query((0, 0), k=2)
# print pointsXY1[indices]
indices = k1.query_ball_point((0, 0), 1)
# print pointsXY1[indices]
indexLists = k1.query_ball_tree(k2, 1)
# for point, indices in izip(pointsXY1, indexLists):
# print '*** %s *** ' % point
# print pointsXY2[indices]
indexPairs = k1.query_pairs(1)
# for index1, index2 in indexPairs:
# print pointsXY1[index1], pointsXY1[index2]
from urllib import urlopen, urlencode
from simplejson import loads
def geocode(address):
url = 'https://maps.googleapis.com/maps/api/geocode/json?' + urlencode([
('address', address),
('sensor', 'false'),
])
json = loads(urlopen(url).read())
result = json['results'][0]
formattedAddress = result['formatted_address']
location = result['geometry']['location']
longitude, latitude = location['lng'], location['lat']
return dict(address=formattedAddress, longitude=longitude, latitude=latitude)
def geocodeLL(address):
valueByKey = geocode(address)
return valueByKey['longitude'], valueByKey['latitude']
addresses = [
'Grand Central',
'Penn Station',
'Wall Street',
'Empire State Building',
'Lincoln Center',
]
pointsLL = np.array([geocodeLL(_) for _ in addresses])
# Longitude, Latitude
from pysal.cg.kdtree import Arc_KDTree
earthRadiusInMeters = 6378100
k3 = Arc_KDTree(pointsLL, radius=earthRadiusInMeters)
distances, indices = k3.query(geocodeLL('Times Square'), 2)
for distance, index in izip(distances, indices):
print '%s is %i meters away' % (addresses[index], distance)
Load the state of Georgia from the provided GADM dataset, transform it to the spherical mercator projection, save the geometry with attributes at three levels of simplification and plot the three simplifications to PNGs.
geometryPath = 'GADM-USA-Georgia.shp.zip'
# Type your solution here and press CTRL-ENTER