lavkach3

1
import sys
2
import json
3
import csv
4
import shapely.wkb
5
import shapely.geometry
6
import shapely.ops
7
import codecs
8
import os
9
import inspect
10
import copy
11
from osgeo import ogr
12
from osgeo import osr
13
from booleano.parser import Grammar, EvaluableParseManager, SymbolTable, Bind
14
from booleano.operations import Variable
15

16

17
class Map:
18
  def __init__(self, name, language):
19
    self.paths = {}
20
    self.name = name
21
    self.language = language
22
    self.width = 0
23
    self.height = 0
24
    self.bbox = []
25

26
  def addPath(self, path, code, name):
27
    self.paths[code] = {"path": path, "name": name}
28

29
  def getJSCode(self):
30
    map = {"paths": self.paths, "width": self.width, "height": self.height, "insets": self.insets, "projection": self.projection}
31
    return "jQuery.fn.vectorMap('addMap', '"+self.name+"_"+self.projection['type']+"',"+json.dumps(map)+');'
32

33

34
class Converter:
35
  def __init__(self, config):
36
    args = {
37
      'buffer_distance': -0.4,
38
      'simplify_tolerance': 0.2,
39
      'longitude0': 0,
40
      'projection': 'mill',
41
      'name': 'world',
42
      'width': 900,
43
      'left': 0,
44
      'top': 0,
45
      'language': 'en',
46
      'precision': 2,
47
      'insets': []
48
    }
49
    args.update(config)
50

51
    self.config = args
52

53
    self.map = Map(args['name'], args.get('language'))
54

55
    if args.get('sources'):
56
      self.sources = args['sources']
57
    else:
58
      self.sources = [{
59
        'input_file': args.get('input_file'),
60
        'where': args.get('where'),
61
        'name_field': args.get('name_field'),
62
        'code_field': args.get('code_field'),
63
        'input_file_encoding': args.get('input_file_encoding')
64
      }]
65

66
    default_source = {
67
      'where': '',
68
      'name_field': 0,
69
      'code_field': 1,
70
      'input_file_encoding': 'iso-8859-1'
71
    }
72

73
    for index in range(len(self.sources)):
74
      for key in default_source:
75
        if self.sources[index].get(key) is None:
76
          self.sources[index][key] = default_source[key]
77

78
    self.features = {}
79
    self.width = args.get('width')
80
    self.left = args.get('left')
81
    self.top = args.get('top')
82
    self.minimal_area = args.get('minimal_area')
83
    self.longitude0 = float(args.get('longitude0'))
84
    self.projection = args.get('projection')
85
    self.precision = args.get('precision')
86
    self.buffer_distance = args.get('buffer_distance')
87
    self.simplify_tolerance = args.get('simplify_tolerance')
88
    self.for_each = args.get('for_each')
89
    self.emulate_longitude0 = args.get('emulate_longitude0')
90
    if args.get('emulate_longitude0') is None and (self.projection == 'merc' or self.projection =='mill') and self.longitude0 != 0:
91
      self.emulate_longitude0 = True
92

93
    if args.get('viewport'):
94
      self.viewport = map(lambda s: float(s), args.get('viewport').split(' '))
95
    else:
96
      self.viewport = False
97

98
    # spatial reference to convert to
99
    self.spatialRef = osr.SpatialReference()
100
    projString = '+proj='+str(self.projection)+' +a=6381372 +b=6381372 +lat_0=0'
101
    if not self.emulate_longitude0:
102
      projString += ' +lon_0='+str(self.longitude0)
103
    self.spatialRef.ImportFromProj4(projString)
104

105
    # handle map insets
106
    if args.get('insets'):
107
      self.insets = args.get('insets')
108
    else:
109
      self.insets = []
110

111

112
  def convert(self, data_source, output_file):
113
    codes = map(lambda g: g.properties[self.config['code_field']], data_source.geometries)
114
    main_codes = copy.copy(codes)
115
    self.map.insets = []
116
    envelope = []
117
    for inset in self.insets:
118
      insetBbox = self.renderMapInset(data_source, inset['codes'], inset['left'], inset['top'], inset['width'])
119
      insetHeight = (insetBbox[3] - insetBbox[1]) * (inset['width'] / (insetBbox[2] - insetBbox[0]))
120
      self.map.insets.append({
121
        "bbox": [{"x": insetBbox[0], "y": -insetBbox[3]}, {"x": insetBbox[2], "y": -insetBbox[1]}],
122
        "left": inset['left'],
123
        "top": inset['top'],
124
        "width": inset['width'],
125
        "height": insetHeight
126
      })
