from kartograph import Kartograph
+from map import projections
+
+__all__ = ['Kartograph', 'projections']
+
+computes a circle cartogram for a given svg map + data file
+class Cartogram: def generate(self, svg_src, attr, csv_src, key, value):
+ regions = self.load_regions_from_svg(svg_src, attr)
+ data = self.load_csv(csv_src, key, value)
+ circles = []
+ for id in regions:
+ cx, cy = regions[id]
+ val = data[id]
+ circles.append(Circle(cx, cy, id, val))
+
+ self.attr = attr
+ self.key = value
+ self.circles = circles
+ self.compute_radii()
+ self.layout(700)
+ self.rescale()
+ self.correct()
+ self.layout(200, True)
+ self.rescale()
+ self.correct()
+ self.layout(100, False)
+ self.rescale()
+ self.correct()
+ self.to_svg() def load_regions_from_svg(self, url, attr):
+ import svg as svgdoc
+ svg = svgdoc.Document.load(url)
+ self.svg = svg
+ g = svg.doc.getElementsByTagName('g')[0]
+ coords = {}
+ for path in g.getElementsByTagName('path'):
+ path_str = path.getAttribte('d')
+ id = path.getAttribte('data-' + attr)
+ poly = restore_poly_from_path_str(path_str)
+ coords[id] = poly.center()
+ return coords def load_csv(self, url, key='id', value='val'):
+ import csv
+ doc = csv.reader(open(url), dialect='excel-tab')
+ head = None
+ data = {}
+ for row in doc:
+ if not head:
+ head = row
+ print head
+ else:
+ id = row[head.index(key)].strip()
+ val = float(row[head.index(value)])
+ data[id] = val
+ return data def compute_radii(self):
+ import sys, math
+ minv = 0
+ maxv = sys.maxint * -1
+ for c in self.circles:
+ minv = min(minv, c.value)
+ maxv = max(maxv, c.value)
+
+ for c in self.circles:
+ c.r = math.pow((c.value - minv) / (maxv - minv), 0.50) * 60
+ c.weight = c.value / maxv def layout(self, steps=100, correct=False):
+ for i in range(steps):if i % 100 == 0: + self.toSVG()
+ self.layout_step(correct) def layout_step(self, correct=False):
+ import math
+ pad = 0
+
+ if correct:
+ for C in self.circles:
+ v = Vector(C.ox - C.x, C.oy - C.y)
+ v.normalize()
+ v.resize(0.5)
+ C._move(v.x, v.y)
+
+ for A in self.circles:
+ for B in self.circles:
+ if A != B:
+ radsq = (A.r + B.r) * (A.r + B.r)
+ d = A.sqdist(B)
+ if radsq + pad > d:move circles away from each other
+ v = Vector(B.x - A.x, B.y - A.y)
+ v.normalize()
+ m = (math.sqrt(radsq) - math.sqrt(d)) * 0.25
+ v.resize(m)
+ A._move(v.x * -1 * B.weight, v.y * -1 * B.weight)
+ B._move(v.x * A.weight, v.y * A.weight)
+
+ for C in self.circles:
+ C.move() def rescale(self):
+ from geometry import BBox, View
+ svg = self.svg
+ svg_view = svg[1][0][0]
+ vh = float(svg_view['h'])
+ vw = float(svg_view['w'])
+
+ bbox = BBox()
+ for c in self.circles:
+ r = c.r
+ bbox.update((c.x + r, c.y + r))
+ bbox.update((c.x + r, c.y - r))
+ bbox.update((c.x - r, c.y + r))
+ bbox.update((c.x - r, c.y - r))
+
+ view = View(bbox, vw, vh)
+ for c in self.circles:
+ c.r *= view.scale
+ x, y = view.project((c.x, c.y))
+ c.x = x
+ c.y = y def correct(self):
+ for A in self.circles:
+ intersects = False
+ for B in self.circles:
+ if A != B:
+ radsq = (A.r + B.r) * (A.r + B.r)
+ d = A.sqdist_o(B)
+ if radsq > d:
+ intersects = True
+ break
+ if not intersects:
+ A.x = A.ox
+ A.y = A.oy def to_svg(self):
+ svg = self.svg
+
+ g = svg.node('g', svg.root, id="cartogram", fill="red", fill_opacity="0.5")
+
+ for circle in self.circles:
+ c = svg.node('circle', g, cx=circle.x, cy=circle.y, r=circle.r)
+ c.setAttribute('data-' + self.attr, circle.id)
+ c.setAttribute('data-' + self.key.lower(), circle.value)
+ g.append(c)
+
+ svg.preview()svg.save('cartogram.svg')
+class Circle: def __init__(self, x, y, id, value):
+ self.x = self.ox = float(x)
+ self.y = self.oy = float(y)
+ self.id = id
+ self.value = float(value)
+ self.dx = 0
+ self.dy = 0 def _move(self, x, y):
+ self.dx += x
+ self.dy += y def move(self):
+ self.x += self.dx
+ self.y += self.dy
+ self.dx = 0
+ self.dy = 0 def __repr__(self):
+ return '<Circle x=%.1f, y=%.1f, id=%s, val=%f >' % (self.x, self.y, self.id, self.value) def sqdist(self, circ):
+ dx = self.x - circ.x
+ dy = self.y - circ.y
+ return dx * dx + dy * dybeen too lazy to code this myself, instead I took code from here +http://www.kokkugia.com/wiki/index.php5?title=Python_vector_class
+ def sqdist_o(self, circ):
+ dx = self.ox - circ.x
+ dy = self.oy - circ.y
+ return dx * dx + dy * dyclass Vector:Class properties
+ def __init__(self, x, y):
+ self.x = float(x)
+ self.y = float(y)represent as a string
+ def __repr__(self):
+ return 'Vector(%s, %s)' % (self.x, self.y)
+
+ '''
+ Class Methods / Behaviours
+ ''' def zero(self):
+ self.x = 0.0
+ self.y = 0.0
+ return self def clone(self):
+ return Vector(self.x, self.y) def normalize(self):
+ from math import sqrt
+ if self.x == 0 and self.y == 0:
+ return self
+ norm = float(1.0 / sqrt(self.x * self.x + self.y * self.y))
+ self.x *= norm
+ self.y *= normself.z *= norm
+ return self def invert(self):
+ self.x = -(self.x)
+ self.y = -(self.y)
+ return self def resize(self, sizeFactor):
+ self.normalize
+ self.scale(sizeFactor)
+ return self def minus(self, t):
+ self.x -= t.x
+ self.y -= t.yself.z -= t.z
+ return self def plus(self, t):
+ self.x += t.x
+ self.y += t.yself.z += t.z
+ return self def roundToInt(self):
+ self.x = int(self.x)
+ self.y = int(self.y)
+ return selfReturns the squared length of this vector.
+ def lengthSquared(self):
+ return float((self.x * self.x) + (self.y * self.y))Returns the length of this vector.
+ def length(self):
+ from math import sqrt
+ return float(sqrt(self.x * self.x + self.y * self.y))Computes the dot product of this vector and vector v2
+ def dot(self, v2):
+ return (self.x * v2.x + self.y * v2.y)Linearly interpolates between vectors v1 and v2 and returns the result point = (1-alpha)v1 + alphav2.
+ def interpolate(self, v2, alpha):
+ self.x = float((1 - alpha) * self.x + alpha * v2.x)
+ self.y = float((1 - alpha) * self.y + alpha * v2.y)
+ return Vector(self.x, self.y)Returns the angle in radians between this vector and the vector parameter; +the return value is constrained to the range [0,PI].
+ def angle(self, v2):
+ from math import acos
+ vDot = self.dot(v2) / (self.length() * v2.length())
+ if vDot < -1.0:
+ vDot = -1.0
+ if vDot > 1.0:
+ vDot = 1.0
+ return float(acos(vDot))Limits this vector to a given size. +NODEBOX USERS: name should change as 'size' and 'scale' are reserved words in Nodebox!
+ def limit(self, size):
+ if (self.length() > size):
+ self.normalize()
+ self.scale(size)Point Methods +Returns the square of the distance between this tuple and tuple t1.
+ def distanceSquared(self, t1):
+ dx = self.x - t1.x
+ dy = self.y - t1.y
+ return (dx * dx + dy * dy)NODEBOX USERS: name should change as 'scale' is reserved word in Nodebox!
+ def scale(self, s):
+ self.x *= s
+ self.y *= s
+ return selfNODEBOX USERS: name should change as 'translate' is reserved word in Nodebox!
+ def translate(self, vec):
+ self.plus(vec) def distance(self, pt):
+ from math import sqrt
+ dx = self.x - pt.x
+ dy = self.y - pt.y
+ return float(sqrt(dx * dx + dy * dy))restores a list of polygons from a SVG path string
+def restore_poly_from_path_str(path_str): contours = path_str.split('Z') # last contour may be empty
+ from Polygon import Polygon as Poly
+ poly = Poly()
+ for c_str in contours:
+ if c_str.strip() != "":
+ pts_str = c_str.strip()[1:].split("L")
+ pts = []
+ for pt_str in pts_str:
+ x, y = map(float, pt_str.split(','))
+ pts.append((x, y))
+ poly.addContour(pts, is_clockwise(pts))
+ return polyreturns true if a given polygon is in clockwise order
+def is_clockwise(pts): s = 0
+ for i in range(len(pts) - 1):
+ if 'x' in pts[i]:
+ x1 = pts[i].x
+ y1 = pts[i].y
+ x2 = pts[i + 1].x
+ y2 = pts[i + 1].y
+ else:
+ x1, y1 = pts[i]
+ x2, y2 = pts[i + 1]
+ s += (x2 - x1) * (y2 + y1)
+ return s >= 0
+
+command line interface for kartograph
+import argparse
+import os.path
+import json
+from errors import KartographErrorclass bcolors:
+ HEADER = '\033[95m'
+ OKBLUE = '\033[94m'
+ OKGREEN = '\033[92m'
+ WARNING = '\033[93m'
+ FAIL = '\033[91m'
+ ENDC = '\033[0m' def disable(self):
+ self.HEADER = ''
+ self.OKBLUE = ''
+ self.OKGREEN = ''
+ self.WARNING = ''
+ self.FAIL = ''
+ self.ENDC = ''
+
+parser = argparse.ArgumentParser(prog='kartograph', description='generating svg maps from shapefiles')subparsers = parser.add_subparsers(help='sub-command help')
+parser_svg = subparsers.add_parser('svg', help='generates svg map')
+parser.add_argument('config', type=argparse.FileType('r'), help='the configuration for the map. accepts json and yaml.')
+parser.add_argument('--output', '-o', metavar='FILE', type=argparse.FileType('w'), help='the file in which the map will be stored')
+parser.add_argument('--verbose', '-v', nargs='?', metavar='', const=True, help='verbose mode')
+parser.add_argument('--format', '-f', metavar='svg|kml', help='output format, if not specified it will be guessed from output filename or default to svg')
+parser.add_argument('--preview', '-p', nargs='?', metavar='', const=True, help='opens the generated svg for preview')
+
+from kartograph import Kartograph
+import time
+import sys
+import osdef parse_config(f):
+ content = f.read()
+ if f.name[-5:].lower() == '.json':
+ try:
+ cfg = json.loads(content)
+ except Exception, e:
+ raise KartographError('parsing of json map configuration failed.\n' + e)
+ else:
+ return cfg
+ elif f.name[-5:].lower() == '.yaml':
+ import yaml
+ try:
+ cfg = yaml.load(content)
+ except Exception, e:
+ raise KartographError('parsing of yaml map configuration failed.\n' + e)
+ else:
+ return cfg
+ else:
+ raise KartographError('supported config formats are .json and .yaml')def render_map(args):
+ cfg = parse_config(args.config)
+ K = Kartograph()
+ if args.format:
+ format = args.format
+ elif args.output:
+ format = os.path.splitext(args.output.name)[1][1:]
+ else:
+ format = 'svg'
+ try:generate the map
+ K.generate(cfg, args.output, verbose=args.verbose, preview=args.preview, format=format)
+
+ except Exception, e:
+ import traceback
+ ignore_path_len = len(__file__) - 7
+ exc = sys.exc_info()
+ for (filename, line, func, code) in traceback.extract_tb(exc[2]):
+ if filename[:len(__file__) - 7] == __file__[:-7]:
+ print ' \033[1;33;40m%s\033[0m, \033[0;37;40min\033[0m %s()\n \033[1;31;40m%d:\033[0m \033[0;37;40m%s\033[0m' % (filename[ignore_path_len:], func, line, code)
+ else:
+ print ' %s, in %s()\n %d: %s' % (filename, func, line, code)
+ print
+ print e
+ exit(-1)
+
+parser.set_defaults(func=render_map)def main():
+ start = time.time()
+
+ try:
+ args = parser.parse_args()
+ except IOError, e:
+ parser.print_help()
+ print '\nIOError:', e
+ except Exception, e:
+ parser.print_help()
+ print '\nError:', e
+ else:
+ args.func(args)
+
+ elapsed = (time.time() - start)
+ print 'execution time: %.3f secs' % elapsed
+ sys.exit(0)
+
+
+if __name__ == "__main__":
+ main()
+
+error classes for kartograph
+Base class for exceptions in this module.
+class KartographError(Exception): def __str__(self):
+ return '\033[0;31;40mKartograph-Error:\033[0m ' + super(KartographError, self).__str__()class KartographOptionParseError(KartographError):
+ passclass KartographShapefileAttributesError(KartographError):
+ passclass KartographLayerSourceError(KartographError):
+ pass
+
+layer filter
+import redef filter_record(filt, record):
+ if isinstance(filt, dict):
+ if 'and' in filt:
+ res = True
+ for sfilt in filt['and']:
+ res = res and filter_record(sfilt, record)
+ elif 'or' in filt:
+ res = False
+ for sfilt in filt['or']:
+ res = res or filter_record(sfilt, record)
+ elif isinstance(filt, list):
+ res = filter_single(filt, record)
+ return resdef filter_single(filt, record):
+ key, comp, val = filt
+ prop = record[key]
+ comp = comp.lower().split(' ')
+
+ if 'in' in comp:
+ res = prop in val
+ elif 'like' in comp:
+ res = re.search('^' + _escape_regex(val).replace('%', '.*') + '$', prop) is not None
+ elif 'matches' in comp:
+ res = re.search(val, prop) is not None
+ elif 'is' in comp or '=' in comp:
+ res = prop == val
+ elif 'greater' in comp or ('>' in comp):
+ res = prop > val
+ elif 'less' in comp or '<' in comp:
+ res = prop < val
+ if 'not' in comp:
+ return not res
+ else:
+ return resdef _escape_regex(s):
+ chars = ('.', '*', '?', '+', '(', ')', '[', ']', '-')
+ for c in chars:
+ s = s.replace(c, '\\' + c)
+ return s
+
+geometry package
+__all__ = ['Feature', 'Geometry', 'SolidGeometry', 'MultiPolygon', 'BBox', 'Point', 'View', 'Line', 'PolyLine', 'create_feature']
+
+from feature import *
+from geometry import Geometry, SolidGeometry
+from polygon import MultiPolygon
+from point import Point
+from bbox import BBox
+from view import View
+from line import Line, PolyLine
+
+from point import Point2D bounding box
+class BBox(object): def __init__(self, width=None, height=None, left=0, top=0):
+ import sys
+ if width == None:
+ self.xmin = sys.maxint
+ self.xmax = sys.maxint * -1
+ else:
+ self.xmin = self.left = left
+ self.xmax = self.right = left + width
+ self.width = width
+ if height == None:
+ self.ymin = sys.maxint
+ self.ymax = sys.maxint * -1
+ else:
+ self.ymin = self.top = top
+ self.ymax = self.bottom = height + top
+ self.height = height def update(self, pt):
+ if not isinstance(pt, Point):
+ pt = Point(pt[0], pt[1])
+ self.xmin = min(self.xmin, pt.x)
+ self.ymin = min(self.ymin, pt.y)
+ self.xmax = max(self.xmax, pt.x)
+ self.ymax = max(self.ymax, pt.y)
+
+ self.left = self.xmin
+ self.top = self.ymin
+ self.right = self.xmax
+ self.bottom = self.ymax
+ self.width = self.xmax - self.xmin
+ self.height = self.ymax - self.yminreturns true if two bounding boxes overlap
+ def intersects(self, bbox): return bbox.left < self.right and bbox.right > self.left and bbox.top < self.bottom and bbox.bottom > self.topcheck if a point is inside the bbox
+ def check_point(self, pt): return pt[0] > self.xmin and pt[0] < self.xmax and pt[1] > self.ymin and pt[1] < self.ymax def __str__(self):
+ return 'BBox(x=%.2f, y=%.2f, w=%.2f, h=%.2f)' % (self.left, self.top, self.width, self.height) def join(self, bbox):
+ self.update(Point(bbox.left, bbox.top))
+ self.update(Point(bbox.right, bbox.bottom)) def inflate(self, amount=0, inflate=False):
+ if inflate:
+ d = min(self.width, self.height)
+ amount += d * inflate
+ self.xmin -= amount
+ self.ymin -= amount
+ self.xmax += amount
+ self.ymax += amount
+
+ self.left = self.xmin
+ self.top = self.ymin
+ self.right = self.xmax
+ self.bottom = self.ymax
+ self.width = self.xmax - self.xmin
+ self.height = self.ymax - self.ymin
+
+DEPRECATED
+base class for all geometry
+class Geometry:project geometry
+ def project(self, proj): raise NotImplementedError('project() is not implemented') def bbox(self):
+ raise NotImplementedError('bbox() is not implemented') def project_view(self, view):
+ raise NotImplementedError('project_view() is not implemented') def to_svg(self, round=0):
+ raise NotImplementedError('to_svg() is not implemented') def crop_to(self, view_bounds):
+ raise NotImplementedError('crop_to() is not implemented') def substract_geom(self, geom):
+ raise NotImplementedError('substract_geom() is not implemented yet') def is_emtpy(self):
+ return False def unify(self, point_store, precision=None):
+ raise NotImplementedError('unify() is not implemented yet')returns a list of point lists
+ def points(self): raise NotImplementedError('points() is not implemented yet') def join(self, geom):
+ raise NotImplementedError('join() is not implemented yet')base class for all solid geometry, e.g. polygons
+class SolidGeometry(Geometry):calculates area for this geometry
+ def area(): raise NotImplementedError('area() is not implemented')calculates centroid for this geometry
+ def centroid(): raise NotImplementedError('centroid() is not implemented') def invalidate(self):
+ self.__area = None
+
+from bbox import BBox
+from geometry import Geometrypolyline
+class PolyLine(Geometry): def __init__(self, contours):
+ self.contours = contoursreturns the bounding box of the line
+ def bbox(self): bbox = BBox()
+ for pts in self.contours:
+ for pt in pts:
+ bbox.update(pt)
+ return bboxprojects the line to a map projection
+ def project(self, proj): contours = []
+ for points in self.contours:
+ pts = []
+ for pt in points:
+ if proj._visible(pt[0], pt[1]):
+ p = proj.project(pt[0], pt[1])
+ if p is not None:
+ pts.append(p)
+ else:
+ if len(pts) > 0:
+ contours.append(pts)
+ pts = []
+ if len(pts) > 0:
+ contours.append(pts)
+ return PolyLine(contours)transforms the line to a new view
+ def project_view(self, view): contours = []
+ for points in self.contours:
+ pts = []
+ for pt in points:
+ p = view.project(pt)
+ pts.append(p)
+ contours.append(pts)
+ return PolyLine(contours)constructs a svg representation of this line
+ def to_svg(self, svg, round): path_str = ""
+ if round is False:
+ fmt = '%f,%f'
+ else:
+ fmt = '%.' + str(round) + 'f'
+ fmt = fmt + ',' + fmt
+
+ for points in self.contours:
+ path_str += "M "
+ pts = []
+ for pt in points:
+ pts.append(fmt % (pt[0], pt[1]))
+ path_str += 'L '.join(pts)
+
+ path = svg.node('path', d=path_str)
+ return path def crop_to(self, geom):skip
+ return self def is_empty(self):
+ return len(self.contours) == 0 def unify(self, point_store, precision):
+ from kartograph.simplify import unify_polygon
+ for i in range(len(self.contours)):
+ self.contours[i] = unify_polygon(self.contours[i], point_store, precision) def points(self):
+ return self.contoursis called after the points of this geometry have been +changed from outside this class
+ def update(self): passsimple line (= list of points)
+class Line(Geometry): def __init__(self, points):
+ self.pts = pointsreturns the bounding box of the line
+ def bbox(self): bbox = BBox()
+ for pt in self.pts:
+ bbox.update(pt)
+ return bboxprojects the line to a map projection
+ def project(self, proj): pts = []
+ for pt in self.pts:
+ p = proj.project(pt[0], pt[1])
+ if p is not None:
+ pts.append(p)
+ return Line(pts)transforms the line to a new view
+ def project_view(self, view): pts = []
+ for pt in self.pts:
+ p = view.project(pt)
+ pts.append(p)
+ return Line(pts)constructs a svg representation of this line
+ def to_svg(self, svg, round): path_str = ""
+ if round is False:
+ fmt = '%f,%f'
+ else:
+ fmt = '%.' + str(round) + 'f'
+ fmt = fmt + ',' + fmt
+
+ for pt in self.pts:
+ if path_str == "":
+ path_str = "M"
+ else:
+ path_str += "L"
+ path_str += fmt % (pt[0], pt[1])
+
+ path = svg.node('path', d=path_str)
+ return path def crop_to(self, geom):skip
+ return self def is_empty(self):
+ return len(self.pts) == 0 def unify(self, point_store, precision):
+ from kartograph.simplify import unify_polygon
+ self.pts = unify_polygon(self.pts, point_store, precision) def points(self):
+ return [self.pts]is called after the points of this geometry have been +changed from outside this class
+ def update(self): pass
+
+base class for points, used by line and bbox
+class Point(): def __init__(self, x, y):
+ self.x = x
+ self.y = yemulate python's container types
+ def project(self, proj):
+ (x, y) = proj.project(self.x, self.y)
+ self.x = x
+ self.y = y def __len__(self):
+ return 2 def __getitem__(self, k):
+ pt = (self.x, self.y)
+ return pt[k] def __setitem__(self, k, value):
+ if k == 0:
+ self.x = value
+ elif k == 1:
+ self.y = value
+ else:
+ raise IndexError def __delitem__(self, key):
+ raise TypeError('deletion not supported')
+
+from geometry import *
+from bbox import BBox
+import utilsDEPRECATED! we use shapely.geometry.MultiPolygon instead
+Complex polygons
+class MultiPolygon(SolidGeometry): def __init__(self, contours):
+ self.__area = None
+ self.__areas = None
+ self.__centroid = None
+ self.apply_contours(contours)constructs a Polygon from contours
+ @staticmethod
+ def fromPoly(poly):
+ contours = []
+ for i in range(len(poly)):
+ pts = poly.contour(i)
+ contours.append(pts)
+ return MultiPolygon(contours)
+
+ def apply_contours(self, contours): self.contours = contours
+ from Polygon import Polygon as GPCPoly
+ poly = GPCPoly()
+ skip = 0
+ for pts_ in contours:
+ pts = []
+ for pt in pts_:
+ if 'deleted' in pt and pt.deleted is True:
+ skip += 1
+ continue
+ pts.append((pt[0], pt[1]))
+ ishole = utils.is_clockwise(pts)
+
+ if len(pts) > 2:
+ poly.addContour(pts, ishole)
+ self.poly = poly def area(self):
+ if self.__area is not None:
+ return self.__area
+ self.__area = self.poly.area()
+ return self.__areareturns array of areas of all sub-polygons areas of holes are < 0
+ def areas(self): if self.__areas is not None:
+ return self.__areas
+ a = []
+ for i in range(len(self.poly)):
+ t = self.poly.area(i)
+ if self.poly.isHole(i):
+ t *= -1
+ a.append(t)
+ self.__areas = a
+ return areturns the center of gravity for this polygon
+ def centroid(self): if self.__centroid is not None:
+ return self.__centroid
+ self.__centroid = self.poly.center()
+ return self.__centroidsmart bounding box
+ def bbox(self, min_area=0): bb = []
+ bbox = BBox()
+ if min_area == 0:
+ bb.append(self.poly.boundingBox())
+ else:
+ areas = self.areas()
+ max_a = max(areas)
+ for i in range(len(self.poly)):
+ if self.poly.isHole(i):
+ continue
+ a = areas[i]
+ if a < max_a * min_area:
+ continue
+ bb.append(self.poly.boundingBox(i))
+ for b in bb:
+ bbox.update((b[0], b[2]))
+ bbox.update((b[1], b[2]))
+ bbox.update((b[0], b[3]))
+ bbox.update((b[1], b[3]))
+ return bboxreturns a new multi-polygon whose contours are +projected to a map projection
+ def project(self, proj): in_contours = self.contours
+ out_contours = []
+ for pts in in_contours:
+ pcont = proj.plot(pts)
+ if pcont != None:
+ out_contours += pcont
+ return MultiPolygon(out_contours)returns a new multi-polygon whose contours are +projected to a new view
+ def project_view(self, view): contours = self.contours
+ out = []
+ for contour in contours:
+ out_c = []
+ for pt in contour:
+ pt = view.project(pt)
+ out_c.append(pt)
+ out.append(out_c)
+ self.contours = out
+ out_poly = MultiPolygon(out)
+ return out_poly def crop_to(self, view_bounds):
+ poly = self.poly & view_bounds.poly
+ return MultiPolygon.fromPoly(poly) def substract_geom(self, geom):
+ if not isinstance(geom, MultiPolygon):
+ raise NotImplementedError('substraction is allowed for polygons only, yet')
+ poly = self.poly - geom.poly
+ return MultiPolygon.fromPoly(poly)constructs a svg representation of this polygon
+ def to_svg(self, svg, round): path_str = ""
+ if round is False:
+ fmt = '%f,%f'
+ else:
+ fmt = '%.' + str(round) + 'f'
+ fmt = fmt + ',' + fmt
+
+ for pts in self.contours:
+ cont_str = ""
+ kept = []
+ for pt in pts:
+ if 'deleted' in pt and pt.deleted is True:
+ continue
+ kept.append(pt)
+
+ if len(kept) <= 3:
+ continue
+ for pt in kept:
+ if cont_str == "":
+ cont_str = "M"
+ else:
+ cont_str += "L"
+ cont_str += fmt % pt
+ cont_str += "Z "
+ path_str += cont_str
+
+ if path_str == "":
+ return None
+
+ path = svg.node('path', d=path_str)
+ return path def is_empty(self):
+ return len(self.contours) == 0 def unify(self, point_store, precision=None):
+ from kartograph.simplify import unify_polygons
+ contours = self.contours
+ contours = unify_polygons(contours, point_store, precision)
+ self.apply_contours(contours) def points(self):
+ return self.contoursis called after the points of this geometry are +changed from outside this class
+ def update(self): self.apply_contours(self.contours) def to_kml(self, round=None):
+ from pykml.factory import KML_ElementMaker as KML
+ poly = KML.Polygon(
+ KML.tesselate("1")
+ )
+ outer = KML.outerBoundaryIs()
+ inner = KML.innerBoundaryIs()
+ has_inner = False
+ for i in range(len(self.poly)):
+ cnt = self.poly[i]
+ coords = ''
+ for p in cnt:
+ coords += ','.join(map(str, p)) + ' '
+ ring = KML.LinearRing(
+ KML.coordinates(coords)
+ )hole = self.poly.isHole(i) +if hole == False:
+ outer.append(ring)else: + inner.append(ring) + has_inner = True
+ poly.append(outer)
+ if has_inner:
+ poly.append(inner)
+
+ return polysplits the geometry into several line segments
+ def to_line_segments(self): self.lines = [] def from_line_segments(self):
+ self.lines = []class Polygon(SolidGeometry): def __init__(self, points):
+ self.__area = None
+ self.__centroid = None
+ self.points = points def area(self):
+ if self.__area is not None:
+ return self.__area
+ a = 0
+ pts = self.points
+ for i in range(len(pts) - 1):
+ p0 = pts[i]
+ p1 = pts[i + 1]
+ a += p0.x * p1.y - p1.x * p0.y
+ self.__area = abs(a) * .5
+ return self.__area def project(self, proj):
+ self.invalidate()
+ self.points = proj.plot(self.points)
+
+geometry utils
+returns true if a given linear ring is in clockwise order
+def is_clockwise(pts): s = 0
+ for i in range(len(pts) - 1):
+ if 'x' in pts[i]:
+ x1 = pts[i].x
+ y1 = pts[i].y
+ x2 = pts[i + 1].x
+ y2 = pts[i + 1].y
+ else:
+ x1, y1 = pts[i]
+ x2, y2 = pts[i + 1]
+ s += (x2 - x1) * (y2 + y1)
+ return s >= 0def bbox_to_polygon(bbox):
+ from shapely.geometry import Polygon
+ s = bbox
+ poly = Polygon([(s.left, s.bottom), (s.left, s.top), (s.right, s.top), (s.right, s.bottom)])
+ return polydef geom_to_bbox(geom, min_area=0):
+ from kartograph.geometry import BBox
+ from shapely.geometry import MultiPolygon
+ if True or min_area == 0 or not isinstance(geom, MultiPolygon):if no minimum area ratio is set or the geometry +is not a multipart geometry, we simply use the +full bbox
+ minx, miny, maxx, maxy = geom.bounds
+ return BBox(width=maxx - minx, height=maxy - miny, left=minx, top=miny)
+ else:for multipart geometry we use only the bbox of +the 'biggest' sub-geometries, depending on min_area
+ bbox = BBox()
+ areas = []
+ bb = []
+ for polygon in geom.geoms:
+ areas.append(polygon.area)
+ max_a = max(areas)
+ for i in range(len(geom.geoms)):
+ a = areas[i]
+ if a < max_a * min_area:ignore this sub polygon since it is too small
+ continue
+ bb.append(geom.geoms[i].bounds)
+ for b in bb:
+ bbox.update((b[0], b[2]))
+ bbox.update((b[1], b[2]))
+ bbox.update((b[0], b[3]))
+ bbox.update((b[1], b[3]))
+ return bboxjoins polygonal features
+def join_features(features, props): from feature import MultiPolygonFeature
+
+ if len(features) == 0:
+ return features
+
+ joined = []
+ polygons = []
+
+ for feat in features:
+ if isinstance(feat, MultiPolygonFeature):
+ polygons.append(feat.geom)
+ else:
+ joined.append(feat) # cannot join this
+
+ polygons = filter(lambda x: x is not None, polygons)
+ if len(polygons) > 0:
+ poly = polygons[0]
+ for poly2 in polygons[1:]:
+ poly = poly.union(poly2)
+ joined.append(MultiPolygonFeature(poly, props))
+ return joined
+
+from shapely.geometry import Polygon, MultiPolygon, LineString, MultiLineString, MultiPoint, Point
+from kartograph.errors import KartographErrortranslates a point to a view
+class View(object): def __init__(self, bbox=None, width=None, height=None, padding=0):
+ self.bbox = bbox
+ self.width = width
+ self.padding = padding
+ self.height = height
+ if bbox:
+ self.scale = min((width - padding * 2) / bbox.width, (height - padding * 2) / bbox.height) def project(self, pt):
+ bbox = self.bbox
+ if not bbox:
+ return pt
+ s = self.scale
+ h = self.height
+ w = self.width
+ px = pt[0]
+ py = pt[1]
+ x = (px - bbox.left) * s + (w - bbox.width * s) * .5
+ y = (py - bbox.top) * s + (h - bbox.height * s) * .5
+ return ((x, y), Point(x, y))[isinstance(pt, Point)] def project_inverse(self, pt):
+ bbox = self.bbox
+ if not bbox:
+ return pt
+ s = self.scale
+ h = self.height
+ w = self.width
+ x = pt[0]
+ y = pt[1]
+ px = (x - (w - bbox.width * s) * .5) / s + bbox.left
+ py = (y - (h - bbox.height * s) * .5) / s + bbox.top
+ return ((px, py), Point(px, py))[isinstance(pt, Point)]converts the given geometry to the view coordinates
+ def project_geometry(self, geometry): geometries = hasattr(geometry, 'geoms') and geometry.geoms or [geometry]
+ res = []at first shift polygons +geometries = [] +for geom in unshifted_geometries: + geometries += self._shift_polygon(geom)
+ for geom in geometries:
+ if isinstance(geom, Polygon):
+ res += self.project_polygon(geom)
+ elif isinstance(geom, LineString):
+ rings = self.project_linear_ring(geom.coords)
+ res += map(LineString, rings)
+ elif isinstance(geom, Point):
+ if self._visible(geom.x, geom.y):
+ res.append(self.project(geom))
+ else:
+ raise KartographError('unknown geometry type %s' % geometry)
+
+ if len(res) > 0:
+ if isinstance(res[0], Polygon):
+ if len(res) > 1:
+ return MultiPolygon(res)
+ else:
+ return res[0]
+ elif isinstance(res[0], LineString):
+ if len(res) > 1:
+ return MultiLineString(res)
+ else:
+ return LineString(res[0])
+ else:
+ if len(res) > 1:
+ return MultiPoint(res)
+ else:
+ return res[0] def project_polygon(self, polygon):
+ ext = self.project_linear_ring(polygon.exterior)
+ if len(ext) == 1:
+ pts_int = []
+ for interior in polygon.interiors:
+ pts_int += self.project_linear_ring(interior)
+ return [Polygon(ext[0], pts_int)]
+ elif len(ext) == 0:
+ return []
+ else:
+ raise KartographError('unhandled case: exterior is split into multiple rings') def project_linear_ring(self, ring):
+ points = []
+ for pt in ring.coords:
+ x, y = self.project(pt)
+ points.append((x, y))
+ return [points] def __str__(self):
+ return 'View(w=%f, h=%f, pad=%f, scale=%f, bbox=%s)' % (self.width, self.height, self.padding, self.scale, self.bbox)
+
+from options import parse_options
+from shapely.geometry import Polygon, LineString, MultiPolygon
+from errors import *
+from copy import deepcopy
+from renderer import SvgRenderer, KmlRenderer
+from map import Mapverbose = FalseThese renderers are currently available. See renderer/svg.py and renderer/kml.py +for more details on those.
+_known_renderer = {
+ 'svg': SvgRenderer,
+ 'kml': KmlRenderer
+}class Kartograph(object): def __init__(self):
+ self.layerCache = {}
+ passGenerates a the map and renders it using the specified output format.
+ def generate(self, opts, outfile=None, format='svg', preview=None): if preview is None:
+ preview = outfile is not NoneCreate a deep copy of the options dictionary so our changes will not be +visible to the calling application.
+ opts = deepcopy(opts)Parse the options dictionary. See options.py for more details.
+ parse_options(opts)Create the map instance. It will do all the hard work for us, so you +definitely should check out map.py for all the fun stuff happending +there..
+ _map = Map(opts, self.layerCache, format=format)Check if the format is handled by a renderer.
+ format = format.lower()
+ if format in _known_renderer:Create a renderer instance and render the map.
+ renderer = _known_renderer[format](_map)
+ renderer.render()If requested, we try to preview the created map now, which means that +we open the created SVG file in Firefox or open the KML in Google Earth. +Of course, preview modes are highly dependent on the operating system, but +we don't care about that now.
+ if preview:
+ renderer.preview()Write the map to a file or return the renderer instance.
+ if outfile is None:
+ return renderer
+ else:
+ renderer.write(outfile)
+ else:
+ raise KartographError('unknown format: %s' % format)Here are some handy methods for debugging Kartograph. It will plot a given shapely +geometry using matplotlib and descartes.
+def _plot_geometry(geom, fill='#ffcccc', stroke='#333333', alpha=1, msg=None):
+ from matplotlib import pyplot
+ from matplotlib.figure import SubplotParams
+ from descartes import PolygonPatch
+
+ if isinstance(geom, (Polygon, MultiPolygon)):
+ b = geom.bounds
+ geoms = hasattr(geom, 'geoms') and geom.geoms or [geom]
+ w, h = (b[2] - b[0], b[3] - b[1])
+ ratio = w / h
+ pad = 0.15
+ fig = pyplot.figure(1, figsize=(5, 5 / ratio), dpi=110, subplotpars=SubplotParams(left=pad, bottom=pad, top=1 - pad, right=1 - pad))
+ ax = fig.add_subplot(111, aspect='equal')
+ for geom in geoms:
+ patch1 = PolygonPatch(geom, linewidth=0.5, fc=fill, ec=stroke, alpha=alpha, zorder=0)
+ ax.add_patch(patch1)
+ p = (b[2] - b[0]) * 0.03 # some padding
+ pyplot.axis([b[0] - p, b[2] + p, b[3] + p, b[1] - p])
+ pyplot.grid(True)
+ if msg:
+ fig.suptitle(msg, y=0.04, fontsize=9)
+ pyplot.show()def _plot_lines(lines):
+ from matplotlib import pyplot def plot_line(ax, line):
+ filtered = []
+ for pt in line:
+ if not pt.deleted:
+ filtered.append(pt)
+ if len(filtered) < 2:
+ return
+ ob = LineString(line)
+ x, y = ob.xy
+ ax.plot(x, y, '-', color='#333333', linewidth=0.5, solid_capstyle='round', zorder=1)
+
+ fig = pyplot.figure(1, figsize=(4, 5.5), dpi=90, subplotpars=SubplotParams(left=0, bottom=0.065, top=1, right=1))
+ ax = fig.add_subplot(111, aspect='equal')
+ for line in lines:
+ plot_line(ax, line)
+ pyplot.grid(False)
+ ax.xaxis.set_visible(False)
+ ax.yaxis.set_visible(False)
+ ax.set_frame_on(False)
+ return (ax, fig)def _debug_show_features(features, message=None):
+ from descartes import PolygonPatch
+ from matplotlib import pyplot
+ from matplotlib.figure import SubplotParams
+
+ fig = pyplot.figure(1, figsize=(9, 5.5), dpi=110, subplotpars=SubplotParams(left=0, bottom=0.065, top=1, right=1))
+ ax = fig.add_subplot(111, aspect='equal')
+ b = (100000, 100000, -100000, -100000)
+ for feat in features:
+ if feat.geom is None:
+ continue
+ c = feat.geom.bounds
+ b = (min(c[0], b[0]), min(c[1], b[1]), max(c[2], b[2]), max(c[3], b[3]))
+ geoms = hasattr(feat.geom, 'geoms') and feat.geom.geoms or [feat.geom]
+ for geom in geoms:
+ patch1 = PolygonPatch(geom, linewidth=0.25, fc='#ddcccc', ec='#000000', alpha=0.75, zorder=0)
+ ax.add_patch(patch1)
+ p = (b[2] - b[0]) * 0.05 # some padding
+ pyplot.axis([b[0] - p, b[2] + p, b[3], b[1] - p])
+ ax.xaxis.set_visible(False)
+ ax.yaxis.set_visible(False)
+ ax.set_frame_on(True)
+ if message:
+ fig.suptitle(message, y=0.04, fontsize=9)
+ pyplot.show()
+
+as of version 2.0 kartograph supports multiple import formats
+__all__ = ['LayerSource', 'ShapefileLayer', 'GraticuleLayer']
+
+from shplayer import ShapefileLayer
+from layersource import LayerSource
+from special import GraticuleLayer, SeaLayer
+from kartograph.errors import *def handle_layer_source(layer, cache={}):
+ if 'src' in layer:
+ src = layer['src']
+ if src in cache:
+ return cache[src]
+ if src[-4:].lower() == ".shp": # shapefile layer
+ src = ShapefileLayer(src)
+ if isinstance(src, LayerSource):
+ cache[layer['src']] = src
+ return src
+ else:
+ raise KartographLayerSourceError('don\'t know how to handle "' + src + '"')
+ elif 'special' in layer:
+ if layer['special'] == 'graticule':
+ return GraticuleLayer()
+ elif layer['special'] == 'sea':
+ return SeaLayer()
+
+base class for layer source data providers (e.g. shapefiles)
+class LayerSource: def get_features(self, attr=None, filter=None, bbox=None):
+ raise NotImplementedError()
+
+shapefile.py
+Provides read and write support for ESRI Shapefiles.
+author: jlawhead
from struct import pack, unpack, calcsize, error
+import os
+import sys
+import time
+import arrayConstants for shape types
+NULL = 0
+POINT = 1
+POLYLINE = 3
+POLYGON = 5
+MULTIPOINT = 8
+POINTZ = 11
+POLYLINEZ = 13
+POLYGONZ = 15
+MULTIPOINTZ = 18
+POINTM = 21
+POLYLINEM = 23
+POLYGONM = 25
+MULTIPOINTM = 28
+MULTIPATCH = 31
+
+PYTHON3 = sys.version_info[0] == 3def b(v):
+ if PYTHON3:
+ if isinstance(v, str):For python 3 encode str to bytes.
+ return v.encode('utf-8')
+ elif isinstance(v, bytes):Already bytes.
+ return v
+ else:Error.
+ raise Exception('Unknown input type')
+ else:For python 2 assume str passed in and return str.
+ return vdef u(v):
+ if PYTHON3:
+ if isinstance(v, bytes):For python 3 decode bytes to str.
+ return v.decode('utf-8')
+ elif isinstance(v, str):Already str.
+ return v
+ else:Error.
+ raise Exception('Unknown input type')
+ else:For python 2 assume str passed in and return str.
+ return vdef is_string(v):
+ if PYTHON3:
+ return isinstance(v, str)
+ else:
+ return isinstance(v, basestring)Converts python tuples to lits of the appropritate type.
+class _Array(array.array):Used to unpack different shapefile header parts.""" +def repr(self): + return str(self.tolist())
+s _Shape: +def init(self, shapeType=None): +Stores the geometry of the different shape types
+ specified in the Shapefile spec. Shape types are
+ usually point, polyline, or polygons. Every shape type
+ except the "Null" type contains points at some level for
+ example verticies in a polygon. If a shape type has
+ multiple shapes containing points within a single
+ geometry record then those shapes are called parts. Parts
+ are designated by their starting index in geometry record's
+ list of shapes."""
+ self.shapeType = shapeType
+ self.points = []
+#DIVIDER
+class _ShapeRecord:
+#DIVIDER
+ def __init__(self, shape=None, record=None):
+ self.shape = shape
+ self.record = record
+#DIVIDER
+class ShapefileException(Exception):
+#DIVIDER
+ pass
+#DIVIDER
+class Reader:
+#DIVIDER
+ def __init__(self, *args, **kwargs):
+ self.shp = None
+ self.shx = None
+ self.dbf = None
+ self.shapeName = "Not specified"
+ self._offsets = []
+ self.shpLength = None
+ self.numRecords = None
+ self.fields = []
+ self.__dbfHdrLength = 0
+#DIVIDER
+ if len(args) > 0:
+ if type(args[0]) is type("stringTest"):
+ self.load(args[0])
+ return
+ if "shp" in kwargs.keys():
+ if hasattr(kwargs["shp"], "read"):
+ self.shp = kwargs["shp"]
+ if hasattr(self.shp, "seek"):
+ self.shp.seek(0)
+ if "shx" in kwargs.keys():
+ if hasattr(kwargs["shx"], "read"):
+ self.shx = kwargs["shx"]
+ if hasattr(self.shx, "seek"):
+ self.shx.seek(0)
+ if "dbf" in kwargs.keys():
+ if hasattr(kwargs["dbf"], "read"):
+ self.dbf = kwargs["dbf"]
+ if hasattr(self.dbf, "seek"):
+ self.dbf.seek(0)
+ if self.shp or self.dbf:
+ self.load()
+ else:
+ raise ShapefileException("Shapefile Reader requires a shapefile or file-like object.")
+#DIVIDER
+ def load(self, shapefile=None):
+#DIVIDER
+ available. If not a ShapefileException is raised."""
+ if not f:
+ raise ShapefileException("Shapefile Reader requires a shapefile or file-like object.")
+ if self.shp and self.shpLength is None:
+ self.load()
+ if self.dbf and len(self.fields) == 0:
+ self.load()
+ return fA shape object of any type.
+ def __restrictIndex(self, i): if not self.shp:
+ raise ShapefileException("Shapefile Reader requires a shapefile or file-like object. (no shp file found")
+ shp = self.shpAn exception to handle shapefile specific problems.
+ shp.seek(24)
+ self.shpLength = unpack(">i", shp.read(4))[0] * 2 shp.seek(32)
+ self.shapeType= unpack("<i", shp.read(4))[0]Reads the three files of a shapefile as a unit or
+ self.bbox = _Array('d', unpack("<4d", shp.read(32)))separately. If one of the three files (.shp, .shx, +.dbf) is missing no exception is thrown until you try +to call a method that depends on that particular file. +The .shx index file is used if available for efficiency +but is not required to read the geometry from the .shp +file. The "shapefile" argument in the constructor is the +name of the file you want to open.
+You can instantiate a Reader without specifying a shapefile +and then specify one later with the load() method.
+Only the shapefile headers are read upon loading. Content +within each file is only accessed when required and as +efficiently as possible. Shapefiles are usually not large +but they can be.
+ self.elevation = _Array('d', unpack("<2d", shp.read(16)))See if a shapefile name was passed as an argument
+ self.measure = _Array('d', unpack("<2d", shp.read(16)))Opens a shapefile from a filename or file-like
+ def __shape(self):object. Normally this method would be called by the +constructor with the file object or file name as an +argument.""" +if shapefile: + (shapeName, ext) = os.path.splitext(shapefile) + self.shapeName = shapeName + try: + self.shp = open("%s.shp" % shapeName, "rb") + except IOError: + raise ShapefileException("Unable to open %s.shp" % shapeName) + try: + self.shx = open("%s.shx" % shapeName, "rb") + except IOError: + raise ShapefileException("Unable to open %s.shx" % shapeName) + try: + self.dbf = open("%s.dbf" % shapeName, "rb") + except IOError: + raise ShapefileException("Unable to open %s.dbf" % shapeName) +if self.shp: + self.shpHeader() +if self.dbf: + self.dbfHeader()
+__getFileObj(self, f): +Checks to see if the requested shapefile file object is
+ f = self.__getFileObj(self.shp)
+ record = _Shape()
+ nParts = nPoints = zmin = zmax = mmin = mmax = None
+ (recNum, recLength) = unpack(">2i", f.read(8))
+ shapeType = unpack("<i", f.read(4))[0]
+ record.shapeType = shapeTypeProvides list-like handling of a record index with a clearer
+ if shapeType == 0:
+ record.points = []error message if the index is out of bounds.""" +if self.numRecords: + rmax = self.numRecords - 1 + if abs(i) > rmax: + raise IndexError("Shape or Record index out of range.") + if i < 0: i = range(self.numRecords)[i] +return i
+__shpHeader(self): +Reads the header information from a .shp or .shx file.
+ elif shapeType in (3,5,8,13,15,18,23,25,28,31):
+ record.bbox = _Array('d', unpack("<4d", f.read(32)))File length (16-bit word * 2 = bytes)
+ if shapeType in (3,5,13,15,23,25,31):
+ nParts = unpack("<i", f.read(4))[0]Shape type
+ if shapeType in (3,5,8,13,15,23,25,31):
+ nPoints = unpack("<i", f.read(4))[0]The shapefile's bounding box (lower left, upper right)
+ if nParts:
+ record.parts = _Array('i', unpack("<%si" % nParts, f.read(nParts * 4)))Elevation
+ if shapeType == 31:
+ record.partTypes = _Array('i', unpack("<%si" % nParts, f.read(nParts * 4)))Measure
+ if nPoints:
+ record.points = [_Array('d', unpack("<2d", f.read(16))) for p in range(nPoints)]Returns the header info and geometry for a single shape.
+ if shapeType in (13,15,18,31):
+ (zmin, zmax) = unpack("<2d", f.read(16))
+ record.z = _Array('d', unpack("<%sd" % nPoints, f.read(nPoints * 8))) if shapeType in (13,18,23,25,28,31):
+ (mmin, mmax) = unpack("<2d", f.read(16))For Null shapes create an empty points list for consistency
+ record.m = []
+ for m in _Array('d', unpack("%sd" % nPoints, f.read(nPoints * 8))):
+ if m > -10e38:
+ record.m.append(m)
+ else:
+ record.m.append(None)All shape types capable of having a bounding box
+ if shapeType in (1,11,21):
+ record.points = [_Array('d', unpack("<2d", f.read(16)))]Shape types with parts
+ if shapeType == 11:
+ record.z = unpack("<d", f.read(8))Shape types with points
+ if shapeType in (11,21):
+ record.m = unpack("<d", f.read(8))
+ return recordRead parts
+ def __shapeIndex(self, i=None):Read part types for Multipatch - 31
+ record file."""
+ shp = self.__getFileObj(self.shp)
+ i = self.__restrictIndex(i)
+ offset = self.__shapeIndex(i)
+ if not offset:
+#DIVIDER
+ shapes = self.shapes()
+ return shapes[i]
+ shp.seek(offset)
+ return self.__shape()
+#DIVIDER
+ def shapes(self):
+#DIVIDER
+ shp = self.__getFileObj(self.shp)
+ shp.seek(100)
+ shapes = []
+ while shp.tell() < self.shpLength:
+ shapes.append(self.__shape())
+ return shapes
+#DIVIDER
+ def __dbfHeaderLength(self):
+#DIVIDER
+ if not self.__dbfHdrLength:
+ if not self.dbf:
+ raise ShapefileException("Shapefile Reader requires a shapefile or file-like object. (no dbf file found)")
+ dbf = self.dbf
+ (self.numRecords, self.__dbfHdrLength) = \
+ unpack("<xxxxLH22x", dbf.read(32))
+ return self.__dbfHdrLength
+#DIVIDER
+ def __dbfHeader(self):
+#DIVIDER
+ if not self.dbf:
+ raise ShapefileException("Shapefile Reader requires a shapefile or file-like object. (no dbf file found)")
+ dbf = self.dbf
+ headerLength = self.__dbfHeaderLength()
+ numFields = (headerLength - 33) // 32
+ for field in range(numFields):
+ fieldDesc = list(unpack("<11sc4xBB14x", dbf.read(32)))
+ name = 0
+ idx = 0
+ if b("\x00") in fieldDesc[name]:
+ idx = fieldDesc[name].index(b("\x00"))
+ else:
+ idx = len(fieldDesc[name]) - 1
+ fieldDesc[name] = fieldDesc[name][:idx]
+ fieldDesc[name] = u(fieldDesc[name])
+ fieldDesc[name] = fieldDesc[name].lstrip()
+ fieldDesc[1] = u(fieldDesc[1])
+ self.fields.append(fieldDesc)
+ terminator = dbf.read(1)
+ assert terminator == b("\r")
+ self.fields.insert(0, ('DeletionFlag', 'C', 1, 0))
+#DIVIDER
+ def __recordFmt(self):
+#DIVIDER
+ if not self.numRecords:
+ self.__dbfHeader()
+ fmt = ''.join(['%ds' % fieldinfo[2] for fieldinfo in self.fields])
+ fmtSize = calcsize(fmt)
+ return (fmt, fmtSize)
+#DIVIDER
+ def __record(self):
+#DIVIDER
+ f = self.__getFileObj(self.dbf)
+ recFmt = self.__recordFmt()
+ recordContents = unpack(recFmt[0], f.read(recFmt[1]))
+ if recordContents[0] != b(' '):
+#DIVIDER
+ return None
+ record = []
+ for (name, typ, size, deci), value in zip(self.fields,
+ recordContents):
+ if name == 'DeletionFlag':
+ continue
+ elif not value.strip():
+ record.append(value)
+ continue
+ elif typ == "N":
+ value = value.replace(b('\0'), b('')).strip()
+ if value == b(''):
+ value = 0
+ elif deci:
+ value = float(value)
+ else:
+ value = int(value)
+ elif typ == b('D'):
+ try:
+ y, m, d = int(value[:4]), int(value[4:6]), int(value[6:8])
+ value = [y, m, d]
+ except:
+ value = value.strip()
+ elif typ == b('L'):
+ value = (value in b('YyTt') and b('T')) or \
+ (value in b('NnFf') and b('F')) or b('?')
+ else:
+ value = u(value)
+ value = value.strip()
+ record.append(value)
+ return record
+#DIVIDER
+ def record(self, i=0):
+#DIVIDER
+ f = self.__getFileObj(self.dbf)
+ if not self.numRecords:
+ self.__dbfHeader()
+ i = self.__restrictIndex(i)
+ recSize = self.__recordFmt()[1]
+ f.seek(0)
+ f.seek(self.__dbfHeaderLength() + (i * recSize))
+ return self.__record()
+#DIVIDER
+ def records(self):
+#DIVIDER
+ if not self.numRecords:
+ self.__dbfHeader()
+ records = []
+ f = self.__getFileObj(self.dbf)
+ f.seek(self.__dbfHeaderLength())
+ for i in range(self.numRecords):
+ r = self.__record()
+ if r:
+ records.append(r)
+ return records
+#DIVIDER
+ def shapeRecord(self, i=0):
+#DIVIDER
+ all records in a shapefile."""
+ shapeRecords = []
+ return [_ShapeRecord(shape=rec[0], record=rec[1]) \
+ for rec in zip(self.shapes(), self.records())]Read points - produces a list of [x,y] values
+class Writer:Read z extremes and values
+ def __init__(self, shapeType=None):
+ self._shapes = []
+ self.fields = []
+ self.records = []
+ self.shapeType = shapeType
+ self.shp = None
+ self.shx = None
+ self.dbf = NoneRead m extremes and values
+ self._offsets = []
+ self._lengths = []Measure values less than -10e38 are nodata values according to the spec
+ self.deletionFlag = 0Read a single point
+ def __getFileObj(self, f):Read a single Z value
+ if not f:
+ raise ShapefileException("No file-like object available.")
+ elif hasattr(f, "write"):
+ return f
+ else:
+ pth = os.path.split(f)[0]
+ if pth and not os.path.exists(pth):
+ os.makedirs(pth)
+ return open(f, "wb")Read a single M value
+ def __shpFileLength(self):Returns the offset in a .shp file for a shape based on information
+ size = 100in the .shx index file.""" +shx = self.shx +if not shx: + return None +if not self._offsets: + # File length (16-bit word * 2 = bytes) - header length + shx.seek(24) + shxRecordLength = (unpack(">i", shx.read(4))[0] * 2) - 100 + numRecords = shxRecordLength // 8 + # Jump to the first record. + shx.seek(100) + for r in range(numRecords): + # Offsets are 16-bit words just like the file length + self._offsets.append(unpack(">i", shx.read(4))[0] * 2) + shx.seek(shx.tell() + 4) +if not i == None: + return self._offsets[i]
+shape(self, i=0): +Returns a shape object for a shape in the the geometry
+ for s in self._shapes:Shx index not available so use the full list.
+ size += 12Returns all shapes in a shapefile.
+ if hasattr(s,'parts'):
+ nParts = len(s.parts)
+ if hasattr(s,'points'):
+ nPoints = len(s.points) if self.shapeType in (3,5,8,13,15,18,23,25,28,31):
+ size += 32Retrieves the header length of a dbf file header.
+ if self.shapeType in (3,5,13,15,23,25,31): size += 4Reads a dbf header. Xbase-related code borrows heavily from ActiveState Python Cookbook Recipe 362715 by Raymond Hettinger
+ size += nParts * 4 if self.shapeType in (3,5,8,13,15,23,25,31):Calculates the size of a .shp geometry record.
+ size += 4 size += 16 * nPointsReads and returns a dbf record row as a list of values.
+ if self.shapeType == 31:
+ size += nParts * 4 if self.shapeType in (13,15,18,31):deleted record
+ size += 16Returns a specific dbf record based on the supplied index.
+ size += 8 * nPoints if self.shapeType in (23,25,31):Returns all records in a dbf file.
+ size += 16 size += 8 * nPointsReturns a combination geometry and attribute record for the
+ if self.shapeType in (1,11,21):
+ size += 16supplied record index.""" +i = self.__restrictIndex(i) +return _ShapeRecord(shape=self.shape(i), + record=self.record(i))
+shapeRecords(self): +Returns a list of combination geometry/attribute records for
+ if self.shapeType == 11:
+ size += 8Provides write support for ESRI Shapefiles.
+ if self.shapeType in (11,21):
+ size += 8 size //= 2
+ return sizeGeometry record offsets and lengths for writing shx file.
+ def __bbox(self, shapes, shapeTypes=[]):
+ x = []
+ y = []
+ for s in shapes:
+ shapeType = self.shapeType
+ if shapeTypes:
+ shapeType = shapeTypes[shapes.index(s)]
+ px, py = list(zip(*s.points))[:2]
+ x.extend(px)
+ y.extend(py)
+ return [min(x), min(y), max(x), max(y)]Use deletion flags in dbf? Default is false (0).
+ def __zbox(self, shapes, shapeTypes=[]):
+ z = []
+ for s in shapes:
+ try:
+ for p in s.points:
+ z.append(p[2])
+ except IndexError:
+ pass
+ if not z: z.append(0)
+ return [min(z), max(z)]Safety handler to verify file-like objects
+ def __mbox(self, shapes, shapeTypes=[]):
+ m = [0]
+ for s in shapes:
+ try:
+ for p in s.points:
+ m.append(p[3])
+ except IndexError:
+ pass
+ return [min(m), max(m)] def bbox(self):Calculates the file length of the shp file.
+ return self.__zbox(self._shapes)Start with header length
+ def mbox(self):Calculate size of all shapes
+ return self.__mbox(self._shapes)Add in record header and shape type fields
+ def __shapefileHeader(self, fileObj, headerType='shp'):nParts and nPoints do not apply to all shapes +if self.shapeType not in (0,1): + nParts = len(s.parts) + nPoints = len(s.points)
+ f = self.__getFileObj(self.dbf)
+ f.seek(0)
+ version = 3
+ year, month, day = time.localtime()[:3]
+ year -= 1900All shape types capable of having a bounding box
+ for field in self.fields:
+ if field[0].startswith("Deletion"):
+ self.fields.remove(field)
+ numRecs = len(self.records)
+ numFields = len(self.fields)
+ headerLength = numFields * 32 + 33
+ recordLength = sum([int(field[2]) for field in self.fields]) + 1
+ header = pack('<BBBBLHH20x', version, year, month, day, numRecs,
+ headerLength, recordLength)
+ f.write(header)Shape types with parts
+ for field in self.fields:
+ name, fieldType, size, decimal = field
+ name = b(name)
+ name = name.replace(b(' '), b('_'))
+ name = name.ljust(11).replace(b(' '), b('\x00'))
+ fieldType = b(fieldType)
+ size = int(size)
+ fld = pack('<11sc4xBB14x', name, fieldType, size, decimal)
+ f.write(fld)Parts count
+ f.write(b('\r'))Parts index array
+ def __shpRecords(self):Shape types with points
+ f = self.__getFileObj(self.shp)
+ f.seek(100)
+ recNum = 1
+ for s in self._shapes:
+ self._offsets.append(f.tell())Points count
+ f.write(pack(">2i", recNum, 0))
+ recNum += 1
+ start = f.tell()Points array
+ f.write(pack("<i", s.shapeType))Calc size of part types for Multipatch (31)
+ if s.shapeType in (3,5,8,13,15,18,23,25,28,31):
+ try:
+ f.write(pack("<4d", *self.__bbox([s])))
+ except error:
+ raise ShapefileException("Falied to write bounding box for record %s. Expected floats." % recNum)Calc z extremes and values
+ if s.shapeType in (3,5,13,15,23,25,31):z extremes
+ f.write(pack("<i", len(s.parts)))z array
+ if s.shapeType in (3,5,8,13,15,23,25,31):Calc m extremes and values
+ f.write(pack("<i", len(s.points)))m extremes
+ if s.shapeType in (3,5,13,15,23,25,31):
+ for p in s.parts:
+ f.write(pack("<i", p))m array
+ if s.shapeType == 31:
+ for pt in s.partTypes:
+ f.write(pack("<i", pt))Calc a single point
+ if s.shapeType in (3,5,8,13,15,23,25,31):
+ try:
+ [f.write(pack("<2d", *p[:2])) for p in s.points]
+ except error:
+ raise ShapefileException("Failed to write points for record %s. Expected floats." % recNum)Calc a single Z value
+ if s.shapeType in (13,15,18,31):
+ try:
+ f.write(pack("<2d", *self.__zbox([s])))
+ except error:
+ raise ShapefileException("Failed to write elevation extremes for record %s. Expected floats." % recNum)
+ try:
+ [f.write(pack("<d", p[2])) for p in s.points]
+ except error:
+ raise ShapefileException("Failed to write elevation values for record %s. Expected floats." % recNum)Calc a single M value
+ if s.shapeType in (23,25,31):
+ try:
+ f.write(pack("<2d", *self.__mbox([s])))
+ except error:
+ raise ShapefileException("Failed to write measure extremes for record %s. Expected floats" % recNum)
+ try:
+ [f.write(pack("<d", p[3])) for p in s.points]
+ except error:
+ raise ShapefileException("Failed to write measure values for record %s. Expected floats" % recNum)Calculate size as 16-bit words
+ if s.shapeType in (1,11,21):
+ try:
+ f.write(pack("<2d", s.points[0][0], s.points[0][1]))
+ except error:
+ raise ShapefileException("Failed to write point for record %s. Expected floats." % recNum) if s.shapeType == 11:
+ try:
+ f.write(pack("<1d", s.points[0][2]))
+ except error:
+ raise ShapefileException("Failed to write elevation value for record %s. Expected floats." % recNum) if s.shapeType in (11,21):
+ try:
+ f.write(pack("<1d", s.points[0][3]))
+ except error:
+ raise ShapefileException("Failed to write measure value for record %s. Expected floats." % recNum) finish = f.tell()
+ length = (finish - start) // 2
+ self._lengths.append(length)Returns the current bounding box for the shapefile which is
+ f.seek(start-4)
+ f.write(pack(">i", length))
+ f.seek(finish)the lower-left and upper-right corners. It does not contain the +elevation or measure extremes.""" +return self.__bbox(self._shapes)
+zbox(self): +Returns the current z extremes for the shapefile.
+ def __shxRecords(self):Returns the current m extremes for the shapefile.
+ f = self.__getFileObj(self.shx)
+ f.seek(100)
+ for i in range(len(self._shapes)):
+ f.write(pack(">i", self._offsets[i] // 2))
+ f.write(pack(">i", self._lengths[i])) def __dbfRecords(self):Writes the specified header type to the specified file-like object.
+ f = self.__getFileObj(self.dbf)
+ for record in self.records:
+ if not self.fields[0][0].startswith("Deletion"):
+ f.write(b(' ')) # deletion flag
+ for (fieldName, fieldType, size, dec), value in zip(self.fields, record):
+ fieldType = fieldType.upper()
+ size = int(size)
+ if fieldType.upper() == "N":
+ value = str(value).rjust(size)
+ elif fieldType == 'L':
+ value = str(value)[0].upper()
+ else:
+ value = str(value)[:size].ljust(size)
+ assert len(value) == size
+ value = b(value)
+ f.write(value)Several of the shapefile formats are so similar that a single generic +method to read or write them is warranted.""" +f = self.__getFileObj(fileObj) +f.seek(0)
+f.write(pack(">6i", 9994,0,0,0,0,0))
+if headerType == 'shp': + f.write(pack(">i", self.__shpFileLength())) +elif headerType == 'shx': + f.write(pack('>i', ((100 + (len(self._shapes) * 8)) // 2)))
+f.write(pack("<2i", 1000, self.shapeType))
+if self.shapeType != 0: + try: + f.write(pack("<4d", *self.bbox())) + except error: + raise ShapefileException("Failed to write shapefile bounding box. Floats required.") +else: + f.write(pack("<4d", 0,0,0,0))
+z = self.zbox()
+m = self.mbox() +try: + f.write(pack("<4d", z[0], z[1], m[0], m[1])) +except error: + raise ShapefileException("Failed to write shapefile elevation and measure values. Floats required.")
+__dbfHeader(self): +Writes the dbf header and field descriptors.
+ def null(self):Remove deletion flag placeholder from fields
+ self._shapes.append(_Shape(NULL))Field descriptors
+ def point(self, x, y, z=0, m=0):Terminator
+ pointShape = _Shape(self.shapeType)
+ pointShape.points.append([x, y, z, m])
+ self._shapes.append(pointShape)Write the shp records
+ def line(self, parts=[], shapeType=POLYLINE): self.poly(parts, shapeType, [])Record number, Content length place holder
+ def poly(self, parts=[], shapeType=POLYGON, partTypes=[]):Shape Type
+ polyShape = _Shape(shapeType)
+ polyShape.parts = []
+ polyShape.points = []
+ for part in parts:
+ polyShape.parts.append(len(polyShape.points))
+ for point in part:All shape types capable of having a bounding box
+ if not isinstance(point, list):
+ point = list(point)Shape types with parts
+ while len(point) < 4:
+ point.append(0)
+ polyShape.points.append(point)
+ if polyShape.shapeType == 31:
+ if not partTypes:
+ for part in parts:
+ partTypes.append(polyShape.shapeType)
+ polyShape.partTypes = partTypes
+ self._shapes.append(polyShape)Number of parts
+ def field(self, name, fieldType="C", size="50", decimal=0):Shape types with multiple points per record
+ self.fields.append((name, fieldType, size, decimal))Number of points
+ def record(self, *recordList, **recordDict):Write part indexes
+ return self._shapesPart types for Multipatch (31)
+ def saveShp(self, target):Write points for multiple-point records
+ if not hasattr(target, "write"):
+ target = os.path.splitext(target)[0] + '.shp'
+ if not self.shapeType:
+ self.shapeType = self._shapes[0].shapeType
+ self.shp = self.__getFileObj(target)
+ self.__shapefileHeader(self.shp, headerType='shp')
+ self.__shpRecords()Write z extremes and values
+ def saveShx(self, target):Write m extremes and values
+ if not hasattr(target, "write"):
+ target = os.path.splitext(target)[0] + '.shx'
+ if not self.shapeType:
+ self.shapeType = self._shapes[0].shapeType
+ self.shx = self.__getFileObj(target)
+ self.__shapefileHeader(self.shx, headerType='shx')
+ self.__shxRecords()Write a single point
+ def saveDbf(self, target):Write a single Z value
+ if not hasattr(target, "write"):
+ target = os.path.splitext(target)[0] + '.dbf'
+ self.dbf = self.__getFileObj(target)
+ self.__dbfHeader()
+ self.__dbfRecords()Write a single M value
+ def save(self, target=None, shp=None, shx=None, dbf=None):Finalize record length as 16-bit words
+ passstart - 4 bytes is the content length field
+ def delete(self, shape=None, part=None, point=None):Writes the shx records.
+ you to update a specific point by shape, part, point of any
+ shape type."""
+#DIVIDER
+ if shape and part and point:
+ try: self._shapes[shape]
+ except IndexError: self._shapes.append([])
+ try: self._shapes[shape][part]
+ except IndexError: self._shapes[shape].append([])
+ try: self._shapes[shape][part][point]
+ except IndexError: self._shapes[shape][part].append([])
+ p = self._shapes[shape][part][point]
+ if x: p[0] = x
+ if y: p[1] = y
+ if z: p[2] = z
+ if m: p[3] = m
+ self._shapes[shape][part][point] = p
+#DIVIDER
+ elif shape and part and not point:
+ try: self._shapes[shape]
+ except IndexError: self._shapes.append([])
+ try: self._shapes[shape][part]
+ except IndexError: self._shapes[shape].append([])
+ points = self._shapes[shape][part]
+ for i in range(len(points)):
+ p = points[i]
+ if x: p[0] = x
+ if y: p[1] = y
+ if z: p[2] = z
+ if m: p[3] = m
+ self._shapes[shape][part][i] = p
+#DIVIDER
+ elif shape and not part and not point:
+ try: self._shapes[shape]
+ except IndexError: self._shapes.append([])
+#DIVIDER
+ if addr:
+ shape, part, point = addr
+ self._shapes[shape][part][point] = [x, y, z, m]
+ else:
+ Writer.point(self, x, y, z, m)
+ if self.autoBalance:
+ self.balance()
+#DIVIDER
+ def validate(self):
+#DIVIDER
+ on which type of record was created to make sure all three files
+ are in synch."""
+ if len(self.records) > len(self._shapes):
+ self.null()
+ elif len(self.records) < len(self._shapes):
+ self.record() def __fieldNorm(self, fieldName):Writes the dbf records.
+ Doctests are contained in the module 'pyshp_usage.py'. This library was developed
+ using Python 2.3. Python 2.4 and above have some excellent improvements in the built-in
+ testing libraries but for now unit testing is done using what's available in
+ 2.3.Creates a null shape.
+Creates a point shape.
+Creates a line shape. This method is just a convienience method
+which wraps 'poly()'.
+Creates a shape that has multiple collections of points (parts)
+including lines, polygons, and even multipoint shapes. If no shape type +is specified it defaults to 'polygon'. If no part types are specified +(which they normally won't be) then all parts default to the shape type.
+Ensure point is list
+Make sure point has z and m values
+Adds a dbf field descriptor to the shapefile.
+Creates a dbf attribute record. You can submit either a sequence of
+field values or keyword arguments of field names and values. Before +adding records you must add fields for the record values using the +fields() method. If the record values exceed the number of fields the +extra ones won't be added. In the case of using keyword arguments to specify +field/value pairs only fields matching the already registered fields +will be added.""" +record = [] +fieldCount = len(self.fields)
+if self.fields[0][0].startswith("Deletion"): fieldCount -= 1 +if recordList: + [record.append(recordList[i]) for i in range(fieldCount)] +elif recordDict: + for field in self.fields: + if field[0] in recordDict: + val = recordDict[field[0]] + if val: + record.append(val) + else: + record.append("") +if record: + self.records.append(record)
+shape(self, i): +return self._shapes[i]
+shapes(self): +Return the current list of shapes.
+Save an shp file.
+Save an shx file.
+Save a dbf file.
+Save the shapefile data to three files or
+three file-like objects. SHP and DBF files can also +be written exclusively using saveShp, saveShx, and saveDbf respectively."""
+if shp: + self.saveShp(shp) +if shx: + self.saveShx(shx) +if dbf: + self.saveDbf(dbf) +elif target: + self.saveShp(target) + self.shp.close() + self.saveShx(target) + self.shx.close() + self.saveDbf(target) + self.dbf.close()
+itor(Writer): +init(self, shapefile=None, shapeType=POINT, autoBalance=1): +self.autoBalance = autoBalance +if not shapefile: + Writer.init(self, shapeType) +elif is_string(shapefile): + base = os.path.splitext(shapefile)[0] + if os.path.isfile("%s.shp" % base): + r = Reader(base) + Writer.init(self, r.shapeType) + self._shapes = r.shapes() + self.fields = r.fields + self.records = r.records()
+select(self, expr): +Select one or more shapes (to be implemented) +TODO: Implement expressions to select shapes.
+Deletes the specified part of any shape by specifying a shape
+number, part number, or point number."""
+if shape and part and point: + del self._shapes[shape][part][point]
+elif shape and part and not point: + del self._shapes[shape][part]
+elif shape and not part and not point: + del self._shapes[shape]
+elif not shape and not part and point: + for s in self._shapes: + if s.shapeType == 1: + del self._shapes[point] + else: + for part in s.parts: + del s[part][point]
+elif not shape and part and point: + for s in self._shapes: + del s[part][point]
+elif not shape and part and not point: + for s in self._shapes: + del s[part]
+point(self, x=None, y=None, z=None, m=None, shape=None, part=None, point=None, addr=None): +Creates/updates a point shape. The arguments allows
+shape, part, point
+shape, part
+shape
+point +part
+An optional method to try and validate the shapefile
+as much as possible before writing it (not implemented)."""
+pass
+balance(self): +Adds a corresponding empty attribute or null geometry record depending
+Normalizes a dbf field name to fit within the spec and the
+expectations of certain ESRI software.""" +if len(fieldName) > 11: fieldName = fieldName[:11] +fieldName = fieldName.upper() +fieldName.replace(' ', '_')
+Testing +(): +rt doctest +est.NORMALIZE_WHITESPACE = 1 +est.testfile("README.txt", verbose=1)
+e == "main__":
+test()
+from layersource import LayerSource
+from kartograph.errors import *
+from kartograph.geometry import BBox, create_featurethis class handles shapefile layers
+class ShapefileLayer(LayerSource):initialize shapefile reader
+ def __init__(self, src): import shapefile
+ if isinstance(src, unicode):
+ src = src.encode('ascii', 'ignore')
+ self.shpSrc = src
+ self.sr = shapefile.Reader(src)
+ self.recs = []
+ self.shapes = {}
+ self.load_records()load shapefile records into memory. note that only the records are loaded and not the shapes.
+ def load_records(self): self.recs = self.sr.records()
+ self.attributes = []
+ for a in self.sr.fields[1:]:
+ self.attributes.append(a[0])
+ i = 0
+ self.attrIndex = {}
+ for attr in self.attributes:
+ self.attrIndex[attr] = i
+ i += 1returns a shape of this shapefile. if requested for the first time, the shape is loaded from shapefile (slow)
+ def get_shape(self, i): if i in self.shapes: # check cache
+ shp = self.shapes[i]
+ else: # load shape from shapefile
+ shp = self.shapes[i] = self.sr.shapeRecord(i).shape
+ return shpreturns a list of features matching to the attr -> value pair
+ def get_features(self, attr=None, filter=None, bbox=None, verbose=False, ignore_holes=False, min_area=False, charset='utf-8'): res = []
+ try_encodings = ('utf-8', 'latin-1', 'iso-8859-2')
+ tried_encodings = [charset]
+ if bbox is not None and not isinstance(bbox, BBox):
+ bbox = BBox(bbox[2] - bbox[0], bbox[3] - bbox[1], bbox[0], bbox[1])
+ ignored = 0
+ for i in range(0, len(self.recs)):
+ drec = {}
+ for j in range(len(self.attributes)):
+ drec[self.attributes[j]] = self.recs[i][j]
+ if filter is None or filter(drec):
+ props = {}
+ for j in range(len(self.attributes)):
+ val = self.recs[i][j]
+ if isinstance(val, str):
+ try:
+ val = val.decode(charset)
+ except:
+ print 'warning: could not decode "%s" to %s' % (val, charset)
+ next_guess = False
+ for enc in try_encodings:
+ if enc not in tried_encodings:
+ next_guess = enc
+ tried_encodings.append(enc)
+ break
+ if next_guess:
+ print 'trying %s now..' % next_guess
+ charset = next_guess
+ j -= 1
+ continue
+ else:
+ raise KartographError('having problems to decode the input data "%s"' % val)
+ if isinstance(val, (str, unicode)):
+ val = val.strip()
+ props[self.attributes[j]] = val
+
+ shp = self.get_shape(i)
+
+ geom = shape2geometry(shp, ignore_holes=ignore_holes, min_area=min_area, bbox=bbox)
+ if geom is None:
+ ignored += 1
+ continue
+
+ feature = create_feature(geom, props)
+ res.append(feature)
+ if bbox is not None and ignored > 0 and verbose:
+ print "-ignoring %d shapes (not in bounds %s )" % (ignored, bbox)
+ return resdef shape2geometry(shp, ignore_holes=False, min_area=False, bbox=False):
+ if bbox:
+ sbbox = BBox(left=shp.bbox[0], top=shp.bbox[1], width=shp.bbox[2] - shp.bbox[0], height=shp.bbox[3] - shp.bbox[1])
+ if not bbox.intersects(sbbox):ignore the shape if it's not within the bbox
+ return None
+
+ if shp.shapeType in (5, 15): # multi-polygon
+ geom = shape2polygon(shp, ignore_holes=ignore_holes, min_area=min_area)
+ elif shp.shapeType == 3: # line
+ geom = points2line(shp)
+ else:
+ raise KartographError('unknown shape type (%d) in shapefile %s' % (shp.shapeType, self.shpSrc))
+ return geomconverts a shapefile polygon to geometry.MultiPolygon
+def shape2polygon(shp, ignore_holes=False, min_area=False):from kartograph.geometry import MultiPolygon
+ from shapely.geometry import Polygon, MultiPolygon
+ from kartograph.geometry.utils import is_clockwise
+ parts = shp.parts[:]
+ parts.append(len(shp.points))
+ exteriors = []
+ holes = []
+ for j in range(len(parts) - 1):
+ pts = shp.points[parts[j]:parts[j + 1]]
+ if shp.shapeType == 15:remove z-coordinate from PolygonZ contours (not supported)
+ for k in range(len(pts)):
+ pts[k] = pts[k][:2]
+ cw = is_clockwise(pts)
+ if cw:
+ exteriors.append(pts)
+ else:
+ holes.append(pts)
+ if ignore_holes:
+ holes = None
+ if len(exteriors) == 1:
+ poly = Polygon(exteriors[0], holes)
+ elif len(exteriors) > 1:use multipolygon, but we need to assign the holes to the right +exteriors
+ from kartograph.geometry import BBox
+ used_holes = set()
+ polygons = []
+ for ext in exteriors:
+ bbox = BBox()
+ my_holes = []
+ for pt in ext:
+ bbox.update(pt)
+ for h in range(len(holes)):
+ if h not in used_holes:
+ hole = holes[h]
+ if bbox.check_point(hole[0]):this is a very weak test but it should be sufficient
+ used_holes.add(h)
+ my_holes.append(hole)
+ polygons.append(Polygon(ext, my_holes))
+ if min_area:compute maximum area
+ max_area = 0
+ for poly in polygons:
+ max_area = max(max_area, poly.area)filter out polygons that are below min_area * max_area
+ polygons = [poly for poly in polygons if poly.area >= min_area * max_area]
+ poly = MultiPolygon(polygons)
+ else:
+ raise KartographError('shapefile import failed - no outer polygon found')
+ return polyconverts a shapefile line to geometry.Line
+def points2line(shp): from kartograph.geometry import PolyLine
+ parts = shp.parts[:]
+ parts.append(len(shp.points))
+ lines = []
+ for j in range(len(parts) - 1):
+ pts = shp.points[parts[j]:parts[j + 1]]
+ lines.append(pts)
+ return PolyLine(lines)
+
+from maplayer import MapLayer
+from geometry.utils import geom_to_bbox
+from geometry import BBox, View
+from shapely.geometry.base import BaseGeometry
+from shapely.geometry import Polygon
+from proj import projections
+from filter import filter_record
+from errors import KartographError
+
+_verbose = Falseclass Map(object): def __init__(me, options, layerCache, verbose=False, format='svg', src_encoding=None):
+ me.options = options
+ me._verbose = verbose
+ me.format = format
+ me.layers = []
+ me.layersById = {}
+ me._bounds_polygons_cache = False
+ me._unprojected_bounds = None
+ if not src_encoding:
+ src_encoding = 'utf-8'
+ me._source_encoding = src_encoding
+
+ for layer_cfg in options['layers']:
+ layer_id = layer_cfg['id']
+ layer = MapLayer(layer_id, layer_cfg, me, layerCache)
+ me.layers.append(layer)
+ me.layersById[layer_id] = layer
+
+ me.proj = me._init_projection()
+ me.bounds_poly = me._init_bounds()
+ me.view = me._get_view()
+ me.view_poly = me._init_view_poly()get features
+ for layer in me.layers:
+ layer.get_features()_debug_show_features(layerFeatures[id], 'original')
+ me._join_layers()_debug_show_features(layerFeatures[id], 'joined')
+ if options['export']['crop-to-view'] and format != 'kml':
+ me._crop_layers_to_view()_debug_show_features(layerFeatures[id], 'cropped to view')
+ me._simplify_layers()_debug_show_features(layerFeatures[id], 'simplified') +self.crop_layers(layers, layerOpts, layerFeatures)
+ me._subtract_layers()instantiates the map projection
+ def _init_projection(self): if self.format in ('kml', 'json'):
+ return projections['ll']() # use no projection for KML
+
+ map_center = self.__get_map_center()
+ opts = self.options
+ projC = projections[opts['proj']['id']]
+ p_opts = {}
+ for prop in opts['proj']:
+ if prop != "id":
+ p_opts[prop] = opts['proj'][prop]
+ if prop == "lon0" and p_opts[prop] == "auto":
+ p_opts[prop] = map_center[0]
+ elif prop == "lat0" and p_opts[prop] == "auto":
+ p_opts[prop] = map_center[1]
+ return projC(**p_opts)used by _init_projection() to determine the center of the +map projection, depending on the bounds config
+ def __get_map_center(self): opts = self.options
+ mode = opts['bounds']['mode']
+ data = opts['bounds']['data']
+
+ lon0 = 0
+
+ if mode == 'bbox':
+ lon0 = data[0] + 0.5 * (data[2] - data[0])
+ lat0 = data[1] + 0.5 * (data[3] - data[1])
+
+ elif mode[:5] == 'point':
+ lon0 = 0
+ lat0 = 0
+ m = 1 / len(data)
+ for (lon, lat) in data:
+ lon0 += m * lon
+ lat0 += m * lat
+
+ elif mode[:4] == 'poly':
+ features = self._get_bounds_polygons()
+ if len(features) > 0:
+ if isinstance(features[0].geom, BaseGeometry):
+ (lon0, lat0) = features[0].geom.representative_point().coords[0]
+ else:
+ lon0 = 0
+ lat0 = 0
+ else:
+ print "unrecognized bound mode", mode
+ return (lon0, lat0)computes the (x,y) bounding box for the map, +given a specific projection
+ def _init_bounds(self): if self.format in ('kml', 'json'):
+ return None # no bounds needed for KML
+
+ from geometry.utils import bbox_to_polygon
+
+ opts = self.options
+ proj = self.proj
+ bnds = opts['bounds']
+ mode = bnds['mode'][:]
+ data = bnds['data']
+
+ if _verbose:
+ print 'using bounds mode', mode
+
+ if mode == "bbox": # catch special case bbox
+ sea = proj.bounding_geometry(data, projected=True)
+ sbbox = geom_to_bbox(sea)
+ sbbox.inflate(sbbox.width * bnds['padding'])
+ return bbox_to_polygon(sbbox)
+
+ bbox = BBox()
+
+ if mode[:5] == "point":
+ for lon, lat in data:
+ pt = proj.project(lon, lat)
+ bbox.update(pt)
+
+ if mode[:4] == "poly":
+ features = self._get_bounds_polygons()
+ ubbox = BBox()
+ if len(features) > 0:
+ for feature in features:
+ ubbox.join(geom_to_bbox(feature.geometry))
+ feature.project(proj)
+ fbbox = geom_to_bbox(feature.geometry, data["min-area"])
+ bbox.join(fbbox)
+ self._unprojected_bounds = ubbox
+ else:
+ raise KartographError('no features found for calculating the map bounds')
+ bbox.inflate(bbox.width * bnds['padding'])
+ return bbox_to_polygon(bbox)for bounds mode "polygons" this helper function +returns a list of all polygons that the map should +be cropped to
+ def _get_bounds_polygons(self): if self._bounds_polygons_cache:
+ return self._bounds_polygons_cache
+
+ opts = self.options
+ features = []
+ data = opts['bounds']['data']
+ id = data['layer']
+
+ if id not in self.layersById:
+ raise KartographError('layer not found "%s"' % id)
+ layer = self.layersById[id]
+
+ if layer.options['filter'] is False:
+ layerFilter = lambda a: True
+ else:
+ layerFilter = lambda rec: filter_record(layer.options['filter'], rec)
+
+ if data['filter']:
+ boundsFilter = lambda rec: filter_record(data['filter'], rec)
+ else:
+ boundsFilter = lambda a: True
+
+ filter = lambda rec: layerFilter(rec) and boundsFilter(rec)
+ features = layer.source.get_features(filter=filter, min_area=data["min-area"], charset=layer.options['charset'])remove features that are too small
+ if layer.options['filter-islands']:
+ features = [feature for feature in features if feature.geometry.area > layer.options['filter-islands']]
+
+ self._bounds_polygons_cache = features
+ return featuresreturns the output view
+ def _get_view(self): if self.format in ('kml', 'json'):
+ return View() # no view transformation needed for KML
+
+ self.src_bbox = bbox = geom_to_bbox(self.bounds_poly)
+ opts = self.options
+ exp = opts["export"]
+ w = exp["width"]
+ h = exp["height"]
+ ratio = exp["ratio"]
+
+ if ratio == "auto":
+ ratio = bbox.width / float(bbox.height)
+
+ if h == "auto":
+ h = w / ratio
+ elif w == "auto":
+ w = h * ratio
+ return View(bbox, w, h - 1)creates a polygon that represents the rectangular view bounds +used for cropping the geometries to not overlap the view
+ def _init_view_poly(self): if self.format in ('kml', 'json'):
+ return None # no view polygon needed for KML
+ w = self.view.width
+ h = self.view.height
+ return Polygon([(0, 0), (0, h), (w, h), (w, 0)])performs polygon simplification
+ def _simplify_layers(self): from simplify import create_point_store, simplify_lines
+
+ point_store = create_point_store() # create a new empty point storecompute topology for all layers
+ for layer in self.layers:
+ if layer.options['simplify'] is not False:
+ for feature in layer.features:
+ feature.compute_topology(point_store, layer.options['unify-precision'])break features into lines
+ for layer in self.layers:
+ if layer.options['simplify'] is not False:
+ for feature in layer.features:
+ feature.break_into_lines()simplify lines
+ total = 0
+ kept = 0
+ for layer in self.layers:
+ if layer.options['simplify'] is not False:
+ for feature in layer.features:
+ lines = feature.break_into_lines()
+ lines = simplify_lines(lines, layer.options['simplify']['method'], layer.options['simplify']['tolerance'])
+ for line in lines:
+ total += len(line)
+ for pt in line:
+ if not pt.deleted:
+ kept += 1
+ feature.restore_geometry(lines, layer.options['filter-islands'])
+ return (total, kept)cuts the layer features to the map view
+ def _crop_layers_to_view(self): for layer in self.layers:out = []
+ for feat in layer.features:
+ if not feat.geometry.is_valid:
+ passprint feat.geometry +_plot_geometry(feat.geometry)
+ feat.crop_to(self.view_poly)if not feat.is_empty(): + out.append(feat) +layer.features = out
+handles crop-to
+ def _crop_layers(self): for layer in self.layers:
+ if layer.options['crop-to'] is not False:
+ cropped_features = []
+ for tocrop in layer.features:
+ cbbox = geom_to_bbox(tocrop.geom)
+ crop_at_layer = layer.options['crop-to']
+ if crop_at_layer not in self.layers:
+ raise KartographError('you want to substract from layer "%s" which cannot be found' % crop_at_layer)
+ for crop_at in self.layersById[crop_at_layer].features:
+ if crop_at.geom.bbox().intersects(cbbox):
+ tocrop.crop_to(crop_at.geom)
+ cropped_features.append(tocrop)
+ layer.features = cropped_featureshandles subtract-from
+ def _subtract_layers(self): for layer in self.layers:
+ if layer.options['subtract-from'] is not False:
+ for feat in layer.features:
+ if feat.geom is None:
+ continue
+ cbbox = geom_to_bbox(feat.geom)
+ for subid in layer.options['subtract-from']:
+ if subid not in self.layers:
+ raise KartographError('you want to subtract from layer "%s" which cannot be found' % subid)
+ for sfeat in self.layersById[subid].features:
+ if sfeat.geom and geom_to_bbox(sfeat.geom).intersects(cbbox):
+ sfeat.subtract_geom(feat.geom)
+ layer.features = []joins features in layers
+ def _join_layers(self): from geometry.utils import join_features
+
+ for layer in self.layers:
+ if layer.options['join'] is not False:
+ unjoined = 0
+ join = layer.options['join']
+ groupBy = join['group-by']
+ groups = join['groups']
+ if not groups:auto populate groups
+ groups = {}
+ for feat in layer.features:
+ fid = feat.props[groupBy]
+ groups[fid] = [fid]
+
+ groupAs = join['group-as']
+ groupFeatures = {}
+ res = []
+ for feat in layer.features:
+ found_in_group = False
+ for g_id in groups:
+ if g_id not in groupFeatures:
+ groupFeatures[g_id] = []
+ if feat.props[groupBy] in groups[g_id] or str(feat.props[groupBy]) in groups[g_id]:
+ groupFeatures[g_id].append(feat)
+ found_in_group = True
+ break
+ if not found_in_group:
+ unjoined += 1
+ res.append(feat)print unjoined,'features were not joined'
+ for g_id in groups:
+ props = {}
+ for feat in groupFeatures[g_id]:
+ fprops = feat.props
+ for key in fprops:
+ if key not in props:
+ props[key] = fprops[key]
+ else:
+ if props[key] != fprops[key]:
+ props[key] = "---"
+
+ if groupAs is not False:
+ props[groupAs] = g_id
+ if g_id in groupFeatures:
+ res += join_features(groupFeatures[g_id], props)
+ layer.features = res
+
+from layersource import handle_layer_source
+from filter import filter_record
+
+
+_verbose = Falseclass MapLayer(object): def __init__(self, id, options, _map, cache):Store layer properties as instance properties
+ self.id = id
+ self.options = options
+ self.map = _map
+ self.cache = cacheMake sure that the layer id is unique within the map.
+ while self.id in self.map.layersById:
+ self.id += "_"Instantiate the layer source which will generate features from the source +geo data such as shapefiles or virtual sources such as graticule lines.
+ self.source = handle_layer_source(self.options, self.cache) def get_features(layer, filter=False, min_area=0): opts = layer.map.options
+ is_projected = FalseLet's see if theres a better bounding box than this..
+ bbox = [-180, -90, 180, 90]Use the clipping mode defined in the map configuration
+ if opts['bounds']['mode'] == "bbox":
+ bbox = opts['bounds']['data']The 'crop' property overrides the clipping settings
+ if 'crop' in opts['bounds'] and opts['bounds']['crop']:If crop is set to "auto", which is the default behaviour, Kartograph +will use the actual bounding geometry to compute the bounding box
+ if opts['bounds']['crop'] == "auto":
+ if layer.map._unprojected_bounds:
+ bbox = layer.map._unprojected_bounds
+ bbox.inflate(inflate=opts['bounds']['padding'] * 2)
+ else:
+ print 'could not compute bounding box for auto-cropping'
+ else:otherwise it will use the user defined bbox in the format +[minLon, minLat, maxLon, maxLat]
+ bbox = opts['bounds']['crop']If the layer has the "src" property, it is a regular map layer source, which +means that there's an exernal file that we load the geometry and meta data from.
+ if 'src' in layer.options:
+ if layer.options['filter'] is False:
+ filter = None
+ else:
+ filter = lambda rec: filter_record(layer.options['filter'], rec)Now we ask the layer source to generate the features that will be displayed +in the map.
+ features = layer.source.get_features(
+ filter=filter,
+ bbox=bbox,
+ ignore_holes='ignore-holes' in layer.options and layer.options['ignore-holes'],
+ charset=layer.options['charset']
+ )
+ if _verbose:
+ print 'loaded %d features from shapefile %s' % (len(features), layer.options['src'])In contrast to regular layers, the geometry for special (or virtual) layers is generated +by Kartograph itself, based on some properties defined in the layer config.
+ elif 'special' in layer.options:The graticule layer generates line features for longitudes and latitudes
+ if layer.options['special'] == "graticule":
+ lats = layer.options['latitudes']
+ lons = layer.options['longitudes']
+ features = layer.source.get_features(lats, lons, layer.map.proj, bbox=bbox)The "sea" layer generates a MultiPolygon that represents the entire boundary +of the map. Especially useful for non-cylindrical map projections.
+ elif layer.options['special'] == "sea":
+ features = layer.source.get_features(layer.map.proj)
+ is_projected = True
+
+ for feature in features:If the features are not projected yet, we project them now.
+ if not is_projected:
+ feature.project(layer.map.proj)Transform features to view coordinates.
+ feature.project_view(layer.map.view)Remove features that don't intersect our clipping polygon
+ if layer.map.view_poly:
+ features = [feature for feature in features
+ if feature.geometry and feature.geometry.intersects(layer.map.view_poly)]
+ layer.features = features
+
+API 2.0 +helper methods for validating options dictionary
+import os.path, proj, errors
+
+
+Error = errors.KartographOptionParseErrordef is_str(s):
+ return isinstance(s, (str, unicode))check out that the option dict is filled correctly
+def parse_options(opts):projection
+ parse_proj(opts)
+ parse_layers(opts)
+ parse_bounds(opts)
+ parse_export(opts)checks projections
+def parse_proj(opts): if 'proj' not in opts:
+ opts['proj'] = {}
+ prj = opts['proj']
+ if 'id' not in prj:
+ if 'bounds' not in opts:
+ prj['id'] = 'robinson'
+ else:
+ prj['id'] = 'laea'
+ if prj['id'] not in proj.projections:
+ raise Error('unknown projection')
+ prjClass = proj.projections[prj['id']]
+ for attr in prjClass.attributes():
+ if attr not in prj:
+ prj[attr] = "auto"
+ else:
+ if prj[attr] != "auto":
+ prj[attr] = float(prj[attr])def parse_layers(opts):
+ if 'layers' not in opts:
+ opts['layers'] = []
+ l_id = 0
+ g_id = 0
+ s_id = 0
+ for layer in opts['layers']:
+ if 'styles' not in layer:
+ layer['styles'] = {}
+ if 'src' not in layer and 'special' not in layer:
+ raise Error('you need to define the source for your layers')
+ if 'src' in layer:
+ if not os.path.exists(layer['src']):
+ raise Error('layer source not found: ' + layer['src'])
+ if 'id' not in layer:
+ layer['id'] = 'layer_' + str(l_id)
+ l_id += 1
+ if 'charset' not in layer:
+ layer['charset'] = 'utf-8'
+ elif 'special' in layer:
+ if layer['special'] == 'graticule':
+ if 'id' not in layer:
+ layer['id'] = 'graticule'
+ if g_id > 0:
+ layer['id'] += '_' + str(g_id)
+ g_id += 1
+ if 'fill' not in layer['styles']:
+ layer['styles']['fill'] = 'None'
+ parse_layer_graticule(layer)
+ elif layer['special'] == 'sea':
+ if 'id' not in layer:
+ layer['id'] = 'sea'
+ if s_id > 0:
+ layer['id'] += '_' + str(s_id)
+ s_id += 1
+
+ parse_layer_attributes(layer)
+ parse_layer_filter(layer)
+ parse_layer_join(layer)
+ parse_layer_simplify(layer)
+ parse_layer_subtract(layer)
+ parse_layer_cropping(layer)def parse_layer_attributes(layer):
+ if 'attributes' not in layer:
+ layer['attributes'] = []
+ return
+ attrs = []
+ for attr in layer['attributes']:
+ if is_str(attr):
+ if isinstance(layer['attributes'], list): # ["ISO_A3", "FIPS"]
+ attrs.append({'src': attr, 'tgt': attr})
+ elif isinstance(layer['attributes'], dict): # { "ISO_A3": "iso" }
+ attrs.append({'src': attr, 'tgt': layer['attributes'][attr]})
+ elif isinstance(attr, dict) and 'src' in attr and 'tgt' in attr:
+ attrs.append(attr)
+ layer['attributes'] = attrsdef parse_layer_filter(layer):
+ if 'filter' not in layer:
+ layer['filter'] = False
+ return
+ return # todo: check valid filter syntax (recursivly, place code in filter.py)
+ filter = layer['filter']
+ if 'type' not in filter:
+ filter['type'] = 'include'
+ if 'attribute' not in filter:
+ raise Error('layer filter must define an attribute to filter on')
+ if 'equals' in filter:
+ if isinstance(filter['equals'], (str, unicode, int, float)):
+ filter['equals'] = [filter['equals']]
+ elif 'greater-than' in filter:
+ try:
+ filter['greater-than'] = float(filter['greater-than'])
+ except ValueError:
+ raise Error('could not convert filter value "greater-than" to float')
+ elif 'less-than' in filter:
+ try:
+ filter['less-than'] = float(filter['less-than'])
+ except ValueError:
+ raise Error('could not convert filter value "less-than" to float')
+ else:
+ raise Error('you must define either "equals", "greater-than" or "less-than" in the filter')def parse_layer_join(layer):
+ if 'join' not in layer:
+ layer['join'] = False
+ return
+ if layer['join'] is False:
+ return
+
+ join = layer['join']
+ if 'group-by' not in join:
+ raise Error('missing attribute "group-by": you need to specify an attribute by which the features should be joined.')
+ if 'groups' not in join:
+ join['groups'] = None
+ if 'group-as' not in join:
+ join['group-as'] = Falsedef parse_layer_simplify(layer):
+ if 'unify-precision' not in layer:
+ layer['unify-precision'] = None
+ if 'simplify' not in layer:
+ layer['simplify'] = False
+ return
+ if layer['simplify'] is False:
+ return
+ if isinstance(layer['simplify'], (int, float, str, unicode)):default to visvalingam-whyatt
+ layer['simplify'] = {"method": "visvalingam-whyatt", "tolerance": float(layer['simplify'])}
+ try:
+ layer['simplify']['tolerance'] = float(layer['simplify']['tolerance'])
+ except ValueError:
+ raise Error('could not convert simplification amount to float')
+ if 'filter-islands' not in layer:
+ layer['filter-islands'] = Falsedef parse_layer_subtract(layer):
+ if 'subtract-from' not in layer:
+ layer['subtract-from'] = False
+ return
+ if isinstance(layer['subtract-from'], (str, unicode)):
+ layer['subtract-from'] = [layer['subtract-from']]def parse_layer_cropping(layer):
+ if 'crop-to' not in layer:
+ layer['crop-to'] = False
+ returndef parse_layer_graticule(layer):
+ if 'latitudes' not in layer:
+ layer['latitudes'] = []
+ elif isinstance(layer['latitudes'], (int, float)):
+ step = layer['latitudes']
+ layer['latitudes'] = [0]
+ for lat in _xfrange(step, 90, step):
+ layer['latitudes'] += [lat, -lat]
+
+ if 'longitudes' not in layer:
+ layer['longitudes'] = []
+ elif isinstance(layer['longitudes'], (int, float)):
+ step = layer['longitudes']
+ layer['longitudes'] = [0]
+ for lon in _xfrange(step, 181, step):
+ if lon == 180:
+ p = [lon]
+ else:
+ p = [lon, -lon]
+ layer['longitudes'] += pdef _xfrange(start, stop, step):
+ while (step > 0 and start < stop) or (step < 0 and start > step):
+ yield start
+ start += stepdef parse_bounds(opts):
+ if 'bounds' not in opts:
+ opts['bounds'] = {}return
+ bounds = opts['bounds']
+ if 'mode' not in bounds:
+ bounds['mode'] = 'bbox'
+
+ if 'data' not in bounds:
+ bounds['data'] = [-180, -90, 180, 90]
+ bounds['mode'] = 'bbox'
+
+ mode = bounds['mode']
+ data = bounds['data']
+
+ if 'crop' not in bounds:
+ bounds['crop'] = 'auto'
+
+ if "padding" not in bounds:
+ bounds["padding"] = 0
+
+ if mode == "bbox":
+ try:
+ if len(data) == 4:
+ for i in range(0, 4):
+ data[i] = float(data[i])
+ else:
+ raise Error('bounds mode bbox requires array with exactly 4 values [lon0,lat0,lon1,lat]')
+ except Error:
+ raise
+ except:
+ raise Error('bounds mode bbox requires array with exactly 4 values [lon0,lat0,lon1,lat]')
+ elif mode == "points":
+ try:
+ for i in range(0, len(data)):
+ pt = data[i]
+ if len(pt) == 2:
+ pt = map(float, pt)
+ else:
+ raise Error('bounds mode points requires array with (lon,lat) tuples')
+ except Error:
+ raise
+ except:
+ raise Error('bounds mode points requires array with (lon,lat) tuples')
+ elif mode in ("polygons", "polygon"):
+ bounds['mode'] = mode = "polygons"
+ if "layer" not in data or not is_str(data["layer"]):
+ raise Error('you must specify a layer for bounds mode ' + mode)
+ if "filter" not in data:
+ data["filter"] = False
+ if "attribute" not in data or not is_str(data["attribute"]):
+ data["attribute"] = None
+ if "values" not in data:
+ if data["attribute"] is None:
+ data["values"] = None
+ else:
+ raise Error('you must specify a list of values to match in bounds mode ' + mode)
+ if is_str(data["values"]):
+ data["values"] = [data["values"]]
+ if "min-area" in data:
+ try:
+ data["min-area"] = float(data["min-area"])
+ except:
+ raise Error('min_area must be an integer or float')
+ else:
+ data['min-area'] = 0def parse_export(opts):
+ if "export" not in opts:
+ opts["export"] = {}
+ exp = opts["export"]
+ if "width" not in exp and "height" not in exp:
+ exp["width"] = 1000
+ exp["height"] = "auto"
+ elif "height" not in exp:
+ exp["height"] = "auto"
+ elif "width" not in exp:
+ exp["width"] = "auto"
+
+ if "ratio" not in exp:
+ exp["ratio"] = "auto"
+ if "round" not in exp:
+ exp["round"] = False
+ else:
+ exp["round"] = int(exp["round"])
+ if "crop-to-view" not in exp:
+ exp['crop-to-view'] = True
+
+kartograph - a svg mapping library
+Copyright (C) 2011 Gregor Aisch
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see <http://www.gnu.org/licenses/>.
+projections = dict()
+
+from base import Proj
+from cylindrical import *
+
+projections['lonlat'] = Equirectangular
+projections['cea'] = CEA
+projections['gallpeters'] = GallPeters
+projections['hobodyer'] = HoboDyer
+projections['behrmann'] = Behrmann
+projections['balthasart'] = Balthasart
+projections['mercator'] = Mercator
+projections['ll'] = LonLat
+
+from pseudocylindrical import *
+
+projections['naturalearth'] = NaturalEarth
+projections['robinson'] = Robinson
+projections['eckert4'] = EckertIV
+projections['sinusoidal'] = Sinusoidal
+projections['mollweide'] = Mollweide
+projections['wagner4'] = WagnerIV
+projections['wagner5'] = WagnerV
+projections['loximuthal'] = Loximuthal
+projections['canters1'] = CantersModifiedSinusoidalI
+projections['goodehomolosine'] = GoodeHomolosine
+projections['hatano'] = Hatano
+projections['aitoff'] = Aitoff
+projections['winkel3'] = Winkel3
+projections['nicolosi'] = Nicolosi
+
+from azimuthal import *
+
+projections['ortho'] = Orthographic
+projections['laea'] = LAEA
+projections['stereo'] = Stereographic
+projections['satellite'] = Satellite
+projections['eda'] = EquidistantAzimuthal
+projections['aitoff'] = Aitoff
+
+from conic import *
+
+projections['lcc'] = LCCclass Proj4(Proj): def __init__(self, projstr):
+ import pyproj
+ self.proj = pyproj.Proj(projstr) def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ x, y = self.proj(lon, lat)
+ return (x, y * -1)class LCC__(Proj4): def __init__(self, lat0=0, lon0=0, lat1=28, lat2=30):
+ Proj4.__init__(self, '+proj=lcc +lat_0=%f +lon_0=%f +lat_1=%f +lat_2=%f' % (lat0, lon0, lat1, lat2)) def _visible(self, lon, lat):
+ return True def _truncate(self, x, y):
+ return (x, y)
+
+
+for pjname in projections:
+ projections[pjname].name = pjname
+
+
+if __name__ == '__main__':
+ import syssome class testing +p = LAEA(52.0,10.0) +x,y = p.project(50,5) +assert (round(x,2),round(y,2)) == (3962799.45, -2999718.85), 'LAEA proj error'
+ from kartograph.geometry import BBox
+
+ print Proj.fromXML(Robinson(lat0=3, lon0=4).toXML(), projections)
+
+ Robinson(lat0=3, lon0=4)
+
+ for pj in projections:
+ Proj = projections[pj]
+ bbox = BBox()
+ try:
+ proj = Proj(lon0=60)
+ print proj.project(0, 0)
+ print proj.world_bounds(bbox)
+ print proj.toXML()
+ except:
+ print 'Error', pj
+ print sys.exc_info()[0]
+ raise
+
+kartograph - a svg mapping library
+Copyright (C) 2011 Gregor Aisch
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see <http://www.gnu.org/licenses/>.
+import math
+from kartograph.errors import KartographError
+from shapely.geometry import Polygon, LineString, Point, MultiPolygon, MultiLineString, MultiPointbase class for projections
+class Proj(object): HALFPI = math.pi * .5
+ QUARTERPI = math.pi * .25
+
+ minLat = -90
+ maxLat = 90
+ minLon = -180
+ maxLon = 180 def _shift_polygon(self, polygon):
+ return [polygon] # no shifting def plot(self, geometry):
+ geometries = hasattr(geometry, 'geoms') and geometry.geoms or [geometry]
+ res = []at first shift polygons +shifted = [] +for geom in geometries: + if isinstance(geom, Polygon): + shifted += self._shift_polygon(geom) + else: + shifted += [geom]
+ for geom in geometries:
+ if isinstance(geom, Polygon):
+ res += self.plot_polygon(geom)
+ elif isinstance(geom, LineString):
+ rings = self.plot_linear_ring(geom.coords)
+ res += map(LineString, rings)
+ elif isinstance(geom, Point):
+ if self._visible(geom.x, geom.y):
+ x, y = self.project(geom.x, geom.y)
+ res.append(Point(x, y))
+ else:
+ passraise KartographError('proj.plot(): unknown geometry type %s' % geom)
+ if len(res) > 0:
+ if isinstance(res[0], Polygon):
+ if len(res) > 1:
+ return MultiPolygon(res)
+ else:
+ return res[0]
+ elif isinstance(res[0], LineString):
+ if len(res) > 1:
+ return MultiLineString(res)
+ else:
+ return LineString(res[0])
+ else:
+ if len(res) > 1:
+ return MultiPoint(res)
+ else:
+ return Point(res[0][0], res[0][1]) def plot_polygon(self, polygon):
+ ext = self.plot_linear_ring(polygon.exterior, truncate=True)
+ if len(ext) == 1:
+ pts_int = []
+ for interior in polygon.interiors:
+ pts_int += self.plot_linear_ring(interior, truncate=True)
+ return [Polygon(ext[0], pts_int)]
+ elif len(ext) == 0:
+ return []
+ else:
+ raise KartographError('unhandled case: exterior is split into multiple rings') def plot_linear_ring(self, ring, truncate=False):
+ ignore = True
+ points = []
+ for (lon, lat) in ring.coords:
+ vis = self._visible(lon, lat)
+ if vis:
+ ignore = False
+ x, y = self.project(lon, lat)
+ if not vis and truncate:
+ points.append(self._truncate(x, y))
+ else:
+ points.append((x, y))
+ if ignore:
+ return []
+ return [points] def ll(self, lon, lat):
+ return (lon, lat) def project(self, lon, lat):
+ assert False, 'Proj is an abstract class' def project_inverse(self, x, y):
+ assert False, 'inverse projection is not supporte by %s' % self.name def _visible(self, lon, lat):
+ assert False, 'Proj is an abstract class' def _truncate(self, x, y):
+ assert False, 'truncation is not implemented' def world_bounds(self, bbox, llbbox=(-180, -90, 180, 90)):
+ sea = self.sea_shape(llbbox)
+ for x, y in sea[0]:
+ bbox.update((x, y))
+ return bboxreturns a WGS84 polygon that represents the limits of this projection +points that lie outside this polygon will not be plotted +this polygon will also be used to render the sea layer in world maps
+defaults to full WGS84 range
+ def bounding_geometry(self, llbbox=(-180, -90, 180, 90), projected=False): from shapely.geometry import Polygon
+ sea = []
+
+ minLon = llbbox[0]
+ maxLon = llbbox[2]
+ minLat = max(self.minLat, llbbox[1])
+ maxLat = min(self.maxLat, llbbox[3]) def xfrange(start, stop, step):
+ if stop > start:
+ while start < stop:
+ yield start
+ start += step
+ else:
+ while stop < start:
+ yield start
+ start -= step
+
+ lat_step = abs((maxLat - minLat) / 180.0)
+ lon_step = abs((maxLon - minLon) / 360.0)
+
+ for lat in xfrange(minLat, maxLat, lat_step):
+ sea.append((minLon, lat))
+ for lon in xfrange(minLon, maxLon, lon_step):
+ sea.append((lon, maxLat))
+ for lat in xfrange(maxLat, minLat, lat_step):
+ sea.append((maxLon, lat))
+ for lon in xfrange(maxLon, minLon, lon_step):
+ sea.append((lon, minLat))
+
+ if projected:
+ sea = [self.project(lon, lat) for (lon, lat) in sea]
+
+ return Polygon(sea) def __str__(self):
+ return 'Proj(' + self.name + ')' def attrs(self):
+ return dict(id=self.name)returns array of attribute names of this projection
+ @staticmethod
+ def attributes(): return [] @staticmethod
+ def fromXML(xml, projections):
+ id = xml['id']
+ if id in projections:
+ ProjClass = projections[id]
+ args = {}
+ for (prop, val) in xml:
+ if prop[0] != "id":
+ args[prop[0]] = float(val)
+ return ProjClass(**args)
+ raise Exception("could not restore projection from xml")
+
+kartograph - a svg mapping library
+Copyright (C) 2011 Gregor Aisch
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see <http://www.gnu.org/licenses/>.
+from base import Proj
+import math
+from math import radians as radclass Conic(Proj): def __init__(self, lat0=0, lon0=0, lat1=0, lat2=0):
+ self.lat0 = lat0
+ self.phi0 = rad(lat0)
+ self.lon0 = lon0
+ self.lam0 = rad(lon0)
+ self.lat1 = lat1
+ self.phi1 = rad(lat1)
+ self.lat2 = lat2
+ self.phi2 = rad(lat2)
+
+ if lon0 != 0.0:
+ self.bounds = self.bounding_geometry() def _visible(self, lon, lat):
+ return True def _truncate(self, x, y):
+ return (x, y) def attrs(self):
+ p = super(Conic, self).attrs()
+ p['lon0'] = self.lon0
+ p['lat0'] = self.lat0
+ p['lat1'] = self.lat1
+ p['lat2'] = self.lat2
+ return pshifts a polygon according to the origin longitude
+ def _shift_polygon(self, polygon): if self.lon0 == 0.0:
+ return [polygon] # no need to shift anything
+
+ from shapely.geometry import Polygonwe need to split and join some polygons
+ poly_coords = []
+ holes = []
+ for (lon, lat) in polygon.exterior.coords:
+ poly_coords.append((lon - self.lon0, lat))
+ for hole in polygon.interiors:
+ hole_coords = []
+ for (lon, lat) in hole.coords:
+ hole_coords.append((lon - self.lon0, lat))
+ holes.append(hole_coords)
+ poly = Polygon(poly_coords, holes)
+
+ polygons = []print "shifted polygons", (time.time() - start) +start = time.time()
+ try:
+ p_in = poly.intersection(self.bounds)
+ polygons += hasattr(p_in, 'geoms') and p_in.geoms or [p_in]
+ except:
+ passprint "computed polygons inside bounds", (time.time() - start) +start = time.time()
+ try:
+ p_out = poly.symmetric_difference(self.bounds)
+ out_geoms = hasattr(p_out, 'geoms') and p_out.geoms or [p_out]
+ except:
+ out_geoms = []
+ passprint "computed polygons outside bounds", (time.time() - start) +start = time.time()
+ for polygon in out_geoms:
+ ext_pts = []
+ int_pts = []
+ s = 0 # at first we compute the avg longitude
+ c = 0
+ for (lon, lat) in polygon.exterior.coords:
+ s += lon
+ c += 1
+ left = s / float(c) < -180 # and use it to decide where to shift the polygon
+ for (lon, lat) in polygon.exterior.coords:
+ ext_pts.append((lon + (-360, 360)[left], lat))
+ for interior in polygon.interiors:
+ pts = []
+ for (lon, lat) in interior.coords:
+ pts.append((lon + (-360, 360)[left], lat))
+ int_pts.append(pts)
+ polygons.append(Polygon(ext_pts, int_pts))print "shifted outside polygons to inside", (time.time() - start)
+ return polygonsLambert Conformal Conic Projection (spherical)
+ @staticmethod
+ def attributes():
+ return ['lon0', 'lat0', 'lat1', 'lat2']
+
+
+class LCC(Conic): def __init__(self, lat0=0, lon0=0, lat1=30, lat2=50):
+ from math import sin, cos, tan, pow, log
+ self.minLat = -60
+ self.maxLat = 85
+ Conic.__init__(self, lat0=lat0, lon0=lon0, lat1=lat1, lat2=lat2)
+ self.n = n = sin(self.phi1)
+ cosphi = cos(self.phi1)
+ secant = abs(self.phi1 - self.phi2) >= 1e-10
+ if secant:
+ n = log(cosphi / cos(self.phi2)) / log(tan(self.QUARTERPI + .5 * self.phi2) / tan(self.QUARTERPI + .5 * self.phi1))
+ self.c = c = cosphi * pow(tan(self.QUARTERPI + .5 * self.phi1), n) / n
+ if abs(abs(self.phi0) - self.HALFPI) < 1e-10:
+ self.rho0 = 0.
+ else:
+ self.rho0 = c * pow(tan(self.QUARTERPI + .5 * self.phi0), -n) def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ phi = rad(lat)
+ lam = rad(lon)
+ n = self.n
+ if abs(abs(phi) - self.HALFPI) < 1e-10:
+ rho = 0.0
+ else:
+ rho = self.c * math.pow(math.tan(self.QUARTERPI + 0.5 * phi), -n)
+ lam_ = (lam - self.lam0) * n
+ x = 1000 * rho * math.sin(lam_)
+ y = 1000 * (self.rho0 - rho * math.cos(lam_))
+
+ return (x, y * -1)
+
+kartograph - a svg mapping library
+Copyright (C) 2011 Gregor Aisch
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see <http://www.gnu.org/licenses/>.
+from base import Proj
+import math
+from math import radians as radclass Cylindrical(Proj): def __init__(self, lon0=0.0, flip=0):
+ self.flip = flip
+ self.lon0 = lon0
+ self.bounds = self.bounding_geometry()shifts a polygon according to the origin longitude
+ def _shift_polygon(self, polygon): if self.lon0 == 0.0:
+ return [polygon] # no need to shift anything
+
+ from shapely.geometry import Polygonwe need to split and join some polygons
+ poly_coords = []
+ holes = []
+ for (lon, lat) in polygon.exterior.coords:
+ poly_coords.append((lon - self.lon0, lat))
+ for hole in polygon.interiors:
+ hole_coords = []
+ for (lon, lat) in hole.coords:
+ hole_coords.append((lon - self.lon0, lat))
+ holes.append(hole_coords)
+ poly = Polygon(poly_coords, holes)
+
+ polygons = []print "shifted polygons", (time.time() - start) +start = time.time()
+ try:
+ p_in = poly.intersection(self.bounds)
+ polygons += hasattr(p_in, 'geoms') and p_in.geoms or [p_in]
+ except:
+ passprint "computed polygons inside bounds", (time.time() - start) +start = time.time()
+ try:
+ p_out = poly.symmetric_difference(self.bounds)
+ out_geoms = hasattr(p_out, 'geoms') and p_out.geoms or [p_out]
+ except:
+ out_geoms = []
+ passprint "computed polygons outside bounds", (time.time() - start) +start = time.time()
+ for polygon in out_geoms:
+ ext_pts = []
+ int_pts = []
+ s = 0 # at first we compute the avg longitude
+ c = 0
+ for (lon, lat) in polygon.exterior.coords:
+ s += lon
+ c += 1
+ left = s / float(c) < -180 # and use it to decide where to shift the polygon
+ for (lon, lat) in polygon.exterior.coords:
+ ext_pts.append((lon + (-360, 360)[left], lat))
+ for interior in polygon.interiors:
+ pts = []
+ for (lon, lat) in interior.coords:
+ pts.append((lon + (-360, 360)[left], lat))
+ int_pts.append(pts)
+ polygons.append(Polygon(ext_pts, int_pts))print "shifted outside polygons to inside", (time.time() - start)
+ return polygons def _visible(self, lon, lat):
+ return True def _truncate(self, x, y):
+ return (x, y) def attrs(self):
+ a = super(Cylindrical, self).attrs()
+ a['lon0'] = self.lon0
+ a['flip'] = self.flip
+ return a def __str__(self):
+ return 'Proj(' + self.name + ', lon0=%s)' % self.lon0Equirectangular Projection, aka lonlat, aka plate carree
+ @staticmethod
+ def attributes():
+ return ['lon0', 'flip']
+
+ def ll(self, lon, lat):
+ if self.flip == 1:
+ return (-lon, -lat)
+ return (lon, lat)
+
+
+class Equirectangular(Cylindrical): def __init__(self, lon0=0.0, lat0=0.0, flip=0):
+ self.lat0 = lat0
+ self.phi0 = rad(lat0 * -1)
+ Cylindrical.__init__(self, lon0=lon0, flip=flip) def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ return (lon * math.cos(self.phi0) * 1000, lat * -1 * 1000)Cylindrical Equal Area Projection
+class CEA(Cylindrical): def __init__(self, lat0=0.0, lon0=0.0, lat1=0.0, flip=0):
+ self.lat0 = lat0
+ self.lat1 = lat1
+ self.phi0 = rad(lat0 * -1)
+ self.phi1 = rad(lat1 * -1)
+ self.lam0 = rad(lon0)
+ Cylindrical.__init__(self, lon0=lon0, flip=flip) def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ lam = rad(lon)
+ phi = rad(lat * -1)
+ x = (lam) * math.cos(self.phi1) * 1000
+ y = math.sin(phi) / math.cos(self.phi1) * 1000
+ return (x, y) def attrs(self):
+ p = super(CEA, self).attrs()
+ p['lat1'] = self.lat1
+ return p @staticmethod
+ def attributes():
+ return ['lon0', 'lat1', 'flip']
+
+ def __str__(self):
+ return 'Proj(' + self.name + ', lon0=%s, lat1=%s)' % (self.lon0, self.lat1)
+
+
+class GallPeters(CEA):
+ def __init__(self, lat0=0.0, lon0=0.0, flip=0):
+ CEA.__init__(self, lon0=lon0, lat0=0, lat1=45, flip=flip)
+
+
+class HoboDyer(CEA):
+ def __init__(self, lat0=0.0, lon0=0.0, flip=0):
+ CEA.__init__(self, lon0=lon0, lat0=lat0, lat1=37.5, flip=flip)
+
+
+class Behrmann(CEA):
+ def __init__(self, lat0=0.0, lon0=0.0, flip=0):
+ CEA.__init__(self, lat1=30, lat0=lat0, lon0=lon0, flip=flip)
+
+
+class Balthasart(CEA):
+ def __init__(self, lat0=0.0, lon0=0.0, flip=0):
+ CEA.__init__(self, lat1=50, lat0=lat0, lon0=lon0, flip=flip)
+
+
+class Mercator(Cylindrical):
+ def __init__(self, lon0=0.0, lat0=0.0, flip=0):
+ Cylindrical.__init__(self, lon0=lon0, flip=flip)
+ self.minLat = -85
+ self.maxLat = 85
+
+ def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ lam = rad(lon)
+ phi = rad(lat * -1)
+ x = lam * 1000
+ y = math.log((1 + math.sin(phi)) / math.cos(phi)) * 1000
+ return (x, y)
+
+
+class LonLat(Cylindrical):
+ def project(self, lon, lat):
+ return (lon, lat)
+
+kartograph - a svg mapping library
+Copyright (C) 2011 Gregor Aisch
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU Affero General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU Affero General Public License for more details.
+
+You should have received a copy of the GNU Affero General Public License
+along with this program. If not, see <http://www.gnu.org/licenses/>.
+from cylindrical import Cylindrical
+import math
+from math import radians as radclass PseudoCylindrical(Cylindrical): def __init__(self, lon0=0.0, flip=0):
+ Cylindrical.__init__(self, lon0=lon0, flip=flip)src: http://www.shadedrelief.com/NE_proj/
+class NaturalEarth(PseudoCylindrical): def __init__(self, lat0=0.0, lon0=0.0, flip=0):
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
+ from math import pi
+ s = self
+ s.A0 = 0.8707
+ s.A1 = -0.131979
+ s.A2 = -0.013791
+ s.A3 = 0.003971
+ s.A4 = -0.001529
+ s.B0 = 1.007226
+ s.B1 = 0.015085
+ s.B2 = -0.044475
+ s.B3 = 0.028874
+ s.B4 = -0.005916
+ s.C0 = s.B0
+ s.C1 = 3 * s.B1
+ s.C2 = 7 * s.B2
+ s.C3 = 9 * s.B3
+ s.C4 = 11 * s.B4
+ s.EPS = 1e-11
+ s.MAX_Y = 0.8707 * 0.52 * pi def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ lplam = rad(lon)
+ lpphi = rad(lat * -1)
+ phi2 = lpphi * lpphi
+ phi4 = phi2 * phi2
+ x = lplam * (self.A0 + phi2 * (self.A1 + phi2 * (self.A2 + phi4 * phi2 * (self.A3 + phi2 * self.A4)))) * 180 + 500
+ y = lpphi * (self.B0 + phi2 * (self.B1 + phi4 * (self.B2 + self.B3 * phi2 + self.B4 * phi4))) * 180 + 270
+ return (x, y)class Robinson(PseudoCylindrical): def __init__(self, lat0=0.0, lon0=0.0, flip=0):
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
+ self.X = [1, -5.67239e-12, -7.15511e-05, 3.11028e-06, 0.9986, -0.000482241, -2.4897e-05, -1.33094e-06, 0.9954, -0.000831031, -4.4861e-05, -9.86588e-07, 0.99, -0.00135363, -5.96598e-05, 3.67749e-06, 0.9822, -0.00167442, -4.4975e-06, -5.72394e-06, 0.973, -0.00214869, -9.03565e-05, 1.88767e-08, 0.96, -0.00305084, -9.00732e-05, 1.64869e-06, 0.9427, -0.00382792, -6.53428e-05, -2.61493e-06, 0.9216, -0.00467747, -0.000104566, 4.8122e-06, 0.8962, -0.00536222, -3.23834e-05, -5.43445e-06, 0.8679, -0.00609364, -0.0001139, 3.32521e-06, 0.835, -0.00698325, -6.40219e-05, 9.34582e-07, 0.7986, -0.00755337, -5.00038e-05, 9.35532e-07, 0.7597, -0.00798325, -3.59716e-05, -2.27604e-06, 0.7186, -0.00851366, -7.0112e-05, -8.63072e-06, 0.6732, -0.00986209, -0.000199572, 1.91978e-05, 0.6213, -0.010418, 8.83948e-05, 6.24031e-06, 0.5722, -0.00906601, 0.000181999, 6.24033e-06, 0.5322, 0., 0., 0.]
+ self.Y = [0, 0.0124, 3.72529e-10, 1.15484e-09, 0.062, 0.0124001, 1.76951e-08, -5.92321e-09, 0.124, 0.0123998, -7.09668e-08, 2.25753e-08, 0.186, 0.0124008, 2.66917e-07, -8.44523e-08, 0.248, 0.0123971, -9.99682e-07, 3.15569e-07, 0.31, 0.0124108, 3.73349e-06, -1.1779e-06, 0.372, 0.0123598, -1.3935e-05, 4.39588e-06, 0.434, 0.0125501, 5.20034e-05, -1.00051e-05, 0.4968, 0.0123198, -9.80735e-05, 9.22397e-06, 0.5571, 0.0120308, 4.02857e-05, -5.2901e-06, 0.6176, 0.0120369, -3.90662e-05, 7.36117e-07, 0.6769, 0.0117015, -2.80246e-05, -8.54283e-07, 0.7346, 0.0113572, -4.08389e-05, -5.18524e-07, 0.7903, 0.0109099, -4.86169e-05, -1.0718e-06, 0.8435, 0.0103433, -6.46934e-05, 5.36384e-09, 0.8936, 0.00969679, -6.46129e-05, -8.54894e-06, 0.9394, 0.00840949, -0.000192847, -4.21023e-06, 0.9761, 0.00616525, -0.000256001, -4.21021e-06, 1., 0., 0., 0]
+ self.NODES = 18
+ self.FXC = 0.8487
+ self.FYC = 1.3523
+ self.C1 = 11.45915590261646417544
+ self.RC1 = 0.08726646259971647884
+ self.ONEEPS = 1.000001
+ self.EPS = 1e-8 def _poly(self, arr, off, z):
+ return arr[off] + z * (arr[off + 1] + z * (arr[off + 2] + z * (arr[off + 3]))) def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ lplam = rad(lon)
+ lpphi = rad(lat * -1)
+
+ phi = abs(lpphi)
+ i = int(phi * self.C1)
+ if i >= self.NODES:
+ i = self.NODES - 1
+ phi = math.degrees(phi - self.RC1 * i)
+ i *= 4
+ x = 1000 * self._poly(self.X, i, phi) * self.FXC * lplam
+ y = 1000 * self._poly(self.Y, i, phi) * self.FYC
+ if lpphi < 0.0:
+ y = -y
+
+ return (x, y)class EckertIV(PseudoCylindrical): def __init__(self, lon0=0.0, lat0=0, flip=0):
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
+
+ self.C_x = .42223820031577120149
+ self.C_y = 1.32650042817700232218
+ self.RC_y = .75386330736002178205
+ self.C_p = 3.57079632679489661922
+ self.RC_p = .28004957675577868795
+ self.EPS = 1e-7
+ self.NITER = 6 def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ lplam = rad(lon)
+ lpphi = rad(lat * -1)
+
+ p = self.C_p * math.sin(lpphi)
+ V = lpphi * lpphi
+ lpphi *= 0.895168 + V * (0.0218849 + V * 0.00826809)
+
+ i = self.NITER
+ while i > 0:
+ c = math.cos(lpphi)
+ s = math.sin(lpphi)
+ V = (lpphi + s * (c + 2.) - p) / (1. + c * (c + 2.) - s * s)
+ lpphi -= V
+ if abs(V) < self.EPS:
+ break
+ i -= 1
+
+ if i == 0:
+ x = self.C_x * lplam
+ y = (self.C_y, - self.C_y)[lpphi < 0]
+ else:
+ x = self.C_x * lplam * (1. + math.cos(lpphi))
+ y = self.C_y * math.sin(lpphi)
+ return (x, y)class Sinusoidal(PseudoCylindrical): def __init__(self, lon0=0.0, lat0=0.0, flip=0):
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip) def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ lam = rad(lon)
+ phi = rad(lat * -1)
+ x = 1032 * lam * math.cos(phi)
+ y = 1032 * phi
+ return (x, y)class Mollweide(PseudoCylindrical): def __init__(self, p=1.5707963267948966, lon0=0.0, lat0=0.0, cx=None, cy=None, cp=None, flip=0):
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
+ self.MAX_ITER = 10
+ self.TOLERANCE = 1e-7
+
+ if p != None:
+ p2 = p + p
+ sp = math.sin(p)
+ r = math.sqrt(math.pi * 2.0 * sp / (p2 + math.sin(p2)))
+ self.cx = 2. * r / math.pi
+ self.cy = r / sp
+ self.cp = p2 + math.sin(p2)
+ elif cx != None and cy != None and cp != None:
+ self.cx = cx
+ self.cy = cy
+ self.cp = cp
+ else:
+ assert False, 'either p or cx,cy,cp must be defined' def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ lam = rad(lon)
+ phi = rad(lat)
+
+ k = self.cp * math.sin(phi)
+ i = self.MAX_ITER
+
+ while i != 0:
+ v = (phi + math.sin(phi) - k) / (1. + math.cos(phi))
+ phi -= v
+ if abs(v) < self.TOLERANCE:
+ break
+ i -= 1
+
+ if i == 0:
+ phi = (self.HALFPI, -self.HALFPI)[phi < 0]
+ else:
+ phi *= 0.5
+
+ x = 1000 * self.cx * lam * math.cos(phi)
+ y = 1000 * self.cy * math.sin(phi)
+ return (x, y * -1)class GoodeHomolosine(PseudoCylindrical): def __init__(self, lon0=0, flip=0):
+ self.lat1 = 41.737
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
+ self.p1 = Mollweide()
+ self.p0 = Sinusoidal() def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)lon = me.clon(lon)
+ if abs(lat) > self.lat1:
+ return self.p1.project(lon, lat)
+ else:
+ return self.p0.project(lon, lat)class WagnerIV(Mollweide): def __init__(self, lon0=0, lat0=0, flip=0):p=math.pi/3
+ Mollweide.__init__(self, p=1.0471975511965976, flip=flip)class WagnerV(Mollweide): def __init__(self, lat0=0, lon0=0, flip=0):
+ Mollweide.__init__(self, cx=0.90977, cy=1.65014, cp=3.00896, flip=flip)class Loximuthal(PseudoCylindrical):
+
+ minLat = -89
+ maxLat = 89 def __init__(self, lon0=0.0, lat0=0.0, flip=0):
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
+ if flip == 1:
+ lat0 = -lat0
+ self.lat0 = lat0
+ self.phi0 = rad(lat0) def project(self, lon, lat):
+ lon, lat = self.ll(lon, lat)
+ lam = rad(lon)
+ phi = rad(lat)
+ if phi == self.phi0:
+ x = lam * math.cos(self.phi0)
+ else:
+ try:
+ x = lam * (phi - self.phi0) / (math.log(math.tan(self.QUARTERPI + phi * 0.5)) - math.log(math.tan(self.QUARTERPI + self.phi0 * 0.5)))
+ except:
+ return None
+ x *= 1000
+ y = 1000 * (phi - self.phi0)
+ return (x, y * -1) def attrs(self):
+ p = super(Loximuthal, self).attrs()
+ p['lat0'] = self.lat0
+ return pCanters, F. (2002) Small-scale Map projection Design. p. 218-219. +Modified Sinusoidal, equal-area.
+implementation borrowed from +http://cartography.oregonstate.edu/temp/AdaptiveProjection/src/projections/Canters1.js
+ @staticmethod
+ def attributes():
+ return ['lon0', 'lat0', 'flip']
+
+
+class CantersModifiedSinusoidalI(PseudoCylindrical): def __init__(self, lon0=0.0, flip=0):
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
+ self.C1 = 1.1966
+ self.C3 = -0.1290
+ self.C3x3 = 3 * self.C3
+ self.C5 = -0.0076
+ self.C5x5 = 5 * self.C5 def project(self, lon, lat):
+ me = self
+ lon, lat = me.ll(lon, lat)
+
+ lon = rad(lon)
+ lat = rad(lat)
+
+ y2 = lat * lat
+ y4 = y2 * y2
+ x = 1000 * lon * math.cos(lat) / (me.C1 + me.C3x3 * y2 + me.C5x5 * y4)
+ y = 1000 * lat * (me.C1 + me.C3 * y2 + me.C5 * y4)
+ return (x, y * -1)class Hatano(PseudoCylindrical): def __init__(me, lon0=0, flip=0):
+ PseudoCylindrical.__init__(me, lon0=lon0, flip=flip)
+ me.NITER = 20
+ me.EPS = 1e-7
+ me.ONETOL = 1.000001
+ me.CN = 2.67595
+ me.CS = 2.43763
+ me.RCN = 0.37369906014686373063
+ me.RCS = 0.41023453108141924738
+ me.FYCN = 1.75859
+ me.FYCS = 1.93052
+ me.RYCN = 0.56863737426006061674
+ me.RYCS = 0.51799515156538134803
+ me.FXC = 0.85
+ me.RXC = 1.17647058823529411764 def project(me, lon, lat):
+ [lon, lat] = me.ll(lon, lat)
+ lam = rad(lon)
+ phi = rad(lat)
+ c = math.sin(phi) * (me.CN, me.CS)[phi < 0.0]
+ for i in range(me.NITER, 0, -1):
+ th1 = (phi + math.sin(phi) - c) / (1.0 + math.cos(phi))
+ phi -= th1
+ if abs(th1) < me.EPS:
+ break
+ phi *= 0.5
+ x = 1000 * me.FXC * lam * math.cos(phi)
+ y = 1000 * math.sin(phi) * (me.FYCN, me.FYCS)[phi < 0.0]
+ return (x, y * -1)Aitoff projection
+implementation taken from +Snyder, Map projections - A working manual
+class Aitoff(PseudoCylindrical): def __init__(self, lon0=0, flip=0):
+ PseudoCylindrical.__init__(self, lon0=lon0, flip=flip)
+ self.winkel = False
+ self.COSPHI1 = 0.636619772367581343 def project(me, lon, lat):
+ [lon, lat] = me.ll(lon, lat)
+ lam = rad(lon)
+ phi = rad(lat)
+ c = 0.5 * lam
+ d = math.acos(math.cos(phi) * math.cos(c))
+ if d != 0:
+ y = 1.0 / math.sin(d)
+ x = 2.0 * d * math.cos(phi) * math.sin(c) * y
+ y *= d * math.sin(phi)
+ else:
+ x = y = 0
+ if me.winkel:
+ x = (x + lam * me.COSPHI1) * 0.5
+ y = (y + phi) * 0.5
+ return (x * 1000, y * -1000)class Winkel3(Aitoff): def __init__(self, lon0=0, flip=0):
+ Aitoff.__init__(self, lon0=lon0, flip=flip)
+ self.winkel = Trueclass Nicolosi(PseudoCylindrical): def __init__(me, lon0=0, flip=0):
+ me.EPS = 1e-10
+ PseudoCylindrical.__init__(me, lon0=lon0, flip=flip)
+ me.r = me.HALFPI * 100
+ sea = []
+ r = me.r
+ for phi in range(0, 361):
+ sea.append((math.cos(rad(phi)) * r, math.sin(rad(phi)) * r))
+ me.sea = sea def _clon(me, lon):
+ lon -= me.lon0
+ if lon < -180:
+ lon += 360
+ elif lon > 180:
+ lon -= 360
+ return lon def _visible(me, lon, lat):lon = me._clon(lon)
+ return lon > -90 and lon < 90 def _truncate(me, x, y):
+ theta = math.atan2(y, x)
+ x1 = me.r * math.cos(theta)
+ y1 = me.r * math.sin(theta)
+ return (x1, y1) def world_bounds(self, bbox, llbbox=(-180, -90, 180, 90)):
+ if llbbox == (-180, -90, 180, 90):
+ d = self.r * 2
+ bbox.update((-d, -d))
+ bbox.update((d, d))
+ else:
+ bbox = super(PseudoCylindrical, self).world_bounds(bbox, llbbox)
+ return bbox def sea_shape(self, llbbox=(-180, -90, 180, 90)):
+ out = []
+ if llbbox == (-180, -90, 180, 90) or llbbox == [-180, -90, 180, 90]:
+ for phi in range(0, 360):
+ x = math.cos(math.radians(phi)) * self.r
+ y = math.sin(math.radians(phi)) * self.r
+ out.append((x, y))
+ out = [out]
+ else:
+ out = super(PseudoCylindrical, self).sea_shape(llbbox)
+ return out def project(me, lon, lat):
+ [lon, lat] = me.ll(lon, lat)
+ lam = rad(lon)
+ phi = rad(lat)
+
+ if abs(lam) < me.EPS:
+ x = 0
+ y = phi
+ elif abs(phi) < me.EPS:
+ x = lam
+ y = 0
+ elif abs(abs(lam) - me.HALFPI) < me.EPS:
+ x = lam * math.cos(phi)
+ y = me.HALFPI * math.sin(phi)
+ elif abs(abs(phi) - me.HALFPI) < me.EPS:
+ x = 0
+ y = phi
+ else:
+ tb = me.HALFPI / lam - lam / me.HALFPI
+ c = phi / me.HALFPI
+ sp = math.sin(phi)
+ d = (1 - c * c) / (sp - c)
+ r2 = tb / d
+ r2 *= r2
+ m = (tb * sp / d - 0.5 * tb) / (1.0 + r2)
+ n = (sp / r2 + 0.5 * d) / (1.0 + 1.0 / r2)
+ x = math.cos(phi)
+ x = math.sqrt(m * m + x * x / (1.0 + r2))
+ x = me.HALFPI * (m + (x, -x)[lam < 0])
+ f = n * n - (sp * sp / r2 + d * sp - 1.0) / (1.0 + 1.0 / r2)
+ if f < 0:
+ y = phi
+ else:
+ y = math.sqrt(f)
+ y = me.HALFPI * (n + (-y, y)[phi < 0])
+ return (x * 100, y * -100) def plot(self, polygon, truncate=True):
+ polygons = self._shift_polygon(polygon)
+ plotted = []
+ for polygon in polygons:
+ points = []
+ ignore = True
+ for (lon, lat) in polygon:
+ vis = self._visible(lon, lat)
+ if vis:
+ ignore = False
+ x, y = self.project(lon, lat)
+ if not vis and truncate:
+ points.append(self._truncate(x, y))
+ else:
+ points.append((x, y))
+ if ignore:
+ continue
+ plotted.append(points)
+ return plotted
+
+class MapRenderer(object): def __init__(self, map):
+ self.map = map def render(self):
+ pass def write(self, filename):
+ raise 'Not implemented yet' def preview(self):
+ raise 'Not implemented yet'
+
+
+from svg import SvgRenderer
+from kml import KmlRenderer
+
+__all__ = ['MapRenderer', 'SvgRenderer', 'KmlRenderer']
+
+from kartograph.renderer import MapRenderer
+from kartograph.errors import KartographError
+from pykml.factory import KML_ElementMaker as KMLclass KmlRenderer(MapRenderer): def render(self):
+ self._init_kml_canvas()
+ self._store_layers_kml() def write(self, filename):
+ outfile = open(filename, 'w')
+ from lxml import etree
+ outfile.write(etree.tostring(self.kml, pretty_print=True)) def preview(self):
+ self.write('tmp.kml')
+ from subprocess import call
+ call(["open", "tmp.kml"]) def _init_kml_canvas(self):
+ self.kml = KML.kml(
+ KML.Document(
+ KML.name('kartograph map')
+ )
+ )store features in kml (projected to WGS84 latlon)
+ def _store_layers_kml(self): from pykml.factory import KML_ElementMaker as KML
+
+ for layer in self.map.layers:
+ if len(layer.features) == 0:
+ continue # ignore empty layers
+ g = KML.Folder(
+ KML.name(id)
+ )
+ for feat in layer.features:
+ g.append(self._render_feature(feat, layer.options['attributes']))
+ self.kml.Document.append(g) def _render_feature(self, feature, attributes=[]):
+ path = self._render_geometry(feature.geometry)
+
+ pm = KML.Placemark(
+ KML.name(unicode(feature.props[attributes[0]['src']])),
+ path
+ )
+
+ xt = KML.ExtendedData()
+
+ for cfg in attributes:
+ if 'src' in cfg:
+ if cfg['src'] not in feature.props:
+ continueraise KartographError(('attribute not found "%s"'%cfg['src']))
+ val = feature.props[cfg['src']]
+ import unicodedata
+ if isinstance(val, str):
+ val = unicode(val, errors='ignore')
+ val = unicodedata.normalize('NFKD', val).encode('ascii', 'ignore')
+ xt.append(KML.Data(
+ KML.value(unicode(val)),
+ name=cfg['tgt']
+ ))
+ elif 'where' in cfg:
+ src = cfg['where']
+ tgt = cfg['set']
+ if len(cfg['equals']) != len(cfg['to']):
+ raise KartographError('attributes: "equals" and "to" arrays must be of same length')
+ for i in range(len(cfg['equals'])):
+ if feature.props[src] == cfg['equals'][i]:
+ xt.append(KML.Data(
+ KML.value(cfg['to'][i]),
+ name=tgt
+ ))
+ pm.append(xt)
+
+ return pm def _render_geometry(self, geometry):
+ from shapely.geometry import Polygon, MultiPolygon
+ if isinstance(geometry, (Polygon, MultiPolygon)):
+ return self._render_polygon(geometry)
+ else:
+ raise KartographError('kml-renderer is not fully implemented yet')renders a Polygon or MultiPolygon as KML node
+ def _render_polygon(self, geometry): geoms = hasattr(geometry, 'geoms') and geometry.geoms or [geometry]
+
+ kml_polys = []
+
+ for geom in geoms:
+ poly = KML.Polygon(
+ KML.tesselate("1")
+ )
+ outer = KML.outerBoundaryIs()
+ coords = ''
+ for pt in geom.exterior.coords:
+ coords += ','.join(map(str, pt)) + ' '
+ outer.append(KML.LinearRing(
+ KML.coordinates(coords)
+ ))
+ poly.append(outer)
+ inner = KML.innerBoundaryIs()
+ for hole in geom.interiors:
+ coords = ''
+ for pt in hole.coords:
+ coords += ','.join(map(str, pt)) + ' '
+ inner.append(KML.LinearRing(
+ KML.coordinates(coords)
+ ))
+ if len(inner) > 0:
+ poly.append(inner)
+ kml_polys.append(poly)
+
+ if len(kml_polys) == 1:
+ return kml_polys[0]
+ multigeometry = KML.MultiGeometry()
+ for p in kml_polys:
+ multigeometry.append(p)
+ return multigeometry
+
+from kartograph.renderer import MapRenderer
+from kartograph.errors import KartographErrorThis script contains everything that is needed by Kartograph to finally +render the processed maps into SVG files.
+The SVG renderer is based on xml.dom.minidom.
+from xml.dom import minidom
+from xml.dom.minidom import parse
+import reclass SvgRenderer(MapRenderer):The render() method prepares a new empty SVG document and +stores all the layer features into SVG groups.
+ def render(self): self._init_svg_doc()
+ self._store_layers_to_svg() def _init_svg_doc(self):Load width and height of the map view +We add two pixels to the height to ensure that +the map fits.
+ w = self.map.view.width
+ h = self.map.view.height + 2SvgDocument is a handy wrapper around xml.dom.minidom. It is defined below.
+ svg = SvgDocument(
+ width='%dpx' % w,
+ height='%dpx' % h,
+ viewBox='0 0 %d %d' % (w, h),
+ enable_background='new 0 0 %d %d' % (w, h),
+ style='stroke-linejoin: round; stroke:#000; fill:#f6f3f0;')
+
+ defs = svg.node('defs', svg.root)
+ style = svg.node('style', defs, type='text/css')
+ css = 'path { fill-rule: evenodd; }\n#context path { fill: #eee; stroke: #bbb; } '
+ svg.cdata(css, style)
+ metadata = svg.node('metadata', svg.root)
+ views = svg.node('views', metadata)
+ view = svg.node('view', views,
+ padding=str(self.map.options['bounds']['padding']), w=w, h=h)
+
+ svg.node('proj', view, **self.map.proj.attrs())
+ svg.node('bbox', view,
+ x=round(self.map.src_bbox.left, 2),
+ y=round(self.map.src_bbox.top, 2),
+ w=round(self.map.src_bbox.width, 2),
+ h=round(self.map.src_bbox.height, 2))
+
+ ll = [-180, -90, 180, 90]
+ if self.map.options['bounds']['mode'] == "bbox":
+ ll = self.map.options['bounds']['data']
+ svg.node('llbbox', view,
+ lon0=ll[0], lon1=ll[2],
+ lat0=ll[1], lat1=ll[3])
+ self.svg = svgstore features in svg
+ def _store_layers_to_svg(self): svg = self.svg
+ for layer in self.map.layers:
+ if len(layer.features) == 0:
+ print "ignoring empty layer", layer.id
+ continue # ignore empty layers
+ g = svg.node('g', svg.root, id=layer.id)
+ for feat in layer.features:
+ node = self._render_feature(feat, layer.options['attributes'])
+ if node is not None:
+ g.appendChild(node)
+ else:
+ print "feature.to_svg is None", feat
+ if 'styles' in layer.options:
+ for prop in layer.options['styles']:
+ g.setAttribute(prop, str(layer.options['styles'][prop])) def _render_feature(self, feature, attributes=[]):
+ node = self._render_geometry(feature.geometry)
+ if node is None:
+ return None
+
+ for cfg in attributes:
+ if 'src' in cfg:
+ tgt = re.sub('(\W|_)+', '-', cfg['tgt'].lower())
+ if cfg['src'] not in feature.props:
+ continueraise KartographError(('attribute not found "%s"'%cfg['src']))
+ val = feature.props[cfg['src']]
+ if isinstance(val, (int, float)):
+ val = str(val)
+ node.setAttribute('data-' + tgt, val)
+ if tgt == "id":
+ node.setAttribute('id', val)
+
+ elif 'where' in cfg:can be used to replace attributes...
+ src = cfg['where']
+ tgt = cfg['set']
+ if len(cfg['equals']) != len(cfg['to']):
+ raise KartographError('attributes: "equals" and "to" arrays must be of same length')
+ for i in range(len(cfg['equals'])):
+ if feature.props[src] == cfg['equals'][i]:
+ node.setAttribute('data-' + tgt, cfg['to'][i])
+
+ if '__color__' in feature.props:
+ node.setAttribute('fill', self.props['__color__'])
+ return node def _render_geometry(self, geometry):
+ from shapely.geometry import Polygon, MultiPolygon
+ if isinstance(geometry, (Polygon, MultiPolygon)):
+ return self._render_polygon(geometry)constructs a svg representation of a polygon
+ def _render_polygon(self, geometry): _round = self.map.options['export']['round']
+ path_str = ""
+ if _round is False:
+ fmt = '%f,%f'
+ else:
+ fmt = '%.' + str(_round) + 'f'
+ fmt = fmt + ',' + fmt
+
+ geoms = hasattr(geometry, 'geoms') and geometry.geoms or [geometry]
+ for polygon in geoms:
+ if polygon is None:
+ continue
+ for ring in [polygon.exterior] + list(polygon.interiors):
+ cont_str = ""
+ kept = []
+ for pt in ring.coords:
+ kept.append(pt)
+ if len(kept) <= 3:
+ continue
+ for pt in kept:
+ if cont_str == "":
+ cont_str = "M"
+ else:
+ cont_str += "L"
+ cont_str += fmt % pt
+ cont_str += "Z "
+ path_str += cont_str
+ if path_str == "":
+ return None
+ path = self.svg.node('path', d=path_str)
+ return path def write(self, filename):
+ self.svg.write(filename) def preview(self):
+ self.svg.preview()SVGDocument is a handy wrapper around xml.dom.minidom which allows us +to quickly build XML structures. It is largely inspired by the SVG class +of the svgfig project, which was +used by one of the earlier versions of Kartograph.
+class SvgDocument(object):Of course, we need to create and XML document with all this +boring SVG header stuff added to it.
+ def __init__(self, **kwargs):
+ imp = minidom.getDOMImplementation('')
+ dt = imp.createDocumentType('svg',
+ '-//W3C//DTD SVG 1.1//EN',
+ 'http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd')
+ self.doc = imp.createDocument('http://www.w3.org/2000/svg', 'svg', dt)
+ self.root = svg = self.doc.getElementsByTagName('svg')[0]
+ svg.setAttribute('xmlns', 'http://www.w3.org/2000/svg')
+ svg.setAttribute('version', '1.1')
+ svg.setAttribute('xmlns:xlink', 'http://www.w3.org/1999/xlink')
+ _add_attrs(self.root, kwargs)This is the magic of SvgDocument. Instead of having to do appendChild() +and addAttribute() for every node we create, we just call svgdoc.node() +which is smart enough to append itself to the parent if we specify one, +and also sets all attributes we pass as keyword arguments.
+ def node(self, name, parent=None, **kwargs):
+ el = self.doc.createElement(name)
+ _add_attrs(el, kwargs)
+ if parent is not None:
+ parent.appendChild(el)
+ return elSometimes we also need a <[CDATA]> block, for instance if we embed +CSS code in the SVG document.
+ def cdata(self, data, parent=None):
+ cd = minidom.CDATASection()
+ cd.data = data
+ if parent is not None:
+ parent.appendChild(cd)
+ return cdHere we finally write the SVG file, and we're brave enough +to try to write it in Unicode.
+ def write(self, outfile):
+ if isinstance(outfile, str):
+ outfile = open(outfile, 'w')
+ raw = self.doc.toxml()
+ try:
+ raw = raw.encode('utf-8')
+ except:
+ print 'warning: could not encode to unicode'
+
+ outfile.write(raw)
+ outfile.close()Don't blame me if you don't have a command-line shortcut to +simply the best free browser of the world.
+ def preview(self):
+ self.write('tmp.svg')
+ from subprocess import call
+ call(["firefox", "tmp.svg"]) def tostring(self):
+ return self.doc.toxml()This is an artifact of an older version of Kartograph, but +maybe we'll need it later. It will load an SVG document from +a file.
+ @staticmethod
+ def load(filename):
+ svg = SvgDocument()
+ dom = parse(filename)
+ svg.doc = dom
+ svg.root = dom.getElementsByTagName('svg')[0]
+ return svg
+
+
+def _add_attrs(node, attrs):
+ for key in attrs:
+ node.setAttribute(key, str(attrs[key]))
+
+__all__ = ['create_point_store', 'unify_polygon', 'unify_polygons', 'simplify_lines']
+
+from unify import *
+from distance import simplify_distance
+from douglas_peucker import simplify_douglas_peucker
+from visvalingam import simplify_visvalingam_whyatt
+
+
+simplification_methods = dict()
+simplification_methods['distance'] = simplify_distance
+simplification_methods['douglas-peucker'] = simplify_douglas_peucker
+simplification_methods['visvalingam-whyatt'] = simplify_visvalingam_whyattsimplifies a set of lines
+def simplify_lines(lines, method, params): simplify = simplification_methods[method]
+ out = []
+ for line in lines:remove duplicate points from line
+ unique = [line[0]]
+ for i in range(1, len(line)):
+ if line[i] != line[i - 1]:
+ unique.append(line[i])
+ simplify(unique, params)
+ out.append(unique)
+ return out
+
+simplifies a line segment using a very simple algorithm that checks the distance +to the last non-deleted point. the algorithm operates on line segments.
+in order to preserve topology of the original polygons the algorithm +- never removes the first or the last point of a line segment +- flags all points as simplified after processing (so it won't be processed twice)
+def simplify_distance(points, dist): dist_sq = dist * dist
+ n = len(points)
+
+ kept = []
+ deleted = 0
+ if n < 4:
+ return points
+
+ for i in range(0, n):
+ pt = points[i]
+ if i == 0 or i == n - 1:never remove first or last point of line
+ pt.simplified = True
+ lpt = pt
+ kept.append(pt)
+ else:
+ d = (pt.x - lpt.x) * (pt.x - lpt.x) + (pt.y - lpt.y) * (pt.y - lpt.y) # compute distance to last point
+ if pt.simplified or d > dist_sq: # if point already handled or distance exceeds threshold..
+ kept.append(pt) # ..keep the point
+ lpt = pt
+ else: # otherwise remove it
+ deleted += 1
+ pt.deleted = True
+ pt.simplified = True
+
+ if len(kept) < 4:
+ for pt in points:
+ pt.deleted = False
+ return points
+
+ return keptprint 'kept %d deleted %d' % (kept, deleted)
+simplifies a line segment using the Douglas-Peucker algorithm.
+taken from +http://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm#Pseudocode
+in order to preserve topology of the original polygons the algorithm +- never removes the first or the last point of a line segment +- flags all points as simplified after processing (so it won't be processed twice)
+def simplify_douglas_peucker(points, epsilon): n = len(points)
+ kept = []
+ if n < 4:
+ return points # skip short lines
+
+ if not points[0].simplified:
+ _douglas_peucker(points, 0, n - 1, epsilon)
+
+ return keptprint 'kept %d deleted %d' % (kept, deleted)
+inner part of Douglas-Peucker algorithm, called recursively
+def _douglas_peucker(points, start, end, epsilon): dmax = 0
+ index = 0Find the point with the maximum distance
+ for i in range(start + 1, end):
+ x1, y1 = points[start]
+ x2, y2 = points[end]
+ if x1 == x2 and y1 == y2:
+ return
+ x3, y3 = points[i]
+ d = _min_distance(x1, y1, x2, y2, x3, y3)
+ if d > dmax:
+ index = i
+ dmax = dIf max distance is greater than epsilon, recursively simplify
+ if dmax >= epsilon and start < index < end:recursivly call
+ _douglas_peucker(points, start, index, epsilon)
+ _douglas_peucker(points, index, end, epsilon)
+ else:remove any point but the first and last
+ for i in range(start, end + 1):
+ points[i].deleted = i == start or i == end
+ points[i].simplified = Truethe perpendicular distance from a point (x3,y3) to the line from (x1,y1) to (x2,y2) +taken from http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/
+def _min_distance(x1, y1, x2, y2, x3, y3): d = _dist(x1, y1, x2, y2)
+ u = (x3 - x1) * (x2 - x1) + (y3 - y1) * (y2 - y1) / (d * d)
+ x = x1 + u * (x2 - x1)
+ y = y1 + u * (y2 - y1)
+ return _dist(x, y, x3, y3)eucledian distance between two points
+def _dist(x1, y1, x2, y2): import math
+ dx = x2 - x1
+ dy = y2 - y1
+ return math.sqrt(dx * dx + dy * dy)
+
+Point class used for polygon simplification
+class MPoint: def __init__(self, x, y):
+ self.x = x
+ self.y = y
+ self.simplified = False
+ self.deleted = False
+ self.keep = False
+ self.features = set() def isDeletable(self):
+ if self.keep or self.simplified or self.three:
+ return False
+ return True def __repr__(self):
+ return 'Pt(%.2f,%.2f)' % (self.x, self.y) def __len__(self):
+ return 2 def __getitem__(self, key):
+ if key == 0:
+ return self.x
+ if key == 1:
+ return self.y
+ raise IndexError() def __contains__(self, key):
+ if key == "deleted":
+ return True
+ return False
+
+the whole point of the unification step is to convert all points into unique MPoint instances
+from mpoint import MPointcreates a new point_store
+def create_point_store(): point_store = {'kept': 0, 'removed': 0}
+ return point_storedef unify_rings(rings, point_store, precision=None, feature=None):
+ out = []
+ for ring in rings:
+ out.append(unify_ring(ring, point_store, precision=precision, feature=feature))
+ return outReplaces duplicate points with MPoint instances
+def unify_ring(ring, point_store, precision=None, feature=None): out_ring = []
+ lptid = ''
+ for pt in ring:
+ if 'deleted' not in pt:
+ pt = MPoint(pt[0], pt[1]) # eventually convert to MPointgenerate hash for point
+ if precision is not None:
+ fmt = '%' + precision + 'f-%' + precision + 'f'
+ else:
+ fmt = '%f-%f'
+ pid = fmt % (pt.x, pt.y)
+ if pid == lptid:
+ continue # skip double points
+ lptid = pid
+ if pid in point_store:load existing point from point store
+ point = point_store[pid]
+ point_store['removed'] += 1
+ else:
+ point = pt
+ point_store['kept'] += 1
+ point_store[pid] = pt
+
+ point.features.add(feature)
+ out_ring.append(point)
+ return out_ring
+
+Visvalingam-Whyatt simplification
+def simplify_visvalingam_whyatt(points, tolerance):mplementation borrowed from @migurski: +ttps://github.com/migurski/Bloch/blob/master/Bloch/init.py#L133
+ if len(points) < 3:
+ return
+ if points[1].simplified:
+ return
+
+ min_area = tolerance ** 2
+
+ pts = range(len(points)) # pts stores an index of all non-deleted points
+
+ while len(pts) > 4:
+ preserved, popped = set(), []
+ areas = []
+
+ for i in range(1, len(pts) - 1):
+ x1, y1 = points[pts[i - 1]]
+ x2, y2 = points[pts[i]]
+ x3, y3 = points[pts[i + 1]]compute and store triangle area
+ areas.append((_tri_area(x1, y1, x2, y2, x3, y3), i))
+
+ areas = sorted(areas)
+
+ if not areas or areas[0][0] > min_area:there's nothing to be done
+ for pt in points:
+ pt.simplified = True
+ breakReduce any segments that makes a triangle whose area is below +the minimum threshold, starting with the smallest and working up. +Mark segments to be preserved until the next iteration.
+ for (area, i) in areas:
+
+ if area > min_area:there won't be any more points to remove.
+ break
+
+ if i - 1 in preserved or i + 1 in preserved:the current segment is too close to a previously-preserved one. +print "-pre", preserved
+ continue
+
+ points[pts[i]].deleted = True
+ popped.append(i)make sure that the adjacent points
+ preserved.add(i - 1)
+ preserved.add(i + 1)
+
+ if len(popped) == 0:no points removed, so break out of loop
+ break
+
+ popped = sorted(popped, reverse=True)
+ for i in popped:remove point from index list
+ pts = pts[:i] + pts[i + 1:]
+
+ for pt in points:
+ pt.simplified = Truecomputes the area of a triangle given by three points +implementation taken from: +http://www.btinternet.com/~se16/hgb/triangle.htm
+def _tri_area(x1, y1, x2, y2, x3, y3): return abs((x2*y1-x1*y2)+(x3*y2-x2*y3)+(x1*y3-x3*y1))/2.0
+
+