jqvmap.py 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573
  1. import sys
  2. import json
  3. import csv
  4. import shapely.wkb
  5. import shapely.geometry
  6. import shapely.ops
  7. import os
  8. import copy
  9. from osgeo import ogr
  10. from osgeo import osr
  11. from booleano.parser import Grammar, EvaluableParseManager, SymbolTable, Bind
  12. from booleano.operations import Variable
  13. class JQVMap:
  14. def __init__(self, name, language):
  15. self.paths = {}
  16. self.name = name
  17. self.language = language
  18. self.width = 0
  19. self.height = 0
  20. self.bbox = []
  21. def addPath(self, path, code, name):
  22. self.paths[code] = {"path": path, "name": name}
  23. def getJSCode(self):
  24. map = {"paths": self.paths, "width": self.width, "height": self.height, "insets": self.insets, "projection": self.projection}
  25. header = "/** JQVMap " + self.projection['type'] + " map for " + self.name + " */"
  26. js = "jQuery.fn.vectorMap('addMap', '" + self.name + "'," + json.dumps(map) + ");"
  27. return header + "\n" + js
  28. class Converter:
  29. def __init__(self, config):
  30. args = {
  31. 'buffer_distance': -0.4,
  32. 'simplify_tolerance': 0.2,
  33. 'longitude0': 0,
  34. 'projection': 'mill',
  35. 'name': 'world',
  36. 'width': 900,
  37. 'left': 0,
  38. 'top': 0,
  39. 'language': 'en',
  40. 'precision': 2,
  41. 'insets': []
  42. }
  43. args.update(config)
  44. self.config = args
  45. self.map = JQVMap(args['name'], args.get('language'))
  46. if args.get('sources'):
  47. self.sources = args['sources']
  48. else:
  49. self.sources = [{
  50. 'input_file': args.get('input_file'),
  51. 'where': args.get('where'),
  52. 'name_field': args.get('name_field'),
  53. 'code_field': args.get('code_field'),
  54. 'input_file_encoding': args.get('input_file_encoding')
  55. }]
  56. default_source = {
  57. 'where': '',
  58. 'name_field': 0,
  59. 'code_field': 1,
  60. 'input_file_encoding': 'iso-8859-1'
  61. }
  62. for index in range(len(self.sources)):
  63. for key in default_source:
  64. if self.sources[index].get(key) is None:
  65. self.sources[index][key] = default_source[key]
  66. self.features = {}
  67. self.width = args.get('width')
  68. self.left = args.get('left')
  69. self.top = args.get('top')
  70. self.minimal_area = args.get('minimal_area')
  71. self.longitude0 = float(args.get('longitude0'))
  72. self.projection = args.get('projection')
  73. self.precision = args.get('precision')
  74. self.buffer_distance = args.get('buffer_distance')
  75. self.simplify_tolerance = args.get('simplify_tolerance')
  76. self.for_each = args.get('for_each')
  77. self.emulate_longitude0 = args.get('emulate_longitude0')
  78. if args.get('emulate_longitude0') is None and (self.projection == 'merc' or self.projection =='mill') and self.longitude0 != 0:
  79. self.emulate_longitude0 = True
  80. if args.get('viewport'):
  81. self.viewport = map(lambda s: float(s), args.get('viewport').split(' '))
  82. else:
  83. self.viewport = False
  84. # spatial reference to convert to
  85. self.spatialRef = osr.SpatialReference()
  86. projString = '+proj='+str(self.projection)+' +a=6381372 +b=6381372 +lat_0=0'
  87. if not self.emulate_longitude0:
  88. projString += ' +lon_0='+str(self.longitude0)
  89. self.spatialRef.ImportFromProj4(projString)
  90. # handle map insets
  91. if args.get('insets'):
  92. self.insets = args.get('insets')
  93. else:
  94. self.insets = []
  95. def convert(self, data_source, output_file):
  96. codes = map(lambda g: g.properties[self.config['code_field']], data_source.geometries)
  97. main_codes = copy.copy(codes)
  98. self.map.insets = []
  99. envelope = []
  100. for inset in self.insets:
  101. insetBbox = self.renderMapInset(data_source, inset['codes'], inset['left'], inset['top'], inset['width'])
  102. insetHeight = (insetBbox[3] - insetBbox[1]) * (inset['width'] / (insetBbox[2] - insetBbox[0]))
  103. self.map.insets.append({
  104. "bbox": [{"x": insetBbox[0], "y": -insetBbox[3]}, {"x": insetBbox[2], "y": -insetBbox[1]}],
  105. "left": inset['left'],
  106. "top": inset['top'],
  107. "width": inset['width'],
  108. "height": insetHeight
  109. })
  110. envelope.append(
  111. shapely.geometry.box(
  112. inset['left'], inset['top'], inset['left'] + inset['width'], inset['top'] + insetHeight
  113. )
  114. )
  115. for code in inset['codes']:
  116. main_codes.remove(code)
  117. insetBbox = self.renderMapInset(data_source, main_codes, self.left, self.top, self.width)
  118. insetHeight = (insetBbox[3] - insetBbox[1]) * (self.width / (insetBbox[2] - insetBbox[0]))
  119. envelope.append( shapely.geometry.box( self.left, self.top, self.left + self.width, self.top + insetHeight ) )
  120. mapBbox = shapely.geometry.MultiPolygon( envelope ).bounds
  121. self.map.width = mapBbox[2] + mapBbox[0]
  122. self.map.height = mapBbox[3] + mapBbox[1]
  123. self.map.insets.append({
  124. "bbox": [{"x": insetBbox[0], "y": -insetBbox[3]}, {"x": insetBbox[2], "y": -insetBbox[1]}],
  125. "left": self.left,
  126. "top": self.top,
  127. "width": self.width,
  128. "height": insetHeight
  129. })
  130. self.map.projection = {"type": self.projection, "centralMeridian": float(self.longitude0)}
  131. open(output_file, 'w').write( self.map.getJSCode() )
  132. if self.for_each is not None:
  133. for code in codes:
  134. childConfig = copy.deepcopy(self.for_each)
  135. for param in ('input_file', 'output_file', 'where', 'name'):
  136. childConfig[param] = childConfig[param].replace('{{code}}', code.lower())
  137. converter = Converter(childConfig)
  138. converter.convert(childConfig['output_file'])
  139. def renderMapInset(self, data_source, codes, left, top, width):
  140. envelope = []
  141. geometries = filter(lambda g: g.properties[self.config['code_field']] in codes, data_source.geometries)
  142. for geometry in geometries:
  143. envelope.append( geometry.geom.envelope )
  144. bbox = shapely.geometry.MultiPolygon( envelope ).bounds
  145. scale = (bbox[2]-bbox[0]) / width
  146. # generate SVG paths
  147. for geometry in geometries:
  148. geom = geometry.geom
  149. if self.buffer_distance:
  150. geom = geom.buffer(self.buffer_distance*scale, 1)
  151. if geom.is_empty:
  152. continue
  153. if self.simplify_tolerance:
  154. geom = geom.simplify(self.simplify_tolerance*scale, preserve_topology=True)
  155. if isinstance(geom, shapely.geometry.multipolygon.MultiPolygon):
  156. polygons = geom.geoms
  157. else:
  158. polygons = [geom]
  159. path = ''
  160. for polygon in polygons:
  161. rings = []
  162. rings.append(polygon.exterior)
  163. rings.extend(polygon.interiors)
  164. for ring in rings:
  165. for pointIndex in range( len(ring.coords) ):
  166. point = ring.coords[pointIndex]
  167. if pointIndex == 0:
  168. path += 'M'+str( round( (point[0]-bbox[0]) / scale + left, self.precision) )
  169. path += ','+str( round( (bbox[3] - point[1]) / scale + top, self.precision) )
  170. else:
  171. path += 'l' + str( round(point[0]/scale - ring.coords[pointIndex-1][0]/scale, self.precision) )
  172. path += ',' + str( round(ring.coords[pointIndex-1][1]/scale - point[1]/scale, self.precision) )
  173. path += 'Z'
  174. self.map.addPath(path, geometry.properties[self.config['code_field']].lower(), geometry.properties[self.config['name_field']])
  175. return bbox
  176. class Geometry:
  177. def __init__(self, geometry, properties):
  178. self.geom = geometry
  179. self.properties = properties
  180. class GeometryProperty(Variable):
  181. operations = set(["equality", "membership"])
  182. def __init__(self, name):
  183. self.name = name
  184. def equals(self, value, context):
  185. return context[self.name] == value
  186. def belongs_to(self, value, context):
  187. return value in context[self.name]
  188. def is_subset(self, value, context):
  189. return set(value).issubset(set(context[self.name]))
  190. def to_python(self, value):
  191. return unicode(value[self.name])
  192. class DataSource:
  193. def __init__(self, config):
  194. default_config = {
  195. "projection": "merc",
  196. "longitude0": 0
  197. }
  198. default_config.update(config)
  199. self.config = default_config
  200. self.spatialRef = osr.SpatialReference()
  201. projString = '+proj='+str(self.config['projection'])+' +a=6381372 +b=6381372 +lat_0=0'
  202. #if 'emulate_longitude0' in self.config and not self.config['emulate_longitude0']:
  203. projString += ' +lon_0='+str(self.config['longitude0'])
  204. self.spatialRef.ImportFromProj4(projString)
  205. def load_data(self):
  206. self.source = ogr.Open( self.config['file_name'], update = 0 )
  207. self.layer = self.source.GetLayer(0)
  208. if 'filter' in self.config and self.config['filter'] is not None:
  209. self.layer.SetAttributeFilter( self.config['filter'].encode('ascii') )
  210. self.layer_dfn = self.layer.GetLayerDefn()
  211. self.fields = []
  212. field_count = self.layer_dfn.GetFieldCount()
  213. for field_index in range(field_count):
  214. field = self.layer_dfn.GetFieldDefn( field_index )
  215. self.fields.append({
  216. 'name': field.GetName(),
  217. 'type': field.GetType(),
  218. 'width': field.GetWidth(),
  219. 'precision': field.GetPrecision()
  220. })
  221. self.geometries = []
  222. for feature in self.layer:
  223. geometry = feature.GetGeometryRef()
  224. geometry.TransformTo( self.spatialRef )
  225. geometry = shapely.wkb.loads( geometry.ExportToWkb() )
  226. if not geometry.is_valid:
  227. geometry = geometry.buffer(0)
  228. properties = {}
  229. for field in self.fields:
  230. properties[field['name']] = feature.GetFieldAsString(field['name']).decode('utf-8')
  231. self.geometries.append( Geometry(geometry, properties) )
  232. self.layer.ResetReading()
  233. self.create_grammar()
  234. def create_grammar(self):
  235. root_table = SymbolTable("root",
  236. map( lambda f: Bind(f['name'], GeometryProperty(f['name'])), self.fields )
  237. )
  238. tokens = {
  239. 'not': 'not',
  240. 'eq': '==',
  241. 'ne': '!=',
  242. 'belongs_to': 'in',
  243. 'is_subset': 'are included in',
  244. 'or': "or",
  245. 'and': 'and'
  246. }
  247. grammar = Grammar(**tokens)
  248. self.parse_manager = EvaluableParseManager(root_table, grammar)
  249. def output(self, output):
  250. if output.get('format') == 'jqvmap':
  251. self.output_jvm(output)
  252. else:
  253. self.output_ogr(output)
  254. def output_ogr(self, output):
  255. driver = ogr.GetDriverByName( 'ESRI Shapefile' )
  256. if os.path.exists( output['file_name'] ):
  257. driver.DeleteDataSource( output['file_name'] )
  258. source = driver.CreateDataSource( output['file_name'] )
  259. layer = source.CreateLayer( self.layer_dfn.GetName(),
  260. geom_type = self.layer_dfn.GetGeomType(),
  261. srs = self.layer.GetSpatialRef() )
  262. for field in self.fields:
  263. fd = ogr.FieldDefn( str(field['name']), field['type'] )
  264. fd.SetWidth( field['width'] )
  265. if 'precision' in field:
  266. fd.SetPrecision( field['precision'] )
  267. layer.CreateField( fd )
  268. for geometry in self.geometries:
  269. if geometry.geom is not None:
  270. feature = ogr.Feature( feature_def = layer.GetLayerDefn() )
  271. for index, field in enumerate(self.fields):
  272. if field['name'] in geometry.properties:
  273. feature.SetField( index, geometry.properties[field['name']].encode('utf-8') )
  274. else:
  275. feature.SetField( index, '' )
  276. feature.SetGeometryDirectly(
  277. ogr.CreateGeometryFromWkb(
  278. shapely.wkb.dumps(
  279. geometry.geom
  280. )
  281. )
  282. )
  283. layer.CreateFeature( feature )
  284. feature.Destroy()
  285. source.Destroy()
  286. def output_jvm(self, output):
  287. params = copy.deepcopy(output['params'])
  288. params.update({
  289. "projection": self.config["projection"],
  290. "longitude0": self.config["longitude0"]
  291. })
  292. converter = Converter(params)
  293. converter.convert(self, output['file_name'])
  294. class PolygonSimplifier:
  295. def __init__(self, geometries):
  296. self.format = '%.8f %.8f'
  297. self.tolerance = 0.05
  298. self.geometries = geometries
  299. connections = {}
  300. counter = 0
  301. for geom in geometries:
  302. counter += 1
  303. polygons = []
  304. if isinstance(geom, shapely.geometry.Polygon):
  305. polygons.append(geom)
  306. else:
  307. for polygon in geom:
  308. polygons.append(polygon)
  309. for polygon in polygons:
  310. if polygon.area > 0:
  311. lines = []
  312. lines.append(polygon.exterior)
  313. for line in polygon.interiors:
  314. lines.append(line)
  315. for line in lines:
  316. for i in range(len(line.coords)-1):
  317. indexFrom = i
  318. indexTo = i+1
  319. pointFrom = self.format % line.coords[indexFrom]
  320. pointTo = self.format % line.coords[indexTo]
  321. if pointFrom == pointTo:
  322. continue
  323. if not (pointFrom in connections):
  324. connections[pointFrom] = {}
  325. connections[pointFrom][pointTo] = 1
  326. if not (pointTo in connections):
  327. connections[pointTo] = {}
  328. connections[pointTo][pointFrom] = 1
  329. self.connections = connections
  330. self.simplifiedLines = {}
  331. self.pivotPoints = {}
  332. def simplifyRing(self, ring):
  333. coords = list(ring.coords)[0:-1]
  334. simpleCoords = []
  335. isPivot = False
  336. pointIndex = 0
  337. while not isPivot and pointIndex < len(coords):
  338. pointStr = self.format % coords[pointIndex]
  339. pointIndex += 1
  340. isPivot = ((len(self.connections[pointStr]) > 2) or (pointStr in self.pivotPoints))
  341. pointIndex = pointIndex - 1
  342. if not isPivot:
  343. simpleRing = shapely.geometry.LineString(coords).simplify(self.tolerance)
  344. if len(simpleRing.coords) <= 2:
  345. return None
  346. else:
  347. self.pivotPoints[self.format % coords[0]] = True
  348. self.pivotPoints[self.format % coords[-1]] = True
  349. simpleLineKey = self.format % coords[0]+':'+self.format % coords[1]+':'+self.format % coords[-1]
  350. self.simplifiedLines[simpleLineKey] = simpleRing.coords
  351. return simpleRing
  352. else:
  353. points = coords[pointIndex:len(coords)]
  354. points.extend(coords[0:pointIndex+1])
  355. iFrom = 0
  356. for i in range(1, len(points)):
  357. pointStr = self.format % points[i]
  358. if ((len(self.connections[pointStr]) > 2) or (pointStr in self.pivotPoints)):
  359. line = points[iFrom:i+1]
  360. lineKey = self.format % line[-1]+':'+self.format % line[-2]+':'+self.format % line[0]
  361. if lineKey in self.simplifiedLines:
  362. simpleLine = self.simplifiedLines[lineKey]
  363. simpleLine = list(reversed(simpleLine))
  364. else:
  365. simpleLine = shapely.geometry.LineString(line).simplify(self.tolerance).coords
  366. lineKey = self.format % line[0]+':'+self.format % line[1]+':'+self.format % line[-1]
  367. self.simplifiedLines[lineKey] = simpleLine
  368. simpleCoords.extend( simpleLine[0:-1] )
  369. iFrom = i
  370. if len(simpleCoords) <= 2:
  371. return None
  372. else:
  373. return shapely.geometry.LineString(simpleCoords)
  374. def simplifyPolygon(self, polygon):
  375. simpleExtRing = self.simplifyRing(polygon.exterior)
  376. if simpleExtRing is None:
  377. return None
  378. simpleIntRings = []
  379. for ring in polygon.interiors:
  380. simpleIntRing = self.simplifyRing(ring)
  381. if simpleIntRing is not None:
  382. simpleIntRings.append(simpleIntRing)
  383. return shapely.geometry.Polygon(simpleExtRing, simpleIntRings)
  384. def simplify(self):
  385. results = []
  386. for geom in self.geometries:
  387. polygons = []
  388. simplePolygons = []
  389. if isinstance(geom, shapely.geometry.Polygon):
  390. polygons.append(geom)
  391. else:
  392. for polygon in geom:
  393. polygons.append(polygon)
  394. for polygon in polygons:
  395. simplePolygon = self.simplifyPolygon(polygon)
  396. if not (simplePolygon is None or simplePolygon._geom is None):
  397. simplePolygons.append(simplePolygon)
  398. if len(simplePolygons) > 0:
  399. results.append(shapely.geometry.MultiPolygon(simplePolygons))
  400. else:
  401. results.append(None)
  402. return results
  403. class Processor:
  404. def __init__(self, config):
  405. self.config = config
  406. def process(self):
  407. self.data_sources = {}
  408. for action in self.config:
  409. getattr(self, action['name'])( action, self.data_sources.get(".") )
  410. def read_data(self, config, data_source):
  411. self.data_sources["."] = DataSource( config )
  412. self.data_sources["."].load_data()
  413. def write_data(self, config, data_source):
  414. data_source.output( config )
  415. def union(self, config, data_source):
  416. groups = {}
  417. geometries = []
  418. for geometry in data_source.geometries:
  419. if geometry.properties[config['by']] in groups:
  420. groups[geometry.properties[config['by']]]['geoms'].append(geometry.geom)
  421. else:
  422. groups[geometry.properties[config['by']]] = {
  423. 'geoms': [geometry.geom],
  424. 'properties': geometry.properties
  425. }
  426. for key in groups:
  427. geometries.append( Geometry(shapely.ops.cascaded_union( groups[key]['geoms'] ), groups[key]['properties']) )
  428. data_source.geometries = geometries
  429. def merge(self, config, data_source):
  430. new_geometries = []
  431. for rule in config['rules']:
  432. expression = data_source.parse_manager.parse( rule['where'] )
  433. geometries = filter(lambda g: expression(g.properties), data_source.geometries)
  434. geometries = map(lambda g: g.geom, geometries)
  435. new_geometries.append( Geometry(shapely.ops.cascaded_union( geometries ), rule['fields']) )
  436. data_source.fields = config['fields']
  437. data_source.geometries = new_geometries
  438. def join_data(self, config, data_source):
  439. field_names = [f['name'] for f in config['fields']]
  440. if 'data' in config:
  441. data_col = config['data']
  442. else:
  443. data_file = open(config['file_name'], 'rb')
  444. data_col = csv.reader(data_file, delimiter='\t', quotechar='"')
  445. data = {}
  446. for row in data_col:
  447. row_dict = dict(zip(field_names, row))
  448. data[row_dict.pop(config['on'])] = row_dict
  449. for geometry in data_source.geometries:
  450. if geometry.properties[config['on']] in data:
  451. geometry.properties.update( data[geometry.properties[config['on']]] )
  452. field_names = map(lambda f: f['name'], data_source.fields)
  453. data_source.fields = data_source.fields + filter(lambda f: f['name'] not in field_names, config['fields'])
  454. def remove(self, config, data_source):
  455. expression = data_source.parse_manager.parse( config['where'] )
  456. data_source.geometries = filter(lambda g: not expression(g.properties), data_source.geometries)
  457. def remove_fields(self, config, data_source):
  458. data_source.fields = filter(lambda f: f.name not in config['fields'], data_source.fields)
  459. def remove_other_fields(self, config, data_source):
  460. data_source.fields = filter(lambda f: f['name'] in config['fields'], data_source.fields)
  461. def buffer(self, config, data_source):
  462. for geometry in data_source.geometries:
  463. geometry.geom = geometry.geom.buffer(config['distance'], config['resolution'])
  464. def simplify_adjancent_polygons(self, config, data_source):
  465. simple_geometries = PolygonSimplifier( map( lambda g: g.geom, data_source.geometries ) ).simplify()
  466. for i in range(len(data_source.geometries)):
  467. data_source.geometries[i].geom = simple_geometries[i]
  468. def intersect_rect(self, config, data_source):
  469. transform = osr.CoordinateTransformation( data_source.layer.GetSpatialRef(), data_source.spatialRef )
  470. point1 = transform.TransformPoint(config['rect'][0], config['rect'][1])
  471. point2 = transform.TransformPoint(config['rect'][2], config['rect'][3])
  472. rect = shapely.geometry.box(point1[0], point1[1], point2[0], point2[1])
  473. for geometry in data_source.geometries:
  474. geometry.geom = geometry.geom.intersection(rect)
  475. def remove_small_polygons(self, config, data_source):
  476. for geometry in data_source.geometries:
  477. if isinstance(geometry.geom, shapely.geometry.multipolygon.MultiPolygon):
  478. polygons = geometry.geom.geoms
  479. else:
  480. polygons = [geometry.geom]
  481. polygons = filter(lambda p: p.area > config['minimal_area'], polygons)
  482. if len(polygons) > 0:
  483. geometry.geom = shapely.geometry.multipolygon.MultiPolygon(polygons)
  484. args = {}
  485. if len(sys.argv) > 1:
  486. paramsJson = open(sys.argv[1], 'r').read()
  487. else:
  488. paramsJson = sys.stdin.read()
  489. paramsJson = json.loads(paramsJson)
  490. processor = Processor(paramsJson)
  491. processor.process()