127
      envelope.append(
128
        shapely.geometry.box(
129
          inset['left'], inset['top'], inset['left'] + inset['width'], inset['top'] + insetHeight
130
        )
131
      )
132
      for code in inset['codes']:
133
        main_codes.remove(code)
134

135
    insetBbox = self.renderMapInset(data_source, main_codes, self.left, self.top, self.width)
136
    insetHeight = (insetBbox[3] - insetBbox[1]) * (self.width / (insetBbox[2] - insetBbox[0]))
137
    envelope.append( shapely.geometry.box( self.left, self.top, self.left + self.width, self.top + insetHeight ) )
138
    mapBbox = shapely.geometry.MultiPolygon( envelope ).bounds
139

140
    self.map.width = mapBbox[2] + mapBbox[0]
141
    self.map.height = mapBbox[3] + mapBbox[1]
142
    self.map.insets.append({
143
      "bbox": [{"x": insetBbox[0], "y": -insetBbox[3]}, {"x": insetBbox[2], "y": -insetBbox[1]}],
144
      "left": self.left,
145
      "top": self.top,
146
      "width": self.width,
147
      "height": insetHeight
148
    })
149
    self.map.projection = {"type": self.projection, "centralMeridian": float(self.longitude0)}
150

151
    open(output_file, 'w').write( self.map.getJSCode() )
152

153
    if self.for_each is not None:
154
      for code in codes:
155
        childConfig = copy.deepcopy(self.for_each)
156
        for param in ('input_file', 'output_file', 'where', 'name'):
157
          childConfig[param] = childConfig[param].replace('{{code}}', code.lower())
158
        converter = Converter(childConfig)
159
        converter.convert(childConfig['output_file'])
160

161
  def renderMapInset(self, data_source, codes, left, top, width):
162
    envelope = []
163
    geometries = filter(lambda g: g.properties[self.config['code_field']] in codes, data_source.geometries)
164
    for geometry in geometries:
165
      envelope.append( geometry.geom.envelope )
166

167
    bbox = shapely.geometry.MultiPolygon( envelope ).bounds
168

169
    scale = (bbox[2]-bbox[0]) / width
170

171
    # generate SVG paths
172
    for geometry in geometries:
173
      geom = geometry.geom
174
      if self.buffer_distance:
175
        geom = geom.buffer(self.buffer_distance*scale, 1)
176
      if geom.is_empty:
177
        continue
178
      if self.simplify_tolerance:
179
        geom = geom.simplify(self.simplify_tolerance*scale, preserve_topology=True)
180
      if isinstance(geom, shapely.geometry.multipolygon.MultiPolygon):
181
        polygons = geom.geoms
182
      else:
183
        polygons = [geom]
184
      path = ''
185
      for polygon in polygons:
186
        rings = []
187
        rings.append(polygon.exterior)
188
        rings.extend(polygon.interiors)
189
        for ring in rings:
190
          for pointIndex in range( len(ring.coords) ):
191
            point = ring.coords[pointIndex]
192
            if pointIndex == 0:
193
              path += 'M'+str( round( (point[0]-bbox[0]) / scale + left, self.precision) )
194
              path += ','+str( round( (bbox[3] - point[1]) / scale + top, self.precision) )
195
            else:
196
              path += 'l' + str( round(point[0]/scale - ring.coords[pointIndex-1][0]/scale, self.precision) )
197
              path += ',' + str( round(ring.coords[pointIndex-1][1]/scale - point[1]/scale, self.precision) )
198
          path += 'Z'
199
      self.map.addPath(path, geometry.properties[self.config['code_field']], geometry.properties[self.config['name_field']])
200
    return bbox
201

202

203
class Geometry:
204
  def __init__(self, geometry, properties):
205
    self.geom = geometry
206
    self.properties = properties
207

208

209
class GeometryProperty(Variable):
210
  operations = set(["equality", "membership"])
211

212
  def __init__(self, name):
213
    self.name = name
214

215
  def equals(self, value, context):
216
    return context[self.name] == value
217

218
  def belongs_to(self, value, context):
219
    return value in context[self.name]
220

221
  def is_subset(self, value, context):
222
    return set(value).issubset(set(context[self.name]))
223

224
  def to_python(self, value):
225
    return unicode(value[self.name])
226

227

228
class DataSource:
229
  def __init__(self, config):
230
    default_config = {
231
      "projection": "merc",
232
      "longitude0": 0
233
    }
234
    default_config.update(config)
235
    self.config = default_config
236

237
    self.spatialRef = osr.SpatialReference()
238
    projString = '+proj='+str(self.config['projection'])+' +a=6381372 +b=6381372 +lat_0=0'
