← Back to branch summary

~awuerl/blitzortung-python/master

« back to all changes in this revision

Viewing changes to blitzortung/geom.py

  • Committer: Andreas Würl
  • Date: 2012-01-29 15:34:23 UTC
  • Revision ID: git-v1:cdca5487c8322e426d0859349fa52643e0a47019
added python code from blitzortung-tracker-tools

Show diffs side-by-side

added added

removed removed

Lines of Context:
blitzortung/geom.py
 
1
# -*- coding: utf8 -*-
 
2
 
 
3
'''
 
4
 
 
5
@author Andreas Würl
 
6
 
 
7
'''
 
8
 
 
9
import subprocess
 
10
import re
 
11
import pyproj
 
12
import numpy
 
13
import math
 
14
import shapely
 
15
 
 
16
from abc import ABCMeta, abstractmethod
 
17
 
 
18
class Point(object):
 
19
 
 
20
  __geod = pyproj.Geod(ellps='WGS84', units='m')
 
21
 
 
22
  __whitespaceRe = re.compile('\s+')
 
23
 
 
24
  def __init__(self, x, y):
 
25
    self.x = x
 
26
    self.y = y
 
27
 
 
28
  def get_x(self):
 
29
    return self.x
 
30
 
 
31
  def get_y(self):
 
32
    return self.y
 
33
 
 
34
  def __invgeod(self, other):
 
35
    return Point.__geod.inv(self.x, self.y, other.x, other.y)
 
36
 
 
37
  def distance(self, other):
 
38
    return self.__invgeod(other)[2]
 
39
 
 
40
  def azimuth(self, other):
 
41
    return self.__invgeod(other)[0]
 
42
 
 
43
class Geometry:
 
44
  '''
 
45
  abstract base class for geometries
 
46
  '''
 
47
  
 
48
  __metaclass__ = ABCMeta
 
49
    
 
50
  DefaultSrid = 4326
 
51
  
 
52
  def __init__(self, srid = DefaultSrid):
 
53
    self.srid = srid
 
54
    
 
55
  def getSrid(self):
 
56
    return self.srid
 
57
  
 
58
  def setSrid(self, srid):
 
59
    self.srid = srid
 
60
 
 
61
  @abstractmethod
 
62
  def getEnv(self):
 
63
    pass
 
64
 
 
65
 
 
66
class Envelope(Geometry):
 
67
  ''' class for definition of coordinate envelopes '''
 
68
 
 
69
  def __init__(self, xmin, xmax, ymin, ymax, srid=Geometry.DefaultSrid):
 
70
    Geometry.__init__(self, srid)
 
71
    self.xmin = xmin
 
72
    self.xmax = xmax
 
73
    self.ymin = ymin
 
74
    self.ymax = ymax
 
75
 
 
76
  def getXMin(self):
 
77
    return self.xmin
 
78
 
 
79
  def getXMax(self):
 
80
    return self.xmax
 
81
 
 
82
  def getYMin(self):
 
83
    return self.ymin
 
84
 
 
85
  def getYMax(self):
 
86
    return self.ymax
 
87
  
 
88
  def getSrid(self):
 
89
    return self.srid
 
90
 
 
91
  def getDeltaY(self):
 
92
    return abs(self.ymax - self.ymin)
 
93
 
 
94
  def getDeltaX(self):
 
95
    return abs(self.xmax - self.xmin)
 
96
 
 
97
  def contains(self, point):
 
98
    if ((point.getX() >= self.xmin) and 
 
99
        (point.getX() <= self.xmax) and
 
100
        (point.getY() >= self.ymin) and
 
101
        (point.getY() <= self.ymax)):
 
102
      return True
 
103
    else:
 
104
      return False
 
105
    
 
106
  def getEnv(self):
 
107
    return shapely.geometry.Polygon(((self.xmin, self.ymin), (self.xmin, self.ymax), (self.xmax, self.ymax), (self.xmax, self.ymin)))
 
108
 
 
109
  def __str__(self):
 
110
    return 'longitude: %.2f .. %.2f, latitude: %.2f .. %.2f' % (self.xmin, self.xmax, self.ymin, self.ymax)
 
111
 
 
112
 
 
113
class Raster(Envelope):
 
114
  ' class for raster characteristics and data '
 
115
  
 
116
  def __init__(self, xmin, xmax, ymin, ymax, xdiv, ydiv, srid=Geometry.DefaultSrid, nodata = 0):
 
117
    Envelope.__init__(self, xmin, xmax, ymin, ymax, srid)
 
118
    self.xdiv = xdiv
 
119
    self.ydiv = ydiv
 
120
    self.nodata = nodata
 
121
    self.data = numpy.zeros((self.getYBinCount(), self.getXBinCount()), dtype=type(self.nodata))
 
122
 
 
123
  def getXDiv(self):
 
124
    return self.xdiv
 
125
  
 
126
  def getYDiv(self):
 
127
    return self.ydiv
 
128
  
 
129
  def getXBinCount(self):
 
130
    return int(math.ceil(1.0 * (self.xmax - self.xmin) / self.xdiv))
 
131
 
 
132
  def getYBinCount(self):
 
133
    return int(math.ceil(1.0 * (self.ymax - self.ymin) / self.ydiv))
 
134
  
 
135
  def set(self, x, y, data):
 
136
    self.data[y][x] = data
 
137
    
 
138
  def get(self, x, y):
 
139
    return self.data[y][x]
 
140
  
 
141
  def getNodataValue(self):
 
142
    return self.nodata
 
143
 
 
144
  def toArcGrid(self):
 
145
    result  = 'NCOLS %d\n' % self.getXBinCount()
 
146
    result += 'NROWS %d\n' % self.getYBinCount()
 
147
    result += 'XLLCORNER %.4f\n' % self.getXMin()
 
148
    result += 'YLLCORNER %.4f\n' % self.getYMin()
 
149
    result += 'CELLSIZE %.4f\n' % self.getXDiv()
 
150
    result += 'NODATA_VALUE %s\n' % str(self.getNodataValue())
 
151
    
 
152
    for row in self.data[::-1]:
 
153
      for cell in row:
 
154
        result += str(cell) + ' '
 
155
      result += '\n'
 
156
      
 
157
    return result.strip()
 
158
 
 
159
  def toMap(self):
 
160
    chars = " .-o*O8"
 
161
    maximum = 0
 
162
    total = 0
 
163
 
 
164
    for row in self.data[::-1]:
 
165
      for cell in row:
 
166
        total += cell
 
167
        if maximum < cell:
 
168
          maximum = cell
 
169
 
 
170
    if maximum > len(chars):
 
171
      divider = float(maximum) / (len(chars) - 1)
 
172
    else:
 
173
      divider = 1
 
174
 
 
175
    result = (self.getXBinCount() + 2) * '-' + '\n'
 
176
    for row in self.data[::-1]:
 
177
      result += "|"
 
178
      for cell in row:
 
179
        index = int(math.floor((cell - 1) / divider + 1))
 
180
        result += chars[index]
 
181
      result += "|\n"
 
182
 
 
183
    result += (self.getXBinCount() + 2) * '-' + '\n'
 
184
    result += 'total count: %d, max per area: %d' %(total, maximum)
 
185
    return result