GRASS Programmer's Manual  6.4.3(2013)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
profile.py
Go to the documentation of this file.
1 """!
2 @package wxplot.profile
3 
4 @brief Profiling using PyPlot
5 
6 Classes:
7  - profile::ProfileFrame
8  - profile::ProfileToolbar
9 
10 (C) 2011-2012 by the GRASS Development Team
11 
12 This program is free software under the GNU General Public License
13 (>=v2). Read the file COPYING that comes with GRASS for details.
14 
15 @author Michael Barton, Arizona State University
16 """
17 
18 import os
19 import sys
20 import math
21 
22 import wx
23 import wx.lib.plot as plot
24 
25 import grass.script as grass
26 
27 try:
28  import numpy
29 except ImportError:
30  msg = _("This module requires the NumPy module, which could not be "
31  "imported. It probably is not installed (it's not part of the "
32  "standard Python distribution). See the Numeric Python site "
33  "(http://numpy.scipy.org) for information on downloading source or "
34  "binaries.")
35  print >> sys.stderr, "wxplot.py: " + msg
36 
37 from wxplot.base import BasePlotFrame, PlotIcons
38 from gui_core.toolbars import BaseToolbar, BaseIcons
39 from wxplot.dialogs import ProfileRasterDialog, PlotStatsFrame
40 from core.gcmd import RunCommand, GWarning, GError, GMessage
41 
42 class ProfileFrame(BasePlotFrame):
43  """!Mainframe for displaying profile of one or more raster maps. Uses wx.lib.plot.
44  """
45  def __init__(self, parent, id = wx.ID_ANY, style = wx.DEFAULT_FRAME_STYLE,
46  size = wx.Size(700, 400),
47  rasterList = [], **kwargs):
48  BasePlotFrame.__init__(self, parent, size = size, **kwargs)
49 
50  self.toolbar = ProfileToolbar(parent = self)
51  self.SetToolBar(self.toolbar)
52  self.SetTitle(_("GRASS Profile Analysis Tool"))
53 
54  #
55  # Init variables
56  #
57  self.rasterList = rasterList
58  self.plottype = 'profile'
59  self.coordstr = '' # string of coordinates for r.profile
60  self.seglist = [] # segment endpoint list
61  self.ppoints = '' # segment endpoints data
62  self.transect_length = 0.0 # total transect length
63  self.ptitle = _('Profile of') # title of window
64  self.colorList = ["blue", "red", "green", "yellow", "magenta", "cyan",
65  "aqua", "black", "grey", "orange", "brown", "purple", "violet",
66  "indigo"]
67 
68  self._initOpts()
69 
70  if len(self.rasterList) > 0: # set raster name(s) from layer manager if a map is selected
71  self.raster = self.InitRasterOpts(self.rasterList, self.plottype)
72  else:
73  self.raster = {}
74 
75  # determine units (axis labels)
76  if self.parent.Map.projinfo['units'] != '':
77  self.xlabel = _('Distance (%s)') % self.parent.Map.projinfo['units']
78  else:
79  self.xlabel = _("Distance along transect")
80  self.ylabel = _("Cell values")
81 
82  def _initOpts(self):
83  """!Initialize plot options
84  """
85  self.InitPlotOpts('profile')
86 
87  def OnDrawTransect(self, event):
88  """!Draws transect to profile in map display
89  """
90  self.mapwin.polycoords = []
91  self.seglist = []
92  self.mapwin.ClearLines(self.mapwin.pdc)
93  self.ppoints = ''
94 
95  self.parent.SetFocus()
96  self.parent.Raise()
97 
98  self.mapwin.mouse['use'] = 'profile'
99  self.mapwin.mouse['box'] = 'line'
100  self.mapwin.pen = wx.Pen(colour = 'Red', width = 2, style = wx.SHORT_DASH)
101  self.mapwin.polypen = wx.Pen(colour = 'dark green', width = 2, style = wx.SHORT_DASH)
102  self.mapwin.SetCursor(self.Parent.cursors["cross"])
103 
104  def OnSelectRaster(self, event):
105  """!Select raster map(s) to profile
106  """
107  dlg = ProfileRasterDialog(parent = self)
108 
109  if dlg.ShowModal() == wx.ID_OK:
110  self.rasterList = dlg.rasterList
111  self.raster = self.InitRasterOpts(self.rasterList, self.plottype)
112 
113  # plot profile
114  if len(self.mapwin.polycoords) > 0 and len(self.rasterList) > 0:
115  self.OnCreateProfile(event = None)
116 
117  dlg.Destroy()
118 
119  def SetupProfile(self):
120  """!Create coordinate string for profiling. Create segment list for
121  transect segment markers.
122  """
123 
124  #
125  # create list of coordinate points for r.profile
126  #
127  dist = 0
128  cumdist = 0
129  self.coordstr = ''
130  lasteast = lastnorth = None
131 
132  region = grass.region()
133  insideRegion = True
134  if len(self.mapwin.polycoords) > 0:
135  for point in self.mapwin.polycoords:
136  if not (region['w'] <= point[0] <= region['e'] and region['s'] <= point[1] <= region['n']):
137  insideRegion = False
138  # build string of coordinate points for r.profile
139  if self.coordstr == '':
140  self.coordstr = '%d,%d' % (point[0], point[1])
141  else:
142  self.coordstr = '%s,%d,%d' % (self.coordstr, point[0], point[1])
143 
144  if not insideRegion:
145  GWarning(message = _("Not all points of profile lie inside computational region."),
146  parent = self)
147 
148  if len(self.rasterList) == 0:
149  return
150 
151  # title of window
152  self.ptitle = _('Profile of')
153 
154  #
155  # create list of coordinates for transect segment markers
156  #
157  if len(self.mapwin.polycoords) > 0:
158  self.seglist = []
159  for point in self.mapwin.polycoords:
160  # get value of raster cell at coordinate point
161  ret = RunCommand('r.what',
162  parent = self,
163  read = True,
164  input = self.rasterList[0],
165  east_north = '%d,%d' % (point[0],point[1]))
166 
167  val = ret.splitlines()[0].split('|')[3]
168  if val == None or val == '*': continue
169  val = float(val)
170 
171  # calculate distance between coordinate points
172  if lasteast and lastnorth:
173  dist = math.sqrt(math.pow((lasteast-point[0]),2) + math.pow((lastnorth-point[1]),2))
174  cumdist += dist
175 
176  #store total transect length
177  self.transect_length = cumdist
178 
179  # build a list of distance,value pairs for each segment of transect
180  self.seglist.append((cumdist,val))
181  lasteast = point[0]
182  lastnorth = point[1]
183 
184  # delete extra first segment point
185  try:
186  self.seglist.pop(0)
187  except:
188  pass
189 
190  #
191  # create datalist of dist/value pairs and y labels for each raster map
192  #
193  self.ylabel = ''
194  i = 0
195 
196  for r in self.raster.iterkeys():
197  self.raster[r]['datalist'] = []
198  datalist = self.CreateDatalist(r, self.coordstr)
199  if len(datalist) > 0:
200  self.raster[r]['datalist'] = datalist
201 
202  # update ylabel to match units if they exist
203  if self.raster[r]['units'] != '':
204  self.ylabel += '%s (%d),' % (r['units'], i)
205  i += 1
206 
207  # update title
208  self.ptitle += ' %s ,' % r.split('@')[0]
209 
210  self.ptitle = self.ptitle.rstrip(',')
211 
212  if self.ylabel == '':
213  self.ylabel = _('Raster values')
214  else:
215  self.ylabel = self.ylabel.rstrip(',')
216 
217  def CreateDatalist(self, raster, coords):
218  """!Build a list of distance, value pairs for points along transect using r.profile
219  """
220  datalist = []
221 
222  # keep total number of transect points to 500 or less to avoid
223  # freezing with large, high resolution maps
224  region = grass.region()
225  curr_res = min(float(region['nsres']),float(region['ewres']))
226  transect_rec = 0
227  if self.transect_length / curr_res > 500:
228  transect_res = self.transect_length / 500
229  else: transect_res = curr_res
230 
231  ret = RunCommand("r.profile",
232  parent = self,
233  input = raster,
234  profile = coords,
235  res = transect_res,
236  null = "nan",
237  quiet = True,
238  read = True)
239 
240  if not ret:
241  return []
242 
243  for line in ret.splitlines():
244  dist, elev = line.strip().split(' ')
245  if dist == None or dist == '' or dist == 'nan' or \
246  elev == None or elev == '' or elev == 'nan':
247  continue
248  dist = float(dist)
249  elev = float(elev)
250  datalist.append((dist,elev))
251 
252  return datalist
253 
254  def OnCreateProfile(self, event):
255  """!Main routine for creating a profile. Uses r.profile to
256  create a list of distance,cell value pairs. This is passed to
257  plot to create a line graph of the profile. If the profile
258  transect is in multiple segments, these are drawn as
259  points. Profile transect is drawn, using methods in mapdisp.py
260  """
261 
262  if len(self.mapwin.polycoords) == 0 or len(self.rasterList) == 0:
263  dlg = wx.MessageDialog(parent = self,
264  message = _('You must draw a transect to profile in the map display window.'),
265  caption = _('Nothing to profile'),
266  style = wx.OK | wx.ICON_INFORMATION | wx.CENTRE)
267  dlg.ShowModal()
268  dlg.Destroy()
269  return
270 
271  self.mapwin.SetCursor(self.parent.cursors["default"])
272  self.SetCursor(self.parent.cursors["default"])
273  self.SetGraphStyle()
274  self.SetupProfile()
275  p = self.CreatePlotList()
276  self.DrawPlot(p)
277 
278  # reset transect
279  self.mapwin.mouse['begin'] = self.mapwin.mouse['end'] = (0.0,0.0)
280  self.mapwin.mouse['use'] = 'pointer'
281  self.mapwin.mouse['box'] = 'point'
282 
283  def CreatePlotList(self):
284  """!Create a plot data list from transect datalist and
285  transect segment endpoint coordinates.
286  """
287  # graph the distance, value pairs for the transect
288  self.plotlist = []
289 
290  # Add segment marker points to plot data list
291  if len(self.seglist) > 0 :
292  self.ppoints = plot.PolyMarker(self.seglist,
293  legend = ' ' + self.properties['marker']['legend'],
294  colour = wx.Color(self.properties['marker']['color'][0],
295  self.properties['marker']['color'][1],
296  self.properties['marker']['color'][2],
297  255),
298  size = self.properties['marker']['size'],
299  fillstyle = self.ptfilldict[self.properties['marker']['fill']],
300  marker = self.properties['marker']['type'])
301  self.plotlist.append(self.ppoints)
302 
303  # Add profile distance/elevation pairs to plot data list for each raster profiled
304  for r in self.rasterList:
305  col = wx.Color(self.raster[r]['pcolor'][0],
306  self.raster[r]['pcolor'][1],
307  self.raster[r]['pcolor'][2],
308  255)
309  self.raster[r]['pline'] = plot.PolyLine(self.raster[r]['datalist'],
310  colour = col,
311  width = self.raster[r]['pwidth'],
312  style = self.linestyledict[self.raster[r]['pstyle']],
313  legend = self.raster[r]['plegend'])
314 
315  self.plotlist.append(self.raster[r]['pline'])
316 
317  if len(self.plotlist) > 0:
318  return self.plotlist
319  else:
320  return None
321 
322  def Update(self):
323  """!Update profile after changing options
324  """
325  self.SetGraphStyle()
326  p = self.CreatePlotList()
327  self.DrawPlot(p)
328 
329  def SaveProfileToFile(self, event):
330  """!Save r.profile data to a csv file
331  """
332  dlg = wx.FileDialog(parent = self,
333  message = _("Choose prefix for file(s) where to save profile values..."),
334  defaultDir = os.getcwd(),
335  wildcard = _("Comma separated value (*.csv)|*.csv"), style = wx.SAVE)
336  pfile = []
337  if dlg.ShowModal() == wx.ID_OK:
338  path = dlg.GetPath()
339  for r in self.rasterList:
340  pfile.append(path + '_' + str(r.replace('@', '_')) + '.csv')
341  if os.path.exists(pfile[-1]):
342  dlgOv = wx.MessageDialog(self,
343  message = _("File <%s> already exists. "
344  "Do you want to overwrite this file?") % pfile[-1],
345  caption = _("Overwrite file?"),
346  style = wx.YES_NO | wx.YES_DEFAULT | wx.ICON_QUESTION)
347  if dlgOv.ShowModal() != wx.ID_YES:
348  pfile.pop()
349  dlgOv.Destroy()
350  continue
351 
352  try:
353  fd = open(pfile[-1], "w")
354  except IOError, e:
355  GError(parent = self,
356  message = _("Unable to open file <%(file)s> for "
357  "writing.\nReason: %(e)s") % {'file': pfile[-1],
358  'e': e})
359  dlg.Destroy()
360  return
361 
362  for datapair in self.raster[r]['datalist']:
363  fd.write('%d,%d\n' % (float(datapair[0]),float(datapair[1])))
364 
365  fd.close()
366 
367  dlg.Destroy()
368  if pfile:
369  message = _("%(l)d files created:\n%(p)s") % {'l': len(pfile),
370  'p':'\n'.join(pfile)}
371  else:
372  message = _("No files generated.")
373 
374  GMessage(parent = self, message = message)
375 
376  def OnStats(self, event):
377  """!Displays regression information in messagebox
378  """
379  message = []
380  title = _('Statistics for Profile(s)')
381 
382  for r in self.raster.iterkeys():
383  try:
384  rast = r.split('@')[0]
385  statstr = 'Profile of %s\n\n' % rast
386 
387  iterable = (i[1] for i in self.raster[r]['datalist'])
388  a = numpy.fromiter(iterable, numpy.float)
389 
390  statstr += 'n: %f\n' % a.size
391  statstr += 'minimum: %f\n' % numpy.amin(a)
392  statstr += 'maximum: %f\n' % numpy.amax(a)
393  statstr += 'range: %f\n' % numpy.ptp(a)
394  statstr += 'mean: %f\n' % numpy.mean(a)
395  statstr += 'standard deviation: %f\n' % numpy.std(a)
396  statstr += 'variance: %f\n' % numpy.var(a)
397  cv = numpy.std(a)/numpy.mean(a)
398  statstr += 'coefficient of variation: %f\n' % cv
399  statstr += 'sum: %f\n' % numpy.sum(a)
400  statstr += 'median: %f\n' % numpy.median(a)
401  statstr += 'distance along transect: %f\n\n' % self.transect_length
402  message.append(statstr)
403  except:
404  pass
405 
406  stats = PlotStatsFrame(self, id = wx.ID_ANY, message = message,
407  title = title)
408 
409  if stats.Show() == wx.ID_CLOSE:
410  stats.Destroy()
411 
412 
413 class ProfileToolbar(BaseToolbar):
414  """!Toolbar for profiling raster map
415  """
416  def __init__(self, parent):
417  BaseToolbar.__init__(self, parent)
418 
419  self.InitToolbar(self._toolbarData())
420 
421  # realize the toolbar
422  self.Realize()
423 
424  def _toolbarData(self):
425  """!Toolbar data"""
426  return self._getToolbarData((('addraster', BaseIcons["addRast"],
427  self.parent.OnSelectRaster),
428  ('transect', PlotIcons["transect"],
429  self.parent.OnDrawTransect),
430  (None, ),
431  ('draw', PlotIcons["draw"],
432  self.parent.OnCreateProfile),
433  ('erase', BaseIcons["erase"],
434  self.parent.OnErase),
435  ('drag', BaseIcons['pan'],
436  self.parent.OnDrag),
437  ('zoom', BaseIcons['zoomIn'],
438  self.parent.OnZoom),
439  ('unzoom', BaseIcons['zoomBack'],
440  self.parent.OnRedraw),
441  (None, ),
442  ('statistics', PlotIcons['statistics'],
443  self.parent.OnStats),
444  ('datasave', PlotIcons["save"],
445  self.parent.SaveProfileToFile),
446  ('image', BaseIcons["saveFile"],
447  self.parent.SaveToFile),
448  ('print', BaseIcons["print"],
449  self.parent.PrintMenu),
450  (None, ),
451  ('settings', PlotIcons["options"],
452  self.parent.PlotOptionsMenu),
453  ('quit', PlotIcons["quit"],
454  self.parent.OnQuit),
455  ))