239
    #if 'emulate_longitude0' in self.config and not self.config['emulate_longitude0']:
240
    projString += ' +lon_0='+str(self.config['longitude0'])
241
    self.spatialRef.ImportFromProj4(projString)
242

243
  def load_data(self):
244
    self.source = ogr.Open( self.config['file_name'], update = 0 )
245
    self.layer = self.source.GetLayer(0)
246
    if 'filter' in self.config and self.config['filter'] is not None:
247
      self.layer.SetAttributeFilter( self.config['filter'].encode('ascii') )
248
    self.layer_dfn = self.layer.GetLayerDefn()
249

250
    self.fields = []
251
    field_count = self.layer_dfn.GetFieldCount()
252
    for field_index in range(field_count):
253
      field = self.layer_dfn.GetFieldDefn( field_index )
254
      self.fields.append({
255
        'name': field.GetName(),
256
        'type': field.GetType(),
257
        'width': field.GetWidth(),
258
        'precision': field.GetPrecision()
259
      })
260

261
    self.geometries = []
262

263
    for feature in self.layer:
264
      geometry = feature.GetGeometryRef()
265
      geometry.TransformTo( self.spatialRef )
266
      geometry = shapely.wkb.loads( geometry.ExportToWkb() )
267
      if not geometry.is_valid:
268
        geometry = geometry.buffer(0)
269
      properties = {}
270
      for field in self.fields:
271
        properties[field['name']] = feature.GetFieldAsString(field['name']).decode('utf-8')
272
      self.geometries.append( Geometry(geometry, properties) )
273

274
    self.layer.ResetReading()
275

276
    self.create_grammar()
277

278
  def create_grammar(self):
279
    root_table = SymbolTable("root",
280
      map( lambda f: Bind(f['name'], GeometryProperty(f['name'])), self.fields )
281
    )
282

283
    tokens = {
284
      'not': 'not',
285
      'eq': '==',
286
      'ne': '!=',
287
      'belongs_to': 'in',
288
      'is_subset': 'are included in',
289
      'or': "or",
290
      'and': 'and'
291
    }
292
    grammar = Grammar(**tokens)
293
    self.parse_manager = EvaluableParseManager(root_table, grammar)
294

295
  def output(self, output):
296
    if output.get('format') == 'jvectormap':
297
      self.output_jvm(output)
298
    else:
299
      self.output_ogr(output)
300

301
  def output_ogr(self, output):
302
    driver = ogr.GetDriverByName( 'ESRI Shapefile' )
303
    if os.path.exists( output['file_name'] ):
304
      driver.DeleteDataSource( output['file_name'] )
305
    source = driver.CreateDataSource( output['file_name'] )
306
    layer = source.CreateLayer( self.layer_dfn.GetName(),
307
                                geom_type = self.layer_dfn.GetGeomType(),
308
                                srs = self.layer.GetSpatialRef() )
309

310
    for field in self.fields:
311
      fd = ogr.FieldDefn( str(field['name']), field['type'] )
312
      fd.SetWidth( field['width'] )
313
      if 'precision' in field:
314
        fd.SetPrecision( field['precision'] )
315
      layer.CreateField( fd )
316

317
    for geometry in self.geometries:
318
      if geometry.geom is not None:
319
        feature = ogr.Feature( feature_def = layer.GetLayerDefn() )
320
        for index, field in enumerate(self.fields):
321
          if field['name'] in geometry.properties:
322
            feature.SetField( index, geometry.properties[field['name']].encode('utf-8') )
323
          else:
324
            feature.SetField( index, '' )
325
        feature.SetGeometryDirectly(
326
          ogr.CreateGeometryFromWkb(
327
            shapely.wkb.dumps(
328
              geometry.geom
329
            )
330
          )
331
        )
332
        layer.CreateFeature( feature )
333
        feature.Destroy()
334

335
    source.Destroy()
336

337
  def output_jvm(self, output):
338
    params = copy.deepcopy(output['params'])
339
    params.update({
340
      "projection": self.config["projection"],
341
      "longitude0": self.config["longitude0"]
342
    })
343
    converter = Converter(params)
344
    converter.convert(self, output['file_name'])
345

346
class PolygonSimplifier:
347
  def __init__(self, geometries):
348
    self.format = '%.8f %.8f'
349
    self.tolerance = 0.05
350
    self.geometries = geometries
351

352
    connections = {}
353
    counter = 0
354
    for geom in geometries:
355
      counter += 1
356
      polygons = []
357

358
      if isinstance(geom, shapely.geometry.Polygon):
359
        polygons.append(geom)
360
      else:
361
        for polygon in geom:
362
          polygons.append(polygon)
363

364
      for polygon in polygons:
365
        if polygon.area > 0:
366
          lines = []
367
          lines.append(polygon.exterior)
368
          for line in polygon.interiors:
369
            lines.append(line)
370

371
          for line in lines:
372
            for i in range(len(line.coords)-1):
373
              indexFrom = i
374
              indexTo = i+1
375
              pointFrom = self.format % line.coords[indexFrom]
376
              pointTo = self.format % line.coords[indexTo]
377
              if pointFrom == pointTo:
378
                continue
379
              if not (pointFrom in connections):
380
                connections[pointFrom] = {}
381
              connections[pointFrom][pointTo] = 1
382
              if not (pointTo in connections):
383
                connections[pointTo] = {}
384
              connections[pointTo][pointFrom] = 1
385
    self.connections = connections
386
    self.simplifiedLines = {}
387
    self.pivotPoints = {}
388

389
  def simplifyRing(self, ring):
390
    coords = list(ring.coords)[0:-1]
391
    simpleCoords = []
392

393
    isPivot = False
394
    pointIndex = 0
395
    while not isPivot and pointIndex < len(coords):
396
      pointStr = self.format % coords[pointIndex]
397
      pointIndex += 1
398
      isPivot = ((len(self.connections[pointStr]) > 2) or (pointStr in self.pivotPoints))
399
    pointIndex = pointIndex - 1
400

401
    if not isPivot:
402
      simpleRing = shapely.geometry.LineString(coords).simplify(self.tolerance)
403
      if len(simpleRing.coords) <= 2:
404
        return None
405
      else:
406
        self.pivotPoints[self.format % coords[0]] = True
407
        self.pivotPoints[self.format % coords[-1]] = True
408
        simpleLineKey = self.format % coords[0]+':'+self.format % coords[1]+':'+self.format % coords[-1]
409
        self.simplifiedLines[simpleLineKey] = simpleRing.coords
410
        return simpleRing
411
    else:
412
      points = coords[pointIndex:len(coords)]
413
      points.extend(coords[0:pointIndex+1])
414
      iFrom = 0
415
      for i in range(1, len(points)):
416
        pointStr = self.format % points[i]
417
        if ((len(self.connections[pointStr]) > 2) or (pointStr in self.pivotPoints)):
418
          line = points[iFrom:i+1]
419
          lineKey = self.format % line[-1]+':'+self.format % line[-2]+':'+self.format % line[0]
420
          if lineKey in self.simplifiedLines:
421
            simpleLine = self.simplifiedLines[lineKey]
422
            simpleLine = list(reversed(simpleLine))
423
          else:
424
            simpleLine = shapely.geometry.LineString(line).simplify(self.tolerance).coords
425
            lineKey = self.format % line[0]+':'+self.format % line[1]+':'+self.format % line[-1]
426
            self.simplifiedLines[lineKey] = simpleLine
427
          simpleCoords.extend( simpleLine[0:-1] )
428
          iFrom = i
429
      if len(simpleCoords) <= 2:
430
        return None
431
      else:
432
        return shapely.geometry.LineString(simpleCoords)
433

434
  def simplifyPolygon(self, polygon):
435
    simpleExtRing = self.simplifyRing(polygon.exterior)
436
    if simpleExtRing is None:
437
      return None
438
    simpleIntRings = []
439
    for ring in polygon.interiors:
440
      simpleIntRing = self.simplifyRing(ring)
441
      if simpleIntRing is not None:
442
        simpleIntRings.append(simpleIntRing)
443
    return shapely.geometry.Polygon(simpleExtRing, simpleIntRings)
444

445
  def simplify(self):
446
    results = []
447
    for geom in self.geometries:
448
      polygons = []
449
      simplePolygons = []
450

451
      if isinstance(geom, shapely.geometry.Polygon):
452
        polygons.append(geom)
453
      else:
454
        for polygon in geom:
455
          polygons.append(polygon)
456

457
      for polygon in polygons:
458
        simplePolygon = self.simplifyPolygon(polygon)
459
        if not (simplePolygon is None or simplePolygon._geom is None):
460
          simplePolygons.append(simplePolygon)
461

462
      if len(simplePolygons) > 0:
463
        results.append(shapely.geometry.MultiPolygon(simplePolygons))
464
      else:
465
        results.append(None)
466
    return results
467

468

469
class Processor:
470
  def __init__(self, config):
471
    self.config = config
472

473
  def process(self):
474
    self.data_sources = {}
475
    for action in self.config:
476
      getattr(self, action['name'])( action, self.data_sources.get(".") )
477

478
  def read_data(self, config, data_source):
479
    self.data_sources["."] = DataSource( config )
480
    self.data_sources["."].load_data()
481

482
  def write_data(self, config, data_source):
483
    data_source.output( config )
484

485
  def union(self, config, data_source):
486
    groups = {}
487
    geometries = []
488
    for geometry in data_source.geometries:
489
      if geometry.properties[config['by']] in groups:
490
        groups[geometry.properties[config['by']]]['geoms'].append(geometry.geom)
491
      else:
492
        groups[geometry.properties[config['by']]] = {
493
          'geoms': [geometry.geom],
494
          'properties': geometry.properties
495
        }
496
    for key in groups:
497
      geometries.append( Geometry(shapely.ops.cascaded_union( groups[key]['geoms'] ), groups[key]['properties']) )
498
    data_source.geometries = geometries
499

500
  def merge(self, config, data_source):
501
    new_geometries = []
502
    for rule in config['rules']:
503
      expression = data_source.parse_manager.parse( rule['where'] )
504
      geometries = filter(lambda g: expression(g.properties), data_source.geometries)
505
      geometries = map(lambda g: g.geom, geometries)
506
      new_geometries.append( Geometry(shapely.ops.cascaded_union( geometries ), rule['fields']) )
507
    data_source.fields = config['fields']
508
    data_source.geometries = new_geometries
509

510
  def join_data(self, config, data_source):
511
    field_names = [f['name'] for f in config['fields']]
512
    if 'data' in config:
513
      data_col = config['data']
514
    else:
515
      data_file = open(config['file_name'], 'rb')
516
      data_col = csv.reader(data_file, delimiter='\t', quotechar='"')
517
    data = {}
518
    for row in data_col:
519
      row_dict = dict(zip(field_names, row))
520
      data[row_dict.pop(config['on'])] = row_dict
521
    for geometry in data_source.geometries:
522
      if geometry.properties[config['on']] in data:
523
        geometry.properties.update( data[geometry.properties[config['on']]] )
524
    field_names = map(lambda f: f['name'], data_source.fields)
525
    data_source.fields = data_source.fields + filter(lambda f: f['name'] not in field_names, config['fields'])
526

527
  def remove(self, config, data_source):
528
    expression = data_source.parse_manager.parse( config['where'] )
529
    data_source.geometries = filter(lambda g: not expression(g.properties), data_source.geometries)
530

531
  def remove_fields(self, config, data_source):
532
    data_source.fields = filter(lambda f: f.name not in config['fields'], data_source.fields)
533

534
  def remove_other_fields(self, config, data_source):
535
    data_source.fields = filter(lambda f: f['name'] in config['fields'], data_source.fields)
536

537
  def buffer(self, config, data_source):
538
    for geometry in data_source.geometries:
539
      geometry.geom = geometry.geom.buffer(config['distance'], config['resolution'])
540

541
  def simplify_adjancent_polygons(self, config, data_source):
542
    simple_geometries = PolygonSimplifier( map( lambda g: g.geom, data_source.geometries ) ).simplify()
543
    for i in range(len(data_source.geometries)):
544
      data_source.geometries[i].geom = simple_geometries[i]
545

546
  def intersect_rect(self, config, data_source):
547
    transform = osr.CoordinateTransformation( data_source.layer.GetSpatialRef(), data_source.spatialRef )
548
    point1 = transform.TransformPoint(config['rect'][0], config['rect'][1])
549
    point2 = transform.TransformPoint(config['rect'][2], config['rect'][3])
550
    rect = shapely.geometry.box(point1[0], point1[1], point2[0], point2[1])
551
    for geometry in data_source.geometries:
552
      geometry.geom = geometry.geom.intersection(rect)
553

554
  def remove_small_polygons(self, config, data_source):
555
    for geometry in data_source.geometries:
556
      if isinstance(geometry.geom, shapely.geometry.multipolygon.MultiPolygon):
557
        polygons = geometry.geom.geoms
558
      else:
559
        polygons = [geometry.geom]
560
      polygons = filter(lambda p: p.area > config['minimal_area'], polygons)
561
      if len(polygons) > 0:
562
        geometry.geom = shapely.geometry.multipolygon.MultiPolygon(polygons)
563

564

565
args = {}
566
if len(sys.argv) > 1:
567
  paramsJson = open(sys.argv[1], 'r').read()
568
else:
569
  paramsJson = sys.stdin.read()
570
paramsJson = json.loads(paramsJson)
571

572
processor = Processor(paramsJson)
573
processor.process()
574